import numpy as np
import matplotlib.pyplot as plt

Restriction and Prolongation

We need a way to transfer data from the current grid to a coarser grid (restriction) or finer grid (prolongation).

The procedure by which we do this depends on the type of grid we are using. We’ll continue to consider a cell-centered finite-difference / finite-volume grid.

Since the finite-volume method is conservative, we want to ensure that the operations that we do here maintain that conservation.

Consider the following 2 grids:

The fine grid (top) is indexed with \(i\) and the coarse grid (bottom) is indexing with \(j\).

Note that the grids are properly nested: two fine cells fit exactly into a single coarse cell.

In general, the jump in resolution can be something other than 2, but we’ll restrict ourselves to jumps by 2 here.

Restriction

Consider the case of bringing the data from the fine grid to the coarse grid. Since two fine cells are contained exactly in a single coarse cell, a conservative aveaging of the fine data down to the coarse grid would be:

\[\phi^c_j = \frac{1}{2} \left ( \phi^f_i + \phi^f_{i+1} \right )\]

This is also second-order accurate at the coarse cell cell-center.

Prolongation

To initialize the two fine cells from a coarse cell, we could do simple direct injection:

\[\phi_i^f = \phi_j^c\]

\[\phi_{i+1}^f = \phi_j^c\]

This is only first order accurate. A second-order accurate method of initializing the fine data starts with reconstructing the coarse data as a line:

\[\phi^{\mathrm{coarse}}(x) = m_x (x - x_j) + \phi^c_j\]

with

\[m_x = \frac{\phi^c_{j+1} - \phi^c_{j-1}}{2\Delta x}\]

Then we integrate this line over the two fine cells to find their values:

\[\phi^f_j = \phi^c_i - \frac{\Delta x}{4} m_x\]

\[\phi^f_{j+1} = \phi^c_i + \frac{\Delta x}{4} m_x\]

Implementation

Let’s create a basic grid class that implements this. We’ll assume homogeneous Dirichlet boundary conditions.

class Grid:
    def __init__(self, nx, ng=1, xmin=0, xmax=1):

        self.xmin = xmin
        self.xmax = xmax
        self.ng = ng
        self.nx = nx

        self.ilo = ng
        self.ihi = ng+nx-1

        # physical coords -- cell-centered
        self.dx = (xmax - xmin)/(nx)
        self.x = xmin + (np.arange(nx+2*ng)-ng+0.5)*self.dx

        # storage for the solution
        self.phi = self.scratch_array()

    def scratch_array(self):
        """return a scratch array dimensioned for our grid """
        return np.zeros((self.nx+2*self.ng), dtype=np.float64)

    def fill_bcs(self):
        """fill the boundary conditions on phi -- we'll assume Dirichlet"""

        self.phi[self.ilo-1] = -self.phi[self.ilo]
        self.phi[self.ihi+1] = -self.phi[self.ihi]

    def restrict(self):
        """conservative restriction"""
        
        # create a coarse array
        ng = self.ng
        nc = self.nx//2

        ilo_c = ng
        ihi_c = ng + nc - 1

        coarse_data = np.zeros((nc + 2*ng), dtype=np.float64)

        coarse_data[ilo_c:ihi_c+1] = 0.5 * (self.phi[self.ilo:self.ihi+1:2] +
                                            self.phi[self.ilo+1:self.ihi+1:2])

        return coarse_data

    def prolong(self):
        """
        prolong the data in the current (coarse) grid to a finer
        (factor of 2 finer) grid using linear reconstruction.
        """

        # allocate an array for the coarsely gridded data
        ng = self.ng
        nf = self.nx * 2

        fine_data = np.zeros((nf + 2*ng), dtype=np.float64)

        ilo_f = ng
        ihi_f = ng + nf - 1

        # slopes for the coarse data
        m_x = self.scratch_array()
        m_x[self.ilo:self.ihi+1] = 0.5 * (self.phi[self.ilo+1:self.ihi+2] -
                                          self.phi[self.ilo-1:self.ihi])

        # fill the left children
        fine_data[ilo_f:ihi_f+1:2] = \
            self.phi[self.ilo:self.ihi+1] - 0.25 * m_x[self.ilo:self.ihi+1]

        # fill the right children
        fine_data[ilo_f+1:ihi_f+1:2] = \
            self.phi[self.ilo:self.ihi+1] + 0.25 * m_x[self.ilo:self.ihi+1]

        return fine_data

Let’s test it out. If we are conservative, we should be able to prolong and then restrict some initial data and get back what we started with.

We’ll create 2 grid objects: a coarse and fine one

nx = 6
coarse = Grid(nx)
fine = Grid(2*nx)

Now initialize the coarse data to just a sequence of integers

coarse.phi[coarse.ilo:coarse.ihi+1] = np.arange(nx)+1
coarse.fill_bcs()
coarse.phi

array([-1.,  1.,  2.,  3.,  4.,  5.,  6., -6.])

Notice that because of this choice of initial data and the homogeneous Dirichlet BC on the right, the data is very steep there. That should not matter for this though—we should still be conservative.

Let’s initialize the fine data

fine.phi[:] = coarse.prolong()
fine.fill_bcs()
fine.phi

array([-0.625,  0.625,  1.375,  1.75 ,  2.25 ,  2.75 ,  3.25 ,  3.75 ,
        4.25 ,  4.75 ,  5.25 ,  7.375,  4.625, -4.625])

Now let’s restrict this and compare to the original coarse data

coarse_new = fine.restrict()
coarse_new

array([0., 1., 2., 3., 4., 5., 6., 0.])

Notice that this is the same as our original data—our prolongation and restriction operators are conservative.