Merge pull request #47 from ypaul21/saftgammamie2

Speeding up SAFT Gamma Mie (and SAFT in general)

Project.toml

@@ -1,7 +1,7 @@
 name = "Clapeyron"
 uuid = "7c7805af-46cc-48c9-995b-ed0ed2dc909a"
 authors = ["Hon Wa Yew <yewhonwa@gmail.com>", "Pierre Walker <pwalker@mit.edu>"]
-version = "0.2.3"
+version = "0.2.4"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -11,6 +11,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 NLSolvers = "337daf1e-9722-11e9-073e-8b9effe078ba"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+PackedVectorsOfVectors = "7713531c-48ef-4bdd-9821-9ff7a8736089"
 PositiveFactorizations = "85a6dd25-e78a-55b7-8502-1745935b8125"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -19,12 +20,14 @@ ThermoState = "e7b6519d-fdf7-4a33-b544-2b37a9c1234a"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
+
 CSV = "0.9"
 DiffResults = "^1.0"
 ForwardDiff = "^0.10"
 LogExpFunctions = "^0.2,^0.3"
 NLSolvers = "^0.1"
 NaNMath = "^0.3"
+PackedVectorsOfVectors = "^0.1.2"
 PositiveFactorizations = "0.2"
 Roots = "^1"
 StaticArrays = "^1"

database/LatticeFluid/SanchezLacombe/SanchezLacombe_like.csv

@@ -0,0 +1,17 @@
+Clapeyron Database File
+SanchezLacombe Like Parameters
+species,epsilon,vol,segment
+carbon dioxide,2276.66,3.638,8.564
+ethane,2444.62,7.865,6.653
+ethylene,2273.34,7.238,6.518
+Xenon,2943.30,11.854,3.644
+Biphenyl,5280.12,14.985,11.708
+Naphthalene,5509.75,12.407,9.600 
+"2,3-Dimethylnaphthalene",5604.69,13.069,11.245
+"2,6-Dimethylnaphthalene",5540.04,13.069,11.244
+Fluorene,5866.55,10.676,11.821
+Phenanthrene,6354.58,13.804,11.378
+Anthracene,6372.58,13.804,11.378
+Pyrene,6732.91,14.733,11.980
+n-Octacosane,5295.78,17.360,22.989
+Benzoic Acid,4920.81,6.196,15.431

database/LatticeFluid/SanchezLacombe/mixing/k0k1l_unlike.csv

@@ -0,0 +1,23 @@
+Clapeyron Database File
+SanchezLacombe k0-k1-l unlike Parameters
+species1,species2,k0,k1,l
+carbon dioxide,biphenyl,0.1840,0.1580,0.0067
+carbon dioxide,naphthalene,0.1985,0.2509,-0.2162
+carbon dioxide,"2,3-m-naphthalene",0.2259,0.0907,-0.2729
+carbon dioxide,"2,6-m-naphthalene",0.2285,0.1858,-0.3110
+carbon dioxide,"fluorene",0.2279,0.1722,-0.4510
+carbon dioxide,"phenanthrene",0.2434,0.0959,-0.2301
+carbon dioxide,"anthracene",0.2471,0.2755,-0.1589
+carbon dioxide,"pyrene",0.2615,0.1518,-0.2040
+carbon dioxide,n-octacosane,0.2711,0.0465,-0.3562
+carbon dioxide,"benzoic acid",0.1191,0.2753,-0.4264
+ethane,biphenyl,0.0740,0.1496,-0.1069
+ethane,naphthalene,0.0678,0.2893,-0.2188
+ethane,anthracene,0.1197,0.2430,-0.1998
+ethane,n-octacosane,0.1249,0.0101,-0.4677
+ethane,"benzoic acid",-0.1455,0.1563,0.1817
+ethylene,fluorene,0.0759,0.1801,-0.3545
+ethylene,phenanthrene,0.2020,-0.1635,-0.0330
+ethylene,pyrene,0.1457,0.1632,-0.3614
+ethylene,benzoic acid,0.0340,0.2542,-0.4361
+Xenon,naphthalene,0.0283,-0.0214,-0.0372

database/properties/molarmass_groups.csv

@@ -1,4 +1,4 @@
-Clapeyron Database File,,,,,,,,,,,,,,,,
+Clapeyron Database File
 GC Single Data
 species,N_a,N_h,Mw
 -CH3,4,3,15.035
