Interactive explainer

Geometric multigrid, hands-on

Gravity couples every cell of a simulation to every other cell. Multigrid computes the potential in a time proportional to the number of cells, by working on a ladder of coarser and coarser copies of the grid. Three experiments show why that works. Everything below runs live in your browser.

A global problem at every step

The gravitational potential \(\Phi\) of a gas with density \(\rho\) obeys the Poisson equation,

\[ \nabla^2 \Phi = 4\pi G \rho . \]

The fluid equations are local: within one time step a cell only exchanges fluxes with its neighbours. The Poisson equation is not. Move one mass and the potential changes everywhere at once. A self-gravitating simulation has to solve it at every step, in lock-step with the gas, so the cost of that solve decides what is affordable.

Gravitational potential of two orbiting masses in AthenaK, with equipotential lines. The multigrid solver is called at every step; the mesh blocks refine around each mass and follow it.

Relaxation only removes the rough part

On a grid of spacing \(h\) the equation says that every value equals the average of its neighbours, minus a source term. In one dimension, with \(f = 4\pi G\rho\):

\[ \Phi_i = \tfrac{1}{2}\left(\Phi_{i-1} + \Phi_{i+1} - h^2 f_i\right) . \]

The simplest solver applies this rule to every cell and repeats. It is called relaxation, and it does converge, slowly. To see why it is slow, look at the error instead of the solution. Below, the error starts as the sum of a smooth, a medium and a rough wave.

After a handful of sweeps the rough wave is gone and the curve is smooth. The smooth wave has hardly changed, and it will take hundreds of sweeps to remove. A sweep only moves information by one cell, so an error that spans the whole grid needs a number of sweeps that grows with the square of the number of cells. Relaxation is a poor solver and an excellent smoother.

The same behaviour on the potential of a binary, in two dimensions: the solution, the residual, and the error that is left for three relaxation methods. After 400 sweeps between three quarters and nine tenths of it is still there.

A smooth error is a rough error on a coarser grid

Once the error is smooth it can be represented on a grid with half as many cells in each direction and nothing is lost. On that grid the same wave spans half as many cells, so relative to the grid it is rougher, and relaxation is effective again. A sweep there is also cheaper: two times in one dimension, eight times in three.

What is passed down is the equation for the error, which has the same form as the original one with the residual \(r = f - \nabla^2\Phi\) as its source. Two operators connect the grids: restriction averages fine values onto the coarse grid, and prolongation interpolates coarse values back. Apply the idea recursively until the grid is so small that solving on it is trivial, and you have a V-cycle.

Each cycle divides the error by about the same factor, whatever the resolution, and its cost is that of a few sweeps on the finest grid: the coarser grids add \(1/8 + 1/64 + \dots\) in three dimensions, one seventh in total. A fixed number of cycles therefore solves the problem, and the work is proportional to the number of cells.

One V-cycle on the three-dimensional potential of a binary, from 257³ points down to 3³ and back. Left: the solution on the finest grid, the correction on the coarser ones. Middle: the residual. Right: the position in the cycle.

Twelve V-cycles on the same problem. Each one makes the residual about ten times smaller.

Gravity in a box

The same algorithm in two dimensions, solving for the potential of three masses in a box whose walls are held at zero. The grid has 128×128 cells and seven coarser levels beneath it. Drag the masses.

Things to try

From nothing. Choose “Relaxation only” and press “Start from zero”. The wells grow outwards from each mass, a cell per sweep. Now choose “Multigrid” and start from zero again.
The error that relaxation leaves. While relaxation is running, switch the view to “Error”. What remains is large and smooth, which is exactly what a coarser grid represents well.
A moving source. Drag a mass under each solver. With relaxation the potential trails behind the mass; with multigrid it follows.
The chart. It plots the error against the work spent, counted in sweeps of the finest grid. A V-cycle costs about seven of them.

Full multigrid

A V-cycle improves a guess. Full multigrid builds the first guess itself: it solves the problem on the coarsest grid, interpolates that solution to the next grid as a starting point, runs a V-cycle there, and continues up to the finest level. By the time the finest grid is reached the guess is already close, and the whole pass costs little more than a V-cycle or two.

Full multigrid on the same problem, stage by stage: the solution, the residual, and the nested sequence of V-cycles from the coarsest grid to the finest.

What a production solver adds

The multigrid self-gravity solver in AthenaK follows the design of Tomida & Stone (2023) for Athena++. It differs from the toy above in a few ways that matter.

Three dimensions, on GPUs. The smoother is red–black Gauss–Seidel: cells are coloured like a checkerboard, and all cells of one colour can be updated at the same time, which is what a GPU needs.
Isolated boundaries. A star or a binary does not sit in a box with walls at zero potential. The potential at the boundary is obtained from a multipole expansion of the mass inside.
Mesh refinement. With refined patches a level may cover only part of the domain. The solver uses the Full Approximation Scheme, in which the coarse levels carry the full solution instead of a correction, so that levels stay consistent where they overlap.

The same V-cycle in the Full Approximation Scheme: the coarse grids now hold the solution itself.

Refinement must not perturb the orbit

The delicate part of mesh refinement is the boundary between a coarse and a fine block. Values there have to be interpolated, and an interpolation that does not conserve mass acts as a small spurious force. A binary is a sensitive test: over a few orbits that force changes the separation. In the comparison below both runs use two levels of refinement that follow the masses.

Left: mass-conserving prolongation; the masses stay on the analytic orbit (dashed). Right: plain interpolation at the coarse–fine boundaries; the orbit drifts. Thin lines are the mesh blocks.