Module methods 

Numerical methods

Everything is written from the published literature; each item below cites its source and links to its implementation. The methods are chosen by name at run time from the [numerics] table of a TOML problem file, but are generic types in the code, so that the inner loops are compiled for the exact scheme in use.

The grid

The domain is a box $\prod_d[a_d,b_d]$ cut into $N_d$ equal cells along each axis, of width $\Delta x_d=(b_d-a_d)/N_d$, in any number of dimensions ([merlin_core::grid::Grid]). Cell $i$ spans $[x_{i-1/2},x_{i+1/2}]$ with centre $x_i=a+(i+\tfrac12)\Delta x$; the unknown stored in it is the cell average of each conserved variable, not a point value. On a cylindrical or spherical grid the first axis is the radius and the cell volume and face areas depend on $r$ (below).

The grid is cut into blocks, and every block stores its cells in a Field padded with ghost cells: extra layers on each side that hold the values a stencil needs from outside the block. A scheme needs as many layers as its stencil reaches (1 for first order and Lax-Wendroff, 2 for MUSCL, 3 for Weno5 and 6 for the eleventh-order rule). Before every operator evaluation the ghost cells are filled, from the neighbouring block’s edge cells where there is one and by the boundary condition at the edge of the domain (Boundary conditions). Each block is then independent of the others, which is what lets blocks run on separate threads.

Two blocks of four cells with one ghost layer. Dashed cells are ghosts; the arrows are the exchange.

Finite volumes and the method of lines

The domain is a uniform grid of cells. Each cell stores the cell average $u_i$, and the integral form of the conservation law becomes one ordinary differential equation per cell (LeVeque 2002, Chapter 4):

\frac{du_i}{dt}=L(u)_i=-\sum_{d}\frac{F_{d,i+1/2}-F_{d,i-1/2}}{\Delta x_d}.

The numerical flux $F_{i+1/2}$ at each face is computed from the states on its two sides: the scheme reconstructs a left and a right value at the face from the cell averages, and a Riemann flux turns the pair into one flux, the same flux for the cell on each side of the face.

The two reconstructed values at a face differ; the Riemann flux resolves the jump.

Because the flux leaving one cell is the flux entering its neighbour, the sum of $u_i\Delta x$ changes only through the domain boundary: conservation to round-off holds by construction, and any time integrator that forms linear combinations of $u+\Delta t\,L(u)$ keeps it. This is the property the conservation tests rely on. The operator is [merlin_finite_volume::solver::operator].

One step

A time step of size $\Delta t$ is, for each Runge-Kutta stage:

1. fill the ghost cells of every block (neighbour copy or boundary condition);
2. for each axis, move it to first place and, in every block, reconstruct the face states, evaluate the Riemann flux at each face and form the flux difference, adding the geometric terms on a curvilinear grid;
3. sum the axes into $L(u)$ (unsplit) or advance one axis at a time (Strang);
4. combine the stage with the old solution by the weights of the integrator.

The time step is set before the step from the largest wave speed, and the last step before an output time is shortened to land on it.

Reconstruction

A flux scheme turns the cell averages around a face into the face flux ([merlin_finite_volume::schemes]).

Scheme	Face states	Order	Reference
FirstOrder	the two adjacent cell values (Godunov)	1	Godunov 1959
LaxWendroff	space-time, linear advection only	2	Lax and Wendroff 1960; LeVeque 2002, §6.1
Muscl	piecewise-linear with a limited slope	2	van Leer 1979; Toro 2009, Ch. 13
Reconstruction	a wide stencil rule: Weno5, or the linear Upwind5 and Upwind11	5 and 11 (1D)	Jiang and Shu 1996; Shu 1998

MUSCL. Each cell is reconstructed as a line $u(x)=u_i+\sigma_i(x-x_i)/\Delta x$, giving $u^L_{i+1/2}=u_i+\sigma_i/2$ and $u^R_{i+1/2}=u_{i+1}-\sigma_{i+1}/2$, which are handed to a Riemann flux. For systems the reconstruction is done in primitive variables (density, velocity, pressure), each limited on its own, and a reconstructed state that is not admissible (negative density or pressure) is replaced by the first-order cell value (Toro 2009, §14.3.1).

