The Finite Volume Method: Theory, Discretization, and General Relativistic Magnetohydrodynamics

Method

The Finite Volume Method

This document provides a comprehensive treatment of the Finite Volume Method (FVM), progressing from foundational continuum mechanics through classical CFD discretization to the full General Relativistic Magnetohydrodynamic (GRMHD) formulation. The hierarchy is organized such that simpler physics emerges as limiting cases of the general theory.

1. Foundations of Continuum Mechanics and Conservation Laws

The mathematical modeling of fluid mechanics and heat transfer rests upon conservation laws. Whether analyzing hypersonic re-entry or slow groundwater percolation, the governing physics remains invariant: mass is conserved, momentum follows Newton’s laws, and energy is transformed but never lost.

The Finite Volume Method (FVM) is distinguished by its direct derivation from integral conservation laws, rather than differential forms. This endows FVM with intrinsic robustness in preserving physical quantities locally and globally, making it the dominant framework in modern Computational Fluid Dynamics (CFD).

1.1 The Reynolds Transport Theorem

The bridge between a system (material volume) and a control volume (fixed in space) is the Reynolds Transport Theorem (RTT). For a scalar property $\phi$ per unit mass:

\frac{D}{Dt} \int_{V_{sys}} \rho \phi \, dV = \frac{\partial}{\partial t} \int_{CV} \rho \phi \, dV + \oint_{CS} \rho \phi (\mathbf{v} \cdot \mathbf{n}) \, dA

FVM starts directly from this integral form, ensuring discrete conservation even on coarse meshes.

1.2 The Navier-Stokes Equations in Integral Form

1.2.1 Conservation of Mass (Continuity)

Setting $\phi = 1$ :

\frac{\partial}{\partial t} \int_{\Omega} \rho \, d\Omega + \oint_{\Gamma} \rho (\mathbf{v} \cdot \mathbf{n}) \, d\Gamma = 0

For incompressible flow ( $\rho = \text{const}$ ), this enforces a divergence-free velocity field.

1.2.2 Conservation of Momentum

Setting $\phi = \mathbf{v}$ :

\frac{\partial}{\partial t} \int_{\Omega} \rho \mathbf{v} \, d\Omega + \oint_{\Gamma} \rho \mathbf{v} (\mathbf{v} \cdot \mathbf{n}) \, d\Gamma = \oint_{\Gamma} \boldsymbol{\tau} \cdot \mathbf{n} \, d\Gamma - \oint_{\Gamma} p \mathbf{n} \, d\Gamma + \int_{\Omega} \rho \mathbf{g} \, d\Omega

Momentum transport is both convective and diffusive.

1.2.3 General Scalar Transport Equation

All conservation laws may be written in unified form:

\underbrace{\frac{\partial}{\partial t} \int_{\Omega} \rho \phi \, d\Omega}_{\text{Transient}} + \underbrace{\oint_{\Gamma} \rho \phi (\mathbf{v} \cdot \mathbf{n}) \, d\Gamma}_{\text{Convection}} = \underbrace{\oint_{\Gamma} \Gamma_\phi \nabla \phi \cdot \mathbf{n} \, d\Gamma}_{\text{Diffusion}} + \underbrace{\int_{\Omega} S_\phi \, d\Omega}_{\text{Source}}

The terms represent:

Transient: Rate of change of $\phi$ inside the control volume. Set to 0 for steady-state.
Convection: Transport of $\phi$ due to fluid velocity. Set to 0 for pure diffusion (e.g., solid conduction).
Diffusion: Transport of $\phi$ due to gradients. Set to 0 for inviscid flows (Euler equations).
Source: Generation or destruction of $\phi$ (e.g., chemical reaction, gravity, heat generation).

1.2.4 Application to Specific Conservation Laws

Conservation Law	$\phi$	$\Gamma_\phi$	$S_\phi$
Mass (Continuity)	1	0	0
$x$ -Momentum	$u$	$\mu$ (Viscosity)	$-\partial p/\partial x + F_x$
$y$ -Momentum	$v$	$\mu$ (Viscosity)	$-\partial p/\partial y + F_y$
Energy	$h$ (Enthalpy)	$k/c_p$	Viscous dissipation, radiation
Species	$c$ (Concentration)	$D$ (Diffusivity)	Reaction rate
Turbulence ( $k$ - $\varepsilon$ )	$k, \varepsilon$	$\mu_t/\sigma_k, \mu_t/\sigma_\varepsilon$	Production/Dissipation

