Numerical methods

FLUME discretises the conservation laws with high-order finite differences on point values: the state ${\mathbf{q}}_i$ is the value at the centre of cell $i$ of a uniform Cartesian block, and the semi-discrete update is the conservative flux difference

$$\frac{d{\mathbf{q}}_i}{dt}=-\sum _d\frac{{\hat{\mathbf{F}}}_{d,i+1/2}-{\hat{\mathbf{F}}}_{d,i-1/2}}{\mathrm{\Delta }{x}_d}+\mathbf{S}({\mathbf{q}}_i),$$

integrated in time by an explicit Runge–Kutta scheme. The numerical flux ${\hat{\mathbf{F}}}_{i+1/2}$ is built so that the difference is a high-order approximation of $\partial \mathbf{F}/\partial x$ (Shu & Osher 1989). Two space operators are available, selected by [numerics] scheme_space:

scheme_space	Face flux	Models
weno	WENO reconstruction of the Lax–Friedrichs-split fluxes, per characteristic field	Euler, MHD
weno-riemann	WENO interpolation of the face states, a Riemann solver and a high-order correction (Chen, Tóth & Gombosi 2016)	Euler (LLF, HLL, HLLC); MHD (LLF, HLL, HLLD)

Both operators store the face flux at the same place, so seams, AMR reflux and the flux difference are shared.

Grid, blocks and ghost cells

The domain is a forest of Cartesian blocks of ni × nj × nk cells, organised by an octree (or quadtree) with Morton ordering. Each block carries ngc ghost cells per side, filled before every residual evaluation from the neighbouring blocks (MPI exchange across ranks), from the physical boundary conditions, or by the coarse–fine interpolation at an AMR seam. A direction flagged null_x/y/z = .true. is inactive: its fluxes are not computed and the solution is constant along it (1-D and 2-D problems are 3-D runs with null directions). The ghost width must cover the stencil: ngc ≥ S for WENO of order $2S-1$ (ngc = 3 for WENO5), and ngc ≥ 2 for weno-riemann.

WENO reconstruction

A WENO scheme of order $2S-1$ ([weno] scheme = weno-u-{3,5,7,9}, $S=2\dots 5$) combines $S$ candidate polynomials $p_k$, each on $S$ points, with nonlinear weights

$${\hat{v}}_{i+1/2}=\sum _{k=0}^{S-1}{\omega }_k{p}_k,{\omega }_k=\frac{{\alpha }_k}{\sum _j{\alpha }_j},{\alpha }_k=\frac{d_k}{(\epsilon +{\beta }_k{)}^m},$$

with the optimal (linear) weights $d_k$, the Jiang–Shu smoothness indicators ${\beta }_k$, $\epsilon ={10}^{-6}$ and the exponent $m=S$ ($S-1$ for $S>4$), the same on the host and on the device. On smooth data ${\omega }_k\to d_k$ and the combination reaches order $2S-1$; the stencils crossing a discontinuity get ${\omega }_k\to 0$. The centred schemes weno-c-* are refused by FLUME.

The absolute $\epsilon$ ties the weights to the units of the data: ${\beta }_k$ scales as the square of the field, so on a field of magnitude ${10}^{-4}$ the ratio ${\beta }_k/\epsilon$ is ${10}^8$ times smaller than on the same field at magnitude 1, the weights collapse to the linear ones and the scheme stops limiting (measured on Sod at an interstellar density, ${10}^{-21}$ kg/m³: the total variation of the density grows by 33%, and by 178% with weno-riemann, whose sensor reads the same weights; NV-6). Two cures: the reference layer, which hands the solver numbers of order 1, or [weno] weights = si, scale-invariant weights after Don, Li, Wang and Wang (2022):

$${\alpha }_k=\frac{d_k}{(\epsilon +{\beta }_k/{\mu }^2{)}^m},$$

