Add heterogeneous reacting surface boundary conditions - #1821

Codecov Report

❌ Patch coverage is 53.92157% with 94 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.55%. Comparing base (ed7a238) to head (13fb4e0).
⚠️ Report is 3 commits behind head on master.

@@            Coverage Diff             @@
##           master    #1821      +/-   ##
==========================================
- Coverage   62.82%   61.55%   -1.27%     
==========================================
  Files          86       86              
  Lines       22394    22721     +327     
  Branches     3305     3329      +24     
==========================================
- Hits        14070    13987      -83     
- Misses       6071     6279     +208     
- Partials     2253     2455     +202     

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

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)

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

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

As a newbie, is there any required of me right now before merging with master branch? Your steps 3-7 will be addressed in a separate PR since I do not want to corrupt this one.

you can leave this here as is, thanks @rocfire11

Ran this on a GPU, which nothing in the test suite does — every reacting-surface golden is CPU-only, so the offload path had never been executed.

It works, and the answer is bit-identical to CPU.

Case: the 2D_ibm_reacting_surface example at cells_per_D = 16, Tstar = 0.05, parallel_io = F, single rank, one GPU.

Worth checking because s_solve_surface, s_blend_ghost_state and s_pick_bath_species are all GPU_ROUTINE(parallelism='[seq]') called from inside the ghost-point parallel loop, and that pattern has two documented traps in this codebase: CCE rejects such a call when it is nested inside another routine seq, and nvfortran 23.11/24.1 segfault when a parameter array from m_thermochem is passed into a declare-target routine. Neither bites here — W_species is copied from molecular_weights into a local before the call, which is the form that works.

Scope. This is one backend of three. It says nothing about CCE OpenACC, Cray or AMD flang OpenMP offload, or --gpu mp on nvfortran — and CCE is the one with the sharpest edges. But it is one lane verified where there were none, and the bit-identical result means the surface solve is not accumulating any device/host divergence.

Happy to run --gpu mp on nvfortran too if that is useful; the AMD and Cray lanes need CI.

@rocfire11 — I pushed two commits to this branch to clear the conflicts and the red CI. Both touch your surface-chemistry code, so please check they do what you intended.

1. 151119c7 — merge resolution

Three of the conflicts needed a judgement call about your code:

m_ibm.fpp. Master extracted s_compute_ghost_point_pressure and s_compute_ghost_point_velocity out of s_ibm_correct_state, taking the v_blow blowing term and the slip/rotation handling with them. I left your reacting-surface block in the loop and kept the Stefan-flow superposition on top of whatever master's helper returns:

! Calculate velocity of ghost cell
call s_compute_ghost_point_velocity(gp, patch_id, radial_vector, vel_IP, pres_IP, vel_g)

if (chemistry .and. patch_ib(patch_id)%inj_species == 0 .and. patch_ib(patch_id)%surface_reaction == 1 &
    & .and. surface_converged) then
    norm(1:3) = gp%levelset_norm
    buf = sqrt(sum(norm**2))
    if (buf > 0._wp) vel_g = vel_g + v_stefan*norm/buf
end if

Master deleted the norm/buf declarations when it moved that code out, so I declared them locally again for this block. Worth your eye: v_blow now gets applied inside master's helper, before your v_stefan is added here. Previously both were superposed in the same place. If the two are ever active together, confirm the ordering is still what you want.

m_particle_cloud.fpp → m_start_up.fpp. Master moved cloud generation to pre_process, where a cloud IB now carries only position, kinematics and radius. Your thermal_bc / Twall / surface_reaction defaults had nowhere to live there, so they moved to s_assign_particle_cloud_ib_defaults in src/simulation/m_start_up.fpp, which is where master now fills in the rest of a cloud patch:

ib_patch%slip = .false.
! Particles are inert surfaces: a cloud IB carries no case-file surface condition, so the thermal,
! reaction and blowing fields must be set here rather than left as whatever patch_ib held.
ib_patch%thermal_bc = 0
ib_patch%Twall = 0._wp
ib_patch%surface_reaction = 0
ib_patch%v_blow = 0._wp

lint_test_suite.py. Master added a gate that asks whether a case's trace carries a "Chemistry" segment. Your f38a1c6f had already made the --only Chemistry label derive from the params instead, which is what lets an auto-registered Example be covered at all. I pointed the gate at case_filter_labels rather than letting it re-derive the rule, so ibm_reacting_surface stays in the chem pre-build instead of being skipped.

2. fcc5b734 — regenerated tests/F52F0D4C/golden.txt

The latest master merge pulled in #1792, which splits s_ibm_correct_state into interpolate-then-apply so a ghost point's image-point stencil can no longer read a cell another ghost point already overwrote. That moves the answer wherever those stencils overlap — on this case, everywhere near the surface. #1792 regenerated the one master golden it moved (127A967A, mibm_cylinder_in_cross_flow); this case is the same class but had no master golden, so nothing caught it.

It is not FP noise — the pre-regeneration failure reproduced with identical values on GNU/CPU locally and on NVHPC 23.11 / 25.11 / 26.1 in CI (var 246, abs 1.59e-03, rel 2.44e-03), so the case is compiler-stable well inside its 1e-3 tolerance.

Magnitude of the change across 119152 entries:

16% of entries move past 1e-3 and 4% past 1e-2, concentrated near the surface with the far field essentially unchanged.

This is the part I most want you to check. I regenerated the golden because the mechanism is understood and reproducible, but I can't tell you whether the new near-surface field is the physically better one — that's your surface chemistry. #1792 fixed a genuine race, so it ought to be an improvement, but please confirm the new solution looks right before this merges rather than taking the regenerated numbers on faith.

Merged master into this branch (e53c1572) to clear the conflict. The only conflict was a comment in toolchain/mfc/test/cases.py: both sides drop the 3D_reacting_mixing_layer test, and I kept master's wording. The second merge brings in #1914 and #1930 with no conflicts. It builds on Frontier on CPU, OpenMP offload and OpenACC. (Merged with Claude Code.)

Lines of Code

Regarding item 1 (m_ibm.fpp, ordering of v_blow and v_stefan):

I checked this carefully and the ordering is fine (both by hand, and by a new AI chat). Before the merge, v_blow and v_stefan were added consecutively to the already-constructed IBM ghost velocity. The merged version performs the same operations, except the v_blow addition now occurs at the end of s_compute_ghost_point_velocity, followed immediately by the existing v_stefan addition in the reacting-surface block. There is no intervening velocity operation, so the result is unchanged. I also checked the Stefan sign: mdot_s is positive for net gas production and the circle levelset_norm points outward, so positive carbon gasification gives outward Stefan flow.

As an additional check, I ran the reacting carbon-cylinder case with the same fixed timestep for 100 steps using the MFC_05 version and the current PR. The initial Silo fields are bit-for-bit identical. After 100 steps, the combined relative differences are $L_1=1.51\times10^{-14}$ and $L_2=6.28\times10^{-14}$. Thus, the two versions produce essentially identical numerical solutions, with differences at the level of accumulated floating-point roundoff.

I will address your point 2 tomorrow.