Combustion: condensed-phase reactive burn (5-/6-equation, Arrhenius kinetics) + IBM burning surfaces

.claude/rules/common-pitfalls.md

@@ -59,6 +59,23 @@ covered in `docs/documentation/contributing.md`.
   Gotcha: ADDING a new file under `toolchain/mfc/params/` needs one reconfigure
   (the custom command's DEPENDS list is globbed at configure time). Under `--case-optimization` the baked-in constants are dropped from the
   namelist, so changing one needs a *rebuild*, not a case edit.
+- Derived-type params (`chem_params`, `lag_params`, `rburn`) are NOT auto-broadcast:
+  `generated_bcast.fpp` covers namelist *scalars* only. Each type needs a hand-written
+  `_emit_<name>` in `toolchain/mfc/params/generators/fortran_gen.py` plus its call site in
+  the `target == "sim"` block, and — if it is read on device — an explicit
+  `$:GPU_UPDATE(device='[name]')` in BOTH `m_global_parameters.fpp` and `m_start_up.fpp`
+  (`GPU_DECLARE` alone does not make it device-resident). Regrouping scalars into a derived
+  type silently drops their broadcast, so every non-root rank keeps the `dflt_real`
+  sentinel; single-rank goldens cannot see this, so pair such a change with a `ppn=2` test
+  and confirm it fails without the emitter.
+- A `patch_ib` member that any `m_ibm` ghost-point code reads must ALSO be set in
+  `s_add_cloud_particle` (`src/simulation/m_particle_cloud.fpp`): `particle_cloud_ibs` is
+  `allocate`d without default initialization, and `s_reduce_ib_patch_array` copies the whole
+  struct into `patch_ib`, overwriting the defaults from
+  `s_assign_default_values_to_user_inputs`. Anything left unset reaches the solver as
+  uninitialized memory, and only where the allocation is not already zero-filled — a
+  garbage `v_blow` failed Frontier AMD with `ICFL is NaN` while every NVIDIA lane and all
+  local CPU/GPU runs passed. A platform-only NaN is the signature of this class.
 - Shared-state pattern: namelist declarations (`#:include 'generated_decls.fpp'`), the
   `eqn_idx`/`sys_size`/`b_size`/`tensor_size` state variables, and the common defaults
   core all live in `src/common/m_global_parameters_common.fpp`. Each per-target

docs/documentation/case.md

@@ -1129,9 +1129,17 @@ Note: For relativistic flow, the conservative and primitive densities are differ
 | `tau_star`        | Real    | Threshold stress for continuum damage model         |
 | `cont_damage_s`   | Real    | Power `s` for continuum damage model                |
 | `alpha_bar`       | Real    | Damage factor (rate) for continuum damage model     |
+| `reactive_burn`   | Logical | Enable condensed-phase reactive burn                |
+| `rburn%%k`         | Real    | Reactive-burn rate coefficient [1/s]                |
+| `rburn%%pign`      | Real    | Reactive-burn ignition pressure threshold [Pa]      |
+| `rburn%%pref`      | Real    | Reactive-burn reference pressure for the drive [Pa] |
+| `rburn%%n`         | Real    | Reactive-burn pressure-drive exponent               |
+| `rburn%%ta`        | Real    | Reactive-burn activation temperature [K] (0 = off)  |
 
 - `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).
 
+- `reactive_burn` converts a "reactant" fluid into a "product" fluid (`num_fluids = 2`, ``chemistry = 'F'``) via a programmed pressure burn `dlambda/dt = rburn%%k (1 - lambda) ((p - rburn%%pign)/rburn%%pref)^rburn%%n`. The two fluids share the same `gamma`/`pi_inf` and differ only in `qv`, so the conversion releases `qv` through the mixture EOS — a reactive-Euler/ZND detonation model on the diffuse-interface framework. It runs on the 5-equation (`model_eqns = 2`) and 6-equation (`model_eqns = 3`) multi-fluid models. Setting `rburn%%ta > 0` multiplies the rate by an Arrhenius factor `exp(-rburn%%ta/T)`, where `T` is the reactant phasic temperature, giving temperature-driven ignition instead of a pure pressure switch.
+
 ### 16. Cylindrical Coordinates
 
 When ``cyl_coord = 'T'`` is set in 3D the following constraints must be met:

docs/module_categories.json

@@ -22,6 +22,7 @@
             "m_viscous",
             "m_hb_function",
             "m_surface_tension",
+            "m_reactive_burn",
             "m_bubbles",
             "m_bubbles_EE",
             "m_bubbles_EL",

examples/1D_propellant_flame/case.py

@@ -0,0 +1,120 @@
+#!/usr/bin/env python3
+# 1D premixed H2/O2/Ar deflagration: a self-propagating premixed flame, used
+# as a stand-in for a solid-propellant burn front. A hot equilibrium-product
+# kernel at the left end ignites a flame that eats into the fresh premix.
+import argparse
+import json
+
+import cantera as ct
+
+parser = argparse.ArgumentParser(prog="1D_propellant_flame")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT", help="MFC toolchain state.")
+parser.add_argument("--length", type=float, default=0.008, help="Domain length [m].")
+parser.add_argument("--m", type=int, default=200, help="Number of grid cells.")
+parser.add_argument("--kernel", type=float, default=1e-3, help="Ignition-kernel width at the left end [m].")
+parser.add_argument("--ar", type=float, default=3.0, help="Ar dilution (H2:2,O2:1,AR:<ar>); higher = calmer/slower flame.")
+parser.add_argument("--tend", type=float, default=1.2e-4, help="Physical end time [s].")
+parser.add_argument("--cfl", type=float, default=0.35, help="Target acoustic CFL number.")
+parser.add_argument("--frames", type=int, default=10, help="Number of saved output frames (raise for a smooth movie).")
+args = parser.parse_args()
+
+ctfile = "h2o2.yaml"
+X = f"H2:2,O2:1,AR:{args.ar}"
+T0, P0 = 300.0, 101325.0
+
+# Fresh premixed reactants (fill the domain).
+fresh = ct.Solution(ctfile)
+fresh.TPX = T0, P0, X
+
+# Ignition kernel: constant-pressure (isobaric) equilibrium products of the
+# same mixture -> a thermodynamically self-consistent hot patch (T, rho, Y
+# all from Cantera) that ignites a flame without an initial pressure jump.
+burned = ct.Solution(ctfile)
+burned.TPX = T0, P0, X
+burned.equilibrate("HP")
+
+L = args.length
+Nx = args.m
+dx = L / Nx
+
+c_max = max(fresh.sound_speed, burned.sound_speed)
+dt = args.cfl * dx / c_max
+
+NT = int(args.tend / dt)
+NS = max(1, NT // args.frames)
+
+case = {
+    "run_time_info": "T",
+    # Domain
+    "x_domain%beg": 0.0,
+    "x_domain%end": L,
+    "m": Nx,
+    "n": 0,
+    "p": 0,
+    "dt": float(dt),
+    "t_step_start": 0,
+    "t_step_stop": NT,
+    "t_step_save": NS,
+    "t_step_print": NS,
+    "parallel_io": "T",
+    # Algorithm
+    "model_eqns": 2,
+    "num_fluids": 1,
+    "num_patches": 2,
+    "mpp_lim": "F",
+    "mixture_err": "F",
+    "time_stepper": 3,
+    "weno_order": 5,
+    "weno_eps": 1e-16,
+    "weno_avg": "F",
+    "mapped_weno": "T",
+    "mp_weno": "T",
+    "riemann_solver": 2,
+    "wave_speeds": 2,
+    "avg_state": 1,
+    # BCs: left reflective wall, right outflow
+    "bc_x%beg": -2,
+    "bc_x%end": -3,
+    # Flame is diffusion-controlled: viscous + species/thermal diffusion on
+    "viscous": "T",
+    "chemistry": "T",
+    "chem_params%diffusion": "T",
+    "chem_params%reactions": "T",
+    "chem_params%transport_model": 2,
+    "cantera_file": ctfile,
+    # Output
+    "format": "silo",
+    "precision": "double",
+    "prim_vars_wrt": "T",
+    "chem_wrt_T": "T",
+    # Patch 1: fresh premix, whole domain
+    "patch_icpp(1)%geometry": 1,
+    "patch_icpp(1)%x_centroid": L / 2,
+    "patch_icpp(1)%length_x": L,
+    "patch_icpp(1)%vel(1)": 0.0,
+    "patch_icpp(1)%pres": fresh.P,
+    "patch_icpp(1)%alpha(1)": 1.0,
+    "patch_icpp(1)%alpha_rho(1)": fresh.density,
+    # Patch 2: hot equilibrium-product ignition kernel, left end
+    "patch_icpp(2)%geometry": 1,
+    "patch_icpp(2)%alter_patch(1)": "T",
+    "patch_icpp(2)%x_centroid": args.kernel / 2,
+    "patch_icpp(2)%length_x": args.kernel,
+    "patch_icpp(2)%vel(1)": 0.0,
+    "patch_icpp(2)%pres": burned.P,
+    "patch_icpp(2)%alpha(1)": 1.0,
+    "patch_icpp(2)%alpha_rho(1)": burned.density,
+    # Fluid EOS (ideal-gas closure is bypassed by chemistry, but gamma/Re must be set)
+    "fluid_pp(1)%gamma": 1.0 / (1.4 - 1.0),
+    "fluid_pp(1)%pi_inf": 0.0,
+    "fluid_pp(1)%Re(1)": 1.0 / fresh.viscosity,
+}
+
+# Species mass fractions per patch + OH output (flame marker)
+for i in range(fresh.n_species):
+    case[f"patch_icpp(1)%Y({i + 1})"] = float(fresh.Y[i])
+    case[f"patch_icpp(2)%Y({i + 1})"] = float(burned.Y[i])
+case[f"chem_wrt_Y({fresh.species_index('OH') + 1})"] = "T"
+
+if __name__ == "__main__":
+    print(json.dumps(case))

examples/1D_reactive_burn/case.py

@@ -0,0 +1,124 @@
+#!/usr/bin/env python3
+# 1D condensed-phase detonation via programmed pressure burn (reactive_burn).
+# A dense stiffened-gas "reactant" is converted to a "product" fluid that shares
+# the same mechanical EOS (gamma, pi_inf) but has a lower reference energy qv --
+# so reactant -> product releases (qv_reactant - qv_product) per unit mass through
+# the mixture EOS, with no explicit energy source. A high-pressure initiator at the
+# left end launches a shock; where the shock raises the pressure above rburn%pign
+# the reactant burns, and the energy release sustains a self-propagating detonation
+# (von Neumann spike + Taylor rarefaction, near-CJ speed). Multi-fluid model,
+# chemistry OFF -- this deliberately bypasses the Cantera num_fluids=1 lock.
+import argparse
+import json
+
+parser = argparse.ArgumentParser(prog="1D_reactive_burn")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT", help="MFC toolchain state.")
+parser.add_argument("--res", type=float, default=1.0, help="Resolution multiplier: m=3000*res.")
+parser.add_argument("--tend", type=float, default=8.0e-6, help="Physical end time [s].")
+parser.add_argument("--frames", type=int, default=80, help="Number of saved output frames.")
+parser.add_argument("--model", type=int, choices=[2, 3], default=2, help="Multi-fluid model: 2 = 5-equation, 3 = 6-equation.")
+args = parser.parse_args()
+
+# 6-eq needs both phases present everywhere (its pressure-relaxation Newton solve is singular
+# on a zero-volume phase); seed a trace of product. 5-eq starts pure reactant.
+eps = 1.0e-6 if args.model == 3 else 0.0
+a1, a2 = 1.0 - eps, eps
+
+# --- Stiffened-gas EOS shared by reactant and product (mechanically identical) ---
+Gamma = 3.0  # physical adiabatic exponent
+Pi = 6.0e8  # stiffening pressure [Pa]
+gamma_p = 1.0 / (Gamma - 1.0)  # MFC gamma parameter
+pi_inf_p = Gamma * Pi / (Gamma - 1.0)  # MFC pi_inf parameter
+Q = 4.0e6  # heat of reaction [J/kg] (reactant qv; product qv = 0)
+
+rho0 = 1600.0  # reactant density [kg/m^3]
+p0 = 1.0e5  # ambient pressure [Pa]
+c0 = (Gamma * (p0 + Pi) / rho0) ** 0.5  # reactant sound speed
+
+# Domain
+L = 0.06
+m = int(3000 * args.res)
+dx = L / m
+x_init = 0.004  # initiator occupies the left 4 mm
+p_init = 5.0e9  # initiator pressure [Pa] -> launches an igniting shock
+
+# CJ-scale detonation speed sets dt.
+D_guess = 8000.0
+cfl = 0.3
+dt = cfl * dx / (D_guess + c0)
+NT = int(args.tend / dt)
+NS = max(1, NT // args.frames)
+
+case = {
+    "run_time_info": "T",
+    "x_domain%beg": 0.0,
+    "x_domain%end": L,
+    "m": m,
+    "n": 0,
+    "p": 0,
+    "dt": float(dt),
+    "t_step_start": 0,
+    "t_step_stop": NT,
+    "t_step_save": NS,
+    "t_step_print": max(1, NT // 20),
+    "parallel_io": "T",
+    # Algorithm
+    "model_eqns": args.model,
+    "num_fluids": 2,
+    "num_patches": 2,
+    "mpp_lim": "T",
+    "mixture_err": "T",
+    "time_stepper": 3,
+    "weno_order": 5,
+    "weno_eps": 1e-16,
+    "mapped_weno": "T",
+    "mp_weno": "T",
+    "weno_avg": "F",
+    "riemann_solver": "hllc",
+    "wave_speeds": "direct",
+    "avg_state": "arithmetic",
+    # Reactant reflects at the initiator wall; open at the far end.
+    "bc_x%beg": -2,
+    "bc_x%end": -3,
+    # Condensed-phase reactive burn
+    "reactive_burn": "T",
+    "rburn%k": 5.0e6,  # rate coefficient [1/s]
+    "rburn%pign": 5.0e8,  # ignition pressure threshold [Pa]
+    "rburn%pref": 1.0e9,  # reference pressure for the pressure drive [Pa]
+    "rburn%n": 1.0,  # pressure exponent
+    # Output
+    "format": "silo",
+    "precision": "double",
+    "prim_vars_wrt": "T",
+    # Patch 1: unreacted reactant fills the domain (fluid 1)
+    "patch_icpp(1)%geometry": 1,
+    "patch_icpp(1)%x_centroid": L / 2,
+    "patch_icpp(1)%length_x": L,
+    "patch_icpp(1)%vel(1)": 0.0,
+    "patch_icpp(1)%pres": p0,
+    "patch_icpp(1)%alpha_rho(1)": rho0 * a1,
+    "patch_icpp(1)%alpha_rho(2)": rho0 * a2,
+    "patch_icpp(1)%alpha(1)": a1,
+    "patch_icpp(1)%alpha(2)": a2,
+    # Patch 2: high-pressure initiator (reactant), left end
+    "patch_icpp(2)%geometry": 1,
+    "patch_icpp(2)%alter_patch(1)": "T",
+    "patch_icpp(2)%x_centroid": x_init / 2,
+    "patch_icpp(2)%length_x": x_init,
+    "patch_icpp(2)%vel(1)": 0.0,
+    "patch_icpp(2)%pres": p_init,
+    "patch_icpp(2)%alpha_rho(1)": rho0 * a1,
+    "patch_icpp(2)%alpha_rho(2)": rho0 * a2,
+    "patch_icpp(2)%alpha(1)": a1,
+    "patch_icpp(2)%alpha(2)": a2,
+    # Fluid EOS: reactant (1) and product (2) share gamma/pi_inf, differ only in qv.
+    "fluid_pp(1)%gamma": gamma_p,
+    "fluid_pp(1)%pi_inf": pi_inf_p,
+    "fluid_pp(1)%qv": Q,
+    "fluid_pp(2)%gamma": gamma_p,
+    "fluid_pp(2)%pi_inf": pi_inf_p,
+    "fluid_pp(2)%qv": 0.0,
+}
+
+if __name__ == "__main__":
+    print(json.dumps(case))