Keep species mass fractions in [0, 1] by construction instead of clipping

[sbryngelson]

Contribution Policy

All contributions are expected to demonstrate:

• A clear understanding of the codebase
• Alignment with product direction
• Thoughtful reasoning behind changes
• Evidence of real-world usage or hands-on experience with the problem

If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.

---

Description

The measurements below come from runs on Tuolumne.

With chemistry = 'T', MFC's operator-split reaction substep (reaction_substeps > 0) clamped negative mass fractions to zero, which hides negative mass fractions and adds mass (sum of rho*Y_k no longer equals rho). This PR replaces the clamp with a scheme that keeps Y in [0, 1] by construction:

1. Advection. Reconstructed face mass fractions are put on the simplex and limited with a Zhang-Shu scaling limiter (one theta shared by all species and both faces). HLLC's species flux is then Y_upwind times the mass flux, so species fluxes sum to the mass flux (Larrouturou 1991) (other Riemann solvers keep the sum consistency but not the upwind form).
2. Diffusion. The diffusion correction velocity carries upwind mass fractions.
3. Chemistry. alpha-QSS with alpha recomputed per stage keeps non-negative compositions non-negative, so the clamp is removed.
4. Fallback. Any cell still out of bounds at the end of an RK stage is cleaned as in PeleC clean_massfrac: rho*Y_k clipped to [0, rho], rho reset to their sum, velocity and specific energy kept. The number of fluid-cell cleanings (per cell per RK stage; IB cells are cleaned silently) is printed, only when the violation exceeds 1e-12 rho or the cell had no positive partial density and was left untouched.

Guarantees and their conditions

Cell averages provably stay non-negative when dt times each cell's outgoing mass flux per unit length (summed over directions) is at most w*rho, with w the reconstruction's Gauss-Lobatto endpoint weight (1/12 WENO5, 1/6 WENO3, 1/2 MUSCL); roughly an advective CFL <= w. At MFC's usual CFL (about 0.5), and for the diffusion terms, the limiter and upwinding are standard practice, not a proof; the fallback cleanup is the safety net. Explicit-diffusion stability itself is the job of chem_params%transport_dt (#1953, not yet on master).

What this does not fix

