6  Equations and numerical methods

6.1 GLM magnetohydrodynamics

The MHD state stores nine conserved fields in structure-of-arrays order (rho, rho*vx, rho*vy, rho*vz, E, Bx, By, Bz, psi). Code magnetic units absorb 4*pi, so magnetic pressure is B^2/2, not B^2/(8*pi). ^{allocate_state}

The local ideal-MHD flux operator has continuity, momentum, and induction terms

\partial_t\rho+\nabla\cdot(\rho\mathbf v)=0, \qquad \partial_t(\rho\mathbf v)+\nabla\cdot\left[ \rho\mathbf v\mathbf v+(p+B^2/2)\mathbf I-\mathbf B\mathbf B\right] =\mathbf S_m, \tag{6.1}

\partial_t\mathbf B+\nabla\cdot(\mathbf v\mathbf B-\mathbf B\mathbf v) +\nabla\psi=0, \qquad p=(\gamma-1)\left(E-\frac{\rho v^2}{2}-\frac{B^2}{2}\right). \tag{6.2}

Here the induction tensor is understood componentwise as the flux v_j B_i - B_j v_i in direction j. The GLM scalar propagates magnetic divergence at the cleaning speed and is damped by the mixed hyperbolic/parabolic cleaning update. The energy flux is (E+p+B^2/2)*v - (v dot B)*B; external gravity, drag, expansion-aware chemistry, and radiation response are coupled by their respective source updates, not inserted as an unqualified isolated-box conservation law. ^glm

The production scheme uses finite-volume Godunov updates with PPM reconstruction and HLLD (Miyoshi and Kusano 2005) with GLM cleaning (Dedner et al. 2002) unless configured otherwise; PLM is used by the tiny installation examples. LLF fallback and positivity handling can change the local numerical diffusion and must be monitored. A portable reference integrator and fused cube integrator coexist. Shared source is not a substitute for backend parity tests.

6.2 Constrained transport (optional)

mhd_solver = "ct" replaces the GLM divergence cleaning by face-staggered constrained transport (Gardiner and Stone 2005); GLM stays the default. The Godunov sweeps are shared with GLM: PLM/PPM reconstruction, the Hancock predictor, and HLLD/HLL/LLF fluxes on cell-centred conserved variables, with the normal magnetic component at each face overwritten by the staggered face value and psi identically zero. The induction update differs:

1. The Godunov magnetic fluxes on the lower faces give face EMFs, and a cell-centred reference EMF E = -v x B is built from the source-updated predicted velocity and the step-start B.
2. Edge EMFs are the Gardiner and Stone CT-contact average of the four adjacent face EMFs, with upwind gradient corrections selected by the sign of the mass flux and anchored on the cell-centred reference EMF.
3. The face fields bxf, byf, bzf advance by the discrete curl of the edge EMFs, so the face divergence is preserved to round-off.
4. The cell-centred B used by every other module (energy, forces, chemistry, radiation, diagnostics) is the average of the two faces.

Compensated face updates. With a strong Float32 guide field the per-step face increments dtdx curl(E) of a weak perturbation are comparable to half an ulp of the face value, and plain summation drops them face by face; because the dropped parts differ from face to face, the exact face divergence then random-walks away from round-off and the perturbation drifts from the Float64 solution (the test_ct.jl guide-field case: components 22 to 35, 198 ulp perturbation, 2000 steps, div_face*dx/Brms 1.15e-7 initially, 7.5e-6 plain, 1.5e-7 compensated). When MHD_RADIATION_MAGNETIC_COMPENSATION=1 the radiation-coupled steps (exactly where GLM compensates its cell B) therefore sum each face with Kahan compensation: the part of an increment that b + inc cannot represent is kept in a per-face remainder (rbx, rby, rbz, three more slots of the CT buffer, checkpointed as ct_residual_1..3) and added to the next increment. The exact sum stored + remainder follows the curl update, so the divergence of the stored faces stays bounded by the remainders (at most half an ulp of each face) instead of growing. The cell-centred B remains the plain average of the stored faces, and plain step! and drag-only steps stay uncompensated, as for GLM. ^{_ct_faceB_comp_k!}

The scheme is a port of ct_mhd.jl from the Solaceonx/GLMMHDTurb sweep-fixes (commits 9d20e52 and 7256715). ^_ct_gs05

