Extending the solver to a new PDE

Once you have a Poisson solver, moving to a different linear elliptic PDE means changing two methods and nothing else. This example adds a reaction term -- a modified Helmholtz problem \(\nabla^2 u - k\,u = f\) -- which is the shape of a screened-Poisson or implicit-diffusion step.

The code block is a region of tests/snippets/poisson_tutorial.cpp, compiled and run with the test suite.

Change the operator and the smoother

// A modified Helmholtz problem:  laplacian(u) - k*u = f.
// Only the operator and the smoother change; the boundary handling, the cycle
// and the I/O are all inherited unchanged.
class ModifiedHelmholtz : public mgrid::LinearMultigrid {
  public:
    ModifiedHelmholtz(const mgrid::Settings &settings, double k)
        : mgrid::LinearMultigrid(settings), k_(k) {
        for (auto edge : mgrid::allBoundaryFlags)
            solution.boundaryConditions.set(edge, mgrid::zeroDirichletCondition);
        solution.propagate_boundary_conditions();

        source_term() = -1.0;
        mark_source_set();
    }

    double differential_operator(mgrid::Level level, int i, int j) override {
        return solution[level].dxx(i, j) + solution[level].dzz(i, j) - k_ * solution[level](i, j);
    }

    void relaxation_updater(mgrid::Level level, int i, int j) override {
        const double hx = solution[level].spacing(0);
        const double hz = solution[level].spacing(1);
        const double cx = 1.0 / (hx * hx);
        const double cz = 1.0 / (hz * hz);
        // The reaction term adds k to the diagonal coefficient.
        solution[level](i, j) =
            ((solution[level](i + 1, j) + solution[level](i - 1, j)) * cx +
             (solution[level](i, j + 1) + solution[level](i, j - 1)) * cz - source[level](i, j)) /
            (2.0 * (cx + cz) + k_);
    }

    std::string filename(std::string root = "") override { return root + "helmholtz"; }

  private:
    double k_;
};

Two changes from the Poisson class:

1. differential_operator gains the \(-k\,u\) term: ... - k_ * solution[level](i, j).
2. relaxation_updater must invert the new diagonal. Discretising \(\nabla^2 u - k u = f\) and solving for \(u_{i,j}\) puts \(k\) into the denominator:

$$ u_{i,j} = \frac{c_x(u_{i+1,j} + u_{i-1,j}) + c_z(u_{i,j+1} + u_{i,j-1}) - f_{i,j}} {2(c_x + c_z) + k}, \qquad c_x = 1/h_x^2,\ c_z = 1/h_z^2. $$

The operator and the smoother must always agree: relaxation_updater is just differential_operator solved for the centre point.

Everything else -- boundary handling, propagate_boundary_conditions(), mark_source_set(), the FMG cycle, write() -- is unchanged. Here the constructor uses all-Dirichlet edges via mgrid::allBoundaryFlags:

for (auto edge : mgrid::allBoundaryFlags)
    solution.boundaryConditions.set(edge, mgrid::zeroDirichletCondition);

Solve it

    mgrid::Settings settings;
    settings.numberOfGrids = 5;
    settings.minimumResolution = 4;

    ModifiedHelmholtz problem(settings, /*k=*/10.0);
    problem.solve();
    mgrid::FDArray &u = problem.get_result();

With \(k > 0\) the operator is more diagonally dominant than plain Poisson, so Gauss--Seidel and the multigrid cycle both converge at least as fast. Far from the boundary the Laplacian term is small and \(-k\,u \approx -1\), so the interior plateau sits near \(1/k\).

What each override is for

Method	Purpose	Must satisfy
differential_operator(level, i, j)	evaluate \(\mathcal{L}u\) at a point	uses solution[level] stencils; may read scratch state (it is non-const)
relaxation_updater(level, i, j)	one in-place point smooth	writes solution[level](i,j); must be differential_operator solved for the centre
filename(root)	output file stem	any string; the examples encode their parameters

Nonlinear problems

For a genuinely nonlinear operator, derive from mgrid::NonlinearMultigrid instead -- same three overrides, but the base runs an FAS cycle. The mosolov example takes the other route common for viscoplastic flow: it stays linear per step and wraps an augmented-Lagrangian iteration around LinearMultigrid.