Stencil rules. Reconstruction turns the $2r+1$ cells around a face into its two face values with one rule, in the conserved, primitive or characteristic variables (variables in [numerics]). The linear rules need an SSP-RK3 or better integrator and a Courant number below a limit found by von Neumann analysis (1.43 for fifth order, 1.04 for eleventh). In more than one dimension the unsplit operator is second order whatever the rule. The derivations, limits and comparisons are on the high-order page.

The five cells and three candidate stencils of Weno5 for the left value at face $i+\tfrac12$ (Jiang and Shu 1996). The nonlinear weights equal the ideal ones, $\tfrac1{10},\tfrac6{10},\tfrac3{10}$, on smooth data and drop a candidate that crosses a discontinuity.

Slope limiters choose $\sigma$ from the one-sided differences $a=u_i-u_{i-1}$ and $b=u_{i+1}-u_i$. All three return zero at a local extremum ($ab\le0$) and otherwise satisfy $\min(|a|,|b|)\le|\sigma|\le2\min(|a|,|b|)$, which is the condition for a total-variation-diminishing scheme (Sweby 1984):

Limiter	$\sigma$ when $ab>0$	Reference
minmod	$\operatorname{sign}(a)\min(\lvert a\rvert,\lvert b\rvert)$	Roe 1986
van_leer	$2ab/(a+b)$	van Leer 1974
mc (monotonized central)	$\operatorname{sign}(a)\min\bigl(2\lvert a\rvert,\,2\lvert b\rvert,\,\tfrac12\lvert a+b\rvert\bigr)$	van Leer 1977

Slope $\sigma$ over the left difference $a$, against the ratio $r=b/a$, for the three limiters. All lie inside the shaded region and return zero for $r\le0$. minmod is the lower edge for $r>1$, mc the upper edge for $r\ge3$.

The limiter that avoids new extrema is also the one that flattens smooth extrema, so a limited scheme is second order in $L_1$ but converges more slowly in $L_2$ and $L_\infty$. The convergence tables show this, and it is expected, not a defect.

Riemann fluxes

A Riemann flux gives the face flux from a left and a right state ([merlin_finite_volume::riemann]).

Flux	Equation	Idea	Reference
upwind	linear advection	the exact Riemann solution	LeVeque 2002, Ch. 4
godunov	scalar	flux of the exact Riemann solution at the face	Godunov 1959
engquist_osher	scalar	split the flux at its minimum; smoother at sonic points	Engquist and Osher 1980
local_lax_friedrichs	any	central average minus diffusion with the larger spectral radius $\alpha$: $F=\tfrac12(f_L+f_R)-\tfrac12\alpha(u_R-u_L)$	Rusanov 1961
roe_no_fix	scalar	upwinding by the Roe speed with no entropy fix; kept to demonstrate the failure	Roe 1986; Harten and Hyman 1983
exact	Euler	flux of the exact Riemann solution at the face	Godunov 1959; Toro 2009, Ch. 4
hll	any (a gas in the problem file)	two waves $S_L,S_R$, one averaged state between them	Harten, Lax and van Leer 1983
hllc	systems with a middle wave (Euler)	HLL with the contact wave restored	Toro, Spruce and Speares 1994

HLL and HLLC take their wave speeds from the system, not from a formula of the flux (Equation::wave_speed_estimates). The default is the signed estimate of Davis, the smallest and largest speed of the two states (Toro 2009, §10.5.1), so HLL works for any equation, a scalar law included; the Euler equations override it with the pressure-based estimates of Toro 2009, §10.5.2, in which the factor by which a shock outruns the sound speed comes from the fluid model. The problem file offers hll for a gas only. HLL is

F=\frac{S_R F_L-S_L F_R+S_LS_R(U_R-U_L)}{S_R-S_L}

between the two speeds. It is robust and satisfies the entropy condition, but with no contact wave a density jump at constant pressure and velocity is smeared. HLLC adds a middle wave at speed

S_*=\frac{p_R-p_L+\rho_Lu_L(S_L-u_L)-\rho_Ru_R(S_R-u_R)}{\rho_L(S_L-u_L)-\rho_R(S_R-u_R)}

and keeps an isolated contact exactly; the validation shows both behaviours. The middle wave and the states beside it are also the system’s: a system implements StarStates to be used with HLLC, and one that does not, a scalar law, is refused at compile time.

The approximate fans in the $(x,t)$ plane, to compare with the exact fan. HLL replaces everything between the outer speeds $S_L$ and $S_R$ by one state; HLLC splits it in two at $S_*$, which is where the contact is. A choice of flux is also a choice about which failures to accept, which is why the limitations page lists the unfixed Roe flux and HLL’s smearing as demonstrations.

