MixtureFugacityTP::solveCubic fails for certain temperatures and pressures

[speth]

Problem description

For certain combinations of T and P, the method MixtureFugacityTP::solveCubic fails, issuing a warning like:

MixtureFugacityTP::solveCubic(T = 271.8109054527264, p = 19614005.835764904): WARNING root didn't converge V = 2174.003159636111

and leaving the phase in an invalid state with pressure and density at inf and -inf (or vice versa), which in turn results in all thermodynamic properties evaluating to NaN.

Steps to reproduce

Conditions where this error occurs are relatively rare, but occur on a set of related points.

import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats

yaml = """
phases:
- name: gasphase
  thermo: Redlich-Kwong
  state: {T: 300.0, P: 1 atm}

species:
- name: H2
  composition: {H: 2}
  thermo:
    model: NASA7
    temperature-ranges: [200.0, 1000.0, 3500.0]
    data:
    - [2.34433112, 7.98052075e-03, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12,
      -917.935173, 0.683010238]
    - [3.3372792, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14,
      -950.158922, -3.20502331]
  equation-of-state:
    model: Redlich-Kwong
    units: {length: m, pressure: Pa, quantity: mol}
    a: 0.14470695468544512
    b: 1.8439567747927248e-05
"""

gas = ct.Solution(yaml=yaml)

TT = np.linspace(250, 350, 2000)
PP = np.linspace(0.1e7, 4.0e7, 2400)
DD = np.zeros((len(TT), len(PP)))

for i, T in enumerate(TT):
    for j, P in enumerate(PP):
        gas.TP = T, P
        DD[i,j] = gas.density

fig, ax = plt.subplots()
T_fail = []
P_fail = []
for i, j in zip(*np.where(~np.isfinite(DD))):
    T_fail.append(TT[i])
    P_fail.append(PP[j])

ax.plot(T_fail, P_fail, 'ko')
ax.plot(TT, res.intercept + res.slope * TT)
res = scipy.stats.linregress(T_fail, P_fail)

ax.set(xlabel='T [K]', ylabel='P [Pa]')
ax.set_ylim(ymin=0)

Behavior

These are points where the solver fails and issues a warning:

The error occurs along a line of points, approximately defined by the line $P=157285712-506093T$. The failures become less common as temperature goes up.

System information

• Cantera version: 3.1.0a2 / main branch at d37a76b
• OS: Windows 11
• Python/MATLAB/other software versions: Python 3.11

Additional context

Originally reported on the Users' Group

[bryanwweber]

How interesting! Is the line the same for other fluids like O2?

[ischoegl]

Is there any connection to #1157?

[wandadars]

That is strange because up in the supercritical region there, the root finder shouldn't have any issues converging like it would potentially do in a subcritical case with a liquid/vapor region to contend with.

[ischoegl]

I added a band-aid fix in #1819 to address some UG comments.

As an aside, I believe we should consider using a 'canned' root solver, e.g. boost's cubic_roots. Based on a cursory google search, there are other alternatives as well, e.g. ACM TOMS Algorithm 954. Unfortunately, there aren't any comments regarding what algorithm is currently used. Edit: Hiding in plain sight: Nickall's method. (Still) tagging @decaluwe and @gkogekar ....