Species bounds do not cure the explicit-diffusion stability blow-up seen in the reacting curtain case; that is chem_params%transport_dt (#1953). The PR does not claim otherwise.

Evidence

Where negatives come from (probe on the curtain restart from step 5440, on the development branch, before this PR). Negatives first appear at the RK stage update (transport); minimum rho*Y was -3.39e-4 at step 0 of the restart. The IB ghost-cell update only carries them along, and the reaction clamp zeroed them so they never showed afterwards.

1D test (UUID 6CA07418, new species-bounds test). On master: min Y -2.9e-10, and sum(rho*Y) vs rho mismatch 3.9e-9. On this branch: min Y 0, mismatch 1.5e-13.

GPU curtain restart from step 5440, with transport_dt (2D, 38 moving reacting particles, 4 ranks): reaches 5.4 us. Min Y is 0 at every saved step, mass mismatch 8.9e-16. The fallback fired 4 times (1 cell each, |d rho| <= 5.1e-12). These runs use this branch merged into the development branch that carries moving reacting particles and transport_dt (not on master; one merge conflict, in s_compute_dt, resolved to the development side).

In the CPU test ibm_reacting_surface (F52F0D4C) the fallback reports on most steps (634 lines, 1 cell each, |d rho| ~4e-12): fluid cells at i = 1 (not IB ghost cells) whose violation is about 1.3e-12 of rho, just above the 1e-12 report threshold; these are roundoff-level.

Cost: on this restart (1407 steps, same step count) the average time per step is 0.309 s against 0.292 s for the unmodified build, about +6%.

GPU curtain restart without transport_dt still fails, and earlier than with the old clamp: about 110 vs 315 steps (split) and about 44 vs 214 (unsplit, CFL 0.7). Bounds hold (min Y >= 0, mass mismatch <= 3e-12, no cleanup lines) up to the failure. Bounds hold to the last saved step before the NaN, so the failure is not a positivity violation; with transport_dt both builds complete. Why this branch fails earlier than the clamp build without transport_dt has not been established. These runs use this branch merged into the development branch that carries moving reacting particles and transport_dt (not on master; one merge conflict, in s_compute_dt, resolved to the development side).

The full-length curtain run with transport_dt on this branch is in progress; I will post the result when it finishes.

Goldens (13 regenerated + 2 new, all chemistry = T)

Max relative change is for the quantity noted; all other tests unchanged. Full CPU suite after regeneration: 711 pass, 1 fail (903FA85A, a harness race at -j48 that passes alone), 50 skipped. After the final review fixes (cleanup counting, vacuum guard, docs, MUSCL test), the Chemistry suite and every regenerated golden were re-run individually and pass; no golden moved.

Test	UUID	Max change	Reason
1D inert shock tubes (6 tests)	6B3AE48F, F8ADA51B, 2BDE2018, 4D7D85B4, 7BBCE9E9, 259B2428	rel ~1e-9 (abs ~4e-5)	Limiter removes trace-species (<= 1e-8) undershoots at the contact
2D ibm_burning_grain	E5B66084	rel 6.8e-2 (abs 9.3e2 on energy)	Limiter at the burning wall, 0-to-nonzero species jump; >99% of the change within 3 cells of the wall
2D IBM Vieille Burn Rate	B317404C	rel 6.8e-2	Same as above
2D ibm_reacting_surface	F52F0D4C	rel 4.66 on a near-zero total-energy cell (abs 1.5e3 on energy of order 1e5; 17% at the worst cell)	Limiter; see below
2D ibm_flameholder	F200F862	~1.5e-11	Roundoff-level (regenerated unnecessarily, harmless)
Spatial Reacting Mixing Layer	56F8C4BC	abs <= 5e-2	Limiter at species interfaces
Reacting Mixing Layer	C4EB58A8	abs <= 5e-2	Limiter at species interfaces
MultiComponent Diffusion	41AD0140	max rel 3.7e-3 (on a 3.7e-6 value); max abs 1.6 (rel 3.8e-7)	Upwinded diffusion correction velocity
Species Bounds (new)	6CA07418	n/a	New in this PR
Species Bounds MUSCL (new)	2E7FF24B	n/a	New in this PR: MUSCL variant (min Y 0, mismatch 1.5e-13)

F52F0D4C bisected by commit. About 99% of the change comes from the limiter. At step 1 it is wall-concentrated (71% of the L1 difference within 3 cells of the wall), then propagates about 10 cells per step (RK3 x WENO5 reach), saturates from step 5 on (energy difference 1.2e3 to 1.5e3), with no NaN, no density sign flips and min rho 0.166. The diffusion, clamp-removal and cleanup changes add <= 1%. Regenerating at master reproduces the existing golden to 5e-7.

References

• Larrouturou, JCP 95:59 (1991), doi:10.1016/0021-9991(91)90253-H. Consistent species fluxes.
• Zhang & Shu, JCP 229:8918 (2010). Positivity-preserving scaling limiter for finite-volume WENO.
• Du & Yang, JCP 469:111548 (2022). Bounds-preserving multi-species reacting flow; clipping does not conserve mass.
• Mott, Oran & van Leer, JCP 164:407 (2000). alpha-QSS.
• PeleC clean_massfrac (Utilities.cpp).

Not verified

OpenACC/NVHPC builds (GPU tested only with OpenMP offload on MI300A); 3D flows and WENO3/WENO7 with species jumps; order of accuracy on smooth flow with the limiter; the full-length curtain run (in progress).

Mixed-precision (--mixed) builds: the build fails on master in M_DERIVED_TYPES independent of this PR, so mixed precision was not compiled or tested.

Type of change

• Bug fix / numerical robustness
• New feature
• Breaking change

Testing

• CPU test suite as above; new 1D tests 6CA07418 (WENO5) and 2E7FF24B (MUSCL) with a check on min Y and mass mismatch.
• GPU (Tuolumne, MI300A) curtain restarts as above.
• Changes CFD results for chemistry = 'T' cases; evidence for each regenerated golden above.

Checklist

• ./mfc.sh format and precheck pass
• Documentation updated (docs/documentation/case.md)

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_species_bounds.fpp	129	+129
src/common/m_chemistry.fpp	404	+9
src/simulation/m_time_steppers.fpp	909	+3
src/common/m_mpi_common.fpp	1496	+2
src/simulation/m_rhs.fpp	1970	+2

Directory	Lines	Diff
common	10489	+11
simulation	28637	+134
total	47659	+145