[ISSUE] low level routine that worked in v7 fails in v8

[longemen3000]

in Julia, using the CoolProp library at v8.0.0, the following test in CoolProp.jl fails:

handle = AbstractState_factory("HEOS", "Water&Ethanol")
pq_inputs = get_input_pair_index("PQ_INPUTS")
t = get_param_index("T")
AbstractState_set_fractions(handle, [0.4, 0.6])
# Update state with P = 101325 Pa, Q = 0 (saturated liquid at 1 bar)
AbstractState_update(handle,pq_inputs,101325, 0)
mole_frac = [0.0, 0.0]
AbstractState_get_mole_fractions(handle, mole_frac)
gas_frac = [0.0, 0.0]
AbstractState_get_mole_fractions_satState(handle, "gas", gas_frac)
liq_frac = [0.0, 0.0]
AbstractState_get_mole_fractions_satState(handle, "liquid", liq_frac)
i = 0 #index of water in mixture
AbstractState_get_fugacity(handle, i)
AbstractState_get_fugacity_coefficient(handle, i)

the error is the following:

ERROR: LoadError: CoolProp: Error: fugacity_coefficient is not well-defined in the two-phase region; evaluate on SatL or SatV instead
Stacktrace:
  [1] error(::String, ::String)
    @ Base ./error.jl:54
  [2] raise(errcode::Base.RefValue{Int64}, message_buffer::Vector{UInt8})
    @ CoolProp ~/src/julia/pkgdev/CoolProp/src/CoolProp.jl:746
  [3] AbstractState_get_fugacity_coefficient(handle::Int64, i::Int64)
    @ CoolProp ~/src/julia/pkgdev/CoolProp/src/CoolProp.jl:982
  [4] top-level scope
    @ ~/src/julia/pkgdev/CoolProp/test/testLow.jl:79
  [5] include(mapexpr::Function, mod::Module, _path::String)
    @ Base ./Base.jl:307
  [6] top-level scope
    @ ~/src/julia/pkgdev/CoolProp/test/runtests.jl:55
  [7] include(mapexpr::Function, mod::Module, _path::String)
    @ Base ./Base.jl:307
  [8] top-level scope
    @ none:6
  [9] eval(m::Module, e::Any)
    @ Core ./boot.jl:489
 [10] exec_options(opts::Base.JLOptions)
    @ Base ./client.jl:283
 [11] _start()
    @ Base ./client.jl:550

The error originates from here:

CoolProp/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp

Lines 3470 to 3501 in b72bc23

	CoolPropDbl HelmholtzEOSMixtureBackend::calc_fugacity_coefficient(std::size_t i) {
	x_N_dependency_flag xN_flag = XN_DEPENDENT;
	if (isTwoPhase()) {
	// phi_i = f_i / (x_i * p). At VLE f_i^L == f_i^V but x_i^L != x_i^V, so
	// phi_i^L != phi_i^V and no convex combination is physically meaningful for
	// the overall two-phase state. Force callers to evaluate on SatL or SatV.
	throw ValueError(format("fugacity_coefficient is not well-defined in the two-phase region; evaluate on SatL or SatV instead"));
	} else if (isHomogeneousPhase()) {
	return exp(MixtureDerivatives::ln_fugacity_coefficient(*this, i, xN_flag));
	} else {
	throw ValueError(format("phase is invalid in calc_fugacity_coefficient"));
	}
	}
	CoolPropDbl HelmholtzEOSMixtureBackend::calc_fugacity(std::size_t i) {
	x_N_dependency_flag xN_flag = XN_DEPENDENT;
	CoolPropDbl localFug;
	if (isTwoPhase()) {
	if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for the two-phase properties"));
	if (std::abs(_Q) < DBL_EPSILON) {
	localFug = SatL->fugacity(i);
	} else if (std::abs(_Q - 1) < DBL_EPSILON) {
	localFug = SatV->fugacity(i);
	} else {
	localFug = _Q * SatV->fugacity(i) + (1 - _Q) * SatL->fugacity(i);
	}
	return static_cast<CoolPropDbl>(localFug);
	} else if (isHomogeneousPhase()) {
	return MixtureDerivatives::fugacity_i(*this, i, xN_flag);
	} else {
	throw ValueError(format("phase is invalid in calc_fugacity"));
	}
	}

The method was updated from 7.2 to 8.0, CoolProp 7.2 had the following definition:

CoolPropDbl HelmholtzEOSMixtureBackend::calc_fugacity_coefficient(std::size_t i) {
    x_N_dependency_flag xN_flag = XN_DEPENDENT;
    return exp(MixtureDerivatives::ln_fugacity_coefficient(*this, i, xN_flag));
}

It seems that the new definition has some safeguards against calculating fugacity coefficients in a mixture, but in this particular case, the vapour fraction was set to 0, so it can be considered a pure liquid state.

related: CoolProp/CoolProp.jl#49

[ibell]

Looks like a bug to me.

Claude: put in the same checks for saturation _Q values and return values for the SatV or SatL like is done for calc_fugacity and error out if not saturated