Hyperbolic Solver Interface

AbstractMesh{Dim}

Abstract supertype for all mesh types in the hyperbolic solver framework. Dim is the spatial dimension (1 or 2).

StructuredMesh1D{FT} <: AbstractMesh{1}

A uniform 1D structured mesh on the interval [xmin, xmax] with ncells cells.

Fields

• xmin::FT: Left boundary coordinate.
• xmax::FT: Right boundary coordinate.
• num_cells::Int: Number of cells.
• dx::FT: Cell width (xmax - xmin) / num_cells.
• cell_centers::Vector{FT}: Coordinates of cell centers.

StructuredMesh2D{FT} <: AbstractMesh{2}

A uniform 2D structured Cartesian mesh on [xmin, xmax] x [ymin, ymax].

Cells are numbered in column-major order: cell (i, j) has global index (j-1)*nx + i. Faces are split into x-faces (vertical, normal in x-direction) and y-faces (horizontal, normal in y-direction).

Fields

• xmin::FT, xmax::FT, ymin::FT, ymax::FT: Domain extents.
• nx::Int, ny::Int: Number of cells in each direction.
• dx::FT, dy::FT: Cell sizes.

StructuredMesh3D{FT} <: AbstractMesh{3}

A uniform 3D structured Cartesian mesh on [xmin, xmax] x [ymin, ymax] x [zmin, zmax].

Cells are numbered in column-major order: cell (i, j, k) has global index (k-1)*nx*ny + (j-1)*nx + i.

Fields

• xmin::FT, xmax::FT, ymin::FT, ymax::FT, zmin::FT, zmax::FT: Domain extents.
• nx::Int, ny::Int, nz::Int: Number of cells in each direction.
• dx::FT, dy::FT, dz::FT: Cell sizes.

ndims_mesh(mesh::AbstractMesh{Dim}) -> Int

Return the spatial dimension of the mesh.

ncells(mesh::AbstractMesh)

Return the number of cells in the mesh.

nfaces(mesh::AbstractMesh)

Return the number of internal faces in the mesh.

cell_center(mesh::AbstractMesh, i::Int)

Return the coordinates of the center of cell i.

cell_volume(mesh::AbstractMesh, i::Int)

Return the volume (or area in 2D, length in 1D) of cell i.

face_area(mesh::AbstractMesh, f::Int)

Return the area (or length in 2D, 1.0 in 1D) of face f.

face_owner(mesh::AbstractMesh, f::Int)

Return the index of the cell that owns face f (the "left" cell).

face_neighbor(mesh::AbstractMesh, f::Int)

Return the index of the neighbor cell across face f (the "right" cell).

AbstractEOS

Abstract supertype for equations of state.

All subtypes must implement:

• pressure(eos, ρ, ε): Compute pressure from density and specific internal energy.
• sound_speed(eos, ρ, P): Compute sound speed from density and pressure.
• internal_energy(eos, ρ, P): Compute specific internal energy from density and pressure.

IdealGasEOS{FT} <: AbstractEOS

Ideal (gamma-law) equation of state: P = (γ - 1) ρ ε, where ε is the specific internal energy. Equivalent to P V = n R T for a calorically perfect gas with constant ratio of specific heats γ = cₚ/cᵥ.

Sound speed: c = √(γ P / ρ).

Fields

• gamma::FT: Adiabatic index. Common values: 5/3 (monatomic), 7/5 (diatomic), 4/3 (relativistic gas or radiation-dominated).

StiffenedGasEOS{FT} <: AbstractEOS

Stiffened gas equation of state: P = (γ - 1) ρ ε - γ P∞.

This generalizes the ideal gas EOS with a stiffening pressure P∞ that models the repulsive molecular interactions in dense fluids (e.g., water). Setting P∞ = 0 recovers the ideal gas law.

Fields

• gamma::FT: Adiabatic index.
• P_inf::FT: Stiffening pressure. Typical values: 0 (ideal gas), ~6.08e8 Pa (water).

pressure(eos::AbstractEOS, ρ, ε)

Compute the thermodynamic pressure given density ρ and specific internal energy ε.

sound_speed(eos::AbstractEOS, ρ, P)

Compute the adiabatic sound speed given density ρ and pressure P.

internal_energy(eos::AbstractEOS, ρ, P)

Compute the specific internal energy given density ρ and pressure P.

total_energy(eos::IdealGasEOS, ρ, v, P) -> E

Compute the total energy density E = P/(γ-1) + ½ρv² for 1D.

total_energy(eos::IdealGasEOS, ρ, vx, vy, P) -> E

Compute the total energy density E = P/(γ-1) + ½ρ(vx² + vy²) for 2D.

total_energy(eos::StiffenedGasEOS, ρ, v, P) -> E

Compute the total energy density for 1D: E = (P + γP∞)/(γ-1) + ½ρv²

total_energy(eos::StiffenedGasEOS, ρ, vx, vy, P) -> E

Compute the total energy density for 2D: E = (P + γP∞)/(γ-1) + ½ρ(vx² + vy²)

total_energy(eos::IdealGasEOS, ρ, vx, vy, vz, P) -> E

Compute the total energy density E = P/(γ-1) + ½ρ(vx² + vy² + vz²) for 3D.

total_energy(eos::StiffenedGasEOS, ρ, vx, vy, vz, P) -> E

Compute the total energy density for 3D stiffened gas: E = (P + γP∞)/(γ-1) + ½ρ(vx² + vy² + vz²)

AbstractConservationLaw{Dim}

Abstract supertype for hyperbolic conservation laws in Dim spatial dimensions.

The general form is: ∂U/∂t + ∇·F(U) = S(U) where U is the vector of conserved variables, F(U) is the flux tensor, and S(U) is the source term.

All subtypes must implement:

• nvariables(law): Number of conserved variables.
• physical_flux(law, u, dir): Compute the physical flux F(U) in direction dir (1=x, 2=y, 3=z).
• max_wave_speed(law, u, dir): Maximum wave speed in direction dir.
• conserved_to_primitive(law, u): Convert conserved variables to primitive variables.
• primitive_to_conserved(law, p): Convert primitive variables to conserved variables.

EulerEquations{Dim, EOS <: AbstractEOS} <: AbstractConservationLaw{Dim}

The compressible Euler equations in Dim spatial dimensions — inviscid conservation of mass, momentum, and total energy for a compressible fluid.

1D (Dim=1)

Conserved variables: U = [ρ, ρv, E] Primitive variables: W = [ρ, v, P] Flux: F = [ρv, ρv² + P, (E + P)v]

2D (Dim=2)

Conserved variables: U = [ρ, ρvx, ρvy, E] Primitive variables: W = [ρ, vx, vy, P] Fluxes: Fx = [ρvx, ρvx² + P, ρvx*vy, (E + P)vx] Fy = [ρvy, ρvx*vy, ρvy² + P, (E + P)vy]

Closed by an equation of state P = P(ρ, ε) (see AbstractEOS).

