Energy equation for reactors using PlasmaPhase #2068

mquiram · 2025-12-10T14:07:11Z

Changes proposed in this pull request

Added Joule heating and elastic power loss terms to energy equation RHS in all reactor classes when thermo: plasma
Added PlasmaPhase methods for cp_mole(), intEnergy_mole(), entropy_mole(), and gibbs_mole() which account for a separate electron temperature. Currently, they are commented out and inherited from IdealGasPhase so that the energy equation can work and all the tests can pass. Not sure if we want to maintain inheritance from IdealGasPhase, skip the test, or maybe change the h = u + pv test for PlasmaPhase to something like h = u + pv + X_e(R(T_e - T)) and similarly change the Gibbs consistency test to account for non-equilibrium electron temperatures. Edit: these methods are now uncommented

Checklist

The pull request includes a clear description of this code change
Commit messages have short titles and reference relevant issues
Build passes (scons build & scons test) and unit tests address code coverage
Style & formatting of contributed code follows contributing guidelines
The pull request is ready for review

codecov · 2025-12-10T14:37:44Z

Codecov Report

❌ Patch coverage is 28.28283% with 71 lines in your changes missing coverage. Please review.
✅ Project coverage is 76.67%. Comparing base (9e50aea) to head (c738296).

Files with missing lines	Patch %	Lines
src/thermo/PlasmaPhase.cpp	48.83%	21 Missing and 1 partial ⚠️
src/equil/MultiPhase.cpp	50.00%	6 Missing ⚠️
src/zeroD/ConstPressureMoleReactor.cpp	0.00%	5 Missing ⚠️
src/zeroD/ConstPressureReactor.cpp	0.00%	5 Missing ⚠️
src/zeroD/IdealGasConstPressureMoleReactor.cpp	0.00%	5 Missing ⚠️
src/zeroD/IdealGasConstPressureReactor.cpp	16.66%	5 Missing ⚠️
src/zeroD/IdealGasMoleReactor.cpp	0.00%	5 Missing ⚠️
src/zeroD/IdealGasReactor.cpp	0.00%	5 Missing ⚠️
src/zeroD/MoleReactor.cpp	0.00%	5 Missing ⚠️
src/zeroD/Reactor.cpp	0.00%	5 Missing ⚠️
... and 2 more

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2068      +/-   ##
==========================================
- Coverage   76.75%   76.67%   -0.09%     
==========================================
  Files         457      457              
  Lines       53699    53779      +80     
  Branches     9121     9134      +13     
==========================================
+ Hits        41219    41235      +16     
- Misses       9375     9442      +67     
+ Partials     3105     3102       -3

☔ View full report in Codecov by Sentry.
📢 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.

BangShiuh

@mquiram Thanks for working on this. I left some comments for the changes in class PlasmaPhase.

BangShiuh · 2025-12-17T13:12:40Z

src/thermo/PlasmaPhase.cpp