2. Finite Volume Discretization Framework

FVM discretizes the control volume, ensuring that flux entering one cell exits its neighbor exactly (“telescoping flux”), guaranteeing global conservation.

2.1 Mesh Topologies and Geometric Definitions

2.1.1 Structured Grids

Implicit connectivity with efficient memory access, but difficult to generate for complex geometries.

2.1.2 Unstructured and Polyhedral Grids

Arbitrary polygons/polyhedra with connectivity stored explicitly:

Each face has Owner / Neighbor
High geometric flexibility
Higher memory and cache cost

2.1.3 Cell-Centered vs Cell-Vertex Formulations

Cell-Centered FVM (CC-FVM): Variables at cell centroids (most commercial solvers)
Cell-Vertex FVM (CV-FVM): Variables at vertices; dual mesh control volumes

2.2 Geometric Parameters and Surface Approximation

Surface integrals are approximated as:

\oint_{\Gamma} \mathbf{J} \cdot \mathbf{n} \, d\Gamma \approx \sum_f \mathbf{J}_f \cdot \mathbf{S}_f

Key mesh quality metrics:

Orthogonality: Alignment of face normal with cell-center connection
Skewness: Deviation of face center from intersection point

3. Convective Transport: High-Resolution Spatial Reconstruction

The key challenge is determining the face value $\phi_f$ .

3.1 Upwind Differencing Scheme (UDS)

\phi_f = \begin{cases} \phi_P & F_f \ge 0 \\ \phi_N & F_f < 0 \end{cases}

First-order accurate
Unconditionally stable
Strong numerical diffusion

3.2 Central Differencing Scheme (CDS)

\phi_f = f_x \phi_P + (1 - f_x)\phi_N

Second-order accurate
Unbounded in convection-dominated flows

3.3 Godunov’s Theorem

No linear, second-order, monotone scheme exists.

3.4 TVD Schemes and Flux Limiters

\phi_f = \phi_{UDS} + \frac{1}{2}\psi(r)\left(\phi_{CDS} - \phi_{UDS}\right)

Gradient ratio:

r = \frac{\phi_P - \phi_{up}}{\phi_N - \phi_P}

Common Limiters

