Adding Cantera fluid model for incompressible reacting flows

[Cristopher-Morales]

Proposed Changes

This pull request aims to add a fluid model called FLUID_CANTERA for coupling Cantera library with SU2. Cantera library is an open-source library for problems involving chemical kinetics, thermodynamics, and transport processes. This is needed for combustion problems using the incompressible solver.

[UPDATE] adding the cantera option to meson is now sufficient to get cantera coupled to su2.

SU2 can be built with Cantera as follows (for example):

    ./meson.py build -Dwarning_level=2 -Denable-autodiff=false -Denable-directdiff=false -Dwith-mpi=enabled -Denable-cgns=true -Denable-tecio=false -Denable-cantera=true --prefix=/path_to_SU2_directory/SU2

• -Denable-tests=true can be used for testing if Cantera is correctly compiled with SU2, running one of the Cantera fluid unit test cases.

Work in progress:

• Improve cantera coupling, only compile the functionalities needed for SU2. Currently, when cantera is compiled, whole functionalities are compiled. Most of the cantera implementations are not needed for incompressible reacting flows.
• Compile cantera with the AD options. For this purpose, a different repository must be cloned and compiled, where the AD implementation of cantera is located.
• Add regression and unit test cases in the SU2 workflow.

Regarding the coupling between SU2 and Cantera, any advice/suggestion how to improve it would be really appreciated.

Thanks in advance!!

Related Work

Related to pull request #2426

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
• My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
• I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
• I have added a test case that demonstrates my contribution, if necessary.
• I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

[Hanquist]

How difficult would it be to extend this to the compressible solver?

[bigfooted]

How difficult would it be to extend this to the compressible solver?

Have a look at all the changed files that have "INC" in the name :-)
Cristopher can give a more detailed answer.

[joshkellyjak]

duplicated code, use helper function

[joshkellyjak]

duplicated from CIncNSSolver

[joshkellyjak]

logic repeated many times, make a config option

[joshkellyjak]

look into performance aspects of this

[joshkellyjak]

may also need jacobian here, check with tape tagging

[joshkellyjak]

comparing signed with unsigned

[joshkellyjak]

newSolution is already a shared_ptr

[joshkellyjak]

are all these mass fractions needed?

[joshkellyjak]

you don't need FlowIndices, nothing in the class uses it

[joshkellyjak]

space

[joshkellyjak]

I benchmarked how MAX_TRANSPORTED_SPECIES affects performance at this PR's head (01ab322). It does cost time, but not for the reason given in the comment. The cause is a full zero-fill of the per-edge scratch storage, and a small, targeted change removes it.

What slows down

case (time per iteration, whole solver)	cap 8	cap 12	cap 16	cap 20	cap 24
1 species eq., 1 rank	7.24 ms	+3.7 %	+8.0 %	+13.2 %	+18.5 %
2 species eq., 1 rank	7.73 ms	+2.6 %	+7.4 %	+13.1 %	+18.8 %
2 species eq., 4 ranks	2.25 ms	+7.8 %	+11.4 %	+17.2 %	+23.3 %
8 species eq., 1 rank	15.55 ms	+1.2 %	+4.2 %	+9.3 %	+11.9 %

The cases are species2_primitiveVenturi, species3_primitiveVenturi_flux_value, and the latter extended to 8 passive species. Each run is 200 iterations; the table shows the median of 5 runs. Final residuals are bit-identical across all caps.

Why

The cost grows with the square of the cap, and by about the same absolute amount whatever the number of species. Retired instructions grow too, so this is extra work, not cache effects. A profile at cap 24 shows memset called from CScalarSolver<CSpeciesVariable>::EdgeFluxResidual taking 17 % of the total run time (2 % at cap 8).

The source is the default constructor of C2DContainer:

C2DContainer() noexcept : Base() {}

Naming Base() in the initializer list value-initializes AccessorImpl. That base has no user-provided constructor, so its static m_data array is zero-filled. As a result, the whole EdgeResidual (4·Size² + 2·Size doubles: 2,176 B at cap 8, 18,816 B at cap 24) is cleared on every edge. This happens before its constructor zeroes the nVar × nVar part that is actually used, so the "loops stop at nVar" design never got the chance to save anything.

