Split-explicit acoustic substepping for low-Mach time integration (with two-tier shock/interface capture)

docs/documentation/case.md

@@ -488,6 +488,9 @@ See @ref equations "Equations" for the mathematical models these parameters cont
 | `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 (\cite Chen22); [2] Velocity (\cite Thornber08)	 |
+| `acoustic_substepping`     | Logical | Enable split-explicit acoustic substepping for low-Mach flows |
+| `n_acoustic_substeps`      | Integer | Number of acoustic substeps per flow step (0 = auto) |
+| `acoustic_div_damp`        | Real    | Dimensionless grad-div divergence damping coefficient for acoustic substepping (default 0.1; stable for <~ 0.5/num_dims) |
 | `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 Toro09)	 |
 | `weno_Re_flux`             | Logical | Compute velocity gradient using scalar divergence theorem	 |

docs/module_categories.json

@@ -4,6 +4,7 @@
         "modules": [
             "m_rhs",
             "m_time_steppers",
+            "m_acoustic_substep",
             "m_weno",
             "m_riemann_solvers",
             "m_riemann_state",

examples/1D_acoustic_lowmach/case.py

@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+"""
+1D weakly-compressible acoustic pulse (M ≈ 0.01) — baseline harness for
+acoustic_substepping (split-explicit low-Mach time integrator).
+
+Physics:
+  - Single-fluid stiffened-gas air (gamma=1.4, pi_inf=0) on a periodic domain.
+  - Background state: p0=101325 Pa, rho0=1.204 kg/m^3.
+  - Sound speed: c0 = sqrt(gamma*p0/rho0) ≈ 343 m/s.
+  - Background velocity: u0 = M*c0 with M=0.01 ≈ 3.43 m/s (low-Mach).
+  - Initial condition: small Gaussian pressure perturbation dp = 1e-3 * p0.
+
+Set acoustic_substepping = 'T' (and build the numerics in later tasks) to
+exercise the split-explicit pathway.  Default here is 'F' (baseline RK3).
+"""
+
+import json
+import math
+
+# Physical parameters
+gamma = 1.4
+p0 = 101325.0  # background pressure [Pa]
+rho0 = 1.204  # background density [kg/m^3]
+c0 = math.sqrt(gamma * p0 / rho0)  # ~343 m/s
+Mach = 0.01
+u0 = Mach * c0  # background velocity ~3.43 m/s
+
+# Domain
+L = 1.0  # domain length [m]
+Nx = 199  # cells (m parameter)
+dx = L / (Nx + 1)
+
+# CFL-based time stepping. In split-explicit (acoustic_substepping) mode the
+# timestep is set from the ADVECTIVE CFL (uses |u|, not |u|+c), so it is ~1/Mach
+# larger than the acoustic CFL; the acoustic waves are resolved by the internal
+# substeps. CFL-based dt is required for split mode (it is how n_substeps is
+# computed) and also works for the baseline (acoustic_substepping='F') path.
+CFL = 0.5
+dt = CFL * dx / c0  # initial dt guess; the solver recomputes from cfl_target each step
+t_end = L / c0  # one acoustic transit
+
+# Gaussian pressure pulse: dp = eps_p * p0 * exp(-((x-0.5)/sigma)^2)
+# sigma = 0.05 * L (narrow pulse, well-resolved on the grid)
+# Written as an analytic IC expression evaluated at cell centers.
+eps_p = 1.0e-3
+sigma = 0.05 * L
+
+# Stiffened-gas EOS: Gamma = 1/(gamma-1), pi_inf = 0 for ideal gas
+Gamma = 1.0 / (gamma - 1.0)
+pi_inf = 0.0
+
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0,
+            "x_domain%end": L,
+            "m": Nx,
+            "n": 0,
+            "p": 0,
+            "dt": dt,
+            "cfl_adap_dt": "T",
+            "cfl_target": CFL,
+            "t_step_start": 0,
+            "n_start": 0,
+            "t_stop": t_end,
+            "t_save": t_end / 10.0,
+            # Simulation Algorithm Parameters
+            "num_patches": 1,
+            "model_eqns": 2,
+            "num_fluids": 1,
+            "alt_soundspeed": "F",
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": "rk3",
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "mp_weno": "F",
+            "weno_avg": "F",
+            "mapped_weno": "F",
+            "null_weights": "F",
+            "riemann_solver": "hllc",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            # Acoustic substepping (split-explicit low-Mach). Set to 'T' to exercise
+            # the split-explicit pathway; CFL-based dt above supports both modes.
+            "acoustic_substepping": "T",
+            # Formatted Database Files Structure Parameters
+            "format": "binary",
+            "precision": "double",
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            # Patch 1: uniform background + small Gaussian pressure pulse
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.5 * L,
+            "patch_icpp(1)%length_x": L,
+            "patch_icpp(1)%vel(1)": u0,
+            "patch_icpp(1)%pres": f"{p0} + {eps_p * p0} * exp(-((x - {0.5 * L}) / {sigma})**2)",
+            "patch_icpp(1)%alpha_rho(1)": rho0,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            # Fluid physical parameters (stiffened-gas air, ideal limit)
+            "fluid_pp(1)%gamma": Gamma,
+            "fluid_pp(1)%pi_inf": pi_inf,
+        }
+    )
+)