Limiter	Formula
Minmod	$\psi(r) = \max(0, \min(1, r))$
Superbee	$\psi(r) = \max(0, \min(2r, 1), \min(r, 2))$
Van Leer	$\psi(r) = (r +

4. Diffusive Transport and Gradient Computation

4.1 Gradient Reconstruction Methods

4.1.1 Green-Gauss Method

(\nabla \phi)_P \approx \frac{1}{V_P} \sum_f \phi_f \mathbf{S}_f

Variants: Cell-based, Node-based

4.1.2 Least Squares Method

Assume linear variation:

\phi_{nb} \approx \phi_P + (\nabla \phi)_P \cdot \mathbf{d}_{nb}

Minimize:

E^2 = \sum_{nb} w_{nb} \left[\phi_{nb} - (\phi_P + \nabla \phi_P \cdot \mathbf{d}_{nb})\right]^2

4.2 Non-Orthogonal Correction

Decompose face area vector:

\mathbf{S}_f = \mathbf{E} + \mathbf{T}

Diffusive flux:

\nabla \phi \cdot \mathbf{S}_f \approx |\mathbf{E}|\frac{\phi_N - \phi_P}{|\mathbf{d}|} + (\nabla \phi)_f \cdot \mathbf{T}

5. Temporal Evolution and Stability Analysis

5.1 Explicit Time Integration

Forward Euler:

\frac{\rho V (\phi_P^{n+1} - \phi_P^n)}{\Delta t} = -\sum_f F_f^n + S_\phi^n V

CFL condition:

\text{CFL} = \frac{u \Delta t}{\Delta x} \le 1

5.1.1 Runge-Kutta Methods

Multi-stage explicit schemes
Large stability region
Favored in compressible flows

5.2 Implicit Time Integration

\frac{\rho V (\phi_P^{n+1} - \phi_P^n)}{\Delta t} = -\sum_f F_f^{n+1} + S_\phi^{n+1} V

Unconditionally stable
Requires matrix solvers (AMG)

6. Pressure-Velocity Coupling Algorithms

6.1 Rhie-Chow Interpolation

Prevents checkerboard pressure on collocated grids:

u_f = \overline{u_f} + \overline{D_f}\left(\overline{\nabla p}_f - (\nabla p)_f\right)

6.2 Segregated Solvers

SIMPLE (Semi-Implicit Method for Pressure-Linked Equations)

Predictor-corrector loop
Requires under-relaxation

PISO (Pressure-Implicit with Splitting of Operators)

Two pressure correctors
Efficient for transient flows

7. Boundary Condition Implementation

General discretized equation:

a_P \phi_P + \sum_{nb} a_{nb} \phi_{nb} = S_u

7.1 Dirichlet (Fixed Value)

Ghost-cell extrapolation:

\phi_G = 2\phi_b - \phi_P

7.2 Neumann (Fixed Flux)

J = q_{\text{fixed}} A

Added directly to source term.

7.3 Robin (Mixed)

Example (convective heat transfer):

-k \frac{\partial T}{\partial n} = h (T_b - T_\infty)

Effective boundary temperature:

T_b = \frac{k T_P + h d T_\infty}{k + h d}

7.4 Periodic Boundary Conditions

Implicit matrix coupling, or
Explicit ghost-cell mapping

8. The Arbitrary Lagrangian-Eulerian (ALE) Formulation

For problems involving moving boundaries (fluid-structure interaction, free surfaces), the formulation must account for control volume movement:

\frac{d}{dt}\int_{\Omega(t)} \rho\phi \, d\Omega + \oint_{\partial\Omega(t)} \rho\phi(\mathbf{v} - \mathbf{v}_g) \cdot \mathbf{n} \, dS = \oint_{\partial\Omega(t)} (\Gamma\nabla\phi) \cdot \mathbf{n} \, dS + \int_{\Omega(t)} S_\phi \, d\Omega

Where $\mathbf{v}_g$ is the mesh velocity:

If $\mathbf{v}_g = 0$ : Eulerian (fixed grid)
If $\mathbf{v}_g = \mathbf{v}$ : Lagrangian (grid moves with fluid)

9. Special Relativistic Magnetohydrodynamics (SRMHD)

When fluid velocities approach the speed of light, standard conservation laws must be modified to account for relativistic effects.

9.1 Flat Spacetime Assumptions

Spacetime is Minkowski (flat): $g_{\mu\nu} = \eta_{\mu\nu} = \text{diag}(-1, 1, 1, 1)$

Implications:

Lapse $\alpha = 1$
Shift $\beta^i = 0$
Determinant $\gamma = 1$
All Christoffel symbols $\Gamma^\lambda_{\mu\nu} = 0$
Geometric source terms vanish

9.2 The Lorentz Factor

W = \frac{1}{\sqrt{1 - v^2/c^2}}

This accounts for time dilation and length contraction.

9.3 Relativistic Velocity Addition

Standard addition fails for relativistic speeds. The correct formula:

\lambda_{\text{lab}} = \frac{v_{\text{fluid}} + \lambda'}{1 + v_{\text{fluid}}\lambda'/c^2}

This guarantees $\lambda_{\text{lab}} < c$ when both component velocities are subluminal.

9.4 Conservation Laws in SRMHD

Mass conservation:

\frac{\partial D}{\partial t} + \nabla \cdot (D \mathbf{v}) = 0

where $D = \rho W$ is the relativistic mass density.

10. Newtonian Magnetohydrodynamics

10.1 Limiting Conditions

Weak Gravity: $\alpha \approx 1$ , $\beta^i \approx 0$
Non-Relativistic Velocities: $v \ll c$ , implying $W \approx 1$
Rest-Mass Dominance: $P \ll \rho c^2$ , $B^2 \ll \rho c^2$

10.2 Recovered Equations

Continuity:

\partial_t \rho + \nabla \cdot (\rho \mathbf{v}) = 0

Momentum (Cauchy):

\partial_t(\rho \mathbf{v}) + \nabla \cdot (\rho \mathbf{v} \otimes \mathbf{v}) = -\nabla P + \mathbf{J} \times \mathbf{B} + \rho \mathbf{g}

Energy:

E_{\text{Newt}} \approx \frac{1}{2}\rho v^2 + \rho \epsilon + \frac{B^2}{2}

Induction:

\partial_t \mathbf{B} = \nabla \times (\mathbf{v} \times \mathbf{B})

10.3 The Magnetic Divergence Constraint

\nabla \cdot \mathbf{B} = 0

This constraint must be maintained numerically to avoid unphysical magnetic monopoles.

11. General Relativistic Magnetohydrodynamics (GRMHD)

GRMHD combines fluid dynamics, electromagnetism, and general relativity. It is used to simulate accretion disks around black holes, neutron star mergers, and relativistic jets.

11.1 The 3+1 Formalism

The numerical evolution of relativistic fields requires breaking four-dimensional covariance using the Arnowitt-Deser-Misner (ADM) formalism, which foliates spacetime into three-dimensional spatial hypersurfaces $\Sigma_t$ evolving in time.

11.1.1 The ADM Metric

ds^2 = -\alpha^2 dt^2 + \gamma_{ij} (dx^i + \beta^i dt)(dx^j + \beta^j dt)

Geometric Variables:

Lapse Function ( $\alpha$ ): Relates coordinate time to proper time, $d\tau = \alpha dt$ . Encapsulates gravitational redshift.
Shift Vector ( $\beta^i$ ): Describes relative motion between spatial coordinates and Eulerian observers. Essential for frame-dragging in rotating spacetimes.
Spatial Metric ( $\gamma_{ij}$ ): Measures distances within the three-dimensional hypersurface.

11.1.2 The Eulerian Observer

The four-velocity of the Eulerian observer (at rest relative to the spatial slice):

n^\mu = \frac{1}{\alpha} (1, -\beta^i), \quad n_\mu = (-\alpha, 0, 0, 0)

The metric determinants are related by: $\sqrt{-g} = \alpha \sqrt{\gamma}$

11.2 The Stress-Energy Tensor

11.2.1 Perfect Fluid

T^{\mu\nu}_{\text{fluid}} = \rho h u^\mu u^\nu + P g^{\mu\nu}

where $h = 1 + \epsilon + P/\rho$ is the specific enthalpy.

11.2.2 Electromagnetic Field (Ideal MHD)

T^{\mu\nu}_{\text{EM}} = b^2 u^\mu u^\nu + \frac{1}{2}b^2 g^{\mu\nu} - b^\mu b^\nu

where $b^\mu$ is the magnetic field four-vector and $b^2 = b^\mu b_\mu$ .

11.2.3 Total Stress-Energy Tensor

T^{\mu\nu} = (\rho h + b^2) u^\mu u^\nu + \left(P + \frac{1}{2}b^2\right)g^{\mu\nu} - b^\mu b^\nu

The magnetic field contributes:

Effective enthalpy ( $b^2$ term)
Isotropic magnetic pressure ( $b^2/2$ )
Anisotropic magnetic tension ( $-b^\mu b^\nu$ )

11.3 The Valencia Formulation (Flux-Conservative Form)

\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}^i}{\partial x^i} = \mathbf{S}

