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

* simulation: fix IBM+chemistry hang by setting consistent species/energy at immersed-boundary points

* examples: enable IBM+chemistry in flame-holder now that the hang is fixed

* examples: add hybrid-rocket fuel slab and 1D propellant flame (solid-fuel analogs)

* simulation: add burning-grain surface fuel injection to IBM (v_blow + inj_species)

* simulation: pressure-coupled IBM burn rate (Vieille's law) for internal-ballistics feedback

* examples: make flameholder a proper reacting bluff-body (fuel-injecting cylinder, wake flame)

* examples: add --tend knob to burning-grain case

* simulation: add condensed-phase reactive burn (multi-fluid programmed detonation)

* simulation: harden reactive_burn preconditions and IBM burn-rate against negative pressure

* test: add reactive_burn regression test (condensed programmed burn)

* examples: fix hybrid slab run time and viscosity (flame now anchors)

* examples: add 3D spherical burning grain with octant symmetry

* example: reactive shock-flame RM roll-up via hardcoded IC 273

* test: golden files for combustion example and reactive-burn tests (generated on master base)

* checker: require rburn_pref > 0 and forbid chemistry + Euler bubbles

Addresses PR review (Copilot):
- reactive_burn divides the pressure drive by rburn_pref; prohibit rburn_pref <= 0
- s_interpolate_image_point selects the bubbles/QBMM branch before chemistry, so the
  Ys_IP species state is not carried when both are on; forbid chemistry + bubbles_euler/qbmm
  (Fortran m_checker + Python case_validator) rather than touch the hot IBM path

* checker: constrain reactive_burn to model_eqns=2 and burn_rate_exp >= 0

Addresses independent review findings on #1670:
- reactive_burn's exact volume-fraction swap and single-pressure read assume the
  5-equation pressure-equilibrium model, but num_fluids=2 alone also admits model_eqns=3;
  prohibit model_eqns /= 2 (m_checker + case_validator.check_reactive_burn + PHYSICS_DOCS)
- Vieille's-law v_blow ~ p^n diverges for a negative exponent as surface pressure -> 0;
  prohibit patch_ib%burn_rate_exp < 0 (case_validator.check_ibm)

* feat(reactive-burn): 6-equation model support + Arrhenius temperature kinetics

Extends the condensed-phase reactive burn beyond the 5-equation model to the 6-equation (Saurel/Pelanti) multi-fluid model, and adds an optional temperature-driven ignition term.

6-equation support: the reactant->product qv release is carried through the phasic-pressure relaxation. m_pressure_relaxation.fpp threaded the formation energy (alpha_rho_k*qv_k) through s_equilibrate_pressure and s_correct_internal_energies so the phasic pressures and the recovered mixture pressure are qv-consistent with s_compute_pressure -- otherwise a nonzero qv inflated the phasic pressure by rho_k*qv_k/gamma_k. The existing 5-eq burn source then works unchanged on model_eqns=3. For qv=0 the new terms are exactly zero, so all 13 existing 6-eq relaxation goldens are bit-identical (verified CPU/OMP/ACC).

Arrhenius kinetics: new rburn_ta (activation temperature [K]); default 0 keeps the pure pressure-driven rate bit-identical. When rburn_ta>0 the rate gains a exp(-rburn_ta/T_r) factor, with T_r=(p+pi_inf)/((gamma-1)*cv*rho) the reactant phasic temperature. Guarded (checker + case_validator) so rburn_ta>0 with cv=0 is rejected instead of silently disabling the factor.

Test: adds a model_eqns=3 reactive-burn regression golden (nested off the 5-eq test), and hardens both reactive goldens by giving the uniform case a small advection velocity so the momentum output is well-resolved rather than near-zero roundoff (the old golden compared ~0 momentum at 1e-12 and failed across recompiles/backends). Validated on CPU, OpenMP-offload, and OpenACC.

* test(reactive-burn): add Arrhenius (rburn_ta>0) regression golden

Nested off the reactive-burn detonation test with rburn_ta=500 and cv=1500, chosen so the exp(-rburn_ta/T_r) factor is O(0.4) at the IC temperature -- exercising the temperature branch (all other reactive goldens use rburn_ta=0). Verified on CPU, OpenMP-offload, and OpenACC.

* refactor(reactive-burn): group rburn params into a derived type; address review panel + CI

Post-review hardening from a 5-reviewer panel, plus the CI-surfaced flame-example fragility.

rburn derived type: the five rburn_* scalars (k, pign, pref, n, ta) become a rburn type instance, matching the sibling chem_params pattern (type in m_derived_types, TYPED_DECLS registration, member _r/_nv, descriptions regex, auto GPU_DECLARE + GPU_UPDATE(device=[rburn]) in m_global_parameters and m_start_up so the type is device-resident). Numerically inert (~1 ULP recompile roundoff); reactive goldens regenerated for the namelist-key change. Verified CPU/OMP/ACC.

Completeness guards (review): rburn%k (>0), rburn%pign (set), rburn%pref (>0), rburn%n (>=0) are now validated in both m_checker and case_validator. Each defaults to the dflt_real sentinel, so an unset value silently produced spurious ignition (pign), NaN via drive**n (n), or a backward reaction (k).

IB inj_species bound (review): m_checker now requires patch_ib%inj_species <= num_species when chemistry + ib; it indexes Ys_IP(1:num_species) at image points, so an out-of-range value was an out-of-bounds write.

Docs (review): reactive_burn requires num_fluids = 2 (checker aborts otherwise), corrected from >= 2 in case.md and the module docstring.

CI: skip the 1D_propellant_flame example (golden removed). Its finite-rate flame front drifts past the 1e-3 Example tolerance across compilers (nvhpc 25.11 golden vs 24.3/GNU/Intel/CCE/AMD) -- same non-portability class as the flame examples already skipped. The reactive_burn detonation goldens remain portable (they passed every CI lane).

* fix(reactive-burn): broadcast the rburn derived type to non-root ranks

The rburn refactor moved the five rburn_* scalars (auto-broadcast) into a derived type, but derived-type members are not auto-broadcast -- they are hand-emitted per type in fortran_gen.py (as for chem_params/lag_params). rburn had no emitter, so on a multi-rank reactive_burn run every non-root rank kept rburn at the dflt_real sentinel: a silent wrong-physics bug the single-rank golden tests cannot catch. Adds _emit_rburn (all members REAL -> mpi_p, under the reactive_burn guard, sim target), mirroring _emit_chem_params. Verified in the regenerated generated_bcast.fpp.

* test(reactive-burn): add a 2-rank MPI regression test for the rburn broadcast

All other reactive goldens run single-rank, so a broken rburn broadcast is invisible to them. This ppn=2 variant of the detonation case burns rank 1's half of the domain with the broadcast rburn; verified to FAIL when the rburn broadcast emitter is removed and PASS with it, so it locks the fix in CI.

* test(ibm): add a golden for the Vieille's-law pressure-coupled burn rate

The ibm_burning_grain/ibm_flameholder Example goldens already cover surface injection (v_blow + the inj_species composition swap), but both leave burn_rate_exp at 0, so the pressure-coupled scaling v_blow*(p/burn_rate_pref)^n shipped with no regression coverage. Registers the same example with --burn_exp 0.5 at the Example suite's 25x25/50-step size, so it differs from the existing golden only by the coupling: turning it on moves the solution by ~5e-3 relative (1.3 absolute in energy), well above the 1e-3 tolerance, so a broken multiplier fails the test. Tolerance pinned to 1e-3 to match the sibling Example golden -- compute_tolerance would otherwise pick 1e-10 via the ib branch, which a 50-step IBM + finite-rate-chemistry run cannot hold across compilers.

* test: skip the 2D_hybrid_slab example golden (not portable across compilers)

The diffusion flame over the fuel slab is positioned by accumulated stiff-kinetics roundoff, so by step 50 the transverse momentum drifts past the 1e-3 Example tolerance: nvhpc lanes pass, Intel and Frontier CCE disagree with the golden by ~2e-3 absolute (up to 44% relative on a small component). Same non-portability class as 1D_propellant_flame and the other flame examples already skipped there, so no single golden passes every lane. The case remains a runnable example; the reactive_burn, shock_flame, ibm_burning_grain and ibm_flameholder goldens are unaffected and pass on all lanes.

* fix(ibm): initialize the burning-surface fields on particle-cloud IB patches

particle_cloud_ibs is allocated without default initialization, and s_add_cloud_particle set only the geometry/motion members. s_reduce_ib_patch_array then copies the whole struct into patch_ib, overwriting the zeroed defaults from s_assign_default_values_to_user_inputs -- so v_blow, inj_species and burn_rate_exp/pref reached the solver as uninitialized memory.

That was latent until this branch added the first read of those members in m_ibm: a garbage v_blow > 0 superimposes a garbage wall-normal velocity on the ghost points, which NaNs the field. It reproduced only where the allocation is not already zero-filled -- Frontier AMD failed 2D_mibm_particle_cloud with 'ICFL is NaN' at step 2 while every NVHPC lane (and local CPU/OMP/ACC) passed. Sets the four fields for the inert particle surfaces; audited that these are exactly the patch_ib members m_ibm reads (airfoil_id/model_id/length_x remain unset but are never read for circle/sphere particles).

* docs: record the derived-type broadcast and patch_ib init traps in common-pitfalls

Both bit this branch, both fail silently, and neither is visible to local testing: a derived-type param regrouping dropped its MPI broadcast (non-root ranks kept the sentinel, invisible to single-rank goldens), and an uninitialized patch_ib member reached the solver through the whole-struct copy from particle_cloud_ibs (garbage v_blow -> ICFL NaN, reproducing only on Frontier AMD where the allocation is not zero-filled). Documents the required steps and the platform-only-NaN signature.

Author: sbryngelson

.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))