examples/1D_bubble_acoustic_lowmach/case.py

@@ -0,0 +1,132 @@
+#!/usr/bin/env python3
+"""
+1D Euler-Euler bubble, acoustically driven, at low Mach -- acoustic-substep validation.
+
+Physics:
+  - Single carrier fluid (model_eqns=2, num_fluids=1) seeded with a dilute population of
+    Euler-Euler bubbles (bubbles_euler='T', adv_n='T', polytropic Rayleigh-Plesset).
+  - Uniform bubbly background (void fraction vf0) on a periodic [0,1] domain with a small
+    Gaussian pressure pulse centered at x=0.5. The pulse drives a gentle bubble response
+    (compression/expansion) -- a MODERATE forcing, the regime the co-subcycle targets.
+  - Low Mach: u0 = M*c with M=0.1 (c=sqrt(gamma)~1.18, normalised p0=rho0=1), so the split
+    integrator takes the advective-CFL outer step and resolves acoustics + co-subcycles the
+    bubble dynamics on the micro-steps.
+
+Primary validation (R(t)/alpha(t) match):
+  Run split (MFC_SPLIT=T -> acoustic_substepping='T') vs standard (MFC_SPLIT=F, the full
+  RK3 + adaptive bubble Strang split) at the SAME resolution and physical t_stop, and
+  compare the acoustically driven bubble response at the pulse-center probe. The
+  co-subcycled equivalent radius R(t) tracks the standard solver in sign, magnitude, and
+  time-trend; the void fraction alpha(t) under-responds because the bubble number density
+  rides the slow advective transport (see docs/documentation/acousticSubstepping.md,
+  "Euler-Euler bubbles"). Conservation (carrier mass + total energy) is exact and the void
+  fraction stays bounded and positive.
+
+  Diagnostics (read from <case>/D/ 1D ASCII): void fraction (prim slot eqn_idx%alf) and the
+  bubble radius moment at the center cell over the save sequence.
+
+Toggle baseline vs split via the MFC_SPLIT env var:
+  MFC_SPLIT=F  -> acoustic_substepping='F' (standard RK3 + adaptive bubble integrator) [default]
+  MFC_SPLIT=T  -> acoustic_substepping='T' (split-explicit, bubble dynamics co-subcycled)
+
+NOTE on scope: this is MODERATE forcing (0.5% pulse). A violent single-bubble collapse
+(strong forcing) drives the adaptive bubble sub-integrator to a stiffness convergence
+failure regardless of the acoustic substep count -- documented out of scope (see
+docs/documentation/acousticSubstepping.md, "Euler-Euler bubbles").
+"""
+
+import json
+import math
+import os
+
+split = os.environ.get("MFC_SPLIT", "F").upper()
+acoustic_substepping = "T" if split == "T" else "F"
+
+gamma = 1.4
+c0 = math.sqrt(gamma)  # carrier sound speed ~1.183 (normalised p0=rho0=1)
+Mach = 0.1
+u0 = Mach * c0  # uniform mean flow
+
+L = 1.0
+Nx = 99  # 100 cells, dx = 1/100
+dx = L / (Nx + 1)
+
+CFL = 0.5
+dt_init = CFL * dx / c0  # acoustic-CFL initial guess (cfl_adap_dt refines it)
+
+vf0 = 1.0e-3  # background void fraction (dilute)
+pulse = 0.005  # 0.5% Gaussian pressure pulse -> MODERATE bubble forcing
+
+print(
+    json.dumps(
+        {
+            "run_time_info": "T",
+            "m": Nx,
+            "n": 0,
+            "p": 0,
+            "x_domain%beg": 0.0,
+            "x_domain%end": 1.0,
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "dt": dt_init,
+            "cfl_adap_dt": "T",
+            "cfl_target": CFL,
+            "n_start": 0,
+            "t_stop": 0.15,
+            "t_save": 0.01,
+            "num_patches": 1,
+            "model_eqns": 2,
+            "num_fluids": 1,
+            "alt_soundspeed": "F",
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "mapped_weno": "F",
+            "null_weights": "F",
+            "mp_weno": "F",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "acoustic_substepping": acoustic_substepping,
+            "n_acoustic_substeps": 20,
+            "acoustic_div_damp": 0.4,
+            "format": "silo",
+            "precision": "double",
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            "fd_order": 1,
+            "probe_wrt": "T",
+            "num_probes": 1,
+            "probe(1)%x": 0.5,
+            # Patch 1: uniform bubbly background + Gaussian pressure pulse
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%length_x": 1.0,
+            "patch_icpp(1)%vel(1)": u0,
+            "patch_icpp(1)%pres": f"1.0 + {pulse} * exp(-((x - 0.5) / 0.15)**2)",
+            "patch_icpp(1)%alpha_rho(1)": (1.0 - vf0) * 1.0,
+            "patch_icpp(1)%alpha(1)": vf0,
+            "patch_icpp(1)%r0": 1.0,
+            "patch_icpp(1)%v0": 0.0,
+            "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
+            "fluid_pp(1)%pi_inf": 0.0,
+            # Euler-Euler bubbles (co-subcycled with the acoustic substep when split)
+            "bubbles_euler": "T",
+            "bubble_model": 2,
+            "polytropic": "T",
+            "adv_n": "T",
+            "polydisperse": "F",
+            "thermal": 3,
+            "nb": 1,
+            "bub_pp%R0ref": 1.0,
+            "bub_pp%p0ref": 1.0,
+            "bub_pp%rho0ref": 1.0,
+            "bub_pp%ss": 0.0,
+            "bub_pp%pv": 0.0,
+            "bub_pp%mu_l": 0.1,
+            "bub_pp%gam_g": 1.4,
+        }
+    )
+)