@@ -41,4 +41,40 @@ O=CH- (aldehyde),3,1,29.018
 -NO2,3,0,46.005
 -SH,2,1,33.07
 -S- (non-ring),1,0,32.06
--S- (ring),1,0,32.06
+-S- (ring),1,0,32.06
+CH3,4,3,15.035
+CH2,3,2,14.027
+CH,2,1,13.019
+C,1,0,12.011
+aCH,2,1,13.019
+aCCH2,4,2,26.03728
+aCCH,3,1,25.02934
+CH2=,3,2,14.027
+CH=,3,2,14.027
+cCH2,3,2,14.027
+COOH,4,1,45.017
+CH3COCH3,10,6,58.08
+COO,3,0,44.009
+H2O,3,1,18.02
+CH3OH,6,4,32.04
+CH4,5,1,16.04
+CO2,3,0,44.01
+OH,2,1,17.008
+CH2OH,5,3,31.03
+CHOH,4,2,30.03
+NH2,3,2,16.024
+NH,2,1,15.015
+N,1,0,14.007
+cNH,2,1,15.015
+cN,1,0,14.007
+C=,1,0,12.011
+aCCH3,5,3,27.04522
+aCOH,3,1,29.018
+cCH,2,1,13.019
+cCHNH,4,2,28.03328
+cCHN,3,1,27.02534
+aCCOaC,4,0,52.0315
+aCCOOH,5,1,57.02814
+aCNHaC,4,1,39.03604
+CH3CO,6,3,43.04462
+[CH2][OCH2],7,4,44.05256

src/Clapeyron.jl

@@ -1,31 +1,36 @@
 module Clapeyron
 using StaticArrays
 using LinearAlgebra
-using NLSolvers,Roots
+import PackedVectorsOfVectors
+const PackedVofV = PackedVectorsOfVectors.PackedVectorOfVectors
+using Roots: Roots
+using NLSolvers
 using DiffResults, ForwardDiff
 include("solvers/Solvers.jl")
 using .Solvers
 using .Solvers: log
 using Unitful
 import LogExpFunctions
 include("constants.jl")
-include("models/basetools.jl")
-include("utils/ParamOptions.jl")
+include("models/basetools.jl") #type hierarchy
+include("utils/ParamOptions.jl") 
+include("utils/vectors.jl")
 include("utils/ClapeyronParam.jl")
 
-include("models/eos/ideal/BasicIdeal.jl") #before macros, because its used there
 
 include("utils/macros.jl")
 using CSV, Tables
 include("utils/database.jl")
 include("utils/misc.jl")
 
+
 include("models/combiningrules.jl")
 
 include("models/eos.jl")
 include("utils/visualisation.jl")
 include("utils/split_model.jl")
 
+include("models/eos/ideal/BasicIdeal.jl") #before macros, because its used there
 include("models/eos/ideal/MonomerIdeal.jl")
 include("models/eos/ideal/ReidIdeal.jl")
 include("models/eos/ideal/WalkerIdeal.jl")
@@ -42,10 +47,13 @@ include("models/eos/SAFT/softSAFT/softSAFT.jl")
 include("models/eos/SAFT/SAFTVRMie/SAFTVRMie.jl")
 include("models/eos/SAFT/SAFTVRMie/variants/SAFTVRQMie.jl")
 include("models/eos/SAFT/SAFTgammaMie/SAFTgammaMie.jl")
+
 include("models/eos/SAFT/CKSAFT/CKSAFT.jl")
 include("models/eos/SAFT/CKSAFT/variants/sCKSAFT.jl")
 include("models/eos/SAFT/BACKSAFT/BACKSAFT.jl")
 
+include("models/eos/SAFT/equations.jl")
+
 include("models/eos/cubic/vdW.jl")
 include("models/eos/cubic/RK/RK.jl")
 include("models/eos/cubic/PR/PR.jl")
@@ -77,6 +85,8 @@ include("models/eos/EmpiricHelmholtz/IAPWS95.jl")
 include("models/eos/EmpiricHelmholtz/PropaneRef.jl")
 include("models/eos/EmpiricHelmholtz/GERG2008/GERG2008.jl")
 
+include("models/eos/LatticeFluid/SanchezLacombe/SanchezLacombe.jl")
+
 include("models/eos/SPUNG/SPUNG.jl")
 
 include("models/eos/cached/CachedEoS.jl")

