A nonlinear PDE (the Bratu problem)

The Poisson and modified-Helmholtz examples are linear: the operator is a fixed matrix. This one is genuinely nonlinear, so it uses mgrid::NonlinearMultigrid, whose base runs a full-approximation-scheme (FAS) cycle rather than a linear correction scheme.

The Bratu problem

\[ \nabla^2 u + \lambda\,e^{u} = f \]

is the standard nonlinear-elliptic multigrid benchmark. Every code block below is a region of tests/snippets/nonlinear_tutorial.cpp, compiled and run as part of the test suite (against a manufactured solution $u^\* = \sin \pi x \sin \pi z$ on the unit square, so the code can check its own accuracy).

Subclass `NonlinearMultigrid`

// The Bratu problem  laplacian(u) + lambda * e^u = f  on the unit square, with
// zero-Dirichlet data on every edge. The e^u term makes it nonlinear, so it
// derives from NonlinearMultigrid: the base runs a full-approximation-scheme
// (FAS) cycle, and this class supplies the discrete operator together with a
// *nonlinear* pointwise smoother.
class Bratu : public mgrid::NonlinearMultigrid {
  public:
    explicit Bratu(const mgrid::Settings &settings) : mgrid::NonlinearMultigrid(settings) {
        for (auto edge : mgrid::allBoundaryFlags)
            solution.boundaryConditions.set(edge, mgrid::zeroDirichletCondition);
        solution.propagate_boundary_conditions();

        // Sample the manufactured right-hand side onto the finest grid.
        auto &f = source_term();
        for (int i = 0; i < f.rows(); ++i)
            for (int j = 0; j < f.cols(); ++j)
                f(i, j) = bratu_rhs(i * f.spacing(0), j * f.spacing(1));
        mark_source_set();
    }

    // L(u) at one point: the 5-point Laplacian (via the FDArray stencils, which
    // stay in bounds on the edges) plus the nonlinear reaction term.
    double differential_operator(mgrid::Level l, int i, int j) override {
        return solution[l].dxx(i, j) + solution[l].dzz(i, j) +
               kLambda * std::exp(solution[l](i, j));
    }

    // One nonlinear Gauss-Seidel update: two Newton steps for w = u(i,j) in
    //   phi(w) = C*w - lambda*e^w - (S - f_ij) = 0,
    // where C is the centre coefficient and S the neighbour sum of the stencil.
    void relaxation_updater(mgrid::Level l, int i, int j) override {
        auto &u = solution[l];
        const double cx = 1.0 / (u.spacing(0) * u.spacing(0));
        const double cz = 1.0 / (u.spacing(1) * u.spacing(1));
        const double c = 2.0 * (cx + cz);
        const double s = (u(i + 1, j) + u(i - 1, j)) * cx + (u(i, j + 1) + u(i, j - 1)) * cz;
        const double rhs_ij = s - source[l](i, j);

        double w = u(i, j);
        for (int step = 0; step < 2; ++step) {
            const double phi = c * w - kLambda * std::exp(w) - rhs_ij;
            const double dphi = c - kLambda * std::exp(w);
            w -= phi / dphi;
        }
        u(i, j) = w;
    }

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

    // Frobenius norm of the nonlinear residual  f - L(u)  on the finest grid.
    double residual_norm() {
        mgrid::FDArray r(solution[finestLevel]); // copy: same shape and geometry
        evaluate_residual(finestLevel, r);
        return mgrid::frobenius_norm(r.field());
    }
};

The two overrides carry the nonlinearity:

differential_operator returns the full nonlinear residual operator $\mathcal{L}(u) = \nabla^2_h u + \lambda e^{u}$ at a point. The FAS cycle evaluates this on both the fine and the restricted-coarse solution to build the tau correction that keeps the coarse problem consistent.
relaxation_updater is a nonlinear Gauss--Seidel step. The pointwise equation

$$ c\,u_{i,j} - \lambda e^{u_{i,j}} = S - f_{i,j}, \qquad c = \tfrac{2}{h_x^2} + \tfrac{2}{h_z^2}, $$

($S$ = the stencil's neighbour sum) is transcendental in $u_{i,j}$, so it is solved with two Newton steps per visit. For a linear operator this reduces to the ordinary weighted-Jacobi/GS update.

source, solution, mark_source_set() and the boundary setup are exactly as in the linear tutorial.

Solve

    mgrid::Settings settings;
    settings.aspectRatio = 1.0; // square cells on the unit square
    settings.numberOfGrids = 5;
    settings.minimumResolution = 4;

    Bratu problem(settings);
    problem.solve(); // FAS full-multigrid cycle

    mgrid::FDArray &u = problem.get_result();

solve() runs one FAS full-multigrid cycle. For $\lambda = 1$ (well below the turning point $\lambda_c \approx 6.81$) the nonlinearity is mild and one FMG pass reaches truncation-error accuracy: the test measures second-order convergence, $8\times 10^{-4}$ at $49\times 49$ down to $9\times 10^{-5}$ at $97\times 97$.

FAS vs. the linear correction scheme

For a linear operator, FAS and the linear coarse-grid correction are algebraically identical and either transfer works for the solution. For a nonlinear operator they are not: FAS must restrict the solution by straight injection ($u_{2h}(i,j) = u_h(2i,2j)$), not full weighting. Full-weighting the solution smooths it and moves the FAS fixed point off the true discrete solution -- a mesh-independent error of a few percent that never converges away. NonlinearMultigrid does this correctly; tests/test_multigrid_nonlinear.cpp pins it with a manufactured Bratu solve at second order.

Harder nonlinearities

For $\lambda$ near $\lambda_c$, or for a stiffer nonlinearity, one FMG pass is not enough -- wrap the solve in an outer loop of FAS V-cycles until the residual stops dropping, or switch to a globalised Newton iteration where each step is a linear LinearMultigrid solve (the pattern the Mosolov example uses for viscoplastic flow).

A nonlinear PDE (the Bratu problem)

Subclass NonlinearMultigrid

Solve

FAS vs. the linear correction scheme

Harder nonlinearities

Subclass `NonlinearMultigrid`