diff --git a/.typos.toml b/.typos.toml
index 72dce1b905..432385eef0 100644
--- a/.typos.toml
+++ b/.typos.toml
@@ -30,4 +30,4 @@ unknwn = "unknwn" # typo for "unknown" - tests unknown param detection
tru = "tru" # typo for "true" in "when_tru" - tests dependency keys
[files]
-extend-exclude = ["docs/documentation/references*", "tests/", "toolchain/cce_simulation_workgroup_256.sh", "build-docs/"]
+extend-exclude = ["docs/documentation/references*", "docs/references.bib", "tests/", "toolchain/cce_simulation_workgroup_256.sh", "build-docs/"]
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9e3fef30d1..f2ad948004 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -720,6 +720,7 @@ if (MFC_DOCUMENTATION)
set(DOXYGEN_HTML_OUTPUT "\"${CMAKE_CURRENT_BINARY_DIR}/${target}\"")
set(DOXYGEN_MATHJAX_CODEFILE "\"${CMAKE_CURRENT_SOURCE_DIR}/docs/config.js\"")
set(DOXYGEN_PROJECT_LOGO "\"${CMAKE_CURRENT_SOURCE_DIR}/docs/res/icon.ico\"")
+ set(DOXYGEN_CITE_BIB_FILES "\"${CMAKE_CURRENT_SOURCE_DIR}/docs/references.bib\"")
set(DOXYGEN_IMAGE_PATH "\"${CMAKE_CURRENT_SOURCE_DIR}/docs/res\"\
\"${CMAKE_CURRENT_SOURCE_DIR}/docs/${target}\"")
@@ -752,6 +753,8 @@ if (MFC_DOCUMENTATION)
"${opt_parameters_dependency}"
"${${target}_SRCs}" "${${target}_DOCs}"
COMMAND "${DOXYGEN_EXECUTABLE}" "${target}-Doxyfile"
+ COMMAND "${Python3_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/docs/postprocess_citations.py"
+ "${CMAKE_CURRENT_BINARY_DIR}/${target}"
WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}"
COMMENT "${target}: Generating documentation"
)
@@ -771,6 +774,7 @@ if (MFC_DOCUMENTATION)
add_custom_target(documentation)
find_package(Doxygen REQUIRED dot REQUIRED)
+ find_package(Python3 REQUIRED COMPONENTS Interpreter)
# > Fetch CSS Theme
ExternalProject_Add(doxygen-awesome-css
diff --git a/docs/Doxyfile.in b/docs/Doxyfile.in
index 26e1a0156b..5a32db918d 100644
--- a/docs/Doxyfile.in
+++ b/docs/Doxyfile.in
@@ -740,7 +740,7 @@ LAYOUT_FILE =
# LATEX_BIB_STYLE. To use this feature you need bibtex and perl available in the
# search path. See also \cite for info how to create references.
-CITE_BIB_FILES =
+CITE_BIB_FILES = @DOXYGEN_CITE_BIB_FILES@
#---------------------------------------------------------------------------
# Configuration options related to warning and progress messages
diff --git a/docs/documentation/case.md b/docs/documentation/case.md
index f9073406d1..d0c6431839 100644
--- a/docs/documentation/case.md
+++ b/docs/documentation/case.md
@@ -381,7 +381,7 @@ The parameters define material's property of compressible fluids that are used i
- `fluid_pp(i)%%Re(1)` and `fluid_pp(i)%%Re(2)` define the shear and volume viscosities of $i$-th fluid, respectively.
When these parameters are undefined, fluids are treated as inviscid.
-Details of implementation of viscosity in MFC can be found in [Coralic (2015)](@ref references).
+Details of implementation of viscosity in MFC can be found in \cite Coralic15.
- `fluid_pp(i)%%cv`, `fluid_pp(i)%%qv`, and `fluid_pp(i)%%qvp` define $c_v$, $q$, and $q'$ as parameters of $i$-th fluid that are used in stiffened gas equation of state.
@@ -419,9 +419,9 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](@
| `ic_eps` | Real | Interface compression threshold (default: 1e-4) |
| `ic_beta` | Real | Interface compression sharpness parameter (default: 1.6) |
| `riemann_solver` | Integer | Riemann solver algorithm: [1] HLL*; [2] HLLC; [3] Exact*; [4] HLLD (only for MHD) |
-| `low_Mach` | Integer | Low Mach number correction for HLLC Riemann solver: [0] None; [1] Pressure (Chen et al. 2022); [2] Velocity (Thornber et al. 2008) |
+| `low_Mach` | Integer | Low Mach number correction for HLLC Riemann solver: [0] None; [1] Pressure (\cite Chen22); [2] Velocity (\cite Thornber08) |
| `avg_state` | Integer | Averaged state evaluation method: [1] Roe average*; [2] Arithmetic mean |
-| `wave_speeds` | Integer | Wave-speed estimation: [1] Direct (Batten et al. 1997); [2] Pressure-velocity* (Toro 1999) |
+| `wave_speeds` | Integer | Wave-speed estimation: [1] Direct (\cite Batten97); [2] Pressure-velocity* (\cite Toro09) |
| `weno_Re_flux` | Logical | Compute velocity gradient using scalar divergence theorem |
| `weno_avg` | Logical | Arithmetic mean of left and right, WENO-reconstructed, cell-boundary values |
| `dt` | Real | Time step size |
@@ -438,7 +438,7 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](@
| `surface_tension` | Logical | Activate surface tension |
| `viscous` | Logical | Activate viscosity |
| `hypoelasticity` | Logical | Activate hypoelasticity* |
-| `igr` | Logical | Enable solution via information geometric regularization (IGR) [Cao (2024)](@ref references) |
+| `igr` | Logical | Enable solution via information geometric regularization (IGR) \cite Cao24 |
| `igr_order` | Integer | Order of reconstruction for IGR [3,5] |
| `alf_factor` | Real | Alpha factor for IGR entropic pressure (default 10) |
| `igr_pres_lim` | Logical | Limit IGR pressure to avoid negative values (default F) |
@@ -452,8 +452,8 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](@
The table lists simulation algorithm parameters.
The parameters are used to specify options in algorithms that are used to integrate the governing equations of the multi-component flow based on the initial condition.
-Models and assumptions that are used to formulate and discretize the governing equations are described in [Bryngelson et al. (2019)](@ref references).
-Details of the simulation algorithms and implementation of the WENO scheme can be found in [Coralic (2015)](@ref references).
+Models and assumptions that are used to formulate and discretize the governing equations are described in \cite Wilfong26 and \cite Bryngelson21.
+Details of the simulation algorithms and implementation of the WENO scheme can be found in \cite Coralic15.
- `bc_[x,y,z]%[beg,end]` specifies the boundary conditions at the beginning and the end of domain boundaries in each coordinate direction by a negative integer from -1 through -16.
See table [Boundary Conditions](#boundary-conditions) for details.
@@ -467,12 +467,12 @@ Tangential velocities require viscosity, `weno_avg = T`, and `bc_[x,y,z]%%beg =
Tangential velocities require viscosity, `weno_avg = T`, and `bc_[x,y,z]%%end = 16` to work properly. Normal velocities require `bc_[x,y,z]%%end = -15` or `\bc_[x,y,z]%%end = -16` to work properly.
- `model_eqns` specifies the choice of the multi-component model that is used to formulate the dynamics of the flow using integers from 1 through 3.
-`model_eqns = 1`, `2`, and `3` correspond to $\Gamma$-$\Pi_\infty$ model ([Johnsen, 2008](@ref references)), 5-equation model ([Allaire et al., 2002](@ref references)), and 6-equation model ([Saurel et al., 2009](@ref references)), respectively.
-The difference of the two models is assessed by ([Schmidmayer et al., 2019](@ref references)).
+`model_eqns = 1`, `2`, and `3` correspond to $\Gamma$-$\Pi_\infty$ model (\cite Johnsen08), 5-equation model (\cite Allaire02), and 6-equation model (\cite Saurel09), respectively.
+The difference of the two models is assessed by (\cite Schmidmayer20).
Note that some code parameters are only compatible with 5-equation model.
- `alt_soundspeed` activates the source term in the advection equations for the volume fractions, $K\nabla\cdot \underline{u}$, that regularizes the speed of sound in the mixture region when the 5-equation model is used.
-The effect and use of the source term are assessed by [Schmidmayer et al., 2019](@ref references).
+The effect and use of the source term are assessed by \cite Schmidmayer20.
- `adv_n` activates the direct computation of number density by the Riemann solver instead of computing number density from the void fraction in the method of classes.
@@ -481,7 +481,7 @@ The effect and use of the source term are assessed by [Schmidmayer et al., 2019]
- `mixture_err` activates correction of solutions to avoid imaginary speed of sound at each grid cell.
- `time_stepper` specifies the order of the Runge-Kutta (RK) time integration scheme that is used for temporal integration in simulation, from the 1st to 5th order by corresponding integer.
-Note that `time_stepper = 3` specifies the total variation diminishing (TVD), third order RK scheme ([Gottlieb and Shu, 1998](@ref references)).
+Note that `time_stepper = 3` specifies the total variation diminishing (TVD), third order RK scheme (\cite Gottlieb98).
- `adap_dt` activates the Strang operator splitting scheme which splits flux and source terms in time marching, and an adaptive time stepping strategy is implemented for the source term. It requires ``bubbles_euler = 'T'``, ``polytropic = 'T'``, ``adv_n = 'T'`` and `time_stepper = 3`. Additionally, it can be used with ``bubbles_lagrange = 'T'`` and `time_stepper = 3`. `adap_dt_tol` and `adap_dt_max_iters` are 1e-4 and 100, respectively, by default.
@@ -490,19 +490,19 @@ Note that `time_stepper = 3` specifies the total variation diminishing (TVD), th
- `weno_eps` specifies the lower bound of the WENO nonlinear weights.
It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$ for other WENO variants.
-- `mapped_weno` activates the WENO-M scheme in place of the default WENO-JS scheme ([Henrick et al., 2005](@ref references)). WENO-M a variant of the WENO scheme that remaps the nonlinear WENO-JS weights by assigning larger weights to non-smooth stencils, reducing dissipation compared to the default WENO-JS scheme, at the expense of higher computational cost. Only one of `mapped_weno`, `wenoz`, and `teno` can be activated.
+- `mapped_weno` activates the WENO-M scheme in place of the default WENO-JS scheme (\cite Henrick05). WENO-M a variant of the WENO scheme that remaps the nonlinear WENO-JS weights by assigning larger weights to non-smooth stencils, reducing dissipation compared to the default WENO-JS scheme, at the expense of higher computational cost. Only one of `mapped_weno`, `wenoz`, and `teno` can be activated.
-- `wenoz` activates the WENO-Z scheme in place of the default WENO-JS scheme ([Borges et al., 2008](@ref references)). WENO-Z is a variant of the WENO scheme that further reduces the dissipation compared to the WENO-M scheme. It has similar computational cost to the WENO-JS scheme.
+- `wenoz` activates the WENO-Z scheme in place of the default WENO-JS scheme (\cite Borges08). WENO-Z is a variant of the WENO scheme that further reduces the dissipation compared to the WENO-M scheme. It has similar computational cost to the WENO-JS scheme.
- `wenoz_q` specifies the power parameter `q` used in the WENO-Z scheme. It controls how aggressively the smoothness coefficients scale the weights. A higher value of `wenoz_q` increases the sensitivity to smoothness, improving stability but worsening numerical dissipation. For WENO3 and WENO5, `q=1` is fixed, so `wenoz_q` must not be set. For WENO7, `wenoz_q` can be set to 2, 3, or 4.
-- `teno` activates the TENO scheme in place of the default WENO-JS scheme ([Fu et al., 2016](@ref references)). TENO is a variant of the ENO scheme that is the least dissipative, but could be less robust for extreme cases. It uses a threshold to identify smooth and non-smooth stencils, and applies optimal weights to the smooth stencils. Only available for `weno_order = 5` and `7`. Requires `teno_CT` to be set. Does not support grid stretching.
+- `teno` activates the TENO scheme in place of the default WENO-JS scheme (\cite Fu16). TENO is a variant of the ENO scheme that is the least dissipative, but could be less robust for extreme cases. It uses a threshold to identify smooth and non-smooth stencils, and applies optimal weights to the smooth stencils. Only available for `weno_order = 5` and `7`. Requires `teno_CT` to be set. Does not support grid stretching.
- `teno_CT` specifies the threshold for the TENO scheme. This dimensionless constant, also known as $C_T$, sets a threshold to identify smooth and non-smooth stencils. Larger values make the scheme more robust but also more dissipative. A recommended value for teno_CT is `1e-6`. When adjusting this parameter, it is recommended to try values like `1e-5` or `1e-7` for TENO5. A smaller value can be used for TENO7.
- `null_weights` activates nullification of the nonlinear WENO weights at the buffer regions outside the domain boundaries when the Riemann extrapolation boundary condition is specified (`bc_[x,y,z]%%beg[end]} = -4`).
-- `mp_weno` activates monotonicity preservation in the WENO reconstruction (MPWENO) such that the values of reconstructed variables do not reside outside the range spanned by WENO stencil ([Balsara and Shu, 2000](@ref references); [Suresh and Huynh, 1997](@ref references)).
+- `mp_weno` activates monotonicity preservation in the WENO reconstruction (MPWENO) such that the values of reconstructed variables do not reside outside the range spanned by WENO stencil (\cite Balsara00; \cite Suresh97).
- `muscl_order` specifies the order of the MUSCL scheme that is used for spatial reconstruction of variables by an integer of 1, or 2, that corresponds to the 1st, and 2nd order respectively. When using `muscl_order = 2`, `muscl_lim` must be defined.
@@ -512,16 +512,16 @@ It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$
- `int_comp` activates interface compression using THINC used in MUSCL Reconstruction, with control parameters (`ic_eps`, and `ic_beta`).
- `riemann_solver` specifies the choice of the Riemann solver that is used in simulation by an integer from 1 through 4.
-`riemann_solver = 1`, `2`, and `3` correspond to HLL, HLLC, and Exact Riemann solver, respectively ([Toro, 2013](@ref references)).
-`riemann_solver = 4` is only for MHD simulations. It resolves 5 of the full seven-wave structure of the MHD equations ([Miyoshi and Kusano, 2005](@ref references)).
+`riemann_solver = 1`, `2`, and `3` correspond to HLL, HLLC, and Exact Riemann solver, respectively (\cite Toro09).
+`riemann_solver = 4` is only for MHD simulations. It resolves 5 of the full seven-wave structure of the MHD equations (\cite Miyoshi05).
-- `low_Mach` specifies the choice of the low Mach number correction scheme for the HLLC Riemann solver. `low_Mach = 0` is default value and does not apply any correction scheme. `low_Mach = 1` and `2` apply the anti-dissipation pressure correction method ([Chen et al., 2022](@ref references)) and the improved velocity reconstruction method ([Thornber et al., 2008](@ref references)). This feature requires `model_eqns = 2` or `3`. `low_Mach = 1` works for `riemann_solver = 1` and `2`, but `low_Mach = 2` only works for `riemann_solver = 2`.
+- `low_Mach` specifies the choice of the low Mach number correction scheme for the HLLC Riemann solver. `low_Mach = 0` is default value and does not apply any correction scheme. `low_Mach = 1` and `2` apply the anti-dissipation pressure correction method (\cite Chen22) and the improved velocity reconstruction method (\cite Thornber08). This feature requires `model_eqns = 2` or `3`. `low_Mach = 1` works for `riemann_solver = 1` and `2`, but `low_Mach = 2` only works for `riemann_solver = 2`.
- `avg_state` specifies the choice of the method to compute averaged variables at the cell-boundaries from the left and the right states in the Riemann solver by an integer of 1 or 2.
`avg_state = 1` and `2` correspond to Roe- and arithmetic averages, respectively.
- `wave_speeds` specifies the choice of the method to compute the left, right, and middle wave speeds in the Riemann solver by an integer of 1 and 2.
-`wave_speeds = 1` and `2` correspond to the direct method ([Batten et al., 1997](@ref references)), and indirect method that approximates the pressures and velocity ([Toro, 2013](@ref references)), respectively.
+`wave_speeds = 1` and `2` correspond to the direct method (\cite Batten97), and indirect method that approximates the pressures and velocity (\cite Toro09), respectively.
- `weno_Re_flux` activates the scalar divergence theorem in computing the velocity gradients using WENO-reconstructed cell boundary values.
If this option is false, velocity gradient is computed using finite difference scheme of order 2 which is independent of the WENO order.
@@ -702,7 +702,7 @@ It also cannot be enabled with `flux_wrt`, `heat_ratio_wrt`, `pres_inf_wrt`, `c_
| `acoustic(i)%%bb_bandwidth` | Real | The bandwidth of each frequency in the broadband wave |
| `acoustic(i)%%bb_lowest_freq` | Real | The lower frequency bound of the broadband wave |
-Details of the transducer acoustic source model can be found in [Maeda and Colonius (2017)](@ref references).
+Details of the transducer acoustic source model can be found in \cite Maeda17.
- `acoustic_source` activates the acoustic source module.
@@ -714,7 +714,7 @@ Details of the transducer acoustic source model can be found in [Maeda and Colon
- `%%loc(j)` specifies the location of the acoustic source in the $j$-th coordinate direction. For planer support, the location defines midpoint of the source plane. For transducer arrays, the location defines the center of the transducer or transducer array (not the focal point; for 3D it's the tip of the spherical cap, for 2D it's the tip of the arc).
-- `%%pulse` specifies the acoustic wave form. `%%pulse = 1`, `2`, `3` and `4` correspond to sinusoidal wave, Gaussian wave, square wave and broadband wave, respectively. The implementation of the broadband wave is based on [Tam (2005)](@ref references)
+- `%%pulse` specifies the acoustic wave form. `%%pulse = 1`, `2`, `3` and `4` correspond to sinusoidal wave, Gaussian wave, square wave and broadband wave, respectively. The implementation of the broadband wave is based on \cite Tam05
- `%%npulse` specifies the number of cycles of the acoustic wave generated. Only applies to `%%pulse = 1 and 3` (sine and square waves), and must be an integer for non-planar waves.
@@ -795,7 +795,7 @@ This table lists the sub-grid bubble model parameters, which can be utilized in
- `bub_pp` specifies simulation parameters for the EE and/or EL bubble model.
-Implementation of the parameters into the model follows [Ando (2010)](@ref references).
+Implementation of the parameters into the model follows \cite Ando10.
#### 9.1 Ensemble-Averaged Bubble Model
@@ -816,10 +816,10 @@ Implementation of the parameters into the model follows [Ando (2010)](@ref refer
This table lists the ensemble-averaged bubble model parameters.
- `polytropic` activates polytropic gas compression in the bubble.
-When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due to heat and mass transfer across the bubble wall with constant heat and mass transfer coefficients based on ([Preston et al., 2007](@ref references)).
+When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due to heat and mass transfer across the bubble wall with constant heat and mass transfer coefficients based on (\cite Preston07).
- `thermal` specifies a model for heat transfer across the bubble interface by an integer from 1 through 3.
-`thermal = 1`, `2`, and `3` correspond to no heat transfer (adiabatic gas compression), isothermal heat transfer, and heat transfer with a constant heat transfer coefficient based on [Preston et al., 2007](@ref references), respectively.
+`thermal = 1`, `2`, and `3` correspond to no heat transfer (adiabatic gas compression), isothermal heat transfer, and heat transfer with a constant heat transfer coefficient based on \cite Preston07, respectively.
- `polydisperse` activates polydispersity in the bubble model through a probability density function (PDF) of the equilibrium bubble radius. Simpson's rule is used for integrating the log-normal PDF of equilibrium bubble radius for polydisperse populations.
@@ -857,15 +857,15 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
- `nBubs_glb` Total number of bubbles. Their initial conditions need to be specified in the ./input/lag_bubbles.dat file. See the example cases for additional information.
-- `solver_approach` Specifies the Euler-Lagrange coupling method: [1] enables a one-way coupling approach, where the bubbles do not influence the Eulerian field. [2] activates the two-way coupling approach based on [Maeda and Colonius (2018)](@ref references), where the effect of the bubbles is added in the Eulerian field as source terms.
+- `solver_approach` Specifies the Euler-Lagrange coupling method: [1] enables a one-way coupling approach, where the bubbles do not influence the Eulerian field. [2] activates the two-way coupling approach based on \cite Maeda18, where the effect of the bubbles is added in the Eulerian field as source terms.
-- `cluster_type` Specifies method to find p_inf (pressure that drives the bubble dynamics): [1] activates the bilinear interpolation of the pressure field, while [2] enables the bubble dynamic closure based on [Maeda and Colonius (2018)](@ref references), the full model is obtained when `pressure_corrector` is true.
+- `cluster_type` Specifies method to find p_inf (pressure that drives the bubble dynamics): [1] activates the bilinear interpolation of the pressure field, while [2] enables the bubble dynamic closure based on \cite Maeda18, the full model is obtained when `pressure_corrector` is true.
-- `smooth_type` Specifies the smoothening method of projecting the lagrangian bubbles in the Eulerian field: [1] activates the gaussian kernel function described in [Maeda and Colonius (2018)](@ref references), while [2] activates the delta kernel function where the effect of the bubble is only seen in the specific bubble location cell.
+- `smooth_type` Specifies the smoothening method of projecting the lagrangian bubbles in the Eulerian field: [1] activates the gaussian kernel function described in \cite Maeda18, while [2] activates the delta kernel function where the effect of the bubble is only seen in the specific bubble location cell.
-- `heatTransfer_model` Activates the heat transfer model at the bubble's interface based on ([Preston et al., 2007](@ref references)).
+- `heatTransfer_model` Activates the heat transfer model at the bubble's interface based on (\cite Preston07).
-- `massTransfer_model` Activates the mass transfer model at the bubble's interface based on ([Preston et al., 2007](@ref references)).
+- `massTransfer_model` Activates the mass transfer model at the bubble's interface based on (\cite Preston07).
### 10. Velocity Field Setup
@@ -899,7 +899,7 @@ The parameters are optionally used to define initial velocity profiles and pertu
$$ u = \mbox{patch\_icpp(1)\%vel(1)} * \tanh( y_{cc} * \mbox{mixlayer\_vel\_coef}) $$
-- `mixlayer_perturb` activates the velocity perturbation for a temporal mixing layer with hyperbolic tangent mean streamwise velocity profile, using an inverter version of the spectrum-based synthetic turbulence generation method proposed by Guo et al. (2023, JFM). This option only works for `p > 0` and `mixlayer_vel_profile = 'T'`.
+- `mixlayer_perturb` activates the velocity perturbation for a temporal mixing layer with hyperbolic tangent mean streamwise velocity profile, using an inverter version of the spectrum-based synthetic turbulence generation method proposed by \cite Guo23. This option only works for `p > 0` and `mixlayer_vel_profile = 'T'`.
### 11. Phase Change Model
| Parameter | Type | Description |
@@ -957,7 +957,7 @@ By convention, positive accelerations in the `x[y,z]` direction are in the posit
- `relativity` only works for `mhd` enabled and activates relativistic MHD (RMHD) simulation.
-- `hyper_cleaning` [Dedner et al., 2002](@ref references) only works with `mhd` in 2D/3D and reduces numerical `div B` errors by propagation and damping. Currently not compatible with HLLD (`riemann_solver = 4`).
+- `hyper_cleaning` \cite Dedner02 only works with `mhd` in 2D/3D and reduces numerical `div B` errors by propagation and damping. Currently not compatible with HLLD (`riemann_solver = 4`).
- `hyper_cleaning_speed` sets the propagation speed of divergence-cleaning waves.
@@ -982,7 +982,7 @@ Note: For relativistic flow, the conservative and primitive densities are differ
| `cont_damage_s` | Real | Power `s` for continuum damage model |
| `alpha_bar` | Real | Damage factor (rate) for continuum damage model |
-- `cont_damage` activates continuum damage model for solid materials. Requires `tau_star`, `cont_damage_s`, and `alpha_bar` to be set (empirically determined) ([Cao et al., 2019](@ref references)).
+- `cont_damage` activates continuum damage model for solid materials. Requires `tau_star`, `cont_damage_s`, and `alpha_bar` to be set (empirically determined) (\cite Cao19).
### 16. Cylindrical Coordinates
@@ -1042,7 +1042,7 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met:
The boundary condition supported by the MFC are listed in table [Boundary Conditions](#boundary-conditions).
Their number (`#`) corresponds to the input value in `input.py` labeled `bc_[x,y,z]%[beg,end]` (see table [Simulation Algorithm Parameters](#5-simulation-algorithm)).
-The entries labeled "Characteristic." are characteristic boundary conditions based on [Thompson (1987)](@ref references) and [Thompson (1990)](@ref references).
+The entries labeled "Characteristic." are characteristic boundary conditions based on \cite Thompson87 and \cite Thompson90.
### Generalized Characteristic Boundary conditions
@@ -1058,7 +1058,7 @@ The entries labeled "Characteristic." are characteristic boundary conditions bas
| `bc_[x,y,z]%alpha_rho_in` | Real Array | Inflow density |
| `bc_[x,y,z]%alpha_in` | Real Array | Inflow void fraction |
-This boundary condition can be used for subsonic inflow (`bc_[x,y,z]%[beg,end]` = -7) and subsonic outflow (`bc_[x,y,z]%[beg,end]` = -8) characteristic boundary conditions. These are based on [Pirozzoli (2013)](@ref references). This enables to provide inflow and outflow conditions outside the computational domain.
+This boundary condition can be used for subsonic inflow (`bc_[x,y,z]%[beg,end]` = -7) and subsonic outflow (`bc_[x,y,z]%[beg,end]` = -8) characteristic boundary conditions. These are based on \cite Pirozzoli13. This enables to provide inflow and outflow conditions outside the computational domain.
### Patch types
@@ -1142,7 +1142,7 @@ The midplane location is [`%%loc(1)`, `%%loc(2)`] and the normal vector is [$\ma
- `%%support = 3` specifies a semi-infinite source plane in 3D simulation.
The midplane location is [`%%loc(1)`, `%%loc(2)`] and the normal vector is [$\mathrm{cos}$(`%%dir`), $\mathrm{sin}$(`%%dir`)]. The length of the source plane is `%%length`, and the plane is perpendicular to the direction of wave propagation (defined by `%%dir`). Note that the plane is infinite in the $z$-direction, so `%%loc(3)` is not required.
-- `%%support = 5` specifies a circular transducer in 2D simulation. The transducer is centered at [`%%loc(1)`, `%%loc(2)`] with a focal length of `%%foc_length` and an aperture of `%%aperture`. The center location is not the focal point; it is the tip of the circular arc (intersection of the arc and the x-axis). The aperture is the length of the projection of the circular arc onto the y-axis. If a semi-circle is desired, set the aperture to double the focal length. Note that this is physically a cylindrical transducer, and due to the complexity of Green's function for 2D wave, no closed-form solution is available for the 2D circular transducer, and an approximate is used (see [Maeda and Colonius (2017)](@ref references) for details). For the mass source term correction factor, the theoretical approximation factor of -0.5 in ($r_{foc}^{-0.5}$) is replaced by an empirically determined factor of -0.85.
+- `%%support = 5` specifies a circular transducer in 2D simulation. The transducer is centered at [`%%loc(1)`, `%%loc(2)`] with a focal length of `%%foc_length` and an aperture of `%%aperture`. The center location is not the focal point; it is the tip of the circular arc (intersection of the arc and the x-axis). The aperture is the length of the projection of the circular arc onto the y-axis. If a semi-circle is desired, set the aperture to double the focal length. Note that this is physically a cylindrical transducer, and due to the complexity of Green's function for 2D wave, no closed-form solution is available for the 2D circular transducer, and an approximate is used (see \cite Maeda17 for details). For the mass source term correction factor, the theoretical approximation factor of -0.5 in ($r_{foc}^{-0.5}$) is replaced by an empirically determined factor of -0.85.
- `%%support = 6` specifies a spherical transducer in 2D axisymmetric simulation. It is identical to `%%support = 5` in terms of simulation parameters. Note that this is physically a spherical 3D transducer, so the equation is exact.
diff --git a/docs/documentation/equations.md b/docs/documentation/equations.md
new file mode 100644
index 0000000000..59c6a5c776
--- /dev/null
+++ b/docs/documentation/equations.md
@@ -0,0 +1,847 @@
+@page equations Equations
+
+# MFC: Comprehensive Equations Reference
+
+This document catalogs every equation solved by MFC, organized by physical model.
+Each section notes the input parameter(s) that activate the corresponding physics module and cross-references the relevant source files.
+
+The models and algorithms described here are detailed in \cite Wilfong26 (MFC 5.0) and \cite Bryngelson21. Foundational references for each model are cited inline; see the \ref citelist "Bibliography" for full details.
+
+---
+
+## 1. Overview
+
+MFC solves the compressible Navier-Stokes equations (or Euler equations when viscosity is off) in a finite volume framework. The general semi-discrete form is:
+
+\f[\frac{\partial \mathbf{q}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{q}) + \mathbf{h}(\mathbf{q})\,\nabla \cdot \mathbf{u} = \mathbf{s}(\mathbf{q})\f]
+
+where:
+- \f$\mathbf{q}\f$ is the conservative variable vector,
+- \f$\mathbf{F}\f$ is the flux tensor,
+- \f$\mathbf{h}(\mathbf{q})\,\nabla \cdot \mathbf{u}\f$ contains non-conservative terms (volume fraction advection),
+- \f$\mathbf{s}(\mathbf{q})\f$ is the source vector (bubbles, body forces, chemistry, etc.).
+
+The parameter `model_eqns` (1, 2, 3, or 4) selects the governing equation set.
+
+**Key source files:** `src/simulation/m_rhs.fpp` (RHS evaluation), `src/common/m_variables_conversion.fpp` (EOS and variable conversion).
+
+---
+
+## 2. Governing PDEs
+
+### 2.1 Five-Equation Model (`model_eqns = 2`)
+
+The primary workhorse model (\cite Allaire02; \cite Wilfong26 Sec. 2.1). The state vector is:
+
+\f[\mathbf{q} = \bigl(\alpha_1 \rho_1,\;\alpha_2 \rho_2,\;\ldots,\;\rho u_1,\;\rho u_2,\;\rho u_3,\;\rho E,\;\alpha_1,\;\alpha_2,\;\ldots\bigr)^T\f]
+
+**Continuity** (one per component):
+
+\f[\frac{\partial (\alpha_i \rho_i)}{\partial t} + \nabla \cdot (\alpha_i \rho_i\,\mathbf{u}) = 0\f]
+
+**Momentum:**
+
+\f[\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot \bigl(\rho\,\mathbf{u} \otimes \mathbf{u} + p\,\mathbf{I} - \boldsymbol{\tau}^v\bigr) = 0\f]
+
+**Energy:**
+
+\f[\frac{\partial (\rho E)}{\partial t} + \nabla \cdot \bigl[(\rho E + p)\,\mathbf{u} - \boldsymbol{\tau}^v \cdot \mathbf{u}\bigr] = 0\f]
+
+**Volume fraction advection:**
+
+\f[\frac{\partial \alpha_i}{\partial t} + \mathbf{u} \cdot \nabla \alpha_i = K\,\nabla \cdot \mathbf{u}\f]
+
+where the \f$K\f$ term enforces interface conditions via the Wood sound speed:
+
+\f[K = \frac{\rho_2 c_2^2 - \rho_1 c_1^2}{\displaystyle\frac{\rho_1 c_1^2}{\alpha_1} + \displaystyle\frac{\rho_2 c_2^2}{\alpha_2}}\f]
+
+Setting `alt_soundspeed = .true.` enables the \f$K\f$ correction (\cite Kapila01, with Wood sound speed). Setting `alt_soundspeed = .false.` uses the Allaire variant without the \f$K\f$ correction, which is conservative but does not strictly obey the second law of thermodynamics.
+
+**Mixture rules:**
+
+\f[1 = \sum_i \alpha_i, \qquad \rho = \sum_i \alpha_i \rho_i, \qquad \rho e = \sum_i \alpha_i \rho_i e_i\f]
+
+### 2.2 Six-Equation Model (`model_eqns = 3`)
+
+Allows pressure disequilibrium between phases (\cite Saurel09; \cite Wilfong26 Sec. 2.1).
+
+**Continuity and momentum:** Same as the five-equation model.
+
+**Separate phasic internal energy:**
+
+\f[\frac{\partial (\alpha_i \rho_i e_i)}{\partial t} + \nabla \cdot (\alpha_i \rho_i e_i\,\mathbf{u}) + \alpha_i p_i\,\nabla \cdot \mathbf{u} = -\mu\,p_I\,(p_2 - p_1) - \alpha_i\,\boldsymbol{\tau}_i^v : \nabla \mathbf{u}\f]
+
+**Volume fraction:**
+
+\f[\frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 = \mu\,(p_1 - p_2)\f]
+
+**Interfacial pressure:**
+
+\f[p_I = \frac{z_2\,p_1 + z_1\,p_2}{z_1 + z_2}, \qquad z_i = \rho_i\,c_i\f]
+
+Infinite pressure relaxation is applied at each Runge-Kutta stage to drive toward pressure equilibrium.
+
+**Mixture speed of sound:**
+
+\f[c^2 = \sum_k Y_k\,c_k^2\f]
+
+With phase change (`relax = .true.`), additional source terms appear in the phasic energy and volume fraction equations:
+- **Pressure relaxation:** \f$\mu\,\delta p\f$ where \f$\delta p = p_1 - p_2\f$
+- **Thermal transfer:** \f$Q = \theta\,(T_2 - T_1)\f$
+- **Mass transfer:** \f$\dot{m} = \nu\,(g_2 - g_1)\f$ (Gibbs free energy difference)
+
+See Section 8 (Phase Change) below for details.
+
+### 2.3 Other Model Variants
+
+- `model_eqns = 1`: **Gamma/pi_inf model** — simplified single-fluid formulation using mixture \f$\gamma\f$ and \f$\pi_\infty\f$ directly without tracking individual volume fractions (\cite Johnsen08).
+- `model_eqns = 4`: **Four-equation model** — reduced model from the six-equation system after full pressure-temperature equilibrium relaxation (Tait-like compressible liquid).
+
+---
+
+## 3. Equations of State
+
+### 3.1 Stiffened Gas EOS (\cite Menikoff89; \cite LeMetayer04; \cite Wilfong26 Sec. 2.2)
+
+The primary closure for each phase:
+
+\f[p_k = (\gamma_k - 1)\,\rho_k\,e_k - \gamma_k\,\pi_{\infty,k}\f]
+
+Equivalently:
+
+\f[e_k = \frac{p_k + \gamma_k\,\pi_{\infty,k}}{(\gamma_k - 1)\,\rho_k}\f]
+
+**Total energy relation:**
+
+\f[\rho E = \Gamma\,p + \Pi_\infty + \frac{1}{2}\rho\,|\mathbf{u}|^2 + q_v\f]
+
+where MFC internally tracks the transformed thermodynamic quantities:
+
+\f[\Gamma_k = \frac{1}{\gamma_k - 1}, \qquad \Pi_{\infty,k} = \frac{\gamma_k\,\pi_{\infty,k}}{\gamma_k - 1}\f]
+
+and the mixture rules are arithmetic averages of these transformed quantities:
+
+\f[\Gamma = \sum_i \frac{\alpha_i}{\gamma_i - 1}, \qquad \Pi_\infty = \sum_i \frac{\alpha_i\,\gamma_i\,\pi_{\infty,i}}{\gamma_i - 1}, \qquad q_v = \sum_i \alpha_i\,\rho_i\,q_{v,i}\f]
+
+The pressure is recovered from the total energy as:
+
+\f[p = \frac{\rho E - \frac{1}{2}\rho\,|\mathbf{u}|^2 - \Pi_\infty - q_v}{\Gamma}\f]
+
+**Phasic speed of sound:**
+
+\f[c_k = \sqrt{\frac{\gamma_k\,(p + \pi_{\infty,k})}{\rho_k}}\f]
+
+**Wood mixture sound speed:**
+
+\f[\frac{1}{\rho\,c^2} = \sum_k \frac{\alpha_k}{\rho_k\,c_k^2}\f]
+
+Input parameters per fluid: `gamma` (\f$\gamma_k\f$), `pi_inf` (\f$\pi_{\infty,k}\f$), `cv` (\f$c_{v,k}\f$), `qv` (\f$q_{v,k}\f$), `qvp` (\f$q'_{v,k}\f$).
+
+### 3.2 Ideal Gas EOS (Chemistry, `chemistry = .true.`)
+
+For reacting gas mixtures:
+
+\f[p = \frac{\rho\,R_u\,T}{W}, \qquad W = \left(\sum_m \frac{Y_m}{W_m}\right)^{-1}\f]
+
+Temperature is obtained from the internal energy by Newton iteration:
+
+\f[e_g - \sum_m e_m(T)\,Y_m = 0\f]
+
+**Species internal energy from enthalpy:**
+
+\f[e_m(T) = \frac{\hat{h}_m(T) - R_u\,T}{W_m}\f]
+
+**NASA polynomial enthalpies** (\cite McBride93):
+
+\f[\frac{\hat{h}_m}{R_u\,T} = \frac{C_0}{T} + \sum_{r=1}^{5} \frac{C_r\,T^{r-1}}{r}\f]
+
+---
+
+## 4. Viscous Stress Tensor (`viscous = .true.`)
+
+**Newtonian viscous stress (no bulk viscosity by default):**
+
+\f[\boldsymbol{\tau}^v = 2\,\eta\left(\mathbf{D} - \frac{1}{3}\,\text{tr}(\mathbf{D})\,\mathbf{I}\right)\f]
+
+where the strain rate tensor is:
+
+\f[\mathbf{D} = \frac{1}{2}\bigl(\nabla \mathbf{u} + (\nabla \mathbf{u})^T\bigr)\f]
+
+**With bulk viscosity:**
+
+\f[\tau_{ij} = \mu\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) + \left(\zeta - \frac{2\mu}{3}\right)\delta_{ij}\,\frac{\partial u_k}{\partial x_k}\f]
+
+**Cartesian components:**
+
+\f[\tau_{xx} = \mu\left(2\,\frac{\partial u}{\partial x} - \frac{2}{3}\nabla\cdot\mathbf{u}\right), \qquad \tau_{xy} = \mu\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)\f]
+
+and similarly for all other components. Cylindrical coordinate formulations include additional \f$1/r\f$ terms.
+
+**Viscosity averaging:**
+
+\f[\frac{1}{\text{Re}_\text{mix}} = \sum_j \frac{\alpha_j}{\text{Re}_j}\f]
+
+Input parameters: `Re_inv` (shear and volume Reynolds numbers per fluid).
+
+---
+
+## 5. Cylindrical Coordinates (`cyl_coord = .true.`) (\cite Wilfong26 Sec. 2.3)
+
+Additional geometric source terms appear with \f$1/r\f$ factors in the continuity, momentum, and energy equations. Key modifications:
+
+- **Radial momentum:** extra \f$p/r\f$ and \f$\tau_{\theta\theta}/r\f$ terms
+- **Viscous stress:** \f$\tau_{yy}\f$ includes \f$v/r\f$ corrections:
+
+\f[\tau_{yy} = \mu\left(\frac{4}{3}\frac{\partial v}{\partial r} - \frac{2}{3}\frac{\partial u}{\partial x} - \frac{2}{3}\frac{v}{r}\right)\f]
+
+- **Axis singularity:** axis placed at cell boundary with spectral filtering in the azimuthal direction
+
+---
+
+## 6. Sub-Grid Bubble Dynamics (\cite Wilfong26 Sec. 4.1)
+
+### 6.1 Euler-Euler Bubbles (`bubbles_euler = .true.`)
+
+**Source:** `src/simulation/m_bubbles_EE.fpp`, `src/simulation/m_bubbles.fpp`
+
+#### 6.1.1 Method of Classes (\cite Commander89; \cite Ando11)
+
+**Modified mixture pressure:**
+
+\f[p = (1 - \alpha)\,p_l + \alpha\left(\frac{R^3\,p_{bw}}{\bar{R}^3} - \frac{\rho\,R^3\,\dot{R}^2}{\bar{R}^3}\right)\f]
+
+**Modified stiffened gas for the liquid phase:**
+
+\f[\Gamma_l\,p_l + \Pi_{\infty,l} = \frac{1}{1 - \alpha}\left(E - \frac{1}{2}\rho\,|\mathbf{u}|^2\right)\f]
+
+**Bubble wall pressure (polytropic):**
+
+\f[p_{bw} = \left(p_0 + \frac{2\sigma}{R_0}\right)\left(\frac{R_0}{R}\right)^{3\gamma} - \frac{4\mu\,\dot{R}}{R} - \frac{2\sigma}{R}\f]
+
+**Void fraction transport:**
+
+\f[\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = \frac{3\,\alpha\,\bar{R}^2\,\dot{R}}{\bar{R}^3}\f]
+
+**Number density conservation:**
+
+\f[\frac{\partial n_\text{bub}}{\partial t} + \nabla \cdot (n_\text{bub}\,\mathbf{u}) = 0\f]
+
+where \f$n = \frac{3}{4\pi}\,\frac{\alpha}{\bar{R}^3}\f$.
+
+**Polydispersity** (`polydisperse = .true.`): Log-normal PDF discretized into \f$N_\text{bin}\f$ equilibrium radii with standard deviation `poly_sigma`, integrated via Simpson's rule.
+
+#### 6.1.2 Rayleigh-Plesset (`bubble_model = 3`) (\cite Rayleigh17; \cite Plesset49)
+
+\f[R\,\ddot{R} + \frac{3}{2}\,\dot{R}^2 = \frac{p_{bw} - p_\infty}{\rho_l}\f]
+
+#### 6.1.3 Keller-Miksis (`bubble_model = 2`) (\cite Keller80)
+
+\f[R\,\ddot{R}\left(1 - \frac{\dot{R}}{c}\right) + \frac{3}{2}\,\dot{R}^2\left(1 - \frac{\dot{R}}{3c}\right) = \frac{p_{bw} - p_\infty}{\rho_l}\left(1 + \frac{\dot{R}}{c}\right) + \frac{R\,\dot{p}_{bw}}{\rho_l\,c}\f]
+
+#### 6.1.4 Gilmore (`bubble_model = 1`) (\cite Gilmore52)
+
+Enthalpy-based formulation with compressibility corrections via the Tait EOS:
+
+\f[R\,\ddot{R}\left(1 - \frac{\dot{R}}{C}\right) + \frac{3}{2}\,\dot{R}^2\left(1 - \frac{\dot{R}}{3C}\right) = H\left(1 + \frac{\dot{R}}{C}\right) + \frac{R\,\dot{H}}{C}\left(1 - \frac{\dot{R}}{C}\right)\f]
+
+where the enthalpy difference is:
+
+\f[H = \frac{n_\text{tait}(1 + B)}{n_\text{tait} - 1}\left[\left(\frac{p_{bw}}{1+B} + 1\right)^{(n_\text{tait}-1)/n_\text{tait}} - \left(\frac{p_\infty}{1+B} + 1\right)^{(n_\text{tait}-1)/n_\text{tait}}\right]\f]
+
+and the local liquid sound speed:
+
+\f[C = \sqrt{n_\text{tait}(1+B)\left(\frac{p_\infty}{1+B} + 1\right)^{(n_\text{tait}-1)/n_\text{tait}} + (n_\text{tait} - 1)\,H}\f]
+
+#### 6.1.5 Non-Polytropic Thermal Model (`polytropic = .false.`) (\cite Preston07)
+
+**Internal bubble pressure ODE:**
+
+\f[\dot{p}_b = \frac{3\gamma_b}{R}\left(-\dot{R}\,p_b + R_v\,T_{bw}\,\dot{m}_v + \frac{\gamma_b - 1}{\gamma_b}\,k_{bw}\left.\frac{\partial T}{\partial r}\right|_R\right)\f]
+
+**Vapor mass flux:**
+
+\f[\dot{m}_v = \frac{D\,\rho_{bw}}{1 - \chi_{vw}}\left.\frac{\partial \chi_v}{\partial r}\right|_R\f]
+
+#### 6.1.6 QBMM Moment Transport (`qbmm = .true.`) (\cite Bryngelson20)
+
+**Population balance equation:**
+
+\f[\frac{\partial f}{\partial t} + \frac{\partial (f\,\dot{R})}{\partial R} + \frac{\partial (f\,\ddot{R})}{\partial \dot{R}} = 0\f]
+
+**Moment transport:**
+
+\f[\frac{\partial (n_\text{bub}\,\mu_i)}{\partial t} + \nabla \cdot (n_\text{bub}\,\mu_i\,\mathbf{u}) = n_\text{bub}\,\dot{\mu}_i\f]
+
+where moments \f$\mu_{i_1,i_2} = \int R^{i_1}\,\dot{R}^{i_2}\,f\,dR\,d\dot{R}\f$.
+
+**CHyQMOM inversion** recovers 4 quadrature nodes \f$(w_j, R_j, \dot{R}_j)\f$ from 6 moments via:
+
+\f[\bar{u} = \frac{\mu_{10}}{\mu_{00}}, \quad \bar{v} = \frac{\mu_{01}}{\mu_{00}}, \quad c_{20} = \frac{\mu_{20}}{\mu_{00}} - \bar{u}^2, \quad c_{11} = \frac{\mu_{11}}{\mu_{00}} - \bar{u}\bar{v}, \quad c_{02} = \frac{\mu_{02}}{\mu_{00}} - \bar{v}^2\f]
+
+### 6.2 Euler-Lagrange Bubbles (`bubbles_lagrange = .true.`) (\cite Maeda18)
+
+**Source:** `src/simulation/m_bubbles_EL.fpp`
+
+Volume-averaged carrier flow equations with bubble source terms:
+
+**Continuity:**
+
+\f[\frac{\partial \rho_l}{\partial t} + \nabla \cdot (\rho_l\,\mathbf{u}_l) = \frac{\rho_l}{1 - \alpha}\left[\frac{\partial \alpha}{\partial t} + \mathbf{u}_l \cdot \nabla \alpha\right]\f]
+
+**Momentum:**
+
+\f[\frac{\partial (\rho_l\,\mathbf{u}_l)}{\partial t} + \nabla \cdot (\rho_l\,\mathbf{u}_l \otimes \mathbf{u}_l + p\,\mathbf{I} - \boldsymbol{\tau}_l) = \frac{\rho_l\,\mathbf{u}_l}{1 - \alpha}\left[\frac{\partial \alpha}{\partial t} + \mathbf{u}_l \cdot \nabla \alpha\right] - \frac{\alpha}{1 - \alpha}\,\nabla \cdot (p\,\mathbf{I} - \boldsymbol{\tau}_l)\f]
+
+**Energy:**
+
+\f[\frac{\partial E_l}{\partial t} + \nabla \cdot \bigl[(E_l + p)\,\mathbf{u}_l - \boldsymbol{\tau}_l \cdot \mathbf{u}_l\bigr] = \frac{E_l}{1 - \alpha}\left[\frac{\partial \alpha}{\partial t} + \mathbf{u}_l \cdot \nabla \alpha\right] - \frac{\alpha}{1 - \alpha}\,\nabla \cdot (p\,\mathbf{u}_l - \boldsymbol{\tau}_l \cdot \mathbf{u}_l)\f]
+
+The left-hand side is the standard conservation law for the liquid phase; the right-hand side source terms capture the effect of the bubbles on the host liquid.
+
+**Void fraction via regularization kernel:**
+
+\f[\alpha(\mathbf{x}) = \sum_n V_n\,\delta_\sigma(\mathbf{x} - \mathbf{x}_n)\f]
+
+where \f$\delta_\sigma\f$ is a Gaussian kernel:
+
+\f[\delta_\sigma(\mathbf{r}) = \frac{1}{(2\pi\sigma^2)^{3/2}}\exp\!\left(-\frac{|\mathbf{r}|^2}{2\sigma^2}\right)\f]
+
+with \f$\sigma = \varepsilon_b \max(\Delta x^{1/3}_\text{cell},\;R_\text{bubble})\f$.
+
+Each bubble is tracked individually with Keller-Miksis dynamics and 4th-order adaptive Runge-Kutta time integration.
+
+---
+
+## 7. Fluid-Structure Interaction
+
+### 7.1 Hypoelastic Model (`hypoelasticity = .true.`) (\cite Rodriguez19; \cite Wilfong26 Sec. 4.1.6)
+
+**Source:** `src/simulation/m_hypoelastic.fpp`
+
+**Cauchy stress decomposition:**
+
+\f[\sigma_{ij} = -p\,\delta_{ij} + \tau_{ij}^{(v)} + \tau_{ij}^{(e)}\f]
+
+**Elastic energy contribution to total energy:**
+
+\f[E = e + \frac{|\mathbf{u}|^2}{2} + \frac{\boldsymbol{\tau}^e : \boldsymbol{\tau}^e}{4\,\rho\,G}\f]
+
+**Elastic stress evolution:**
+
+\f[\frac{\partial (\rho\,\boldsymbol{\tau}^e)}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\tau}^e \otimes \mathbf{u}) = \mathbf{S}^e\f]
+
+**Source term:**
+
+\f[\mathbf{S}^e = \rho\bigl(\mathbf{l} \cdot \boldsymbol{\tau}^e + \boldsymbol{\tau}^e \cdot \mathbf{l}^T - \boldsymbol{\tau}^e\,\text{tr}(\mathbf{D}) + 2G\,\mathbf{D}^d\bigr)\f]
+
+where \f$\mathbf{l} = \nabla \mathbf{u}\f$ is the velocity gradient and \f$\mathbf{D}^d = \mathbf{D} - \frac{1}{3}\text{tr}(\mathbf{D})\,\mathbf{I}\f$ is the deviatoric strain rate.
+
+**Lie objective temporal derivative (Kelvin-Voigt):**
+
+\f[\hat{\boldsymbol{\tau}}^e = \frac{D\boldsymbol{\tau}^e}{Dt} - \mathbf{l} \cdot \boldsymbol{\tau}^e - \boldsymbol{\tau}^e \cdot \mathbf{l}^T + \boldsymbol{\tau}^e\,\text{tr}(\mathbf{D}) = 2G\,\mathbf{D}^d\f]
+
+This adds 6 additional transport equations in 3D (symmetric stress tensor: \f$\tau_{xx}^e, \tau_{xy}^e, \tau_{yy}^e, \tau_{xz}^e, \tau_{yz}^e, \tau_{zz}^e\f$).
+
+### 7.2 Hyperelastic Model (`hyperelasticity = .true.`) (\cite Kamrin12; \cite Wilfong26 Sec. 4.1.6)
+
+**Source:** `src/simulation/m_hyperelastic.fpp`
+
+**Reference map evolution:**
+
+\f[\frac{\partial (\rho\,\boldsymbol{\xi})}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\xi} \otimes \mathbf{u}) = 0\f]
+
+**Deformation gradient from reference map:**
+
+\f[\mathbf{F} = (\nabla \boldsymbol{\xi})^{-1}\f]
+
+**Left Cauchy-Green tensor:**
+
+\f[\mathbf{b} = \mathbf{F}\,\mathbf{F}^T\f]
+
+**Neo-Hookean Cauchy stress:**
+
+\f[\boldsymbol{\tau}^e = \frac{G}{J}\left(\mathbf{b} - \frac{\text{tr}(\mathbf{b})}{3}\,\mathbf{I}\right)\f]
+
+where \f$J = \det(\mathbf{F})\f$.
+
+**Hyperelastic energy:**
+
+\f[e^e = \frac{G}{2}\bigl(I_{\mathbf{b}} - 3\bigr), \qquad I_{\mathbf{b}} = \text{tr}(\mathbf{b})\f]
+
+---
+
+## 8. Phase Change (`relax = .true.`) (\cite Wilfong26 Sec. 4.1.3)
+
+**Source:** `src/common/m_phase_change.fpp`
+
+### 8.1 pT-Relaxation (`relax_model = 5`) (\cite Saurel08)
+
+\f$N\f$-fluid pressure-temperature equilibrium. The equilibrium condition is:
+
+\f[f(p) = \sum_i \alpha_i - 1 = 0\f]
+
+**Temperature from energy conservation:**
+
+\f[T = \frac{\rho e + p - \sum_i (\alpha_i \rho_i)\,q_{v,i}}{\sum_i (\alpha_i \rho_i)\,c_{v,i}\,\gamma_i}\f]
+
+**Newton residual:**
+
+\f[g(p) = \sum_i \frac{(\gamma_i - 1)\,(\alpha_i \rho_i)\,c_{v,i}}{(p + \pi_{\infty,i})} \cdot \frac{\rho e + p - \sum_j (\alpha_j \rho_j)\,q_{v,j}}{\sum_j (\alpha_j \rho_j)\,c_{v,j}\,\gamma_j}\f]
+
+Solved via Newton's method for the equilibrium pressure.
+
+### 8.2 pTg-Relaxation (`relax_model = 6`) (\cite Zein10)
+
+Two coupled equations for \f$(\alpha_1 \rho_1,\;p)\f$:
+
+**Gibbs free energy equilibrium (Clausius-Clapeyron):**
+
+\f[F_1 = T\left[(c_{v,l}\gamma_l - c_{v,v}\gamma_v)(1 - \ln T) - (q'_l - q'_v) + c_{v,l}(\gamma_l - 1)\ln(p + \pi_{\infty,l}) - c_{v,v}(\gamma_v - 1)\ln(p + \pi_{\infty,v})\right] + q_{v,l} - q_{v,v} = 0\f]
+
+**Energy conservation constraint:**
+
+\f[F_2 = \rho e + p + m_l\,(q_{v,v} - q_{v,l}) - m_T\,q_{v,v} - m_{qD} + \frac{m_l\,(c_{v,v}\,\gamma_v - c_{v,l}\,\gamma_l) - m_T\,c_{v,v}\,\gamma_v - m_{cpD}}{m_l\left(\frac{c_{v,l}\,(\gamma_l - 1)}{p + \pi_{\infty,l}} - \frac{c_{v,v}\,(\gamma_v - 1)}{p + \pi_{\infty,v}}\right) + \frac{m_T\,c_{v,v}\,(\gamma_v - 1)}{p + \pi_{\infty,v}} + m_{cvgp}} = 0\f]
+
+where \f$m_T\f$ is the total mass, \f$m_l = \alpha_l \rho_l\f$ is the liquid partial density, and \f$m_{qD}\f$, \f$m_{cpD}\f$, \f$m_{cvgp}\f$ are auxiliary thermodynamic sums over additional fluids (beyond the phase-changing pair).
+
+Solved via 2D Newton-Raphson.
+
+---
+
+## 9. Chemistry and Combustion (`chemistry = .true.`) (\cite Wilfong26 Sec. 4.1.7)
+
+**Source:** `src/common/m_chemistry.fpp`
+
+**Species transport:**
+
+\f[\frac{\partial (\rho_g\,Y_m)}{\partial t} + \frac{\partial (\rho_g\,u_i\,Y_m)}{\partial x_i} = W_m\,\dot{\omega}_m\f]
+
+**Net production rate:**
+
+\f[\dot{\omega}_m = \sum_n (\nu''_{mn} - \nu'_{mn})\,\mathcal{R}_n\f]
+
+**Reaction rate (law of mass action):**
+
+\f[\mathcal{R}_n = k_n(T)\left[\prod_j \left(\frac{\rho_g\,Y_j}{W_j}\right)^{\nu'_{jn}} - \frac{1}{K_n}\prod_k \left(\frac{\rho_g\,Y_k}{W_k}\right)^{\nu''_{kn}}\right]\f]
+
+**Arrhenius rate:**
+
+\f[k_n(T) = A_n\,T^{b_n}\exp\!\left(-\frac{T_{a,n}}{T}\right)\f]
+
+**Molecular diffusion** (`transport_model`):
+- **Mixture-average:** Species-specific diffusion coefficients \f$D_m^\text{mix}\f$, mass flux: \f$\dot{m}_k = \rho\,D_k^\text{mix}\,(W_k / W_\text{mix})\,\partial X_k / \partial x\f$
+- **Unity Lewis number:** \f$D_m = \lambda / (\rho\,c_p)\f$
+
+Enthalpy flux with diffusion:
+
+\f[q_\text{diff} = \lambda\,\frac{\partial T}{\partial x} + \sum_k h_k\,\dot{m}_k\f]
+
+Reaction mechanisms are code-generated via Pyrometheus (\cite Cisneros26), which provides symbolic abstractions for thermochemistry that enable portable GPU computation and automatic differentiation of chemical source terms.
+
+---
+
+## 10. Surface Tension (`surface_tension = .true.`) (\cite Schmidmayer17; \cite Wilfong26 Sec. 4.1.8)
+
+**Source:** `src/simulation/m_surface_tension.fpp`, `src/simulation/include/inline_capillary.fpp`
+
+**Color function advection:**
+
+\f[\frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla c = 0\f]
+
+**Capillary stress tensor (CSF model):**
+
+\f[\boldsymbol{\Omega} = -\sigma\left(\|\nabla c\|\,\mathbf{I} - \frac{\nabla c \otimes \nabla c}{\|\nabla c\|}\right)\f]
+
+In component form, with \f$\hat{w}_i = (\partial c / \partial x_i) / \|\nabla c\|\f$:
+
+\f[\Omega_{xx} = -\sigma\,(\hat{w}_y^2 + \hat{w}_z^2)\,\|\nabla c\|, \qquad \Omega_{xy} = \sigma\,\hat{w}_x\,\hat{w}_y\,\|\nabla c\|\f]
+
+The capillary stress divergence is added to the momentum and energy equations. The total energy equation becomes:
+
+\f[\frac{\partial (\rho E + \varepsilon_0)}{\partial t} + \nabla \cdot \bigl[(\rho E + \varepsilon_0 + p)\,\mathbf{u} + (\boldsymbol{\Omega} - \boldsymbol{\tau}^v) \cdot \mathbf{u}\bigr] = 0\f]
+
+**Capillary mixture energy:**
+
+\f[\varepsilon_0 = \sigma\,\|\nabla c\|\f]
+
+---
+
+## 11. Magnetohydrodynamics
+
+### 11.1 Ideal MHD (`mhd = .true.`) (\cite Wilfong26 Sec. 4.1.9)
+
+**Continuity:**
+
+\f[\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho\,\mathbf{u}) = 0\f]
+
+**Momentum:**
+
+\f[\frac{\partial (\rho\,\mathbf{u})}{\partial t} + \nabla \cdot \left[\rho\,\mathbf{u} \otimes \mathbf{u} + \left(p + \frac{|\mathbf{B}|^2}{2}\right)\mathbf{I} - \mathbf{B} \otimes \mathbf{B}\right] = 0\f]
+
+**Energy:**
+
+\f[\frac{\partial \mathcal{E}}{\partial t} + \nabla \cdot \left[\left(\mathcal{E} + p + \frac{|\mathbf{B}|^2}{2}\right)\mathbf{u} - (\mathbf{u} \cdot \mathbf{B})\,\mathbf{B}\right] = 0\f]
+
+**Induction:**
+
+\f[\frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u}) = 0\f]
+
+**Total energy:**
+
+\f[\mathcal{E} = \rho\,e + \frac{1}{2}\rho\,|\mathbf{u}|^2 + \frac{|\mathbf{B}|^2}{2}\f]
+
+**Fast magnetosonic speed:**
+
+\f[c_f = \sqrt{\frac{1}{2}\left(c_s^2 + v_A^2 + \sqrt{(c_s^2 + v_A^2)^2 - 4\,c_s^2\,v_A^2\cos^2\theta}\right)}\f]
+
+**Alfven speed:**
+
+\f[v_A = \sqrt{\frac{|\mathbf{B}|^2}{\rho}}\f]
+
+Uses the HLLD Riemann solver (`riemann_solver = 4`). Hyperbolic divergence cleaning (`hyper_cleaning = .true.`) via the GLM method (\cite Dedner02).
+
+### 11.2 Relativistic MHD (`relativity = .true.`) (\cite Wilfong26 Sec. 4.1.10)
+
+**Conserved variables:**
+
+\f[\mathbf{U} = (D,\;\mathbf{m},\;\tau,\;\mathbf{B})^T\f]
+
+where:
+
+\f[D = \Gamma\,\rho, \qquad \mathbf{m} = \Gamma^2\rho h\,\mathbf{u} + |\mathbf{B}|^2\mathbf{u} - (\mathbf{u} \cdot \mathbf{B})\,\mathbf{B}\f]
+
+\f[\tau = \Gamma^2\rho h - p + \frac{|\mathbf{B}|^2}{2} + \frac{|\mathbf{u}|^2|\mathbf{B}|^2 - (\mathbf{B} \cdot \mathbf{u})^2}{2} - \Gamma\,\rho\f]
+
+Primitive recovery uses Newton-Raphson on the nonlinear conserved-to-primitive relation.
+
+---
+
+## 12. Information Geometric Regularization (`igr = .true.`) (\cite Wilfong25a)
+
+**Source:** `src/simulation/m_igr.fpp`
+
+**Modified momentum with entropic pressure** \f$\Sigma\f$**:**
+
+\f[\frac{\partial (\rho\,\mathbf{u})}{\partial t} + \nabla \cdot \bigl[\rho\,\mathbf{u} \otimes \mathbf{u} + (p + \Sigma)\,\mathbf{I} - \boldsymbol{\tau}\bigr] = 0\f]
+
+**Elliptic PDE for** \f$\Sigma\f$**:**
+
+\f[\alpha\left[\text{tr}(\nabla \mathbf{u})^2 + \text{tr}^2(\nabla \mathbf{u})\right] = \frac{\Sigma}{\rho} - \alpha\,\nabla \cdot \left(\frac{\nabla \Sigma}{\rho}\right)\f]
+
+where \f$\alpha \sim \Delta x^2\f$ (regularization strength proportional to mesh spacing squared):
+
+\f[\alpha_\text{IGR} = \alpha_\text{factor} \cdot \max(\Delta x,\;\Delta y,\;\Delta z)^2\f]
+
+**RHS strain-rate source (3D):**
+
+\f[\text{RHS} = \alpha\left[2\left(\frac{\partial u}{\partial y}\frac{\partial v}{\partial x} + \frac{\partial u}{\partial z}\frac{\partial w}{\partial x} + \frac{\partial v}{\partial z}\frac{\partial w}{\partial y}\right) + \left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial w}{\partial z}\right)^2 + (\nabla \cdot \mathbf{u})^2\right]\f]
+
+**Iterative solver:** Jacobi (`igr_iter_solver = 1`) or Gauss-Seidel (`igr_iter_solver = 2`), up to `num_igr_iters` iterations (default 5).
+
+Uses Lax-Friedrichs flux (replaces WENO + Riemann solver).
+
+---
+
+## 13. Body Forces (`bf_x`, `bf_y`, `bf_z`)
+
+**Source:** `src/simulation/m_body_forces.fpp`
+
+**Time-dependent acceleration:**
+
+\f[a_i(t) = g_i + k_i\sin(\omega_i\,t - \phi_i)\f]
+
+**Momentum source:**
+
+\f[\frac{\partial (\rho\,u_i)}{\partial t}\bigg|_\text{bf} = \rho\,a_i(t)\f]
+
+**Energy source:**
+
+\f[\frac{\partial E}{\partial t}\bigg|_\text{bf} = \rho\,\mathbf{u} \cdot \mathbf{a}(t)\f]
+
+---
+
+## 14. Acoustic Sources (`acoustic_source = .true.`)
+
+**Source:** `src/simulation/m_acoustic_src.fpp`
+
+Source terms added to the RHS of the governing equations.
+
+Formulation follows \cite Maeda17.
+
+**Discrete delta function (spatial support):**
+
+\f[\delta_h(r) = \frac{1}{(2\pi\sigma^2)^{d/2}}\exp\!\left(-\frac{r^2}{2\sigma^2}\right)\f]
+
+**Forcing form (added to conservative variables):**
+
+\f[\mathbf{s}_\text{ac} = \Omega_\Gamma\,f(t)\left[\frac{1}{c_0},\;\cos\theta,\;\sin\theta,\;\frac{c_0^2}{\gamma - 1}\right]^T\f]
+
+**Temporal profiles:**
+- **Sine** (`pulse = 1`): \f$S(t) = M\sin(\omega(t - t_\text{delay}))\f$
+- **Gaussian** (`pulse = 2`): \f$S(t) = M\exp\!\bigl(-\frac{1}{2}((t - t_\text{delay})/\sigma_t)^2\bigr)\f$
+- **Square** (`pulse = 3`): \f$S(t) = M\,\text{sign}(\sin(\omega(t - t_\text{delay})))\f$
+- **Broadband** (`pulse = 4`): superposition of multiple frequencies across a bandwidth
+
+**Spatial supports:** planar, spherical transducer, cylindrical transducer, transducer array (arcuate, annular, circular).
+
+---
+
+## 15. Numerical Methods
+
+### 15.1 Spatial Reconstruction
+
+#### WENO (`weno_order = 3, 5, 7`)
+
+**Source:** `src/simulation/m_weno.fpp`
+
+Weighted sum of candidate polynomials at cell interfaces:
+
+\f[f_{i+1/2} = \sum_r \omega_r\,f_{i+1/2}^{(r)}\f]
+
+**WENO-JS** (\cite Jiang96, default):
+
+\f[\alpha_r = \frac{d_r}{(\beta_r + \varepsilon)^2}, \qquad \omega_r = \frac{\alpha_r}{\sum_s \alpha_s}\f]
+
+where \f$d_r\f$ are ideal weights, \f$\beta_r\f$ are smoothness indicators, and \f$\varepsilon\f$ is a small regularization parameter (`weno_eps`).
+
+**WENO-M** (`mapped_weno = .true.`): \cite Henrick05 mapped weights for improved accuracy at critical points:
+
+\f[\omega_M^{(r)} = \frac{d_r\bigl(1 + d_r - 3\omega_0^{(r)} + (\omega_0^{(r)})^2\bigr)\,\omega_0^{(r)}}{d_r^2 + \omega_0^{(r)}(1 - 2d_r)}, \qquad \omega^{(r)} = \frac{\omega_M^{(r)}}{\sum_s \omega_M^{(s)}}\f]
+
+**WENO-Z** (`wenoz = .true.`): \cite Borges08 improved weights with global smoothness measure:
+
+\f[\alpha_r = d_r\left(1 + \left(\frac{\tau}{\beta_r + \varepsilon}\right)^q\right), \qquad \tau = |\beta_0 - \beta_{k-1}|\f]
+
+The parameter \f$q\f$ controls the convergence rate at critical points (typically \f$q = 1\f$ for fifth-order reconstruction, as used in MFC).
+
+**TENO** (`teno = .true.`): \cite Fu16 targeted ENO with smoothness threshold \f$C_T\f$ (`teno_CT`):
+
+\f[\gamma_r = 1 + \frac{\tau}{\beta_r}, \qquad \xi_r = \frac{\gamma_r}{\sum_s \gamma_s}\f]
+
+If \f$\xi_r < C_T\f$, set \f$\alpha_r = 0\f$ (stencil excluded).
+
+Primitive variable reconstruction is used to avoid spurious oscillations at interfaces.
+
+#### MUSCL (`muscl_order = 2`)
+
+**Source:** `src/simulation/m_muscl.fpp`
+
+\f[q_L(j) = q(j) - \frac{1}{2}\,\phi(r)\,\Delta q, \qquad q_R(j) = q(j) + \frac{1}{2}\,\phi(r)\,\Delta q\f]
+
+**Five slope limiters** (`muscl_lim`):
+1. **Minmod:** \f$\phi = \text{sign}(a)\,\min(|a|,\;|b|)\f$ if \f$ab > 0\f$, else \f$0\f$
+2. **MC (Monotonized Central):** \f$\phi = \text{sign}(a)\,\min(2|a|,\;2|b|,\;\frac{1}{2}(|a|+|b|))\f$ if \f$ab > 0\f$
+3. **OSPRE (Van Albada):** \f$\phi = (a^2 b + a b^2)/(a^2 + b^2)\f$
+4. **Van Leer:** \f$\phi = 2ab/(a + b)\f$ if \f$ab > 0\f$
+5. **Superbee:** \f$\phi = \max(\min(2|a|,|b|),\;\min(|a|,2|b|))\f$ if \f$ab > 0\f$
+
+where \f$a\f$ and \f$b\f$ are the left and right slope differences.
+
+**THINC interface compression** (`int_comp = .true.`): applies hyperbolic tangent compression near material interfaces:
+
+\f[q_\text{THINC} = q_\min + \frac{q_\max}{2}\left(1 + \text{sign}(s)\,\frac{\tanh(\beta) + A}{1 + A\,\tanh(\beta)}\right)\f]
+
+where \f$A = \frac{\exp(\text{sign}(s)\,\beta\,(2C - 1))\,/\,\cosh(\beta) - 1}{\tanh(\beta)}\f$ and \f$\beta\f$ controls compression steepness.
+
+#### IGR Reconstruction
+
+5th-order or 3rd-order polynomial interpolation without WENO nonlinear weights, using Lax-Friedrichs numerical flux. Stencil coefficients:
+
+**5th-order:** \f$\text{coeff}_L = [-3/60,\;27/60,\;47/60,\;-13/60,\;2/60]\f$
+
+**3rd-order:** \f$\text{coeff}_L = [2/6,\;5/6,\;-1/6]\f$
+
+### 15.2 Riemann Solvers
+
+**Source:** `src/simulation/m_riemann_solvers.fpp`
+
+#### HLL (`riemann_solver = 1`) (\cite Harten83)
+
+\f[\mathbf{F}^\text{HLL} = \frac{S_R\,\mathbf{F}_L - S_L\,\mathbf{F}_R + S_L\,S_R\,(\mathbf{q}_R - \mathbf{q}_L)}{S_R - S_L}\f]
+
+with wave speed estimates \f$S_L = \min(u_L - c_L,\;u_R - c_R)\f$, \f$S_R = \max(u_R + c_R,\;u_L + c_L)\f$.
+
+#### HLLC (`riemann_solver = 2`) (\cite Toro94)
+
+Four-state solver resolving the contact discontinuity. Star-state satisfies:
+
+\f[p_L^* = p_R^* = p^*, \qquad u_L^* = u_R^* = u^*\f]
+
+#### Exact (`riemann_solver = 3`)
+
+Iterative exact Riemann solver.
+
+#### HLLD (`riemann_solver = 4`, MHD only)
+
+Seven-state solver for ideal MHD resolving fast magnetosonic, Alfven, and contact waves (\cite Miyoshi05). The Riemann fan is divided by outer wave speeds \f$S_L\f$, \f$S_R\f$, Alfven speeds \f$S_L^*\f$, \f$S_R^*\f$, and a middle contact \f$S_M\f$:
+
+\f[S_M = \frac{(S_R - u_R)\rho_R u_R - (S_L - u_L)\rho_L u_L - p_{T,R} + p_{T,L}}{(S_R - u_R)\rho_R - (S_L - u_L)\rho_L}\f]
+
+\f[S_L^* = S_M - \frac{|B_x|}{\sqrt{\rho_L^*}}, \qquad S_R^* = S_M + \frac{|B_x|}{\sqrt{\rho_R^*}}\f]
+
+where \f$p_T = p + |\mathbf{B}|^2/2\f$ is the total (thermal + magnetic) pressure. Continuity of normal velocity and total pressure is enforced across the Riemann fan.
+
+### 15.3 Time Integration
+
+**Source:** `src/simulation/m_time_steppers.fpp`
+
+#### TVD Runge-Kutta (`time_stepper = 1, 2, 3`) (\cite Gottlieb98)
+
+**RK1 (Forward Euler):**
+
+\f[\mathbf{q}^{n+1} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)\f]
+
+**RK2:**
+
+\f[\mathbf{q}^{(1)} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)\f]
+
+\f[\mathbf{q}^{n+1} = \frac{1}{2}\mathbf{q}^n + \frac{1}{2}\mathbf{q}^{(1)} + \frac{1}{2}\Delta t\,\mathbf{L}(\mathbf{q}^{(1)})\f]
+
+**RK3 (SSP):**
+
+\f[\mathbf{q}^{(1)} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)\f]
+
+\f[\mathbf{q}^{(2)} = \frac{3}{4}\mathbf{q}^n + \frac{1}{4}\mathbf{q}^{(1)} + \frac{1}{4}\Delta t\,\mathbf{L}(\mathbf{q}^{(1)})\f]
+
+\f[\mathbf{q}^{n+1} = \frac{1}{3}\mathbf{q}^n + \frac{2}{3}\mathbf{q}^{(2)} + \frac{2}{3}\Delta t\,\mathbf{L}(\mathbf{q}^{(2)})\f]
+
+#### Adaptive Time Stepping (`adap_dt = .true.`)
+
+Embedded RK pairs for error estimation with Hairer-Norsett-Wanner algorithm for initial step size.
+
+#### Strang Splitting (\cite Strang68) (for stiff bubble dynamics)
+
+For equations of the form \f$\partial \mathbf{q}/\partial t = -\nabla \cdot \mathbf{F}(\mathbf{q}) + \mathbf{s}(\mathbf{q})\f$, the Strang splitting scheme integrates three sub-equations per time step:
+
+\f[\frac{\partial \mathbf{q}'}{\partial t} = \mathbf{s}(\mathbf{q}'), \quad t \in [0,\,\Delta t/2], \quad \mathbf{q}'(0) = \mathbf{q}^n\f]
+
+\f[\frac{\partial \mathbf{q}''}{\partial t} = -\nabla \cdot \mathbf{F}(\mathbf{q}''), \quad t \in [0,\,\Delta t], \quad \mathbf{q}''(0) = \mathbf{q}'(\Delta t/2)\f]
+
+\f[\frac{\partial \mathbf{q}'''}{\partial t} = \mathbf{s}(\mathbf{q}'''), \quad t \in [0,\,\Delta t/2], \quad \mathbf{q}^{n+1} = \mathbf{q}'''(\Delta t/2)\f]
+
+The stiff bubble source terms are integrated using an adaptive embedded Runge-Kutta scheme with error control.
+
+### 15.4 CFL Condition
+
+**Inviscid:**
+
+\f[\Delta t = \text{CFL} \cdot \min_{i,j,k}\left(\frac{\Delta x_\text{cell}}{|u| + |v| + |w| + c}\right)\f]
+
+**Viscous:**
+
+\f[\Delta t_v = \text{CFL} \cdot \min\left(\frac{\Delta x^2}{\nu}\right)\f]
+
+where \f$c\f$ is the mixture sound speed and \f$\nu\f$ is the kinematic viscosity.
+
+### 15.5 Finite Differences (`fd_order = 1, 2, 4`)
+
+**Source:** `src/common/m_finite_differences.fpp`
+
+Used for viscous fluxes and velocity gradients.
+
+**1st-order:**
+
+\f[f'_i = \frac{f_{i+1} - f_i}{x_{i+1} - x_i}\f]
+
+**2nd-order centered:**
+
+\f[f'_i = \frac{f_{i+1} - f_{i-1}}{x_{i+1} - x_{i-1}}\f]
+
+**4th-order centered:**
+
+\f[f'_i = \frac{-f_{i+2} + 8f_{i+1} - 8f_{i-1} + f_{i-2}}{-x_{i+2} + 8x_{i+1} - 8x_{i-1} + x_{i-2}}\f]
+
+**Boundary-adjusted one-sided stencils (3rd-order):**
+
+\f[f'_i\big|_\text{left} = \frac{-3f_i + 4f_{i+1} - f_{i+2}}{x_{i+2} - x_i}\f]
+
+\f[f'_i\big|_\text{right} = \frac{3f_i - 4f_{i-1} + f_{i-2}}{x_i - x_{i-2}}\f]
+
+---
+
+## 16. Boundary Conditions
+
+**Source:** `src/simulation/m_cbc.fpp`, `src/simulation/m_compute_cbc.fpp`, `src/common/m_constants.fpp`
+
+### 16.1 Simple BCs
+
+| BC Code | Type |
+|---------|------|
+| `-1` | Periodic |
+| `-2` | Reflective |
+| `-3` | Ghost cell extrapolation |
+| `-4` | Riemann extrapolation |
+| `-14` | Axis symmetry |
+| `-15` | Slip wall |
+| `-16` | No-slip wall |
+
+### 16.2 Characteristic BCs (\cite Thompson87, \cite Thompson90; `bc_x%beg = -5` to `-12`)
+
+**Characteristic decomposition:**
+
+\f[\frac{\partial \mathbf{q}_c}{\partial t} + \mathbf{R}_x\,\boldsymbol{\Lambda}_x\,\mathbf{L}_x\,\frac{\partial \mathbf{q}_c}{\partial x} = 0\f]
+
+**Characteristic amplitudes** \f$\mathcal{L}_1\f$ through \f$\mathcal{L}_5\f$ define wave interactions at boundaries:
+
+\f[\mathcal{L}_1 = (u - c)\left(\frac{\partial p}{\partial x} - \rho\,c\,\frac{\partial u}{\partial x}\right), \qquad \mathcal{L}_2 = u\left(c^2\,\frac{\partial \rho}{\partial x} - \frac{\partial p}{\partial x}\right)\f]
+
+\f[\mathcal{L}_3 = u\,\frac{\partial v}{\partial x}, \qquad \mathcal{L}_4 = u\,\frac{\partial w}{\partial x}, \qquad \mathcal{L}_5 = (u + c)\left(\frac{\partial p}{\partial x} + \rho\,c\,\frac{\partial u}{\partial x}\right)\f]
+
+For non-reflecting boundaries, the incoming wave amplitudes are set to zero.
+
+**GRCBC (incoming wave from ghost point):**
+
+\f[\mathcal{L} = -\boldsymbol{\Lambda}^o\,\mathbf{L}_x\,\frac{\partial \mathbf{q}_c}{\partial x} + n_x\,\lambda^i\,\mathbf{L}_x\,\frac{\mathbf{P}\,(\mathbf{q}_p^{(g)} - \mathbf{q}_p)}{\Delta}\f]
+
+**8 types:** slip wall (`-5`), non-reflecting buffer (`-6`), non-reflecting sub. inflow (`-7`), non-reflecting sub. outflow (`-8`), force-free sub. outflow (`-9`), constant-pressure sub. outflow (`-10`), supersonic inflow (`-11`), supersonic outflow (`-12`).
+
+### 16.3 Immersed Boundary Method (`ib = .true.`) (\cite Tseng03; \cite Mittal05; \cite Wilfong26 Sec. 4.1.1)
+
+**Source:** `src/simulation/m_ibm.fpp`
+
+Ghost-cell IBM with level set field \f$\phi\f$.
+
+**Image point:**
+
+\f[\mathbf{x}_{ip} = \mathbf{x}_{gp} + 2\,\phi(\mathbf{x}_{gp})\,\hat{\mathbf{n}}\f]
+
+**Velocity boundary conditions:**
+- **No-slip:** \f$\mathbf{u}_{gp} = \mathbf{0}\f$
+- **Slip:** \f$\mathbf{u}_{gp} = \mathbf{u}_{ip} - (\hat{\mathbf{n}} \cdot \mathbf{u}_{ip})\,\hat{\mathbf{n}}\f$
+
+**Neumann BC** for pressure, density, and volume fractions.
+
+Supports STL/OBJ geometry import with ray tracing for inside/outside determination.
+
+---
+
+## 17. Low Mach Number Corrections (\cite Wilfong26 Sec. 4.2.4)
+
+**Chen correction** (`low_Mach = 1`, \cite Chen22): anti-dissipation pressure correction (APC) added to the HLLC flux:
+
+\f[p_d = \frac{\rho_L\,\rho_R\,(S_R - u_R)(S_L - u_L)}{\rho_R(S_R - u_R) - \rho_L(S_L - u_L)}\f]
+
+\f[\text{APC}_{p\mathbf{u}} = (z - 1)\,p_d\,\hat{\mathbf{n}}, \qquad \text{APC}_{pE} = (z - 1)\,p_d\,S_*\f]
+
+The corrected HLLC flux in the star region becomes \f$\mathbf{F}^{*} + \text{APC}\f$.
+
+**Thornber correction** (`low_Mach = 2`, \cite Thornber08): modifies the reconstructed velocities at cell interfaces:
+
+\f[u'_L = \frac{u_L + u_R}{2} + z\,\frac{u_L - u_R}{2}, \qquad u'_R = \frac{u_L + u_R}{2} - z\,\frac{u_L - u_R}{2}\f]
+
+\f[z = \min\!\left(\max\!\left(\frac{|u_L|}{c_L},\;\frac{|u_R|}{c_R}\right),\;1\right)\f]
+
+This reduces numerical dissipation at low Mach numbers while recovering the standard scheme at high Mach numbers.
+
+Additionally, the **artificial Mach number** technique (`pi_fac`) modifies \f$\pi_\infty\f$ to artificially increase the local Mach number.
+
+---
+
+## 18. Flux Limiting
+
+**Volume fraction limiting:** enforces \f$0 \le \alpha_k \le 1\f$ and rescales as:
+
+\f[\alpha_k \leftarrow \frac{\alpha_k}{\sum_j \alpha_j}\f]
+
+**Advective flux limiting** based on local volume fraction gradient \f$\chi\f$ to prevent spurious oscillations near material interfaces.
+
diff --git a/docs/documentation/readme.md b/docs/documentation/readme.md
index 1d70332830..d72395cae0 100644
--- a/docs/documentation/readme.md
+++ b/docs/documentation/readme.md
@@ -11,6 +11,7 @@ Welcome to the Multi-component Flow Code (MFC) documentation.
## Reference
+- @ref equations "Equations" - Comprehensive equations reference
- @ref parameters "Case Parameters" - All ~3,400 parameters
- @ref cli-reference "CLI Reference" - Command line options
- @ref case_constraints "Case Creator Guide" - Feature compatibility
diff --git a/docs/documentation/references.md b/docs/documentation/references.md
index 184dd82b0b..7a9f1cfe6c 100644
--- a/docs/documentation/references.md
+++ b/docs/documentation/references.md
@@ -2,74 +2,6 @@
# References
-- Allaire, G., Clerc, S., and Kokh, S. (2002). A five-equation model for the simulation of interfaces between compressible fluids. Journal of Computational Physics, 181(2):577–616.
+All references cited throughout the MFC documentation are managed via BibTeX and rendered automatically by Doxygen.
-- Ando, K. (2010). Effects of polydispersity in bubbly flows. PhD thesis, California Institute of Technology.
-
-- Balsara, D. S. and Shu, C.-W. (2000). Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. Journal of Computational Physics, 160(2):405–452.
-
-- Batten, P., Clarke, N., Lambert, C., and Causon, D. M. (1997). On the choice of wavespeeds for the hllc riemann solver. SIAM Journal on Scientific Computing, 18(6):1553–1570.
-
-- Bryngelson, S. H., Schmidmayer, K., Coralic, V., Meng, J. C., Maeda, K., and Colonius, T. (2019). Mfc: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver. arXiv preprint arXiv:1907.10512.
-
-- Cao A. and. Schäfer F. (2024). Information Geometric Regularization of the Barotropic Euler Equations. arXiv preprint arXiv:2308.14127.
-
-- Chen, S. S., Li, J. P., Li, Z., Yuan, W., & Gao, Z. H. (2022). Anti-dissipation pressure correction under low Mach numbers for Godunov-type schemes. Journal of Computational Physics, 456, 111027.
-
-- Childs, H., Brugger, E., Whitlock, B., Meredith, J., Ahern, S., Pugmire, D., Biagas, K., Miller, M., Harrison, C., Weber, G. H., Krishnan, H., Fogal, T., Sanderson, A., Garth, C., Bethel, E. W., Camp, D., R¨ubel, O., Durant, M., Favre, J. M., and Navr´atil, P. (2012). VisIt: An End-User Tool For Visualizing and Analyzing Very Large Data. In High Performance Visualization–Enabling Extreme-Scale Scientific Insight, pages 357–372.
-
-- Coralic, V. (2015). Simulation of shock-induced bubble collapse with application to vascular injury in shockwave lithotripsy. PhD thesis, California Institute of Technology.
-
-- Coralic, V. and Colonius, T. (2014). Finite-volume weno scheme for viscous compressible multicomponent flows. Journal of computational physics, 274:95–121.
-
-- Gottlieb, S. and Shu, C.-W. (1998). Total variation diminishing runge-kutta schemes. Mathematics of computation of the American Mathematical Society, 67(221):73–85.
-
-- Guo, H., Jiang, P., Ye, L., & Zhu, Y. (2023). An efficient and low-divergence method for generating inhomogeneous and anisotropic turbulence with arbitrary spectra. Journal of Fluid Mechanics, 970, A2.
-
-- Henrick, A. K., Aslam, T. D., and Powers, J. M. (2005). Mapped weighted essentially nonoscillatory schemes: achieving optimal order near critical points. Journal of Computational Physics, 207(2):542–567.
-
-- Borges, R., Carmona, M., Costa, B., and Don, W. S. (2008). An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws. Journal of computational physics, 227(6):3191–3211.
-
-- Fu, L., Hu, X. Y., and Adams, N. A. (2016). A family of high-order targeted ENO schemes for compressible-fluid simulations. Journal of Computational Physics, 305:333–359.
-
-- Johnsen, E. (2008). Numerical simulations of non-spherical bubble collapse: With applications to shockwave lithotripsy. PhD thesis, California Institute of Technology.
-
-- Maeda, K. and Colonius, T. (2017). A source term approach for generation of one-way acoustic waves in the euler and navier–stokes equations. Wave Motion, 75:36–49.
-
-- Maeda, K. and Colonius, T. (2018). Eulerian–lagrangian method for simulation of cloud cavitation. Journal of computational physics, 371:994–1017.
-
-- Meng, J. C. C. (2016). Numerical simulations of droplet aerobreakup. PhD thesis, California Institute of Technology.
-
-- Michalke, A. (1964). On the inviscid instability of the hyperbolictangent velocity profile. Journal of Fluid Mechanics, 19(4), 543-556.
-
-- Pirozzoli, S., and Colonius, T. (2013). Generalized characteristic relaxation boundary conditions for unsteady compressible flow simulations. Journal of Computational Physics, 248:109-126.
-
-- Preston, A., Colonius, T., and Brennen, C. (2007). A reduced-order model of diffusive effects on the dynamics of bubbles. Physics of Fluids, 19(12):123302.
-
-- Saurel, R., Petitpas, F., and Berry, R. A. (2009). Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures. journal of Computational Physics, 228(5):1678–1712
-
-- Schmidmayer, K., Bryngelson, S. H., and Colonius, T. (2019). An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics. arXiv preprint arXiv:1903.08242.
-
-- Suresh, A. and Huynh, H. (1997). Accurate monotonicity-preserving schemes with runge–kutta time stepping. Journal of Computational Physics, 136(1):83–99.
-
-- Tam, C. K., Ju, H., Jones, M. G., Watson, W. R., and Parrott, T. L. (2005). A computational and experimental study of slit resonators. Journal of Sound and Vibration, 284(3-5), 947-984.
-
-- Thompson, K. W. (1987). Time dependent boundary conditions for hyperbolic systems. Journal of computational physics, 68(1):1–24.
-
-- Thompson, K. W. (1990). Time-dependent boundary conditions for hyperbolic systems, ii. Journal of computational physics, 89(2):439–461.
-
-- Thornber, B., Mosedale, A., Drikakis, D., Youngs, D., & Williams, R. J. (2008). An improved reconstruction method for compressible flows with low Mach number features. Journal of computational Physics, 227(10), 4873-4894.
-
-- Titarev, V. A. and Toro, E. F. (2004). Finite-volume weno schemes for three-dimensional conservation laws. Journal of Computational Physics, 201(1):238–260.
-
-- Tiwari, A., Freund, J. B., and Pantano, C. (2013). A diffuse interface model with immiscibility preservation. Journal of computational physics, 252:290–309.
-
-- Toro, E. F. (2013). Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Science & Business Media.
-
-- Miyoishi, T., and Kusano, K. (2005). A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics. Journal of Computational Physics, 208(1), 315-344.
-
-- Dedner, A., Kemm, F., Kröner, D., Munz, C. D., Schnitzer, T., & Wesenberg, M. (2002). Hyperbolic divergence cleaning for the MHD equations. Journal of Computational Physics, 175(2), 645-673.
-
-- Cao, S., Zhang, Y., Liao, D., Zhong, P., and Wang, K. G. (2019). Shock-induced damage and dynamic fracture in cylindrical bodies submerged in liquid. International Journal of Solids and Structures, 169:55–71. Elsevier.
-
-- Xu, W., Gao, Y., Deng, Y., Liu, J., and Liu, C. (2019). An explicit expression for the calculation of the Rortex vector. Physics of Fluids, 31(9)..
\ No newline at end of file
+See the \ref citelist "Bibliography" for the full list of cited works.
diff --git a/docs/documentation/visualization.md b/docs/documentation/visualization.md
index ac40d10b1c..10c1fca271 100644
--- a/docs/documentation/visualization.md
+++ b/docs/documentation/visualization.md
@@ -40,7 +40,7 @@ For many cases, this step will require resizing the render view window.
VisIt is an alternative open-source interactive parallel visualization and graphical analysis tool for viewing scientific data.
Versions of VisIt after 2.6.0 have been confirmed to work with the MFC databases for some parallel environments.
Nevertheless, installation and configuration of VisIt can be environment-dependent and are left to the user.
-Further remarks on parallel flow visualization, analysis, and processing of the MFC database using VisIt can also be found in [Coralic (2015)](@ref references) and [Meng (2016)](@ref references).
+Further remarks on parallel flow visualization, analysis, and processing of the MFC database using VisIt can also be found in \cite Coralic15 and \cite Meng16.
The user can launch VisIt and open the index files under `/silo_hdf5/root`.
Once the database is loaded, flow field variables contained in the database can be added to the plot.
diff --git a/docs/postprocess_citations.py b/docs/postprocess_citations.py
new file mode 100644
index 0000000000..a5c19df5c3
--- /dev/null
+++ b/docs/postprocess_citations.py
@@ -0,0 +1,212 @@
+"""Post-process Doxygen HTML to add author names before citation numbers.
+
+Transforms bare "[N]" citation links into "Author et al. [N]" style.
+Parses docs/references.bib for author information and walks all HTML
+files in the Doxygen output directory to insert author text.
+
+Usage: python3 postprocess_citations.py
+"""
+
+import re
+import sys
+from pathlib import Path
+
+
+def parse_bib_authors(bib_path: Path) -> dict[str, tuple[str, int, str | None]]:
+ """Parse references.bib to extract first-author surname and author count.
+
+ Returns:
+ Dict mapping lowercase bib key ->
+ (first_author_surname, author_count, second_author_surname | None)
+ """
+ text = bib_path.read_text(encoding="utf-8")
+ entries: dict[str, tuple[str, int, str | None]] = {}
+
+ # Find all bib entries: @type{Key,
+ for m in re.finditer(r"@\w+\{(\w+),", text):
+ key = m.group(1)
+ # Find the author field for this entry
+ start = m.end()
+ # Look for author = {..} or author = "..."
+ author_match = re.search(
+ r"author\s*=\s*\{([^{}]*(?:\{[^}]*\}[^{}]*)*)\}",
+ text[start:start + 2000],
+ )
+ if not author_match:
+ continue
+
+ author_str = author_match.group(1)
+ # Split by " and " to get individual authors
+ authors = re.split(r"\s+and\s+", author_str)
+ count = len(authors)
+
+ # Extract surname of first author
+ first = authors[0].strip()
+ surname = _extract_surname(first)
+
+ # Also get second author surname for 2-author papers
+ second_surname = None
+ if count == 2:
+ second_surname = _extract_surname(authors[1].strip())
+
+ entries[key.lower()] = (surname, count, second_surname)
+
+ return entries
+
+
+def _extract_surname(author: str) -> str:
+ """Extract surname from a BibTeX author name.
+
+ Handles:
+ - "G. Allaire" -> "Allaire"
+ - "O. {Le M\\'etayer}" -> "Le Métayer"
+ - "{Lord Rayleigh}" -> "Lord Rayleigh"
+ - "B. {van Leer}" -> "van Leer"
+ - "M. {Rodriguez Jr.}" -> "Rodriguez Jr."
+ """
+ # If the whole name is braced, return it (e.g., {Lord Rayleigh})
+ if author.startswith("{") and author.endswith("}"):
+ return _clean_latex(author[1:-1])
+
+ # Check for braced surname at the end: "O. {Le Métayer}"
+ brace_match = re.search(r"\{([^}]+)\}\s*$", author)
+ if brace_match:
+ return _clean_latex(brace_match.group(1))
+
+ # Standard case: last word is surname
+ parts = author.split()
+ if parts:
+ return _clean_latex(parts[-1])
+ return author
+
+
+def _clean_latex(s: str) -> str:
+ """Remove LaTeX formatting from a string."""
+ s = s.replace("\\'", "") # accent commands
+ s = s.replace('\\"', "")
+ s = s.replace("\\", "")
+ s = s.replace("{", "").replace("}", "")
+ return s.strip()
+
+
+# Pattern matching Doxygen citation links
+CITE_RE = re.compile(
+ r'\[(\d+)\]'
+)
+
+
+def process_html(html_path: Path, authors: dict) -> bool:
+ """Process a single HTML file, inserting author names before citations.
+
+ Returns True if the file was modified.
+ """
+ text = html_path.read_text(encoding="utf-8")
+ if "citelist.html#CITEREF_" not in text:
+ return False
+
+ def replace_cite(match: re.Match) -> str:
+ key = match.group(1)
+ full_link = match.group(0)
+
+ info = authors.get(key)
+ if not info:
+ return full_link
+
+ surname, count, second_surname = info
+
+ # Build author prefix
+ if count == 1:
+ prefix = surname
+ elif count == 2 and second_surname:
+ prefix = f"{surname} and {second_surname}"
+ else:
+ prefix = f"{surname} et al."
+
+ # Dedup: only skip if this exact prefix was already inserted
+ # (prevents double-insertion on re-runs of the script)
+ preceding = text[max(0, match.start() - len(prefix) - 5):match.start()]
+ if preceding.rstrip().endswith(prefix):
+ return full_link
+
+ return f"{prefix} {full_link}"
+
+ new_text = CITE_RE.sub(replace_cite, text)
+ if new_text != text:
+ html_path.write_text(new_text, encoding="utf-8")
+ return True
+ return False
+
+
+def check_bare_citations(html_dir: Path, authors: dict) -> list[str]:
+ """Find citation links that lack an author-name prefix.
+
+ Returns a list of warning strings for each bare citation found.
+ """
+ bare = []
+ for html_file in sorted(html_dir.rglob("*.html")):
+ if html_file.name == "citelist.html":
+ continue
+ text = html_file.read_text(encoding="utf-8")
+ for m in CITE_RE.finditer(text):
+ key = m.group(1)
+ if key not in authors:
+ bare.append(f" {html_file.name}: [{m.group(2)}] (CITEREF_{key}) — no bib entry")
+ continue
+ # Check if an author name precedes the link
+ before = text[max(0, m.start() - 60):m.start()].rstrip()
+ surname, _count, second_surname = authors[key]
+ has_prefix = (
+ before.endswith(surname) or
+ "et al." in before[-20:] or
+ (second_surname and before.endswith(second_surname))
+ )
+ if not has_prefix:
+ bare.append(f" {html_file.name}: [{m.group(2)}] (CITEREF_{key}) — missing author prefix")
+ return bare
+
+
+def main():
+ check_only = "--check" in sys.argv
+ args = [a for a in sys.argv[1:] if a != "--check"]
+
+ if len(args) != 1:
+ print(f"Usage: {sys.argv[0]} [--check] ", file=sys.stderr)
+ sys.exit(1)
+
+ html_dir = Path(args[0])
+ if not html_dir.is_dir():
+ print(f"Error: {html_dir} is not a directory", file=sys.stderr)
+ sys.exit(1)
+
+ # Find references.bib relative to this script
+ script_dir = Path(__file__).resolve().parent
+ bib_path = script_dir / "references.bib"
+ if not bib_path.exists():
+ print(f"Error: {bib_path} not found", file=sys.stderr)
+ sys.exit(1)
+
+ authors = parse_bib_authors(bib_path)
+ print(f" Parsed {len(authors)} bib entries")
+
+ if not check_only:
+ modified = 0
+ html_files = list(html_dir.rglob("*.html"))
+ for html_file in html_files:
+ if html_file.name == "citelist.html":
+ continue
+ if process_html(html_file, authors):
+ modified += 1
+ print(f" Updated citations in {modified}/{len(html_files)} HTML files")
+
+ # Verify all citations have author prefixes
+ bare = check_bare_citations(html_dir, authors)
+ if bare:
+ print("Bare citations found (missing author prefix):")
+ for b in bare:
+ print(b)
+ sys.exit(1)
+ print(" All citations have author prefixes")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/docs/references.bib b/docs/references.bib
new file mode 100644
index 0000000000..3aabb266bf
--- /dev/null
+++ b/docs/references.bib
@@ -0,0 +1,661 @@
+% ==========================================================================
+% MFC Documentation References
+% Used by Doxygen's \cite command (CITE_BIB_FILES in Doxyfile.in)
+% ==========================================================================
+
+% --- Governing equations / multiphase models ---
+
+@article{Allaire02,
+ author = {G. Allaire and S. Clerc and S. Kokh},
+ title = {A five-equation model for the simulation of interfaces between compressible fluids},
+ journal = {Journal of Computational Physics},
+ volume = {181},
+ number = {2},
+ pages = {577--616},
+ year = {2002},
+ doi = {10.1006/jcph.2002.7143}
+}
+
+@article{Kapila01,
+ author = {A. K. Kapila and R. Menikoff and J. B. Bdzil and S. F. Son and D. S. Stewart},
+ title = {Two-phase modeling of deflagration-to-detonation transition in granular materials: {R}educed equations},
+ journal = {Physics of Fluids},
+ volume = {13},
+ number = {10},
+ pages = {3002--3024},
+ year = {2001},
+ doi = {10.1063/1.1398042}
+}
+
+@article{Saurel09,
+ author = {R. Saurel and F. Petitpas and R. A. Berry},
+ title = {Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures},
+ journal = {Journal of Computational Physics},
+ volume = {228},
+ number = {5},
+ pages = {1678--1712},
+ year = {2009},
+ doi = {10.1016/j.jcp.2008.11.002}
+}
+
+@article{Saurel08,
+ author = {R. Saurel and F. Petitpas and R. Abgrall},
+ title = {Modelling phase transition in metastable liquids: application to cavitating and flashing flows},
+ journal = {Journal of Fluid Mechanics},
+ volume = {607},
+ pages = {313--350},
+ year = {2008},
+ doi = {10.1017/S0022112008002061}
+}
+
+@article{Zein10,
+ author = {A. Zein and M. Hantke and G. Warnecke},
+ title = {Modeling phase transition for compressible two-phase flows applied to metastable liquids},
+ journal = {Journal of Computational Physics},
+ volume = {229},
+ number = {8},
+ pages = {2964--2998},
+ year = {2010},
+ doi = {10.1016/j.jcp.2009.12.026}
+}
+
+% --- Equations of state ---
+
+@article{Menikoff89,
+ author = {R. Menikoff and B. J. Plohr},
+ title = {The {R}iemann problem for fluid flow of real materials},
+ journal = {Reviews of Modern Physics},
+ volume = {61},
+ number = {1},
+ pages = {75--130},
+ year = {1989},
+ doi = {10.1103/RevModPhys.61.75}
+}
+
+@article{LeMetayer04,
+ author = {O. {Le M\'etayer} and J. Massoni and R. Saurel},
+ title = {Elaborating equations of state of a liquid and its vapor for two-phase flow models},
+ journal = {International Journal of Thermal Sciences},
+ volume = {43},
+ number = {3},
+ pages = {265--276},
+ year = {2004},
+ doi = {10.1016/j.ijthermalsci.2003.09.002}
+}
+
+@techreport{McBride93,
+ author = {B. J. McBride and S. Gordon and M. A. Reno},
+ title = {Coefficients for calculating thermodynamic and transport properties of individual species},
+ institution = {NASA},
+ number = {TM-4513},
+ year = {1993}
+}
+
+% --- Bubble dynamics ---
+
+@article{Rayleigh17,
+ author = {{Lord Rayleigh}},
+ title = {On the pressure developed in a liquid during the collapse of a spherical cavity},
+ journal = {Philosophical Magazine},
+ volume = {34},
+ number = {200},
+ pages = {94--98},
+ year = {1917},
+ doi = {10.1080/14786440808635681}
+}
+
+@article{Plesset49,
+ author = {M. S. Plesset},
+ title = {The dynamics of cavitation bubbles},
+ journal = {Journal of Applied Mechanics},
+ volume = {16},
+ pages = {277--282},
+ year = {1949},
+ doi = {10.1115/1.4009975}
+}
+
+@article{Keller80,
+ author = {J. B. Keller and M. J. Miksis},
+ title = {Bubble oscillations of large amplitude},
+ journal = {Journal of the Acoustical Society of America},
+ volume = {68},
+ number = {2},
+ pages = {628--633},
+ year = {1980},
+ doi = {10.1121/1.384720}
+}
+
+@techreport{Gilmore52,
+ author = {F. R. Gilmore},
+ title = {The growth or collapse of a spherical bubble in a viscous compressible liquid},
+ institution = {California Institute of Technology},
+ number = {Report 26-4},
+ year = {1952}
+}
+
+@article{Commander89,
+ author = {K. W. Commander and A. Prosperetti},
+ title = {Linear pressure waves in bubbly liquids: {C}omparison between theory and experiments},
+ journal = {Journal of the Acoustical Society of America},
+ volume = {85},
+ number = {2},
+ pages = {732--746},
+ year = {1989},
+ doi = {10.1121/1.397599}
+}
+
+@article{Preston07,
+ author = {A. T. Preston and T. Colonius and C. E. Brennen},
+ title = {A reduced-order model of diffusive effects on the dynamics of bubbles},
+ journal = {Physics of Fluids},
+ volume = {19},
+ number = {12},
+ pages = {123302},
+ year = {2007},
+ doi = {10.1063/1.2825018}
+}
+
+@phdthesis{Ando10,
+ author = {K. Ando},
+ title = {Effects of polydispersity in bubbly flows},
+ school = {California Institute of Technology},
+ year = {2010}
+}
+
+@article{Ando11,
+ author = {K. Ando and T. Colonius and C. E. Brennen},
+ title = {Numerical simulation of shock propagation in a polydisperse bubbly liquid},
+ journal = {International Journal of Multiphase Flow},
+ volume = {37},
+ number = {6},
+ pages = {596--608},
+ year = {2011},
+ doi = {10.1016/j.ijmultiphaseflow.2011.03.007}
+}
+
+@article{Bryngelson20,
+ author = {S. H. Bryngelson and T. Colonius and R. O. Fox},
+ title = {{QBMMlib}: {A} library of quadrature-based moment methods},
+ journal = {SoftwareX},
+ volume = {12},
+ pages = {100615},
+ year = {2020},
+ doi = {10.1016/j.softx.2020.100615}
+}
+
+% --- Euler-Lagrange bubbles ---
+
+@article{Maeda18,
+ author = {K. Maeda and T. Colonius},
+ title = {{Eulerian--Lagrangian} method for simulation of cloud cavitation},
+ journal = {Journal of Computational Physics},
+ volume = {371},
+ pages = {994--1017},
+ year = {2018},
+ doi = {10.1016/j.jcp.2018.05.029}
+}
+
+% --- Elasticity ---
+
+@article{Rodriguez19,
+ author = {M. Rodriguez and E. Johnsen},
+ title = {A high-order accurate five-equations compressible multiphase approach for viscoelastic fluids and solids with relaxation and elasticity},
+ journal = {Journal of Computational Physics},
+ volume = {379},
+ pages = {70--90},
+ year = {2019},
+ doi = {10.1016/j.jcp.2018.10.035}
+}
+
+@article{Kamrin12,
+ author = {K. Kamrin and C. H. Rycroft and J.-C. Nave},
+ title = {Reference map technique for finite-strain elasticity and fluid-solid interaction},
+ journal = {Journal of the Mechanics and Physics of Solids},
+ volume = {60},
+ number = {11},
+ pages = {1952--1969},
+ year = {2012},
+ doi = {10.1016/j.jmps.2012.06.003}
+}
+
+@article{Cao19,
+ author = {S. Cao and Y. Zhang and D. Liao and P. Zhong and K. G. Wang},
+ title = {Shock-induced damage and dynamic fracture in cylindrical bodies submerged in liquid},
+ journal = {International Journal of Solids and Structures},
+ volume = {169},
+ pages = {55--71},
+ year = {2019},
+ doi = {10.1016/j.ijsolstr.2019.04.002}
+}
+
+% --- Surface tension ---
+
+@article{Schmidmayer17,
+ author = {K. Schmidmayer and F. Petitpas and E. Daniel and N. Favrie and S. L. Gavrilyuk},
+ title = {A model and numerical method for compressible flows with capillary effects},
+ journal = {Journal of Computational Physics},
+ volume = {334},
+ pages = {468--496},
+ year = {2017},
+ doi = {10.1016/j.jcp.2017.01.001}
+}
+
+% --- Chemistry ---
+
+@article{Cisneros26,
+ author = {E. Cisneros-Garibay and H. {Le Berre} and D. Adam and S. H. Bryngelson and J. B. Freund},
+ title = {Pyrometheus: {S}ymbolic abstractions for {XPU} and automatically differentiated computation of combustion kinetics and thermodynamics},
+ journal = {Computer Physics Communications},
+ volume = {320},
+ pages = {109987},
+ year = {2026},
+ doi = {10.1016/j.cpc.2025.109987}
+}
+
+% --- MHD ---
+
+@article{Miyoshi05,
+ author = {T. Miyoshi and K. Kusano},
+ title = {A multi-state {HLL} approximate {R}iemann solver for ideal magnetohydrodynamics},
+ journal = {Journal of Computational Physics},
+ volume = {208},
+ number = {1},
+ pages = {315--344},
+ year = {2005},
+ doi = {10.1016/j.jcp.2005.02.017}
+}
+
+@article{Dedner02,
+ author = {A. Dedner and F. Kemm and D. Kr{\"o}ner and C.-D. Munz and T. Schnitzer and M. Wesenberg},
+ title = {Hyperbolic divergence cleaning for the {MHD} equations},
+ journal = {Journal of Computational Physics},
+ volume = {175},
+ number = {2},
+ pages = {645--673},
+ year = {2002},
+ doi = {10.1006/jcph.2001.6961}
+}
+
+% --- IGR ---
+
+@inproceedings{Wilfong25a,
+ author = {B. Wilfong and A. Radhakrishnan and H. {Le Berre} and D. J. Vickers and T. Prathi and N. Tselepidis and B. Dorschner and R. Budiardja and B. Cornille and S. Abbott and F. Sch{\"a}fer and S. H. Bryngelson},
+ title = {Simulating many-engine spacecraft: {E}xceeding 1 quadrillion degrees of freedom via information geometric regularization},
+ booktitle = {SC'25: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis},
+ pages = {14--24},
+ year = {2025},
+ doi = {10.1145/3712285.3771783},
+ note = {{*}Equal contribution}
+}
+
+@article{Cao24,
+ author = {A. Cao and F. Sch{\"a}fer},
+ title = {Information geometric regularization of the barotropic {E}uler equations},
+ journal = {arXiv preprint arXiv:2308.14127},
+ year = {2024},
+ doi = {10.48550/arXiv.2308.14127}
+}
+
+% --- Acoustic sources ---
+
+@article{Maeda17,
+ author = {K. Maeda and T. Colonius},
+ title = {A source term approach for generation of one-way acoustic waves in the {E}uler and {N}avier--{S}tokes equations},
+ journal = {Wave Motion},
+ volume = {75},
+ pages = {36--49},
+ year = {2017},
+ doi = {10.1016/j.wavemoti.2017.08.004}
+}
+
+% --- Riemann solvers ---
+
+@article{Harten83,
+ author = {A. Harten and P. D. Lax and B. {van Leer}},
+ title = {On upstream differencing and {G}odunov-type schemes for hyperbolic conservation laws},
+ journal = {SIAM Review},
+ volume = {25},
+ number = {1},
+ pages = {35--61},
+ year = {1983},
+ doi = {10.1137/1025002}
+}
+
+@article{Toro94,
+ author = {E. F. Toro and M. Spruce and W. Speares},
+ title = {Restoration of the contact surface in the {HLL-Riemann} solver},
+ journal = {Shock Waves},
+ volume = {4},
+ pages = {25--34},
+ year = {1994},
+ doi = {10.1007/BF01414629}
+}
+
+@book{Toro09,
+ author = {E. F. Toro},
+ title = {Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction},
+ publisher = {Springer},
+ edition = {3rd},
+ year = {2009},
+ doi = {10.1007/b79761}
+}
+
+@article{Batten97,
+ author = {P. Batten and N. Clarke and C. Lambert and D. M. Causon},
+ title = {On the choice of wavespeeds for the {HLLC} {R}iemann solver},
+ journal = {SIAM Journal on Scientific Computing},
+ volume = {18},
+ number = {6},
+ pages = {1553--1570},
+ year = {1997},
+ doi = {10.1137/S1064827593260140}
+}
+
+% --- WENO / spatial reconstruction ---
+
+@article{Jiang96,
+ author = {G.-S. Jiang and C.-W. Shu},
+ title = {Efficient implementation of weighted {ENO} schemes},
+ journal = {Journal of Computational Physics},
+ volume = {126},
+ number = {2},
+ pages = {202--228},
+ year = {1996},
+ doi = {10.1006/jcph.1996.0130}
+}
+
+@article{Henrick05,
+ author = {A. K. Henrick and T. D. Aslam and J. M. Powers},
+ title = {Mapped weighted essentially non-oscillatory schemes: {A}chieving optimal order near critical points},
+ journal = {Journal of Computational Physics},
+ volume = {207},
+ number = {2},
+ pages = {542--567},
+ year = {2005},
+ doi = {10.1016/j.jcp.2005.01.023}
+}
+
+@article{Borges08,
+ author = {R. Borges and M. Carmona and B. Costa and W. S. Don},
+ title = {An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws},
+ journal = {Journal of Computational Physics},
+ volume = {227},
+ number = {6},
+ pages = {3191--3211},
+ year = {2008},
+ doi = {10.1016/j.jcp.2007.11.038}
+}
+
+@article{Fu16,
+ author = {L. Fu and X. Y. Hu and N. A. Adams},
+ title = {A family of high-order targeted {ENO} schemes for compressible-fluid simulations},
+ journal = {Journal of Computational Physics},
+ volume = {305},
+ pages = {333--359},
+ year = {2016},
+ doi = {10.1016/j.jcp.2015.10.037}
+}
+
+@article{Balsara00,
+ author = {D. S. Balsara and C.-W. Shu},
+ title = {Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy},
+ journal = {Journal of Computational Physics},
+ volume = {160},
+ number = {2},
+ pages = {405--452},
+ year = {2000},
+ doi = {10.1006/jcph.2000.6443}
+}
+
+@article{Suresh97,
+ author = {A. Suresh and H. Huynh},
+ title = {Accurate monotonicity-preserving schemes with {R}unge--{K}utta time stepping},
+ journal = {Journal of Computational Physics},
+ volume = {136},
+ number = {1},
+ pages = {83--99},
+ year = {1997},
+ doi = {10.1006/jcph.1997.5745}
+}
+
+@article{Titarev04,
+ author = {V. A. Titarev and E. F. Toro},
+ title = {Finite-volume {WENO} schemes for three-dimensional conservation laws},
+ journal = {Journal of Computational Physics},
+ volume = {201},
+ number = {1},
+ pages = {238--260},
+ year = {2004},
+ doi = {10.1016/j.jcp.2004.05.015}
+}
+
+% --- Time integration ---
+
+@article{Gottlieb98,
+ author = {S. Gottlieb and C.-W. Shu},
+ title = {Total variation diminishing {R}unge--{K}utta schemes},
+ journal = {Mathematics of Computation},
+ volume = {67},
+ number = {221},
+ pages = {73--85},
+ year = {1998},
+ doi = {10.1090/S0025-5718-98-00913-2}
+}
+
+@article{Strang68,
+ author = {G. Strang},
+ title = {On the construction and comparison of difference schemes},
+ journal = {SIAM Journal on Numerical Analysis},
+ volume = {5},
+ number = {3},
+ pages = {506--517},
+ year = {1968},
+ doi = {10.1137/0705041}
+}
+
+% --- Boundary conditions ---
+
+@article{Thompson87,
+ author = {K. W. Thompson},
+ title = {Time dependent boundary conditions for hyperbolic systems},
+ journal = {Journal of Computational Physics},
+ volume = {68},
+ number = {1},
+ pages = {1--24},
+ year = {1987},
+ doi = {10.1016/0021-9991(87)90041-6}
+}
+
+@article{Thompson90,
+ author = {K. W. Thompson},
+ title = {Time-dependent boundary conditions for hyperbolic systems, {II}},
+ journal = {Journal of Computational Physics},
+ volume = {89},
+ number = {2},
+ pages = {439--461},
+ year = {1990},
+ doi = {10.1016/0021-9991(90)90152-Q}
+}
+
+@article{Pirozzoli13,
+ author = {S. Pirozzoli and T. Colonius},
+ title = {Generalized characteristic relaxation boundary conditions for unsteady compressible flow simulations},
+ journal = {Journal of Computational Physics},
+ volume = {248},
+ pages = {109--126},
+ year = {2013},
+ doi = {10.1016/j.jcp.2013.04.021}
+}
+
+% --- Immersed boundary ---
+
+@article{Tseng03,
+ author = {Y. H. Tseng and J. H. Ferziger},
+ title = {A ghost-cell immersed boundary method for flow in complex geometry},
+ journal = {Journal of Computational Physics},
+ volume = {192},
+ number = {2},
+ pages = {593--623},
+ year = {2003},
+ doi = {10.1016/j.jcp.2003.07.024}
+}
+
+@article{Mittal05,
+ author = {R. Mittal and G. Iaccarino},
+ title = {Immersed boundary methods},
+ journal = {Annual Review of Fluid Mechanics},
+ volume = {37},
+ pages = {239--261},
+ year = {2005},
+ doi = {10.1146/annurev.fluid.37.061903.175743}
+}
+
+% --- Low Mach corrections ---
+
+@article{Thornber08,
+ author = {B. Thornber and A. Mosedale and D. Drikakis and D. Youngs and R. J. R. Williams},
+ title = {An improved reconstruction method for compressible flows with low {M}ach number features},
+ journal = {Journal of Computational Physics},
+ volume = {227},
+ number = {10},
+ pages = {4873--4894},
+ year = {2008},
+ doi = {10.1016/j.jcp.2008.01.036}
+}
+
+@article{Chen22,
+ author = {S. S. Chen and J. P. Li and Z. Li and W. Yuan and Z. H. Gao},
+ title = {Anti-dissipation pressure correction under low {M}ach numbers for {G}odunov-type schemes},
+ journal = {Journal of Computational Physics},
+ volume = {456},
+ pages = {111027},
+ year = {2022},
+ doi = {10.1016/j.jcp.2022.111027}
+}
+
+% --- MFC papers ---
+
+@article{Wilfong26,
+ author = {B. Wilfong and H. {Le Berre} and A. Radhakrishnan and A. Gupta and D. J. Vickers and D. Vaca-Revelo and D. Adam and H. Yu and H. Lee and J. R. Chreim and M. {Carcana Barbosa} and Y. Zhang and E. Cisneros-Garibay and A. Gnanaskandan and M. {Rodriguez Jr.} and R. D. Budiardja and S. Abbott and T. Colonius and S. H. Bryngelson},
+ title = {{MFC 5.0}: {A}n exascale many-physics flow solver},
+ journal = {Computer Physics Communications},
+ volume = {322},
+ pages = {110055},
+ year = {2026},
+ doi = {10.1016/j.cpc.2026.110055}
+}
+
+@article{Bryngelson21,
+ author = {S. H. Bryngelson and K. Schmidmayer and V. Coralic and J. C. Meng and K. Maeda and T. Colonius},
+ title = {{MFC}: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver},
+ journal = {Computer Physics Communications},
+ volume = {266},
+ pages = {107396},
+ year = {2021},
+ doi = {10.1016/j.cpc.2020.107396}
+}
+
+@article{Schmidmayer20,
+ author = {K. Schmidmayer and S. H. Bryngelson and T. Colonius},
+ title = {An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics},
+ journal = {Journal of Computational Physics},
+ volume = {402},
+ pages = {109080},
+ year = {2020},
+ doi = {10.1016/j.jcp.2019.109080}
+}
+
+@article{Coralic14,
+ author = {V. Coralic and T. Colonius},
+ title = {Finite-volume {WENO} scheme for viscous compressible multicomponent flows},
+ journal = {Journal of Computational Physics},
+ volume = {274},
+ pages = {95--121},
+ year = {2014},
+ doi = {10.1016/j.jcp.2014.06.003}
+}
+
+% --- Theses ---
+
+@phdthesis{Coralic15,
+ author = {V. Coralic},
+ title = {Simulation of shock-induced bubble collapse with application to vascular injury in shockwave lithotripsy},
+ school = {California Institute of Technology},
+ year = {2015}
+}
+
+@phdthesis{Johnsen08,
+ author = {E. Johnsen},
+ title = {Numerical simulations of non-spherical bubble collapse: with applications to shockwave lithotripsy},
+ school = {California Institute of Technology},
+ year = {2008}
+}
+
+@phdthesis{Meng16,
+ author = {J. C. C. Meng},
+ title = {Numerical simulations of droplet aerobreakup},
+ school = {California Institute of Technology},
+ year = {2016}
+}
+
+% --- Other references ---
+
+@article{Tiwari13,
+ author = {A. Tiwari and J. B. Freund and C. Pantano},
+ title = {A diffuse interface model with immiscibility preservation},
+ journal = {Journal of Computational Physics},
+ volume = {252},
+ pages = {290--309},
+ year = {2013},
+ doi = {10.1016/j.jcp.2013.06.021}
+}
+
+@article{Michalke64,
+ author = {A. Michalke},
+ title = {On the inviscid instability of the hyperbolic-tangent velocity profile},
+ journal = {Journal of Fluid Mechanics},
+ volume = {19},
+ number = {4},
+ pages = {543--556},
+ year = {1964},
+ doi = {10.1017/S0022112064000908}
+}
+
+@article{Tam05,
+ author = {C. K. Tam and H. Ju and M. G. Jones and W. R. Watson and T. L. Parrott},
+ title = {A computational and experimental study of slit resonators},
+ journal = {Journal of Sound and Vibration},
+ volume = {284},
+ number = {3--5},
+ pages = {947--984},
+ year = {2005},
+ doi = {10.1016/j.jsv.2004.07.013}
+}
+
+@article{Guo23,
+ author = {H. Guo and P. Jiang and L. Ye and Y. Zhu},
+ title = {An efficient and low-divergence method for generating inhomogeneous and anisotropic turbulence with arbitrary spectra},
+ journal = {Journal of Fluid Mechanics},
+ volume = {970},
+ pages = {A2},
+ year = {2023},
+ doi = {10.1017/jfm.2023.548}
+}
+
+@inproceedings{Childs12,
+ author = {H. Childs and E. Brugger and B. Whitlock and J. Meredith and S. Ahern and D. Pugmire and K. Biagas and M. Miller and C. Harrison and G. H. Weber and H. Krishnan and T. Fogal and A. Sanderson and C. Garth and E. W. Bethel and D. Camp and O. R{\"u}bel and M. Durant and J. M. Favre and P. Navr{\'a}til},
+ title = {{VisIt}: An end-user tool for visualizing and analyzing very large data},
+ booktitle = {High Performance Visualization---Enabling Extreme-Scale Scientific Insight},
+ pages = {357--372},
+ year = {2012}
+}
+
+@article{Xu2019,
+ author = {W. Xu and Y. Gao and Y. Deng and J. Liu and C. Liu},
+ title = {An explicit expression for the calculation of the {R}ortex vector},
+ journal = {Physics of Fluids},
+ volume = {31},
+ number = {9},
+ year = {2019},
+ doi = {10.1063/1.5116374}
+}
diff --git a/toolchain/mfc/lint_docs.py b/toolchain/mfc/lint_docs.py
index 4b46950905..f6318b8f03 100644
--- a/toolchain/mfc/lint_docs.py
+++ b/toolchain/mfc/lint_docs.py
@@ -1,4 +1,4 @@
-"""Check that file paths referenced in documentation still exist."""
+"""Check that file paths, cite keys, and parameters in docs still exist."""
import re
import sys
@@ -10,6 +10,7 @@
"docs/documentation/gpuParallelization.md",
"docs/documentation/running.md",
"docs/documentation/case.md",
+ "docs/documentation/equations.md",
".github/copilot-instructions.md",
]
@@ -19,8 +20,23 @@
# Skip paths with placeholders, globs, or patterns (not real file paths)
SKIP_RE = re.compile(r"[<>*()\[\]{}]|/\.\.\.|%|\$")
+# Match \cite Key references in markdown (Doxygen-style)
+CITE_RE = re.compile(r"\\cite\s+(\w+)")
+
+# Match backtick-quoted parameter names like `param = value`, `param > 0`, or `param`
+PARAM_RE = re.compile(r"`([a-z_][a-z0-9_%]*)(?:\s*[=><][^`]*)?\s*`")
+
+# Parameter-like names to skip (not actual MFC parameters)
+PARAM_SKIP = re.compile(
+ r"^(src/|toolchain/|\.github|docs/|examples/|tests/)" # file paths
+ r"|^\.(?:true|false)\.$" # Fortran logicals
+ r"|^\d" # numeric values
+ r"|^[A-Z]" # constants/types (uppercase start)
+)
+
def check_docs(repo_root: Path) -> list[str]:
+ """Check that file paths referenced in documentation still exist."""
errors = []
for doc in DOCS:
doc_path = repo_root / doc
@@ -38,12 +54,98 @@ def check_docs(repo_root: Path) -> list[str]:
return errors
+def parse_bib_keys(bib_path: Path) -> set[str]:
+ """Parse references.bib and return the set of valid citation keys (lowercase)."""
+ text = bib_path.read_text(encoding="utf-8")
+ return {m.group(1).lower() for m in re.finditer(r"@\w+\{(\w+),", text)}
+
+
+def check_cite_keys(repo_root: Path) -> list[str]:
+ """Check that \\cite keys in doc files exist in references.bib."""
+ bib_path = repo_root / "docs" / "references.bib"
+ if not bib_path.exists():
+ return []
+
+ valid_keys = parse_bib_keys(bib_path)
+ errors = []
+
+ doc_dir = repo_root / "docs" / "documentation"
+ if not doc_dir.exists():
+ return []
+
+ for md_file in sorted(doc_dir.glob("*.md")):
+ text = md_file.read_text(encoding="utf-8")
+ rel = md_file.relative_to(repo_root)
+ for match in CITE_RE.finditer(text):
+ key = match.group(1)
+ if key.lower() not in valid_keys:
+ errors.append(f" {rel} uses \\cite {key} but no bib entry found")
+
+ return errors
+
+
+def check_param_refs(repo_root: Path) -> list[str]:
+ """Check that parameter names in equations.md exist in the MFC registry."""
+ eq_path = repo_root / "docs" / "documentation" / "equations.md"
+ if not eq_path.exists():
+ return []
+
+ # Import REGISTRY from the toolchain
+ toolchain_dir = str(repo_root / "toolchain")
+ if toolchain_dir not in sys.path:
+ sys.path.insert(0, toolchain_dir)
+ try:
+ from mfc.params import REGISTRY # pylint: disable=import-outside-toplevel
+ except ImportError:
+ print(" Warning: could not import REGISTRY, skipping parameter check")
+ return []
+
+ valid_params = set(REGISTRY.all_params.keys())
+ # Build set of sub-parameter suffixes (the part after %)
+ sub_params = {p.split("%")[-1] for p in valid_params if "%" in p}
+ text = eq_path.read_text(encoding="utf-8")
+ errors = []
+ seen = set()
+
+ for match in PARAM_RE.finditer(text):
+ param = match.group(1)
+ if param in seen:
+ continue
+ seen.add(param)
+
+ # Skip non-parameter identifiers
+ if PARAM_SKIP.search(param):
+ continue
+ # Skip single-character names (too ambiguous: m, n, p are grid dims)
+ if len(param) <= 1:
+ continue
+ # Skip names containing % (struct members like x_domain%beg are valid
+ # but the base name before % is what matters)
+ base = param.split("%")[0] if "%" in param else param
+ # Check for indexed parameters: strip trailing _N (e.g., patch_icpp(1)%alpha(1))
+ # In docs these appear as e.g. `patch_icpp(i)%vel(j)` — skip indexed refs
+ if "(" in param or ")" in param:
+ continue
+
+ if base not in valid_params and base not in sub_params:
+ # Check if it's a known parameter family prefix
+ if not any(p.startswith(base) for p in valid_params):
+ errors.append(f" equations.md references parameter '{param}' not in REGISTRY")
+
+ return errors
+
+
def main():
repo_root = Path(__file__).resolve().parents[2]
- errors = check_docs(repo_root)
- if errors:
+
+ all_errors = []
+ all_errors.extend(check_docs(repo_root))
+ all_errors.extend(check_cite_keys(repo_root))
+ all_errors.extend(check_param_refs(repo_root))
+
+ if all_errors:
print("Doc reference check failed:")
- for e in errors:
+ for e in all_errors:
print(e)
sys.exit(1)
diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py
index 4ef006a304..c08249c3d8 100644
--- a/toolchain/mfc/params/definitions.py
+++ b/toolchain/mfc/params/definitions.py
@@ -781,11 +781,6 @@ def get_value_label(param_name: str, value: int) -> str:
"requires": ["sigma"],
}
},
- "mhd": {
- "when_true": {
- "recommends": ["hyper_cleaning"],
- }
- },
"relativity": {
"when_true": {
"requires": ["mhd"],