Add heterogeneous reacting surface boundary conditions

[rocfire11]

Contribution Policy

We do not accept pull requests generated primarily by AI without genuine understanding or real-world usage context.

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.

---

Acknowledgement

• [ X] I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

---

PR template credit: junegunn

Summary

This PR adds heterogeneous reacting surface boundary conditions for immersed
boundaries. The implementation was developed for reacting carbon-particle
simulations in which heterogeneous surface chemistry is coupled to the
compressible gas-phase species equations.

The surface treatment supports:

• zero-normal-gradient temperature (thermal_bc = 0)
• prescribed surface temperature (thermal_bc = 1)
• coupled reacting-surface energy balance (thermal_bc = 2)
• heterogeneous species-flux boundary conditions
• Cantera-based heterogeneous surface mechanisms through
  surface_cantera_file and surface_phase

For a reacting surface, the species boundary condition balances diffusive
transport, Stefan mass flux, and heterogeneous surface production. One species
equation is replaced by the mass-fraction closure. When thermal_bc = 2, the
surface temperature is included as an additional Newton unknown and the
conductive and heterogeneous reaction heat fluxes are balanced.

Motivation

The immediate application is heterogeneous oxidation/gasification of carbon
particles using an immersed-boundary representation. The implementation is
intended to remain general with respect to the number of gas species and the
Cantera surface mechanism rather than hard-coding a particular carbon
mechanism.

Testing

The implementation has been tested using an 11-species reduced GRI-based gas-phase
mechanism together with a compatible heterogeneous carbon surface mechanism, including:

• prescribed surface temperature
• coupled surface species/energy solution
• MPI CPU/GPU simulations on HiPerGator

Both prescribed-temperature and coupled-energy carbon cases run successfully
and produce the expected surface reaction products, oxygen consumption,
thermal field, and reacting wake.

The branch was rebased from current MFC master before the surface changes were
introduced, and ./mfc.sh precheck and the simulation build pass.

[codecov]

Codecov Report

❌ Patch coverage is 56.97674% with 74 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.03%. Comparing base (e7139bc) to head (94ac09d).
⚠️ Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simulation/m_ibm.fpp	56.84%	55 Missing and 8 partials ⚠️
src/simulation/m_checker.fpp	15.38%	7 Missing and 4 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1821      +/-   ##
==========================================
- Coverage   61.26%   61.03%   -0.23%     
==========================================
  Files          84       84              
  Lines       22330    22529     +199     
  Branches     3265     3286      +21     
==========================================
+ Hits        13680    13751      +71     
- Misses       6207     6307     +100     
- Partials     2443     2471      +28     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_ibm.fpp	1543	+290
src/simulation/m_checker.fpp	83	+13
src/common/m_derived_types.fpp	472	+3
src/pre_process/m_global_parameters.fpp	453	+3
src/simulation/m_global_parameters.fpp	774	+3
src/simulation/m_particle_cloud.fpp	414	+3
src/pre_process/m_mpi_proxy.fpp	132	+2
src/simulation/m_mpi_proxy.fpp	525	+2

Directory	Lines	Diff
common	10356	+3
pre_process	4518	+5
simulation	28397	+311
total	46659	+319

[sbryngelson]

requires me to merge #1852 before this works on amd compilers

will make 'ready to review' when that merges. work on this pr can continue if needed (or ideally a follow up pr later)

[sbryngelson]

Review of 94ac09d. format and precheck are clean, so none of this is lint.

Worth saying first, because it shapes how to read the rest: the parts most likely to be silently wrong are right. Species indexing agrees end to end - gas_index is built from the main mechanism's species_names, which is the ordering m_thermochem, Ys_IP and eqn_idx%species%beg all use. The Stefan closure mdot_s = sum(W_k omega_k) correctly recovers the carbon mass crossing from the solid (reaction 5: 2*28.01 - 31.998 = 24.02 = 2 W_C). The Y_k * sum_BG diffusive-flux correction is the standard zero-sum form, and dropping species num_species for sum(Y) - 1 is legitimate because sum_k R_k vanishes at the root. norm = gp%levelset_norm points into the fluid, so the blowing sign is right. The GPU private clause is complete. molecular_weights is copied into a kernel local before being passed down, which is exactly what the nvfortran 23.11/24.1 trap requires. And s_add_cloud_particle sets all three new patch_ib members - the single thing our pitfalls file says is most often missed.

