libmultigrid

A finite-difference multigrid solver for elliptic PDEs on rectangular domains.

The library discretises a second-order elliptic PDE on a rectangular grid and drives it to its discrete solution with a geometric multigrid cycle: red-black Gauss-Seidel smoothing, full-weighting restriction and bilinear interpolation for transferring residuals and corrections between levels, and V-, W- or full-multigrid (FMG) cycling. Solutions are written to netCDF. You describe a particular PDE by subclassing mgrid::LinearMultigrid (or mgrid::NonlinearMultigrid) and supplying a discrete operator and a point smoother; the cycle machinery, grid hierarchy, boundary handling and output are provided.

Originally written in 2010-2011 against Blitz++/Boost/netCDF; this version is a C++23 modernisation with a hand-written array layer and a containerised build. See the design notes below.

Build & test

Prerequisites: Docker or Podman, plus the dvt CLI that drives the devcontainer.

dvt up                       # first time: builds the devcontainer
dvt run just build           # configure + compile
dvt run just test            # build, then run the full ctest suite
dvt run just lint            # clang-format + clang-tidy quality gate

dvt run <cmd> executes <cmd> inside the running devcontainer through the workspace user's login shell, so the Homebrew Clang toolchain, Catch2 and the just recipes are on PATH without entering the container; dvt ssh opens an interactive shell there if you'd rather run the recipes directly.

just build configures a Ninja RelWithDebInfo tree under build/ and compiles the library, examples and tests; just test builds then runs ctest; just clean empties the build tree.

just lint is the static-analysis gate: just format-check fails if any tracked C++ source is not clang-format-clean (the just format fix-up set), and just tidy runs clang-tidy (--warnings-as-errors='*', config in .clang-tidy) over src/ and include/multigrid/.

.github/workflows/ci.yml runs just test and just lint in the devcontainer on every push and pull request.

The C++ toolchain (Homebrew Clang 23, C++23) is provided by the cpp-devtools devcontainer feature. The library dependencies -- netcdf-cxx (the netcdf-cxx4 C++ interface) and catch2 -- are installed from the Brewfile by the container's post-create step. No brew install / apt-get install of Boost, Blitz++, cmake or netCDF is needed.

The documentation toolchain is separate and runs through pixi:

just docs                    # pixi run -e docs docs  (mkdocs build --strict)

Visualising the examples

The worked examples' netCDF output is turned into figures by the mgviz Python package (viz/), which runs in its own pixi environment:

dvt run just viz-data        # build + run poisson/mosolov/vorticity into viz/_data/
dvt run just viz             # render docs/src/assets/figures/**  (pixi -e viz)
dvt run just viz-notebooks   # open the marimo notebooks in viz/notebooks/

dvt run targets the one running devcontainer, so all three recipes execute there: just viz-data runs the C++ examples and needs the compiler toolchain, while just viz and just viz-notebooks only need Python and run through pixi run -e viz (the viz environment has no compiler dependencies, so it resolves in the default container fine). .devcontainer/viz is an optional leaner, Python-only container for editor users -- open it with your editor's "Reopen in Container" picker; there is no dvt flag to select it. Data is shared through the workspace bind mount -- viz/_data/ is git-ignored; the committed PNGs under docs/src/assets/figures/ are the source of truth.

Using the library -- subclassing mgrid::LinearMultigrid

To solve a new PDE, subclass mgrid::LinearMultigrid (linear, constant-coefficient operators) or mgrid::NonlinearMultigrid, and override three virtuals:

Override	Signature	Purpose
differential_operator	double differential_operator(mgrid::Level level, int i, int j)	The discrete operator L(u) evaluated at grid point (i, j) on level.
relaxation_updater	void relaxation_updater(mgrid::Level level, int i, int j)	One in-place smoother step: writes the relaxed value straight into solution[level](i, j).
filename	std::string filename(std::string root = "")	The output file stem (without .nc); root is an optional prefix/path.

Inside these, solution[level] and source[level] are mgrid::FDArrays. They provide indexed access (i, j), finite-difference stencils that fall back to one-sided forms near boundaries (dxx, dzz, dxz, dx, dz, ...), grid spacings spacing(0) / spacing(1), and rows() / cols(). finestLevel and coarsestLevel name the ends of the hierarchy.

In the constructor you set boundary conditions on solution and the source term:

#include "multigrid/multigrid.hpp"