Time integration

All integrators are strong-stability-preserving Runge-Kutta methods in Shu-Osher form (Shu and Osher 1988; Gottlieb, Shu and Tadmor 2001): $u^{(k)}=a_ku^n+b_k\bigl(u^{(k-1)}+\Delta t\,L(u^{(k-1)})\bigr)$ with non-negative weights and $a_k+b_k=1$. Each stage is a convex combination of forward-Euler steps, so any property of forward Euler under a time-step limit (no new extrema, positivity, TVD) carries over to the whole step.

Name	Stages	Order	Reference
forward_euler	1	1
ssp_rk2	2	2	Shu and Osher 1988 (Heun)
ssp_rk3	3	3	Shu and Osher 1988, Eq. 2.18

They are one table of weights over a single stage loop ([merlin_core::time::integrator::TimeIntegrator]). The time step is $\Delta t=\mathrm{CFL}\,\Delta x/\max|\lambda|$, recomputed every step from the solution, with $\mathrm{CFL}\le1$ for the first-order and Lax-Wendroff schemes and $\le\tfrac12$ for MUSCL. Output times are met exactly by shortening the last step before them, and the schedule tests confirm that a shortened step does not break conservation.

Level set

merlin_levelset advects $\phi$ in $\phi_t+u\cdot\nabla\phi=0$ with the same stencil rules applied as one-sided derivatives, chosen by the sign of $u$, on point values at the cell centres. It runs on the same blocks, executor and time integrators (third order or better). Two reinitialisations rebuild a signed distance function: a pseudo-time equation (Sussman, Smereka and Osher 1994) and fast sweeping (Zhao 2005); both are on the high-order page.

The derivative $\phi_x$ at cell $i$ for Weno5: six points, biased towards the side the information comes from, which is chosen by the sign of $u$. Upwind5 uses the same points, Upwind11 twelve.

Dimensions

The schemes and Riemann fluxes are written for one direction. The operator applies them along each axis by moving that axis into first place (Equation::along, swap_axes), so for a gas the momentum along the axis plays the part of the x-momentum, and moves the fluxes back. The 1D results are therefore unchanged to the last bit, and the tests compare a 1D problem embedded in 2D along either axis with the 1D run bit for bit. Two ways of combining the axes are available ([merlin_finite_volume::solver::DimensionalMethod]):

• Unsplit: $L=L_x+L_y+L_z$ advanced together, with $\sum_d|a_d|\Delta t/\Delta x_d\le\mathrm{CFL}_{\max}$.
• Strang split: one axis at a time in the order $X(\Delta t/2)\,Y(\Delta t)\,X(\Delta t/2)$ (Strang 1968; LeVeque 2002, §17.1), which needs only the per-axis limit but adds a splitting error.

Cylindrical and spherical coordinates

A radial cell of each geometry, per unit angle (cylinder, also per unit height) or unit solid angle (sphere). The flux difference is weighted by the face areas $A_\pm$ and divided by the volume $V$; subtracting $P$ first is what makes the update vanish exactly for a gas at rest.

In the radial direction the flux difference is weighted by the face areas $A$ and divided by the cell volume $V$ of the geometry ([merlin_core::geometry::Geometry]), per unit angle:

L_i=-\frac{A_+\,(F_+-P_i)-A_-\,(F_--P_i)}{V_i}.

$P$ is the pressure part of the flux, the vector $(0,p,0)$ for a gas (Equation::pressure_vector). Subtracting it before weighting and adding the cell value back is the same as the usual source term $\alpha p/r$, written so that it cancels exactly for a gas at rest: for a uniform $p$ every $F_\pm-P_i$ vanishes in the momentum component and the update is exactly zero, not merely small. The sweep along $z$ is the Cartesian one, since face area and volume scale together. The axis $r=0$ has zero face area and a reflecting boundary, and integrals use the physical cell volume, including the factor $2\pi$ or $4\pi$.

Because the cell next to the axis has $A_+/V=2/\Delta r$ (cylinder) or $3/\Delta r$ (sphere), against $1/\Delta x$ for a flat cell, the time step multiplies the radial speed by that geometric factor, keeping the Courant number of every cell at most the nominal one.