$\mu$ a magnitude of the reconstructed field with its units, so that ${\beta }_k/{\mu }^2$ does not depend on them. Which magnitude matters. The mean of $|v|$ over the stencil (the descaler of Don et al.) fails on the characteristic fields of FLUME: projected on the eigenvectors of the face state, a field such as the shear wave $\rho (v-v_{\text{Roe}})$, or one carrying only a small wave, is near zero by construction, its $\mu$ is its own variation, $\epsilon {\mu }^2\to 0$ and the weights turn fully nonlinear on smooth data (issue #49: the isentropic vortex fell from order 5.79 to 4.10). FLUME takes instead the magnitude of the state that the field projects,

$${\mu }_k=⟨\sum _v|l_{kv}|\frac{1}{2}(|f_v|+{\alpha }_{\text{max}}{\hat{q}}_v){⟩}_{\text{stencil}}$$

for the split fluxes ($⟨\sum _v|l_{kv}|{\hat{q}}_v⟩$ for the interpolated fields of weno-riemann), $l_k$ the left eigenvector (the identity for conservative variables), ${\alpha }_{\text{max}}$ the largest speed of the stencil (MHD: $max(|u_n|+c_f)$, never the cleaning speed $c_h$) and ${\hat{q}}_v\ge |q_v|$ a reference magnitude of each component: $|\rho |$, $|E|$, $\parallel \rho \mathbf{u}{\parallel }_1+\sqrt{\rho p}$ for every momentum component, $\parallel \mathbf{B}{\parallel }_1+\sqrt{p}$ for every field component, $|\psi |+c_h{\hat{q}}_B$ (GLM) or $|\psi |+{\hat{q}}_B$ (EGLM). The reference magnitudes matter: with $|q_v|$, a field built on components that vanish (the Alfvén rows of a 2-D problem touch only $\rho w$ and $B_z$) or are exponentially small (the field in the far field of the magnetised vortex, $\sim {10}^{-21}$) gets a descaler of its own size, the weights treat round-off as data, and the vortex went unstable at 256² (issue #49, MV-7). ${\mu }_k$ bounds $|f_k^{\pm }|$ cell by cell, never vanishes on a physical state, and at unit scale acts like the absolute $\epsilon$ (on the isentropic vortex, $L_1=7.22\cdot {10}^{-8}$ at 256², order 4.97, against $8.38\cdot {10}^{-8}$ with js). Every term scales exactly like the field ($\sqrt{\rho p}$ and $\sqrt{p}$ by exact powers of two), and a power of two is exact in floating point, so multiplying the data by $2^n$ multiplies the reconstruction by $2^n$ bit for bit (unit tests test_flume_weno_weights on the host, test_flume_weno_interpolation_fnl on the device; the whole solver: the scaling oracle). $\mu$ is floored at the smallest normal number. Primitive variables have no projection: a velocity or a field component crossing zero has no magnitude of its own, so si with weno-riemann on primitive variables is refused. One property of js is lost at round-off: the smoothness indicator of a uniform field is a quadratic form, round-off rather than zero, and the si descaler varies from face to face, so a uniform field in a non-uniform state (the $B_n$ of a 1-D MHD problem) is reconstructed with face values that differ in the last bits (RJ2a: GLM equals the run without cleaning to ${10}^{-13}$, not bitwise).

The default js is the first formula, bitwise unchanged: the kernel receives $1/\mu$, which is exactly 1 for js, so the choice costs no branch in the kernels.

The WENO kernel takes its coefficient tables as arguments, so one kernel serves two purposes:

• reconstruction (weno): from point values of a flux, the face value whose difference approximates the derivative;
• interpolation (weno-riemann): from point values of the state, the state at the face. The interpolation tables (for WENO5: candidates exact to degree 2, combination $(3,-20,90,60,-5)/128$, linear weights $1/16,10/16,5/16$) are generated with exact rational arithmetic by riemann-flux/weno_interpolation_tables.py.

scheme_space = weno: flux splitting

At each face $i+1/2$ of direction $d$, over the stencil of cells $i-S+1\dots i+S$:

1. the eigenvectors are evaluated at the face (Euler: the Roe average of cells $i$, $i+1$; MHD: their arithmetic average);
2. the state and the physical flux of every stencil cell are projected on the left eigenvectors, $w_k={\mathbf{l}}_k\cdot \mathbf{q}$, $g_k={\mathbf{l}}_k\cdot \mathbf{F}$ (reconstruction_variables = characteristic), or kept component by component (conservative);
3. each field is split with a local Lax–Friedrichs speed, $g_k^{\pm }=\frac{1}{2}(g_k\pm {\alpha }_k{w}_k)$, where ${\alpha }_k=\underset{\text{stencil}}{max}|{\lambda }_k|$ per characteristic field (conservative: one speed $max(|u_n|+c)$, and with GLM at least $c_h$);
4. $g^{+}$ is reconstructed from the left-biased stencil and $g^{-}$ from the right-biased one, and the face flux is projected back, $\hat{\mathbf{F}}=\sum _k{\mathbf{r}}_k({\hat{g}}_k^{+}+{\hat{g}}_k^{-})$.

For MHD the characteristic decomposition is block-diagonal: the 7×7 wave core plus the $(B_n,\psi )$ pair with speeds $\mp c_h$ (without cleaning the $B_n$ row has speed 0 and its face flux is exactly zero). The core is the Roe–Balsara set in the normalisation of Stone et al. (2008, appendix B), with ${\alpha }_f^2=(a^2-c_s^2)/(c_f^2-c_s^2)$ and ${\alpha }_s^2=(c_f^2-a^2)/(c_f^2-c_s^2)$. At high β with $\mathbf{B}$ near the face normal one of the two numerators is a difference of $O(a^2)$ speeds far smaller than them, so FLUME evaluates the smaller one without subtracting, $2a^2{c}_t^2/((c_f^2-c_s^2)\pm (v_{A,n}^2+c_t^2-a^2))$ with $c_t^2=B_t^2/\rho$, and the other as its complement. Subtracted directly, the normalisation carried an absolute error of $\sim {\epsilon }^{1/2}$ that the characteristic WENO turned into a ${10}^{-6}$ difference of $\mathbf{B}$ between compilers (#79).

scheme_space = weno-riemann: interpolation, Riemann flux and correction

The face flux of the hybrid scheme (Chen, Tóth & Gombosi 2016) is

$${\hat{\mathbf{F}}}_{i+1/2}={\mathbf{F}}^{RS}+s[c_1{\mathbf{F}}^{RS}+c_2({\mathbf{f}}_i+{\mathbf{f}}_{i+1})+c_3({\mathbf{f}}_{i-1}+{\mathbf{f}}_{i+2})-{\mathbf{F}}^{RS}],$$

with ${\mathbf{F}}^{RS}$ the Riemann flux of the two interpolated face states, ${\mathbf{f}}_m$ the physical flux of cell $m$ and $s\in \{0,1\}$ a per-face sensor:

1. Face states. The fields of the stencil are interpolated to the face from the left and from the right (reconstruction_variables: characteristic, the recommended default for Euler and MHD, with the eigenvectors of the face; or primitive, $(\rho ,u,v,w,p)$, plus $\mathbf{B}$ and $\psi$ for MHD). A face state with non-positive density or pressure is replaced by the adjacent cell's state.

2. Riemann flux (riemann_solver):

   • llf: ${\mathbf{F}}^{RS}=\frac{1}{2}({\mathbf{f}}_L+{\mathbf{f}}_R)-\frac{1}{2}\alpha ({\mathbf{q}}_R-{\mathbf{q}}_L)$, $\alpha =max(|u_n|+a)$ of the two states;
   • hll (Harten, Lax & van Leer 1983) with the Einfeldt speeds $S_L=min(u_{n,L}-a_L,{\tilde{u}}_n-\tilde{a})$, $S_R=max(u_{n,R}+a_R,{\tilde{u}}_n+\tilde{a})$ (Roe averages): ${\mathbf{F}}^{RS}=(S_R{\mathbf{f}}_L-S_L{\mathbf{f}}_R+S_L{S}_R({\mathbf{q}}_R-{\mathbf{q}}_L))/(S_R-S_L)$ for $S_L<0<S_R$, the upwind physical flux otherwise;
   • hllc (Toro, Spruce & Speares 1994) with the same outer speeds and the contact speed of Batten et al. (1997), $S_M={\displaystyle \frac{p_R-p_L+{\rho }_L{u}_{n,L}(S_L-u_{n,L})-{\rho }_R{u}_{n,R}(S_R-u_{n,R})}{{\rho }_L(S_L-u_{n,L})-{\rho }_R(S_R-u_{n,R})}}$; it resolves an isolated contact exactly.

3. Correction (flux_correction), coefficients derived with sympy so that the flux difference approximates the derivative of point values:

   flux_correction	$c_1$	$c_2$	$c_3$	Order of the flux difference
   6th	$64/45$	$-13/60$	$1/180$	6 (the WENO5 interpolation limits the scheme to 5)
   4th	$4/3$	$-1/6$	$0$	4
   none	$1$	$0$	$0$	2

4. Sensor (flux_correction_sensor): with weno, $s=1$ where every interpolated field keeps its nonlinear weights close to the linear ones, $\underset{k}{min}{\omega }_k/d_k\ge 0.2$ on both sides of the face, and $s=0$ elsewhere: the correction, which amplifies jumps by $64/45$, is switched off at discontinuities. none keeps it everywhere.

For MHD (riemann_solver = llf | hll | hlld) the solvers of adam_flume_mhd_riemann_library work in the frame of the direction, so the x, y and z runs of a rotated problem stay bitwise equal: LLF, HLL, and HLLD (Miyoshi & Kusano 2005) with its degenerate cases and a per-face fallback to HLL when a star state is not admissible or the wave speeds are out of order; the fallbacks are counted and logged per stage (HLLD fallbacks to HLL: N faces). With GLM the $(B_n,\psi )$ subsystem is solved exactly at the face and both states take the resulting ${\tilde{B}}_n$; with EGLM the same in field units, the cleaning energy ${\psi }^2/2$ kept out of the solver states and added to the energy flux (see EGLM). On Ryu–Jones 2a (256 cells) HLLD with characteristic interpolation reaches $L_1=3.76\cdot {10}^{-2}$ (splitting: $3.84\cdot {10}^{-2}$), with primitive interpolation $5.13\cdot {10}^{-2}$; HLL is worse with either ($4.32$, $5.98\cdot {10}^{-2}$). The MHD verification suite runs on weno-riemann through the --numerics option of its scripts (issue #47, M3-P3c): the linear waves and the circularly polarised Alfvén wave converge at order 4.96–5.00, Brio–Wu and Ryu–Jones 4d need no floor and no HLLD fallback, Orszag–Tang keeps its symmetry and conserves to ${10}^{-16}$, the field loop loses as much magnetic energy as with the splitting (primitive) or less (characteristic). Both interpolations agree on smooth waves (the WENO weights are linear there); at discontinuities characteristic interpolation is more accurate (the recommended default, issue #47 D-5), primitive keeps the 7×7 eigenvector projection out of the face kernel.

Measured (np 2, WENO5, 6th-order correction, CPU and FNL identical to the printed digits; details in the verification gallery): on the isentropic vortex the HLLC variant reaches order 6.5 between 128² and 256² with $L_1(\rho )=8.0\cdot {10}^{-8}$ (splitting: $8.4\cdot {10}^{-8}$); on Sod and Shu–Osher its $L_1$ error is 9% and 7% below the splitting scheme's, on Lax 1.6% above.

Positivity limiter

[numerics] positivity_limiter = cell (issue #47, D-9; both space schemes) makes each stage keep the density and the pressure positive, so the floors are no longer needed on hard problems. Per cell, the backbone is the forward-Euler update with the first-order Lax–Friedrichs fluxes of the cell states (Euler: Rusanov; MHD: the Wu speed, which adds $|{\mathbf{B}}_L-{\mathbf{B}}_R|/(\sqrt{{\rho }_L}+\sqrt{{\rho }_R})$ to the fastest wave, and with EGLM at least $c_h$) and, for EGLM, the second-order nonconservative sources; with the Godunov–Powell (EGLM) sources it is admissible under the CFL bound (Wu 2018, Wu & Shu 2018). Each face carries the antidiffusive difference $\pm \mathrm{\Delta }t/\mathrm{\Delta }x(F-F^{LF})$ of the high-order flux $F$ (and EGLM the difference of the high- and second-order sources). The cell factor $\mathrm{\Lambda }$ is the largest value in $[0,1]$ for which every corner of the box $[0,\mathrm{\Lambda }{]}^{\text{faces}}$ keeps $\rho$ and the internal energy above their floors; both are concave in the conservative state, so each corner has a closed-form factor (Zhang & Shu 2012), and $2^{2D}$ corners (16 in 2-D, 64 in 3-D) bound every combination. The floor of a quantity is $\kappa$ times its backbone value, $\kappa =0.1$: positivity alone leaves a state admissible but not usable (on the planar Sedov blast with HLLC the former absolute floor ${10}^{-13}$ let one stage drain the centre cell to $\rho ={10}^{-13}$ with its energy kept, the sound speed grew by ${10}^6$ and the next stage, whose $\mathrm{\Delta }t$ is that of the step's first state, had no admissible backbone), so a stage may not take a cell below a tenth of its first-order update. The floor is relative only: an absolute one would make the solution depend on the units of the input (issue #49), and below a backbone value of ${10}^{-12}$ it was the absolute one that acted. A smooth update differs from the backbone by $O(\mathrm{\Delta }x)$, so the relative floor does not act there. A face takes the smaller factor of its two cells, $F\leftarrow F^{LF}+\theta (F-F^{LF})$ (Xu 2014; Christlieb et al. 2015), so the update is conservative and each forward-Euler step of size $\mathrm{\Delta }t$ is admissible; an SSP Runge–Kutta stage is a convex combination of such steps, which is why the limiter requires an SSP scheme. A face whose two cells keep $\mathrm{\Lambda }=1$ is not touched, so smooth runs are bitwise unchanged. The cell factors are exchanged like a field (intra-realm copies and MPI); a physical-boundary face takes the interior cell's factor. The log reports, per stage, the limited faces and the cells whose backbone is inadmissible (zero in every verified case; they would fall back to the floors). A cell with a non-finite high-order face flux takes $\mathrm{\Lambda }=0$ and its faces the backbone flux itself (a NaN would otherwise pass every comparison of the corner factors and survive the blend); such cells are logged on a line of their own, and none occurs in the verified cases.

On the Balsara–Spicer blast (β = 2.5·10⁻⁴, EGLM) the splitting scheme fails at step 3 without the limiter and reaches $t=0.01$ with it (at most 4% of the faces limited per stage); weno-riemann HLLD passes with or without it.

At a 2:1 AMR seam the end-of-step reflux (below) would replace the coarse face flux with the restricted fine one, a correction that no cell factor bounds: on the blast across a refined box it drives coarse cells beside the seam to $\rho e<0$ (issue #50). A limited run therefore synchronises the seam flux at every stage, before the update:

• the fine seam faces take the backbone $F^{LF}(q_C,q_f)$, with the coarse cell $q_C$ as outer state;
• the coarse cell's factor ${\mathrm{\Lambda }}_C$ uses, on the seam face, the means of the four fine backbone and high-order fluxes; the coarse backbone update with the mean of $F^{LF}(q_C,q_f)$ is a convex combination of Lax–Friedrichs updates of $q_C$, so it is admissible;
• both sides blend with ${\theta }_s=min({\mathrm{\Lambda }}_C,min{\mathrm{\Lambda }}_f)$, so each cell stays in its corner box and the coarse flux is exactly the mean of the fine ones. The stage is conservative, and the reflux only corrects round-off.

With self-gravity the register also carries the seam gravity work (#81 P5b), so the end-of-step correction is not round-off: on a star whose surface a 2:1 seam cuts, in a ${10}^{-4}$ atmosphere, it drove one coarse atmosphere cell beside the seam to $\rho e<0$ at step 688, and the run stopped non-finite (#88). A limited run therefore gives each coarse cell the largest fraction $\theta \in [0,1]$ of its end-of-step correction that keeps $\rho$ and $\rho e$ at least $\kappa$ times their values before it (reflux_positivity_theta; $\rho e$ is concave in $\theta$, so a bisection finds the bound). A cell whose full correction is admissible receives it exactly, so a run that never clips is unchanged; a clipped cell gives up conservation by the part it did not receive, and each step that clips logs reflux positivity: +N coarse cells limited. On that star one cell was limited on five steps, and the run reached its end. reflux_cadence = stage needs no such bound: its correction acts inside the stages, before the limiter.

The coarse states, the fine means and the factors travel over the ranks in seam skins (four reductions of seam size per stage); on the GPU the seam cells are gathered to the host, which recomputes their factors and returns them with the seam face fluxes.

Dissipative fluxes (Navier–Stokes)

With a viscosity or a conductivity in [physics] (input), Euler becomes the compressible Navier–Stokes system (issue #65, P2), and so does each MHD model (P3), which also takes an Ohmic resistivity. Each term is the divergence of a flux, added to the inviscid one in direction $d$:

$$G_d^{\rho u_i}=-{\tau }_{di},G_d^E=-u_i{\tau }_{di}-k\frac{\partial T}{\partial x_d},\tau =\mu (\mathrm{\nabla }\mathbf{u}+\mathrm{\nabla }{\mathbf{u}}^T-\frac{2}{3}(\mathrm{\nabla }\cdot \mathbf{u})I),$$

with no bulk viscosity (Stokes hypothesis), $T=p/(\rho R)$, and the laws $\mu ={\mu }_0(T/T_{\mathrm{ref}}{)}^{\omega }$, $k=k_0(T/T_{\mathrm{ref}}{)}^{{\omega }_k}$ (${\omega }_k=\omega$ when $k$ comes from the Prandtl number, otherwise 0; the constant law is $\omega =0$). The pointwise flux is one pure routine shared by both backends (adam_flume_dissipation_library).

Ohmic resistivity (MHD). With the magnetic diffusivity $\eta$ and $\mathbf{J}=\mathrm{\nabla }\times \mathbf{B}$ (code units), Ohm's law $\mathbf{E}=-\mathbf{u}\times \mathbf{B}+\eta \mathbf{J}$ in ${\partial }_t\mathbf{B}=-\mathrm{\nabla }\times \mathbf{E}$ and the Poynting flux $\mathbf{E}\times \mathbf{B}$ give

$$G_d^{B_i}=-\eta (\frac{\partial B_i}{\partial x_d}-\frac{\partial B_d}{\partial x_i}),G_d^E=\eta (\mathbf{J}\times \mathbf{B}{)}_d=\eta B_j(\frac{\partial B_d}{\partial x_j}-\frac{\partial B_j}{\partial x_d}),$$

with $G_d^E=\sum _i{B}_i{G}_d^{B_i}$, so the field's energy loss is exactly the Ohmic heating $\eta |\mathbf{J}{|}^2$ that the internal energy receives. The field flux is the curl form, ${ϵ}_{idk}\eta J_k$, not $-\eta {\partial }_d{B}_i$: the two differ by $\eta \mathrm{\nabla }(\mathrm{\nabla }\cdot \mathbf{B})$, which is not zero under GLM, and the curl form keeps the resistive term out of the divergence constraint. (The plan in #65 printed the field flux with the opposite sign, an anti-diffusion; VV-5 measures the decay rate $\eta k^2$.) The Ohmic fluxes have their own kernels, with the same stencils on $W=\mathbf{B}$, called only when the resistivity is on; the viscous and heat kernels serve the MHD models unchanged ($(u,v,w,T)$ sit at the same auxiliary indexes). The psi equation and the EGLM sources are untouched.

Where it enters. The dissipative face fluxes are computed by their own kernels, after the inviscid face fluxes and before the seam accumulation and the flux difference, in both space schemes. Two consequences: the AMR reflux and the inter-realm flux register see the total flux, so conservation across seams needs no new machinery; and the WENO face kernels are untouched (their FNL register count does not change). A run without coefficients never calls the kernels and stays bitwise identical to the ideal solver. The kernels read $(u,v,w,T)$ from the auxiliary variables, whose indexes are common to every model, and the host picks the kernel of the order ([numerics] dissipative_order), never a branch inside a loop.

Face derivatives. The scheme stores point values, so a face flux of order $p$ does not by itself give a flux difference of order $p$: the conservative numerical flux is $h=G-\frac{\mathrm{\Delta }{x}^2}{24}{G}^{″}+O(\mathrm{\Delta }{x}^4)$ (Shu & Osher 1989). Face $i+\frac{1}{2}$ along $d$, cells $m=-1\dots 2$ counted from the left one:

	dissipative_order = 2	dissipative_order = 4 (default)
face state	$(W_0+W_1)/2$	$(9(W_0+W_1)-(W_{-1}+W_2))/16$
normal derivative	$(W_1-W_0)/\mathrm{\Delta }x$	$(27(W_1-W_0)-(W_2-W_{-1}))/(24\mathrm{\Delta }x)$
tangential derivative	mean of the 2nd-order central differences of cells 0 and 1	the 4th-order central differences $D_m$ of cells $-1\dots 2$, interpolated as the state
flux	$G(W_f,\mathrm{\nabla }{W}_f)$	$G(W_f,\mathrm{\nabla }{W}_f)-(G_2-G_1-G_0+G_{-1})/48$

where $G_m$ is the flux of cell $m$ from 2nd-order central gradients, so the correction term is $\mathrm{\Delta }{x}^2{G}^{″}/24$ to second order. The same form serves the linear and the nonlinear terms ($u_i{\tau }_{di}$, $\mu (T)$). Across a strong temperature jump the four-point face temperature undershoots, possibly below zero, where a power law $(T/T_{\mathrm{ref}}{)}^{\omega }$ is NaN (a wall 300 times hotter than the gas did it in the input-contract check DC); where it is not positive, the face takes the mean of its two cells for the laws, which a smooth positive $T$ never triggers. For constant $\mu$ it is a 4th-order $u_{xx}$ on six points, not the compact five-point one; its Fourier symbol is non-positive at every wavenumber ($-4.67\nu /\mathrm{\Delta }{x}^2$ at the grid Nyquist mode against $-5.33$ for the compact operator), so it damps every mode, and its error is $O({\theta }^6)$ in the wavenumber $\theta$. The 4th-order fluxes read cells $-2\dots 3$ along $d$ and $\pm 2$ along the tangents, so they need [grid] ngc $\ge 3$ (fatal otherwise); the tangential differences read the edge ghosts, whose fill is exact on every path (the ghost probe GP, verification). A null direction has a zero step, so its derivatives are zero exactly and its ghosts are never read.

Walls. wall-noslip and wall-isothermal (boundary conditions) are mirrors, second-order at the wall: a 4th-order interior with them converges at order 2 overall (VV-3).

Time step. The explicit unsplit integration adds the diffusive limit to the cell sum of the CFL bound,

$$\mathrm{\Delta }t=\frac{\mathrm{CFL}}{\underset{\text{cells}}{max}[\sum _d\frac{|u_d|+a}{\mathrm{\Delta }{x}_d}+2\nu \sum _d\frac{1}{\mathrm{\Delta }{x}_d^2}]},\nu =max(\frac{4}{3}\frac{\mu }{\rho },\gamma \frac{k}{\rho c_p},\eta ),$$

so under AMR the diffusive part tightens four times per level, against twice for the hyperbolic one. With a power law and an isothermal wall the laws in $\nu$ are evaluated at $max(T,T_{w,max}^2/T)$, the geometric-mirror ghost temperature of a cell beside the hottest wall: a wall much hotter than the gas gives its face a coefficient above every cell's, which the cell bound alone misses (a wall 330 times hotter than the gas blew up at CFL 0.5). The first step logs both maxima, the minimum cell Reynolds number $\sum _d(|u_d|+a)/\mathrm{\Delta }{x}_d/(\nu \sum _d\mathrm{\Delta }{x}_d^{-2})$ ($(|u|+a)\mathrm{\Delta }x/\nu$ on a uniform 1-D grid), and whether the run is hyperbolic- or diffusion-limited.

Limits. Explicit only: a diffusion-limited run pays $\mathrm{\Delta }t\propto \mathrm{\Delta }{x}^2$ (super-time-stepping is a follow-up). The positivity limiter is refused with dissipative terms (D-M4-5), and so are immersed solids, whose walls are inviscid. A null direction freezes the momentum normal to it (the CHASE semantics of the flux difference), so a reduced-dimension run cannot carry a velocity along a null axis: a shear layer or a Couette flow needs its velocity axis active, if thin.

Across seams. The dissipative fluxes enter the AMR reflux and the inter-realm register with the inviscid ones, so a 2:1 face conserves to round-off (VV-7). A 2:1 face is second order, for them as for the inviscid fluxes: the mean restriction of the coarse ghosts and the reflux of point-value fluxes each cap the order at 2 (issue #68 F1); without the reflux the dissipative flux at the face is first order (VV-8).

Time integration

The library Runge–Kutta schemes ([runge_kutta] scheme, listed in the input reference) integrate the semi-discrete system: strong-stability-preserving (SSP) and low-storage schemes. Each stage refills the ghost cells, recomputes the auxiliary variables, applies the positivity floors (MHD) and evaluates the residual. The time step follows the CFL bound, with the diffusive part of the dissipative fluxes when they are on; [time] it_max and time_max stop the run.

Adaptive mesh refinement

The block tree is refined at initialisation ([initial_conditions] amr_iterations passes) by the markers of the [amr] section: a geometric box, a variable gradient, or the surface of an immersed solid. Neighbouring blocks differ by at most one level (2:1 balance). At a coarse–fine face:

• the ghost cells of the fine block are interpolated from the coarse one (the library seam ghost fill), and the coarse ghosts are restricted from the fine cells. The interpolant (tricubic by default) works variable by variable, so next to a strong shock it can return an inadmissible state: at low β the internal energy is a small residual of $E-|\mathbf{B}{|}^2/2$, and on the Balsara–Spicer blast across a 2:1 seam the fine ghosts reach $\rho e<0$, which makes the high-order flux non-finite (issue #50). After every ghost fill, an inadmissible face ghost is therefore blended toward the adjacent interior cell, $g\leftarrow q_{\text{in}}+t(g-q_{\text{in}})$, with the largest $t$ that keeps $\rho$ and $\rho e$ above the limiter's floors of $q_{\text{in}}$ (a non-finite ghost takes $q_{\text{in}}$). An admissible ghost is never touched, so the blend changes nothing else and costs one scan of the face ghosts; the log reports the ghosts blended, per rank;
• conservation is restored by the Berger–Colella reflux: the fluxes of every Runge–Kutta stage through the seam are accumulated in a flux register, weighted by the stage's coefficient in the step, and at the end of the step the coarse cells next to the seam are corrected by the difference between the fine fluxes and their own ([numerics] reflux = .true.; with reflux_cadence = stage the correction is applied at every stage instead). Under GLM $\mathbf{B}$ and $\psi$ are refluxed; under EGLM $\mathbf{B}$ is not, since the sources $-\mathbf{u}\mathrm{\nabla }\cdot \mathbf{B}$ make it a non-conserved integral (G3b). With reflux the volume integrals of a periodic box stay constant to round-off, without it they drift at ${10}^{-5}$ (conservation example).

With markers the tree is an octree (ratio = 8, and with a null axis nk ≥ 4) or a quadtree (ratio = 4, any nk, nk = 1 for a true 2-D run): a quadtree seam is 2:1 in x and y and 1:1 in z, and the ghost restriction, the coarse-fine interpolation, the reflux and the limiter's seam synchronisation work per axis (issue #46, verified by MV-15). For MHD, a cell-centred $\mathbf{B}$ cannot be kept divergence-free across a 2:1 face (Tóth & Roe 2002): the seams inject a truncation-level $\mathrm{\nabla }\cdot \mathbf{B}$ that GLM transports and damps; the div(B) history reports it (seam_max_divb).

Runtime regridding (#74). With [amr] frequency = n > 0, after every $n$-th step and before its output, the realm regrids: up to iters sweeps, each one refreshing the ghost cells, combining the markers into one set of flags (a block is refined when any marker asks for it, coarsened only when all agree, never below the base level), and calling the library regrid (2:1 balance, prolongation, restriction, redistribution, maps). The distance function of the immersed solids is recomputed on the new grid, the new state is checked (density and pressure positive, fatal otherwise, with the cell), and the forest rebuilds the 2:1 flux register, so the reflux follows seams that appear and disappear. The time step stays global (no subcycling). The FNL backend regrids by a host round trip: each sweep refreshes the ghosts on the device, copies the state to the host, marks and regrids there with the same code, and copies the state, the block coordinates, the maps and the distance function back (a map the new grid no longer has is freed on the device); its log line adds the wall time. The state and distance-function copies move the blocks in use only: the device array leads with the block index, so a kernel packs those blocks, transposed to the host layout, into a scratch buffer that one contiguous copy moves (#75; copying the whole capacity cost about 2.5 s per regrid on the development box, two ranks of 17195 blocks, whatever the grid, against 0.02 s now on a 52-block grid). The Löhner marker (Löhner 1987, mode 4) is the scale-free second-derivative estimator

$$E=\sqrt{\frac{\sum _d{(u_{+}-2u+u_{-})}^2}{\sum _d{(|u_{+}-u|+|u-u_{-}|+ϵ(|u_{+}|+2|u|+|u_{-}|)+\delta )}^2}}\in [0,1],$$

in the dimension-sum form of FLASH and PLUTO ($u_{\pm }$ the neighbours along $d$): about 0 where the variable is smooth on the grid, about 1 at a jump. The epsilon term filters ripples small relative to the variable; on a variable that vanishes in quiet regions it filters nothing, and round-off noise reads $E\approx 1$ (the GLM residue of $|\mathbf{B}|\sim {10}^{-6}$ outside a field loop refined the whole box). The absolute floor $\delta$ (default 0, the FLASH/PLUTO form) makes the variations well below it noise. A block takes the largest $E$ over its cells and buffer ghost layers; it is refined above refine_tol, coarsened below derefine_tol, kept between.

Keep shocks off the 2:1 seams. A captured shock that crosses a 2:1 seam into finer cells leaves a defect behind it (about 2 % of the jump in density, with an entropy defect of the same size) that does not shrink with the resolution and that no choice of seam operator removes: restriction, coarse-to-fine fill and reflux move it by at most ±30 % (A3, E0). The evidence points to the start-up error of a captured shock whose numerical structure is squeezed into finer cells. Refine ahead of the shocks instead: a gradient or Löhner marker on density or pressure with a buffer of a few cells and a regrid frequency short enough that a shock cannot outrun the refined region between two regrids, so shocks stay inside the finest level and never cross a seam (the standard AMR practice; E0 measured the defect, not this remedy).

Regridding transfers. When a block is derefined, the parent takes the mean of its children. When one is refined, [amr] regrid_prolongation chooses how the children are filled. FLUME's default is conservative: each child is the parent plus monotonized-central slopes (van Leer 1977) times a quarter of a cell, one per refined axis, scaled by the largest $\varphi \le 1$ that keeps every child within the parent's and its face neighbours' range (Barth & Jespersen 1989). The children average to the parent, so a regrid changes no conserved integral and a refine followed by a derefine returns the state exactly; a positive variable stays positive; linear data are reproduced. The library's linear prolongation (tensor interpolation, weights 1/4 and 3/4) is second order too but not conservative (RG). Both read the parent's first ghost layer.

Immersed boundary (Euler)

Static rigid solids ([solids], [solid_N], e.g. definition = analytical_circle) are imposed as inviscid walls:

1. a signed distance function $\varphi$ marks the fluid and the solid cells;
2. the fluid state is extended into the solid by n_eikonal sweeps of an eikonal extrapolation along the surface normal $\mathrm{\nabla }\varphi$ (Jacobi sweeps: every increment is computed from the old state, so the result does not depend on the loop order or the thread count, and the CPU and FNL backends agree);
3. the extended state is inverted into the mirror state of a slip wall (the normal velocity reversed);
4. the fluid cells cut by the surface use the cut spacing in the flux difference, and the solid cells are masked in the Runge–Kutta stages. The cut spacing of a fluid cell $c$ with a solid neighbour $s$ is $\mathrm{\Delta }/2-{\varphi }_c\mathrm{\Delta }/({\varphi }_s-{\varphi }_c)$; the denominator is positive (${\varphi }_s>0>{\varphi }_c$), so it carries no guard (the absolute ${10}^{-12}$ inherited from CHASE was a length, removed with issue #49).

The solid surface can drive the AMR (the refined ring of the shock-over-a-cylinder example). mhd-ideal with solids is refused.

Multi-realm runs

A forest manifest glues several realms, each with its own input file, through inter-realm seams (mirror coupling of blocks of the same resolution). The seam ghosts are filled from the peer realm at every Runge–Kutta stage (stage_coincident) or once per step (end_of_step); the union of a split problem is bitwise equal to the single-realm run (verification multirealm/).

Output and diagnostics

Checkpoints are XH5F files (one HDF5 file per rank plus an XDMF description readable by ParaView), written every it_save steps, with optional slices, residual and conservation histories, and restart files. MHD adds the div(B) history <basename>-divb_history.dat (it time max_divb l1_divb seam_max_divb, plus the D_flux columns with [diagnostics] divb_flux) and the divb_tol monitor, on the normalised D_flux max h|D_flux| / max|B| (#68 G1).

References

• Barth T. J., Jespersen D. C. (1989), The design and application of upwind schemes on unstructured meshes, AIAA paper 89-0366.
• Batten P. et al. (1997), On the choice of wavespeeds for the HLLC Riemann solver, SIAM J. Sci. Comput. 18, 1553–1570.
• Berger M. J., Colella P. (1989), Local adaptive mesh refinement for shock hydrodynamics, J. Comput. Phys. 82, 64–84.
• Chen Y., Tóth G., Gombosi T. I. (2016), A fifth-order finite difference scheme for hyperbolic equations on block-adaptive curvilinear grids, J. Comput. Phys. 305, 604–621.
• Christlieb A. J., Liu Y., Tang Q., Xu Z. (2015), Positivity-preserving finite difference weighted ENO schemes with constrained transport for ideal magnetohydrodynamic equations, SIAM J. Sci. Comput. 37, A1825–A1845, doi:10.1137/140971208.
• Einfeldt B. et al. (1991), On Godunov-type methods near low densities, J. Comput. Phys. 92, 273–295.
• Don W. S., Li R., Wang B.-S., Wang Y. H. (2022), A novel and robust scale-invariant WENO scheme for hyperbolic conservation laws, J. Comput. Phys. 448, 110724.
• Jiang G.-S., Shu C.-W. (1996), Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126, 202–228.
• Becker R. (1922), Stoßwelle und Detonation, Z. Phys. 8, 321–362.
• Löhner R. (1987), An adaptive finite element scheme for transient problems in CFD, Comput. Methods Appl. Mech. Engrg. 61, 323–338.
• Miyoshi T., Kusano K. (2005), A multi-state HLL approximate Riemann solver for ideal MHD, J. Comput. Phys. 208, 315–344.
• Shu C.-W., Osher S. (1989), Efficient implementation of essentially non-oscillatory shock-capturing schemes II, J. Comput. Phys. 83, 32–78.
• Toro E. F., Spruce M., Speares W. (1994), Restoration of the contact surface in the HLL-Riemann solver, Shock Waves 4, 25–34.
• Tóth G., Roe P. L. (2002), Divergence- and curl-preserving prolongation and restriction formulas, J. Comput. Phys. 180, 736–750.
• van Leer B. (1977), Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection, J. Comput. Phys. 23, 276–299.
• Wu K. (2018), doi:10.1137/18M1168017; Wu K., Shu C.-W. (2018), doi:10.1137/18M1168042: positivity of the first-order and high-order schemes for ideal MHD with the Godunov–Powell source.
• Xu Z. (2014), Parametrized maximum principle preserving flux limiters for high order schemes solving hyperbolic conservation laws, Math. Comp. 83, 2213–2238, doi:10.1090/S0025-5718-2013-02788-3.
• Zhang X., Shu C.-W. (2012), Positivity-preserving high order finite difference WENO schemes for compressible Euler equations, J. Comput. Phys. 231, 2245–2258, doi:10.1016/j.jcp.2011.11.020.