The findings below are about the guards around that, not the physics.

1. Ys_g = 2*Ys_s - Ys_IP goes negative, and your own example reaches it

src/simulation/m_ibm.fpp:316. The mirrored extrapolation is unclamped, so whenever the surface consumes a species faster than diffusion resupplies it, Ys_s < Ys_IP/2 and the ghost mass fraction is negative. It lands in a conserved variable at line 461.

Worked through with your mechanism, your Twall = 1200 K, Cantera mixture-averaged diffusivities, on the 25x25 grid the Example suite caps this case to (d ~ 2.8e-4 m):

species	k (m/s)	Da	Ys_s/Ys_IP	Ys_g/Ys_IP
OH	57.2	48	0.020	-0.959
O	118.1	97	0.010	-0.980
O2	0.032	0.041	0.961	+0.921

R1 and R2 have Ea = 0 and k ~ sqrt(T), so the radical channels are diffusion-limited at any surface temperature - this is not a high-temperature corner. sum(Ys_g) is still exactly 1 because the products compensate, so no global check catches it. The alpha-QSS clip at m_chemistry.fpp:271,288 hides the symptom, and only when reaction_substeps > 0.

Fix: clamp the mirror, or drop to first order where the mirror would leave the physical range. The second keeps the surface value exact:

if (2._wp*Ys_s(q) - Ys_IP(q) < 0._wp) then
    Ys_g(q) = Ys_s(q)          ! first-order ghost; the mirror would be unphysical
else
    Ys_g(q) = 2._wp*Ys_s(q) - Ys_IP(q)
end if

Renormalising afterwards is not enough on its own - the sum is already 1.

2. T_g = 2*T_s - T_IP has the same problem, and it bites thermal_bc = 1 too

src/simulation/m_ibm.fpp:296 (inert) and :317 (reacting). T_g feeds alpha_rho_IP(1) = alpha_IP(1)*pres_IP*mw_g/(gas_constant*T_g) and get_mixture_energy_mass(T_g, ...). A 1200 K wall under a 2500 K flame - the carbon-combustion case - gives T_g = -100 K: negative ghost density, and NASA polynomials evaluated at negative T (the a6/T term). The Newton solve clamps T_s internally to [200, 5000]; nothing clamps the extrapolation. Same one-line treatment as above, with a floor rather than zero.

3. thermal_bc and Twall are silently ignored without chemistry

The whole block at m_ibm.fpp:288 sits inside if (chemistry .and. patch_ib(patch_id)%inj_species == 0), but m_checker.fpp:104 was loosened from if (ib .and. chemistry) to if (ib), and docs/documentation/case.md:366,409 documents them as general IB thermal boundary conditions. So an isothermal cylinder in a plain Navier-Stokes IB run validates, documents as working, and is silently adiabatic. thermal_bc is read nowhere else in src/ - I grepped. Either gate the parameters on chemistry or lift the thermal branch out of the conditional; the docs should match whichever you pick.

4. The generator mis-emits sticking and reversible surface reactions