11.3.1 Conserved Variables

\mathbf{U} = \sqrt{\gamma} \begin{bmatrix} D \\ S_j \\ \tau \\ B^k \end{bmatrix}

where:

Relativistic Mass Density:

D = \rho W

Momentum Density:

S_j = (\rho h + b^2) W^2 v_j - \alpha b^0 b_j

Energy Density (minus rest mass):

\tau = (\rho h + b^2) W^2 - \left(P + \frac{1}{2}b^2\right) - \alpha^2 (b^0)^2 - D

11.3.2 Flux Vectors

\mathbf{F}^i = \sqrt{\gamma} \begin{bmatrix} D (v^i - \beta^i / \alpha) \\ S_j (v^i - \beta^i / \alpha) + P_{\text{tot}} \delta^i_j - b^i b_j / W^2 \\ \tau (v^i - \beta^i / \alpha) + P_{\text{tot}} v^i - \alpha b^0 b^i / W^2 \\ B^i (v^i - \beta^i / \alpha) - B^j (v^j - \beta^j / \alpha) \end{bmatrix}

where $P_{\text{tot}} = P + b^2/2$ is the total pressure.

11.3.3 Geometric Source Terms

\mathbf{S} = \sqrt{\gamma} \begin{bmatrix} 0 \\ \frac{1}{2} \alpha T^{\mu\nu} \partial_j g_{\mu\nu} + S_k \partial_j \beta^k - (\tau + D) \partial_j \alpha \\ \alpha T^{\mu\nu} \nabla_\nu n_\mu \\ 0 \end{bmatrix}