Fields

• eos::EOS: Equation of state.

References

• Toro, E. F. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed., Springer.
• LeVeque, R. J. (2002). Finite Volume Methods for Hyperbolic Problems, Cambridge University Press.

IdealMHDEquations{Dim, EOS <: AbstractEOS} <: AbstractConservationLaw{Dim}

The ideal magnetohydrodynamics (MHD) equations in Dim spatial dimensions — conservation of mass, momentum, energy, and magnetic flux for an electrically conducting inviscid fluid in the infinite-conductivity (ideal) limit.

Variables (8 components in all dimensions)

• Primitive: W = [ρ, vx, vy, vz, P, Bx, By, Bz]
• Conserved: U = [ρ, ρvx, ρvy, ρvz, E, Bx, By, Bz]

Total energy: E = P/(γ-1) + ½ρ|v|² + ½|B|²

The system must satisfy the solenoidal constraint ∇·B = 0, which is maintained to machine precision via constrained transport (see CTData2D).

1D

Bx is constant (its x-flux is zero). The 1D system still has 8 variables but the Bx equation is trivially satisfied.

2D

The ∇·B = 0 constraint is maintained via a staggered (Yee-type) update of face-centered magnetic fields using the corner EMF.

Fields

• eos::EOS: Equation of state.

References

• Powell, K. G. et al. (1999). "A solution-adaptive upwind scheme for ideal magnetohydrodynamics." J. Comput. Phys., 154(2), 284–309.
• Tóth, G. (2000). "The ∇·B = 0 constraint in shock-capturing magnetohydrodynamics codes." J. Comput. Phys., 161(2), 605–652.

NavierStokesEquations{Dim, EOS <: AbstractEOS} <: AbstractConservationLaw{Dim}

The compressible Navier-Stokes equations in Dim spatial dimensions.

Wraps EulerEquations and adds viscous transport coefficients. All inviscid interface methods (physicalflux, wavespeeds, etc.) are delegated to the underlying Euler equations.

Fields

• euler::EulerEquations{Dim, EOS}: The inviscid part.
• mu::Float64: Dynamic viscosity μ.
• Pr::Float64: Prandtl number.

nvariables(law::AbstractConservationLaw) -> Int

Return the number of conserved variables for the conservation law.

physical_flux(law::AbstractConservationLaw, u::SVector, dir::Int) -> SVector

Compute the physical flux vector F(U) in direction dir (1=x, 2=y, 3=z). The input u should be a vector of primitive variables.

max_wave_speed(law::AbstractConservationLaw, u::SVector, dir::Int)

Return the maximum absolute wave speed in direction dir. The input u should be a vector of primitive variables.

wave_speeds(law::EulerEquations{1}, w::SVector{3}, ::Int) -> (λ_min, λ_max)

Return the minimum and maximum wave speeds from primitive variables.

wave_speeds(law::EulerEquations{2}, w::SVector{4}, dir::Int) -> (λ_min, λ_max)

Return the minimum and maximum wave speeds in direction dir from primitive variables.

wave_speeds(law::IdealMHDEquations, w::SVector{8}, dir::Int) -> (λ_min, λ_max)

Return the fastest left-going and right-going wave speeds (fast magnetosonic).

wave_speeds(law::ShallowWaterEquations{1}, w::SVector{2}, ::Int) -> (lambda_min, lambda_max)

Return the minimum and maximum wave speeds from primitive variables.

wave_speeds(law::ShallowWaterEquations{2}, w::SVector{3}, dir::Int) -> (lambda_min, lambda_max)

Return the minimum and maximum wave speeds in direction dir from primitive variables.

wave_speeds(law::SRHydroEquations{1}, w::SVector{3}, dir::Int) -> (λ_min, λ_max)

Return the fastest left-going and right-going wave speeds for 1D SR hydro.

wave_speeds(law::SRHydroEquations{2}, w::SVector{4}, dir::Int) -> (λ_min, λ_max)

Return the fastest left-going and right-going wave speeds for 2D SR hydro.

wave_speeds(law::TwoFluidEquations{1}, w::SVector{6}, ::Int) -> (λ_min, λ_max)

Return the envelope of wave speeds across both species: (min(v_i - c_i, v_e - c_e), max(v_i + c_i, v_e + c_e)).

wave_speeds(law::TwoFluidEquations{2}, w::SVector{8}, dir::Int) -> (λ_min, λ_max)

Return the envelope of wave speeds in direction dir across both species.

wave_speeds(law::SRMHDEquations, w::SVector{8}, dir::Int) -> (λ_min, λ_max)

Return the fastest left-going and right-going wave speeds.

wave_speeds(law::GRMHDEquations, w::SVector{8}, dir::Int) -> (lam_min, lam_max)

Return the fastest left-going and right-going wave speeds in flat-space form.

wave_speeds(law::EulerEquations{3}, w::SVector{5}, dir::Int) -> (λ_min, λ_max)

Return the minimum and maximum wave speeds in direction dir from primitive variables.

wave_speeds(law::IdealMHDEquations{3}, w::SVector{8}, dir::Int) -> (λ_min, λ_max)

Return the fastest left-going and right-going wave speeds for 3D MHD.

wave_speeds(law::ReactiveEulerEquations{1}, w, ::Int) -> (λ_min, λ_max)

Wave speed bounds — delegate to Euler.

wave_speeds(law::ReactiveEulerEquations{2}, w, dir::Int) -> (λ_min, λ_max)

Wave speed bounds — delegate to Euler.

conserved_to_primitive(law::AbstractConservationLaw, u::SVector) -> SVector

Convert a vector of conserved variables to primitive variables.

primitive_to_conserved(law::AbstractConservationLaw, w::SVector) -> SVector

Convert a vector of primitive variables to conserved variables.

fast_magnetosonic_speed(law::IdealMHDEquations, w::SVector{8}, dir::Int) -> cf

Compute the fast magnetosonic speed in direction dir: cf² = ½(a² + b² + √((a² + b²)² - 4 a² bn²)) where a² = γP/ρ, b² = |B|²/ρ, bn = Bn/√ρ.

fast_magnetosonic_speed(law::IdealMHDEquations{3}, w::SVector{8}, dir::Int) -> cf

Compute the fast magnetosonic speed in direction dir (1=x, 2=y, 3=z).

AbstractRiemannSolver

Abstract supertype for approximate Riemann solvers used at cell interfaces.

All subtypes must support calling via solve_riemann(solver, law, wL, wR, dir).

LaxFriedrichsSolver <: AbstractRiemannSolver

The Lax-Friedrichs (Rusanov) approximate Riemann solver.

The numerical flux is: F* = ½(F(UL) + F(UR)) - ½ λ_max (UR - UL) where λ_max = max(|vL| + cL, |vR| + cR).

This is the most diffusive but most robust Riemann solver. First-order accurate at discontinuities; suitable as a fallback or baseline solver.

References

