From 1a50fbad45f96eb0b28d79f4f12d4b44ef241adb Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Thu, 12 Feb 2026 11:41:50 -0500 Subject: [PATCH 1/6] Add comprehensive equations documentation Catalog every equation solved by MFC organized by physical model, including governing PDEs, equations of state, viscous stress, bubble dynamics, fluid-structure interaction, phase change, chemistry, surface tension, MHD, IGR, and all numerical methods. Each section references the activating input parameter, source file paths, and paper citations. Co-Authored-By: Claude Opus 4.6 --- docs/equations.md | 895 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 895 insertions(+) create mode 100644 docs/equations.md diff --git a/docs/equations.md b/docs/equations.md new file mode 100644 index 0000000000..401ff80629 --- /dev/null +++ b/docs/equations.md @@ -0,0 +1,895 @@ +# 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, cross-references the relevant source files, and cites the originating papers. + +**References:** +- **[bryngelson-CPC-20]** Bryngelson et al., "MFC: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver," *Computer Physics Communications* **266**, 107396 (2021). +- **[radhakrishnan-CPC-24]** Radhakrishnan et al., "Method for portable, scalable, and performant GPU-accelerated simulation of multiphase compressible flow," *Computer Physics Communications* **302**, 109238 (2024). +- **[wilfong-mfc-26]** Wilfong et al., "MFC 5.0: An exascale many-physics flow solver," *Computer Physics Communications* **322**, 110055 (2026). +- **[rodriguez-CAV-21]** Rodriguez et al., "Acoustically-induced bubble growth and phase change dynamics near compliant surfaces," *11th International Symposium on Cavitation* (2021). +- **[cisneros-cpc-25]** Cisneros-Garibay et al., "Pyrometheus: Symbolic abstractions for XPU and automatically differentiated computation of combustion kinetics and thermodynamics," *Computer Physics Communications* **320**, 109987 (2026). +- **[wilfong-igr-SC25]** Wilfong et al., "Information Geometric Regularization of the Barotropic Euler Equations," SC'25 (2025). +- **[wilfong-arxiv-25]** Wilfong et al., arxiv preprint (2025). + +--- + +## 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: + +$$\frac{\partial \mathbf{q}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{q}) + \mathbf{h}(\mathbf{q})\,\nabla \cdot \mathbf{u} = \mathbf{s}(\mathbf{q})$$ + +where: +- $\mathbf{q}$ is the conservative variable vector, +- $\mathbf{F}$ is the flux tensor, +- $\mathbf{h}(\mathbf{q})\,\nabla \cdot \mathbf{u}$ contains non-conservative terms (volume fraction advection), +- $\mathbf{s}(\mathbf{q})$ 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 (Allaire et al., 2002). The state vector is: + +$$\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$$ + +**Continuity** (one per component): + +$$\frac{\partial (\alpha_i \rho_i)}{\partial t} + \nabla \cdot (\alpha_i \rho_i\,\mathbf{u}) = 0$$ + +**Momentum:** + +$$\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot \bigl(\rho\,\mathbf{u} \otimes \mathbf{u} + p\,\mathbf{I} - \boldsymbol{\tau}^v\bigr) = 0$$ + +**Energy:** + +$$\frac{\partial (\rho E)}{\partial t} + \nabla \cdot \bigl[(\rho E + p)\,\mathbf{u} - \boldsymbol{\tau}^v \cdot \mathbf{u}\bigr] = 0$$ + +**Volume fraction advection:** + +$$\frac{\partial \alpha_i}{\partial t} + \mathbf{u} \cdot \nabla \alpha_i = K\,\nabla \cdot \mathbf{u}$$ + +where the $K$ term enforces interface conditions via the Wood sound speed: + +$$K = \frac{\rho_2 c_2^2 - \rho_1 c_1^2}{\dfrac{\rho_1 c_1^2}{\alpha_1} + \dfrac{\rho_2 c_2^2}{\alpha_2}}$$ + +Setting `alt_soundspeed = .true.` enables the $K$ correction (Kapila model with Wood sound speed). Setting `alt_soundspeed = .false.` uses the Allaire variant without the $K$ correction, which is conservative but does not strictly obey the second law of thermodynamics. + +**Mixture rules:** + +$$1 = \sum_i \alpha_i, \qquad \rho = \sum_i \alpha_i \rho_i, \qquad \rho e = \sum_i \alpha_i \rho_i e_i$$ + +**Ref:** [wilfong-mfc-26] Section 2.1.1, [bryngelson-CPC-20] Section 3.1, [radhakrishnan-CPC-24] Section 2.1 Eq. 2. + +### 2.2 Six-Equation Model (`model_eqns = 3`) + +Allows pressure disequilibrium between phases (Saurel et al., 2009). + +**Continuity and momentum:** Same as the five-equation model. + +**Separate phasic internal energy:** + +$$\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}$$ + +**Volume fraction:** + +$$\frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 = \mu\,(p_1 - p_2)$$ + +**Interfacial pressure:** + +$$p_I = \frac{z_2\,p_1 + z_1\,p_2}{z_1 + z_2}, \qquad z_i = \rho_i\,c_i$$ + +Infinite pressure relaxation is applied at each Runge-Kutta stage to drive toward pressure equilibrium. + +**Mixture speed of sound:** + +$$c^2 = \sum_k Y_k\,c_k^2$$ + +With phase change (`relax = .true.`), additional source terms appear in the phasic energy and volume fraction equations: +- **Pressure relaxation:** $\mu\,\delta p$ where $\delta p = p_1 - p_2$ +- **Thermal transfer:** $Q = \theta\,(T_2 - T_1)$ +- **Mass transfer:** $\dot{m} = \nu\,(g_2 - g_1)$ (Gibbs free energy difference) + +See [Section 8](#8-phase-change-relax--true) for details. + +**Ref:** [wilfong-mfc-26] Section 2.1.2, [bryngelson-CPC-20] Section 3.2, [rodriguez-CAV-21] Section 2.1, `src/simulation/m_pressure_relaxation.fpp`. + +### 2.3 Other Model Variants + +- `model_eqns = 1`: **Gamma/pi_inf model** — simplified single-fluid formulation using mixture $\gamma$ and $\pi_\infty$ directly without tracking individual volume fractions (Johnsen, 2008). +- `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 + +The primary closure for each phase: + +$$p_k = (\gamma_k - 1)\,\rho_k\,e_k - \gamma_k\,\pi_{\infty,k}$$ + +Equivalently: + +$$e_k = \frac{p_k + \gamma_k\,\pi_{\infty,k}}{(\gamma_k - 1)\,\rho_k}$$ + +**Total energy relation:** + +$$\rho E = \Gamma\,p + \Pi_\infty + \frac{1}{2}\rho\,|\mathbf{u}|^2 + q_v$$ + +where MFC internally tracks the transformed thermodynamic quantities: + +$$\Gamma_k = \frac{1}{\gamma_k - 1}, \qquad \Pi_{\infty,k} = \frac{\gamma_k\,\pi_{\infty,k}}{\gamma_k - 1}$$ + +and the mixture rules are arithmetic averages of these transformed quantities: + +$$\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}$$ + +The pressure is recovered from the total energy as: + +$$p = \frac{\rho E - \frac{1}{2}\rho\,|\mathbf{u}|^2 - \Pi_\infty - q_v}{\Gamma}$$ + +**Phasic speed of sound:** + +$$c_k = \sqrt{\frac{\gamma_k\,(p + \pi_{\infty,k})}{\rho_k}}$$ + +**Wood mixture sound speed:** + +$$\frac{1}{\rho\,c^2} = \sum_k \frac{\alpha_k}{\rho_k\,c_k^2}$$ + +Input parameters per fluid: `gamma` ($\gamma_k$), `pi_inf` ($\pi_{\infty,k}$), `cv` ($c_{v,k}$), `qv` ($q_{v,k}$), `qvp` ($q'_{v,k}$). + +**Ref:** [wilfong-mfc-26] Section 2.2, [radhakrishnan-CPC-24] Section 2.1 Eq. 1, `src/common/m_variables_conversion.fpp`. + +### 3.2 Ideal Gas EOS (Chemistry, `chemistry = .true.`) + +For reacting gas mixtures: + +$$p = \frac{\rho\,R_u\,T}{W}, \qquad W = \left(\sum_m \frac{Y_m}{W_m}\right)^{-1}$$ + +Temperature is obtained from the internal energy by Newton iteration: + +$$e_g - \sum_m e_m(T)\,Y_m = 0$$ + +**Species internal energy from enthalpy:** + +$$e_m(T) = \frac{\hat{h}_m(T) - R_u\,T}{W_m}$$ + +**NASA polynomial enthalpies:** + +$$\frac{\hat{h}_m}{R_u\,T} = \frac{C_0}{T} + \sum_{r=1}^{5} \frac{C_r\,T^{r-1}}{r}$$ + +**Ref:** [wilfong-mfc-26] Section 4.1.7 Eqs. 18-22. + +--- + +## 4. Viscous Stress Tensor (`viscous = .true.`) + +**Newtonian viscous stress (no bulk viscosity by default):** + +$$\boldsymbol{\tau}^v = 2\,\eta\left(\mathbf{D} - \frac{1}{3}\,\text{tr}(\mathbf{D})\,\mathbf{I}\right)$$ + +where the strain rate tensor is: + +$$\mathbf{D} = \frac{1}{2}\bigl(\nabla \mathbf{u} + (\nabla \mathbf{u})^T\bigr)$$ + +**With bulk viscosity:** + +$$\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}$$ + +**Cartesian components:** + +$$\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)$$ + +and similarly for all other components. Cylindrical coordinate formulations include additional $1/r$ terms. + +**Viscosity averaging:** + +$$\frac{1}{\text{Re}_\text{mix}} = \sum_j \frac{\alpha_j}{\text{Re}_j}$$ + +Input parameters: `Re_inv` (shear and volume Reynolds numbers per fluid). + +**Ref:** [wilfong-mfc-26] Eqs. 1-2, `src/simulation/m_viscous.fpp`. + +--- + +## 5. Cylindrical Coordinates (`cyl_coord = .true.`) + +Additional geometric source terms appear with $1/r$ factors in the continuity, momentum, and energy equations. Key modifications: + +- **Radial momentum:** extra $p/r$ and $\tau_{\theta\theta}/r$ terms +- **Viscous stress:** $\tau_{yy}$ includes $v/r$ corrections: + +$$\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)$$ + +- **Axis singularity:** axis placed at cell boundary with spectral filtering in the azimuthal direction + +**Ref:** [wilfong-mfc-26] Section 2.3, [bryngelson-CPC-20] Section 4.3. + +--- + +## 6. Sub-Grid Bubble Dynamics + +### 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 + +**Modified mixture pressure:** + +$$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)$$ + +**Modified stiffened gas for the liquid phase:** + +$$\Gamma_l\,p_l + \Pi_{\infty,l} = \frac{1}{1 - \alpha}\left(E - \frac{1}{2}\rho\,|\mathbf{u}|^2\right)$$ + +**Bubble wall pressure (polytropic):** + +$$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}$$ + +**Void fraction transport:** + +$$\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = \frac{3\,\alpha\,\bar{R}^2\,\dot{R}}{\bar{R}^3}$$ + +**Number density conservation:** + +$$\frac{\partial n_\text{bub}}{\partial t} + \nabla \cdot (n_\text{bub}\,\mathbf{u}) = 0$$ + +where $n = \frac{3}{4\pi}\,\frac{\alpha}{\bar{R}^3}$. + +**Polydispersity** (`polydisperse = .true.`): Log-normal PDF discretized into $N_\text{bin}$ equilibrium radii with standard deviation `poly_sigma`, integrated via Simpson's rule. + +#### 6.1.2 Rayleigh-Plesset (`bubble_model = 3`) + +$$R\,\ddot{R} + \frac{3}{2}\,\dot{R}^2 = \frac{p_{bw} - p_\infty}{\rho_l}$$ + +#### 6.1.3 Keller-Miksis (`bubble_model = 2`) + +$$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}$$ + +#### 6.1.4 Gilmore (`bubble_model = 1`) + +Enthalpy-based formulation with compressibility corrections via the Tait EOS: + +$$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)$$ + +where the enthalpy difference is: + +$$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]$$ + +and the local liquid sound speed: + +$$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}$$ + +#### 6.1.5 Non-Polytropic Thermal Model (`polytropic = .false.`) + +**Internal bubble pressure ODE:** + +$$\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)$$ + +**Vapor mass flux:** + +$$\dot{m}_v = \frac{D\,\rho_{bw}}{1 - \chi_{vw}}\left.\frac{\partial \chi_v}{\partial r}\right|_R$$ + +#### 6.1.6 QBMM Moment Transport (`qbmm = .true.`) + +**Population balance equation:** + +$$\frac{\partial f}{\partial t} + \frac{\partial (f\,\dot{R})}{\partial R} + \frac{\partial (f\,\ddot{R})}{\partial \dot{R}} = 0$$ + +**Moment transport:** + +$$\frac{\partial (n_\text{bub}\,\mu_i)}{\partial t} + \nabla \cdot (n_\text{bub}\,\mu_i\,\mathbf{u}) = n_\text{bub}\,\dot{\mu}_i$$ + +where moments $\mu_{i_1,i_2} = \int R^{i_1}\,\dot{R}^{i_2}\,f\,dR\,d\dot{R}$. + +**CHyQMOM inversion** recovers 4 quadrature nodes $(w_j, R_j, \dot{R}_j)$ from 6 moments via: + +$$\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$$ + +**Ref:** [wilfong-mfc-26] Section 4.1.4, `src/simulation/m_qbmm.fpp`. + +### 6.2 Euler-Lagrange Bubbles (`bubbles_lagrange = .true.`) + +**Source:** `src/simulation/m_bubbles_EL.fpp` + +Volume-averaged carrier flow equations with bubble source terms (Eq. 12 of [wilfong-mfc-26]): + +**Continuity:** + +$$\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]$$ + +**Momentum:** + +$$\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)$$ + +**Energy:** + +$$\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)$$ + +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:** + +$$\alpha(\mathbf{x}) = \sum_n V_n\,\delta_\sigma(\mathbf{x} - \mathbf{x}_n)$$ + +where $\delta_\sigma$ is a Gaussian kernel: + +$$\delta_\sigma(\mathbf{r}) = \frac{1}{(2\pi\sigma^2)^{3/2}}\exp\!\left(-\frac{|\mathbf{r}|^2}{2\sigma^2}\right)$$ + +with $\sigma = \varepsilon_b \max(\Delta x^{1/3}_\text{cell},\;R_\text{bubble})$. + +Each bubble is tracked individually with Keller-Miksis dynamics and 4th-order adaptive Runge-Kutta time integration. + +**Ref:** [wilfong-mfc-26] Section 4.1.5 Eqs. 13-14. + +--- + +## 7. Fluid-Structure Interaction + +### 7.1 Hypoelastic Model (`hypoelasticity = .true.`) + +**Source:** `src/simulation/m_hypoelastic.fpp` + +**Cauchy stress decomposition:** + +$$\sigma_{ij} = -p\,\delta_{ij} + \tau_{ij}^{(v)} + \tau_{ij}^{(e)}$$ + +**Elastic energy contribution to total energy:** + +$$E = e + \frac{|\mathbf{u}|^2}{2} + \frac{\boldsymbol{\tau}^e : \boldsymbol{\tau}^e}{4\,\rho\,G}$$ + +**Elastic stress evolution:** + +$$\frac{\partial (\rho\,\boldsymbol{\tau}^e)}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\tau}^e \otimes \mathbf{u}) = \mathbf{S}^e$$ + +**Source term:** + +$$\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)$$ + +where $\mathbf{l} = \nabla \mathbf{u}$ is the velocity gradient and $\mathbf{D}^d = \mathbf{D} - \frac{1}{3}\text{tr}(\mathbf{D})\,\mathbf{I}$ is the deviatoric strain rate. + +**Lie objective temporal derivative (Kelvin-Voigt):** + +$$\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$$ + +This adds 6 additional transport equations in 3D (symmetric stress tensor: $\tau_{xx}^e, \tau_{xy}^e, \tau_{yy}^e, \tau_{xz}^e, \tau_{yz}^e, \tau_{zz}^e$). + +**Ref:** [wilfong-mfc-26] Section 4.1.6, [radhakrishnan-CPC-24] Section 2.2 Eqs. 3-5. + +### 7.2 Hyperelastic Model (`hyperelasticity = .true.`) + +**Source:** `src/simulation/m_hyperelastic.fpp` + +**Reference map evolution:** + +$$\frac{\partial (\rho\,\boldsymbol{\xi})}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\xi} \otimes \mathbf{u}) = 0$$ + +**Deformation gradient from reference map:** + +$$\mathbf{F} = (\nabla \boldsymbol{\xi})^{-1}$$ + +**Left Cauchy-Green tensor:** + +$$\mathbf{b} = \mathbf{F}\,\mathbf{F}^T$$ + +**Neo-Hookean Cauchy stress:** + +$$\boldsymbol{\tau}^e = \frac{G}{J}\left(\mathbf{b} - \frac{\text{tr}(\mathbf{b})}{3}\,\mathbf{I}\right)$$ + +where $J = \det(\mathbf{F})$. + +**Hyperelastic energy:** + +$$e^e = \frac{G}{2}\bigl(I_{\mathbf{b}} - 3\bigr), \qquad I_{\mathbf{b}} = \text{tr}(\mathbf{b})$$ + +**Ref:** [wilfong-mfc-26] Section 4.1.6. + +--- + +## 8. Phase Change (`relax = .true.`) + +**Source:** `src/common/m_phase_change.fpp` + +### 8.1 pT-Relaxation (`relax_model = 5`) + +$N$-fluid pressure-temperature equilibrium. The equilibrium condition is: + +$$f(p) = \sum_i \alpha_i - 1 = 0$$ + +**Temperature from energy conservation:** + +$$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}$$ + +**Newton residual:** + +$$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}$$ + +Solved via Newton's method for the equilibrium pressure. + +### 8.2 pTg-Relaxation (`relax_model = 6`) + +Two coupled equations for $(\alpha_1 \rho_1,\;p)$: + +**Gibbs free energy equilibrium (Clausius-Clapeyron):** + +$$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$$ + +**Energy conservation constraint:** + +$$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$$ + +where $m_T$ is the total mass, $m_l = \alpha_l \rho_l$ is the liquid partial density, and $m_{qD}$, $m_{cpD}$, $m_{cvgp}$ are auxiliary thermodynamic sums over additional fluids (beyond the phase-changing pair). + +Solved via 2D Newton-Raphson. + +**Ref:** [wilfong-mfc-26] Section 4.1.3, [rodriguez-CAV-21] Section 2. + +--- + +## 9. Chemistry and Combustion (`chemistry = .true.`) + +**Source:** `src/common/m_chemistry.fpp` + +**Species transport:** + +$$\frac{\partial (\rho_g\,Y_m)}{\partial t} + \frac{\partial (\rho_g\,u_i\,Y_m)}{\partial x_i} = W_m\,\dot{\omega}_m$$ + +**Net production rate:** + +$$\dot{\omega}_m = \sum_n (\nu''_{mn} - \nu'_{mn})\,\mathcal{R}_n$$ + +**Reaction rate (law of mass action):** + +$$\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]$$ + +**Arrhenius rate:** + +$$k_n(T) = A_n\,T^{b_n}\exp\!\left(-\frac{T_{a,n}}{T}\right)$$ + +**Molecular diffusion** (`transport_model`): +- **Mixture-average:** Species-specific diffusion coefficients $D_m^\text{mix}$, mass flux: $\dot{m}_k = \rho\,D_k^\text{mix}\,(W_k / W_\text{mix})\,\partial X_k / \partial x$ +- **Unity Lewis number:** $D_m = \lambda / (\rho\,c_p)$ + +Enthalpy flux with diffusion: + +$$q_\text{diff} = \lambda\,\frac{\partial T}{\partial x} + \sum_k h_k\,\dot{m}_k$$ + +Reaction mechanisms are code-generated via Pyrometheus, which provides symbolic abstractions for thermochemistry that enable portable GPU computation and automatic differentiation of chemical source terms. + +**Ref:** [wilfong-mfc-26] Section 4.1.7 Eqs. 15-22, [cisneros-cpc-25] Section 5. + +--- + +## 10. Surface Tension (`surface_tension = .true.`) + +**Source:** `src/simulation/m_surface_tension.fpp`, `src/simulation/include/inline_capillary.fpp` + +**Color function advection:** + +$$\frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla c = 0$$ + +**Capillary stress tensor (CSF model):** + +$$\boldsymbol{\Omega} = -\sigma\left(\|\nabla c\|\,\mathbf{I} - \frac{\nabla c \otimes \nabla c}{\|\nabla c\|}\right)$$ + +In component form, with $\hat{w}_i = (\partial c / \partial x_i) / \|\nabla c\|$: + +$$\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\|$$ + +The capillary stress divergence is added to the momentum and energy equations. The total energy equation becomes: + +$$\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$$ + +**Capillary mixture energy:** + +$$\varepsilon_0 = \sigma\,\|\nabla c\|$$ + +**Ref:** [wilfong-mfc-26] Section 4.1.8. + +--- + +## 11. Magnetohydrodynamics + +### 11.1 Ideal MHD (`mhd = .true.`) + +**Continuity:** + +$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho\,\mathbf{u}) = 0$$ + +**Momentum:** + +$$\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$$ + +**Energy:** + +$$\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$$ + +**Induction:** + +$$\frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u}) = 0$$ + +**Total energy:** + +$$\mathcal{E} = \rho\,e + \frac{1}{2}\rho\,|\mathbf{u}|^2 + \frac{|\mathbf{B}|^2}{2}$$ + +**Fast magnetosonic speed:** + +$$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)}$$ + +**Alfven speed:** + +$$v_A = \sqrt{\frac{|\mathbf{B}|^2}{\rho}}$$ + +Uses the HLLD Riemann solver (`riemann_solver = 3`). Hyperbolic divergence cleaning (`hyper_cleaning = .true.`) via the GLM method (Dedner et al., 2002). + +**Ref:** [wilfong-mfc-26] Section 4.1.9. + +### 11.2 Relativistic MHD (`relativity = .true.`) + +**Conserved variables:** + +$$\mathbf{U} = (D,\;\mathbf{m},\;\tau,\;\mathbf{B})^T$$ + +where: + +$$D = \Gamma\,\rho, \qquad \mathbf{m} = \Gamma^2\rho h\,\mathbf{u} + |\mathbf{B}|^2\mathbf{u} - (\mathbf{u} \cdot \mathbf{B})\,\mathbf{B}$$ + +$$\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$$ + +Primitive recovery uses Newton-Raphson on the nonlinear conserved-to-primitive relation. + +**Ref:** [wilfong-mfc-26] Section 4.1.10. + +--- + +## 12. Information Geometric Regularization (`igr = .true.`) + +**Source:** `src/simulation/m_igr.fpp` + +**Modified momentum with entropic pressure $\Sigma$:** + +$$\frac{\partial (\rho\,\mathbf{u})}{\partial t} + \nabla \cdot \bigl[\rho\,\mathbf{u} \otimes \mathbf{u} + (p + \Sigma)\,\mathbf{I} - \boldsymbol{\tau}\bigr] = 0$$ + +**Elliptic PDE for $\Sigma$:** + +$$\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)$$ + +where $\alpha \sim \Delta x^2$ (regularization strength proportional to mesh spacing squared): + +$$\alpha_\text{IGR} = \alpha_\text{factor} \cdot \max(\Delta x,\;\Delta y,\;\Delta z)^2$$ + +**RHS strain-rate source (3D):** + +$$\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]$$ + +**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). + +**Ref:** [wilfong-arxiv-25] Eqs. 6-9. + +--- + +## 13. Body Forces (`bodyforces = .true.`) + +**Source:** `src/simulation/m_body_forces.fpp` + +**Time-dependent acceleration:** + +$$a_i(t) = g_i + k_i\sin(\omega_i\,t - \phi_i)$$ + +**Momentum source:** + +$$\frac{\partial (\rho\,u_i)}{\partial t}\bigg|_\text{bf} = \rho\,a_i(t)$$ + +**Energy source:** + +$$\frac{\partial E}{\partial t}\bigg|_\text{bf} = \rho\,\mathbf{u} \cdot \mathbf{a}(t)$$ + +--- + +## 14. Acoustic Sources (`acoustic_source = .true.`) + +**Source:** `src/simulation/m_acoustic_src.fpp` + +Source terms added to the RHS of the governing equations. + +**Discrete delta function (spatial support):** + +$$\delta_h(r) = \frac{1}{(2\pi\sigma^2)^{d/2}}\exp\!\left(-\frac{r^2}{2\sigma^2}\right)$$ + +**Forcing form (added to conservative variables):** + +$$\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$$ + +**Temporal profiles:** +- **Sine** (`pulse = 1`): $S(t) = M\sin(\omega(t - t_\text{delay}))$ +- **Gaussian** (`pulse = 2`): $S(t) = M\exp\!\bigl(-\tfrac{1}{2}((t - t_\text{delay})/\sigma_t)^2\bigr)$ +- **Square** (`pulse = 3`): $S(t) = M\,\text{sign}(\sin(\omega(t - t_\text{delay})))$ +- **Broadband** (`pulse = 4`): superposition of multiple frequencies across a bandwidth + +**Spatial supports:** planar, spherical transducer, cylindrical transducer, transducer array (arcuate, annular, circular). + +**Ref:** [bryngelson-CPC-20] Section 4.4 Eqs. 35-38. + +--- + +## 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_{i+1/2} = \sum_r \omega_r\,f_{i+1/2}^{(r)}$$ + +**WENO-JS** (default): + +$$\alpha_r = \frac{d_r}{(\beta_r + \varepsilon)^2}, \qquad \omega_r = \frac{\alpha_r}{\sum_s \alpha_s}$$ + +where $d_r$ are ideal weights, $\beta_r$ are smoothness indicators, and $\varepsilon$ is a small regularization parameter (`weno_eps`). + +**WENO-M** (`mapped_weno = .true.`): Henrick et al. (2005) mapped weights for improved accuracy at critical points: + +$$\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)}}$$ + +**WENO-Z** (`wenoz = .true.`): Borges et al. (2008) improved weights with global smoothness measure: + +$$\alpha_r = d_r\left(1 + \left(\frac{\tau}{\beta_r + \varepsilon}\right)^q\right), \qquad \tau = |\beta_0 - \beta_{k-1}|$$ + +The parameter $q$ controls the convergence rate at critical points (typically $q = 1$ for fifth-order reconstruction, as used in MFC). + +**TENO** (`teno = .true.`): Fu et al. (2016) targeted ENO with smoothness threshold $C_T$ (`teno_CT`): + +$$\gamma_r = 1 + \frac{\tau}{\beta_r}, \qquad \xi_r = \frac{\gamma_r}{\sum_s \gamma_s}$$ + +If $\xi_r < C_T$, set $\alpha_r = 0$ (stencil excluded). + +Primitive variable reconstruction is used to avoid spurious oscillations at interfaces. + +**Ref:** [wilfong-mfc-26] Sections 3.2.1, 4.2.2. + +#### MUSCL (`muscl_order = 2`) + +**Source:** `src/simulation/m_muscl.fpp` + +$$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$$ + +**Five slope limiters** (`muscl_lim`): +1. **Minmod:** $\phi = \text{sign}(a)\,\min(|a|,\;|b|)$ if $ab > 0$, else $0$ +2. **MC (Monotonized Central):** $\phi = \text{sign}(a)\,\min(2|a|,\;2|b|,\;\tfrac{1}{2}(|a|+|b|))$ if $ab > 0$ +3. **OSPRE (Van Albada):** $\phi = (a^2 b + a b^2)/(a^2 + b^2)$ +4. **Van Leer:** $\phi = 2ab/(a + b)$ if $ab > 0$ +5. **Superbee:** $\phi = \max(\min(2|a|,|b|),\;\min(|a|,2|b|))$ if $ab > 0$ + +where $a$ and $b$ are the left and right slope differences. + +**THINC interface compression** (`interface_compression > 0`): applies hyperbolic tangent compression near material interfaces: + +$$q_\text{THINC} = q_\min + \frac{q_\max}{2}\left(1 + \text{sign}(s)\,\frac{\tanh(\beta) + A}{1 + A\,\tanh(\beta)}\right)$$ + +where $A = \exp(\text{sign}(s)\,\beta\,(2C - 1)) / \cosh(\beta) - 1) / \tanh(\beta)$ and $\beta$ 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:** $\text{coeff}_L = [-3/60,\;27/60,\;47/60,\;-13/60,\;2/60]$ + +**3rd-order:** $\text{coeff}_L = [2/6,\;5/6,\;-1/6]$ + +**Ref:** [wilfong-arxiv-25] Algorithm 1. + +### 15.2 Riemann Solvers + +**Source:** `src/simulation/m_riemann_solvers.fpp` + +#### HLL (`riemann_solver = 1`) + +$$\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}$$ + +with wave speed estimates $S_L = \min(u_L - c_L,\;u_R - c_R)$, $S_R = \max(u_R + c_R,\;u_L + c_L)$. + +#### HLLC (`riemann_solver = 2`) + +Four-state solver resolving the contact discontinuity. Star-state satisfies: + +$$p_L^* = p_R^* = p^*, \qquad u_L^* = u_R^* = u^*$$ + +#### HLLD (`riemann_solver = 3`, MHD only) + +Seven-state solver for ideal MHD resolving fast magnetosonic, Alfven, and contact waves (Miyoshi and Kusano, 2005). The Riemann fan is divided by outer wave speeds $S_L$, $S_R$, Alfven speeds $S_L^*$, $S_R^*$, and a middle contact $S_M$: + +$$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}$$ + +$$S_L^* = S_M - \frac{|B_x|}{\sqrt{\rho_L^*}}, \qquad S_R^* = S_M + \frac{|B_x|}{\sqrt{\rho_R^*}}$$ + +where $p_T = p + |\mathbf{B}|^2/2$ is the total (thermal + magnetic) pressure. Continuity of normal velocity and total pressure is enforced across the Riemann fan. + +#### Exact (`riemann_solver = 4`) + +Iterative exact Riemann solver. + +### 15.3 Time Integration + +**Source:** `src/simulation/m_time_steppers.fpp` + +#### TVD Runge-Kutta (`time_stepper = 1, 2, 3`) + +**RK1 (Forward Euler):** + +$$\mathbf{q}^{n+1} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)$$ + +**RK2:** + +$$\mathbf{q}^{(1)} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)$$ + +$$\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)})$$ + +**RK3 (SSP):** + +$$\mathbf{q}^{(1)} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)$$ + +$$\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)})$$ + +$$\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)})$$ + +#### Adaptive Time Stepping (`adapt_dt = .true.`) + +Embedded RK pairs for error estimation with Hairer-Norsett-Wanner algorithm for initial step size. + +#### Strang Splitting (for stiff bubble dynamics) + +For equations of the form $\partial \mathbf{q}/\partial t = -\nabla \cdot \mathbf{F}(\mathbf{q}) + \mathbf{s}(\mathbf{q})$, the Strang splitting scheme integrates three sub-equations per time step: + +$$\frac{\partial \mathbf{q}'}{\partial t} = \mathbf{s}(\mathbf{q}'), \quad t \in [0,\,\Delta t/2], \quad \mathbf{q}'(0) = \mathbf{q}^n$$ + +$$\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)$$ + +$$\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)$$ + +The stiff bubble source terms are integrated using an adaptive embedded Runge-Kutta scheme with error control. + +### 15.4 CFL Condition + +**Inviscid:** + +$$\Delta t = \text{CFL} \cdot \min_{i,j,k}\left(\frac{\Delta x_\text{cell}}{|u| + |v| + |w| + c}\right)$$ + +**Viscous:** + +$$\Delta t_v = \text{CFL} \cdot \min\left(\frac{\Delta x^2}{\nu}\right)$$ + +where $c$ is the mixture sound speed and $\nu$ 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'_i = \frac{f_{i+1} - f_i}{x_{i+1} - x_i}$$ + +**2nd-order centered:** + +$$f'_i = \frac{f_{i+1} - f_{i-1}}{x_{i+1} - x_{i-1}}$$ + +**4th-order centered:** + +$$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}}$$ + +**Boundary-adjusted one-sided stencils (3rd-order):** + +$$f'_i\big|_\text{left} = \frac{-3f_i + 4f_{i+1} - f_{i+2}}{x_{i+2} - x_i}$$ + +$$f'_i\big|_\text{right} = \frac{3f_i - 4f_{i-1} + f_{i-2}}{x_i - x_{i-2}}$$ + +--- + +## 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 (Thompson/LODI, `bc_x%beg = -5` to `-12`) + +**Characteristic decomposition:** + +$$\frac{\partial \mathbf{q}_c}{\partial t} + \mathbf{R}_x\,\boldsymbol{\Lambda}_x\,\mathbf{L}_x\,\frac{\partial \mathbf{q}_c}{\partial x} = 0$$ + +**Characteristic amplitudes** $\mathcal{L}_1$ through $\mathcal{L}_5$ define wave interactions at boundaries: + +$$\mathcal{L}_1 = (u - c)\left(\frac{\partial p}{\partial x} - \rho\,c\,\frac{\partial u}{\partial x}\right), \qquad \mathcal{L}_2 = a\left(\frac{\partial \rho}{\partial x} + \frac{2}{\partial x}\frac{\partial p}{\partial x}\right)$$ + +$$\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)$$ + +For non-reflecting boundaries, the incoming wave amplitudes are set to zero. + +**GRCBC (incoming wave from ghost point):** + +$$\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}$$ + +**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`). + +**Ref:** [wilfong-mfc-26] Section 4.2.1. + +### 16.3 Immersed Boundary Method (`ib = .true.`) + +**Source:** `src/simulation/m_ibm.fpp` + +Ghost-cell IBM with level set field $\phi$. + +**Image point:** + +$$\mathbf{x}_{ip} = \mathbf{x}_{gp} + 2\,\phi(\mathbf{x}_{gp})\,\hat{\mathbf{n}}$$ + +**Velocity boundary conditions:** +- **No-slip:** $\mathbf{u}_{gp} = \mathbf{0}$ +- **Slip:** $\mathbf{u}_{gp} = \mathbf{u}_{ip} - (\hat{\mathbf{n}} \cdot \mathbf{u}_{ip})\,\hat{\mathbf{n}}$ + +**Neumann BC** for pressure, density, and volume fractions. + +Supports STL/OBJ geometry import with ray tracing for inside/outside determination. + +**Ref:** [wilfong-mfc-26] Section 4.1.1. + +--- + +## 17. Low Mach Number Corrections + +**Thornber correction** (`low_Mach = 1`): modifies the reconstructed velocities at cell interfaces: + +$$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}$$ + +$$z = \min\!\left(\max\!\left(\frac{|u_L|}{c_L},\;\frac{|u_R|}{c_R}\right),\;1\right)$$ + +This reduces numerical dissipation at low Mach numbers while recovering the standard scheme at high Mach numbers. + +**Chen correction** (`low_Mach = 2`): anti-dissipation pressure correction (APC) added to the HLLC flux: + +$$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)}$$ + +$$\text{APC}_{p\mathbf{u}} = (z - 1)\,p_d\,\hat{\mathbf{n}}, \qquad \text{APC}_{pE} = (z - 1)\,p_d\,S_*$$ + +The corrected HLLC flux in the star region becomes $\mathbf{F}^{*} + \text{APC}$. + +Additionally, the **artificial Mach number** technique (`pi_fac`) modifies $\pi_\infty$ to artificially increase the local Mach number. + +**Ref:** [wilfong-mfc-26] Section 4.2.4. + +--- + +## 18. Flux Limiting + +**Volume fraction limiting:** enforces $0 \le \alpha_k \le 1$ and rescales as: + +$$\alpha_k \leftarrow \frac{\alpha_k}{\sum_j \alpha_j}$$ + +**Advective flux limiting** based on local volume fraction gradient $\chi$ to prevent spurious oscillations near material interfaces. + +**Ref:** [bryngelson-CPC-20] Section 4.2. From 2845cfc29aaf54f86ab4eb15bc9b7f0093578eb0 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Thu, 12 Feb 2026 11:45:51 -0500 Subject: [PATCH 2/6] Integrate equations reference into Doxygen documentation Move docs/equations.md to docs/documentation/equations.md so it is picked up by the GEN_DOCS(documentation) CMake target. Add @page directive for Doxygen and link it from the main documentation hub under the Reference section. Co-Authored-By: Claude Opus 4.6 --- docs/{ => documentation}/equations.md | 2 ++ docs/documentation/readme.md | 1 + 2 files changed, 3 insertions(+) rename docs/{ => documentation}/equations.md (99%) diff --git a/docs/equations.md b/docs/documentation/equations.md similarity index 99% rename from docs/equations.md rename to docs/documentation/equations.md index 401ff80629..19cae537ed 100644 --- a/docs/equations.md +++ b/docs/documentation/equations.md @@ -1,3 +1,5 @@ +@page equations Equations + # MFC: Comprehensive Equations Reference This document catalogs every equation solved by MFC, organized by physical model. 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 From 5298d5ed8a516b5bcb3fc68d3f4fd2480a0c432b Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Thu, 12 Feb 2026 12:14:21 -0500 Subject: [PATCH 3/6] Fix Doxygen rendering: use native math syntax and clean up refs Convert all math delimiters from GitHub markdown dollar signs to Doxygen native (backslash-f) so that LaTeX commands like dot are passed verbatim to MathJax instead of being intercepted by Doxygen command parser (which was interpreting dot as a Graphviz directive and truncating the page). Also replace per-section paper citations with a single link to the papers page and fix bold markers around inline formulas. Co-Authored-By: Claude Opus 4.6 --- docs/documentation/equations.md | 447 ++++++++++++++------------------ 1 file changed, 197 insertions(+), 250 deletions(-) diff --git a/docs/documentation/equations.md b/docs/documentation/equations.md index 19cae537ed..1035457eab 100644 --- a/docs/documentation/equations.md +++ b/docs/documentation/equations.md @@ -3,16 +3,9 @@ # 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, cross-references the relevant source files, and cites the originating papers. +Each section notes the input parameter(s) that activate the corresponding physics module and cross-references the relevant source files. -**References:** -- **[bryngelson-CPC-20]** Bryngelson et al., "MFC: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver," *Computer Physics Communications* **266**, 107396 (2021). -- **[radhakrishnan-CPC-24]** Radhakrishnan et al., "Method for portable, scalable, and performant GPU-accelerated simulation of multiphase compressible flow," *Computer Physics Communications* **302**, 109238 (2024). -- **[wilfong-mfc-26]** Wilfong et al., "MFC 5.0: An exascale many-physics flow solver," *Computer Physics Communications* **322**, 110055 (2026). -- **[rodriguez-CAV-21]** Rodriguez et al., "Acoustically-induced bubble growth and phase change dynamics near compliant surfaces," *11th International Symposium on Cavitation* (2021). -- **[cisneros-cpc-25]** Cisneros-Garibay et al., "Pyrometheus: Symbolic abstractions for XPU and automatically differentiated computation of combustion kinetics and thermodynamics," *Computer Physics Communications* **320**, 109987 (2026). -- **[wilfong-igr-SC25]** Wilfong et al., "Information Geometric Regularization of the Barotropic Euler Equations," SC'25 (2025). -- **[wilfong-arxiv-25]** Wilfong et al., arxiv preprint (2025). +For full citations, see @ref papers. --- @@ -20,13 +13,13 @@ Each section notes the input parameter(s) that activate the corresponding physic 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: -$$\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[\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: -- $\mathbf{q}$ is the conservative variable vector, -- $\mathbf{F}$ is the flux tensor, -- $\mathbf{h}(\mathbf{q})\,\nabla \cdot \mathbf{u}$ contains non-conservative terms (volume fraction advection), -- $\mathbf{s}(\mathbf{q})$ is the source vector (bubbles, body forces, chemistry, etc.). +- \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. @@ -40,35 +33,33 @@ The parameter `model_eqns` (1, 2, 3, or 4) selects the governing equation set. The primary workhorse model (Allaire et al., 2002). The state vector is: -$$\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[\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): -$$\frac{\partial (\alpha_i \rho_i)}{\partial t} + \nabla \cdot (\alpha_i \rho_i\,\mathbf{u}) = 0$$ +\f[\frac{\partial (\alpha_i \rho_i)}{\partial t} + \nabla \cdot (\alpha_i \rho_i\,\mathbf{u}) = 0\f] **Momentum:** -$$\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[\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:** -$$\frac{\partial (\rho E)}{\partial t} + \nabla \cdot \bigl[(\rho E + p)\,\mathbf{u} - \boldsymbol{\tau}^v \cdot \mathbf{u}\bigr] = 0$$ +\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:** -$$\frac{\partial \alpha_i}{\partial t} + \mathbf{u} \cdot \nabla \alpha_i = K\,\nabla \cdot \mathbf{u}$$ +\f[\frac{\partial \alpha_i}{\partial t} + \mathbf{u} \cdot \nabla \alpha_i = K\,\nabla \cdot \mathbf{u}\f] -where the $K$ term enforces interface conditions via the Wood sound speed: +where the \f$K\f$ term enforces interface conditions via the Wood sound speed: -$$K = \frac{\rho_2 c_2^2 - \rho_1 c_1^2}{\dfrac{\rho_1 c_1^2}{\alpha_1} + \dfrac{\rho_2 c_2^2}{\alpha_2}}$$ +\f[K = \frac{\rho_2 c_2^2 - \rho_1 c_1^2}{\dfrac{\rho_1 c_1^2}{\alpha_1} + \dfrac{\rho_2 c_2^2}{\alpha_2}}\f] -Setting `alt_soundspeed = .true.` enables the $K$ correction (Kapila model with Wood sound speed). Setting `alt_soundspeed = .false.` uses the Allaire variant without the $K$ correction, which is conservative but does not strictly obey the second law of thermodynamics. +Setting `alt_soundspeed = .true.` enables the \f$K\f$ correction (Kapila model 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:** -$$1 = \sum_i \alpha_i, \qquad \rho = \sum_i \alpha_i \rho_i, \qquad \rho e = \sum_i \alpha_i \rho_i e_i$$ - -**Ref:** [wilfong-mfc-26] Section 2.1.1, [bryngelson-CPC-20] Section 3.1, [radhakrishnan-CPC-24] Section 2.1 Eq. 2. +\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`) @@ -78,34 +69,32 @@ Allows pressure disequilibrium between phases (Saurel et al., 2009). **Separate phasic internal energy:** -$$\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[\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:** -$$\frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 = \mu\,(p_1 - p_2)$$ +\f[\frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 = \mu\,(p_1 - p_2)\f] **Interfacial pressure:** -$$p_I = \frac{z_2\,p_1 + z_1\,p_2}{z_1 + z_2}, \qquad z_i = \rho_i\,c_i$$ +\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:** -$$c^2 = \sum_k Y_k\,c_k^2$$ +\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:** $\mu\,\delta p$ where $\delta p = p_1 - p_2$ -- **Thermal transfer:** $Q = \theta\,(T_2 - T_1)$ -- **Mass transfer:** $\dot{m} = \nu\,(g_2 - g_1)$ (Gibbs free energy difference) - -See [Section 8](#8-phase-change-relax--true) for details. +- **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) -**Ref:** [wilfong-mfc-26] Section 2.1.2, [bryngelson-CPC-20] Section 3.2, [rodriguez-CAV-21] Section 2.1, `src/simulation/m_pressure_relaxation.fpp`. +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 $\gamma$ and $\pi_\infty$ directly without tracking individual volume fractions (Johnsen, 2008). +- `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 (Johnsen, 2008). - `model_eqns = 4`: **Four-equation model** — reduced model from the six-equation system after full pressure-temperature equilibrium relaxation (Tait-like compressible liquid). --- @@ -116,59 +105,55 @@ See [Section 8](#8-phase-change-relax--true) for details. The primary closure for each phase: -$$p_k = (\gamma_k - 1)\,\rho_k\,e_k - \gamma_k\,\pi_{\infty,k}$$ +\f[p_k = (\gamma_k - 1)\,\rho_k\,e_k - \gamma_k\,\pi_{\infty,k}\f] Equivalently: -$$e_k = \frac{p_k + \gamma_k\,\pi_{\infty,k}}{(\gamma_k - 1)\,\rho_k}$$ +\f[e_k = \frac{p_k + \gamma_k\,\pi_{\infty,k}}{(\gamma_k - 1)\,\rho_k}\f] **Total energy relation:** -$$\rho E = \Gamma\,p + \Pi_\infty + \frac{1}{2}\rho\,|\mathbf{u}|^2 + q_v$$ +\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: -$$\Gamma_k = \frac{1}{\gamma_k - 1}, \qquad \Pi_{\infty,k} = \frac{\gamma_k\,\pi_{\infty,k}}{\gamma_k - 1}$$ +\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: -$$\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[\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: -$$p = \frac{\rho E - \frac{1}{2}\rho\,|\mathbf{u}|^2 - \Pi_\infty - q_v}{\Gamma}$$ +\f[p = \frac{\rho E - \frac{1}{2}\rho\,|\mathbf{u}|^2 - \Pi_\infty - q_v}{\Gamma}\f] **Phasic speed of sound:** -$$c_k = \sqrt{\frac{\gamma_k\,(p + \pi_{\infty,k})}{\rho_k}}$$ +\f[c_k = \sqrt{\frac{\gamma_k\,(p + \pi_{\infty,k})}{\rho_k}}\f] **Wood mixture sound speed:** -$$\frac{1}{\rho\,c^2} = \sum_k \frac{\alpha_k}{\rho_k\,c_k^2}$$ - -Input parameters per fluid: `gamma` ($\gamma_k$), `pi_inf` ($\pi_{\infty,k}$), `cv` ($c_{v,k}$), `qv` ($q_{v,k}$), `qvp` ($q'_{v,k}$). +\f[\frac{1}{\rho\,c^2} = \sum_k \frac{\alpha_k}{\rho_k\,c_k^2}\f] -**Ref:** [wilfong-mfc-26] Section 2.2, [radhakrishnan-CPC-24] Section 2.1 Eq. 1, `src/common/m_variables_conversion.fpp`. +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: -$$p = \frac{\rho\,R_u\,T}{W}, \qquad W = \left(\sum_m \frac{Y_m}{W_m}\right)^{-1}$$ +\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: -$$e_g - \sum_m e_m(T)\,Y_m = 0$$ +\f[e_g - \sum_m e_m(T)\,Y_m = 0\f] **Species internal energy from enthalpy:** -$$e_m(T) = \frac{\hat{h}_m(T) - R_u\,T}{W_m}$$ +\f[e_m(T) = \frac{\hat{h}_m(T) - R_u\,T}{W_m}\f] **NASA polynomial enthalpies:** -$$\frac{\hat{h}_m}{R_u\,T} = \frac{C_0}{T} + \sum_{r=1}^{5} \frac{C_r\,T^{r-1}}{r}$$ - -**Ref:** [wilfong-mfc-26] Section 4.1.7 Eqs. 18-22. +\f[\frac{\hat{h}_m}{R_u\,T} = \frac{C_0}{T} + \sum_{r=1}^{5} \frac{C_r\,T^{r-1}}{r}\f] --- @@ -176,45 +161,41 @@ $$\frac{\hat{h}_m}{R_u\,T} = \frac{C_0}{T} + \sum_{r=1}^{5} \frac{C_r\,T^{r-1}}{ **Newtonian viscous stress (no bulk viscosity by default):** -$$\boldsymbol{\tau}^v = 2\,\eta\left(\mathbf{D} - \frac{1}{3}\,\text{tr}(\mathbf{D})\,\mathbf{I}\right)$$ +\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: -$$\mathbf{D} = \frac{1}{2}\bigl(\nabla \mathbf{u} + (\nabla \mathbf{u})^T\bigr)$$ +\f[\mathbf{D} = \frac{1}{2}\bigl(\nabla \mathbf{u} + (\nabla \mathbf{u})^T\bigr)\f] **With bulk viscosity:** -$$\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[\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:** -$$\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[\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 $1/r$ terms. +and similarly for all other components. Cylindrical coordinate formulations include additional \f$1/r\f$ terms. **Viscosity averaging:** -$$\frac{1}{\text{Re}_\text{mix}} = \sum_j \frac{\alpha_j}{\text{Re}_j}$$ +\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). -**Ref:** [wilfong-mfc-26] Eqs. 1-2, `src/simulation/m_viscous.fpp`. - --- ## 5. Cylindrical Coordinates (`cyl_coord = .true.`) -Additional geometric source terms appear with $1/r$ factors in the continuity, momentum, and energy equations. Key modifications: +Additional geometric source terms appear with \f$1/r\f$ factors in the continuity, momentum, and energy equations. Key modifications: -- **Radial momentum:** extra $p/r$ and $\tau_{\theta\theta}/r$ terms -- **Viscous stress:** $\tau_{yy}$ includes $v/r$ corrections: +- **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: -$$\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[\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 -**Ref:** [wilfong-mfc-26] Section 2.3, [bryngelson-CPC-20] Section 4.3. - --- ## 6. Sub-Grid Bubble Dynamics @@ -227,112 +208,108 @@ $$\tau_{yy} = \mu\left(\frac{4}{3}\frac{\partial v}{\partial r} - \frac{2}{3}\fr **Modified mixture pressure:** -$$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[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:** -$$\Gamma_l\,p_l + \Pi_{\infty,l} = \frac{1}{1 - \alpha}\left(E - \frac{1}{2}\rho\,|\mathbf{u}|^2\right)$$ +\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):** -$$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[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:** -$$\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = \frac{3\,\alpha\,\bar{R}^2\,\dot{R}}{\bar{R}^3}$$ +\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:** -$$\frac{\partial n_\text{bub}}{\partial t} + \nabla \cdot (n_\text{bub}\,\mathbf{u}) = 0$$ +\f[\frac{\partial n_\text{bub}}{\partial t} + \nabla \cdot (n_\text{bub}\,\mathbf{u}) = 0\f] -where $n = \frac{3}{4\pi}\,\frac{\alpha}{\bar{R}^3}$. +where \f$n = \frac{3}{4\pi}\,\frac{\alpha}{\bar{R}^3}\f$. -**Polydispersity** (`polydisperse = .true.`): Log-normal PDF discretized into $N_\text{bin}$ equilibrium radii with standard deviation `poly_sigma`, integrated via Simpson's rule. +**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`) -$$R\,\ddot{R} + \frac{3}{2}\,\dot{R}^2 = \frac{p_{bw} - p_\infty}{\rho_l}$$ +\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`) -$$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[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`) Enthalpy-based formulation with compressibility corrections via the Tait EOS: -$$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[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: -$$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[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: -$$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[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.`) **Internal bubble pressure ODE:** -$$\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[\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:** -$$\dot{m}_v = \frac{D\,\rho_{bw}}{1 - \chi_{vw}}\left.\frac{\partial \chi_v}{\partial r}\right|_R$$ +\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.`) **Population balance equation:** -$$\frac{\partial f}{\partial t} + \frac{\partial (f\,\dot{R})}{\partial R} + \frac{\partial (f\,\ddot{R})}{\partial \dot{R}} = 0$$ +\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:** -$$\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[\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 $\mu_{i_1,i_2} = \int R^{i_1}\,\dot{R}^{i_2}\,f\,dR\,d\dot{R}$. +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 $(w_j, R_j, \dot{R}_j)$ from 6 moments via: +**CHyQMOM inversion** recovers 4 quadrature nodes \f$(w_j, R_j, \dot{R}_j)\f$ from 6 moments via: -$$\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$$ - -**Ref:** [wilfong-mfc-26] Section 4.1.4, `src/simulation/m_qbmm.fpp`. +\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.`) **Source:** `src/simulation/m_bubbles_EL.fpp` -Volume-averaged carrier flow equations with bubble source terms (Eq. 12 of [wilfong-mfc-26]): +Volume-averaged carrier flow equations with bubble source terms: **Continuity:** -$$\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[\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:** -$$\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[\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:** -$$\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[\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:** -$$\alpha(\mathbf{x}) = \sum_n V_n\,\delta_\sigma(\mathbf{x} - \mathbf{x}_n)$$ +\f[\alpha(\mathbf{x}) = \sum_n V_n\,\delta_\sigma(\mathbf{x} - \mathbf{x}_n)\f] -where $\delta_\sigma$ is a Gaussian kernel: +where \f$\delta_\sigma\f$ is a Gaussian kernel: -$$\delta_\sigma(\mathbf{r}) = \frac{1}{(2\pi\sigma^2)^{3/2}}\exp\!\left(-\frac{|\mathbf{r}|^2}{2\sigma^2}\right)$$ +\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 $\sigma = \varepsilon_b \max(\Delta x^{1/3}_\text{cell},\;R_\text{bubble})$. +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. -**Ref:** [wilfong-mfc-26] Section 4.1.5 Eqs. 13-14. - --- ## 7. Fluid-Structure Interaction @@ -343,29 +320,27 @@ Each bubble is tracked individually with Keller-Miksis dynamics and 4th-order ad **Cauchy stress decomposition:** -$$\sigma_{ij} = -p\,\delta_{ij} + \tau_{ij}^{(v)} + \tau_{ij}^{(e)}$$ +\f[\sigma_{ij} = -p\,\delta_{ij} + \tau_{ij}^{(v)} + \tau_{ij}^{(e)}\f] **Elastic energy contribution to total energy:** -$$E = e + \frac{|\mathbf{u}|^2}{2} + \frac{\boldsymbol{\tau}^e : \boldsymbol{\tau}^e}{4\,\rho\,G}$$ +\f[E = e + \frac{|\mathbf{u}|^2}{2} + \frac{\boldsymbol{\tau}^e : \boldsymbol{\tau}^e}{4\,\rho\,G}\f] **Elastic stress evolution:** -$$\frac{\partial (\rho\,\boldsymbol{\tau}^e)}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\tau}^e \otimes \mathbf{u}) = \mathbf{S}^e$$ +\f[\frac{\partial (\rho\,\boldsymbol{\tau}^e)}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\tau}^e \otimes \mathbf{u}) = \mathbf{S}^e\f] **Source term:** -$$\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[\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 $\mathbf{l} = \nabla \mathbf{u}$ is the velocity gradient and $\mathbf{D}^d = \mathbf{D} - \frac{1}{3}\text{tr}(\mathbf{D})\,\mathbf{I}$ is the deviatoric strain rate. +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):** -$$\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$$ - -This adds 6 additional transport equations in 3D (symmetric stress tensor: $\tau_{xx}^e, \tau_{xy}^e, \tau_{yy}^e, \tau_{xz}^e, \tau_{yz}^e, \tau_{zz}^e$). +\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] -**Ref:** [wilfong-mfc-26] Section 4.1.6, [radhakrishnan-CPC-24] Section 2.2 Eqs. 3-5. +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.`) @@ -373,27 +348,25 @@ This adds 6 additional transport equations in 3D (symmetric stress tensor: $\tau **Reference map evolution:** -$$\frac{\partial (\rho\,\boldsymbol{\xi})}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\xi} \otimes \mathbf{u}) = 0$$ +\f[\frac{\partial (\rho\,\boldsymbol{\xi})}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\xi} \otimes \mathbf{u}) = 0\f] **Deformation gradient from reference map:** -$$\mathbf{F} = (\nabla \boldsymbol{\xi})^{-1}$$ +\f[\mathbf{F} = (\nabla \boldsymbol{\xi})^{-1}\f] **Left Cauchy-Green tensor:** -$$\mathbf{b} = \mathbf{F}\,\mathbf{F}^T$$ +\f[\mathbf{b} = \mathbf{F}\,\mathbf{F}^T\f] **Neo-Hookean Cauchy stress:** -$$\boldsymbol{\tau}^e = \frac{G}{J}\left(\mathbf{b} - \frac{\text{tr}(\mathbf{b})}{3}\,\mathbf{I}\right)$$ +\f[\boldsymbol{\tau}^e = \frac{G}{J}\left(\mathbf{b} - \frac{\text{tr}(\mathbf{b})}{3}\,\mathbf{I}\right)\f] -where $J = \det(\mathbf{F})$. +where \f$J = \det(\mathbf{F})\f$. **Hyperelastic energy:** -$$e^e = \frac{G}{2}\bigl(I_{\mathbf{b}} - 3\bigr), \qquad I_{\mathbf{b}} = \text{tr}(\mathbf{b})$$ - -**Ref:** [wilfong-mfc-26] Section 4.1.6. +\f[e^e = \frac{G}{2}\bigl(I_{\mathbf{b}} - 3\bigr), \qquad I_{\mathbf{b}} = \text{tr}(\mathbf{b})\f] --- @@ -403,38 +376,36 @@ $$e^e = \frac{G}{2}\bigl(I_{\mathbf{b}} - 3\bigr), \qquad I_{\mathbf{b}} = \text ### 8.1 pT-Relaxation (`relax_model = 5`) -$N$-fluid pressure-temperature equilibrium. The equilibrium condition is: +\f$N\f$-fluid pressure-temperature equilibrium. The equilibrium condition is: -$$f(p) = \sum_i \alpha_i - 1 = 0$$ +\f[f(p) = \sum_i \alpha_i - 1 = 0\f] **Temperature from energy conservation:** -$$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[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:** -$$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[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`) -Two coupled equations for $(\alpha_1 \rho_1,\;p)$: +Two coupled equations for \f$(\alpha_1 \rho_1,\;p)\f$: **Gibbs free energy equilibrium (Clausius-Clapeyron):** -$$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[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_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[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 $m_T$ is the total mass, $m_l = \alpha_l \rho_l$ is the liquid partial density, and $m_{qD}$, $m_{cpD}$, $m_{cvgp}$ are auxiliary thermodynamic sums over additional fluids (beyond the phase-changing pair). +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. -**Ref:** [wilfong-mfc-26] Section 4.1.3, [rodriguez-CAV-21] Section 2. - --- ## 9. Chemistry and Combustion (`chemistry = .true.`) @@ -443,32 +414,30 @@ Solved via 2D Newton-Raphson. **Species transport:** -$$\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[\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:** -$$\dot{\omega}_m = \sum_n (\nu''_{mn} - \nu'_{mn})\,\mathcal{R}_n$$ +\f[\dot{\omega}_m = \sum_n (\nu''_{mn} - \nu'_{mn})\,\mathcal{R}_n\f] **Reaction rate (law of mass action):** -$$\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[\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:** -$$k_n(T) = A_n\,T^{b_n}\exp\!\left(-\frac{T_{a,n}}{T}\right)$$ +\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 $D_m^\text{mix}$, mass flux: $\dot{m}_k = \rho\,D_k^\text{mix}\,(W_k / W_\text{mix})\,\partial X_k / \partial x$ -- **Unity Lewis number:** $D_m = \lambda / (\rho\,c_p)$ +- **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: -$$q_\text{diff} = \lambda\,\frac{\partial T}{\partial x} + \sum_k h_k\,\dot{m}_k$$ +\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, which provides symbolic abstractions for thermochemistry that enable portable GPU computation and automatic differentiation of chemical source terms. -**Ref:** [wilfong-mfc-26] Section 4.1.7 Eqs. 15-22, [cisneros-cpc-25] Section 5. - --- ## 10. Surface Tension (`surface_tension = .true.`) @@ -477,25 +446,23 @@ Reaction mechanisms are code-generated via Pyrometheus, which provides symbolic **Color function advection:** -$$\frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla c = 0$$ +\f[\frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla c = 0\f] **Capillary stress tensor (CSF model):** -$$\boldsymbol{\Omega} = -\sigma\left(\|\nabla c\|\,\mathbf{I} - \frac{\nabla c \otimes \nabla c}{\|\nabla c\|}\right)$$ +\f[\boldsymbol{\Omega} = -\sigma\left(\|\nabla c\|\,\mathbf{I} - \frac{\nabla c \otimes \nabla c}{\|\nabla c\|}\right)\f] -In component form, with $\hat{w}_i = (\partial c / \partial x_i) / \|\nabla c\|$: +In component form, with \f$\hat{w}_i = (\partial c / \partial x_i) / \|\nabla c\|\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[\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: -$$\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[\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:** -$$\varepsilon_0 = \sigma\,\|\nabla c\|$$ - -**Ref:** [wilfong-mfc-26] Section 4.1.8. +\f[\varepsilon_0 = \sigma\,\|\nabla c\|\f] --- @@ -505,80 +472,74 @@ $$\varepsilon_0 = \sigma\,\|\nabla c\|$$ **Continuity:** -$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho\,\mathbf{u}) = 0$$ +\f[\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho\,\mathbf{u}) = 0\f] **Momentum:** -$$\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[\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:** -$$\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[\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:** -$$\frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u}) = 0$$ +\f[\frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u}) = 0\f] **Total energy:** -$$\mathcal{E} = \rho\,e + \frac{1}{2}\rho\,|\mathbf{u}|^2 + \frac{|\mathbf{B}|^2}{2}$$ +\f[\mathcal{E} = \rho\,e + \frac{1}{2}\rho\,|\mathbf{u}|^2 + \frac{|\mathbf{B}|^2}{2}\f] **Fast magnetosonic speed:** -$$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[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:** -$$v_A = \sqrt{\frac{|\mathbf{B}|^2}{\rho}}$$ +\f[v_A = \sqrt{\frac{|\mathbf{B}|^2}{\rho}}\f] Uses the HLLD Riemann solver (`riemann_solver = 3`). Hyperbolic divergence cleaning (`hyper_cleaning = .true.`) via the GLM method (Dedner et al., 2002). -**Ref:** [wilfong-mfc-26] Section 4.1.9. - ### 11.2 Relativistic MHD (`relativity = .true.`) **Conserved variables:** -$$\mathbf{U} = (D,\;\mathbf{m},\;\tau,\;\mathbf{B})^T$$ +\f[\mathbf{U} = (D,\;\mathbf{m},\;\tau,\;\mathbf{B})^T\f] where: -$$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[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] -$$\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[\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. -**Ref:** [wilfong-mfc-26] Section 4.1.10. - --- ## 12. Information Geometric Regularization (`igr = .true.`) **Source:** `src/simulation/m_igr.fpp` -**Modified momentum with entropic pressure $\Sigma$:** +**Modified momentum with entropic pressure** \f$\Sigma\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[\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 $\Sigma$:** +**Elliptic PDE for** \f$\Sigma\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[\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 $\alpha \sim \Delta x^2$ (regularization strength proportional to mesh spacing squared): +where \f$\alpha \sim \Delta x^2\f$ (regularization strength proportional to mesh spacing squared): -$$\alpha_\text{IGR} = \alpha_\text{factor} \cdot \max(\Delta x,\;\Delta y,\;\Delta z)^2$$ +\f[\alpha_\text{IGR} = \alpha_\text{factor} \cdot \max(\Delta x,\;\Delta y,\;\Delta z)^2\f] **RHS strain-rate source (3D):** -$$\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[\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). -**Ref:** [wilfong-arxiv-25] Eqs. 6-9. - --- ## 13. Body Forces (`bodyforces = .true.`) @@ -587,15 +548,15 @@ Uses Lax-Friedrichs flux (replaces WENO + Riemann solver). **Time-dependent acceleration:** -$$a_i(t) = g_i + k_i\sin(\omega_i\,t - \phi_i)$$ +\f[a_i(t) = g_i + k_i\sin(\omega_i\,t - \phi_i)\f] **Momentum source:** -$$\frac{\partial (\rho\,u_i)}{\partial t}\bigg|_\text{bf} = \rho\,a_i(t)$$ +\f[\frac{\partial (\rho\,u_i)}{\partial t}\bigg|_\text{bf} = \rho\,a_i(t)\f] **Energy source:** -$$\frac{\partial E}{\partial t}\bigg|_\text{bf} = \rho\,\mathbf{u} \cdot \mathbf{a}(t)$$ +\f[\frac{\partial E}{\partial t}\bigg|_\text{bf} = \rho\,\mathbf{u} \cdot \mathbf{a}(t)\f] --- @@ -607,22 +568,20 @@ Source terms added to the RHS of the governing equations. **Discrete delta function (spatial support):** -$$\delta_h(r) = \frac{1}{(2\pi\sigma^2)^{d/2}}\exp\!\left(-\frac{r^2}{2\sigma^2}\right)$$ +\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):** -$$\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[\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`): $S(t) = M\sin(\omega(t - t_\text{delay}))$ -- **Gaussian** (`pulse = 2`): $S(t) = M\exp\!\bigl(-\tfrac{1}{2}((t - t_\text{delay})/\sigma_t)^2\bigr)$ -- **Square** (`pulse = 3`): $S(t) = M\,\text{sign}(\sin(\omega(t - t_\text{delay})))$ +- **Sine** (`pulse = 1`): \f$S(t) = M\sin(\omega(t - t_\text{delay}))\f$ +- **Gaussian** (`pulse = 2`): \f$S(t) = M\exp\!\bigl(-\tfrac{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). -**Ref:** [bryngelson-CPC-20] Section 4.4 Eqs. 35-38. - --- ## 15. Numerical Methods @@ -635,64 +594,60 @@ $$\mathbf{s}_\text{ac} = \Omega_\Gamma\,f(t)\left[\frac{1}{c_0},\;\cos\theta,\;\ Weighted sum of candidate polynomials at cell interfaces: -$$f_{i+1/2} = \sum_r \omega_r\,f_{i+1/2}^{(r)}$$ +\f[f_{i+1/2} = \sum_r \omega_r\,f_{i+1/2}^{(r)}\f] **WENO-JS** (default): -$$\alpha_r = \frac{d_r}{(\beta_r + \varepsilon)^2}, \qquad \omega_r = \frac{\alpha_r}{\sum_s \alpha_s}$$ +\f[\alpha_r = \frac{d_r}{(\beta_r + \varepsilon)^2}, \qquad \omega_r = \frac{\alpha_r}{\sum_s \alpha_s}\f] -where $d_r$ are ideal weights, $\beta_r$ are smoothness indicators, and $\varepsilon$ is a small regularization parameter (`weno_eps`). +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.`): Henrick et al. (2005) mapped weights for improved accuracy at critical points: -$$\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[\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.`): Borges et al. (2008) improved weights with global smoothness measure: -$$\alpha_r = d_r\left(1 + \left(\frac{\tau}{\beta_r + \varepsilon}\right)^q\right), \qquad \tau = |\beta_0 - \beta_{k-1}|$$ +\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 $q$ controls the convergence rate at critical points (typically $q = 1$ for fifth-order reconstruction, as used in MFC). +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.`): Fu et al. (2016) targeted ENO with smoothness threshold $C_T$ (`teno_CT`): +**TENO** (`teno = .true.`): Fu et al. (2016) targeted ENO with smoothness threshold \f$C_T\f$ (`teno_CT`): -$$\gamma_r = 1 + \frac{\tau}{\beta_r}, \qquad \xi_r = \frac{\gamma_r}{\sum_s \gamma_s}$$ +\f[\gamma_r = 1 + \frac{\tau}{\beta_r}, \qquad \xi_r = \frac{\gamma_r}{\sum_s \gamma_s}\f] -If $\xi_r < C_T$, set $\alpha_r = 0$ (stencil excluded). +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. -**Ref:** [wilfong-mfc-26] Sections 3.2.1, 4.2.2. - #### MUSCL (`muscl_order = 2`) **Source:** `src/simulation/m_muscl.fpp` -$$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[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:** $\phi = \text{sign}(a)\,\min(|a|,\;|b|)$ if $ab > 0$, else $0$ -2. **MC (Monotonized Central):** $\phi = \text{sign}(a)\,\min(2|a|,\;2|b|,\;\tfrac{1}{2}(|a|+|b|))$ if $ab > 0$ -3. **OSPRE (Van Albada):** $\phi = (a^2 b + a b^2)/(a^2 + b^2)$ -4. **Van Leer:** $\phi = 2ab/(a + b)$ if $ab > 0$ -5. **Superbee:** $\phi = \max(\min(2|a|,|b|),\;\min(|a|,2|b|))$ if $ab > 0$ +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|,\;\tfrac{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 $a$ and $b$ are the left and right slope differences. +where \f$a\f$ and \f$b\f$ are the left and right slope differences. **THINC interface compression** (`interface_compression > 0`): applies hyperbolic tangent compression near material interfaces: -$$q_\text{THINC} = q_\min + \frac{q_\max}{2}\left(1 + \text{sign}(s)\,\frac{\tanh(\beta) + A}{1 + A\,\tanh(\beta)}\right)$$ +\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 $A = \exp(\text{sign}(s)\,\beta\,(2C - 1)) / \cosh(\beta) - 1) / \tanh(\beta)$ and $\beta$ controls compression steepness. +where \f$A = \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:** $\text{coeff}_L = [-3/60,\;27/60,\;47/60,\;-13/60,\;2/60]$ +**5th-order:** \f$\text{coeff}_L = [-3/60,\;27/60,\;47/60,\;-13/60,\;2/60]\f$ -**3rd-order:** $\text{coeff}_L = [2/6,\;5/6,\;-1/6]$ - -**Ref:** [wilfong-arxiv-25] Algorithm 1. +**3rd-order:** \f$\text{coeff}_L = [2/6,\;5/6,\;-1/6]\f$ ### 15.2 Riemann Solvers @@ -700,25 +655,25 @@ where $A = \exp(\text{sign}(s)\,\beta\,(2C - 1)) / \cosh(\beta) - 1) / \tanh(\be #### HLL (`riemann_solver = 1`) -$$\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[\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 $S_L = \min(u_L - c_L,\;u_R - c_R)$, $S_R = \max(u_R + c_R,\;u_L + c_L)$. +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`) Four-state solver resolving the contact discontinuity. Star-state satisfies: -$$p_L^* = p_R^* = p^*, \qquad u_L^* = u_R^* = u^*$$ +\f[p_L^* = p_R^* = p^*, \qquad u_L^* = u_R^* = u^*\f] #### HLLD (`riemann_solver = 3`, MHD only) -Seven-state solver for ideal MHD resolving fast magnetosonic, Alfven, and contact waves (Miyoshi and Kusano, 2005). The Riemann fan is divided by outer wave speeds $S_L$, $S_R$, Alfven speeds $S_L^*$, $S_R^*$, and a middle contact $S_M$: +Seven-state solver for ideal MHD resolving fast magnetosonic, Alfven, and contact waves (Miyoshi and Kusano, 2005). 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$: -$$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[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] -$$S_L^* = S_M - \frac{|B_x|}{\sqrt{\rho_L^*}}, \qquad S_R^* = S_M + \frac{|B_x|}{\sqrt{\rho_R^*}}$$ +\f[S_L^* = S_M - \frac{|B_x|}{\sqrt{\rho_L^*}}, \qquad S_R^* = S_M + \frac{|B_x|}{\sqrt{\rho_R^*}}\f] -where $p_T = p + |\mathbf{B}|^2/2$ is the total (thermal + magnetic) pressure. Continuity of normal velocity and total pressure is enforced across the Riemann fan. +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. #### Exact (`riemann_solver = 4`) @@ -732,21 +687,21 @@ Iterative exact Riemann solver. **RK1 (Forward Euler):** -$$\mathbf{q}^{n+1} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)$$ +\f[\mathbf{q}^{n+1} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)\f] **RK2:** -$$\mathbf{q}^{(1)} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)$$ +\f[\mathbf{q}^{(1)} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)\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[\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):** -$$\mathbf{q}^{(1)} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)$$ +\f[\mathbf{q}^{(1)} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)\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[\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] -$$\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[\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 (`adapt_dt = .true.`) @@ -754,13 +709,13 @@ Embedded RK pairs for error estimation with Hairer-Norsett-Wanner algorithm for #### Strang Splitting (for stiff bubble dynamics) -For equations of the form $\partial \mathbf{q}/\partial t = -\nabla \cdot \mathbf{F}(\mathbf{q}) + \mathbf{s}(\mathbf{q})$, the Strang splitting scheme integrates three sub-equations per time step: +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: -$$\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[\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] -$$\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[\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] -$$\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[\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. @@ -768,13 +723,13 @@ The stiff bubble source terms are integrated using an adaptive embedded Runge-Ku **Inviscid:** -$$\Delta t = \text{CFL} \cdot \min_{i,j,k}\left(\frac{\Delta x_\text{cell}}{|u| + |v| + |w| + c}\right)$$ +\f[\Delta t = \text{CFL} \cdot \min_{i,j,k}\left(\frac{\Delta x_\text{cell}}{|u| + |v| + |w| + c}\right)\f] **Viscous:** -$$\Delta t_v = \text{CFL} \cdot \min\left(\frac{\Delta x^2}{\nu}\right)$$ +\f[\Delta t_v = \text{CFL} \cdot \min\left(\frac{\Delta x^2}{\nu}\right)\f] -where $c$ is the mixture sound speed and $\nu$ is the kinematic viscosity. +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`) @@ -784,21 +739,21 @@ Used for viscous fluxes and velocity gradients. **1st-order:** -$$f'_i = \frac{f_{i+1} - f_i}{x_{i+1} - x_i}$$ +\f[f'_i = \frac{f_{i+1} - f_i}{x_{i+1} - x_i}\f] **2nd-order centered:** -$$f'_i = \frac{f_{i+1} - f_{i-1}}{x_{i+1} - x_{i-1}}$$ +\f[f'_i = \frac{f_{i+1} - f_{i-1}}{x_{i+1} - x_{i-1}}\f] **4th-order centered:** -$$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[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'_i\big|_\text{left} = \frac{-3f_i + 4f_{i+1} - f_{i+2}}{x_{i+2} - x_i}$$ +\f[f'_i\big|_\text{left} = \frac{-3f_i + 4f_{i+1} - f_{i+2}}{x_{i+2} - x_i}\f] -$$f'_i\big|_\text{right} = \frac{3f_i - 4f_{i-1} + f_{i-2}}{x_i - x_{i-2}}$$ +\f[f'_i\big|_\text{right} = \frac{3f_i - 4f_{i-1} + f_{i-2}}{x_i - x_{i-2}}\f] --- @@ -822,76 +777,68 @@ $$f'_i\big|_\text{right} = \frac{3f_i - 4f_{i-1} + f_{i-2}}{x_i - x_{i-2}}$$ **Characteristic decomposition:** -$$\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[\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** $\mathcal{L}_1$ through $\mathcal{L}_5$ define wave interactions at boundaries: +**Characteristic amplitudes** \f$\mathcal{L}_1\f$ through \f$\mathcal{L}_5\f$ define wave interactions at boundaries: -$$\mathcal{L}_1 = (u - c)\left(\frac{\partial p}{\partial x} - \rho\,c\,\frac{\partial u}{\partial x}\right), \qquad \mathcal{L}_2 = a\left(\frac{\partial \rho}{\partial x} + \frac{2}{\partial x}\frac{\partial p}{\partial x}\right)$$ +\f[\mathcal{L}_1 = (u - c)\left(\frac{\partial p}{\partial x} - \rho\,c\,\frac{\partial u}{\partial x}\right), \qquad \mathcal{L}_2 = a\left(\frac{\partial \rho}{\partial x} + \frac{2}{\partial x}\frac{\partial p}{\partial x}\right)\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[\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):** -$$\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[\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`). -**Ref:** [wilfong-mfc-26] Section 4.2.1. - ### 16.3 Immersed Boundary Method (`ib = .true.`) **Source:** `src/simulation/m_ibm.fpp` -Ghost-cell IBM with level set field $\phi$. +Ghost-cell IBM with level set field \f$\phi\f$. **Image point:** -$$\mathbf{x}_{ip} = \mathbf{x}_{gp} + 2\,\phi(\mathbf{x}_{gp})\,\hat{\mathbf{n}}$$ +\f[\mathbf{x}_{ip} = \mathbf{x}_{gp} + 2\,\phi(\mathbf{x}_{gp})\,\hat{\mathbf{n}}\f] **Velocity boundary conditions:** -- **No-slip:** $\mathbf{u}_{gp} = \mathbf{0}$ -- **Slip:** $\mathbf{u}_{gp} = \mathbf{u}_{ip} - (\hat{\mathbf{n}} \cdot \mathbf{u}_{ip})\,\hat{\mathbf{n}}$ +- **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. -**Ref:** [wilfong-mfc-26] Section 4.1.1. - --- ## 17. Low Mach Number Corrections **Thornber correction** (`low_Mach = 1`): modifies the reconstructed velocities at cell interfaces: -$$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[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] -$$z = \min\!\left(\max\!\left(\frac{|u_L|}{c_L},\;\frac{|u_R|}{c_R}\right),\;1\right)$$ +\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. **Chen correction** (`low_Mach = 2`): anti-dissipation pressure correction (APC) added to the HLLC flux: -$$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[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] -$$\text{APC}_{p\mathbf{u}} = (z - 1)\,p_d\,\hat{\mathbf{n}}, \qquad \text{APC}_{pE} = (z - 1)\,p_d\,S_*$$ +\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 $\mathbf{F}^{*} + \text{APC}$. +The corrected HLLC flux in the star region becomes \f$\mathbf{F}^{*} + \text{APC}\f$. -Additionally, the **artificial Mach number** technique (`pi_fac`) modifies $\pi_\infty$ to artificially increase the local Mach number. - -**Ref:** [wilfong-mfc-26] Section 4.2.4. +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 $0 \le \alpha_k \le 1$ and rescales as: - -$$\alpha_k \leftarrow \frac{\alpha_k}{\sum_j \alpha_j}$$ +**Volume fraction limiting:** enforces \f$0 \le \alpha_k \le 1\f$ and rescales as: -**Advective flux limiting** based on local volume fraction gradient $\chi$ to prevent spurious oscillations near material interfaces. +\f[\alpha_k \leftarrow \frac{\alpha_k}{\sum_j \alpha_j}\f] -**Ref:** [bryngelson-CPC-20] Section 4.2. +**Advective flux limiting** based on local volume fraction gradient \f$\chi\f$ to prevent spurious oscillations near material interfaces. From 3a8f9e3ffc691709b5fe69bebb6a810a8b2be239 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Thu, 12 Feb 2026 14:05:51 -0500 Subject: [PATCH 4/6] Switch documentation citations to Doxygen BibTeX cite system Replace manual HTML anchor references and inline text citations across all docs with Doxygen native cite command backed by a central references.bib file. This auto-generates a Bibliography page and provides clickable citation links throughout the documentation. Co-Authored-By: Claude Opus 4.6 --- .typos.toml | 2 +- CMakeLists.txt | 1 + docs/Doxyfile.in | 2 +- docs/documentation/case.md | 68 ++-- docs/documentation/equations.md | 79 ++-- docs/documentation/references.md | 72 +--- docs/documentation/visualization.md | 2 +- docs/references.bib | 592 ++++++++++++++++++++++++++++ 8 files changed, 673 insertions(+), 145 deletions(-) create mode 100644 docs/references.bib 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..f803027c55 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}\"") 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..1e2a5498c7 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 Toro13) | | `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 Bryngelson19. +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 Schmidmayer19). 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 Schmidmayer19. - `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 Toro13). +`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 Toro13), 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 index 1035457eab..7fb2635f20 100644 --- a/docs/documentation/equations.md +++ b/docs/documentation/equations.md @@ -5,7 +5,7 @@ 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. -For full citations, see @ref papers. +For full citations of MFC papers, see @ref papers. Foundational references for each model are cited inline; see the \ref citelist "Bibliography" for full details. --- @@ -31,7 +31,7 @@ The parameter `model_eqns` (1, 2, 3, or 4) selects the governing equation set. ### 2.1 Five-Equation Model (`model_eqns = 2`) -The primary workhorse model (Allaire et al., 2002). The state vector is: +The primary workhorse model (\cite Allaire02). 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] @@ -53,9 +53,9 @@ The primary workhorse model (Allaire et al., 2002). The state vector is: 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}{\dfrac{\rho_1 c_1^2}{\alpha_1} + \dfrac{\rho_2 c_2^2}{\alpha_2}}\f] +\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 (Kapila model 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. +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:** @@ -63,7 +63,7 @@ Setting `alt_soundspeed = .true.` enables the \f$K\f$ correction (Kapila model w ### 2.2 Six-Equation Model (`model_eqns = 3`) -Allows pressure disequilibrium between phases (Saurel et al., 2009). +Allows pressure disequilibrium between phases (\cite Saurel09). **Continuity and momentum:** Same as the five-equation model. @@ -94,14 +94,14 @@ 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 (Johnsen, 2008). +- `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 +### 3.1 Stiffened Gas EOS (\cite Menikoff89; \cite LeMetayer04) The primary closure for each phase: @@ -151,7 +151,7 @@ Temperature is obtained from the internal energy by Newton iteration: \f[e_m(T) = \frac{\hat{h}_m(T) - R_u\,T}{W_m}\f] -**NASA polynomial enthalpies:** +**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] @@ -204,7 +204,7 @@ Additional geometric source terms appear with \f$1/r\f$ factors in the continuit **Source:** `src/simulation/m_bubbles_EE.fpp`, `src/simulation/m_bubbles.fpp` -#### 6.1.1 Method of Classes +#### 6.1.1 Method of Classes (\cite Commander89; \cite Ando11) **Modified mixture pressure:** @@ -230,15 +230,15 @@ 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`) +#### 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`) +#### 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`) +#### 6.1.4 Gilmore (`bubble_model = 1`) (\cite Gilmore52) Enthalpy-based formulation with compressibility corrections via the Tait EOS: @@ -252,7 +252,7 @@ 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.`) +#### 6.1.5 Non-Polytropic Thermal Model (`polytropic = .false.`) (\cite Preston07) **Internal bubble pressure ODE:** @@ -262,7 +262,7 @@ and the local liquid sound speed: \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.`) +#### 6.1.6 QBMM Moment Transport (`qbmm = .true.`) (\cite Bryngelson23) **Population balance equation:** @@ -278,7 +278,7 @@ where moments \f$\mu_{i_1,i_2} = \int R^{i_1}\,\dot{R}^{i_2}\,f\,dR\,d\dot{R}\f$ \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.`) +### 6.2 Euler-Lagrange Bubbles (`bubbles_lagrange = .true.`) (\cite Maeda18) **Source:** `src/simulation/m_bubbles_EL.fpp` @@ -314,7 +314,7 @@ Each bubble is tracked individually with Keller-Miksis dynamics and 4th-order ad ## 7. Fluid-Structure Interaction -### 7.1 Hypoelastic Model (`hypoelasticity = .true.`) +### 7.1 Hypoelastic Model (`hypoelasticity = .true.`) (\cite Rodriguez19) **Source:** `src/simulation/m_hypoelastic.fpp` @@ -342,7 +342,7 @@ where \f$\mathbf{l} = \nabla \mathbf{u}\f$ is the velocity gradient and \f$\math 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.`) +### 7.2 Hyperelastic Model (`hyperelasticity = .true.`) (\cite Kamrin12) **Source:** `src/simulation/m_hyperelastic.fpp` @@ -374,7 +374,7 @@ where \f$J = \det(\mathbf{F})\f$. **Source:** `src/common/m_phase_change.fpp` -### 8.1 pT-Relaxation (`relax_model = 5`) +### 8.1 pT-Relaxation (`relax_model = 5`) (\cite Saurel08) \f$N\f$-fluid pressure-temperature equilibrium. The equilibrium condition is: @@ -390,7 +390,7 @@ where \f$J = \det(\mathbf{F})\f$. Solved via Newton's method for the equilibrium pressure. -### 8.2 pTg-Relaxation (`relax_model = 6`) +### 8.2 pTg-Relaxation (`relax_model = 6`) (\cite Zein10) Two coupled equations for \f$(\alpha_1 \rho_1,\;p)\f$: @@ -436,11 +436,11 @@ 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, which provides symbolic abstractions for thermochemistry that enable portable GPU computation and automatic differentiation of chemical source terms. +Reaction mechanisms are code-generated via Pyrometheus (\cite Cisneros25), which provides symbolic abstractions for thermochemistry that enable portable GPU computation and automatic differentiation of chemical source terms. --- -## 10. Surface Tension (`surface_tension = .true.`) +## 10. Surface Tension (`surface_tension = .true.`) (\cite Schmidmayer17) **Source:** `src/simulation/m_surface_tension.fpp`, `src/simulation/include/inline_capillary.fpp` @@ -498,7 +498,7 @@ The capillary stress divergence is added to the momentum and energy equations. T \f[v_A = \sqrt{\frac{|\mathbf{B}|^2}{\rho}}\f] -Uses the HLLD Riemann solver (`riemann_solver = 3`). Hyperbolic divergence cleaning (`hyper_cleaning = .true.`) via the GLM method (Dedner et al., 2002). +Uses the HLLD Riemann solver (`riemann_solver = 3`). Hyperbolic divergence cleaning (`hyper_cleaning = .true.`) via the GLM method (\cite Dedner02). ### 11.2 Relativistic MHD (`relativity = .true.`) @@ -516,7 +516,7 @@ Primitive recovery uses Newton-Raphson on the nonlinear conserved-to-primitive r --- -## 12. Information Geometric Regularization (`igr = .true.`) +## 12. Information Geometric Regularization (`igr = .true.`) (\cite Wilfong25a) **Source:** `src/simulation/m_igr.fpp` @@ -566,6 +566,8 @@ Uses Lax-Friedrichs flux (replaces WENO + Riemann solver). 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] @@ -576,7 +578,7 @@ Source terms added to the RHS of the governing equations. **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(-\tfrac{1}{2}((t - t_\text{delay})/\sigma_t)^2\bigr)\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 @@ -596,23 +598,23 @@ 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** (default): +**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.`): Henrick et al. (2005) mapped weights for improved accuracy at critical points: +**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.`): Borges et al. (2008) improved weights with global smoothness measure: +**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.`): Fu et al. (2016) targeted ENO with smoothness threshold \f$C_T\f$ (`teno_CT`): +**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] @@ -628,7 +630,7 @@ Primitive variable reconstruction is used to avoid spurious oscillations at inte **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|,\;\tfrac{1}{2}(|a|+|b|))\f$ if \f$ab > 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$ @@ -653,13 +655,13 @@ where \f$A = \exp(\text{sign}(s)\,\beta\,(2C - 1)) / \cosh(\beta) - 1) / \tanh(\ **Source:** `src/simulation/m_riemann_solvers.fpp` -#### HLL (`riemann_solver = 1`) +#### 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`) +#### HLLC (`riemann_solver = 2`) (\cite Toro94) Four-state solver resolving the contact discontinuity. Star-state satisfies: @@ -667,7 +669,7 @@ Four-state solver resolving the contact discontinuity. Star-state satisfies: #### HLLD (`riemann_solver = 3`, MHD only) -Seven-state solver for ideal MHD resolving fast magnetosonic, Alfven, and contact waves (Miyoshi and Kusano, 2005). 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$: +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] @@ -683,7 +685,7 @@ Iterative exact Riemann solver. **Source:** `src/simulation/m_time_steppers.fpp` -#### TVD Runge-Kutta (`time_stepper = 1, 2, 3`) +#### TVD Runge-Kutta (`time_stepper = 1, 2, 3`) (\cite Gottlieb98) **RK1 (Forward Euler):** @@ -707,7 +709,7 @@ Iterative exact Riemann solver. Embedded RK pairs for error estimation with Hairer-Norsett-Wanner algorithm for initial step size. -#### Strang Splitting (for stiff bubble dynamics) +#### 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: @@ -773,7 +775,7 @@ Used for viscous fluxes and velocity gradients. | `-15` | Slip wall | | `-16` | No-slip wall | -### 16.2 Characteristic BCs (Thompson/LODI, `bc_x%beg = -5` to `-12`) +### 16.2 Characteristic BCs (\cite Thompson87, \cite Thompson90; `bc_x%beg = -5` to `-12`) **Characteristic decomposition:** @@ -793,7 +795,7 @@ For non-reflecting boundaries, the incoming wave amplitudes are set to zero. **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.`) +### 16.3 Immersed Boundary Method (`ib = .true.`) (\cite Tseng03; \cite Mittal05) **Source:** `src/simulation/m_ibm.fpp` @@ -815,7 +817,7 @@ Supports STL/OBJ geometry import with ray tracing for inside/outside determinati ## 17. Low Mach Number Corrections -**Thornber correction** (`low_Mach = 1`): modifies the reconstructed velocities at cell interfaces: +**Thornber correction** (`low_Mach = 1`, \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] @@ -823,7 +825,7 @@ Supports STL/OBJ geometry import with ray tracing for inside/outside determinati This reduces numerical dissipation at low Mach numbers while recovering the standard scheme at high Mach numbers. -**Chen correction** (`low_Mach = 2`): anti-dissipation pressure correction (APC) added to the HLLC flux: +**Chen correction** (`low_Mach = 2`, \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] @@ -842,3 +844,4 @@ Additionally, the **artificial Mach number** technique (`pi_fac`) modifies \f$\p \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/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/references.bib b/docs/references.bib new file mode 100644 index 0000000000..abbd1fc4e9 --- /dev/null +++ b/docs/references.bib @@ -0,0 +1,592 @@ +% ========================================================================== +% 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} +} + +@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} +} + +@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} +} + +@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} +} + +@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} +} + +% --- 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} +} + +@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} +} + +@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} +} + +@article{Plesset49, + author = {M. S. Plesset}, + title = {The dynamics of cavitation bubbles}, + journal = {Journal of Applied Mechanics}, + volume = {16}, + pages = {277--282}, + year = {1949} +} + +@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} +} + +@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} +} + +@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} +} + +@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} +} + +@article{Bryngelson23, + author = {S. H. Bryngelson and R. O. Fox and T. Colonius}, + title = {{QBMMlib}: {A} library of quadrature-based moment methods}, + journal = {SoftwareX}, + volume = {21}, + pages = {101283}, + year = {2023} +} + +% --- 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} +} + +% --- 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} +} + +@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} +} + +@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} +} + +% --- 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} +} + +% --- Chemistry --- + +@article{Cisneros25, + author = {E. Cisneros-Garibay and T. Passot and A. Brau and T. L. Jackson}, + 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 = {2025} +} + +% --- 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} +} + +@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} +} + +% --- IGR --- + +@inproceedings{Wilfong25a, + author = {B. Wilfong and A. Jain and S. H. Bryngelson}, + title = {Information geometric regularization of the barotropic {E}uler equations}, + booktitle = {SC'25: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis}, + year = {2025} +} + +@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} +} + +% --- 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} +} + +% --- 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} +} + +@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} +} + +@book{Toro13, + author = {E. F. Toro}, + title = {Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction}, + publisher = {Springer}, + edition = {3rd}, + year = {2013} +} + +@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} +} + +% --- 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} +} + +@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} +} + +@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} +} + +@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} +} + +@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} +} + +@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} +} + +@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} +} + +% --- 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} +} + +@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} +} + +% --- 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} +} + +@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} +} + +@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} +} + +% --- 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} +} + +@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} +} + +% --- 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} +} + +@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} +} + +% --- MFC papers --- + +@article{Bryngelson19, + 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 = {arXiv preprint arXiv:1907.10512}, + year = {2019} +} + +@article{Schmidmayer19, + 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 = {arXiv preprint arXiv:1903.08242}, + year = {2019} +} + +@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} +} + +% --- 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} +} + +@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} +} + +@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} +} + +@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} +} + +@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} +} From 793fa7aba4abdf21871d5c9f6d560d0b6070c5af Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Thu, 12 Feb 2026 17:34:40 -0500 Subject: [PATCH 5/6] Add author-name citations, CI drift checks, and bib corrections Post-process Doxygen HTML to transform bare [N] citation links into "Author et al. [N]" style, integrated into the CMake doc build. Extend lint_docs.py with cite key validation against references.bib and parameter name validation against the REGISTRY. Fix bib metadata (DOIs, published venues, author lists) and align keys with publication years. Remove stale hyper_cleaning recommendation from mhd dependency. Co-Authored-By: Claude Opus 4.6 --- CMakeLists.txt | 3 + docs/documentation/case.md | 12 +- docs/documentation/equations.md | 56 ++++---- docs/postprocess_citations.py | 212 ++++++++++++++++++++++++++++ docs/references.bib | 201 +++++++++++++++++--------- toolchain/mfc/lint_docs.py | 110 ++++++++++++++- toolchain/mfc/params/definitions.py | 6 +- 7 files changed, 491 insertions(+), 109 deletions(-) create mode 100644 docs/postprocess_citations.py diff --git a/CMakeLists.txt b/CMakeLists.txt index f803027c55..f2ad948004 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -753,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" ) @@ -772,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/documentation/case.md b/docs/documentation/case.md index 1e2a5498c7..d0c6431839 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -421,7 +421,7 @@ Details of implementation of viscosity in MFC can be found in \cite Coralic15. | `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 (\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 (\cite Batten97); [2] Pressure-velocity* (\cite Toro13) | +| `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 | @@ -452,7 +452,7 @@ Details of implementation of viscosity in MFC can be found in \cite Coralic15. 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 \cite Bryngelson19. +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. @@ -468,11 +468,11 @@ Tangential velocities require viscosity, `weno_avg = T`, and `bc_[x,y,z]%%end = - `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 (\cite Johnsen08), 5-equation model (\cite Allaire02), and 6-equation model (\cite Saurel09), respectively. -The difference of the two models is assessed by (\cite Schmidmayer19). +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 \cite Schmidmayer19. +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. @@ -512,7 +512,7 @@ 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 (\cite Toro13). +`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 (\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`. @@ -521,7 +521,7 @@ It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$ `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 (\cite Batten97), and indirect method that approximates the pressures and velocity (\cite Toro13), 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. diff --git a/docs/documentation/equations.md b/docs/documentation/equations.md index 7fb2635f20..dff5462e20 100644 --- a/docs/documentation/equations.md +++ b/docs/documentation/equations.md @@ -5,7 +5,7 @@ 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. -For full citations of MFC papers, see @ref papers. Foundational references for each model are cited inline; see the \ref citelist "Bibliography" for full details. +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. --- @@ -31,7 +31,7 @@ The parameter `model_eqns` (1, 2, 3, or 4) selects the governing equation set. ### 2.1 Five-Equation Model (`model_eqns = 2`) -The primary workhorse model (\cite Allaire02). The state vector is: +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] @@ -63,7 +63,7 @@ Setting `alt_soundspeed = .true.` enables the \f$K\f$ correction (\cite Kapila01 ### 2.2 Six-Equation Model (`model_eqns = 3`) -Allows pressure disequilibrium between phases (\cite Saurel09). +Allows pressure disequilibrium between phases (\cite Saurel09; \cite Wilfong26 Sec. 2.1). **Continuity and momentum:** Same as the five-equation model. @@ -101,7 +101,7 @@ See Section 8 (Phase Change) below for details. ## 3. Equations of State -### 3.1 Stiffened Gas EOS (\cite Menikoff89; \cite LeMetayer04) +### 3.1 Stiffened Gas EOS (\cite Menikoff89; \cite LeMetayer04; \cite Wilfong26 Sec. 2.2) The primary closure for each phase: @@ -151,7 +151,7 @@ Temperature is obtained from the internal energy by Newton iteration: \f[e_m(T) = \frac{\hat{h}_m(T) - R_u\,T}{W_m}\f] -**NASA polynomial enthalpies** (\cite McBride93)**:** +**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] @@ -185,7 +185,7 @@ Input parameters: `Re_inv` (shear and volume Reynolds numbers per fluid). --- -## 5. Cylindrical Coordinates (`cyl_coord = .true.`) +## 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: @@ -198,7 +198,7 @@ Additional geometric source terms appear with \f$1/r\f$ factors in the continuit --- -## 6. Sub-Grid Bubble Dynamics +## 6. Sub-Grid Bubble Dynamics (\cite Wilfong26 Sec. 4.1) ### 6.1 Euler-Euler Bubbles (`bubbles_euler = .true.`) @@ -262,7 +262,7 @@ and the local liquid sound speed: \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 Bryngelson23) +#### 6.1.6 QBMM Moment Transport (`qbmm = .true.`) (\cite Bryngelson20) **Population balance equation:** @@ -314,7 +314,7 @@ Each bubble is tracked individually with Keller-Miksis dynamics and 4th-order ad ## 7. Fluid-Structure Interaction -### 7.1 Hypoelastic Model (`hypoelasticity = .true.`) (\cite Rodriguez19) +### 7.1 Hypoelastic Model (`hypoelasticity = .true.`) (\cite Rodriguez19; \cite Wilfong26 Sec. 4.1.6) **Source:** `src/simulation/m_hypoelastic.fpp` @@ -342,7 +342,7 @@ where \f$\mathbf{l} = \nabla \mathbf{u}\f$ is the velocity gradient and \f$\math 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) +### 7.2 Hyperelastic Model (`hyperelasticity = .true.`) (\cite Kamrin12; \cite Wilfong26 Sec. 4.1.6) **Source:** `src/simulation/m_hyperelastic.fpp` @@ -370,7 +370,7 @@ where \f$J = \det(\mathbf{F})\f$. --- -## 8. Phase Change (`relax = .true.`) +## 8. Phase Change (`relax = .true.`) (\cite Wilfong26 Sec. 4.1.3) **Source:** `src/common/m_phase_change.fpp` @@ -408,7 +408,7 @@ Solved via 2D Newton-Raphson. --- -## 9. Chemistry and Combustion (`chemistry = .true.`) +## 9. Chemistry and Combustion (`chemistry = .true.`) (\cite Wilfong26 Sec. 4.1.7) **Source:** `src/common/m_chemistry.fpp` @@ -436,11 +436,11 @@ 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 Cisneros25), which provides symbolic abstractions for thermochemistry that enable portable GPU computation and automatic differentiation of chemical source terms. +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) +## 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` @@ -468,7 +468,7 @@ The capillary stress divergence is added to the momentum and energy equations. T ## 11. Magnetohydrodynamics -### 11.1 Ideal MHD (`mhd = .true.`) +### 11.1 Ideal MHD (`mhd = .true.`) (\cite Wilfong26 Sec. 4.1.9) **Continuity:** @@ -498,9 +498,9 @@ The capillary stress divergence is added to the momentum and energy equations. T \f[v_A = \sqrt{\frac{|\mathbf{B}|^2}{\rho}}\f] -Uses the HLLD Riemann solver (`riemann_solver = 3`). Hyperbolic divergence cleaning (`hyper_cleaning = .true.`) via the GLM method (\cite Dedner02). +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.`) +### 11.2 Relativistic MHD (`relativity = .true.`) (\cite Wilfong26 Sec. 4.1.10) **Conserved variables:** @@ -542,7 +542,7 @@ Uses Lax-Friedrichs flux (replaces WENO + Riemann solver). --- -## 13. Body Forces (`bodyforces = .true.`) +## 13. Body Forces (`bf_x`, `bf_y`, `bf_z`) **Source:** `src/simulation/m_body_forces.fpp` @@ -641,7 +641,7 @@ where \f$a\f$ and \f$b\f$ are the left and right slope differences. \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 = \exp(\text{sign}(s)\,\beta\,(2C - 1)) / \cosh(\beta) - 1) / \tanh(\beta)\f$ and \f$\beta\f$ controls compression steepness. +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 @@ -667,7 +667,11 @@ Four-state solver resolving the contact discontinuity. Star-state satisfies: \f[p_L^* = p_R^* = p^*, \qquad u_L^* = u_R^* = u^*\f] -#### HLLD (`riemann_solver = 3`, MHD only) +#### 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$: @@ -677,10 +681,6 @@ Seven-state solver for ideal MHD resolving fast magnetosonic, Alfven, and contac 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. -#### Exact (`riemann_solver = 4`) - -Iterative exact Riemann solver. - ### 15.3 Time Integration **Source:** `src/simulation/m_time_steppers.fpp` @@ -705,7 +705,7 @@ Iterative exact Riemann solver. \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 (`adapt_dt = .true.`) +#### Adaptive Time Stepping (`adap_dt = .true.`) Embedded RK pairs for error estimation with Hairer-Norsett-Wanner algorithm for initial step size. @@ -783,7 +783,7 @@ Used for viscous fluxes and velocity gradients. **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 = a\left(\frac{\partial \rho}{\partial x} + \frac{2}{\partial x}\frac{\partial p}{\partial x}\right)\f] +\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] @@ -795,7 +795,7 @@ For non-reflecting boundaries, the incoming wave amplitudes are set to zero. **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) +### 16.3 Immersed Boundary Method (`ib = .true.`) (\cite Tseng03; \cite Mittal05; \cite Wilfong26 Sec. 4.1.1) **Source:** `src/simulation/m_ibm.fpp` @@ -815,7 +815,7 @@ Supports STL/OBJ geometry import with ray tracing for inside/outside determinati --- -## 17. Low Mach Number Corrections +## 17. Low Mach Number Corrections (\cite Wilfong26 Sec. 4.2.4) **Thornber correction** (`low_Mach = 1`, \cite Thornber08): modifies the reconstructed velocities at cell interfaces: diff --git a/docs/postprocess_citations.py b/docs/postprocess_citations.py new file mode 100644 index 0000000000..a016bf0247 --- /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]] = {} + + # 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 index abbd1fc4e9..3aabb266bf 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -12,7 +12,8 @@ @article{Allaire02 volume = {181}, number = {2}, pages = {577--616}, - year = {2002} + year = {2002}, + doi = {10.1006/jcph.2002.7143} } @article{Kapila01, @@ -22,7 +23,8 @@ @article{Kapila01 volume = {13}, number = {10}, pages = {3002--3024}, - year = {2001} + year = {2001}, + doi = {10.1063/1.1398042} } @article{Saurel09, @@ -32,7 +34,8 @@ @article{Saurel09 volume = {228}, number = {5}, pages = {1678--1712}, - year = {2009} + year = {2009}, + doi = {10.1016/j.jcp.2008.11.002} } @article{Saurel08, @@ -41,7 +44,8 @@ @article{Saurel08 journal = {Journal of Fluid Mechanics}, volume = {607}, pages = {313--350}, - year = {2008} + year = {2008}, + doi = {10.1017/S0022112008002061} } @article{Zein10, @@ -51,7 +55,8 @@ @article{Zein10 volume = {229}, number = {8}, pages = {2964--2998}, - year = {2010} + year = {2010}, + doi = {10.1016/j.jcp.2009.12.026} } % --- Equations of state --- @@ -63,7 +68,8 @@ @article{Menikoff89 volume = {61}, number = {1}, pages = {75--130}, - year = {1989} + year = {1989}, + doi = {10.1103/RevModPhys.61.75} } @article{LeMetayer04, @@ -73,7 +79,8 @@ @article{LeMetayer04 volume = {43}, number = {3}, pages = {265--276}, - year = {2004} + year = {2004}, + doi = {10.1016/j.ijthermalsci.2003.09.002} } @techreport{McBride93, @@ -93,7 +100,8 @@ @article{Rayleigh17 volume = {34}, number = {200}, pages = {94--98}, - year = {1917} + year = {1917}, + doi = {10.1080/14786440808635681} } @article{Plesset49, @@ -102,7 +110,8 @@ @article{Plesset49 journal = {Journal of Applied Mechanics}, volume = {16}, pages = {277--282}, - year = {1949} + year = {1949}, + doi = {10.1115/1.4009975} } @article{Keller80, @@ -112,7 +121,8 @@ @article{Keller80 volume = {68}, number = {2}, pages = {628--633}, - year = {1980} + year = {1980}, + doi = {10.1121/1.384720} } @techreport{Gilmore52, @@ -130,7 +140,8 @@ @article{Commander89 volume = {85}, number = {2}, pages = {732--746}, - year = {1989} + year = {1989}, + doi = {10.1121/1.397599} } @article{Preston07, @@ -140,7 +151,8 @@ @article{Preston07 volume = {19}, number = {12}, pages = {123302}, - year = {2007} + year = {2007}, + doi = {10.1063/1.2825018} } @phdthesis{Ando10, @@ -157,16 +169,18 @@ @article{Ando11 volume = {37}, number = {6}, pages = {596--608}, - year = {2011} + year = {2011}, + doi = {10.1016/j.ijmultiphaseflow.2011.03.007} } -@article{Bryngelson23, - author = {S. H. Bryngelson and R. O. Fox and T. Colonius}, +@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 = {21}, - pages = {101283}, - year = {2023} + volume = {12}, + pages = {100615}, + year = {2020}, + doi = {10.1016/j.softx.2020.100615} } % --- Euler-Lagrange bubbles --- @@ -177,7 +191,8 @@ @article{Maeda18 journal = {Journal of Computational Physics}, volume = {371}, pages = {994--1017}, - year = {2018} + year = {2018}, + doi = {10.1016/j.jcp.2018.05.029} } % --- Elasticity --- @@ -188,7 +203,8 @@ @article{Rodriguez19 journal = {Journal of Computational Physics}, volume = {379}, pages = {70--90}, - year = {2019} + year = {2019}, + doi = {10.1016/j.jcp.2018.10.035} } @article{Kamrin12, @@ -198,7 +214,8 @@ @article{Kamrin12 volume = {60}, number = {11}, pages = {1952--1969}, - year = {2012} + year = {2012}, + doi = {10.1016/j.jmps.2012.06.003} } @article{Cao19, @@ -207,7 +224,8 @@ @article{Cao19 journal = {International Journal of Solids and Structures}, volume = {169}, pages = {55--71}, - year = {2019} + year = {2019}, + doi = {10.1016/j.ijsolstr.2019.04.002} } % --- Surface tension --- @@ -218,18 +236,20 @@ @article{Schmidmayer17 journal = {Journal of Computational Physics}, volume = {334}, pages = {468--496}, - year = {2017} + year = {2017}, + doi = {10.1016/j.jcp.2017.01.001} } % --- Chemistry --- -@article{Cisneros25, - author = {E. Cisneros-Garibay and T. Passot and A. Brau and T. L. Jackson}, +@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 = {2025} + year = {2026}, + doi = {10.1016/j.cpc.2025.109987} } % --- MHD --- @@ -241,7 +261,8 @@ @article{Miyoshi05 volume = {208}, number = {1}, pages = {315--344}, - year = {2005} + year = {2005}, + doi = {10.1016/j.jcp.2005.02.017} } @article{Dedner02, @@ -251,23 +272,28 @@ @article{Dedner02 volume = {175}, number = {2}, pages = {645--673}, - year = {2002} + year = {2002}, + doi = {10.1006/jcph.2001.6961} } % --- IGR --- @inproceedings{Wilfong25a, - author = {B. Wilfong and A. Jain and S. H. Bryngelson}, - title = {Information geometric regularization of the barotropic {E}uler equations}, + 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}, - year = {2025} + 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} + year = {2024}, + doi = {10.48550/arXiv.2308.14127} } % --- Acoustic sources --- @@ -278,7 +304,8 @@ @article{Maeda17 journal = {Wave Motion}, volume = {75}, pages = {36--49}, - year = {2017} + year = {2017}, + doi = {10.1016/j.wavemoti.2017.08.004} } % --- Riemann solvers --- @@ -290,7 +317,8 @@ @article{Harten83 volume = {25}, number = {1}, pages = {35--61}, - year = {1983} + year = {1983}, + doi = {10.1137/1025002} } @article{Toro94, @@ -299,15 +327,17 @@ @article{Toro94 journal = {Shock Waves}, volume = {4}, pages = {25--34}, - year = {1994} + year = {1994}, + doi = {10.1007/BF01414629} } -@book{Toro13, +@book{Toro09, author = {E. F. Toro}, title = {Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction}, publisher = {Springer}, edition = {3rd}, - year = {2013} + year = {2009}, + doi = {10.1007/b79761} } @article{Batten97, @@ -317,7 +347,8 @@ @article{Batten97 volume = {18}, number = {6}, pages = {1553--1570}, - year = {1997} + year = {1997}, + doi = {10.1137/S1064827593260140} } % --- WENO / spatial reconstruction --- @@ -329,7 +360,8 @@ @article{Jiang96 volume = {126}, number = {2}, pages = {202--228}, - year = {1996} + year = {1996}, + doi = {10.1006/jcph.1996.0130} } @article{Henrick05, @@ -339,7 +371,8 @@ @article{Henrick05 volume = {207}, number = {2}, pages = {542--567}, - year = {2005} + year = {2005}, + doi = {10.1016/j.jcp.2005.01.023} } @article{Borges08, @@ -349,7 +382,8 @@ @article{Borges08 volume = {227}, number = {6}, pages = {3191--3211}, - year = {2008} + year = {2008}, + doi = {10.1016/j.jcp.2007.11.038} } @article{Fu16, @@ -358,7 +392,8 @@ @article{Fu16 journal = {Journal of Computational Physics}, volume = {305}, pages = {333--359}, - year = {2016} + year = {2016}, + doi = {10.1016/j.jcp.2015.10.037} } @article{Balsara00, @@ -368,7 +403,8 @@ @article{Balsara00 volume = {160}, number = {2}, pages = {405--452}, - year = {2000} + year = {2000}, + doi = {10.1006/jcph.2000.6443} } @article{Suresh97, @@ -378,7 +414,8 @@ @article{Suresh97 volume = {136}, number = {1}, pages = {83--99}, - year = {1997} + year = {1997}, + doi = {10.1006/jcph.1997.5745} } @article{Titarev04, @@ -388,7 +425,8 @@ @article{Titarev04 volume = {201}, number = {1}, pages = {238--260}, - year = {2004} + year = {2004}, + doi = {10.1016/j.jcp.2004.05.015} } % --- Time integration --- @@ -400,7 +438,8 @@ @article{Gottlieb98 volume = {67}, number = {221}, pages = {73--85}, - year = {1998} + year = {1998}, + doi = {10.1090/S0025-5718-98-00913-2} } @article{Strang68, @@ -410,7 +449,8 @@ @article{Strang68 volume = {5}, number = {3}, pages = {506--517}, - year = {1968} + year = {1968}, + doi = {10.1137/0705041} } % --- Boundary conditions --- @@ -422,7 +462,8 @@ @article{Thompson87 volume = {68}, number = {1}, pages = {1--24}, - year = {1987} + year = {1987}, + doi = {10.1016/0021-9991(87)90041-6} } @article{Thompson90, @@ -432,7 +473,8 @@ @article{Thompson90 volume = {89}, number = {2}, pages = {439--461}, - year = {1990} + year = {1990}, + doi = {10.1016/0021-9991(90)90152-Q} } @article{Pirozzoli13, @@ -441,7 +483,8 @@ @article{Pirozzoli13 journal = {Journal of Computational Physics}, volume = {248}, pages = {109--126}, - year = {2013} + year = {2013}, + doi = {10.1016/j.jcp.2013.04.021} } % --- Immersed boundary --- @@ -453,7 +496,8 @@ @article{Tseng03 volume = {192}, number = {2}, pages = {593--623}, - year = {2003} + year = {2003}, + doi = {10.1016/j.jcp.2003.07.024} } @article{Mittal05, @@ -462,7 +506,8 @@ @article{Mittal05 journal = {Annual Review of Fluid Mechanics}, volume = {37}, pages = {239--261}, - year = {2005} + year = {2005}, + doi = {10.1146/annurev.fluid.37.061903.175743} } % --- Low Mach corrections --- @@ -474,7 +519,8 @@ @article{Thornber08 volume = {227}, number = {10}, pages = {4873--4894}, - year = {2008} + year = {2008}, + doi = {10.1016/j.jcp.2008.01.036} } @article{Chen22, @@ -483,23 +529,40 @@ @article{Chen22 journal = {Journal of Computational Physics}, volume = {456}, pages = {111027}, - year = {2022} + year = {2022}, + doi = {10.1016/j.jcp.2022.111027} } % --- MFC papers --- -@article{Bryngelson19, +@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 = {arXiv preprint arXiv:1907.10512}, - year = {2019} + journal = {Computer Physics Communications}, + volume = {266}, + pages = {107396}, + year = {2021}, + doi = {10.1016/j.cpc.2020.107396} } -@article{Schmidmayer19, +@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 = {arXiv preprint arXiv:1903.08242}, - year = {2019} + journal = {Journal of Computational Physics}, + volume = {402}, + pages = {109080}, + year = {2020}, + doi = {10.1016/j.jcp.2019.109080} } @article{Coralic14, @@ -508,7 +571,8 @@ @article{Coralic14 journal = {Journal of Computational Physics}, volume = {274}, pages = {95--121}, - year = {2014} + year = {2014}, + doi = {10.1016/j.jcp.2014.06.003} } % --- Theses --- @@ -542,7 +606,8 @@ @article{Tiwari13 journal = {Journal of Computational Physics}, volume = {252}, pages = {290--309}, - year = {2013} + year = {2013}, + doi = {10.1016/j.jcp.2013.06.021} } @article{Michalke64, @@ -552,7 +617,8 @@ @article{Michalke64 volume = {19}, number = {4}, pages = {543--556}, - year = {1964} + year = {1964}, + doi = {10.1017/S0022112064000908} } @article{Tam05, @@ -562,7 +628,8 @@ @article{Tam05 volume = {284}, number = {3--5}, pages = {947--984}, - year = {2005} + year = {2005}, + doi = {10.1016/j.jsv.2004.07.013} } @article{Guo23, @@ -571,7 +638,8 @@ @article{Guo23 journal = {Journal of Fluid Mechanics}, volume = {970}, pages = {A2}, - year = {2023} + year = {2023}, + doi = {10.1017/jfm.2023.548} } @inproceedings{Childs12, @@ -588,5 +656,6 @@ @article{Xu2019 journal = {Physics of Fluids}, volume = {31}, number = {9}, - year = {2019} + year = {2019}, + doi = {10.1063/1.5116374} } diff --git a/toolchain/mfc/lint_docs.py b/toolchain/mfc/lint_docs.py index 4b46950905..6c060a7ddb 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 assignments like `param = value` or `param` +PARAM_RE = re.compile(r"`([a-z_][a-z0-9_%]*?)(?:\s*=\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..09c110225d 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -781,11 +781,7 @@ def get_value_label(param_name: str, value: int) -> str: "requires": ["sigma"], } }, - "mhd": { - "when_true": { - "recommends": ["hyper_cleaning"], - } - }, + "mhd": {}, "relativity": { "when_true": { "requires": ["mhd"], From 04f57ed2799bab7f4d3f053dcd82436fa18f4fd9 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Thu, 12 Feb 2026 17:48:41 -0500 Subject: [PATCH 6/6] Fix low_Mach index swap, int_comp param name, and reviewer feedback MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Swap low_Mach = 1 (Chen pressure correction) and low_Mach = 2 (Thornber velocity correction) to match source code and case.md. Fix interface_compression → int_comp (actual parameter name). Remove empty mhd dependency dict. Fix type annotations and unused variable in postprocess_citations.py. Extend PARAM_RE to catch comparison operators (>, <) in addition to assignment (=). Co-Authored-By: Claude Opus 4.6 --- docs/documentation/equations.md | 18 +++++++++--------- docs/postprocess_citations.py | 4 ++-- toolchain/mfc/lint_docs.py | 4 ++-- toolchain/mfc/params/definitions.py | 1 - 4 files changed, 13 insertions(+), 14 deletions(-) diff --git a/docs/documentation/equations.md b/docs/documentation/equations.md index dff5462e20..59c6a5c776 100644 --- a/docs/documentation/equations.md +++ b/docs/documentation/equations.md @@ -637,7 +637,7 @@ Primitive variable reconstruction is used to avoid spurious oscillations at inte where \f$a\f$ and \f$b\f$ are the left and right slope differences. -**THINC interface compression** (`interface_compression > 0`): applies hyperbolic tangent compression near material interfaces: +**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] @@ -817,21 +817,21 @@ Supports STL/OBJ geometry import with ray tracing for inside/outside determinati ## 17. Low Mach Number Corrections (\cite Wilfong26 Sec. 4.2.4) -**Thornber correction** (`low_Mach = 1`, \cite Thornber08): modifies the reconstructed velocities at cell interfaces: +**Chen correction** (`low_Mach = 1`, \cite Chen22): anti-dissipation pressure correction (APC) added to the HLLC flux: -\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[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[z = \min\!\left(\max\!\left(\frac{|u_L|}{c_L},\;\frac{|u_R|}{c_R}\right),\;1\right)\f] +\f[\text{APC}_{p\mathbf{u}} = (z - 1)\,p_d\,\hat{\mathbf{n}}, \qquad \text{APC}_{pE} = (z - 1)\,p_d\,S_*\f] -This reduces numerical dissipation at low Mach numbers while recovering the standard scheme at high Mach numbers. +The corrected HLLC flux in the star region becomes \f$\mathbf{F}^{*} + \text{APC}\f$. -**Chen correction** (`low_Mach = 2`, \cite Chen22): anti-dissipation pressure correction (APC) added to the HLLC flux: +**Thornber correction** (`low_Mach = 2`, \cite Thornber08): modifies the reconstructed velocities at cell interfaces: -\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[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[\text{APC}_{p\mathbf{u}} = (z - 1)\,p_d\,\hat{\mathbf{n}}, \qquad \text{APC}_{pE} = (z - 1)\,p_d\,S_*\f] +\f[z = \min\!\left(\max\!\left(\frac{|u_L|}{c_L},\;\frac{|u_R|}{c_R}\right),\;1\right)\f] -The corrected HLLC flux in the star region becomes \f$\mathbf{F}^{*} + \text{APC}\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. diff --git a/docs/postprocess_citations.py b/docs/postprocess_citations.py index a016bf0247..a5c19df5c3 100644 --- a/docs/postprocess_citations.py +++ b/docs/postprocess_citations.py @@ -20,7 +20,7 @@ def parse_bib_authors(bib_path: Path) -> dict[str, tuple[str, int, str | None]]: (first_author_surname, author_count, second_author_surname | None) """ text = bib_path.read_text(encoding="utf-8") - entries: dict[str, tuple[str, int]] = {} + entries: dict[str, tuple[str, int, str | None]] = {} # Find all bib entries: @type{Key, for m in re.finditer(r"@\w+\{(\w+),", text): @@ -154,7 +154,7 @@ def check_bare_citations(html_dir: Path, authors: dict) -> list[str]: 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] + surname, _count, second_surname = authors[key] has_prefix = ( before.endswith(surname) or "et al." in before[-20:] or diff --git a/toolchain/mfc/lint_docs.py b/toolchain/mfc/lint_docs.py index 6c060a7ddb..f6318b8f03 100644 --- a/toolchain/mfc/lint_docs.py +++ b/toolchain/mfc/lint_docs.py @@ -23,8 +23,8 @@ # Match \cite Key references in markdown (Doxygen-style) CITE_RE = re.compile(r"\\cite\s+(\w+)") -# Match backtick-quoted parameter assignments like `param = value` or `param` -PARAM_RE = re.compile(r"`([a-z_][a-z0-9_%]*?)(?:\s*=\s*[^`]*)?\s*`") +# 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( diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 09c110225d..c08249c3d8 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -781,7 +781,6 @@ def get_value_label(param_name: str, value: int) -> str: "requires": ["sigma"], } }, - "mhd": {}, "relativity": { "when_true": { "requires": ["mhd"],