+        }
+        return m_eedfSolver->getElectronMobility();
+    } else {
+        throw NotImplementedError("PlasmaPhase::electronMobility",


Electron mobility can be supported without Boltzmann-two-term if we provide the cross-section data and the EEDF. The equation is different according to the bolsig document. I am fine with the current implementation. But I would make another PR to improve this so that the plasma reactor can be more flexible. In addition, IonGasTransport.::getMobilities can use the electron mobility. We can cache the electron mobility value in PlasmaPhase so we don't need to call getElectronMobility() when it stays the same value.

BangShiuh · 2025-12-17T13:51:15Z

src/thermo/PlasmaPhase.cpp

+    // sigma = e * n_e * mu_e   [S/m];   q_J = sigma * E^2   [W/m^3]
+    const double mu_e = electronMobility();    // m^2 / (V·s)
+    if (!(mu_e > 0.0) || !std::isfinite(mu_e)) {
+        return 0.0;
+    }
+    const double ne = concentration(m_electronSpeciesIndex) * Avogadro; // m^-3
+    if (!(ne > 0.0) || !std::isfinite(ne)) {
+        return 0.0;
+    }
+    const double E  = electricField(); // V/m
+    if (!std::isfinite(E) || E == 0.0) {
+        return 0.0;
+    }
+    const double sigma = ElectronCharge * ne * mu_e; // S/m
+    return sigma * E * E; // W/m^3
+}


I don't think checking std::isfinite for each term is needed. The error should be caught when calling the function. We can check if any term is or smaller than zero and return zero..

You're right. This, along with lines 577-579 is a sort of debugging artifact, and I agree that it's now unnecessary looking back at it.

src/thermo/PlasmaPhase.cpp

… in sample

…ng single temperature

… tests g_eq_sum_gk_Xk

speth

Thanks, @mquiram, I think this is a good step towards making PlasmaPhase work more smoothly with other parts of Cantera.

I think the new equilibrate behavior makes sense, but I'd like to suggest making the implementation a little more generic. Instead of introducing the setLockTeToT method on PlasmaPhase, we could instead introduce a pair of virtual methods on the ThermoPhase class that are always called at the start / end of the equilibrium solver, for example beginEquilibrate and endEquilibrate. The equilibrium "lock" class would then just call this method for every phase, with the special PlasmaPhase behavior being implemented in its overrides of these methods. While there aren't currently any other uses for this, I could imagine other phase types implementing some special behavior during an equilibrium solve, for example skipping some dependent property lookups or adjusting some internal tolerances, so this feature might come in use again.

Likewise, I think we want to avoid having to do conditional dynamic_casts to PlasmaPhase within the reactor classes, and introduce a virtual method on the base ThermoPhase class that gets called no matter what (and returns 0 for most phase types). I'm not sure what the best name for this method would be -- something like intrinsicHeating, with Joule heating being one example.

speth · 2025-12-24T03:10:50Z

src/zeroD/IdealGasConstPressureMoleReactor.cpp

-    if (m_thermo->type() != "ideal-gas") {
+    /* if (m_thermo->type() != "ideal-gas") {
        throw CanteraError("IdealGasConstPressureMoleReactor::initialize",
                           "Incompatible phase type '{}' provided", m_thermo->type());
-    }
+    } */


Rather than disabling this check, it would be better to make it allow both ideal-gas and plasma as the type() to avoid letting other ThermoPhase models through.

speth · 2025-12-24T03:26:32Z

src/thermo/PlasmaPhase.cpp

+        if (!m_eedfSolver) {
+            throw CanteraError("PlasmaPhase::electronMobility",
+                "EEDF solver is not initialized.");
+        }


I think this case is impossible, given that m_eedfSolver is initialized in the PlasmaPhase constructor.

speth · 2025-12-24T03:33:04Z

src/thermo/PlasmaPhase.cpp

+    if (m_electronEnergyDist.size() != m_nPoints
+        || m_electronEnergyDistDiff.size() != m_nPoints) {
+        writelog("PlasmaPhase::elasticPowerLoss: "
+             "EEDF not initialized, returning 0.\n");
+        return 0.0;
+    }


I think this needs to be a bit more firm about when it does/doesn't work. If the Boltzmann-two-term mode is active, then an empty EEDF should trigger an update. And otherwise, calling this method should be an presumably be an error (like for electronMobility).

speth · 2025-12-24T03:33:42Z

src/thermo/PlasmaPhase.cpp

        rate += concentration(m_targetSpeciesIndices[i]) *
                m_elasticElectronEnergyLossCoefficients[i];
    }
+    const double q_elastic = Avogadro * Avogadro * ElectronCharge * concentration(m_electronSpeciesIndex) * rate;


Please try to stick to our style recommendation of 88 characters per line maximum.

speth · 2026-01-04T23:18:28Z

test/thermo_consistency/consistency.cpp

-    EXPECT_NEAR(g, phase->mean_X(gk), atol);
+    if (phase->type() == "plasma") {
+        // PlasmaPhase is a two-temperature, non-equilibrium model. The mixture Gibbs
+        // free energy is defined via g = h - T*s using the gas temperature T, while
+        // the electron chemical potential is defined using Te and pe. As a result,
+        // the single-temperature identity g = sum_k X_k * mu_k is not enforced for
+        // PlasmaPhase, so this particular consistency check is skipped.
+        GTEST_SKIP() << "PlasmaPhase: g = sum_k X_k * mu_k identity not enforced";
+    } else {
+        // Standard single-temperature identity
+        EXPECT_NEAR(g, phase->mean_X(gk), atol);
+    }


Should the Gibbs free energy be calculated from the chemical potentials instead, using the species-specific temperature as is done for $h_k$ and $s_k$, so that there isn't a discrepancy with using the gas temperature for the portion of the mixture $g$ that comes from the electrons? That may make other tests fail, but it may be preferable to resolve those rather than fail to satisfy $g = \sum \mu_k X_k$.

speth · 2026-01-04T23:26:49Z

src/equil/MultiPhase.cpp

 }

