All mach projection method

cmake/MFCTargets.cmake

@@ -70,7 +70,12 @@ exit 0
         #                          six-equation relaxation kernel silently drops the
         #                          internal-energy update, leaving alpha and
         #                          alpha_rho correct and int_en zero.
-        foreach(_no_inline_file m_start_up m_cbc m_pressure_relaxation)
+        #   m_multigrid            likewise: inlining its device helpers into the
+        #                          PCG and preconditioner kernels drops their
+        #                          stores, and PCG diverges.
+        #   m_multigrid_exchange   likewise: inlining drops the receive offsets
+        #                          of its message layout, and MPI truncates.
+        foreach(_no_inline_file m_start_up m_cbc m_pressure_relaxation m_multigrid m_multigrid_exchange)
             set_source_files_properties(
                 "${CMAKE_BINARY_DIR}/fypp/${ARGS_TARGET}/${_no_inline_file}.fpp.f90"
                 TARGET_DIRECTORY ${ARGS_TARGET}

docs/documentation/equations.md

@@ -410,7 +410,7 @@ using a single thermal-equilibrium mixture temperature \f$T\f$ (the same value c
 
 The flux is direction-split and face-centered: a two-point difference of \f$T\f$ across each face, exact for this term since it has no cross-derivatives (an axis-cell correction applies in cylindrical coordinates; see the limitation below). The closure is implemented for the stiffened-gas and ideal-gas equations of state only, and for `model_eqns = 2` (5-equation) or `model_eqns = 3` (6-equation) — both carry the volume fractions \f$\alpha_i\f$ that weight \f$k\f$, which `model_eqns = 1` (gamma law) does not. Heat conduction is independent of `viscous`: it can be active in an otherwise inviscid simulation. It is not supported together with `igr` or `chemistry`; see @ref sec-fluid-materials in the case documentation for the full set of input constraints.
 
-A thermal diffusion CFL limit (`TCFL`) is added to the adaptive time-step candidates alongside `ICFL`/`VCFL`/`CCFL`.
+A thermal diffusion CFL limit (`TCFL`), \f$k\,\Delta t\sum_d \Delta x_d^{-2}/(\rho c_v) \le\f$ `cfl_target`, is added to the adaptive time-step candidates alongside `ICFL`/`VCFL`/`CCFL`.
 
 **Known limitation:** in cylindrical coordinates, the cell adjacent to the axis carries a non-converging \f$\sim\f$1-3% error in the conduction term. It is inherited from the two-point face-gradient every MFC diffusive flux uses, applied across the coordinate singularity at the axis — the same pattern as the two-point gradient in the chemistry diffusion flux (`src/common/m_chemistry.fpp`).
 
@@ -1016,6 +1016,38 @@ Used for viscous fluxes and velocity gradients.
 
 \f[f'_i\big|_\text{right} = \frac{3f_i - 4f_{i-1} + f_{i-2}}{x_i - x_{i-2}}\f]
 
+### 15.6 All-Mach Pressure Projection (`proj_method = .true.`)
+
+**Source:** `src/simulation/m_projection.fpp`, `src/simulation/m_projection_helpers.fpp` (projection); `src/simulation/m_multigrid.fpp`, `src/simulation/m_multigrid_levels.fpp`, `src/simulation/m_multigrid_exchange.fpp`, `src/simulation/m_multigrid_bottom.fpp` (multigrid pressure solve)
+
+Following \cite Fuster18, the five-equation model is advanced with the pressure implicit, so the time step is limited by the flow speed rather than the sound speed. Each Runge-Kutta stage, of effective step \f$\tau\f$, has four parts.
+
+**Advection.** Partial densities, volume fractions and momentum are transported by a face velocity \f$u_f\f$ kept between stages and upwinded on its sign; the momentum flux is the summed mass flux times the upwind velocity. The pressure is advected, \f$\partial_t p = -\mathbf{u}_f\cdot\nabla p\f$, giving \f$p_\text{adv}\f$. Energy is not transported.
+
+**Face predictor.** From the advected (star) state,
+
+\f[u^*_f = \frac{1}{2}\left[\left(\frac{\rho u}{\rho}\right)_L + \left(\frac{\rho u}{\rho}\right)_R\right] + \tau\,a_f\f]
+
+where \f$a_f\f$ is the body-force and, with `surface_tension_model = 2`, capillary acceleration on the face.
+
+**Pressure equation.**
+
+\f[p - \rho c^2\tau^2\,\nabla\cdot\left(\frac{1}{\rho_f}\nabla p\right) = p_\text{adv} - \rho c^2\tau\,\nabla\cdot\mathbf{u}^*_f\f]
+
+discretized as a symmetric positive-definite seven-point system with face conductances \f$A_f/(\rho_f d_f)\f$, \f$\rho_f\f$ the arithmetic face density. It is solved by flexible preconditioned CG, with one V- or K-cycle (\cite Notay08) of geometric multigrid per iteration; see @ref sec-projection-solve-tuning.
+
+**Correction.** The faces take the solved pressure gradient,
+
+\f[u_f = u^*_f - \frac{\tau}{\rho_f}\frac{p_R - p_L}{d_f}\f]
+
+so \f$\nabla\cdot u_f\f$ matches the solved pressure exactly. Each cell's momentum gains \f$\rho\f$ times the mean of its faces' net acceleration (body force and capillary force less the pressure gradient, zero on walls), so a hydrostatic balance on the faces leaves the cells at rest. The energy is rebuilt from the solved pressure, \f$\rho E = \Gamma p + \Pi_\infty + \rho q_v + \frac{1}{2}\rho|\mathbf{u}|^2\f$. With `proj_lag_corr` the transported variables are then moved by \f$u_f - u^\text{old}_f\f$, the change of the face velocity they were transported with.
+
+Viscous stresses, capillary stresses and heat conduction are explicit; conduction enters the pressure equation as \f$\partial_t p \mathrel{+}= -\nabla\cdot\mathbf{q}/\Gamma\f$.
+
+**Time step.** `ICFL` becomes the advective \f$\Delta t\max_d |u_d|/\Delta x_d\f$ (reported as `AdvCFL`), and `proj_max_acfl` optionally caps the acoustic CFL (`AcCFL`); the viscous, capillary and thermal limits are unchanged.
+
+**Limitations:** `model_eqns = 2`, stiffened- or ideal-gas fluids, Cartesian coordinates. A mixed interface cell compresses with the advected mixture's modulus rather than Wood's, so gas bubbles barely expand or contract with pressure.
+
 ---
 
 ## 16. Boundary Conditions

docs/module_categories.json

@@ -14,7 +14,13 @@
             "m_riemann_solver_hypo_hlld",
             "m_muscl",
             "m_variables_conversion",
-            "m_thinc"
+            "m_thinc",
+            "m_projection",
+            "m_projection_helpers",
+            "m_multigrid",
+            "m_multigrid_levels",
+            "m_multigrid_exchange",
+            "m_multigrid_bottom"
         ]
     },
     {

docs/postprocess_citations.py

@@ -107,7 +107,7 @@ def process_html(html_path: Path, authors: dict) -> bool:
         return False
 
     def replace_cite(match: re.Match) -> str:
-        key = match.group(1)
+        key = match.group(1).lower()
         full_link = match.group(0)
 
         info = authors.get(key)
@@ -150,7 +150,7 @@ def check_bare_citations(html_dir: Path, authors: dict) -> list[str]:
             continue
         text = html_file.read_text(encoding="utf-8")
         for m in CITE_RE.finditer(text):
-            key = m.group(1)
+            key = m.group(1).lower()
             if key not in authors:
                 bare.append(f"  {html_file.name}: [{m.group(2)}] (CITEREF_{key}) — no bib entry")
                 continue

docs/references.bib

@@ -722,3 +722,24 @@ @phdthesis{Spratt24
   year   = {2024},
   doi    = {10.7907/g34e-6p65}
 }
+
+@article{Fuster18,
+  author  = {D. Fuster and S. Popinet},
+  title   = {An all-{M}ach method for the simulation of bubble dynamics problems in the presence of surface tension},
+  journal = {Journal of Computational Physics},
+  volume  = {374},
+  pages   = {752--768},
+  year    = {2018},
+  doi     = {10.1016/j.jcp.2018.07.055}
+}
+
+@article{Notay08,
+  author  = {Y. Notay and P. S. Vassilevski},
+  title   = {Recursive {K}rylov-based multigrid cycles},
+  journal = {Numerical Linear Algebra with Applications},
+  volume  = {15},
+  number  = {5},
+  pages   = {473--487},
+  year    = {2008},
+  doi     = {10.1002/nla.542}
+}

examples/2D_couette_projection/README.md

@@ -0,0 +1,7 @@
+# Plane Couette Flow Startup (2D, all-Mach pressure projection)
+
+Fluid at rest between a fixed lower wall and an upper wall that starts moving at U, periodic along the gap. The fluid is a
+stiffened gas at Mach 0.01. `result.png` shows the validation result: the velocity profile across the gap at several times (symbols)
+against the exact transient (lines).
+
+<img src='result.png' height='MAX_HEIGHT'/>

examples/2D_couette_projection/case.py

@@ -0,0 +1,85 @@
+#!/usr/bin/env python3
+"""
+Startup of plane Couette flow at Mach 0.01 with the all-Mach pressure projection.
+
+Fluid at rest between a fixed lower wall and an upper wall that starts moving at U (bc_y%ve1), periodic in x. Viscosity carries
+the wall's momentum across the gap toward the linear profile u = U y/H (README.md compares the exact transient). rho = U = H =
+1 and Re = U H/nu = 100. A stiffened gas with c = U/Mach makes the flow nearly incompressible, so the projection steps at the
+flow's and viscosity's pace while --explicit (HLLC, low_Mach = 2) must resolve the sound speed, 100 times faster.
+"""
+
+import argparse
+import json
+import math
+import sys
+
+parser = argparse.ArgumentParser(description="Plane Couette flow startup, all-Mach pressure projection")
+parser.add_argument("--explicit", action="store_true", help="explicit HLLC (low_Mach = 2) at the acoustic limit instead of the projection")
+args, _ = parser.parse_known_args()
+
+U, H, rho, Re, mach, gamma, p0 = 1.0, 1.0, 1.0, 100.0, 0.01, 1.4, 1.0
+ny, nx, tstop, saves, cfl = 64, 32, 30.0, 30, 0.5  # cells across and along the gap, end time (0.3 H^2/nu), outputs, CFL
+nu = U * H / Re
+c = U / mach
+pi_inf = rho * c**2 / gamma - p0  # stiffened gas with sound speed c at p0
+
+# Fixed step: the advective and viscous (4/3 nu dt sum 1/dx^2) limits, or the acoustic one with --explicit
+dy = H / ny
+dt = cfl * min(dy / U, dy**2 / (2 * 4 / 3 * nu)) if not args.explicit else cfl * dy / (U + c)
+nt = int(math.ceil(tstop / dt))
+dt = tstop / nt
+print(f"dt = {dt:.3e}, {nt} steps", file=sys.stderr)
+
+print(
+    json.dumps(
+        {
+            "run_time_info": "T",
+            "x_domain%beg": 0.0,
+            "x_domain%end": nx * dy,
+            "y_domain%beg": 0.0,
+            "y_domain%end": H,
+            "m": nx - 1,
+            "n": ny - 1,
+            "p": 0,
+            "dt": dt,
+            "t_step_start": 0,
+            "t_step_stop": nt,
+            "t_step_save": nt // saves,
+            "num_patches": 1,
+            "model_eqns": 2,
+            "num_fluids": 1,
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "mapped_weno": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "bc_y%beg": -16,
+            "bc_y%end": -16,
+            "bc_y%ve1": U,
+            "viscous": "T",
+            **({"low_Mach": 2} if args.explicit else {"proj_method": "T"}),
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "T",
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.5 * nx * dy,
+            "patch_icpp(1)%y_centroid": 0.5 * H,
+            "patch_icpp(1)%length_x": nx * dy,
+            "patch_icpp(1)%length_y": H,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": p0,
+            "patch_icpp(1)%alpha_rho(1)": rho,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            "fluid_pp(1)%eos": "stiffened_gas",
+            "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
+            "fluid_pp(1)%pi_inf": gamma * pi_inf / (gamma - 1.0),
+            "fluid_pp(1)%Re(1)": 1.0 / (rho * nu),
+        }
+    )
+)

examples/2D_lid_driven_cavity_projection/README.md

@@ -0,0 +1,19 @@
+# Lid-Driven Cavity (2D, all-Mach pressure projection)
+
+A unit square of a single ideal gas with no-slip walls; the top wall moves at 0.5. With ρ = 1 and μ = 1e-4 the Reynolds number is
+5000, and the lid's Mach number, against the gas's sound speed of 2.6, is about 0.2. The case runs to t = 200 on 200 × 200 cells;
+`--explicit` runs the HLLC solver (`low_Mach = 2`) instead.
+
+```shell
+./mfc.sh run examples/2D_lid_driven_cavity_projection/case.py --case-optimization
+```
+
+`result.png` shows the validation result: the u velocity along the vertical centerline and the v velocity along the horizontal
+centerline at t = 200, against the Re = 5000 data of Ghia, Ghia & Shin (J. Computational Phys. 48:387-411, 1982), normalized by the
+lid speed.
+
+<img src='result.png' height='MAX_HEIGHT'/>
+
+Both centerline profiles agree with the reference to within 0.04 of the lid speed (on Ghia et al.'s 257 × 257 grid against
+200 × 200 here). At t = 200 the flow is close to steady: its velocities change by under 1% of the lid speed over the last 10
+time units.

examples/2D_lid_driven_cavity_projection/case.py

@@ -0,0 +1,82 @@
+#!/usr/bin/env python3
+"""
+2D lid-driven cavity at Re = 5000 (lid speed 0.5, rho = 1, Re(1) = 1/mu = 1e4) with the all-Mach pressure projection: a single
+ideal gas, no-slip walls, the top one moving. --explicit runs the explicit HLLC solver (low_Mach = 2) at the acoustic limit
+instead.
+"""
+
+import argparse
+import json
+
+parser = argparse.ArgumentParser(description="2D lid-driven cavity, all-Mach pressure projection")
+parser.add_argument("--explicit", action="store_true", help="explicit HLLC (low_Mach = 2) at the acoustic limit instead of the projection")
+args, _ = parser.parse_known_args()
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0,
+            "x_domain%end": 1.0,
+            "y_domain%beg": 0.0,
+            "y_domain%end": 1.0,
+            "m": 199,
+            "n": 199,
+            "p": 0,
+            "cfl_adap_dt": "T",
+            "cfl_target": 0.8,
+            "n_start": 0,
+            "t_stop": 200.0,
+            "t_save": 1.0,
+            # Simulation Algorithm Parameters
+            "num_patches": 1,
+            "model_eqns": "5eq",
+            "alt_soundspeed": "F",
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            # "mixture_err": "T",
+            "time_stepper": "rk3",
+            "weno_order": 5,
+            "weno_eps": 1e-16,
+            "mapped_weno": "T",
+            # "weno_Re_flux": "T",
+            "mp_weno": "T",
+            "weno_avg": "T",
+            "riemann_solver": "hllc",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
+            "bc_x%beg": -16,
+            "bc_x%end": -16,
+            "bc_y%beg": -16,
+            "bc_y%end": -16,
+            "bc_y%ve1": 0.5,
+            "viscous": "T",
+            **({"low_Mach": 2} if args.explicit else {"proj_method": "T", "proj_max_acfl": 200}),
+            # Formatted Database Files Structure Parameters
+            "format": "silo",
+            "precision": "double",
+            "prim_vars_wrt": "T",
+            "omega_wrt(3)": "T",
+            "fd_order": 4,
+            "parallel_io": "T",
+            # Patch 1: Base
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%y_centroid": 0.5,
+            "patch_icpp(1)%length_x": 1.0,
+            "patch_icpp(1)%length_y": 1.0,
+            "patch_icpp(1)%vel(1)": 0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": 5,
+            "patch_icpp(1)%alpha_rho(1)": 1.0,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0 / (1.4 - 1.0),
+            "fluid_pp(1)%eos": "ideal_gas",
+            "fluid_pp(1)%Re(1)": 1e4,
+        }
+    )
+)