• Rusanov, V. V. (1962). "The calculation of the interaction of non-stationary shock waves with barriers." Zh. Vychisl. Mat. Mat. Fiz., 1(2), 267–279.
• Lax, P. D. (1954). "Weak solutions of nonlinear hyperbolic equations and their numerical computation." Comm. Pure Appl. Math., 7(1), 159–193.

HLLSolver <: AbstractRiemannSolver

The Harten-Lax-van Leer (HLL) approximate Riemann solver.

Uses estimates of the fastest left-going and right-going wave speeds SL, SR to compute a two-wave flux:

If SL ≥ 0: F* = FL If SR ≤ 0: F* = FR Otherwise: F* = (SR*FL - SL*FR + SL*SR*(UR - UL)) / (SR - SL)

Resolves rarefactions and shocks but smears contact discontinuities (see HLLC). Wave speed estimates use the Davis bounds.

References

• Harten, A., Lax, P. D., & van Leer, B. (1983). "On upstream differencing and Godunov-type schemes for hyperbolic conservation laws." SIAM Review, 25(1), 35–61.
• Davis, S. F. (1988). "Simplified second-order Godunov-type methods." SIAM J. Sci. Stat. Comput., 9(3), 445–473.

HLLCSolver <: AbstractRiemannSolver

The HLLC (Harten-Lax-van Leer-Contact) approximate Riemann solver.

A 3-wave model that resolves:

• Left-going wave (speed SL)
• Contact discontinuity (speed S*)
• Right-going wave (speed SR)

Exactly resolves isolated contact discontinuities and shear waves, which the 2-wave HLL solver smears. The standard choice for the Euler equations.

References

• Toro, E. F., Spruce, M., & Speares, W. (1994). "Restoration of the contact surface in the HLL-Riemann solver." Shock Waves, 4(1), 25–34.
• Toro, E. F. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed., Springer. Chapter 10.

HLLDSolver <: AbstractRiemannSolver

The HLLD (Harten-Lax-van Leer-Discontinuities) approximate Riemann solver for ideal MHD equations.

Resolves 5 waves: fast magnetosonic shocks (SL, SR), Alfvén/rotational discontinuities (SL, SR), and the contact/entropy wave (SM). This is the standard Riemann solver for ideal MHD, analogous to HLLC for hydrodynamics.

References

• Miyoshi, T. & Kusano, K. (2005). "A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics." J. Comput. Phys., 208(1), 315–344.

solve_riemann(solver::AbstractRiemannSolver, law, wL, wR, dir) -> SVector

Compute the numerical flux at a cell interface given left and right primitive states wL, wR, and direction dir (1=x, 2=y, 3=z).

Returns the numerical flux vector.

NoReconstruction

First-order (piecewise constant) scheme with no reconstruction.

CellCenteredMUSCL{L <: AbstractLimiter}

MUSCL (Monotonic Upstream-centered Scheme for Conservation Laws) reconstruction for cell-centered data on structured meshes.

Constructs a piecewise-linear approximation in each cell:

\[u(x) = \bar{u}_i + \sigma_i \, (x - x_i)\]

where $\bar{u}_i$ is the cell average, $x_i$ the cell centre, and $\sigma_i$ is a limited slope that enforces TVD (or TVB) properties. The limiter bounds the ratio of consecutive slopes to prevent spurious oscillations near discontinuities while retaining second-order accuracy in smooth regions.

Fields

• limiter::L: Slope limiter to use for reconstruction.

References

• van Leer, B. (1979). "Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method." J. Comput. Phys., 32(1), 101–136.

reconstruct_interface(recon::CellCenteredMUSCL, uLL, uL, uR, uRR) -> (uL_face, uR_face)

Reconstruct interface states at the face between cell L and cell R.

• uLL: value two cells to the left
• uL: value one cell to the left (owner)
• uR: value one cell to the right (neighbor)
• uRR: value two cells to the right

Returns (uL_face, uR_face) — the left and right states at the interface.

reconstruct_interface(recon::PPMReconstruction, wLL, wL, wR, wRR) -> (wL_face, wR_face)

Reconstruct interface states at the face between cell L and cell R using PPM reconstruction on primitive variables.

• wLL: value two cells to the left
• wL: value one cell to the left (owner)
• wR: value one cell to the right (neighbor)
• wRR: value two cells to the right

Returns (wL_face, wR_face) — the left and right states at the interface.

reconstruct_interface(recon::WENO3, wLL, wL, wR, wRR) -> (wL_face, wR_face)

Reconstruct interface states at the face between cell L and cell R using WENO-3 reconstruction on primitive variables.

• wLL: value two cells to the left
• wL: value one cell to the left (owner)
• wR: value one cell to the right (neighbor)
• wRR: value two cells to the right

This follows the same interface as CellCenteredMUSCL.

reconstruct_interface(cw::CharacteristicWENO{<:WENO3}, wLL, wL, wR, wRR) -> (wL_face, wR_face)

Component-wise fallback when called without law information.

The actual characteristic projection is performed in the 1D/2D dispatch functions (reconstruct_interface_1d, _reconstruct_face_2d, etc.) where the conservation law is available.

WENO3{FT}

Weighted Essentially Non-Oscillatory reconstruction of order 3.

Uses a 3-point stencil with 2 candidate linear polynomials and nonlinear weights that adapt to solution smoothness. Compatible with the existing 2-ghost-cell padding.

Fields

• epsilon::FT: Small parameter to avoid division by zero in weights (default 1e-6).

WENO5{FT}

Fifth-order Weighted Essentially Non-Oscillatory (WENO) reconstruction with WENO-Z nonlinear weights.

Constructs a convex combination of three candidate quadratic polynomials on a 5-point stencil {i-2, …, i+2}. The nonlinear weights adapt so that stencils crossing discontinuities are assigned near-zero weight, achieving fifth-order accuracy in smooth regions while remaining essentially non-oscillatory at shocks.

Requires 3 ghost cells on each side of the domain.

Fields

• epsilon::FT: Small parameter to avoid division by zero (default 1e-6).
• p::Int: Exponent for weight computation (default 2).

References

• Jiang, G.-S. & Shu, C.-W. (1996). "Efficient implementation of weighted ENO schemes." J. Comput. Phys., 126(1), 202–228.
• Borges, R. et al. (2008). "An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws." J. Comput. Phys., 227(6), 3191–3211. (WENO-Z weights)

CharacteristicWENO{R}

Wrapper that applies characteristic decomposition before WENO reconstruction and projects back afterwards.

The inner reconstruction recon must be a WENO-type scheme (e.g., WENO3 or WENO5).

Fields

• recon::R: Underlying WENO reconstruction scheme.

nghost(recon) -> Int

Return the number of ghost cells required on each side of the domain for the given reconstruction scheme.

left_eigenvectors(law, w, dir) -> L

Return the matrix of left eigenvectors (rows are left eigenvectors) of the flux Jacobian ∂F/∂U in direction dir, evaluated at the primitive state w.

