Evaluation of CPA at too high a density yields no error, just NaN

[ianhbell]

If you evaluate the CPA model for a density greater than $1/b$, $b$ being the covolume, the returned value for pressure is NaN because the term $\ln(1-b\rho)$ is NaN because $b\rho > 1$, but you should instead get an error that the density is above the maximum.

[pw0908]

Are you using the internal function pressure? We do partly need to be able to evaluate this function as, within our iterative solvers, they need an output of NaN to realise that they are exploring the unphysical region. Throwing an error would cause the solvers to fail the instant they step out of bounds (which one might expect to happen, especially when solving for the liquid root.

@longemen3000 thoughts?

[ianhbell]

I see your point, but it isn't obvious what is going on. Sometimes this is handled with an argument that says what to do when the argument is out of range, either to error or give wrong value.

[longemen3000]

To be fair, pressure is not an internal function, but it is a function that works in the same fundamental variables as the Helmholtz energy (V,T,n). this is not a problem specific to cubics. (another (old) related error is https://github.com/ClapeyronThermo/Clapeyron.jl/issues/24).

One design decision we took with Clapeyron some time ago was to make the eos function itself not error, but return NaN on failure. We even have our own definitions of log, sqrt and ^, because the default julia functions error:

julia> Clapeyron.log(-1)
NaN

julia> log(-1)
ERROR: DomainError with -1.0:
log was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(f::Symbol, x::Float64)
   @ Base.Math .\math.jl:33
 [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
   @ Base.Math .\special\log.jl:301
 [3] log
   @ .\special\log.jl:267 [inlined]
 [4] log(x::Int64)
   @ Base.Math .\math.jl:1578
 [5] top-level scope
   @ REPL[6]:1

One problem we encountered early on is that when errors of this style appeared, they were really deep in the call chain. (for example, one iteration of bubble_pressure could fall in invalid volume areas and throw errors). As we use external packages for (some) solvers, those solvers are capable of handling NaN data, but not errors in the function itself.

Coming back on the maximum density, some equations of state (like QCPR and Cubic + Chain) have a T-dependent lower bound volume, (we don't have an adequate api for those cases, as our lb_volume function does not accept temperature, that is an improvement that can be done). In particular, Cubic + Chain (and their modern version, Cubic + Chain + Assoc) are only defined on low reduced temperatures (the maximum density becomes negative on temperatures over 2-3 Tc).

Other equations of state are just empirically defined and can be evaluated outside the domain of maximum density (multiparameter EoS), but the value is still required (in our particular case, we use some empiric rules at construction time to define a sufficiently high density). There are some cases where the theoretical lower bound volume does not correspond to a valid fluid (or solid) volume, like softSAFT. this is due to the nature of the MBWR eos used for their LJ Helmholtz , if we replace that model for the Thol Helmholtz energy for LJ fluids (Clapeyron.softSAFT2016), The spurious zone vanishes.

Finally, if we consider that the maximum density is the maximum fluid density, then there are EoS that have a solid phase defined when rho > rho_max: https://github.com/ClapeyronThermo/Clapeyron.jl/blob/2a3d6eae261393bde743f57fd358d9ee5b955665/src/models/AnalyticalSLV/AnalyticalSLV.jl#L153-L179

On the checks themselves, for a V,T state, there are two invalid ranges, when rho > rho_max and when the values of density are inside the spinodal. The check on the maximum density, while not free, is cheap enough to not add significant penalty (< 1 eos evaluation). the spinodal check requires at least a dpdrho evaluation (and for multiparameter EoS that is not enough, as there may be invalid areas where dpdrho > 0). on multicomponent models, there is the additional stability checks that require additional eos evaluations at higher derivatives.

What we could do, is document better those pitfalls. We we are aware of those, but not everybody starting to use a thermodynamic library.

[ianhbell]

Yeah this comes down to a philosophical discussion about what to do about invalid values. I think that a high-level function like pressure should probably throw an error, while lower-level functions should be permitted to work in NaN. Or perhaps even better, accept the approach of rust of railroading - return a variant from the intermediate function that is either a correct value or error information. But as you get deep into call stacks, exception handling still remains a cleaner option in my opinion.

Also, you are right that the multicomponent models have big problems with bubbles of stability inside the spinodal. This is a big problem that has been covered in some depth by people like Oivind Wilhelmsen's group, but as of yet we don't have a great solution. It also doesn't help that some EOS don't have a well-formed spinodal, so you can't even use the spinodal as your determination in all cases.