Nonsense `ValueError: math domain error` and `ZeroDivisionError: float division by zero` thrown when `hot_start` is provided to `flash`

[alexchandel]

Attempting a CEOSGas + GibbsExcessLiquid[NRTL] flash with a hot_start to save time results a nonsensical ValueError: math domain error being thrown by stability_test_Michelsen in flash_vl.py. If no hot_start is provided, the error does not occur.

The traceback is provided below, and example code after. Note that the solver in flash_utils.py also malfunctions and throws an additional ZeroDivisionError: float division by zero.

The ValueError is similar to what I see if I call flash() with a moderately high temperature, like #flasher.flash(zs=zs, P=pres, T=800) in the example code provided. I suspect thermo is calculating dH_dT_P() at the guess, adding ∆H/dH_dT_P() to the temperature, obtaining a high temperature (due to skipping latent heats / vaporizations at intermediate temperatures), and flashing at this temperature. Then the flash crashes because of the noncondensable hydrogen.

$ ./thermo_error.py 
Traceback (most recent call last):
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_vl.py", line 812, in flash_TPV_HSGUA
    res, flash_convergence = self.solve_PT_HSGUA_NP_guess_bisect(zs, fixed_val, spec_val,
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_vl.py", line 939, in solve_PT_HSGUA_NP_guess_bisect
    sln_val = secant(to_solve, guess, xtol=self.TPV_HSGUA_BISECT_XTOL, ytol=ytol,
  File "/usr/local/lib/python3.9/site-packages/fluids/numerics/__init__.py", line 3161, in secant
    q1 = func(p1, *args, **kwargs)
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_vl.py", line 933, in to_solve
    res = self.flash(**kwargs)
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_base.py", line 275, in flash
    g, ls, ss, betas, flash_convergence = self.flash_TPV(T=T, P=P, V=V, zs=zs, solution=solution, hot_start=hot_start)
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_vl.py", line 746, in flash_TPV
    return self.flash_TP_stability_test(T, P, zs, self.liquid, self.gas, solution=solution)
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_vl.py", line 640, in flash_TP_stability_test
    stable, (trial_zs, appearing_zs, V_over_F, stab_guess_name, stab_guess_number, stab_sum_zs_test, stab_lnK_2_tot) = self.stability_test_Michelsen(
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_vl.py", line 549, in stability_test_Michelsen
    lnK = log(Ks[k])
ValueError: math domain error

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/alex/Desktop/./thermo_error.py", line 93, in <module>
    example()
  File "/Users/alex/Desktop/./thermo_error.py", line 89, in example
    print(add_heat(flasher, hrxn=48709, zs=zs, T0=temp, P=pres, extra_cp=extra_cp))
  File "/Users/alex/Desktop/./thermo_error.py", line 76, in add_heat
    eqs = flasher.flash(zs=zs, P=P, H=current_val + hrxn_vl, hot_start=eqs)
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_base.py", line 382, in flash
    raise e
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_base.py", line 376, in flash
    g, ls, ss, betas, flash_convergence = self.flash_TPV_HSGUA(fixed_var_val, spec_val, fixed_var, spec, iter_var, zs=zs, solution=solution, hot_start=hot_start)
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_vl.py", line 817, in flash_TPV_HSGUA
    g, ls, ss, betas, flash_convergence = self.solve_PT_HSGUA_NP_guess_newton_2P(zs, fixed_val, spec_val,
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_vl.py", line 874, in solve_PT_HSGUA_NP_guess_newton_2P
    sln = nonlin_spec_NP(guess, fixed_val, spec_val, zs, [xs, ys], [1.0-VF, VF],
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_utils.py", line 993, in nonlin_spec_NP
    sln = root(f_jac_numpy, guesses, tol=tol, jac=(True if jac else None), method=method, **solve_kwargs)
  File "/usr/local/lib/python3.9/site-packages/fluids/numerics/__init__.py", line 4093, in root
    return sp_root(*args, **kwargs)
  File "/usr/local/lib/python3.9/site-packages/scipy/optimize/_root.py", line 234, in root
    sol = _root_hybr(fun, x0, args=args, jac=jac, **options)
  File "/usr/local/lib/python3.9/site-packages/scipy/optimize/_minpack_py.py", line 226, in _root_hybr
    shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,))
  File "/usr/local/lib/python3.9/site-packages/scipy/optimize/_minpack_py.py", line 24, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
  File "/usr/local/lib/python3.9/site-packages/scipy/optimize/_optimize.py", line 76, in __call__
    self._compute_if_needed(x, *args)
  File "/usr/local/lib/python3.9/site-packages/scipy/optimize/_optimize.py", line 70, in _compute_if_needed
    fg = self.fun(x, *args)
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_utils.py", line 969, in f_jac_numpy
    ans = to_solve(flows_guess)
  File "/usr/local/lib/python3.9/site-packages/thermo/flash/flash_utils.py", line 895, in to_solve
    gi = trunc_log(xs[i]/xs_ref[i]) + lnphis[i] - lnphis_ref[i]
ZeroDivisionError: float division by zero

The following self-contained example reproduces the above error.

#!/usr/bin/env python

from thermo import ChemicalConstantsPackage, PRMIX, CEOSGas, FlashVL, \
    EquilibriumState, GibbsExcessLiquid
from thermo.interaction_parameters import IPDB
from thermo.nrtl import NRTL

def get_nrtl(zs: list[float], T: float = 298.15) -> NRTL:
    """Build NRTL."""
    tau_as = (
        (0,         -5.1549,    0,          0,      0),
        (5.8547,    0,          -0.40926,   0,      0),
        (0,         -0.39036,   0,          0,      0),
        (0,         0,          0,          0,      0),
        (0,         0,          0,          0,      0),
    )
    tau_bs = (
        (0,         2270.62,    284.966,    0,      0),
        (229.497,   0,          1479.46,    0,      0),
        (-216.256,  447.003,    0,          0,      0),
        (0,         0,          0,          0,      0),
        (0,         0,          0,          0,      0),
    )
    alpha_cs = (
        (0,     0.2,    0.3,    0,      0),
        (0.2,   0,      0.46,   0,      0),
        (0.3,   0.46,   0,      0,      0),
        (0,     0,      0,      0,      0),
        (0,     0,      0,      0,      0),
    )
    return NRTL(T=T, xs=zs, tau_as=tau_as, tau_bs=tau_bs, alpha_cs=alpha_cs)

def get_nrtl_flasher(
    names: list[str],
    zs: list[float] = None,
    T: float = 298.15,
    P: float = 101325,
) -> tuple[NRTL, FlashVL]:
    """Build NRTL flasher."""
    if zs is None:
        zs = [0.2, 0.2, 0.2, 0.2, 0.2]
    consts, props = ChemicalConstantsPackage.from_IDs(names)
    kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', consts.CASs, 'kij')
    eos_kwargs = {'Pcs': consts.Pcs, 'Tcs': consts.Tcs, 'omegas': consts.omegas, 'kijs': kijs}
    gas = CEOSGas(PRMIX, eos_kwargs=eos_kwargs, HeatCapacityGases=props.HeatCapacityGases)

    nrtl = get_nrtl(zs=zs, T=T)
    liquid = GibbsExcessLiquid(
        VaporPressures=props.VaporPressures,
        HeatCapacityGases=props.HeatCapacityGases,
        VolumeLiquids=props.VolumeLiquids,
        GibbsExcessModel=nrtl,
        equilibrium_basis='Psat', caloric_basis='Psat',
        T=T, P=P, zs=zs)

    flasher = FlashVL(consts, props, liquid=liquid, gas=gas)
    return flasher

def add_heat(flasher: FlashVL,
             zs: list[float],
             hrxn: float,
             T0: float,
             P: float = 101325,
             extra_cp: float = 0,) -> EquilibriumState:
    """Flash with Hrxn."""
    eqs: EquilibriumState = flasher.flash(T=T0, P=P, zs=zs)
    current_val = eqs.H()
    cp0 = eqs.dH_dT_P()

    hrxn_vl = hrxn * (cp0 / (cp0 + extra_cp))

    eqs = flasher.flash(zs=zs, P=P, H=current_val + hrxn_vl, hot_start=eqs)
    return eqs

def example():
    """Run."""
    stoich = {'nitrobenzene': -1, 'water': +2, 'aniline': +1, 'hydrogen': -3, 'methane': 0}
    flasher = get_nrtl_flasher(list(stoich.keys()))
    zs = [0, 0.81, 0.16, 0.02, 0.01]
    temp = 273.15+50
    pres = 1750000
    extra_cp = 0.06  # J/K / moles of fluid

    print(add_heat(flasher, hrxn=48709, zs=zs, T0=temp, P=pres, extra_cp=extra_cp))
    print(flasher.flash(zs=zs, P=pres, T=800))

if __name__ == '__main__':
    example()

[alexchandel]

The error appears to be related to the presence of noncondensables. Removing the hydrogen results in no error, either for add_heat or flasher.flash(zs=zs, P=pres, T=800). But noncondensables and inerts are necessary.

The issue is specifically in thermo's liquid phase calculations.

[alexchandel]

This probably related to #105. There is no documented way to add substances that don't participate in equilibrium (like hydrogen).

For example, removing the hot_start option, the 323.15K case calculates the H2 fraction at 2.88144893296719e-304. Much hotter and it probably underfllows to 0. This needs to be handled correctly, especially given #105.

[alexchandel]

I've opened #107 for the underlying liquid issue. However, the flasher also needs to be fixed, as during exception handling it divides by zero, and should not crash if the liquid phase does not exist.

[alexchandel]

This issue was present in thermo==0.2.14 chemicals==1.0.17, and likely earlier and later. But #107 and PH flashing with a hot_start work in thermo==0.2.18 chemicals==1.0.19, although TP with hot_start still does not (see #108).