toolchain/mfc/run/input.py:236-247. append_reaction_rate guards only on hasattr(rate, "pre_exponential_factor"). Checked against the Cantera 3.1.0 in the project venv: ct.StickingArrheniusRate(0.1, 0, 0).pre_exponential_factor returns 0.1, so a sticking reaction passes the guard and its dimensionless sticking probability is emitted as an Arrhenius A - wrong by 1e4 to 1e12. Separately reaction.reversible is never consulted, so a mechanism written with <=> (Cantera's default for interface reactions) silently loses its reverse branch, and rate.coverage_dependencies is dropped.

Your own mechanism uses irreversible => with explicit orders and no sticking, so none of this affects your results - it is a trap for the next user. The generator already raises cleanly for surface-site species and unsupported thermo; three more raises in the same style would close it.

5. Two documented device-routine portability traps

s_surface_species_residual (m_ibm.fpp:1705) and s_surface_energy_residual (:1740) are GPU_ROUTINE(parallelism='[seq]') and call get_species_mass_diffusivities_mixavg / get_mixture_thermal_conductivity_mixavg / transitively get_species_enthalpies_rt from inside. .claude/rules/common-pitfalls.md records this exactly: calling get_species_* from inside a GPU_ROUTINE rather than from the kernel gave CCE OpenMP a runtime Memory access fault by GPU node-N on the first step, while every other backend ran. The build stays clean and only a case that reaches the path shows it, so NVIDIA and CPU lanes passing is not evidence. Every existing call site (m_chemistry.fpp:409-416) evaluates these in the GPU_PARALLEL_LOOP body and passes the arrays down - same restructuring here.

(I checked and discarded the related ftn-7066 concern: num_species is an integer, parameter in the pyrometheus output, so those bounds are constants, not device globals.)

6. Newton non-convergence is a silent BC switch, and it will make the golden flaky

m_ibm.fpp:1875,1898 return converged = .false., consumed at :321 and :454, and the ghost point silently reverts to an inert zero-flux surface - no counter, no warning, no output field. Beyond the diagnostic problem, converged is a discontinuous branch driven by a strict norm_trial < norm_R with no Armijo slack, so a one-ulp difference between compilers flips a ghost point between two O(1)-different boundary conditions. tests/F52F0D4C is an Example-tolerance golden on top of GRI-11 kinetics; we have already had to skip 1D_propellant_flame, 2D_hybrid_slab and 1D_flamelet for subtler drift than this. Expect it to be flaky across lanes.

Related: the Example test is not step-capped. cases.py:3177 caps t_step_stop only if "t_step_stop" in case, and this example uses cfl_adap_dt with t_stop/t_save, so the golden is the full run - ~1e3 steps of GRI-11 with a 30-iteration Newton (12 residual evaluations each, each O(Ns^2)) per ghost point per RK stage. Several adaptive-dt examples are already in casesToSkip for this. A small step-capped dedicated test would serve better than the full-run Example golden.

7. Smaller

• No Python-side validation: case_validator.py is untouched, so ./mfc.sh validate accepts a bad thermal_bc and a real run does a full pre_process before aborting. The IB block at lines 906-944 already validates airfoil_id/model_id/burn_rate_exp; these belong there with PHYSICS_DOCS entries.
• Build slug: build.py:310 correctly keys the gas mechanism on .source - the identical argument applies to surface_cantera_file/surface_phase, which are not hashed at all, so two cases differing only in surface mechanism share a binary.
• input.py:118-122 has a broad except Exception: if the local surface file exists but fails to parse, the error prints dim and a same-named file in MFC_MECHANISMS_DIR loads instead - silently a different mechanism.
• d = abs(gp%levelset) has no floor and 1/d is in both residuals. The NaN is contained (every line-search comparison is false for NaN, so it exhausts and returns converged = .false.) but that lands in finding 6. max(d, small) is clearer than relying on NaN comparison semantics.
• dx = fd_eps_Y*max(abs(Ys_s(j)), 1._wp) (:1841) - mass fractions are <= 1 so the max is always 1 and the scaling is dead. min was probably intended.
• The pivot test pivot_value <= epsilon(1._wp) (:1935) is absolute on a matrix with O(1e7) entries; it will never fire. A relative test against the column norm would.
• Convention: the scalars use second-order mirroring while the Stefan velocity is added once to vel_g (:426), matching the first-order v_blow convention. If reconstruction sees the wall-normal velocity as ~v_stefan/2, the convected mass flux is half what the species BC assumed. May well be deliberate, but the two halves of one surface condition using different ghost conventions deserves a comment.
• get_mixture_molecular_weight(Ys_IP, ...) and Xs_IP are recomputed inside every residual evaluation though Ys_IP never changes; the transport coefficients could be frozen at the outer Newton level too.

Suggested order

1 and 2 are the ones that produce wrong answers rather than crashes, and both are a few lines. 3 is a one-line gate. 4 is three raises matching the ones already there. 5 is the restructuring m_cbc/m_ibm already use. 6 and 7 argue together for a cheap step-capped test instead of the full-run Example golden.

Separately: your Frontier AMD CPU build failures are not yours to fix - they are the 11-species mechanism against the dimension(10) literal under the USING_AMD guard, which #1852 raises. I cherry-picked #1852 onto 94ac09d locally and the tree builds clean on MI210 with F52F0D4C passing, so that lane clears when #1852 lands. The gpu-omp [2/2] lane was a bad node (syscheck failed three times on frontier10212).