Initial conditions are generated as a cell field. Once all initial-condition and import stages are complete, the driver builds divergence-free faces from it once with the minimal change of the cell field (projection onto the face-solenoidal subspace with the tan(k dx/2) symbol), keeps the thermal pressure fixed, and logs CT_INITIAL_FACES cell_b_rel_l2_change, magnetic_energy_ratio, and divb_face_dx_over_brms. Because the PMF generators make the field solenoidal under another discrete symbol (centred sin 2θ, or continuum θ), this changes the field by O((k dx)^2) near the cutoff. The measured rms (brms_measured) reports the post-projection cell field; the reference-band normalization values remain those of the generated field. ^{ct_faces_from_cells!}

Terminal and hybrid drag. The overdamped paths advance the faces. The cell EMF \mathbf v\times\mathbf B is computed exactly as for GLM (cell-centred \mathbf v from the Lorentz force and pressure response, cell \mathbf B = face average); each edge EMF is the mean of the four cells sharing the edge, and the faces advance by the discrete curl of the edge EMFs. The SSPRK2 terminal_ct_induction!, the pressure-implicit step and the radiation terminal induction keep the GLM stages, subcycling, cfl and eta; the order-2/order-4 dissipation is the Laplacian or bilaplacian of every face component with the same stencil, which commutes with the face divergence on a periodic grid. After each stage cell \mathbf B is the face average and E changes by the magnetic-energy change (as GLM). The first-order terminal velocity steps and the implicit B diffusion (a Fourier multiplier applied to each face component) do the same; the B projection of the terminal paths is skipped because the face divergence is already round-off. Face sums are Kahan-compensated through ct.rbx/rby/rbz, because plain Float32 summation of sub-ulp increments over 1000 Euler steps raised \Delta x\,\max|\nabla\cdot B|/B_{rms} to 1.6e-5 (1.4e-7 compensated). Genuinely cell-wise options are errors for CT: MHD_TERMINAL_INDUCTION=ssprk2 (cell SSPRK2 with a GLM-style projection) and, with MHD_TERMINAL_PRESSURE_IMPLICIT=1, any induction other than ct. ^{terminal_ct_induction!}

Limits. CT is periodic-only and implemented on the ref integrator only; the cube integrators and the GLM projection are rejected (terminal and hybrid drag are supported, see below). Metal allows 31 buffer arguments per kernel and counts each scalar argument as one, so the packed reference kernel (step_ref_packed_kernel!, generated from the same body as the CPU GLM step_ref_kernel!) takes the grid/BC integers and the 13 scalars as tuples (and, for CT, the CT arrays as one packed buffer); it runs on Metal for both solvers. GLM :ref on a GPU backend launches this kernel; the CPU keeps the historic GLM kernel, whose results are bit-for-bit unchanged. CUDA is NOT_RUN. CT allocates 21 extra N^3 arrays (three face fields, 15 flux/EMF workspace arrays and three face remainders, in one buffer) on top of GLM’s 18. From a state at rest the face mass fluxes are rounding noise and the upwind choice of the edge EMF follows their sign, so CT trajectories on different devices separate at first faster than with a definite flow (Metal vs CPU Float32, 32^3 PMF field, 20 steps: 1.2e-4 relative from rest, 1.3e-6 with a uniform flow). The explicit step is not limited by the GLM cleaning speed. divb_face_dx_over_brms in history.tsv and performance.tsv is the exact face divergence (empty/NaN for GLM); the centred divBdx_over_brms is a different diagnostic and is not zero for CT.

Positivity protection of the radiation-coupled step. GLM protects its cube radiation step: the trial state is checked (density at least max(smallr, llf_dmin), finite state, internal energy at least pfl/(gamma-1), trial wave speed at most positivity_wave_growth_limit * ch), a rejected trial restores the pre-step state and is retried as a complete minmod-PLM/HLL step, and a second rejection either raises RADIATION_CHECKPOINT_CFL_FALLBACK (checkpoint fallback) or halves the step recursively (positivity_max_subdivisions). CT has the same protection on ref (^{step_radiation_godunov!}), with identical counters (fallback_steps, plm_hll_fallback_steps, subdivision_steps) and warnings. Differences from GLM: the face fields and their Kahan remainders are updated in place by ct_update!, so a snapshot of them (6 extra N^3 arrays in the radiation workspace, ct_snapshot: 15% of the 39 state and CT arrays, plus one device copy per radiation step) is taken before the trial and restored exactly with the cell fields (GLM resets its remainders; CT does not lose them); the minmod-PLM/HLL retry passes its reconstruction and Riemann codes to the kernel without changing s.recon/s.riemann; and the internal-energy and finiteness checks run on the stored trial state (the ref kernel has no inadmissible-cell flag). The GLM ref path is deliberately unchanged and has no such protection (it is the CPU default, not the production GPU path). CT callers must pass a positive ch (the wave-speed gate is growth * ch; CT dynamics do not use it). Runner: CT passes ch = smax, so MHD_RADIATION_POSITIVITY_WAVE_GROWTH_LIMIT is multiplied by MHD_GLM_CH_FAC for CT to reproduce GLM’s absolute wave-speed limit. Whether the PLM/HLL trial is rejected can differ between GLM and CT in low-beta gas, where the thermal energy is a small difference of large energies and the schemes treat B differently (comparing CT and GLM runs).