examples/1D_material_interface_lowmach/case.py

@@ -0,0 +1,128 @@
+#!/usr/bin/env python3
+"""
+1D material-interface advection at low Mach — two-tier acoustic-substep validation
+(split-explicit low-Mach integrator).
+
+Physics:
+  - Two ideal gases (model_eqns=2, num_fluids=2) on a periodic [0,1] domain.
+  - A heavy-fluid slab (fluid 2, rho=3) in a light ambient (fluid 1, rho=1) forms a
+    crisp material interface (sharp alpha / density contrast). The interface is RESOLVED
+    over a few cells (tanh, half-width w ~ 4 dx): a centered transport scheme advects a
+    resolved profile without Gibbs ringing, whereas a 1-cell-sharp contact would ring
+    (the substep's convective mass/energy transport is centered, not upwinded — see the
+    task report; the robust tier protects the acoustic pressure flux, not mass transport).
+  - Pressure and velocity are UNIFORM (p0, u0), so the exact solution is the slab advecting
+    at u0 with p kept flat; any pressure deviation at the interface is numerical ringing.
+  - Low Mach: u0 = M*c with M=0.05 (c ~ O(1)), so the split integrator takes the
+    advective-CFL outer step and resolves acoustics by the substep loop.
+
+Two-tier check: the alpha-jump criterion flags the interface band, so the robust
+WENO+HLLC-acoustic tier acts there while the uniform far field uses the cheap centered tier.
+
+Diagnostics (read from <case>/D/ ASCII files, 1D):
+  1. Non-oscillatory: pressure (prim slot E) stays ~p0 across the interface (no ringing).
+  2. Per-species mass conservation: sum_x of each partial density (cons slots 1,2) is
+     invariant in time to ~machine precision (periodic BCs, pure advection).
+
+Toggle baseline (full HLLC) vs split via the MFC_SPLIT env var:
+  MFC_SPLIT=F  -> acoustic_substepping='F' (standard RK3 + full HLLC)  [default]
+  MFC_SPLIT=T  -> acoustic_substepping='T' (split-explicit two-tier)
+Run both at the SAME resolution for the baseline-vs-split comparison.
+"""
+
+import json
+import math
+import os
+
+split = os.environ.get("MFC_SPLIT", "F").upper()
+acoustic_substepping = "T" if split == "T" else "F"
+
+# Two distinct ideal gases (exercises the mixture gamma/pi_inf sums)
+gamma1 = 1.4
+gamma2 = 2.4
+eps = 1.0e-6  # minority volume fraction (keeps mixture EOS well-defined)
+
+p0 = 1.0
+rho1 = 1.0  # light ambient fluid
+rho2 = 3.0  # heavy slab fluid (material density contrast)
+
+c1 = math.sqrt(gamma1 * p0 / rho1)  # light-fluid sound speed ~1.183
+Mach = 0.05
+u0 = Mach * c1  # uniform advection velocity
+
+L = 1.0
+Nx = 199  # 200 cells, dx = 1/200
+dx = L / (Nx + 1)
+w = 4.0 * dx  # interface half-width (resolved over a few cells)
+
+CFL = 0.5
+dt_init = CFL * dx / c1
+
+# Advect the slab ~0.25 of the domain so the interface is well resolved; short run.
+T_stop = 0.25 / u0
+t_save = T_stop / 5.0
+
+# Resolved slab indicator phi(x) ~ 1 inside [0.4,0.6], ~0 outside (tanh shoulders).
+phi = f"(0.5*(tanh((x - 0.4)/{w}) - tanh((x - 0.6)/{w})))"
+a2 = f"({eps} + {1.0 - 2.0 * eps}*{phi})"  # fluid-2 volume fraction
+a1 = f"(1.0 - {a2})"
+
+print(
+    json.dumps(
+        {
+            "run_time_info": "T",
+            "x_domain%beg": 0.0,
+            "x_domain%end": L,
+            "m": Nx,
+            "n": 0,
+            "p": 0,
+            "dt": dt_init,
+            "cfl_adap_dt": "T",
+            "cfl_target": CFL,
+            "t_step_start": 0,
+            "n_start": 0,
+            "t_stop": T_stop,
+            "t_save": t_save,
+            "num_patches": 1,
+            "model_eqns": 2,
+            "num_fluids": 2,
+            "alt_soundspeed": "F",
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": "rk3",
+            "weno_order": 5,
+            "weno_eps": 1.0e-6,
+            "mp_weno": "F",
+            "weno_avg": "F",
+            "mapped_weno": "F",
+            "null_weights": "F",
+            "riemann_solver": "hllc",
+            "wave_speeds": "direct",
+            "avg_state": "arithmetic",
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "acoustic_substepping": acoustic_substepping,
+            "n_acoustic_substeps": 0,
+            "acoustic_div_damp": 0.1,
+            "format": "binary",
+            "precision": "double",
+            "prim_vars_wrt": "T",
+            "cons_vars_wrt": "T",
+            "parallel_io": "F",
+            # Patch 1: resolved heavy slab in light ambient (analytic tanh profiles)
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.5 * L,
+            "patch_icpp(1)%length_x": L,
+            "patch_icpp(1)%vel(1)": u0,
+            "patch_icpp(1)%pres": p0,
+            "patch_icpp(1)%alpha_rho(1)": f"{rho1}*{a1}",
+            "patch_icpp(1)%alpha_rho(2)": f"{rho2}*{a2}",
+            "patch_icpp(1)%alpha(1)": a1,
+            "patch_icpp(1)%alpha(2)": a2,
+            "fluid_pp(1)%gamma": 1.0 / (gamma1 - 1.0),
+            "fluid_pp(1)%pi_inf": 0.0,
+            "fluid_pp(2)%gamma": 1.0 / (gamma2 - 1.0),
+            "fluid_pp(2)%pi_inf": 0.0,
+        }
+    )
+)