improve stability of crit_pure (#501)

* add modification to crit_pure

* support multicomponent BACKSAFT

* change serialized model to eos_repr created model

* improve x0_sat_pure_virial when b constant is less than lb_volume

Author: longemen3000

database/SAFT/BACKSAFT/BACKSAFT_like.csv

@@ -1,47 +1,47 @@
 Clapeyron Database File
 BACKSAFT Like Parameters
-species,DIPPR Number,Mw,segment,c,vol,epsilon,alpha,source
-argon,,16,1.0,1.0,12.08,149.2,1.0,
-nitrogen,,30,1.0,1.0,14.05,128.9,1.031,
-oxygen,,44,1.0,1.0,11.8,155.9,1.017,
-carbon dioxide,,44,1.348,10.0,10.95,324.3,1.067,
-methane,,44,1.0,1.0,15.93,191.2,1.012,
-ethane,,44,1.07,10.0,21.43,308.4,1.028,
-propane~|~n-propane,,44,1.24,10.0,25.46,387.4,1.051,
-butane~|~n-butane,,44,1.34,10.0,29.88,457.9,1.071,
-pentane~|~n-pentane,,44,1.446,10.0,33.8,517.6,1.083,
-hexane~|~n-hexane,,44,1.573,10.0,36.64,570.7,1.097,
-heptane~|~n-heptane,,44,1.682,10.0,39.65,617.8,1.108,
-octane~|~n-octane,,44,1.814,10.0,42.05,659.7,1.117,
-nonane~|~n-nonane,,44,1.905,10.0,45.06,698.0,1.123,
-decane~|~n-decane,,44,2.018,10.0,47.37,732.0,1.13,
-dodecane~|~n-dodecane,,44,2.255,10.0,55.1,787.5,1.137,
-tetradecane~|~n-tetradecane,,44,2.434,10.0,60.0,841.4,1.148,
-hexadecane~|~n-hexadecane,,44,2.524,10.0,65.61,893.9,1.157,
-octadecane~|~n-octadecane,,44,2.786,10.0,71.14,926.5,1.162,
-eicosane~|~n-eicosane,,44,3.039,10.0,84.15,950.9,1.163,
-isobutane~|~2-methylpropane,,44,1.319,10.0,30.69,435.8,1.066,
-isopentane~|~methylbutane~|~2-methylbutane,,44,1.411,10.0,34.33,501.0,1.073,
-"neopentane~|~2,2-dimethylpropane",,44,1.353,10.0,36.55,463.9,1.06,
-isohexane~|~2-methylpentane,,44,1.541,10.0,37.41,553.3,1.088,
-"2,2-dimethylbutane",,44,1.433,10.0,39.68,533.3,1.074,
-"2,3-dimethylbutane",,44,1.474,10.0,38.29,550.3,1.081,
-3-methylhexane~|~3-methyl-hexane,,44,1.64,10.0,39.68,606.0,1.1,
-3-ethylpentane,,44,1.623,10.0,40.02,609.0,1.096,
-"2,2-dimethylpentane",,44,1.563,10.0,41.93,581.1,1.09,
-"2,2,3-trimethylbutane~|~2,2,3-trimethyl-butane",,44,1.488,10.0,43.04,585.0,1.08,
-3-methylheptane,,44,1.742,10.0,43.01,647.4,1.11,
-"3,3-dimethylhexane",,44,1.645,10.0,44.35,635.5,1.096,
-2-methyl-3-ethylpentane~|~3-ethyl-2-methylpentane,,44,1.661,10.0,43.73,642.4,1.097,
-"2,2,3-trimethylpentane",,44,1.6,10.0,45.4,630.5,1.086,
-"2,3,3-trimethylpentane",,44,1.575,10.0,45.92,640.6,1.087,
-cyclopropane,,44,1.171,10.0,21.99,416.0,1.052,
-ethylcyclopentane,,44,1.55,10.0,38.48,632.4,1.085,
-ethene~|~ethylene~|~1-ethylene,,44,1.039,10.0,20.03,283.4,1.034,
-propene~|~propylene,,44,1.202,10.0,24.22,379.5,1.047,
-pentene~|~1-pentene,,44,1.424,10.0,32.21,510.3,1.082,
-ethine~|~ethyne~|~acetylene,,44,1.254,10.0,14.18,326.6,1.064,
-propyne~|~1-propyne~|~methylacetylene,,44,1.355,10.0,18.63,435.5,1.076,
-benzene,,44,1.403,10.0,28.7,615.8,1.078,
-toluene,,44,1.503,10.0,33.34,655.9,1.076,
-"pxylene~|~p-xylene~|~1,4-dimethylbenzene",,44,1.576,10.0,37.15,703.4,1.101,
+species,segment,c,vol,epsilon,alpha,source
+argon,1.0,1.0,12.08,149.2,1.0,10.1016/S0378-3812(01)00521-0
+nitrogen,1.0,1.0,14.05,128.9,1.031,10.1016/S0378-3812(01)00521-0
+oxygen,1.0,1.0,11.8,155.9,1.017,10.1016/S0378-3812(01)00521-0
+carbon dioxide,1.348,10.0,10.95,324.3,1.067,10.1016/S0378-3812(01)00521-0
+methane,1.0,1.0,15.93,191.2,1.012,10.1016/S0378-3812(01)00521-0
+ethane,1.07,10.0,21.43,308.4,1.028,10.1016/S0378-3812(01)00521-0
+propane~|~n-propane,1.24,10.0,25.46,387.4,1.051,10.1016/S0378-3812(01)00521-0
+butane~|~n-butane,1.34,10.0,29.88,457.9,1.071,10.1016/S0378-3812(01)00521-0
+pentane~|~n-pentane,1.446,10.0,33.8,517.6,1.083,10.1016/S0378-3812(01)00521-0
+hexane~|~n-hexane,1.573,10.0,36.64,570.7,1.097,10.1016/S0378-3812(01)00521-0
+heptane~|~n-heptane,1.682,10.0,39.65,617.8,1.108,10.1016/S0378-3812(01)00521-0
+octane~|~n-octane,1.814,10.0,42.05,659.7,1.117,10.1016/S0378-3812(01)00521-0
+nonane~|~n-nonane,1.905,10.0,45.06,698.0,1.123,10.1016/S0378-3812(01)00521-0
+decane~|~n-decane,2.018,10.0,47.37,732.0,1.13,10.1016/S0378-3812(01)00521-0
+dodecane~|~n-dodecane,2.255,10.0,55.1,787.5,1.137,10.1016/S0378-3812(01)00521-0
+tetradecane~|~n-tetradecane,2.434,10.0,60.0,841.4,1.148,10.1016/S0378-3812(01)00521-0
+hexadecane~|~n-hexadecane,2.524,10.0,65.61,893.9,1.157,10.1016/S0378-3812(01)00521-0
+octadecane~|~n-octadecane,2.786,10.0,71.14,926.5,1.162,10.1016/S0378-3812(01)00521-0
+eicosane~|~n-eicosane,3.039,10.0,84.15,950.9,1.163,10.1016/S0378-3812(01)00521-0
+isobutane~|~2-methylpropane,1.319,10.0,30.69,435.8,1.066,10.1016/S0378-3812(01)00521-0
+isopentane~|~methylbutane~|~2-methylbutane,1.411,10.0,34.33,501.0,1.073,10.1016/S0378-3812(01)00521-0
+"neopentane~|~2,2-dimethylpropane",1.353,10.0,36.55,463.9,1.06,10.1016/S0378-3812(01)00521-0
+isohexane~|~2-methylpentane,1.541,10.0,37.41,553.3,1.088,10.1016/S0378-3812(01)00521-0
+"2,2-dimethylbutane",1.433,10.0,39.68,533.3,1.074,10.1016/S0378-3812(01)00521-0
+"2,3-dimethylbutane",1.474,10.0,38.29,550.3,1.081,10.1016/S0378-3812(01)00521-0
+3-methylhexane~|~3-methyl-hexane,1.64,10.0,39.68,606.0,1.1,10.1016/S0378-3812(01)00521-0
+3-ethylpentane,1.623,10.0,40.02,609.0,1.096,10.1016/S0378-3812(01)00521-0
+"2,2-dimethylpentane",1.563,10.0,41.93,581.1,1.09,10.1016/S0378-3812(01)00521-0
+"2,2,3-trimethylbutane~|~2,2,3-trimethyl-butane",1.488,10.0,43.04,585.0,1.08,10.1016/S0378-3812(01)00521-0
+3-methylheptane,1.742,10.0,43.01,647.4,1.11,10.1016/S0378-3812(01)00521-0
+"3,3-dimethylhexane",1.645,10.0,44.35,635.5,1.096,10.1016/S0378-3812(01)00521-0
+2-methyl-3-ethylpentane~|~3-ethyl-2-methylpentane,1.661,10.0,43.73,642.4,1.097,10.1016/S0378-3812(01)00521-0
+"2,2,3-trimethylpentane",1.6,10.0,45.4,630.5,1.086,10.1016/S0378-3812(01)00521-0
+"2,3,3-trimethylpentane",1.575,10.0,45.92,640.6,1.087,10.1016/S0378-3812(01)00521-0
+cyclopropane,1.171,10.0,21.99,416.0,1.052,10.1016/S0378-3812(01)00521-0
+ethylcyclopentane,1.55,10.0,38.48,632.4,1.085,10.1016/S0378-3812(01)00521-0
+ethene~|~ethylene~|~1-ethylene,1.039,10.0,20.03,283.4,1.034,10.1016/S0378-3812(01)00521-0
+propene~|~propylene,1.202,10.0,24.22,379.5,1.047,10.1016/S0378-3812(01)00521-0
+pentene~|~1-pentene,1.424,10.0,32.21,510.3,1.082,10.1016/S0378-3812(01)00521-0
+ethine~|~ethyne~|~acetylene,1.254,10.0,14.18,326.6,1.064,10.1016/S0378-3812(01)00521-0
+propyne~|~1-propyne~|~methylacetylene,1.355,10.0,18.63,435.5,1.076,10.1016/S0378-3812(01)00521-0
+benzene,1.403,10.0,28.7,615.8,1.078,10.1016/S0378-3812(01)00521-0
+toluene,1.503,10.0,33.34,655.9,1.076,10.1016/S0378-3812(01)00521-0
+"pxylene~|~p-xylene~|~1,4-dimethylbenzene",1.576,10.0,37.15,703.4,1.101,10.1016/S0378-3812(01)00521-0

docs/src/api/parameters.md

@@ -94,10 +94,8 @@ Clapeyron.xlogx
 ## Model Exporting
 
 ```@docs
+Clapeyron.eos_repr
 Clapeyron.export_model
-Clapeyron.FluidCorrelation
-Clapeyron.GammaPhi
-Clapeyron.ECS
 ```
 
 ## Model Citing

src/base/eosshow.jl

@@ -115,4 +115,79 @@ function show_reference_state(io::IO,ref,model::EoSModel,space)
     end
 end
 
-export eosshow
+"""
+    eos_repr(io::IO,model,newlines = true)
+    eos_repr(model;newlines = true)::String
+
+Given a `model::EoSModel`, returns a copy-pastable text representation of that model, designed to be parseable julia code.
+By default, `eos_repr` inserts some newlines for easier human reading of the output, that can be disabled via the `newlines` keyword argument.
+
+## Examples:
+
+```julia-repr
+julia> model = PCSAFT("methane")
+PCSAFT{BasicIdeal, Float64} with 1 component:
+ "methane"
+Contains parameters: Mw, segment, sigma, epsilon, epsilon_assoc, bondvol
+
+julia> push!(model.references,"test")
+3-element Vector{String}:
+ "10.1021/ie0003887"
+ "10.1021/ie010954d"
+ "test"
+
+julia> s = eos_repr(model);
+
+julia> model2 = eval(Meta.parse(s)) #exactly the same model
+PCSAFT{BasicIdeal, Float64} with 1 component:
+ "methane"
+Contains parameters: Mw, segment, sigma, epsilon, epsilon_assoc, bondvol
+
+julia> model2.references
+3-element Vector{String}:
+ "10.1021/ie0003887"
+ "10.1021/ie010954d"
+ "test"
+```
+"""
+function eos_repr(io::IO,model;inside = false,newlines = true)
+    M = typeof(model)
+    n = fieldnames(M)
+    newlines && print(io,"\n")
+    print(io,M)
+    print(io,"(")
+    if !inside && newlines
+        print(io,"\n")
+    end
+    k = 0
+    nf = length(n)
+    for i in n
+        k += 1
+        f = getfield(model,i)
+        F = typeof(f)
+        if F.name.module == Clapeyron || F isa EoSModel
+            eos_repr(io,f,inside = true)
+        elseif f isa PackedVofV
+            print(io,"Clapeyron.PackedVofV(")
+            show(io,f.p)
+            print(io,", ")
+            show(io,f.v)
+            print(io,")")
+        else
+            show(io,f)
+        end
+        k != nf && print(io,", ")
+    end
+    print(io,")")
+    
+
+    return nothing
+end
+
+function eos_repr(model;newlines = false)
+    io = IOBuffer()
+    eos_repr(io,model;newlines)
+    return String(take!(io))
+end
+
+export eosshow, eos_repr

src/methods/initial_guess.jl

@@ -267,6 +267,10 @@ function x0_sat_pure_virial(model,T)
     =#
 
     a,b = vdw_coefficients(vl,pl0,vv_virial,pv_eos,T)
+    if b < lb_v
+        b = lb_v
+        a = RT*(b - B)
+    end
     T̃ = RT*b/a
     T̃max = 0.2962962962962963 # Ωb/Ωa = 8/27
     T̃min = 0.1*T̃max
@@ -381,8 +385,13 @@ function x0_sat_pure_near0(model, T, vl0 = volume(model,zero(T),T,phase = :l);B
     ares = a_res(model, vl0, T, z)
     lnϕ_liq0 = ares - 1 + log(RT/vl0)
     p = exp(lnϕ_liq0)
-    vv = volume_virial(B,p,T)
     pB = -0.25*RT/B
+    if lnϕ_liq0 < log(eps(eltype(T)))
+        p = oneunit(p)*0.5*pB
+        vl0 = volume(model,p,T,z,vol0 = vl0,phase = :l)
+    end
+    vv = volume_virial(B,p,T)
+    
     if refine_vl && pB/p > 10
         vl = volume(model,p,T,z,vol0 = vl0,phase = :l)
     else

src/methods/property_solvers/singlecomponent/crit_pure.jl

@@ -19,44 +19,67 @@ function crit_pure(model::EoSModel)
     end
 end
 
+
+function crit_x_to_v(lbv,x)
+    lo = 0.001
+    hi = 1.0
+    η = lo + (hi - lo)/(exp(-x) + 1)
+    return lbv/η
+end
+
+function crit_v_to_x(lbv,v)
+    η = lbv/v
+    lo = 0.001
+    hi = 1.0
+    return -log((hi - lo)/(η - lo) - 1)
+end
+
 function crit_pure(model::EoSModel,x0,z = SA[1.0];options = NEqOptions())
     check_arraysize(model,z)
     #f! = (F,x) -> obj_crit(model, F, x[1]*T̄, exp10(x[2]))
+    zp = primalval(z)
+    primalmodel = primalval(model)
     if x0 === nothing
-        x0 = x0_crit_pure(model,z)
+        x0 = x0_crit_pure(primalmodel,zp)
     end
     x01,x02 = x0
-    T̄  = T_scale(model,z)*one(x01*one(x02))
-    #if type !== nothing
-    #    T̄ = T̄*oneunit(type)
-    #end
+    T̄  = T_scale(primalmodel,zp)*one(x01*one(x02))
+    Tc0 = x01*T̄
+    lbv0 = lb_volume(primalmodel,Tc0,zp)
     _1 = oneunit(T̄)
-    #x0 = SVector(_1*x01,_1*x02)
-    x0 = vec2(primalval(x01),primalval(x02*log(10)),primalval(_1))
-    primalmodel = primalval(model)
+    Tx0 = primalval(x01)
+    vc0 = exp10(primalval(x02))
+    x0 = vec2(Tx0,crit_v_to_x(lbv0,vc0),_1)
     zz = z/sum(z)
     f!(F,x) = ObjCritPure(primalmodel,F,primalval(T̄),x,zz)
     solver_res = Solvers.nlsolve(f!, x0, TrustRegion(Newton(), NLSolvers.NWI()), options)
-    #display(solver_res)
     r  = Solvers.x_sol(solver_res)
+    !all(<(solver_res.options.f_abstol),solver_res.info.best_residual) && (r .= NaN)
     T_c = r[1]*T̄
-    V_c = exp(r[2])
+    lbv = lb_volume(primalmodel,T_c,z)
+    V_c = crit_x_to_v(lbv,r[2])
     p_c = pressure(model, V_c, T_c, zz)
+    if p_c < 0
+        p_c *= NaN
+        V_c *= NaN
+        T_c *= NaN
+    end
     crit = (T_c, p_c, V_c)
     return crit_pure_ad(model,crit,z)
 end
 
 function crit_pure_ad(model,crit,z)
-    if has_dual(model)
+    if has_dual(model) || has_dual(z)
         T_c_primal, p_c_primal, V_c_primal = crit
         T̄  = T_scale(model)
-        x = SVector(T_c_primal/T̄,log(V_c_primal))
-        f(zz) = __ObjCritPure(model,T̄,z,zz)
+        lbv = lb_volume(model,T̄,z)
+        x = SVector(T_c_primal/T̄,crit_v_to_x(lbv,V_c_primal))
+        f(zz) = __ObjCritPure(model,T̄,z,zz,lbv)
         F,J = Solvers.J2(f,x)
         ∂x = J\F
         r = x .- ∂x
         T_c = r[1]*T̄
-        V_c = exp(r[2])
+        V_c = crit_x_to_v(lbv,r[2])
         P_c = pressure(model,V_c,T_c,z)
         return (T_c,P_c,V_c)
     else
@@ -71,9 +94,10 @@ function ObjCritPure(model::T,F,T̄,x,z) where T
     return F
 end
 
-function __ObjCritPure(model::T,T̄,x,z) where T
+function __ObjCritPure(model::T,T̄,x,z,) where T
     T_c = x[1]*T̄
-    V_c = exp(x[2])
+    lbv = lb_volume(model,T_c,z)
+    V_c = crit_x_to_v(lbv,x[2])
     RT = Rgas(model)*T_c
     ∂²A∂V²_scale = V_c*V_c/RT
     ∂³A∂V³_scale = ∂²A∂V²_scale*V_c

src/methods/property_solvers/singlecomponent/saturation/ChemPotV.jl

@@ -10,9 +10,10 @@
                         atol = 1e-8,
                         rtol = 1e-12,
                         max_iters = 10^4)
-Default `saturation_pressure` Saturation method used by `Clapeyron.jl`. It uses equality of Chemical Potentials with a volume basis. If no volumes are provided, it will use  [`x0_sat_pure`](@ref).
-If those initial guesses fail and the specification is near critical point, it will try one more time, using Corresponding States instead.
-when `crit_retry` is true, if the initial solve fail, it will try to obtain a better estimate by calculating the critical point.
+Default `saturation_pressure` Saturation method used by `Clapeyron.jl`.
+It uses equality of Chemical Potentials with a volume basis.
+If no volumes are provided, it will use  [`x0_sat_pure`](@ref).
+When `crit_retry` is true, if the initial solve fail, it will try to obtain a better estimate by calculating the critical point.
 `f_limit`, `atol`, `rtol`, `max_iters` are passed to the non linear system solver.
 """
 struct ChemPotVSaturation{T,C} <: SaturationMethod