Skip to content

Periodic multigrid and barotropic vorticity

The core solver was written for a bounded rectangle with a boundary condition on every edge (see discretisation and boundary conditions). Making an axis periodic instead -- so that the last grid line is a neighbour of the first -- turns the same multigrid machinery into a spectral-quality Poisson solver on a torus, which is what the barotropic-vorticity example is built on.

Periodic multigrid

The Axis helper

Every place the discretisation used to branch on "is this the first/last row?" now goes through a small value type, one per direction, carried by FDBase:

Member Bounded axis Periodic axis
n, h point count and spacing point count and spacing
periodic false true
lo(), hi() 1, n - 1 0, n
edge_lo(), edge_hi() true at an actual edge always false
shift(i, d) i + d (i + d + n) % n -- wraps

hi() is one past the last swept index: a sweep runs for (i = lo(); i < hi(); ++i), so a bounded axis visits the interior 1 .. n - 2 and a periodic axis visits every line 0 .. n - 1.

The stencil methods, the red--black sweep and the grid-transfer operators all read neighbours through shift and ask edge_lo/hi whether a one-sided form is needed. On a periodic axis edge_lo/hi is never true and shift wraps, so the scattered if (periodic) special cases collapse into one code path.

Wrap-around stencils

An interior central difference such as \((\,f_{i+1,j} - f_{i-1,j})/2h_x\) is used at every line of a periodic axis; the indices i-1 and i+1 are taken mod n, so the point at i = 0 reads i = n-1 on one side and i = 1 on the other. No one-sided boundary stencil is ever assembled for that axis, and the discrete Laplacian is the plain 5-point operator everywhere.

Periodic restriction and interpolation

Full-weighting restriction is the 9-point average

\[ r_{2h}(I,J) = \frac{1}{16}\!\!\sum_{|p|,|q|\le 1} w_{p}\,w_{q}\, r_h(2I+p,\,2J+q), \qquad w_0 = 2,\ w_{\pm 1} = 1, \]

and bilinear interpolation is its transpose. On a bounded grid the stencil is truncated at an edge (the missing neighbours are dropped and the weights no longer sum the same way). On a periodic axis nothing is dropped: the missing neighbour is simply the one across the wrap, shift supplies it, and the divisor stays \(16\). The transfers are therefore exact adjoints on the torus, which is what keeps the coarse-grid correction from injecting a spurious mean.

The hierarchy

Stack still builds \(n_x = n_z\) at every level. A bounded axis coarsens \(n \to 2(n-1)+1\) (endpoints shared between levels); a periodic axis has no shared endpoint, so it refines by straight doubling, \(n \to 2n\), and the finest level is minimumResolution * 2^(numberOfGrids - 1) points per axis. The coarsest resolution must be even: the red--black sweep colours a point by the parity of \(i + j\), and only an even n makes that parity continuous across the wrap (\(i = n\) and \(i = 0\) are the same line and must share a colour).

The constant null space

With a periodic (or all-Neumann) axis the discrete Laplacian annihilates the constant field: \(\mathcal{L}_h\,\mathbf{1} = 0\). The problem \(\mathcal{L}_h\psi = \zeta\) is then solvable only when \(\langle \zeta \rangle = 0\), and its solution is fixed only up to an additive constant. The solver handles both ends of this:

  • the right-hand side is made mean-zero once, with project_out_constant, before the cycle starts (PeriodicPoisson::invert);
  • project_out_constant is applied to the solution after every relaxation sweep (MultigridBase::relax, guarded by fully_periodic()), which stops the unconstrained constant mode from random-walking away during the iteration and pins \(\psi\) mean-zero on return.

Because one full-multigrid pass only reaches the discretisation-error level, PeriodicPoisson::invert wraps it in an iterative-refinement loop, driving the finest-grid residual down to \({\sim}10^{-11}\) so the discrete energy identity the time integrator relies on holds to round-off.

Barotropic vorticity

The example integrates the barotropic vorticity equation on a beta-plane over a doubly-periodic square:

\[ \frac{\partial \zeta}{\partial t} = -J(\psi, \zeta) - \beta\,\frac{\partial \psi}{\partial x} + D(\zeta), \qquad \nabla^2 \psi = \zeta . \]

Here \(\zeta\) is the (scalar) vertical vorticity, \(\psi\) the streamfunction, and the non-divergent velocity is

\[ u = -\frac{\partial \psi}{\partial z}, \qquad w = \phantom{-}\frac{\partial \psi}{\partial x}. \]

\(\beta = \mathrm{d}f/\mathrm{d}y\) is the planetary-vorticity gradient; the \(-\beta\,\partial_x\psi\) term is what makes Rossby waves propagate. \(D\) is optional scale-selective dissipation,

\[ D(\zeta) = \nu\,(-1)^{p+1}\,\nabla^{2p}\zeta , \]

with \(p\) = dissipationOrder: \(p = 1\) is ordinary Laplacian friction, \(p = 2\) is biharmonic hyperviscosity (the sign makes it a sink at every wavenumber). Each time step recovers \(\psi\) from \(\zeta\) with the periodic multigrid Poisson solve above.

