Primitive-Variable Recovery

Relativistic finite-volume schemes update conserved variables but need density, pressure, velocity, and magnetic fields to reconstruct fluxes. Recovering these primitives is a nonlinear inversion that can become ill-conditioned at high Lorentz factor, high magnetization, low internal energy, or near EOS boundaries. A robust implementation distinguishes valid difficult states from conserved states that no physical primitive can produce, reports fallbacks and repairs, and validates round trips across the intended parameter domain.

Helpful background. Formulations of General Relativity supplies the spacetime split, Ordinary Differential Equations supplies nonlinear root solving, Relativistic Matter and Compact Objects supplies EOS and material regimes, Computational Relativity supplies evidence standards, and Fluid and Kinetic Theory supplies fluid variables and physical admissibility.

The inverse map

In pure hydrodynamics, conserved variables $D$, $S_i$, and $\tau$ are nonlinear functions of $\rho$, internal energy $\epsilon$, pressure $p$, and velocity. In GRMHD, magnetic terms couple to the Lorentz factor and enthalpy. Many algorithms introduce

$$Z=\rho hW^2$$

and reduce the inversion to one or two scalar equations for $Z$, $W$, pressure, or temperature, using the EOS at each iteration. Other schemes solve a higher-dimensional system. The implementation must state its unknowns, residual scaling, brackets, stopping condition, derivative evaluation, maximum iterations, and fallback chain.

Noble et al. 2006, §§3–5 and §5.1, Tables 2–3 compares several GRMHD inversions over millions of trial states. Its success rates depend on the sampled distribution and EOS; they are not universal guarantees.

A scalar hydrodynamic recovery

For a nonmagnetic $\Gamma$-law gas, use undensitized $D,S_i,\tau$ as defined on Conservative Relativistic Hydrodynamics. If the code stores $\sqrt\gamma$ times these fields, divide by $\sqrt\gamma$ first. With $S^2=\gamma^{ij}S_iS_j$, a trial pressure determines

$$Q=\tau+D+p,\qquad v^2=\frac{S^2}{Q^2},\qquad W=(1-v^2)^{-1/2},\qquad \rho=\frac{D}{W}.$$

Using $h=1+\Gamma\epsilon$ and $p=(\Gamma-1)\rho\epsilon$ reduces the inverse problem to

$$f(p)=p-\frac{\Gamma-1}{\Gamma} \left(\frac{Q}{W^2}-\frac{D}{W}\right)=0.$$

Only evaluate it where $p\geq0$ and $Q\gt \sqrt{S^2}$; those inequalities define an admissible search region, not by themselves a root bracket. Establish a sign-changing bracket before a safeguarded solve. Once $p$ is known, $v_i=S_i/Q$ and $\epsilon=p/[(\Gamma-1)\rho]$. For the forward-test state $(D,S_x,\tau)=(1.25,1.40625,0.89375)$ and $\Gamma=5/3$, inserting $p=0.2$ gives $W=1.25$, $\rho=1$, $v^x=0.6$, and $f=0$.

In ideal MHD, with Eulerian $B^i$ and $\sqrt{4\pi}$ absorbed, the forward map instead contains

$$\begin{aligned} S_i&=(\rho hW^2+B^2)v_i-(B_jv^j)B_i,\\ \tau&=\rho hW^2+B^2-p-D -\frac12\left[\frac{B^2}{W^2}+(B_iv^i)^2\right]. \end{aligned}$$

Here $B^2=\gamma_{ij}B^iB^j$. Setting $B^i=0$ recovers hydrodynamics; setting $v^i=0$ gives $\tau=\rho\epsilon+B^2/2$. The latter limit exposes a fundamental loss of information. With $\rho=1$, $\epsilon=10^{-10}$, and $B^2/2=5000$, storing $\tau$ in binary64 and subtracting $5000$ gives approximately $1.00044417\times10^{-10}$: a relative thermal error $4.4\times10^{-4}$ before any nonlinear solve. A smaller solver tolerance cannot restore discarded digits. Kastaun et al. 2021, §§II.3 and V.2–V.3 separates recovery accuracy from this amplification of errors in the evolved state.

Physical admissibility before iteration

Require $D\gt 0$, $\rho\gt 0$, $p\geq0$, $W\geq1$, and subluminal velocity, plus EOS-domain and causality restrictions. Magnetized conserved variables obey additional inequalities. Screen obviously invalid inputs separately from solver nonconvergence. Silently flooring an invalid state and counting it as a successful recovery hides a failure of the evolution update.

Use analytic or rigorous brackets where possible. Newton iteration is fast near a root but can leave the physical domain; safeguarded Newton or Brent steps preserve a bracket. Scale residuals relative to the conserved magnitudes, and check both nonlinear residual and primitive change. Kastaun et al. 2021, §II including §II.H and §§III–V, Figures 3–5 develops robust recovery with explicit admissibility and error controls.

After convergence, reconstruct conserved variables from the recovered primitives using an independently coded forward map. Define a normalized round-trip error for every component. A small root residual can coexist with a forward mismatch when EOS derivatives, metric lowering, or magnetic contractions are inconsistent.

Stress-test the intended domain

Construct a reproducible sample spanning log density and temperature, Lorentz factor, magnetization, field–velocity angle, metric condition, and EOS regions. Include boundary cases deliberately but report the sampling measure; a uniform grid in one parametrization can nearly omit the difficult regime. Separate valid samples from injected invalid conserved states.

For each method report success fraction, iterations, forward residual, primitive error against known inputs, fallback rate, repair magnitude, and timing. Refine nonlinear tolerance and working precision. Compare at least two independent algorithms, such as a bracketed one-dimensional method and a multidimensional Newton solve. The benchmark suite should include shock-tube states, TOV surfaces, magnetized torus cells, and atmosphere interfaces.