examples/2D_oscillating_bubble_projection/README.md

@@ -0,0 +1,26 @@
+# Oscillating Bubble (2D, all-Mach pressure projection)
+
+A validation of surface tension. A gas bubble of radius R in a liquid (density ratio 1000) starts at rest as an ellipse of
+semi-axes R(1 ± 0.05), the n = 2 shape mode, and oscillates about its circle. For small amplitude ε the axes follow
+R ± εR cos(ωt) e^(−βt), with Lamb's frequency of a 2D (cylindrical) interface,
+
+ω² = n(n² − 1) σ / ((ρ_g + ρ_l) R³),
+
+and the viscous decay of the irrotational mode in the liquid outside the bubble, β = 2n(n + 1) ν_l / R² (the dissipation of the
+potential flow r⁻ⁿ cos nθ over twice its kinetic energy). One quadrant is simulated, 32 cells per radius, with symmetry planes
+along both axes. The liquid is a stiffened gas with sound speed 100 and the bubble an ideal gas (sound speed ~370) against flow
+speeds of ~0.1, so the projection steps at the capillary limit; `--explicit` runs the HLLC solver at the acoustic one. The
+reference pressure, 100, keeps the bubble's breathing mode near ω ≈ 14, well above the shape mode, so the two do not exchange
+energy.
+
+```shell
+./mfc.sh run examples/2D_oscillating_bubble_projection/case.py -n 4
+```
+
+`result.png` shows the validation result: the deviation of the bubble's x and y semi-axes from the equivalent radius over five
+periods (symbols) against the linear theory (line).
+
+<img src='result.png' height='MAX_HEIGHT'/>
+
+The measured frequency, from the zero crossings of the axis difference, is 2.454 against Lamb's 2.448 (+0.23%), and the decay
+rate fitted to the peaks is 0.125 against 0.120. The bubble's area changes by less than 1e-4 over the run.