Back to the main page Back to the documentation 1. Overview: Multigrid solver 2. Example: Poisson equation 3. Methods to help defining the equation 4. More examples

This is a multigrid solver for very general PDEs: it can handle any type (combination of types), scalar, vectors whatever. The only restriction is that the solver type class is a ring type (i.e. it implements all the arithmetic operations +,-,* and multiplication/division with a scalar). What the solver does is to first make a guess for the solution and then iterate over the grid and update the solution by some root finding method (we use Newtons method). To improve the convergence we apply the multigrid technique: the solution is restricted down to a coarser grid and we repeat the steps above. We do this recursively down to the coarsest level (typically $2$ cells per dimension). Then we restrict the solution up and use this to correct the solution at the finer level and do this recursively up to the finest grid. One such step is called a V-cycle. We repeat this until convergence. Between every sweep of the grid where boundary cells are modified a communication is performed. For more info about how the multigrid solver works see for example Section 3.2 of Winther et al. 2015.

Its hairy doing this in parallel if we allow any number of gridcells. The solution is that we place a few restrictions on the user: the number of tasks $N_{\rm Tasks}$ must divide the gridsize $N_{\rm grid}$ and $\frac{N_{\rm grid}}{N_{\rm Tasks}}$ must be divisible by $2^M$ where $M$ is the number of levels we need. These restrictions (which are not very restrictive) is enough to ensure that there are no coarse cells split between tasks and we don't have to do any communication when restricting down and prolonging up between levels. The coarsest level we can have is one where there are $1$ cell per MPI tasks. If we have a large number of tasks then the convergence will not be as fast as it would be if we could go all the way down. If this is an issue then use fewer tasks if possible (and compansate with more threads, i.e. use a hybrid OpenMP + MPI).

Figure 1: Sketch of how the multigrid solver works.

Below is a simple example for how to solve the PDE $$\nabla^2 \phi(x) = S(x)$$ We write this as $\mathcal{L} = \nabla^2\phi - S$ thus we want to find $\phi$ such that $\mathcal{L} = 0$. The equation can be discretized (here in 1D for simplicity) as $$\mathcal{L}_i = \frac{\phi_{i+1} + \phi_{i-1} - 2\phi_{i}}{h^2} - S_i$$ We also need the derivative of this function $$\frac{d\mathcal{L}}{d\phi_i} = -\frac{2}{h^2}$$ In general dimensions this would be a sum over the $xx$, $yy$, etc. derivatives and the derivative of $\mathcal{L}$ would be $\frac{-2{\rm NDIM}}{h^2}$. These two functions are all we need to implement in the solver.

Set up the solver:

// Solver parameters
using SolverType     = double;  // The type of phi 
const int N          = 128;     // Number of cells
const int Nlevels    = -1;      // Number of multigrid level (-1 let the solver pick this)
const int NDIM       = 2;       // Dimension we are working in
const int nleft      = 1;       // Extra slices (we need atleast 1 for derivatives at the boundary)
const int nright     = 1;       //  --- || ---
const bool periodic  = true;    // Periodic box?
const bool verbose   = true;    

MultiGridSolver<NDIM,SolverType> mg(N, Nlevels, verbose, periodic, nleft, nright);

Define the equation and the convergence criterion you want:

// The right hand side S: this is here just D^2 of phi = Prod_i sin(2 pi x_i)
auto source = [=](std::array<double,NDIM> & x) -> SolverType {
    SolverType sol = (-4.0 * M_PI * M_PI);
    for(int idim = 0; idim < NDIM; idim++)
      sol *= sin(2.0 * M_PI * x[idim]);
    return sol;
}

// Define the equation: we here compute and return L and dLdy
// This function gets called for every cell and index is the index of the given cell
MultiGridFunction<NDIM,SolverType> Equation = [&](MultiGridSolver<NDIM,SolverType> *sol, int level, IndexInt index){
  auto index_list = sol->get_neighbor_gridindex(level, index);
  auto coordinate = sol->get_Coordinate(level, index);
  auto L  = sol->get_Laplacian(level, index_list) - source(coordinate);
  auto dL = sol->get_derivLaplacian(level, index_list);
  return std::pair<SolverType,SolverType>{L, dL};
};

// Set the convergence criterion: here we just use the fiducial choice 
// (ConvergenceCriterionResidual is rms-residual < epsilon)
MultiGridConvCrit ConvergenceCriterion = [&](double rms_residual, double rms_residual_ini, int step_number){
  return mg.ConvergenceCriterionResidual(rms_residual, rms_residual_ini, step_number);
};

Set the initial guess, run the solver and fetch the solution:

// Set the initial guess
mg.set_initial_guess( SolverType(0.0) );

// Run the solver
mg.solve(Equation, ConvergenceCriterion);

// Get the solution
auto sol = mg.get_grid(0);

In this example its very simple as we could use built in functions to define the equation. With the given source the solution should just be $\phi(x_1,\ldots) = \prod_{i=1}^{\rm NDIM} \sin(2\pi x_i)$.

We provide some functions to make it more easy to define the equation. If this is not enough then you will have to do it yourself. There are much more information in the headerfile defining the solver.

// The solution in a given cell
T get_Field(int level, IndexInt index);

// Position of the cell in the global grid
std::array<double,NDIM>   get_Coordinate(int level, IndexInt index);

// The closest 2NDIM cells (plus the cell itself). This is what the methods below require
std::array<IndexInt,2*NDIM+1> get_neighbor_gridindex(int level, IndexInt index);

// The gradent Df and its derivative
std::array<T,NDIM>  get_Gradient     (int level, const std::array<IndexInt,2*NDIM+1> & nbor_index_list);
std::array<T,NDIM>  get_derivGradient(int level, const std::array<IndexInt,2*NDIM+1> & nbor_index_list);

// The Laplacian operator D^2f and its derivative
T get_Laplacian     (int level, const std::array<IndexInt,2*NDIM+1> & nbor_index_list);
T get_derivLaplacian(int level, const std::array<IndexInt,2*NDIM+1> & nbor_index_list);

// If you need even more cells: get all 3^NDIM cells around a given cell
std::array<IndexInt, FML::power(3,NDIM) > get_cube_gridindex(int level, IndexInt index);

In the source code on github we provide a lot more examples like: Poisson equation: $$\nabla^2\phi = \delta$$ Density reconstruction equation for the displacement field $\Psi$: $$\nabla\cdot \Psi + \frac{f}{b}\nabla\cdot[(\Psi\cdot r)r] = -\frac{\delta}{b}$$ The symmetron model equation: $$\nabla^2\phi = \frac{1}{2}\frac{B^2}{L^2}\left[ (1+\delta)\left(\frac{a_{\rm ssb}}{a}\right)^3\phi - \phi + \phi^3\right]$$ Hu-Sawicky $f(R)$ gravity equation: $$\nabla^2 f_R = A\delta + B \frac{1}{f_R^{1/(n+1)}}$$ Quasi-linear continuity equation for the velocity field: $$fH\delta + \nabla\cdot((1+\delta)\vec{v}) = 0$$ DGP equation: $$\nabla^2 \phi + A[(\nabla^2\phi)^2 - (\nabla_i\nabla_j\phi)^2] = B\delta$$ And the qubic variant: $$(\nabla^2 \phi)^3 + p (\nabla^2 \phi) + q = 0$$ And examples of how to using your own types with the solver.