src/methods/initial_guess.jl

@@ -264,8 +264,7 @@ function x0_sat_pure(model,T,z=SA[1.0])
     end
     p = -0.25*R̄*T/B
     vl = x0_volume(model,p,T,z,phase=:l)
-
-    vl = volume_compress(model,p,T,z,V0=vl)
+    vl = volume_compress(model,p,T,SA[1.0],V0=vl)
     #=the basis is that p = RT/v-b - a/v2
     we have a (p,v,T) pair
     and B = 2nd virial coefficient = b-a/RT
@@ -286,10 +285,12 @@ function x0_sat_pure(model,T,z=SA[1.0])
     _b = γ - B - vl
     Δ = _b*_b - 4*_c
     if isnan(vl) | (Δ < 0)
+        #fails on two ocassions:
+        #near critical point, or too low.
         #old strategy
         x0l = 4*lb_v
         x0v = -2*B + 2*lb_v
-        return [log10(x0l),log10(x0v)]   
+        return [log10(x0l),log10(x0v)]
     end
     Δsqrt = sqrt(Δ)
     b1 = 0.5*(-_b + Δsqrt)
@@ -323,7 +324,14 @@ function x0_sat_pure(model,T,z=SA[1.0])
         x0v = -2*B #gas volume as high as possible
         return [log10(x0l),log10(x0v)]   
     end
-
+    Vl0,Vv0 = vdw_x0_xat_pure(T,Tc,Pc,Vc)
+    x0l = min(Vl0,vl)
+    x0v = min(1e4*one(Vv0),Vv0) #cutoff volume
+    return [log10(x0l),log10(x0v)] 
+end
+
+function vdw_x0_xat_pure(T,T_c,P_c,V_c)
+    Tr = T/T_c
     Trm1 = 1.0-Tr
     Trmid = sqrt(Trm1)
     if Tr >= 0.7 
@@ -343,25 +351,19 @@ function x0_sat_pure(model,T,z=SA[1.0])
         #Eq. 31 valid in 0.25 < Tr < 1
         mean_c = 1.0 + 0.4*Trm1 + 0.161*Trm1*Trm1     
         c_v = 2*mean_c - c_l
-
     end
     #volumes predicted by vdW
-    Vl0 = (1/c_l)*Vc
-    Vv0 = (1/c_v)*Vc 
-    x0l = min(Vl0,vl)
-    x0v = min(1e4*one(Vv0),Vv0) #cutoff volume
-    return [log10(x0l),log10(x0v)] 
+    Vl0 = (1/c_l)*V_c
+    Vv0 = (1/c_v)*V_c
+    return (Vl0,Vv0)
 end
-
 function scale_sat_pure(model::EoSModel,z=SA[1.0])
     p    = 1/p_scale(model,z)
     μ    = 1/R̄/T_scale(model,z)
     return p,μ
 end
 
-
 #=x0_crit_pure=#