The left eigenvector matrix satisfies L * R = I where R is the right eigenvector matrix.

right_eigenvectors(law, w, dir) -> R

Return the matrix of right eigenvectors (columns are right eigenvectors) of the flux Jacobian ∂F/∂U in direction dir, evaluated at the primitive state w.

AbstractHyperbolicBC

Abstract supertype for ghost-cell boundary conditions in the hyperbolic solver.

TransmissiveBC <: AbstractHyperbolicBC

Zero-gradient (outflow) boundary condition. Ghost cell values are copied from the nearest interior cells (extrapolation of order 0).

ReflectiveBC <: AbstractHyperbolicBC

Reflective (slip wall) boundary condition. The normal velocity component is negated in the ghost cells while density, pressure, and tangential velocities are copied.

NoSlipBC <: AbstractHyperbolicBC

No-slip wall boundary condition. ALL velocity components are negated in the ghost cells, placing zero velocity at the wall face midpoint.

Unlike ReflectiveBC (slip wall, which only negates the wall-normal velocity), NoSlipBC enforces zero tangential velocity as well, appropriate for viscous (Navier-Stokes) flows.

DirichletHyperbolicBC{N, FT} <: AbstractHyperbolicBC

Fixed-state boundary condition for hyperbolic problems. The ghost cells are set to maintain the prescribed primitive state at the boundary.

Fields

• state::SVector{N, FT}: Prescribed primitive state.

InflowBC{N, FT} <: AbstractHyperbolicBC

Prescribes all primitive variables at the boundary.

Fields

• state::SVector{N, FT}: Prescribed primitive state [ρ, v, P] (1D) or [ρ, vx, vy, P] (2D).

PeriodicHyperbolicBC <: AbstractHyperbolicBC

Periodic boundary condition: the left ghost cells are filled from the right interior cells and vice versa.

HyperbolicProblem{Law, Mesh, RS, Rec, BCL, BCR, IC, FT}

A cell-centered finite volume problem for hyperbolic conservation laws.

Fields

• law: The conservation law (e.g., EulerEquations).
• mesh: The computational mesh.
• riemann_solver: The approximate Riemann solver.
• reconstruction: Reconstruction scheme (e.g., CellCenteredMUSCL).
• bc_left: Left boundary condition.
• bc_right: Right boundary condition.
• initial_condition: Function (x) -> SVector{N} of primitive variables, or a vector of conserved SVectors.
• initial_time::FT: Start time.
• final_time::FT: End time.
• cfl::FT: CFL number for time step control.

HyperbolicProblem2D{Law, Mesh, RS, Rec, BC_L, BC_R, BC_B, BC_T, IC, FT}

A cell-centered finite volume problem for 2D hyperbolic conservation laws on structured Cartesian meshes.

Fields

• law: The conservation law (e.g., EulerEquations{2}).
• mesh: The 2D computational mesh (StructuredMesh2D).
• riemann_solver: The approximate Riemann solver.
• reconstruction: Reconstruction scheme.
• bc_left, bc_right, bc_bottom, bc_top: Boundary conditions on each side.
• initial_condition: Function (x, y) -> SVector{N} of primitive variables.
• initial_time::FT, final_time::FT: Time span.
• cfl::FT: CFL number for time step control.

HyperbolicProblem3D{Law, Mesh, RS, Rec, BC_L, BC_R, BC_Bo, BC_T, BC_F, BC_Ba, IC, FT}

A cell-centered finite volume problem for 3D hyperbolic conservation laws on structured Cartesian meshes.

Fields

• law: The conservation law (e.g., EulerEquations{3}).
• mesh: The 3D computational mesh (StructuredMesh3D).
• riemann_solver: The approximate Riemann solver.
• reconstruction: Reconstruction scheme.
• bc_left, bc_right: Boundary conditions in x-direction.
• bc_bottom, bc_top: Boundary conditions in y-direction.
• bc_front, bc_back: Boundary conditions in z-direction.
• initial_condition: Function (x, y, z) -> SVector{N} of primitive variables.
• initial_time::FT, final_time::FT: Time span.
• cfl::FT: CFL number for time step control.

mhd_stage_limiter(cache) -> Function

Create a stage limiter function for MHD/CT constrained transport.

The returned function enforces cell-centered B = avg(face-centered B) after each RK stage. Pass it to the ODE algorithm:

ode_prob = ODEProblem(prob)
limiter = mhd_stage_limiter(ode_prob.p)
sol = solve(ode_prob, SSPRK33(; stage_limiter! = limiter); adaptive = false, dt = dt0)

compute_dt(prob::HyperbolicProblem, U::AbstractVector, t) -> FT

Compute the time step from the CFL condition: Δt = cfl * Δx / max(|λ|)

compute_dt(prob::HyperbolicProblem{<:NavierStokesEquations{1}}, U, t) -> dt

Compute the time step accounting for both hyperbolic and viscous stability:

• dt_hyp = CFL · dx / max|λ|
• dt_visc = 0.5 · ρ_min · dx² / μ
• dt = CFL · min(dt_hyp, dt_visc)

compute_dt_2d(prob::HyperbolicProblem2D, U::AbstractMatrix, t) -> FT

Compute the time step from the 2D CFL condition: Δt = cfl / (max(|λx|)/Δx + max(|λy|)/Δy)

compute_dt_2d(prob::HyperbolicProblem2D{<:NavierStokesEquations{2}}, U, t) -> dt

Compute the time step accounting for both hyperbolic and viscous stability:

• Hyperbolic: dt_hyp = cfl / (max(|λx|)/dx + max(|λy|)/dy)
• Viscous: dt_visc = 0.5 · ρ_min / (μ · (1/dx² + 1/dy²))

compute_dt_2d(prob::HyperbolicProblem2D{<:GRMHDEquations{2}}, U, t, md) -> dt

Compute the time step using metric-corrected wave speeds: lambda_coord = alpha * lambda_flat - beta

compute_dt_3d(prob::HyperbolicProblem3D, U::AbstractArray{T,3}, t) -> FT

Compute the time step from the 3D CFL condition: dt = cfl / (max(|lx|)/dx + max(|ly|)/dy + max(|lz|)/dz)

thermal_conductivity(ns::NavierStokesEquations) -> κ

Compute the thermal conductivity κ = μ γ / (Pr (γ - 1)).

viscous_flux_1d(ns::NavierStokesEquations{1}, wL::SVector{3}, wR::SVector{3}, dx) -> SVector{3}

Compute the viscous flux at a 1D face between left and right primitive states.

Stress tensor: τ_xx = (4/3) μ ∂v/∂x Heat flux: q_x = -κ ∂T/∂x where T = P/ρ (non-dimensional temperature)

Returns SVector(0, τ_xx, v_face * τ_xx - q_x).

