Skip to content

Fix MultiPhaseEquil constant HP/SP solve interaction with temperature bounds#2126

Draft
speth wants to merge 2 commits into
Cantera:mainfrom
speth:fix-268-multiphase-HP
Draft

Fix MultiPhaseEquil constant HP/SP solve interaction with temperature bounds#2126
speth wants to merge 2 commits into
Cantera:mainfrom
speth:fix-268-multiphase-HP

Conversation

@speth
Copy link
Copy Markdown
Member

@speth speth commented Jun 1, 2026

Changes proposed in this pull request

Fix two issues affecting the "Gibbs" MultiPhaseEquil solver:

In MultiPhase::equilibrate_MultiPhaseEquil, the HP temperature search constructs MultiPhaseEquil e(this, strt) and then runs the equilibration inside a try block. But the constructor was placed outside the try block. When the T-search guesses a high temperature where liquid water's thermo is invalid while H2O(L) still has moles, the constructor throws with start=false. Because construction was outside the try, the catch — which exists precisely to retry with strt=true (recompute the composition) — never ran, so the error escaped to the user. The SP case had the identical latent bug, which is also fixed here.

Fixing the previous error exposed that the solver was returning an incorrect solution that did not satisfy element conservation. The start=true recovery path (added in 2219f65 to fix #160) redistributes an excluded species' atoms by adding b_missing[m] moles directly to component species m. That's only correct when the component formula matrix B is the identity — i.e. monatomic components, as in the Issue #160 Pb/O/C test. For the polyatomic components used in #268 (CO₂, H₂O, N₂…) it effectively adds B·b_missing instead of b_missing. This is fixed by computing b_missing after computeN() (which can reorder elements) and solving B·δ = b_missing for the correct per-component increment.

If applicable, fill in the issue number this pull request is fixing

If applicable, provide an example illustrating new features this pull request is introducing

The following equilibrium problem can now be successfully solved:

import cantera as ct

liquid = ct.Solution('water.yaml', 'liquid_water')
gas = ct.Solution('gri30.yaml')
gas.TPX = 300, 101325, 'O2:1.0, N2:3.76, CH4:0.5'
mix = ct.Mixture([(gas, 1.0), (liquid, 2.0)])
mix.T = 300
mix.P = 101325
mix.equilibrate('HP', solver='gibbs')
mix()

which prints:

************ Phase gri30 ************
Moles:  2.5703389291228995

  gri30:

       temperature   362.81 K
          pressure   1.0133e+05 Pa
           density   0.73081 kg/m^3
  mean mol. weight   21.757 kg/kmol
   phase of matter   gas

                          1 kg             1 kmol     
                     ---------------   ---------------
          enthalpy       -8.1844e+06       -1.7807e+08  J
   internal energy        -8.323e+06       -1.8109e+08  J
           entropy            9331.9        2.0304e+05  J/K
    Gibbs function        -1.157e+07       -2.5173e+08  J
 heat capacity c_p              1511             32874  J/K
 heat capacity c_v            1128.8             24560  J/K

                      mass frac. Y      mole frac. X     chem. pot. / RT
                     ---------------   ---------------   ---------------
                O2        4.4984e-11        3.0587e-11           -48.948
               H2O           0.56711           0.68491           -103.33
               CO2          0.074806          0.036982           -159.54
                N2           0.35808           0.27811           -24.376
     [  +49 minor]        1.5717e-18        1.0713e-18  

************ Phase liquid_water ************
Moles:  0.42966107086394384

  liquid_water:

       temperature   362.81 K
          pressure   1.0132e+05 Pa
           density   1000 kg/m^3
  mean mol. weight   18.015 kg/kmol
   phase of matter   unspecified

                          1 kg             1 kmol     
                     ---------------   ---------------
          enthalpy       -1.5596e+07       -2.8096e+08  J
   internal energy       -1.5596e+07       -2.8096e+08  J
           entropy            4703.2             84728  J/K
    Gibbs function       -1.7302e+07        -3.117e+08  J
 heat capacity c_p            4222.4             76067  J/K
 heat capacity c_v            4222.4             76067  J/K

                      mass frac. Y      mole frac. X     chem. pot. / RT
                     ---------------   ---------------   ---------------
            H2O(L)                 1                 1           -103.33

AI Statement (required)

  • Extensive use of generative AI. Significant portions of code or documentation were generated with AI, including logic and implementation decisions. All generated code and documentation were reviewed and understood by the contributor. Implemented using Claude Code (Opus 4.8).

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
  • AI Statement is included
  • The pull request is ready for review

speth and others added 2 commits June 1, 2026 13:09
When the guessed temperature is outside the bounds of a phase, construction
of the MultiPhaseEquil object raises an exception. Therefore, this needs to
be done within the 'try' block to force a new temperature guess.

Fixes Cantera#268

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The previous handling of this redistribution (introduced in Cantera#2116 / 2219f65)
did not correctly handle polyatomic species.

Partially addresses Cantera#268

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
@codecov
Copy link
Copy Markdown

codecov Bot commented Jun 1, 2026

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 77.79%. Comparing base (9f82790) to head (53e34ab).

Additional details and impacted files
@@           Coverage Diff           @@
##             main    #2126   +/-   ##
=======================================
  Coverage   77.78%   77.79%           
=======================================
  Files         452      452           
  Lines       53602    53608    +6     
  Branches     8958     8961    +3     
=======================================
+ Hits        41693    41702    +9     
+ Misses       8883     8881    -2     
+ Partials     3026     3025    -1     

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Failure in Gibbs multiphase equilibrium solver at constant enthalpy+pressure Multi-Phase Equilibria

1 participant