Updates to make HP queries return values instead of failing for higher pressures

[wandadars]

I was running a test case for a 1D counterflow diffusion flame at very high pressures using the Peng-Robinson equation of state, and the initial condition that used the HP equilibrate call return EoS failures. Codex 5.3 identified this as a possible reason for why the real/ideal gas mixed PREOS might be failing, and it did resolve the EoS that was happening for high pressure H-P equilibrate calls.

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. Examples: Output from agentic coding
  tools and/or substantial refactoring by LLMs (web-based or local).

For additional information on Cantera's AI policy, see
https://github.com/Cantera/cantera/blob/main/CONTRIBUTING.md#restrictions-on-generative-ai-usage -->

[codecov]

Codecov Report

❌ Patch coverage is 81.00000% with 19 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.66%. Comparing base (05de7fe) to head (3dafa16).
⚠️ Report is 33 commits behind head on main.

Files with missing lines	Patch %	Lines
src/thermo/PengRobinson.cpp	78.04%	5 Missing and 4 partials ⚠️
src/thermo/RedlichKwongMFTP.cpp	84.31%	2 Missing and 6 partials ⚠️
src/thermo/MixtureFugacityTP.cpp	75.00%	0 Missing and 2 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2087      +/-   ##
==========================================
+ Coverage   77.53%   77.66%   +0.12%     
==========================================
  Files         451      451              
  Lines       52810    53126     +316     
  Branches     8828     8883      +55     
==========================================
+ Hits        40947    41261     +314     
+ Misses       8888     8878      -10     
- Partials     2975     2987      +12     

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

[speth]

I think this may be partially on the right track, but I'd suggest focusing on a narrower definition of the problem. The equilibrium solver here is acting as a way of "finding" conditions that cause trouble for the equation of state solver, but isn't itself part of the problem.

The equations of the P-R model do not like the case where $a=b=0$. There are a lot of places where this results in 0/0, such that all the resulting thermodynamic properties are nan. For example:

import cantera as ct

def partial_lookup_peng_robinson_phase():
    """Build a mixed-parameter Peng-Robinson phase used by the original test."""
    species = ct.Species.list_from_file("test/data/co2_PR_example.yaml")
    phase_species = []
    for sp in species:
        data = dict(sp.input_data)
        if data["name"] in {"CH4", "O2"}:
            data.pop("equation-of-state", None)
        else:
            data["equation-of-state"] = {
                "model": "Peng-Robinson",
                "a": 0,
                "b": 0,
                "acentric-factor": 0.0,
            }
        phase_species.append(ct.Species.from_dict(data))

    return ct.Solution(
        thermo="Peng-Robinson",
        kinetics="bulk",
        species=phase_species,
        reactions=[],
        name="partial-lookup-pr",
    )

gas = partial_lookup_peng_robinson_phase()
gas.TPX = 500, 1e6, 'CO2:0.333, H2O:0.666'
gas()

which outputs:

  partial-lookup-pr:

       temperature   500 K
          pressure   1e+06 Pa
           density   6.4177 kg/m^3
  mean mol. weight   26.68 kg/kmol
   phase of matter   unspecified

                          1 kg             1 kmol     
                     ---------------   ---------------
          enthalpy   nan               nan              J
   internal energy   nan               nan              J
           entropy   nan               nan              J/K
    Gibbs function   nan               nan              J
 heat capacity c_p   nan               nan              J/K
 heat capacity c_v   nan               nan              J/K

                      mass frac. Y      mole frac. X     chem. pot. / RT
                     ---------------   ---------------   ---------------
               CO2           0.54984           0.33333   nan            
               H2O           0.45016           0.66667   nan            
     [   +5 minor]                 0                 0  

So it's pretty clear why the equilibrium solver can't do anything with such a mixture.

If you want to pick out specific cases where the cubic solver is causing trouble, I'd put a try/catch around the place in the equilibrium solver where it sets the mixture state and print out the full state that it's trying to set if that fails. You probably also want to do any exploration using a specific equilibrium solver (vcs or gibbs), as the automatic switching behavior may have the effect of hiding some errors that arise in the cubic equation of state solver.

[wandadars]

Good catch @speth . We have checks that watch for a=0, b=0 and use ideal-gas fallbacks for those quantities. Also those fallback were moved up to the base class since this is a common issue in both PREOS and RK cubic equations of state.

[speth]

Thanks for the updates here, @wandadars. I had a few thoughts about the contents of this PR.

• While there are lots of singularities as $b \to 0$, is there a reason to restrict evaluation of the full non-ideal properties in the case where $a=0$? I don't think there is any division by $a$ in any of these expressions. Though this may not really matter, as the only likely reason for either coefficient to be zero is if you have a phase with only "ideal" species.
• I'd like to explore a way of handling this case without introducing all of these new *_ideal methods. I think we can do so by having the structure of the actual property calculation methods be (1) evaluate the ideal contribution (2) check whether $b=0$ and return if so (3) otherwise, add the non-ideal terms. This would require adjusting where a couple of the existing terms are added in for some of the properties, but I don't think that's too difficult.
• I think this idea of only refining the physically-meaningful roots makes sense, but the exceptional case here -- only one real root which is non-physical -- is not exercised by the test suite. Can you come up with a case that would trigger this?
• Likewise, even with the changes to PengRobinson::densityCalc and RedlichKwongMFTP::densityCalc removed, the test suite still passes. Can you identify some conditions where this problematic behavior occurs? Or is that all cleared up by correctly handling the pure ideal case? Perhaps some phases very close to the ideal limit would be good ones to test.

[wandadars]

Thanks @speth . I like the idea of that ideal gas + departure gating in the cubic EoS. That's much better. I think the gating based on the value of b effectively prevents some of that nonphysical behavior that we were seeing.