viscous_flux_x_2d(ns::NavierStokesEquations{2}, wL::SVector{4}, wR::SVector{4},
                  dvx_dy, dvy_dy, dx) -> SVector{4}

Compute the viscous flux at a 2D x-direction face.

Normal gradients are computed from the direct neighbors; cross-derivatives ∂vx/∂y and ∂vy/∂y must be provided (4-cell average).

Returns SVector(0, τ_xx, τ_xy, vx*τ_xx + vy*τ_xy - q_x).

viscous_flux_y_2d(ns::NavierStokesEquations{2}, wB::SVector{4}, wT::SVector{4},
                  dvx_dx, dvy_dx, dy) -> SVector{4}

Compute the viscous flux at a 2D y-direction face.

Normal gradients are computed from the direct neighbors; cross-derivatives ∂vx/∂x and ∂vy/∂x must be provided (4-cell average).

Returns SVector(0, τ_yx, τ_yy, vx*τ_yx + vy*τ_yy - q_y).

CTData2D{FT}

Constrained transport data for 2D MHD on a structured mesh.

Stores face-centered magnetic field components and corner EMFs.

Fields

• Bx_face::Matrix{FT}: x-component of B at x-faces, size (nx+1) × ny.
• By_face::Matrix{FT}: y-component of B at y-faces, size nx × (ny+1).
• emf_z::Matrix{FT}: z-component of EMF at corners, size (nx+1) × (ny+1).

CTData3D{FT}

Constrained transport data for 3D MHD on a structured mesh.

Stores face-centered magnetic field components and edge-centered EMFs.

Fields

• Bx_face::Array{FT,3}: x-component of B at x-faces, size (nx+1) x ny x nz.
• By_face::Array{FT,3}: y-component of B at y-faces, size nx x (ny+1) x nz.
• Bz_face::Array{FT,3}: z-component of B at z-faces, size nx x ny x (nz+1).
• emf_x::Array{FT,3}: x-component of EMF at x-edges, size nx x (ny+1) x (nz+1).
• emf_y::Array{FT,3}: y-component of EMF at y-edges, size (nx+1) x ny x (nz+1).
• emf_z::Array{FT,3}: z-component of EMF at z-edges, size (nx+1) x (ny+1) x nz.

initialize_ct!(ct::CTData2D, prob::HyperbolicProblem2D, mesh)

Initialize face-centered B fields from the initial condition. Face values are set by evaluating the initial condition at the face center and extracting the appropriate B component.

initialize_ct_from_potential!(ct::CTData2D, Az_func, mesh)

Initialize face-centered B fields from a vector potential function Az(x, y). This guarantees ∇·B = 0 to machine precision via Stokes' theorem:

Bx_face[i,j] = (Az(x, y_top) - Az(x, y_bottom)) / dy
By_face[i,j] = -(Az(x_right, y) - Az(x_left, y)) / dx

face_to_cell_B!(U::AbstractMatrix, ct::CTData2D, nx::Int, ny::Int)

Update cell-centered B in the conserved variable array U from face-centered values. Cell-centered B is the arithmetic mean of the two adjacent face values.

Interior cell (ix, iy) maps to U[ix+2, iy+2] in the padded array.

compute_emf_2d!(ct::CTData2D, Fx_faces::AbstractMatrix, Fy_faces::AbstractMatrix,
                nx::Int, ny::Int)

Compute the corner EMF emf_z from the x-direction and y-direction numerical fluxes.

Arguments

• Fx_faces: x-fluxes at x-faces, size (nx+1) × ny. Fx_faces[i,j] is the flux at the face between cells (i-1,j) and (i,j).
• Fy_faces: y-fluxes at y-faces, size nx × (ny+1). Fy_faces[i,j] is the flux at the face between cells (i,j-1) and (i,j).

The EMF at x-faces is Ez_x = -Fx[7] (negative of the By induction flux). The EMF at y-faces is Ez_y = +Fy[6] (positive of the Bx induction flux).

Corner values are the arithmetic average of the 4 surrounding face EMFs.

ct_update!(ct::CTData2D, dt, dx, dy, nx::Int, ny::Int)

Update face-centered B-fields using the EMF at cell corners. This guarantees ∇·B = 0 to machine precision.

Bx_face[i,j] -= dt/dy * (emf_z[i,j+1] - emf_z[i,j]) By_face[i,j] += dt/dx * (emf_z[i+1,j] - emf_z[i,j])

compute_divB(ct::CTData2D, dx, dy, nx::Int, ny::Int) -> Matrix{FT}

Compute the cell-centered ∇·B from face-centered B values: ∇·B[i,j] = (Bx_face[i+1,j] - Bx_face[i,j]) / dx + (By_face[i,j+1] - By_face[i,j]) / dy

Should be zero to machine precision when CT is used.

max_divB(ct::CTData2D, dx, dy, nx::Int, ny::Int) -> FT

Compute the maximum |∇·B| over all cells.

max_divB(ct::CTData2D, mesh::StructuredMesh2D) -> FT

Convenience method that extracts mesh parameters automatically.

l2_divB(ct::CTData2D, dx, dy, nx::Int, ny::Int) -> FT

Compute the L2 norm of ∇·B.

AbstractIMEXScheme

Abstract supertype for IMEX Runge-Kutta time integration schemes.

All subtypes should have a corresponding imex_tableau method returning the Butcher tableaux.

IMEX_SSP3_433 <: AbstractIMEXScheme

Third-order IMEX scheme with 4 explicit stages and 3 implicit stages, based on the SSP3(4,3,3) method of Pareschi & Russo (2005).

This is a good default for moderately stiff problems. The explicit part is a 4-stage, 3rd-order SSP scheme. The implicit part is a 3rd-order L-stable DIRK scheme.

IMEX_ARS222 <: AbstractIMEXScheme

Second-order IMEX scheme with 2 implicit stages and 2 explicit stages, from Ascher, Ruuth & Spiteri (1997).

A simple, robust choice for mildly stiff problems.

IMEX_Midpoint <: AbstractIMEXScheme

First-order IMEX method based on the implicit-explicit midpoint rule. Primarily useful for testing the IMEX infrastructure.

2 stages:

• Stage 1: explicit evaluation at t^n
• Stage 2: implicit solve at t^{n+1}

imex_tableau(scheme::AbstractIMEXScheme) -> NamedTuple

Return the Butcher tableaux for the IMEX scheme as a named tuple: (A_ex, b_ex, c_ex, A_im, b_im, c_im, s)

where s is the number of stages, and:

• A_ex: Explicit RK matrix (s x s, strictly lower triangular).
• b_ex: Explicit weights (length s).
• c_ex: Explicit abscissae (length s).
• A_im: Implicit RK matrix (s x s, lower triangular with non-zero diagonal).
• b_im: Implicit weights (length s).
• c_im: Implicit abscissae (length s).

imex_nstages(scheme::AbstractIMEXScheme) -> Int

Return the number of stages for the IMEX scheme.

AbstractStiffSource