The Arakawa Jacobian

The advection term \(J(\psi,\zeta) = \psi_x\zeta_z - \psi_z\zeta_x\) is discretised with the energy- and enstrophy-conserving 9-point form of Arakawa (1966),

\[ J = \tfrac{1}{3}\bigl(J^{++} + J^{+\times} + J^{\times+}\bigr), \]

the arithmetic mean of the three second-order finite-difference Jacobians: the product of two centred first differences (\(J^{++}\)), and the two forms (\(J^{+\times}\), \(J^{\times+}\)) that difference one factor and average the other across the diagonal.

On the continuous torus, integration by parts gives three identities: \(\langle J(\psi,\zeta)\rangle = 0\), \(\langle \psi\,J(\psi,\zeta)\rangle = 0\) (kinetic energy is advected, not created), and \(\langle \zeta\,J(\psi,\zeta)\rangle = 0\) (so is enstrophy). A generic discrete Jacobian preserves at most one of the last two. Arakawa's key observation is that each of \(J^{++}\), \(J^{+\times}\), \(J^{\times+}\) fails the discrete summation-by-parts identities in a different way, and the symmetric average of all three satisfies both

\[ \sum_{i,j} \psi_{i,j}\,J(\psi,\zeta)_{i,j} = 0, \qquad \sum_{i,j} \zeta_{i,j}\,J(\psi,\zeta)_{i,j} = 0 \]

exactly on a periodic grid -- the discrete operator is antisymmetric enough that the semi-discrete system conserves domain-mean energy \(-\tfrac12\langle \psi\,\zeta\rangle\) and enstrophy \(\tfrac12\langle \zeta^2\rangle\) to machine precision (up to the time-stepping error and any explicit \(D\)). This suppresses the systematic enstrophy build-up at the grid scale that sinks a naive scheme in 2-D turbulence.

A. Arakawa, Computational design for long-term numerical integration of the equations of fluid motion, J. Comput. Phys. 1 (1966) 119.

Time stepping

Time integration uses the three-stage, third-order strong-stability-preserving Runge--Kutta scheme (SSP-RK3) in Shu--Osher form. Writing \(L(u) = -J(\psi,u) - \beta\psi_x + D(u)\) for the full spatial operator:

\[ \begin{aligned} u^{(1)} &= u^{n} + \Delta t\,L(u^{n}), \\[2pt] u^{(2)} &= \tfrac{3}{4}u^{n} + \tfrac{1}{4}u^{(1)} + \tfrac{1}{4}\Delta t\,L(u^{(1)}), \\[2pt] u^{n+1} &= \tfrac{1}{3}u^{n} + \tfrac{2}{3}u^{(2)} + \tfrac{2}{3}\Delta t\,L(u^{(2)}). \end{aligned} \]

Each stage is a forward Euler step followed by a convex average with the old state, so any stability bound that forward Euler satisfies is preserved by the full step -- useful here because the Arakawa Jacobian is only neutrally stable. The scheme is third-order accurate in \(\Delta t\) and needs three RHS evaluations (three Poisson solves) per step.

The step size is set by an advective CFL condition,

\[ \Delta t = C\,\frac{\min(h_x, h_z)}{\max |\nabla \psi|}, \]

with a Courant factor \(C \in (0,1]\); BarotropicVorticity::cfl_dt evaluates this from the streamfunction left by the last RHS call (falling back to \(C\,\min(h_x,h_z)\) for a fluid at rest).

Decaying turbulence

With \(\beta = 0\) and only weak hyperviscosity, a random initial vorticity field relaxes into freely decaying two-dimensional turbulence. The dynamics are dominated by the forward enstrophy cascade: enstrophy moves to small scales and is removed there by \(D\), while energy stays almost constant (it can only cascade to large scales). As the flow evolves, like-signed vorticity merges into a sparse population of long-lived coherent vortices separated by quiescent filamented background.

Two diagnostics track that transition:

  • the vorticity kurtosis \(\langle\zeta^4\rangle/\langle\zeta^2\rangle^2\) starts near the Gaussian value \(3\) for the random field and rises well above it as the field becomes intermittent (vorticity concentrated in a few strong cores);
  • the energy spectrum \(E(k)\) steepens. The Batchelor/Kraichnan enstrophy- cascade prediction is \(E(k) \sim k^{-3}\); the emergence of coherent vortices steepens the observed mid-range slope past the Batchelor \(-3\). This example settles near \(-3.5\) over the \([8,\,N/4]\) fit band (the steeper \(-4\) to \(-5\) regime lives in the dissipation range above \(k \approx 25\), outside the band); the regression assertion only requires the slope to fall in \([-5,\,-3]\).

J. C. McWilliams, The emergence of isolated coherent vortices in turbulent flow, J. Fluid Mech. 146 (1984) 21.

The example's spectral regression test runs this scenario and asserts that, after the vortices have formed, the vorticity kurtosis exceeds \(5\), the fitted mid-range spectral slope lies in \([-5, -3]\), energy is conserved to a few percent, and enstrophy has measurably decayed.