During an evolution, output where and when recovery fails. Spatial clustering often identifies a surface, refinement boundary, or unphysical update. Global failure counts alone hide repeated repairs in a scientifically important ejecta region. The validity claim is restricted to the sampled EOS, magnetization, Lorentz factor, and metric range.

Specified round-trip benchmark

Use flat spacetime and a $\Gamma=4/3$ ideal-gas EOS, $p=(\Gamma-1)\rho\epsilon$ and $h=1+\epsilon+p/\rho$. Let

$$\rho\in[10^{-12},1], \qquad \epsilon\in[10^{-10},10^2], \qquad W\in\{1,2,10,100\}.$$

At sample levels $n=5,9,17$, use $n$ equally spaced points in each of $\log_{10}\rho$ and $\log_{10}\epsilon$. Cross them with magnetizations $\sigma=b^2/(\rho h)\in\{0,10^{-4},1,10^4\}$ and field–velocity cosines $\cos\theta\in\{-1,0,1\}$. Here $b^\mu$ is the comoving magnetic four-vector, $b^2=b^\mu b_\mu$, and $\sqrt{4\pi}$ is absorbed into the magnetic field. Choose $v^i=\sqrt{1-W^{-2}}\,\delta^i{}_x$ and an Eulerian $B^i$ in the $x$–$y$ plane; set its magnitude from $b^2=B^2/W^2+(B_iv^i)^2$. The parameter box is the complete algebraic domain, so there is no spatial boundary condition. Its EOS limits and the physical conditions $\rho\gt 0$, $\epsilon\gt 0$, and $W\geq1$ are the state-space boundaries.

Generate each conserved state with an independently coded binary128 forward map, retain the generating primitives as the reference, and recover with both binary64 and binary128 arithmetic at nonlinear tolerances $10^{-8},10^{-10},10^{-12}$. For every valid state add one negative control by replacing $D$ with $-|D|$ while retaining the other conserved components; it must be classified invalid before iteration.

At $W=1$ there is no velocity direction, and at $B=0$ there is no magnetic direction. Treat the angular labels as redundant in those cases, use the $x$ axis only as an arbitrary reference, and count each distinct state once in success fractions. Density is reported relative to one fixed reference density so its logarithmic sampling coordinate is dimensionless.

Define the primitive error by

$$E_\infty=\max( \left|\ln\frac{\rho_{\rm rec}}{\rho}\right|, \left|\ln\frac{\epsilon_{\rm rec}}{\epsilon}\right|, \frac{\|v^i_{\rm rec}-v^i\|_2}{\max(1,\|v^i\|_2)} )$$

and let $R_\infty$ be the maximum conserved-variable round-trip difference, componentwise normalized by $\max(|U_A|,D)$. Preserve two references: the original high-precision primitives, and a high-precision recovery of the same rounded input supplied to the binary64 solver. The first measures input-rounding plus recovery error; the second isolates recovery error. Screen the rounded state again for physical admissibility because rounding can move a boundary state outside the physical domain.

For the $n=17$, tolerance-$10^{-12}$ study, require all deliberately negative-$D$ controls to be rejected and every repair to be reported. Require $R_\infty\leq10^{-12}$ for unrepaired recoveries. Apply the $E_\infty\leq10^{-9}$ primitive target only to the subset for which perturbing each input by its rounding uncertainty, then recovering at high precision, changes the primitives by at most $10^{-10}$ in the same norm. On the remaining cold, highly boosted, or magnetically dominated states, report the attainable componentwise accuracy and compare it with that conditioning estimate; explicitly flag a requested accuracy that the input cannot support. Escalating recovery arithmetic alone does not repair an already rounded conservative state. Tightening the sample grid, input precision, and nonlinear tolerance independently must distinguish new failure classes from these information limits.

Observed discrepancy: not available—specification only.

Environment and failure criteria

Pin variable definitions, EOS implementation and version, metric convention, algorithm and brackets, tolerances, derivatives, precision, fallback order, floors and repair policy, sample seed and distribution, and code revision. A clean interface could be primitive-test --config extreme-states.toml --output run/, emitting inputs, recovered states, round-trip residuals, iteration traces, classifications, failure maps, sweeps, and environment metadata. This specifies a rerun record rather than claiming software is supplied.

Reject a method if it returns superluminal or EOS-invalid states, labels repaired inputs as converged, lacks tolerance convergence, or fails a forward round trip. Other failure modes include cancellation in thermal energy, an EOS derivative on the wrong branch, loss of a root bracket, magnetic contractions formed with the wrong metric, and a fallback that changes total mass or energy without recording it.

Diagnostic

A Newton solve reaches a residual of $10^{-12}$ and returns positive pressure, but reconstructing conserved energy differs from the input by $10^{-3}$. A bracketed fallback gives a larger root residual, $10^{-9}$, but a $10^{-10}$ round-trip error. Which solution is supported?

Solution

The bracketed result is supported. The Newton residual pertains to its implemented scalar equation, which may use an inconsistent EOS derivative or contraction; the independent forward map rejects it. Verify residual scaling and the forward implementation, tighten the fallback tolerance if needed, and require physicality plus round-trip agreement, not the smallest internal residual alone.

References

• Kastaun, W., Kalinani, J. V., and Ciolfi, R. (2021). “Robust recovery of primitive variables in relativistic ideal magnetohydrodynamics.” Physical Review D 103, 023018. arXiv.
• Noble, S. C., Gammie, C. F., McKinney, J. C., and Del Zanna, L. (2006). “Primitive variable solvers for conservative general relativistic magnetohydrodynamics.” Astrophysical Journal 641, 626–637. arXiv.