Abstract supertype for stiff source terms treated implicitly in IMEX time integration.

All subtypes must implement:

• evaluate_stiff_source(src, law, w, u) -> SVector
• stiff_source_jacobian(src, law, w, u) -> SMatrix

where w is the primitive state and u is the conserved state.

NullSource <: AbstractStiffSource

A no-op stiff source that returns zero. Useful for testing the IMEX framework in purely explicit mode.

CoolingSource{FT, F} <: AbstractStiffSource

Optically thin radiative cooling source term for the energy equation.

The cooling rate depends on local density and temperature: S_E = -ρ² Λ(T)

where Λ(T) is the cooling function.

Fields

• cooling_function::F: Function T -> Λ(T) returning the cooling rate coefficient.
• mu_mol::FT: Mean molecular weight (for temperature computation).

ResistiveSource{FT} <: AbstractStiffSource

Resistive (Ohmic) source term for ideal MHD equations.

In the resistive MHD framework, the magnetic field evolves as: ∂B/∂t = ∇×(v×B) + η ∇²B

The diffusive part η ∇²B is stiff when η is large. In a spatially discretised form, this gives a source: S = η / dx² * (B_{i-1} - 2B_i + B_{i+1})

For the IMEX framework, we treat the local contribution as: S_i = -2η / dx² * B_i (implicit part, stiff) and the neighbour contributions explicitly.

Fields

• eta::FT: Magnetic diffusivity (resistivity / μ₀).

evaluate_stiff_source(src::AbstractStiffSource, law, w, u) -> SVector

Evaluate the stiff source term S(U) at the given state.

Arguments

• src: The stiff source term.
• law: The conservation law.
• w: Primitive variable vector.
• u: Conserved variable vector.

Returns

Source vector in conserved variable space (same length as u).

stiff_source_jacobian(src::AbstractStiffSource, law, w, u) -> SMatrix

Compute the Jacobian ∂S/∂U of the stiff source term with respect to the conserved variables.

Arguments

• src: The stiff source term.
• law: The conservation law.
• w: Primitive variable vector.
• u: Conserved variable vector.

Returns

Jacobian matrix ∂S/∂U as an SMatrix.

AMRBlock{N, FT, Dim}

A single block in the AMR hierarchy.

Fields

• id::Int: Unique block identifier.
• level::Int: Refinement level (0 = coarsest).
• origin::NTuple{Dim, FT}: Physical coordinates of the block's lower-left corner.
• dims::NTuple{Dim, Int}: Number of cells in each direction (nx, ny) or (nx, ny, nz).
• dx::NTuple{Dim, FT}: Cell sizes in each direction.
• U::Array{SVector{N, FT}, Dim}: Conserved variable data (interior only, no ghosts).
• parent_id::Int: ID of the parent block (-1 if root).
• child_ids::Vector{Int}: IDs of child blocks (empty if leaf).
• active::Bool: Whether this block is a leaf (active for computation).
• neighbors::Dict{Symbol, Int}: Neighbor block IDs keyed by face (:left, :right, etc.).

AMRGrid{N, FT, Dim, Law, Criterion}

The AMR grid managing the block hierarchy.

Fields

• blocks::Dict{Int, AMRBlock{N, FT, Dim}}: All blocks keyed by ID.
• max_level::Int: Maximum allowed refinement level.
• refinement_ratio::Int: Refinement factor (always 2).
• block_size::NTuple{Dim, Int}: Number of cells per block in each direction.
• law::Law: The conservation law.
• criterion::Criterion: The refinement criterion.
• next_id::Base.RefValue{Int}: Counter for assigning block IDs.

AbstractRefinementCriterion

Abstract supertype for refinement criteria that decide when to refine or coarsen.

GradientRefinement{FT} <: AbstractRefinementCriterion

Refinement based on the maximum gradient of a specified variable. Refines when max(|grad(var)|) * dx > refine_threshold. Coarsens when max(|grad(var)|) * dx < coarsen_threshold.

Fields

• variable_index::Int: Index of the variable to monitor (in conserved vector).
• refine_threshold::FT: Threshold for refinement.
• coarsen_threshold::FT: Threshold for coarsening.

CurrentSheetRefinement{FT} <: AbstractRefinementCriterion

Refinement based on current sheet detection for MHD problems. Refines when |J| = |curl(B)| * dx > refine_threshold. Coarsens when |J| * dx < coarsen_threshold.

Designed for 2D/3D MHD where current sheets indicate magnetic reconnection or sharp field gradients.

Fields

• refine_threshold::FT: Threshold for refinement.
• coarsen_threshold::FT: Threshold for coarsening.

AMRProblem{Grid, RS, Rec, BCs, FT}

An AMR problem definition wrapping the grid, solver settings, and boundary conditions.

Fields

• grid::Grid: The AMR grid hierarchy.
• riemann_solver::RS: The Riemann solver.
• reconstruction::Rec: The reconstruction scheme.
• boundary_conditions::BCs: Boundary conditions (applied at domain boundaries).
• initial_time::FT, final_time::FT: Time span.
• cfl::FT: CFL number.
• regrid_interval::Int: Number of coarse steps between regrids.

compute_dt_amr(prob::AMRProblem, level::Int) -> FT

Compute the time step for a given level based on the CFL condition. The time step is the minimum over all active blocks at this level.

is_leaf(block::AMRBlock) -> Bool

Check if the block is a leaf (has no children).

active_blocks(grid::AMRGrid) -> Vector{AMRBlock}

Return all active (leaf) blocks in the grid, sorted by level.

blocks_at_level(grid::AMRGrid, level::Int) -> Vector{AMRBlock}

Return all active blocks at a specific refinement level.

max_active_level(grid::AMRGrid) -> Int

Return the maximum refinement level among active blocks.

refine_block!(grid::AMRGrid, block_id::Int) -> Vector{Int}

Refine a block by splitting it into 2^Dim child blocks. Returns the IDs of the new child blocks. The parent block is deactivated.

coarsen_block!(grid::AMRGrid, block_id::Int) -> Bool

Coarsen a block by removing all its children and reactivating it. Returns true if coarsening was performed. All children must be active leaves for coarsening to proceed.

regrid!(grid::AMRGrid)

Evaluate refinement criteria on all active blocks and perform refinement/coarsening.

prolongate!(parent::AMRBlock, children::Vector{AMRBlock}, law)

Prolongate (interpolate) the coarse parent solution onto the fine child blocks. Uses conservative linear interpolation with slope limiting.

restrict!(parent::AMRBlock, children::Vector{AMRBlock})

Restrict (average) the fine child solutions onto the coarse parent block. Uses simple volume-weighted averaging (all fine cells have equal volume within a coarse cell for uniform refinement).

FluxRegister{N, FT, Dim}

Accumulates flux corrections at fine/coarse boundaries.

Fields