The source terms arise from Christoffel symbols and represent the exchange of momentum and energy between matter and the gravitational field.

11.4 The Simplification Hierarchy

To Recover…	Set…	Result
SRMHD	$\alpha=1$ , $\beta^i=0$ , $\gamma_{ij}=\delta_{ij}$	Minkowski (Flat) Space
Newtonian MHD	$v \ll c$ , $P \ll \rho c^2$ , $W \to 1$	Standard Ideal MHD
Navier-Stokes	$B=0$ , add viscosity	Compressible NS
Euler Equations	Viscosity $\to 0$	Compressible Euler
Heat Equation	$v=0$ , $S_{\text{geom}}=0$	Pure Diffusion

12. Constrained Transport for the Magnetic Field

12.1 The Divergence Constraint Problem

Numerical truncation errors introduce $\nabla \cdot \mathbf{B} \neq 0$ , creating unphysical magnetic monopoles that cause simulation instability.

12.2 Staggered Grid Variable Placement

Cell Centers $(i,j,k)$ : Volume-averaged quantities ( $D$ , $P$ , $\tau$ )
Face Centers $(i+1/2, j, k)$ : Area-averaged magnetic flux ( $B^x$ on $x$ -face)
Edge Centers $(i+1/2, j+1/2, k)$ : Line-averaged EMF ( $\mathcal{E}^z$ on $z$ -edge)

12.3 The Flux Update via Stokes’ Theorem

The evolution of magnetic flux $\Phi$ through a face is governed by Faraday’s Law:

\frac{\partial \Phi_B}{\partial t} = -\oint_{\partial\text{Face}} \mathbf{E} \cdot d\mathbf{l}

For the $x$ -face:

\frac{\partial B^x_{i+1/2,j,k}}{\partial t} = -\left( \frac{\mathcal{E}^z_{i+1/2,j+1/2,k} - \mathcal{E}^z_{i+1/2,j-1/2,k}}{\Delta y} - \frac{\mathcal{E}^y_{i+1/2,j,k+1/2} - \mathcal{E}^y_{i+1/2,j,k-1/2}}{\Delta z} \right)

12.4 Computing the Relativistic EMF

In curved spacetime:

\mathcal{E}_i = \epsilon_{ijk}(\alpha v^j - \beta^j)B^k

where:

$\alpha v^j$ : Fluid velocity scaled by time dilation
$\beta^j$ : Coordinate grid sliding speed (shift vector)

13. Riemann Solvers for Relativistic MHD

13.1 The Relativistic Wave Fan

When two cells with different states touch, they generate a fan of 7 waves:

Fast Magnetosonic Waves (Left & Right)
Alfven Waves (Left & Right)
Slow Magnetosonic Waves (Left & Right)
Contact Discontinuity (Center)

13.2 The HLL Solver

Assumes only two waves separating Left and Right states:

\mathbf{F}_{\text{HLL}} = \frac{\lambda_R \mathbf{F}_L - \lambda_L \mathbf{F}_R + \lambda_L \lambda_R (\mathbf{U}_R - \mathbf{U}_L)}{\lambda_R - \lambda_L}

Very diffusive (lumps intermediate waves)
Extremely stable
Preserves positivity of density and pressure

13.3 The HLLD Solver (5-Wave Model)

Approximates the Riemann fan with five waves:

Fast Shocks ( $\lambda_L$ , $\lambda_R$ ): Outermost waves at relativistic fast magnetosonic speed
Alfven Waves ( $\lambda_L^*$ , $\lambda_R^*$ ): Rotational discontinuities where magnetic field direction changes
Contact Discontinuity ( $\lambda_M$ ): Central wave where density changes but pressure is continuous

The HLLD solver provides significantly lower dissipation than HLL while maintaining stability.

13.4 Causality Enforcement

Wave speeds are calculated using relativistic eigenvalues that guarantee subluminal propagation:

\lambda_{\text{lab}} = \frac{v + \lambda'}{1 + v\lambda'/c^2}

14. Primitive Variable Recovery (Con2Prim)

14.1 The Inversion Problem

The evolution updates conserved variables $\mathbf{U} = (D, S_j, \tau, B^k)$ , but flux calculation requires primitive variables $\mathbf{P} = (\rho, v^i, P, B^k)$ .

The inverse mapping requires solving transcendental equations for the Lorentz factor:

f(W) = W - \frac{1}{\sqrt{1 - \frac{S^2}{(\rho h W^2)^2}}} = 0

Newton-Raphson iteration is standard.

14.2 Floors and Fail-Safes

If the solver fails or returns negative pressure/density:

Variables reset to “atmosphere” values (small positive density/pressure)
Momentum rescaled to maintain subluminal velocity ( $v < 1$ )

15. Equations of State

15.1 Gamma-Law (Ideal Gas)

P = (\Gamma - 1)\rho \epsilon

$\Gamma = 4/3$ : Radiation-dominated gas
$\Gamma = 5/3$ : Non-relativistic gas

15.2 Hybrid EoS

Combines “cold” nuclear physics (tabulated) with “thermal” ideal gas for shock heating. Used in neutron star mergers.

15.3 Tabulated EoS

Large lookup tables including neutrino cooling and nuclear composition. Required for microphysical accuracy.

16. Time Integration for GRMHD

16.1 Explicit Runge-Kutta

Standard choices: SSP-RK3 (Strong Stability Preserving) or RK4.

Efficient
Limited by CFL condition

16.2 IMEX Schemes

For resistive GRMHD with stiff source terms:

Stiff terms treated implicitly
Flux terms remain explicit

GRMHD simulations span vast dynamic ranges (event horizon to jet lobes). AMR creates hierarchies of nested grids.

Challenge: Maintaining $\nabla \cdot \mathbf{B} = 0$ across refinement boundaries requires divergence-preserving prolongation operators.

18. Advanced Topics

18.1 WENO on Unstructured Meshes

Arbitrarily high order
Nonlinear stencil weighting
Shock-capturing with low dissipation

18.2 Machine Learning Accelerated FVM

ML replaces turbulence closures
Conservation laws enforced by FVM structure
Neural networks as constitutive models

18.3 The Israel-Stewart Formulation (Viscous Relativistic Fluids)

Standard viscosity violates causality in relativity. The Israel-Stewart formulation elevates the viscous stress tensor to a dynamic variable:

\tau_\pi \dot{\pi}^{\mu\nu} + \pi^{\mu\nu} = -2\eta \sigma^{\mu\nu}

Taking $\tau_\pi \to 0$ recovers the relativistic Navier-Stokes structure.

References

3+1 formalism and bases of numerical relativity - Observatoire de Paris, https://people-lux.obspm.fr/gourgoulhon/pdf/form3p1.pdf
Introduction to GRMHD - Einstein Toolkit, https://einsteintoolkit.github.io/et2021uiuc/lectures/21-VassiliosMewes/slides.pdf
General-relativistic Resistive Magnetohydrodynamics - Research Explorer, https://pure.uva.nl/ws/files/45365323/General_relativistic_Resistive_Magnetohydrodynamics.pdf
A five-wave HLL Riemann solver for relativistic MHD - arXiv, https://arxiv.org/abs/0811.1483
The Black Hole Accretion Code: adaptive mesh refinement and constrained transport, https://d-nb.info/1363398083/34
IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes, https://par.nsf.gov/servlets/purl/10293408
Theories of Relativistic Dissipative Fluid Dynamics - MDPI, https://www.mdpi.com/1099-4300/26/3/189