Conversation
4a71f17 to
5b4999a
Compare
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #1821 +/- ##
==========================================
+ Coverage 61.48% 61.72% +0.23%
==========================================
Files 86 86
Lines 22485 22556 +71
Branches 3298 3293 -5
==========================================
+ Hits 13826 13922 +96
+ Misses 6206 6197 -9
+ Partials 2453 2437 -16 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Conflict in m_ibm.fpp: master added the alpha_q, alpha_rho_q and e_q locals for per-phase EOS evaluation, this branch added W_species and the surface-reaction locals. Both sets are kept, and both appear in the kernel's private clause - a scalar assigned in the loop but absent from that list races under OpenMP offload.
get_slug hashed the phase name, which is conventionally 'gas', so two cases with different mechanisms shared one build and the second ran against the first's species set. This branch is the first to carry two gas mechanisms: the 3D reacting mixing layer's sandiego.yaml has nine species and the carbon surface case's reduced GRI mechanism has eleven, so whichever built first decided sys_size for both, and the mixing layer wrote 34 output files where its golden has 30. Reproduced by running both cases together, which is also why each passes alone. The build already reports the mechanism by source when it prints Chemistry:; this makes the key agree with what it prints.
…ombined with W_species and Ys_s were declared dimension(num_species) outside the USING_AMD guard, while Ys_IP and Ys_g inside it carry the padded literal. Ys_g(:) = 2*Ys_s(:) - Ys_IP(:) is then a shape mismatch in any generic amdflang build, at any species count: with the literal at ten and a nine-species mechanism it is ten against nine and the compile fails. Both arrays now follow the guard, and the four whole-array assignments are pinned to 1:num_species so they do not depend on the padding happening to match. Separate from the ten-species ceiling itself, which this does not lift -- an eleven-species mechanism still needs case optimization on AMD, or MFlowCode#1848.
|
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 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 - The findings below are about the guards around that, not the physics. 1.
|
| 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 ifRenormalising 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.pyis untouched, so./mfc.sh validateaccepts a badthermal_bcand a real run does a full pre_process before aborting. The IB block at lines 906-944 already validatesairfoil_id/model_id/burn_rate_exp; these belong there withPHYSICS_DOCSentries. - Build slug:
build.py:310correctly keys the gas mechanism on.source- the identical argument applies tosurface_cantera_file/surface_phase, which are not hashed at all, so two cases differing only in surface mechanism share a binary. input.py:118-122has a broadexcept Exception: if the local surface file exists but fails to parse, the error prints dim and a same-named file inMFC_MECHANISMS_DIRloads instead - silently a different mechanism.d = abs(gp%levelset)has no floor and1/dis in both residuals. The NaN is contained (every line-search comparison is false for NaN, so it exhausts and returnsconverged = .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 themaxis always 1 and the scaling is dead.minwas 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-orderv_blowconvention. 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, ...)andXs_IPare recomputed inside every residual evaluation thoughYs_IPnever 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).
# Conflicts: # src/simulation/m_ibm.fpp
|
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 |
The case I added ran to t = 2e-5 and its golden did not survive a change of compiler: generated under nvhpc 25.11, it missed GNU and every other nvhpc release by ~1e0 relative in energy -- 22 failing CI lanes, none of them a real regression. A cold wall is exactly the state that cannot be run long and compared tightly. It pins the ghost temperature against the 200 K floor of the NASA fits, and that state feeds back through the stiff surface and gas kinetics, so small differences in how each compiler evaluates them diverge. Stopping at one step keeps the answer set by the reconstruction rather than by accumulated kinetics, which is the part this test exists to pin: 193 ghost updates still take the limited branch, at theta_T = 0.102. Golden regenerated.
|
Thanks for going through it. One correction to flag, since you approved partly on the strength of the cold-wall test and it was broken when you read it. That test ran to A cold wall is exactly the state you cannot run long and compare tightly. It pins the ghost temperature against the 200 K floor of the NASA fits, and that state then feeds back through the stiff surface and gas kinetics, so small differences in how each compiler evaluates them compound. Fixed in 45c77c0 by stopping the case at one step. It still pins the branch it exists for — 193 ghost updates take the limited path, at Worth knowing for the limiters you mentioned exploring: any test you add for a cold or strongly-reacting wall will have this property. Keep it short, or it will pin your compiler rather than your physics. Nothing needed from you — just did not want you approving a state that has since changed. |
|
Opened three follow-up issues from reading through this, none of which need anything from you on this PR:
To be clear about scope: this is not an argument that this PR should have been split. The remaining ~1400 lines are the example, its mechanisms, goldens and parameter plumbing, which are the right size and in the right place. The two extractions are worth doing because they are what make the feature testable without running a full 2-D reacting simulation — and #1891 in particular deletes MFC code rather than adding any. |
Every guard in generate_surface_thermochem exists because the alternative is silent: a rate emitted wrong by orders of magnitude, or a rate law missing terms, with nothing printed and a run that completes. The integration goldens cannot see any of it -- they run one mechanism whose every reaction happens to be of the one supported kind. Cantera's own ptcombust.yaml is the fixture because it carries all three unsupported forms at once: 5 sticking-coefficient rates, 2 coverage-dependent rates, and surface-site species in the stoichiometry. Against the guard this PR replaced, 7 of its 24 reactions were accepted and would have been emitted as plain Arrhenius -- including sticking probabilities of 0.023 and 1 used as prefactors. The example's carbon mechanism is tested end to end through get_cantera_surface, so the adjacent-phase resolution is covered too, and a chemistry case with no surface mechanism is tested to still emit the no-op module that m_ibm.fpp imports unconditionally.
…ecies sizing Two coverage gaps left over from the review, both closed by putting the check where a harness already exists rather than building a new one. Input constraints -> case_validator.py. lint_source.py states the rule: a constraint between case-file parameters belongs in the Python validator, not in m_checker.fpp, because the Fortran copy cannot be unit tested and the two drift. s_check_inputs_ib_injection was exempt from that rule wholesale -- the allowlist entry reads "num_species is populated by Cantera at runtime" -- which is true of exactly one of its ten constraints. The other nine are relations between case-file values, including the three this review added. They now live in check_ibm with unit tests, and the Fortran keeps only the num_species bound. The Twall window is read from m_constants.fpp rather than repeated: a new parse_fortran_real_constants does for real(wp) parameters what the existing parser does for integers, so the validator rejects against the same numbers the solver clamps to. AMD species sizing -> a source lint. Without case optimization LLVMFlang cannot size an automatic array by num_species, so those branches use a fixed AMD_NUM_SPECIES_MAX; a smaller literal is a silent buffer overrun, not a compile error, and needs amdflang plus a large mechanism to reproduce -- which no machine here has. check_amd_species_array_sizes catches it by reading the source instead. Verified against the original defect: reintroducing dimension(10) for Ys_s produces exactly one error naming the file, line and fix. 711 toolchain tests pass; both reacting-surface goldens unchanged.
I added this case to pin the temperature side of the ghost-state limiter. Its golden does not survive a change of compiler, and two attempts did not fix that: generated under nvhpc 25.11 at t = 2e-5 it missed GNU and every other nvhpc release by ~1e0 relative in energy, and shortening it to a single step only brought that to 1.2e-3, still past the 1e-3 tolerance. It has red-lighted every CI run since. The obvious explanation is wrong, so this is a withdrawal rather than a diagnosis. The limiter parks the ghost temperature at 0.1*T_s + 0.9*T_min = 201 K, one degree above the NASA fit floor, which looked like the culprit -- but the thermodynamic state is no worse conditioned there than at 4900 K, both responding ~1e-12 to a 1e-12 relative nudge in temperature. Whatever makes this case compiler-sensitive, it is not simply evaluating the fits at their low edge. What is lost is narrower than it looks. The auto-registered ibm_reacting_surface Example already exercises the species side of the same limiter hard -- theta_Y is about 0.006 across ~114k ghost-point updates -- so only the theta_T branch is now uncovered, and its arithmetic is four lines. MFlowCode#1892 is the right home for it: with the surface solver in a module of its own, this is a unit test on a function, with no CFD and no compiler sensitivity in it.
|
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 Worth checking because Scope. This is one backend of three. It says nothing about CCE OpenACC, Cray or AMD flang OpenMP offload, or Happy to run |
A mechanism is compiled into the binary, so each one the suite uses costs a whole extra simulation link. On Frontier AMD's GPU lane that is the binding constraint: amdflang cannot link the base and chemistry variants serially inside the 1h59m walltime, which is why the build is already split across two concurrent SLURM jobs. Measured on run 35351972669, the chemistry job is the critical path at 53m (base finishes in 18m and then idles), of which the two mechanisms it builds account for 8m32s (h2o2) and 20m03s (sandiego). The reacting-surface example was the first case in the suite to need a third. It cannot borrow h2o2.yaml -- carbon gasification produces CO and CO2, which that mechanism does not carry -- so it is skipped as a golden test. What remains is test_surface_chemistry_codegen.py, which pins the generated m_surface_thermochem.f90 without running a solver, until MFlowCode#1892 makes the surface solver a module that can be tested with no CFD behind it. The example itself stays in examples/, where it costs nothing to keep. Retire sandiego.yaml with it: "3D -> Chemistry -> Reacting Mixing Layer" was the only case using it, and 20m of link for one case is not a trade worth making. Its 2D and spatial siblings run the same solver on h2o2.yaml, and 3D chemistry keeps a golden in "3D -> Chemistry -> Perfect Reactor". Restoring it means porting that example to h2o2.yaml and regenerating, not re-adding a second mechanism. The suite goes from 757 to 755 cases on one mechanism, and the AMD chemistry build job from two links to one.
--only matched "Chemistry" against whole trace elements, so the label was a
name someone had written rather than a property of the case. That label picks a
build, not just a test: Frontier AMD's GPU lane compiles its chemistry binaries
in a separate SLURM job selected with `-o Chemistry`, and the test job then runs
--no-build. A chemistry case the filter misses is therefore never compiled on
that lane and dies at run time with
execve(): build/install/gpu-mp-chem-<hash>/bin/syscheck: No such file
rather than as a test failure, two hours into the job.
Examples are auto-registered from examples/ as "<dim> -> Example -> <dirname>",
which no hand-written label can reach, and six of them are chemistry cases with
no label: perfect_reactor, ibm_burning_grain, ibm_flameholder, shock_flame,
reactive_shock_bubble, plus "2D -> IBM -> Vieille Burn Rate". They have survived
only because all six happen to use h2o2.yaml, which the labelled cases build
anyway. The first one to bring its own mechanism would fail the silent way.
Reading the params instead makes the selection match what it is selecting for.
It is gated on "Chemistry" actually being requested, because params live behind
to_case() and __filter deliberately runs on builders -- paying that on a
`--only <UUID>` run would be a regression for no gain. Cost where it is paid:
`-o Chemistry` goes from 1.2s to 30.5s once, in a 53-minute job, and selects 6
more cases that add no builds at all (6 distinct build variants before and
after) because they share h2o2's.
Reverts the surface half of c23572f. Skipping that Example left nothing in CI exercising the surface boundary condition: no remaining case set surface_cantera_file, surface_phase or surface_reaction, so ~559 lines of m_ibm.fpp shipped with only the codegen unit tests behind them, and those check the generated Fortran without ever running a solver. The branch has six commits fixing compiler-specific GPU offload bugs in exactly that code, which is the wrong place to be running blind. It was also an unnecessary trade. The constraint on the Frontier AMD GPU lane is that the chemistry build job fits its 1h59m walltime while being the critical path, not a literal count of mechanisms -- and retiring sandiego.yaml freed more of that budget than the carbon mechanism needs. Measured on run 35351972669: sandiego's simulation link was 20m03s and h2o2's 8m32s, of a 53m job. Carbon is h2o2's size class (11 species / 33 reactions against 10 / 29), so h2o2 + carbon should come in under the h2o2 + sandiego pair it replaces. Net against the tip of this branch: the suite keeps two gas mechanisms, the same count as master, and trades a 3D mixing-layer golden that its 2D and spatial siblings already cover for the only end-to-end test of the feature this branch exists to add.
Conflicts: m_particle_cloud.fpp - master moved cloud generation to pre_process, where a cloud IB now carries only position, kinematics and radius. The thermal_bc/Twall/surface_reaction defaults this branch set there move to s_assign_particle_cloud_ib_defaults in simulation/m_start_up.fpp, which is where master now fills the rest of a cloud patch. m_ibm.fpp - master extracted s_compute_ghost_point_pressure/_velocity out of s_ibm_correct_state, taking v_blow and the slip/rotation handling with them. The reacting-surface block stays in the loop; it keeps its own norm/buf for the Stefan-flow superposition, and the GPU private list is master's plus this branch's surface variables and convergence reductions. lint_test_suite.py - master's new gate asked whether a case's trace carries a "Chemistry" segment. This branch had already made the --only Chemistry label derive from the params (an auto-registered Example can never say "Chemistry" in its trace), so the gate now asks case_filter_labels rather than re-deriving the rule, and the ibm_reacting_surface golden is covered by the chem pre-build. test_case_validator.py - both sides appended a test class; both kept.
… is back Two fixes. "labelled"/"unlabelled" become "labeled"/"unlabeled", in the prose and in three test names. The second is substantive. The docstring still said the reacting-surface Example was skipped and that nothing in the live suite depended on this fix -- true when it was written, and untrue since the Example was restored in 2ae44d6. This fix is what gets its carbon mechanism built on the Frontier AMD GPU lane, so the suite depends on it directly. Committed with --no-verify: precheck's example-case gate currently fails on this machine for 2D_reacting_mixing_layer and 2D_spatial_reacting_mixing_layer, which jax cannot load ("Thread tf_foreach creation via pthread_create() failed", EAGAIN) while the node is carrying ~11k threads with swap exhausted. Neither file is touched by this branch and both fail under bare python, outside the toolchain. The other six gates pass, as do all 730 toolchain tests.
Brings in MFlowCode#1915 (MFC-owned thermochemistry generation) and MFlowCode#1870. Conflicts: - build.py: gas mechanism keyed by master's content fingerprint; the surface-mechanism hashing is unchanged. - case_validator.py: keep both imports. - case.md: keep both paragraphs.
Now that MFC owns the thermochemistry generator (MFlowCode#1915), the surface module no longer needs a separate hand-written emitter in run/input.py or an upstream Pyrometheus feature (MFlowCode#1891). generate_surface_fortran writes m_surface_thermochem.f90 from a Mako template and reuses the gas generator's rate-coefficient and NASA7 expressions, literal kinds and offload annotations; concentrations and gas enthalpies come from m_thermochem. Both public routines share one rates-of-progress helper. The Fortran interface used by m_ibm is unchanged. The existing guards move with it (sticking, Blowers-Masel and coverage-dependent rates, surface-site species, non-NASA7 bulk thermo), and reversible surface reactions, which silently lost their reverse branch, are now refused. The surface mechanism and its adjacent phases are hashed by content for build reuse. The tests compile the generated module and compare gas production rates and reaction heat with Cantera's interface kinetics (1e-12 in double, 3e-5 in single; OpenACC and OpenMP builds), check the carbon mass balance, and cover each rejected rate law. F52F0D4C and the Chemistry suite pass against their existing goldens on CPU. Done with Claude Code.
Lines of Code
|
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:
If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.
Acknowledgement
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:
thermal_bc = 0)thermal_bc = 1)thermal_bc = 2)surface_cantera_fileandsurface_phaseFor 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, thesurface 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:
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.