• coarse_flux::Dict{Symbol, Array{SVector{N,FT}, Dim-1}}: Coarse flux at each face, keyed by direction (:left, :right, etc.).
• fine_flux_sum::Dict{Symbol, Array{SVector{N,FT}, Dim-1}}: Sum of fine fluxes at each coarse face location.
• n_fine::Int: Number of fine faces per coarse face (2 in 2D, 4 in 3D).

reset_flux_register!(reg::FluxRegister)

Zero out all stored fluxes.

accumulate_fine_flux!(reg, face, fine_face_idx, coarse_face_idx, flux)

Add a fine flux contribution at a coarse face location.

store_coarse_flux!(reg, face, coarse_face_idx, flux)

Store the coarse flux at a face location.

apply_flux_correction!(U_coarse, reg, dt, dx, face, law, nx, ny)

Apply the Berger-Colella flux correction to the coarse solution at the interface.

The correction is: dU = dt/dx * (Ffineavg - Fcoarse) where Ffineavg = sum(Ffine) / n_fine.

This is added to the coarse cells adjacent to the interface.

apply_flux_correction_3d!(U, reg, dt, dx, dy, dz, nx, ny, nz, face)

Apply Berger-Colella flux-register corrections to a padded 3D coarse-grid state array U along the coarse-fine interface given by face.

The correction uses the stored difference between the coarse-face flux and the averaged fine-face flux:

\[\Delta U = \frac{\Delta t}{\Delta x_f} (F_{\text{fine,avg}} - F_{\text{coarse}})\]

with the directional spacing dx, dy, or dz chosen from face.

solve_hyperbolic(prob::HyperbolicProblem; method=:ssprk3) -> (x, U_final, t_final)

Solve the 1D hyperbolic problem using explicit time integration.

Keyword Arguments

• method::Symbol: Time integration method. :euler for forward Euler, :ssprk3 for 3rd-order strong stability preserving Runge-Kutta (default).

Returns

• x::Vector: Cell center coordinates.
• U_final::Vector{SVector{N}}: Final conserved variable vectors at cell centers.
• t_final::Real: Final time reached.

solve_hyperbolic(prob::HyperbolicProblem2D{<:IdealMHDEquations{2}}; method=:ssprk3)
    -> (coords, U_final, t_final, ct)

Solve the 2D MHD problem using constrained transport for ∇·B = 0.

Returns

• coords: Cell center coordinates (x, y) matrix.
• U_final: Final conserved variable matrix (nx × ny).
• t_final: Final time reached.
• ct: Final CTData2D (for inspecting ∇·B, face-centered B, etc.).

solve_hyperbolic(prob::HyperbolicProblem2D{<:SRMHDEquations{2}}; method=:ssprk3, vector_potential=nothing)

Solve the 2D SRMHD problem using constrained transport for ∇·B = 0. Same structure as the IdealMHDEquations{2} solver.

solve_hyperbolic(prob::HyperbolicProblem2D{<:GRMHDEquations{2}};
                 method=:ssprk3, vector_potential=nothing)
    -> (coords, U_final, t_final, ct)

Solve the 2D GRMHD problem using the Valencia formulation with:

• Constrained transport for divergence-free B
• Metric-corrected fluxes (Valencia flux: alphaF - betaU)
• Geometric source terms from spacetime curvature

Returns

• coords: Cell center coordinates (x, y) matrix.
• U_final: Final conserved variable matrix (nx x ny, undensitized).
• t_final: Final time reached.
• ct: Final CTData2D for inspecting div(B) and face-centered B.

solve_hyperbolic(prob::HyperbolicProblem3D; method=:ssprk3)
    -> (coords, U_final, t_final)

Solve the 3D hyperbolic problem using explicit time integration.

Returns

• coords::Array{Tuple{Float64,Float64,Float64},3}: Cell center coordinates.
• U_final::Array{SVector{N},3}: Final conserved variable array (nx x ny x nz, interior only).
• t_final::Real: Final time reached.

solve_hyperbolic(prob::HyperbolicProblem3D{<:IdealMHDEquations{3}}; method=:ssprk3)
    -> (coords, U_final, t_final, ct)

Solve the 3D MHD problem using constrained transport for div(B) = 0.

Returns

• coords: Cell center coordinates (x, y, z) array.
• U_final: Final conserved variable array (nx x ny x nz).
• t_final: Final time reached.
• ct: Final CTData3D (for inspecting div(B), face-centered B, etc.).

solve_hyperbolic(prob::UnstructuredHyperbolicProblem; method=:ssprk3)
    -> (centroids, U_final, t_final)

Solve the hyperbolic conservation law on an unstructured mesh.

Returns

• centroids::Vector{Tuple{FT,FT}}: Triangle centroid coordinates.
• U_final::Vector{SVector{N,FT}}: Final conserved variable vector per cell.
• t_final::Real: Final time reached.

solve_hyperbolic_imex(prob::HyperbolicProblem, stiff_source::AbstractStiffSource;
                      scheme=IMEX_SSP3_433(), newton_tol=1e-10, newton_maxiter=5)

Solve a 1D hyperbolic problem with stiff source terms using IMEX Runge-Kutta time integration.

Arguments

• prob::HyperbolicProblem: The hyperbolic problem (provides mesh, law, BCs, etc.).
• stiff_source::AbstractStiffSource: The stiff source term treated implicitly.

Keyword Arguments

• scheme::AbstractIMEXScheme: The IMEX RK scheme (default: IMEX_SSP3_433()).
• newton_tol::Real: Tolerance for implicit Newton solves (default: 1e-10).
• newton_maxiter::Int: Maximum Newton iterations per implicit stage (default: 5).

Returns

• x::Vector: Cell center coordinates.
• U_final::Vector{SVector{N}}: Final conserved variable vectors at cell centers.
• t_final::Real: Final time reached.

solve_hyperbolic_imex(prob::HyperbolicProblem2D, stiff_source::AbstractStiffSource;
                      scheme=IMEX_SSP3_433(), newton_tol=1e-10, newton_maxiter=5,
                      parallel=false)

Solve a 2D hyperbolic problem with stiff source terms using IMEX Runge-Kutta time integration.

Keyword Arguments

• parallel::Bool: Use multi-threaded flux and implicit solve (default: false).

Returns

• coords: Cell center coordinates (nx x ny matrix of tuples).
• U_final::Matrix{SVector{N}}: Final conserved variable matrix (nx x ny).
• t_final::Real: Final time reached.

solve_amr(prob::AMRProblem; method=:subcycling) -> (grid, t_final)

Solve an AMR problem with Berger-Oliger subcycling time integration.

Returns

• grid: The final AMR grid with solution data.
• t_final: The final time reached.

AbstractMetric{Dim}

Abstract supertype for spacetime metrics in Dim spatial dimensions.

All subtypes must implement:

• lapse(m, x, y): Lapse function alpha at position (x, y).
• shift(m, x, y): Shift vector SVector{Dim} at position (x, y).
• spatial_metric(m, x, y): Spatial metric tensor SMatrix{Dim,Dim}.
• sqrt_gamma(m, x, y): Square root of the determinant of the spatial metric.
• inv_spatial_metric(m, x, y): Inverse spatial metric tensor SMatrix{Dim,Dim}.