class Poisson : public mgrid::LinearMultigrid {
public:
    explicit Poisson(const mgrid::Settings& settings) : mgrid::LinearMultigrid(settings) {
        // Zero-Neumann on left/top, zero-Dirichlet on right/bottom.
        solution.boundaryConditions.set(mgrid::leftBoundary,   mgrid::zeroNeumannCondition);
        solution.boundaryConditions.set(mgrid::rightBoundary,  mgrid::zeroDirichletCondition);
        solution.boundaryConditions.set(mgrid::topBoundary,    mgrid::zeroNeumannCondition);
        solution.boundaryConditions.set(mgrid::bottomBoundary, mgrid::zeroDirichletCondition);
        solution.propagate_boundary_conditions();   // down-sample BCs onto every level

        source_term() = -1.0;                       // constant RHS: nabla^2 u = -1
        mark_source_set();                          // (or fill source_term() element-wise first)
    }

    void solve() override { multigrid(); }

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

    void relaxation_updater(mgrid::Level level, int i, int j) override {
        const double xx = 1.0 / (solution[level].spacing(0) * solution[level].spacing(0));
        const double zz = 1.0 / (solution[level].spacing(1) * solution[level].spacing(1));
        solution[level](i, j) =
            ((solution[level](i + 1, j) + solution[level](i - 1, j)) * xx
             + (solution[level](i, j + 1) + solution[level](i, j - 1)) * zz
             - source[level](i, j)) / (2 * (xx + zz));
    }

    std::string filename(std::string root = "") override {
        return std::format("{}A{:.1f}", root, aspect);
    }
};

Driving it:

mgrid::Settings settings;
settings.aspectRatio = 2.0;
Poisson problem(settings);
problem.solve();                       // or problem.multigrid();
mgrid::FDArray& u = problem.get_result();   // finest-grid solution
problem.write(1, "poisson_");          // writes poisson_A2.0.nc

write's first argument selects how much to emit: 1 = solution only, 2 = solution + gradient magnitude, 3 = also the base-10 log of the residual.

examples/poisson/ is the model for a linear problem and examples/mosolov/ for a nonlinear one (viscoplastic channel flow, solved as a sequence of linear problems with a source term recomputed each step). The tutorials build compiled, tested walkthroughs from these.

Design notes

• Array layer. Fields are a hand-written mgrid::Field2D -- row-major storage over a std::vector<double> -- and mgrid::FDArray composes one Field2D with grid geometry and the finite-difference stencils. mgrid::FDVecArray is a struct of two FDArray components. No Blitz++, no Boost.
• Boundary conditions are copied onto every coarser level (stride-2 down-sample) once, at Stack construction. Mutating them afterwards has no effect on the coarse levels until you call solution.propagate_boundary_conditions() again. update_boundaries applies each condition per node; domain corners are owned by the left/right edges. The 4th-order one-sided Neumann stencil degrades to 3rd-, 2nd- or 1st-order on grids too small to fit it (the coarse multigrid levels).
• Console output uses std::format with std::cout.
• netCDF output goes through mgrid::NcWriter, a PIMPL class whose header (multigrid/netcdf_writer.hpp) pulls in no netCDF headers -- <netcdf> stays confined to src/. The Mosolov example writes a bespoke file: the velocity and strain_rate variables and the bingham_number / aspect_ratio / total_flux global attributes, the dimensionless quantities that define the viscoplastic channel-flow regime.

Documentation

Full API reference, multigrid theory notes and tutorials: https://jesserobertson.github.io/multigrid/ (built with just docs; deployed to GitHub Pages by .github/workflows/docs.yml on every push to master).

Licence

I'm releasing this code for academic use under the Community Research and Academic Programming License. You can read the terms of this license here. Specifically, this licence is designed to achieve the following:

Most open source licenses (1) require source and modifications to be shared with binaries, and (2) absolve authors of legal liability.

An open source license for academics has additional needs: (1) it should require that source and modifications used to validate scientific claims be released with those claims; and (2) more importantly, it should absolve authors of shame, embarrassment and ridicule for ugly code.

The Author reserves all rights to the Program, except for any rights granted under any additional licenses attached to the Program.

Basically I'm pretty happy to let you use this code for non-commercial/academic/research use, provided you cite my work when/if you publish your work. However, if this is too restrictive for you then drop me a line and we can have a chat.