CalebBell / thermo

Thermodynamics and Phase Equilibrium component of Chemical Engineering Design Library (ChEDL)
MIT License
594 stars 114 forks source link

`'FAILED from hot start TP'` when `hot_start` is the exact solution #108

Closed alexchandel closed 2 years ago

alexchandel commented 2 years ago

When calculating a PT flash, and providing ~an accurate hot_start guess that is the correct P and within 1K of the true solution~ the exact solution, thermo's flasher prints FAILED from hot start TP and seems to ignore the provided guess. The error is emitted here on L742 with no indication of why the accurate guess was rejected.

If the hot_start is ~close but not exact~ the exact solution, why does it fail? Note this does not occur for a PH flash, just PT.

The following code replicates this:

#!/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 = [1.0/len(names)] * len(names)
    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)

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

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))

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

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  # this guess is accurate
    pres = 1750000

    hot_start = flasher.flash(zs=zs, P=pres, T=temp)
    print((hot_start.T, hot_start.H()))

    # use exact solution:
    res = flasher.flash(zs=zs, P=pres, T=temp, hot_start=hot_start)  # prints 'FAILED from hot start TP' !!!
    print(res)
    print((res.T, res.H()))

if __name__ == '__main__':
    example()
CalebBell commented 2 years ago

Hi Alex,

I have added a warning to the documentation saying not to use hot_start if you prefer convergence. https://thermo.readthedocs.io/thermo.flash.html#thermo.flash.Flash.flash Maybe one day this won't be the case but it is today. These algorithms are very complex and at this time I am not doing development on the hot_start feature.

Caleb

alexchandel commented 2 years ago

@CalebBell What exactly is going on to cause this issue? Why does providing a near-solution result in "FAILED"?

I need performance, and would like to use a close initial guess if it's faster than searching the entire space from no guess.

alexchandel commented 2 years ago

@CalebBell I've narrowed this issue, and it's almost certainly a bug now. The TP flash fails when given the exact solution, but the PH does not. Could you please reopen?