Swirl. For a gas with an azimuthal velocity $w$ the pressure vector becomes $P=(0,\,p+\rho w^2,\,0,\ldots)$: the centrifugal force $\rho w^2/r$ is the same kind of term as $p/r$. The momentum $\rho w$ is not conserved, but $r\rho w$ is, and its flux through a face of radius $r$ is $r$ times the flux of $\rho w$; so the radial sweep updates that component with the face weight $r^2$ and the volume $\int r^2\,dr$ (those of a spherical shell) in place of $r$ and $r\,dr$, and no source term: $\Delta(\rho w)_i=-\bigl(a_+F_+-a_-F_-\bigr)/W_i$ with $a=r^2$ and $W_i=(r_+^3-r_-^3)/3$. The sum of $W_i(\rho w)_i$ is then conserved to round-off, which the tests confirm as a drift of exactly 0. The sweep along $z$ is unchanged, because the weight $W_i/V_i$ of a cell cancels. The symmetry axis reverses $\rho u$ and $\rho w$ (Reflecting::also_flipping), and an outer wall, which is free-slip, reverses only $\rho u$. The equation tells the operator which component is the angular momentum (Equation::angular_momentum_component); for a gas it is the third momentum, so the swirl is the gas with three velocity components on a two-dimensional cylindrical grid.

Boundary conditions

Ghost cells are filled before every operator evaluation ([merlin_core::boundary]): periodic, outflow (zero gradient), reflecting wall (the normal momentum changes sign), and a fixed value or state. Each axis has its own pair, so a problem can be periodic in one direction and walled in another.

The two ghost layers next to the left edge of a domain, for each condition. R is the reflection: it flips the sign of the momentum normal to the wall and leaves the other variables alone.

Blocks and threads

The grid is cut into a regular layout of blocks, each owning its own field with ghost cells. Before a stencil is applied, each block’s ghost cells are filled from its neighbours’ edge cells, one axis at a time over the full padded extent, so the corner and edge ghosts need no special case ([merlin_core::blocks]).

A kernel reads one block and its ghosts and writes only that block, so the blocks can be processed by different threads without locks. They run on a rayon pool wrapped in [merlin_core::parallel::Executor], the only module that names rayon. Block results are combined in block order (or by taking a maximum), never in the order in which threads finish, so the answer is identical to the last bit for any number of threads and any block layout. This is a design property that the tests compare exactly, not an approximate agreement; see Verification method.

Floating point

f64 throughout, no unsafe, and no allocation inside the time loop (scratch buffers are made once, outside it).

References

• Godunov, S. K. (1959), A difference scheme for numerical solution of discontinuous solution of hydrodynamic equations, Mat. Sbornik 47.
• Lax, P. D. and Wendroff, B. (1960), Systems of conservation laws, Comm. Pure Appl. Math. 13.
• Rusanov, V. V. (1961), The calculation of the interaction of non-stationary shock waves with barriers, Zh. Vychisl. Mat. Mat. Fiz. 1.
• Strang, G. (1968), On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5.
• van Leer, B. (1974, 1977, 1979), Towards the ultimate conservative difference scheme II, III, V, J. Comput. Phys. 14, 23, 32.
• Engquist, B. and Osher, S. (1980), Stable and entropy satisfying approximations for transonic flow calculations, Math. Comp. 34.
• Harten, A., Lax, P. D. and van Leer, B. (1983), On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Review 25.
• Harten, A. and Hyman, J. M. (1983), Self adjusting grid methods for one-dimensional hyperbolic conservation laws, J. Comput. Phys. 50.
• Sweby, P. K. (1984), High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal. 21.
• Sussman, M., Smereka, P. and Osher, S. (1994), A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114.
• Roe, P. L. (1986), Characteristic-based schemes for the Euler equations, Annu. Rev. Fluid Mech. 18.
• Shu, C.-W. and Osher, S. (1988), Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77.
• Toro, E. F., Spruce, M. and Speares, W. (1994), Restoration of the contact surface in the HLL Riemann solver, Shock Waves 4.
• Jiang, G.-S. and Shu, C.-W. (1996), Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126.
• Gottlieb, S., Shu, C.-W. and Tadmor, E. (2001), Strong stability-preserving high-order time discretization methods, SIAM Review 43.
• Zhao, H. (2005), A fast sweeping method for Eikonal equations, Math. Comp. 74.
• LeVeque, R. J. (2002), Finite Volume Methods for Hyperbolic Problems, Cambridge.
• Toro, E. F. (2009), Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed., Springer.