Stack frame size and stack probes are not the cause. The frame is set up once per call, not per edge. A test build with the cap-24 frame but no zero-fill ran at cap-8 speed.

Proposed fix

Add an opt-in constructor that skips the zero-fill, and use it only for EdgeResidual. Every other container keeps its current behaviour.

+/*!
+ * \brief Tag to construct a static-size C2DContainer without initializing its data.
+ */
+struct C2DUninitialized {};

   C2DContainer() noexcept : Base() {}
 
+  /*!
+   * \brief Static-size ctor that leaves the data uninitialized, unlike the default ctor which
+   *        value-initializes (zeroes) it. For hot-path temporaries that write before they read.
+   */
+  explicit C2DContainer(C2DUninitialized) noexcept {
+    static_assert(StaticRows != DynamicSize && StaticCols != DynamicSize, "Requires a static size.");
+  }

-  FORCEINLINE explicit EdgeResidual(size_t nEqn) : nVar(nEqn) {
+  FORCEINLINE explicit EdgeResidual(size_t nEqn)
+      : flux_i(C2DUninitialized{}),
+        flux_j(C2DUninitialized{}),
+        jac_ii(C2DUninitialized{}),
+        jac_ij(C2DUninitialized{}),
+        jac_ji(C2DUninitialized{}),
+        jac_jj(C2DUninitialized{}),
+        nVar(nEqn) {

This is safe because every read of an EdgeResidual is bounded by nVar, and the existing constructor loop still zeroes that region before anything accumulates into it. The full zero-fill was protecting nothing.

With the fix, compared with the current cap 8 (5 new runs per build, interleaved):

case	cap 8 + fix	cap 12 + fix	cap 24 + fix
1 species eq., 1 rank	−1.9 %	−2.9 %	−0.8 %
2 species eq., 1 rank	−1.7 %	−1.3 %	+0.3 %
2 species eq., 4 ranks	−2.0 %	−1.9 %	+1.9 %
8 species eq., 1 rank	−0.8 %	−0.9 %	−0.1 %

So with the fix, cap 12 is faster than today's cap 8, and cap 24 is within about 2 % of it, which is within run-to-run noise for the serial cases.

Correctness checks:

• The memset no longer appears in the edge loop's disassembly, and binary size is unchanged.
• On the three benchmark cases, final residuals are bit-identical to the original cap 8, on 1 and 4 ranks.
• 13 regression cases, run for 15 iterations each, give byte-identical history and restart files compared with the unpatched build. They cover compressible and incompressible SA/SST, MUSCL, species with turbulence, flamelet (steady and unsteady), axisymmetric, multizone, and 3D ONERA M6.
• Not tested: discrete-adjoint and OpenMP builds.

Related: the effective cap is 20

CSysMatrix::MAXNVAR = 20 sizes the matrix's stack buffers, and CSysMatrix::Initialize errors out for nVar > 20. A cap of 21–24 therefore compiles, but implicit runs with more than 20 species would fail at run time. It would be good to tie the two together, for example with static_assert(CSysMatrix<su2double>::MAXNVAR >= MAX_TRANSPORTED_SPECIES) (which needs MAXNVAR to be accessible), or to raise MAXNVAR together with the cap.

Suggestions

1. Apply the EdgeResidual fix in this PR or just before it. With it, the "larger values slow the kernel down" note no longer holds, and there is no performance reason to prefer 12 over 20.
2. Keep MAX_TRANSPORTED_SPECIES ≤ CSysMatrix::MAXNVAR, or raise both together.
3. Static dispatch for 1–4 equations (DispatchScheme<CScalarFlux_Species, 1, 2, 3, 4, Dynamic>) also hides the cost for small cases. It does nothing for 5+ species and triples the compile time of CSpeciesSolver.cpp (4.2 s → 13.3 s), so the fix above is the better option.

_{Machine: Apple M2 Pro, macOS 15.7, Apple clang 17.0.0, Open MPI 5.0.8, -O3 -march=native release build, OpenMP off. Time per iteration is the mean ITER_TIME over iterations 20–198. No hardware L1-miss counters were available on macOS. Profiles come from Instruments Time Profiler; instruction and cycle counts come from /usr/bin/time -l.}