-
 # x0_crit_pure(model::LJSAFT) = [1.5, log10(π/6*model.params.segment[model.components[1]]*model.params.b[model.components[1]]/0.3)]
 
 """
@@ -444,6 +446,11 @@ function T_scale(model::SAFTModel,z=SA[1.0])
     return prod(ϵ)^(1/length(ϵ))
 end
 
+function T_scale(model::SAFTgammaMieModel,z=SA[1.0])
+    return T_scale(model.vrmodel,z)
+end
+
+
 function T_scale(model::LJSAFTModel,z=SA[1.0])
     T̃ = model.params.T_tilde.diagvalues
     return prod(T̃)^(1/length(T̃))
@@ -514,16 +521,9 @@ function p_scale(model::LJSAFTModel,z=SA[1.0])
 end
 
 function p_scale(model::SAFTgammaMieModel,z=SA[1.0])
-    vk  = model.groups.n_flattenedgroups
-    seg = model.params.segment.values
-    S   = model.params.shapefactor.values
-    σ   = model.params.sigma.values
-    m̄  = sum(vk[1][k]*seg[k]*S[k] for k in @groups(model.icomponents[1]))
-
-    σ̄3 = sum(vk[1][k]*vk[1][l]*
-                     seg[k]*seg[l]*
-                     S[k]*S[l]*
-                     σ[k,l]^3 for k in @groups(model.icomponents[1]) for l in @groups(model.icomponents[1]))/m̄^2
+    V = zero(first(z))
+    T = zero(first(z))
+    σ̄3 = @f(σ3x)
     ϵ̄ = T_scale(model,z)
     val    = σ̄3*N_A/R̄/ϵ̄
     return 1/val

src/methods/property_solvers/singlecomponent.jl

@@ -11,8 +11,13 @@ function try_sat_pure(model,V0,f!,T,result,error_val,method = LineSearch(Newton(
     end
 
     (P_sat,V_l,V_v) = result[]
+    _,dpdvl = p∂p∂V(model,V_l,T,SA[1.0])
+    _,dpdvv = p∂p∂V(model,V_v,T,SA[1.0])
+    (dpdvl > 0) | (dpdvv > 0) && return false
     ε = abs(V_l-V_v)/(eps(typeof(V_l-V_v)))
-    if (ε > 8)
+    ups = (eps(typeof(V_l-V_v)))
+    #@show abs(V_l-V_v),ups
+    if (ε > 5e9)
         #if ΔV > ε then Vl and Vv are different values
         return true
     end
@@ -47,6 +52,9 @@ function sat_pure(model::EoSModel, T, V0 = x0_sat_pure(model,T))
         return result[]
     end
     (T_c, p_c, V_c) = crit_pure(model)
+    if abs(T_c-T) < eps(typeof(T))
+        return (p_c,V_c,V_c)
+    end
     if T_c < T
         @error "initial temperature $T greater than critical temperature $T_c. returning NaN"
     else
@@ -60,34 +68,32 @@ function sat_pure(model::EoSModel, T, V0 = x0_sat_pure(model,T))
     return res0
 end
 
-#=
-    x0_sat_pure_crit(model,T,T_c,P_c,V_c)
 
-based on:
-DOI: 10.1007/s10910-007-9272-4
-Journal of Mathematical Chemistry, Vol. 43, No. 4, May 2008 (© 2007)
-The van der Waals equation: analytical and approximate solutions
+ #   x0_sat_pure_crit(model,T,T_c,P_c,V_c)
 
-we only use the upper part of the approximation
-but we use the calculated Tc, Vc and Pc. that shifts the 
-result enough to be useful as a starting point
-=#
-function x0_sat_pure_crit(model,T,T_c,P_c,V_c)
-    _1 = one(T)
-    Tr = T/T_c
-    Trm1 = _1-Tr
-    Trmid = sqrt(Trm1)
-    P_sat0 = (_1 -4.0*(Trm1)+4.8(Trm1*Trm1))*P_c
-    #c = 1/Vr
-    c_l = 1.0+2.0*Trmid + 0.4*Trm1 - 0.52*Trmid*Trm1 +0.115*Trm1*Trm1
-    c_v = 1.0-2.0*Trmid + 0.4*Trm1 + 0.52*Trmid*Trm1 +0.207*Trm1*Trm1
-    Vl0 = (1/c_l)*V_c
-    Vv0 = (1/c_v)*V_c
-    #Vl = volume_compress(model,P_sat0,T,V0=Vl0)
-    #Vv = volume_compress(model,P_sat0,T,V0=Vv0)
+
+
+#with the critical point, we can perform a
+#corresponding states approximation with the
+#propane reference equation of state
+function x0_sat_pure_crit(model,T,T_c,P_c,V_c) 
+    h = V_c*5000
+    T0 = 369.89*T/T_c
+    Vl0 = (1.0/_propaneref_rholsat(T0))*h
+    Vv0 = (1.0/_propaneref_rhovsat(T0))*h
+    _1 = SA[1.0]
+    #μ_l = only(VT_chemical_potential(model,Vl0,T,_1))
+    #μ_v = only(VT_chemical_potential(model,Vv0,T,_1))
+    #@show (μ_l < μ_v,T/T_c)
+    #if μ_l < μ_v 
+      #@show μ_l,μ_v
+    #end
+    # _,dpdvv = p∂p∂V(model,Vv0,T,SA[1.0])
+    # @show dpdvv*Vv0
+    # _,dpdvv = p∂p∂V(model,2*Vv0,T,SA[1.0])
+    # @show dpdvv*Vv0
     return [log10(Vl0),log10(Vv0)]
 end
-
 function sat_pure(model::EoSModel,V0,f!,T,method =LineSearch(Newton()) )  
     r = Solvers.nlsolve(f!, V0 ,method )
     Vsol = Solvers.x_sol(r)