6.3 Magnetic-force decomposition

The centered Maxwell-stress diagnostic uses

\mathbf F_P=-\nabla(B^2/2),\qquad \mathbf F_T=\nabla\cdot(\mathbf B\mathbf B),\qquad \mathbf F_L=\mathbf F_P+\mathbf F_T. \tag{6.3}

The divergence form equals (B dot grad) B only when magnetic divergence is negligible. The force used by the radiation coupling and diagnostic split uses the same centered stencil. ^{terminal_magnetic_force_split!}

6.4 Expansion, units, and clocks

The run evolves a=1/(1+z) and a dimensionless supercomoving clock tau with

\frac{da}{d\tau}= \sqrt{a^3\left(\Omega_m+\Omega_\Lambda a^3+\Omega_k a+\Omega_r/a\right)}, \qquad \Omega_k=1-\Omega_m-\Omega_\Lambda-\Omega_r. \tag{6.4}

Radiation density is included in early starts. The driver converts this clock to physical source times and MHD code time; it does not equate an MHD step to a physical second. legacy and jedamzik clock scalings are distinct model choices, with the facade using the latter for cosmological runs. ^{dadtau function _hydro_time_per_tau}

For a comoving field amplitude in nG specified at z_ref, the initializer uses

B_{\rm phys}=10^{-9}B_{0,\rm nG}(1+z_{\rm ref})^2\ {\rm G},\qquad B_{\rm code}=\frac{B_{\rm phys}}{v_{\rm unit}\sqrt{4\pi\rho_b(z_{\rm ref})}}. \tag{6.5}

3 pG = 0.003 nG. The physical density follows the selected baryon cosmology. MHD_CHEM_VUNIT is in cm/s; it is not a box velocity in km/s. The pressure initializer with pressure="auto" uses the same CMB convention as the chemistry, T_CMB=2.725(1+z) K, including particle-number contributions from helium and electrons. ^{function _pmf_code_brms_from_comoving_ng function _cosmological_ionized_p0}

The TOML box is ckpc/h. The facade divides it by h for the driver’s comoving ckpc length. The diagnostic spectrum has a separate numeric box parameter; see wavenumber units before plotting it.

6.5 Photon-baryon response

The evolved-moment closure retains photon density and three photon velocity components. In its linear Fourier response subproblem, consistent code-time derivatives obey

\begin{aligned} \delta_b'&=-i\mathbf k\cdot\mathbf v_b,\\ \mathbf v_b'&=-i\mathbf k c_s^2\delta_b+ \Gamma(\mathbf v_\gamma-\mathbf v_b)+\mathbf F_{\rm MHD},\\ \delta_\gamma'&=-\tfrac43 i\mathbf k\cdot\mathbf v_\gamma,\\ \mathbf v_\gamma'&=-\tfrac14 i c^2\mathbf k\delta_\gamma+ R\Gamma(\mathbf v_b-\mathbf v_\gamma)-S_\gamma(k)\mathbf v_\gamma, \end{aligned} \tag{6.6}

where R=3*rho_b*c^2/(4*rho_gamma_energy) and Gamma is the converted Compton-drag rate. Equal and opposite scattering momentum transfers recover the tight-coupling sound speed c^2/[3(1+R)]. An algebraic nonlocal stress closure replaces independent quadrupole fields, approaching diffusion damping at small k*lambda_gamma and phase mixing toward free streaming. ^{_photon_stress_rates}

moments uses a second-order L-stable SDIRK response; moments_exponential is the exponential variant chosen by the facade. They should not be described as identical integrators. Legacy angular and bridge models remain comparison options, not higher-order photon hierarchies.

The full MHD update remains active. Spectral predictors and correctors adjust both density and velocity relative to local drag, avoiding double-counting its effect. Default photon ICs are adiabatic and comoving with gas; controls are MHD_RADIATION_PHOTON_ADIABATIC and MHD_RADIATION_PHOTON_COMOVING. ^{step_radiation_godunov!}

The positive MHD_RADIATION_FREE_STREAM_TOL optimization can omit the spectral gas correction when its bound is small, while still evolving photons. Compare against zero tolerance before interpreting its accuracy or speed.