MinkowskiMetric{Dim} <: AbstractMetric{Dim}

Flat Minkowski spacetime in Dim spatial dimensions.

Lapse alpha = 1, shift beta = 0, spatial metric gammaij = deltaij. Using this metric in GRMHDEquations should recover SRMHDEquations exactly.

SchwarzschildMetric{FT} <: AbstractMetric{2}

Schwarzschild black hole spacetime in Kerr-Schild (horizon-penetrating) coordinates.

The metric is parameterized by the black hole mass M. Coordinates are Cartesian (x, y) with r = sqrt(x^2 + y^2).

A floor radius r_min prevents singularities at r = 0 in numerical simulations.

Fields

• M::FT: Black hole mass.
• r_min::FT: Minimum radius floor (default 0.1 * M).

KerrMetric{FT} <: AbstractMetric{2}

Kerr black hole spacetime in Kerr-Schild (horizon-penetrating) coordinates, restricted to the equatorial plane.

Fields

• M::FT: Black hole mass.
• a::FT: Spin parameter (dimensionless, |a| <= M).
• r_min::FT: Minimum radius floor.

lapse(m::AbstractMetric, x, y) -> Real

Compute the lapse function alpha at position (x, y). The lapse measures the rate of proper time flow relative to coordinate time.

shift(m::AbstractMetric{Dim}, x, y) -> SVector{Dim}

Compute the shift vector beta^i at position (x, y). The shift relates spatial coordinates between neighboring time slices.

spatial_metric(m::AbstractMetric{Dim}, x, y) -> SMatrix{Dim,Dim}

Compute the spatial metric tensor gamma_ij at position (x, y). This is the induced metric on the spatial hypersurface.

sqrt_gamma(m::AbstractMetric, x, y) -> Real

Compute sqrt(det(gamma_ij)) at position (x, y). This factor is used to densitize conserved variables in the Valencia formulation.

inv_spatial_metric(m::AbstractMetric{Dim}, x, y) -> SMatrix{Dim,Dim}

Compute the inverse spatial metric tensor gamma^ij at position (x, y). Satisfies gamma^ik gammakj = delta^ij.

MetricData2D{FT}

Precomputed metric quantities at cell centers for a 2D structured mesh.

All arrays are nx x ny and correspond to interior cells (ix, iy) for ix = 1:nx, iy = 1:ny.

Fields

• alpha::Matrix{FT}: Lapse function.
• beta_x::Matrix{FT}: x-component of shift vector.
• beta_y::Matrix{FT}: y-component of shift vector.
• gamma_xx::Matrix{FT}: (1,1) component of spatial metric.
• gamma_xy::Matrix{FT}: (1,2) component of spatial metric.
• gamma_yy::Matrix{FT}: (2,2) component of spatial metric.
• gammaI_xx::Matrix{FT}: (1,1) component of inverse spatial metric.
• gammaI_xy::Matrix{FT}: (1,2) component of inverse spatial metric.
• gammaI_yy::Matrix{FT}: (2,2) component of inverse spatial metric.
• sqrtg::Matrix{FT}: Square root of spatial metric determinant.

precompute_metric(metric::AbstractMetric{2}, mesh::StructuredMesh2D) -> MetricData2D

Precompute all metric quantities at cell centers for a 2D structured mesh.

This should be called once before the time integration loop for stationary spacetimes.

SRMHDEquations{Dim, EOS <: AbstractEOS} <: AbstractConservationLaw{Dim}

The special relativistic magnetohydrodynamics equations in Dim spatial dimensions.

Variables (8 components in all dimensions)

• Primitive: W = [ρ, vx, vy, vz, P, Bx, By, Bz]
• Conserved: U = [D, Sx, Sy, Sz, τ, Bx, By, Bz]

The conserved↔primitive conversion requires iterative root-finding (Con2Prim).

Fields

• eos::EOS: Equation of state.
• con2prim_tol::Float64: Tolerance for con2prim convergence (default 1e-12).
• con2prim_maxiter::Int: Maximum con2prim iterations (default 50).

GRMHDEquations{Dim, EOS <: AbstractEOS, M <: AbstractMetric{Dim}} <: AbstractConservationLaw{Dim}

The general relativistic magnetohydrodynamics equations in Dim spatial dimensions using the Valencia formulation.

Variables (8 components)

• Primitive: W = [rho, vx, vy, vz, P, Bx, By, Bz]
• Conserved (densitized): U = sqrt(gamma) * [D, Sx, Sy, Sz, tau, Bx, By, Bz]

The physical_flux method returns the flat-space-like flux f^i. The solver applies the Valencia correction alpha * F_riemann - beta * U at each face.

Fields

• eos::EOS: Equation of state.
• metric::M: Spacetime metric.
• con2prim_tol::Float64: Tolerance for con2prim convergence.
• con2prim_maxiter::Int: Maximum con2prim iterations.

lorentz_factor(vx, vy, vz) -> W

Compute the Lorentz factor W = 1/√(1 − v²).

srmhd_con2prim(eos, u::SVector{8}, tol, maxiter) -> (SVector{8}, Con2PrimResult)

Recover primitive variables [ρ, vx, vy, vz, P, Bx, By, Bz] from conserved variables [D, Sx, Sy, Sz, τ, Bx, By, Bz] using the Palenzuela et al. (2015) Newton-Raphson method.

Iterates on auxiliary variable ξ = ρhW² until the energy equation residual f(ξ) = ξ + B² − P_tot − (τ + D) converges to zero.

Con2PrimResult{FT}

Result of a conservative-to-primitive recovery.

Fields

• converged::Bool: Whether the iteration converged.
• iterations::Int: Number of iterations performed.
• residual::FT: Final residual.

grmhd_con2prim(law::GRMHDEquations, u_tilde::SVector{8}, x, y)
    -> (SVector{8}, Con2PrimResult)

Recover primitive variables from densitized conserved variables at position (x, y) using the metric from law.metric.

Steps:

1. Undensitize: u = u_tilde / sqrt(gamma)
2. Compute S^2 = gamma^ij S_i S_j
3. Run the modified Palenzuela Newton-Raphson iteration
4. Return primitives [rho, vx, vy, vz, P, Bx, By, Bz]

grmhd_source_terms(law::GRMHDEquations{2}, w::SVector{8}, u_densitized::SVector{8},
                    md::MetricData2D, mesh, ix::Int, iy::Int) -> SVector{8}

Compute the geometric source terms for the GRMHD equations at cell (ix, iy).

Returns S = (0, Sx_src, Sy_src, Sz_src, tau_src, 0, 0, 0).

The source terms are computed using the stress-energy tensor contracted with metric derivatives (Christoffel symbols). We use centered finite differences for the metric derivatives.