+struct TeLockGuard {
+    std::vector<std::pair<PlasmaPhase*, double>> items;


No need to use the std:: prefix for vector or pair -- we have using directives for these and other common types in ct_defs.h.

speth · 2026-01-05T01:44:42Z

src/equil/MultiPhase.cpp

        throw CanteraError("MultiPhase::equilibrate_MultiPhaseEquil",
                           "No convergence for T");
    } else if (XY == SP) {
+        TeLockGuard te_guard(*this); // also single-T during SP


You may want to move this up to the more general MultiPhase::equilibrate method, as it will then be activated for both the MultiPhaseEquil (aka "Gibbs") solver and the VCS solver. You should probably also implement something on this front for the ChemEquil solver (aka "element potential").

speth · 2026-01-05T14:22:22Z

src/thermo/EEDFTwoTermApproximation.cpp

-            throw CanteraError("WeaklyIonizedGas::converge", "Convergence failed");
+            throw CanteraError("WeaklyIonizedGas::converge", "Convergence failed; returning last iterate.\n");
+            return;


No need for the unreachable return after throw. Also, watch the line length.

speth · 2026-01-05T14:34:00Z

src/thermo/PlasmaPhase.cpp

+double PlasmaPhase::cp_mole() const
+{
+    return IdealGasPhase::cp_mole();
+}


If this doesn't do anything other than what the base class does, there's no need for an explicit override. I do think the class docstring could use a bit more documentation about how this and the definitions for the other properties result in a consistent set of thermodynamic properties (or even what the definition of $c_p$ is when there isn't a single temperature)

speth · 2026-01-05T14:36:05Z

src/thermo/PlasmaPhase.cpp

+
+double PlasmaPhase::intEnergy_mole() const
+{
+    std::vector<double> ubar(m_kk);


Might be worth having a temporary m_work array as a member of the class to avoid repeatedly allocating this vector.

ischoegl added the Reactors label Dec 11, 2025

BangShiuh suggested changes Dec 17, 2025

View reviewed changes

mquiram added 4 commits January 4, 2026 17:19

fixed energy equation RHS for reactors using PlasmaPhase

6a46169

addressed redundant finite checks in PlasmaPhase, turned on energy eq…

257f2f5

… in sample

enable HP and SP equilibrate when thermo: plasma by temporarilyassumi…

f87cca1

…ng single temperature

implemented 2T thermo consistency test for h_eq_u_plus_Pv and skipped…

c738296

… tests g_eq_sum_gk_Xk

speth force-pushed the plasma-energy-eq branch from 74dc550 to c738296 Compare January 4, 2026 22:19

speth requested changes Jan 5, 2026

View reviewed changes

speth mentioned this pull request Jan 5, 2026

Cantera does not compute a multi-temperature ideal gas law for now #1872

Open

Uh oh!

Energy equation for reactors using PlasmaPhase #2068

Are you sure you want to change the base?

Energy equation for reactors using PlasmaPhase #2068

Conversation

mquiram commented Dec 10, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov bot commented Dec 10, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

BangShiuh left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

speth left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

mquiram commented Dec 10, 2025 •

edited

Loading

codecov bot commented Dec 10, 2025 •

edited

Loading