6.6 Gravity and chemistry

The periodic Poisson solve removes the homogeneous source and combines gas and lattice-DM perturbations with the configured baryon/DM weights. The default midpoint density source and particle updates are time-centered. Lattice displacements reduce cancellation when tiny density perturbations are evolved; they are not an imported cosmological particle realization. ^{gravity _assemble_total_perturbation_k!}

The full mixing chemistry couples hydrogen/helium abundances, thermal energy, and recombination response. The facade uses symmetric half-steps around MHD and an expansion-aware source clock. Centered hydrogen storage protects small perturbations relative to a changing homogeneous ionization reference; it must be reconstructed before interpreting an abundance array as x_HII. ^{chem_hydrogen_storage}

Gaussian smoothing approximates mixing. Optional nonlocal response tables add a separately modeled, support-limited response. Neither path solves full three-dimensional Lyman-alpha transfer; validation limits are part of the model specification.

6.7 Baryon–DM streaming and planar studies

A particle velocity is not an MHD velocity. The periodic particle solver stores u_c=dx/d\tau with d\tau=H_0dt/a^2, so its physical peculiar velocity is v_c=(H_0 L/a)u_c, where L is the physical comoving box length, independent of the spectral coordinate box. In Jedamzik scaling the gas conversion is v_b=v_{\rm unit}\sqrt{a_{\rm ref}/a}\,u_b. The gas momentum decay u_b\propto a^{-1/2} and constant force-free particle u_c both give physical v\propto a^{-1}. Their gravitational impulses agree after these conversions. ^{_particle_velocity_cms_per_code} baryon-DM common physical clocks

The optional [streaming] components vx_kms, vy_kms, and vz_kms specify physical peculiar baryon minus DM velocity at cosmology.zstart. The built-in initial state uses the DM rest frame: gas receives this uniform boost, preserving its internal energy, and DM begins at rest. No normalization epoch is inferred. Photons follow the existing comoving-gas initialization unless MHD_RADIATION_PHOTON_COMOVING=false is explicitly selected. That choice changes the photon rest frame and must be included in an experiment’s specification. Streaming requires cosmology, full gravity with gas kicks, and Jedamzik scaling; a boost is rejected with an imported correlated state to avoid adding it twice. ^{boost_gas_velocity!}

Tseliakhovich & Hirata (2010) describe the coupling of coherent relative flows to short density modes. The supplied boost is a local controlled background, not a draw from their cosmological distribution. The built-in PMF gas/DM initialization is uniform in density and has no primordial adiabatic transfer realization. Its growing DM perturbations are PMF-driven. A cosmological suppression measurement needs correlated density, velocity, thermal and photon initial conditions and matched zero-stream controls. It also needs phase-sensitive density growth at common epochs and separate timestep/grid convergence. A plot of streaming RMS cannot establish that measurement.

tools/reference/streaming.py supplies an independent REFERENCE_ONLY moving-background ODE in the DM frame, with matter-only expansion and prescribed adiabatic sound speed. Its analytic growing-mode, transverse-stream and conjugate phase tests check conventions; they do not certify the nonlinear coupled solver. test_streaming.py

magnetic.invariant_axis=1|2|3 retains only magnetic modes with zero wavenumber along that axis before normalization. All three magnetic and velocity components remain active. This defines planar (2.5D) initial conditions evolved by the existing 3D periodic cube kernels. It still allocates and updates N cubed cells; it is not an N-by-N-by-1 accelerated solver. The planar spectrum is a different ensemble from isotropic 3D Batchelor statistics, so do not interpret its 3D shell estimates as an isotropic cosmological spectrum. Uniform streaming, uniform chemistry, and the periodic lattice respect this symmetry analytically; tests check numerical invariance. Imported states require their own symmetry checks. ^{pmf_batchelor_field} streaming_checks.jl

Dedner, A., F. Kemm, D. Kröner, C.-D. Munz, T. Schnitzer, et al. 2002. “Hyperbolic Divergence Cleaning for the MHD Equations.” Journal of Computational Physics 175: 645–73. https://doi.org/10.1006/jcph.2001.6961.

Gardiner, Thomas A., and James M. Stone. 2005. “An Unsplit Godunov Method for Ideal MHD via Constrained Transport.” Journal of Computational Physics 205: 509–39. https://doi.org/10.1016/j.jcp.2004.11.016.

Miyoshi, Takahiro, and Kanya Kusano. 2005. “A Multi-State HLL Approximate Riemann Solver for Ideal Magnetohydrodynamics.” Journal of Computational Physics 208: 315–44. https://doi.org/10.1016/j.jcp.2005.02.017.