diff --git a/database/properties/critical.csv b/database/properties/critical.csv index 161d0af5b..b8223e4a3 100644 --- a/database/properties/critical.csv +++ b/database/properties/critical.csv @@ -12,7 +12,7 @@ n-Heptane,142825,540.2,2.72e6,0.346 n-Octane,111659,568.7,2.47e6,0.396 n-Nonane,111842,594.6,2.31e6,0.446 n-Decane,124185,617.7,2.09e6,0.488 -n-Undecane,1120214,639,1.95e6,0.53 +undecane,1120214,639,1.95e6,0.53 n-Dodecane,112403,658,1.82e6,0.577 n-Tridecane,629505,675,1.68e6,0.617 n-Tetradecane,629594,693,1.57e6,0.643 diff --git a/docs/src/api/parameters.md b/docs/src/api/parameters.md index 61e786ae5..dc5c2036f 100644 --- a/docs/src/api/parameters.md +++ b/docs/src/api/parameters.md @@ -27,4 +27,5 @@ Clapeyron.SiteParam ## Model Splitting ```@docs Clapeyron.split_model +Clapeyron.is_splittable ``` diff --git a/src/Clapeyron.jl b/src/Clapeyron.jl index 13f7080e3..46c2c8bd1 100644 --- a/src/Clapeyron.jl +++ b/src/Clapeyron.jl @@ -48,11 +48,18 @@ include("models/eos/SAFT/CKSAFT/variants/sCKSAFT.jl") include("models/eos/SAFT/BACKSAFT/BACKSAFT.jl") -include("models/eos/cubic/alphas/alphas.jl") include("models/eos/cubic/vdW.jl") include("models/eos/cubic/RK/RK.jl") include("models/eos/cubic/PR/PR.jl") + +include("models/eos/cubic/alphas/alphas.jl") + +include("models/eos/cubic/mixing/mixing.jl") + +include("models/eos/cubic/RK/variants/SRK.jl") +include("models/eos/cubic/PR/variants/PR78.jl") + include("models/eos/cubic/equations.jl") include("models/eos/EmpiricHelmholtz/IAPWS95.jl") diff --git a/src/methods/differentials.jl b/src/methods/differentials.jl index dd8d7390d..9259f83c1 100644 --- a/src/methods/differentials.jl +++ b/src/methods/differentials.jl @@ -154,7 +154,7 @@ returns the second order volume (`V`) and temperature (`T`) derivatives of the t ``` [ ∂²f/∂V² ∂²f/∂V∂T - ∂²f/∂V∂T ∂²f/∂V²] + ∂²f/∂V∂T ∂²f/∂T²] ``` use this instead of the ∂2f if you only need second order information. ∂2f also gives zeroth and first order derivative information, but due to a bug in the used AD, it allocates more than necessary. diff --git a/src/models/eos/SPUNG/SPUNG.jl b/src/models/eos/SPUNG/SPUNG.jl index 1eb643b36..3731545d5 100644 --- a/src/models/eos/SPUNG/SPUNG.jl +++ b/src/models/eos/SPUNG/SPUNG.jl @@ -40,8 +40,8 @@ end function shape_factors(model::SPUNG{<:ABCubicModel},V,T,z=SA[1.0]) n = sum(z) x = z * (1/n) - a,b = cubic_ab(model.shape_model,T,x) - a0,b0 = cubic_ab(model.shape_ref,T,SA[1.0]) + a,b = cubic_ab(model.shape_model,V,T,x) + a0,b0 = cubic_ab(model.shape_ref,V,T,SA[1.0]) h = b/b0 fh = n*a/a0 f = fh/h diff --git a/src/models/eos/cubic/PR/PR.jl b/src/models/eos/cubic/PR/PR.jl index d8226fba7..34d9ca0b0 100644 --- a/src/models/eos/cubic/PR/PR.jl +++ b/src/models/eos/cubic/PR/PR.jl @@ -9,11 +9,11 @@ struct PRParam <: EoSParam Mw::SingleParam{Float64} end -struct PR{T <: IdealModel,α,γ} <: RKModel +struct PR{T <: IdealModel,α,γ} <:PRModel components::Array{String,1} icomponents::UnitRange{Int} alpha::α - activity::γ + mixing::γ params::PRParam idealmodel::T absolutetolerance::Float64 @@ -34,16 +34,16 @@ end Base.length(model::PR) = Base.length(model.icomponents) -molecular_weight(model::PR,z=SA[1.0]) = group_molecular_weight(model.groups,mw(model),z) +molecular_weight(model::PR,z=SA[1.0]) = comp_molecular_weight(mw(model),z) export PR function PR(components::Vector{String}; idealmodel=BasicIdeal, - alpha = SoaveAlpha, - activity = nothing, + alpha = PRAlpha, + mixing = vdW1fRule, userlocations=String[], ideal_userlocations=String[], alpha_userlocations = String[], - activity_userlocations = String[], + mixing_userlocations = String[], verbose=false) params = getparams(components, ["properties/critical.csv", "properties/molarmass.csv","SAFT/PCSAFT/PCSAFT_unlike.csv"]; userlocations=userlocations, verbose=verbose) k = params["k"] @@ -59,11 +59,11 @@ function PR(components::Vector{String}; idealmodel=BasicIdeal, init_idealmodel = init_model(idealmodel,components,ideal_userlocations,verbose) init_alpha = init_model(alpha,components,alpha_userlocations,verbose) - init_activity = init_model(activity,components,activity_userlocations,verbose) + init_mixing = init_model(mixing,components,mixing_userlocations,verbose) icomponents = 1:length(components) packagedparams = PRParam(a, b, params["Tc"],_pc,Mw) references = String[] - model = PR(components,icomponents,init_alpha,init_activity,packagedparams,init_idealmodel,1e-12,references) + model = PR(components,icomponents,init_alpha,init_mixing,packagedparams,init_idealmodel,1e-12,references) return model end @@ -72,24 +72,19 @@ function ab_consts(::Type{<:PRModel}) return 0.457235,0.077796 end -function cubic_ab(model::PR{<:Any,SoaveAlpha},T,z=SA[1.0],n=sum(z)) +function cubic_ab(model::PR{<:Any,<:Any,<:Any},V,T,z=SA[1.0],n=sum(z)) + invn2 = (one(n)/n)^2 a = model.params.a.values b = model.params.b.values - α = model.alpha - ω = α.params.acentricfactor.values - Tc = model.params.Tc.values - invn = one(n)/n - αx = @. (1+(0.37464+1.54226*ω-0.26992*ω^2)*(1-√(T/Tc))) * z * invn - āᾱ = dot(αx,Symmetric(a),αx) - b̄ = dot(z,Symmetric(b),z) * invn*invn - return āᾱ ,b̄ + α = @f(α_function,model.alpha) + ā,b̄ = @f(mixing_rule,model.mixing,α,a,b) + return ā ,b̄ end - function cubic_abp(model::PRModel, V, T, z) n = sum(z) v = V/n - āᾱ ,b̄ = cubic_ab(model,T,z,n) + āᾱ ,b̄ = cubic_ab(model,V,T,z,n) _1 = one(b̄) denom = evalpoly(v,(-b̄*b̄,2*b̄,_1)) p = R̄*T/(v-b̄) - āᾱ /denom @@ -97,7 +92,7 @@ function cubic_abp(model::PRModel, V, T, z) end function cubic_poly(model::PRModel,p,T,z) - a,b = cubic_ab(model,T,z) + a,b = cubic_ab(model,p,T,z) RT⁻¹ = 1/(R̄*T) A = a*p*RT⁻¹*RT⁻¹ B = b*p*RT⁻¹ @@ -105,6 +100,7 @@ function cubic_poly(model::PRModel,p,T,z) k₁ = -B*(3*B+2.0) + A k₂ = B-1.0 k₃ = 1.0 + return [k₀,k₁,k₂,k₃] end #= (-B2-2(B2+B)+A) @@ -113,7 +109,7 @@ end =# function a_res(model::PRModel, V, T, z) n = sum(z) - a,b = cubic_ab(model,T,z,n) + a,b = cubic_ab(model,V,T,z,n) Δ1 = 1+√2 Δ2 = 1-√2 ΔPRΔ = 2*√2 diff --git a/src/models/eos/cubic/PR/variants/PR78.jl b/src/models/eos/cubic/PR/variants/PR78.jl new file mode 100644 index 000000000..3cf5a8337 --- /dev/null +++ b/src/models/eos/cubic/PR/variants/PR78.jl @@ -0,0 +1,19 @@ +#just a function, no struct +function PR78(components::Vector{String}; idealmodel=BasicIdeal, + mixing = vdW1fRule, + userlocations=String[], + ideal_userlocations=String[], + alpha_userlocations = String[], + mixing_userlocations = String[], + verbose=false) + + return PR(components; + idealmodel = idealmodel, + alpha = PR78Alpha, + mixing=mixing, + ideal_userlocations = ideal_userlocations, + alpha_userlocations = alpha_userlocations, + mixing_userlocations = mixing_userlocations, + verbose = verbose) +end +export PR78 \ No newline at end of file diff --git a/src/models/eos/cubic/RK/RK.jl b/src/models/eos/cubic/RK/RK.jl index 6946ad618..c12b898b9 100644 --- a/src/models/eos/cubic/RK/RK.jl +++ b/src/models/eos/cubic/RK/RK.jl @@ -12,7 +12,7 @@ struct RK{T <: IdealModel,α,γ} <: RKModel components::Array{String,1} icomponents::UnitRange{Int} alpha::α - activity::γ + mixing::γ params::RKParam idealmodel::T absolutetolerance::Float64 @@ -37,12 +37,12 @@ molecular_weight(model::RK,z=SA[1.0]) = comp_molecular_weight(mw(model),z) export RK function RK(components::Vector{String}; idealmodel=BasicIdeal, - alpha = nothing, - activity = nothing, + alpha = RKAlpha, + mixing = vdW1fRule, userlocations=String[], ideal_userlocations=String[], alpha_userlocations = String[], - activity_userlocations = String[], + mixing_userlocations = String[], verbose=false) params = getparams(components, ["properties/critical.csv", "properties/molarmass.csv","SAFT/PCSAFT/PCSAFT_unlike.csv"]; userlocations=userlocations, verbose=verbose) k = params["k"] @@ -58,11 +58,11 @@ function RK(components::Vector{String}; idealmodel=BasicIdeal, init_idealmodel = init_model(idealmodel,components,ideal_userlocations,verbose) init_alpha = init_model(alpha,components,alpha_userlocations,verbose) - init_activity = init_model(activity,components,activity_userlocations,verbose) + init_mixing = init_model(mixing,components,mixing_userlocations,verbose) icomponents = 1:length(components) packagedparams = RKParam(a, b, params["Tc"],_pc,Mw) references = String[] - model = RK(components,icomponents,init_alpha,init_activity,packagedparams,init_idealmodel,1e-12,references) + model = RK(components,icomponents,init_alpha,init_mixing,packagedparams,init_idealmodel,1e-12,references) return model end @@ -71,20 +71,18 @@ function ab_consts(::Type{<:RKModel}) Ωb = (2^(1/3)-1)/3 return Ωa,Ωb end -function cubic_ab(model::RK{<:Any,Nothing},T,z=SA[1.0],n=sum(z)) +function cubic_ab(model::RK{<:Any,<:Any,<:Any},V,T,z=SA[1.0],n=sum(z)) invn2 = (one(n)/n)^2 a = model.params.a.values b = model.params.b.values - Tc = model.params.Tc.values - T̄c = sum(sqrt.(Tc*Tc')) - āᾱ = dot(z,Symmetric(a),z) * invn2*sqrt(T̄c/T) - b̄ = dot(z,Symmetric(b),z) * invn2 - return āᾱ ,b̄ + α = @f(α_function,model.alpha) + ā,b̄ = @f(mixing_rule,model.mixing,α,a,b) + return ā ,b̄ end function cubic_abp(model::RKModel, V, T, z) n = sum(z) - a,b = cubic_ab(model,T,z,n) + a,b = cubic_ab(model,V,T,z,n) v = V/n p = R̄*T/(v-b) - a/((v+b)*v) return a,b,p @@ -92,7 +90,7 @@ end function cubic_poly(model::RKModel,p,T,z) n = sum(z) - a,b = cubic_ab(model,T,z,n) + a,b = cubic_ab(model,p,T,z,n) RT⁻¹ = 1/(R̄*T) A = a*p* RT⁻¹* RT⁻¹ B = b*p* RT⁻¹ @@ -103,7 +101,7 @@ end function a_res(model::RKModel, V, T, z) n = sum(z) - ā,b̄ = cubic_ab(model,T,z,n) + ā,b̄ = cubic_ab(model,V,T,z,n) ρ = n/V RT⁻¹ = 1/(R̄*T) return -log(1-b̄*ρ) - ā*RT⁻¹*log(b̄*ρ+1)/b̄ @@ -112,4 +110,4 @@ end cubic_zc(::RKModel) = 1/3 -include("variants/SRK.jl") +# include("variants/SRK.jl") diff --git a/src/models/eos/cubic/RK/variants/SRK.jl b/src/models/eos/cubic/RK/variants/SRK.jl index 7b20d14a2..855cc21d9 100644 --- a/src/models/eos/cubic/RK/variants/SRK.jl +++ b/src/models/eos/cubic/RK/variants/SRK.jl @@ -1,32 +1,19 @@ -function cubic_ab(model::RK{<:Any,SoaveAlpha},T,z=SA[1.0],n=sum(z)) - a = model.params.a.values - b = model.params.b.values - αmodel = model.alpha - ω = αmodel.params.acentricfactor.values - Tc = model.params.Tc.values - _1 = one(T+n) - invn = one(n)/n - αx = @. (_1+(0.480+1.547*ω-0.176*ω^2)*(1-√(T/Tc))) * z * invn - āᾱ =dot(αx, Symmetric(a), αx) - b̄ = dot(z, Symmetric(b), z) * invn * invn - return āᾱ ,b̄ -end #just a function, no struct function SRK(components::Vector{String}; idealmodel=BasicIdeal, - activity = nothing, + mixing = vdW1fRule, userlocations=String[], ideal_userlocations=String[], alpha_userlocations = String[], - activity_userlocations = String[], + mixing_userlocations = String[], verbose=false) return RK(components; idealmodel = idealmodel, alpha = SoaveAlpha, - activity=activity, + mixing=mixing, ideal_userlocations = ideal_userlocations, alpha_userlocations = alpha_userlocations, - activity_userlocations = activity_userlocations, + mixing_userlocations = mixing_userlocations, verbose = verbose) end export SRK \ No newline at end of file diff --git a/src/models/eos/cubic/alphas/BM.jl b/src/models/eos/cubic/alphas/BM.jl new file mode 100644 index 000000000..5536c8509 --- /dev/null +++ b/src/models/eos/cubic/alphas/BM.jl @@ -0,0 +1,34 @@ +abstract type BMAlphaModel <: AlphaModel end + +struct BMAlphaParam <: EoSParam + acentricfactor::SingleParam{Float64} +end + +@newmodelsimple BMAlpha BMAlphaModel BMAlphaParam + +export BMAlpha +function BMAlpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false) + params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose) + acentricfactor = SingleParam(params["w"],"acentric factor") + packagedparams = BMAlphaParam(acentricfactor) + model = BMAlpha(packagedparams, verbose=verbose) + return model +end + +function α_function(model::RKModel,V,T,z,alpha_model::BMAlphaModel) + Tc = model.params.Tc.values + Tr = @. T/Tc + ω = alpha_model.params.acentricfactor.values + m = @. 0.480+1.547*ω-0.176*ω^2 + α = @. (Tr>1)*(exp((1-2/(2+m))*(1-Tr^(1+m/2))))^2+(Tr<=1)*(1+m*(1-√(Tr)))^2 + return α +end + +function α_function(model::PRModel,V,T,z,alpha_model::BMAlphaModel) + Tc = model.params.Tc.values + Tr = @. T/Tc + ω = alpha_model.params.acentricfactor.values + m = @. 0.37464+1.54226*ω-0.26992*ω^2 + α = @. (Tr>1)*(exp((1-2/(2+m))*(1-Tr^(1+m/2))))^2+(Tr<=1)*(1+m*(1-√(Tr)))^2 + return α +end \ No newline at end of file diff --git a/src/models/eos/cubic/alphas/PR78Alpha.jl b/src/models/eos/cubic/alphas/PR78Alpha.jl new file mode 100644 index 000000000..136eef201 --- /dev/null +++ b/src/models/eos/cubic/alphas/PR78Alpha.jl @@ -0,0 +1,24 @@ +abstract type PR78AlphaModel <: AlphaModel end + +struct PR78AlphaParam <: EoSParam + acentricfactor::SingleParam{Float64} +end + +@newmodelsimple PR78Alpha PR78AlphaModel PR78AlphaParam + +export PR78Alpha +function PR78Alpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false) + params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose) + acentricfactor = SingleParam(params["w"],"acentric factor") + packagedparams = PR78AlphaParam(acentricfactor) + model = PR78Alpha(packagedparams, verbose=verbose) + return model +end + +function α_function(model::CubicModel,V,T,z,alpha_model::PR78AlphaModel) + Tc = model.params.Tc.values + ω = alpha_model.params.acentricfactor.values + m = @. (ω<=0.491)*(0.37464+1.54226*ω-0.26992*ω^2)+(ω>0.491)*(0.379642+1.487503*ω-0.164423*ω^2-0.016666*ω^3) + α = @. (1+m*(1-√(T/Tc)))^2 + return α +end \ No newline at end of file diff --git a/src/models/eos/cubic/alphas/PRAlpha.jl b/src/models/eos/cubic/alphas/PRAlpha.jl new file mode 100644 index 000000000..b633826e0 --- /dev/null +++ b/src/models/eos/cubic/alphas/PRAlpha.jl @@ -0,0 +1,23 @@ +abstract type PRAlphaModel <: AlphaModel end + +struct PRAlphaParam <: EoSParam + acentricfactor::SingleParam{Float64} +end + +@newmodelsimple PRAlpha PRAlphaModel PRAlphaParam + +export PRAlpha +function PRAlpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false) + params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose) + acentricfactor = SingleParam(params["w"],"acentric factor") + packagedparams = PRAlphaParam(acentricfactor) + model = PRAlpha(packagedparams, verbose=verbose) + return model +end + +function α_function(model::CubicModel,V,T,z,alpha_model::PRAlphaModel) + Tc = model.params.Tc.values + ω = alpha_model.params.acentricfactor.values + α = @. (1+(0.37464+1.54226*ω-0.26992*ω^2)*(1-√(T/Tc)))^2 + return α +end \ No newline at end of file diff --git a/src/models/eos/cubic/alphas/RKAlpha.jl b/src/models/eos/cubic/alphas/RKAlpha.jl new file mode 100644 index 000000000..739ac9e49 --- /dev/null +++ b/src/models/eos/cubic/alphas/RKAlpha.jl @@ -0,0 +1,22 @@ +abstract type RKAlphaModel <: AlphaModel end + +struct RKAlphaParam <: EoSParam +end + + +@newmodelsimple RKAlpha RKAlphaModel RKAlphaParam +is_splittable(::RKAlpha) = false + +export RKAlpha +function RKAlpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false) + params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose) + packagedparams = RKAlphaParam() + model = RKAlpha(packagedparams, verbose=verbose) + return model +end + +function α_function(model::CubicModel,V,T,z,alpha_model::RKAlphaModel) + Tc = model.params.Tc.values + α = @. 1 /√(T/Tc) + return α +end \ No newline at end of file diff --git a/src/models/eos/cubic/alphas/alphas.jl b/src/models/eos/cubic/alphas/alphas.jl index de4d16c37..f3fefe58b 100644 --- a/src/models/eos/cubic/alphas/alphas.jl +++ b/src/models/eos/cubic/alphas/alphas.jl @@ -4,9 +4,15 @@ end function init_model(model::Type{<:AlphaModel},components,userlocations,verbose) verbose && @info("""Now creating alpha model: - $idealmodel""") + $model""") return model(components;userlocations,verbose) end +has_sites(::Type{<:AlphaModel})=false +has_groups(::Type{<:AlphaModel})=false -include("soave.jl") \ No newline at end of file +include("RKAlpha.jl") +include("PRAlpha.jl") +include("PR78Alpha.jl") +include("soave.jl") +include("BM.jl") \ No newline at end of file diff --git a/src/models/eos/cubic/alphas/soave.jl b/src/models/eos/cubic/alphas/soave.jl index 1313644bf..5df33b8b5 100644 --- a/src/models/eos/cubic/alphas/soave.jl +++ b/src/models/eos/cubic/alphas/soave.jl @@ -4,7 +4,6 @@ struct SoaveAlphaParam <: EoSParam acentricfactor::SingleParam{Float64} end -abstract type PRModel <: ABCubicModel end @newmodelsimple SoaveAlpha SoaveAlphaModel SoaveAlphaParam export SoaveAlpha @@ -16,4 +15,9 @@ function SoaveAlpha(components::Vector{String}; userlocations::Vector{String}=St return model end - \ No newline at end of file +function α_function(model::CubicModel,V,T,z,alpha_model::SoaveAlpha) + Tc = model.params.Tc.values + ω = alpha_model.params.acentricfactor.values + α = @. (1+(0.480+1.547*ω-0.176*ω^2)*(1-√(T/Tc)))^2 + return α +end \ No newline at end of file diff --git a/src/models/eos/cubic/mixing/Kay.jl b/src/models/eos/cubic/mixing/Kay.jl new file mode 100644 index 000000000..e706650ae --- /dev/null +++ b/src/models/eos/cubic/mixing/Kay.jl @@ -0,0 +1,26 @@ +abstract type KayRuleModel <: MixingRule end + +struct KayRuleParam <: EoSParam +end + +@newmodelsimple KayRule KayRuleModel KayRuleParam + +export KayRule + +function KayRule(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false) + #params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose) + #acentricfactor = SingleParam(params["w"],"acentric factor") + packagedparams = KayRuleParam() + model = KayRule(packagedparams, verbose=verbose) + return model +end + +function mixing_rule(model::ABCubicModel,V,T,z,mixing_model::KayRuleModel,α,a,b) + n = sum(z) + invn2 = (one(n)/n)^2 + b̄ = (dot(z,Symmetric(b.^(1/3)),z) * invn2).^3 + ā = √(dot(z,Symmetric(a .* sqrt.(α*α') ./ b).^2,z)* invn2) * b̄ + return ā,b̄ +end + +is_splittable(::KayRule) = false diff --git a/src/models/eos/cubic/mixing/mixing.jl b/src/models/eos/cubic/mixing/mixing.jl new file mode 100644 index 000000000..ee39fcbfc --- /dev/null +++ b/src/models/eos/cubic/mixing/mixing.jl @@ -0,0 +1,15 @@ +abstract type MixingRule <:EoSModel end +function init_model(model::MixingRule,components,userlocations,verbose) + return model +end + +function init_model(model::Type{<:MixingRule},components,userlocations,verbose) + verbose && @info("""Now creating mixing model: + $model""") + return model(components;userlocations,verbose) +end + +has_sites(::Type{<:MixingRule})=false + +include("vdW1f.jl") +include("Kay.jl") \ No newline at end of file diff --git a/src/models/eos/cubic/mixing/vdW1f.jl b/src/models/eos/cubic/mixing/vdW1f.jl new file mode 100644 index 000000000..5e74a2ec2 --- /dev/null +++ b/src/models/eos/cubic/mixing/vdW1f.jl @@ -0,0 +1,25 @@ +abstract type vdW1fRuleModel <: MixingRule end + +struct vdW1fRuleParam <: EoSParam +end + +@newmodelsimple vdW1fRule vdW1fRuleModel vdW1fRuleParam + +export vdW1fRule +function vdW1fRule(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false) + params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose) + acentricfactor = SingleParam(params["w"],"acentric factor") + packagedparams = vdW1fRuleParam() + model = vdW1fRule(packagedparams, verbose=verbose) + return model +end + +function mixing_rule(model::ABCubicModel,V,T,z,mixing_model::vdW1fRuleModel,α,a,b) + n = sum(z) + invn2 = (one(n)/n)^2 + b̄ = dot(z,Symmetric(b),z) * invn2 + ā = dot(z,Symmetric(a .* sqrt.(α*α')),z) * invn2 + return ā,b̄ +end + +is_splittable(::vdW1fRule) = false diff --git a/src/models/eos/ideal/BasicIdeal.jl b/src/models/eos/ideal/BasicIdeal.jl index f23173ad2..940bf8e4f 100644 --- a/src/models/eos/ideal/BasicIdeal.jl +++ b/src/models/eos/ideal/BasicIdeal.jl @@ -15,6 +15,7 @@ end function BasicIdeal(; userlocations::Array{String,1}=String[], verbose=false) return BasicIdeal(BasicIdealParam()) end +is_splittable(::BasicIdeal) = false function a_ideal(model::BasicIdeal, V, T, z) N = ∑(z) diff --git a/src/utils/macros.jl b/src/utils/macros.jl index a94ff1cfe..bd600243e 100644 --- a/src/utils/macros.jl +++ b/src/utils/macros.jl @@ -160,7 +160,7 @@ macro newmodelgc(name, parent, paramstype) built_by_macro(::Type{<:$name}) = true function Base.show(io::IO, mime::MIME"text/plain", model::$name) - return eosshow(io, mime, model) + return gc_eosshow(io, mime, model) end function Base.show(io::IO, model::$name) @@ -300,16 +300,26 @@ function (::Type{model})(params::EoSParam, absolutetolerance::Float64 = 1e-12, verbose::Bool = false) where model <:EoSModel - arbparam = arbitraryparam(params) - components = arbparam.components + if has_sites(model) + arbparam = arbitraryparam(params) + components = arbparam.components sites = SiteParam(components) return model(params,sites,idealmodel;ideal_userlocations,references,absolutetolerance,verbose) end #With sites out of the way, this is a simplemodel, no need to initialize the ideal model + #if there isnt any params, just put empty values. + if length(fieldnames(typeof(params))) > 0 + arbparam = arbitraryparam(params) + components = arbparam.components icomponents = 1:length(components) + else + components = String[] + icomponents =1:0 + end return model(components,icomponents,params,absolutetolerance,references) + end function init_model(idealmodel::IdealModel,components,userlocations,verbose) diff --git a/src/utils/split_model.jl b/src/utils/split_model.jl index ba4536ed2..e7789ff38 100644 --- a/src/utils/split_model.jl +++ b/src/utils/split_model.jl @@ -18,6 +18,17 @@ julia> split_model(gerg2) """ function split_model end + +""" + is_splittable(model)::Bool + +Trait to determine if a `EoSModel` should be splitted by itself or can be simply filled into a vector. +This is useful in the case of models without any parameters, as those models are impossible by definition to split, because they don't have any underlying data. + +The Default is `is_splittable(model) = true`. +""" +is_splittable(model) = true + function split_model(param::SingleParam{T}, splitter =split_model(1:length(param.components))) where T function generator(I) @@ -149,36 +160,31 @@ function auto_split_model(Base.@nospecialize(model::EoSModel)) allfields[:icomponents] = [1:1 for _ in 1:len_comps] end - #generate a vector of all params - if hasfield(typeof(model),:params) - allfields[:params] = split_model(model.params,splitter) - end - - if hasfield(typeof(model),:idealmodel) - if model.idealmodel == BasicIdeal() - allfields[:idealmodel] = fill(BasicIdeal(),len) + #process all model fields + modelfields = filter(x->getproperty(model,x) isa EoSModel,fieldnames(M)) + for modelkey in modelfields + modelx = getproperty(model,modelkey) + if is_splittable(modelx) + allfields[modelkey]= split_model(modelx) else - allfields[:idealmodel] = split_model(model.idealmodel) + allfields[modelkey] = fill(modelx,len_comps) end end - if hasfield(typeof(model),:references) - allfields[:references] = fill(model.references,len_comps) - end - if hasfield(typeof(model),:alpha) - if isnothing(model.alpha) - allfields[:alpha] = fill(nothing,len) - else - allfields[:alpha] = split_model(model.alpha) - end + #process all empty (Missing,Nothing) fields + emptyfields = filter(x->getproperty(model,x) isa Union{Nothing,Missing},fieldnames(M)) + + for emptykey in emptyfields + allfields[emptykey] = fill(model.emptykey,len_comps) + end + + + if hasfield(typeof(model),:params) + allfields[:params] = split_model(model.params,splitter) end - if hasfield(typeof(model),:activity) - if isnothing(model.activity) - allfields[:activity] = fill(nothing,len) - else - allfields[:activity] = split_model(model.activity) - end + if hasfield(typeof(model),:references) + allfields[:references] = fill(model.references,len_comps) end if hasfield(typeof(model),:absolutetolerance) @@ -187,9 +193,11 @@ function auto_split_model(Base.@nospecialize(model::EoSModel)) if hasfield(typeof(model),:sites) allfields[:sites] = split_model(model.sites,splitter) end + + return [M((allfields[k][i] for k in fieldnames(M))...) for i in 1:len] catch e - throw(e) + @show model return simple_split_model(model) end end diff --git a/src/utils/visualisation.jl b/src/utils/visualisation.jl index 9334c7166..eadbefda7 100644 --- a/src/utils/visualisation.jl +++ b/src/utils/visualisation.jl @@ -1,6 +1,6 @@ -function eosshow(io::IO, ::MIME"text/plain", model::EoSModel) +function eosshow(io::IO, ::MIME"text/plain", Base.@nospecialize(model::EoSModel)) print(io, typeof(model)) length(model) == 1 && println(io, " with 1 component:") length(model) > 1 && println(io, " with ", length(model), " components:") @@ -8,16 +8,18 @@ function eosshow(io::IO, ::MIME"text/plain", model::EoSModel) print(io, " \"", model.components[i], "\"") println(io) end - print(io,"Contains parameters: ") + paramnames = fieldnames(typeof(model.params)) + len_params = length(paramnames) + !iszero(len_params) && print(io,"Contains parameters: ") firstloop = true - for fieldname in fieldnames(typeof(model.params)) + for fieldname in paramnames firstloop == false && print(io, ", ") print(io, fieldname) firstloop = false end end -function eosshow(io::IO, model::EoSModel) +function eosshow(io::IO, Base.@nospecialize(model::EoSModel)) print(io, typeof(model)) firstloop = true print(io, "(") @@ -29,16 +31,16 @@ function eosshow(io::IO, model::EoSModel) print(io, ")") end -function Base.show(io::IO, ::MIME"text/plain", model::GCSAFTModel) +function gc_eosshow(io::IO, ::MIME"text/plain", Base.@nospecialize(model::EoSModel)) print(io, typeof(model)) length(model) == 1 && println(io, " with 1 component:") length(model) > 1 && println(io, " with ", length(model), " components:") - for i in model.icomponents + for i in 1:length(model) print(io, " \"", model.components[i], "\": ") firstloop = true - for k in 1:length(model.allcomponentgroups[i]) + for k in 1:length(model.groups.groups[i]) firstloop == false && print(io, ", ") - print(io, "\"", model.allcomponentgroups[i][k], "\" => ", model.allcomponentngroups[i][k]) + print(io, "\"", model.groups.groups[i][k], "\" => ", model.groups.n_groups[i][k]) firstloop = false end println(io) @@ -52,15 +54,5 @@ function Base.show(io::IO, ::MIME"text/plain", model::GCSAFTModel) end end -function Base.show(io::IO, model::GCSAFTModel) - print(io, typeof(model)) - firstloop = true - for i in model.icomponents - firstloop == false && print(io, ", ") - print(io, "\"", model.components[i], "\"") - firstloop = false - end - print(io, ")") -end export eosshow \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 3e4e6a7f6..9438f18e9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,8 @@ using Clapeyron, Test, Unitful @testset "All tests" begin include("test_database.jl") + include("test_solvers.jl") + include("test_differentials.jl") include("test_models.jl") include("test_methods.jl") end diff --git a/test/test_differentials.jl b/test/test_differentials.jl new file mode 100644 index 000000000..e64b01267 --- /dev/null +++ b/test/test_differentials.jl @@ -0,0 +1,79 @@ +struct TestModel <: EoSModel end + +function Clapeyron.eos(model::TestModel,V,T,z) + z_part = 1*z[1] + 2*z[2]+3*z[3] + v_part = log(V) + T_part = T^6 + vt_part = V*T*V*T + return vt_part+T_part+v_part+z_part +end + + +∂f∂T_analytical(model::TestModel,V,T,z) =6*T^5 +2*T*V*V +∂f∂V_analytical(model::TestModel,V,T,z) =1/V +2*V*T*T +∂2f∂T2_analytical(model::TestModel,V,T,z) =5*6*T^4 +2*V*V +∂2f∂V2_analytical(model::TestModel,V,T,z) =-1/(V*V) +2*T*T +∂2f∂V∂T_analytical(model::TestModel,V,T,z) =4*V*T +∂3f∂V3_analytical(model::TestModel,V,T,z) =2/(V*V*V) +∂3f∂V2∂T_analytical(model::TestModel,V,T,z) = 4*T +∂3f∂V∂T2_analytical(model::TestModel,V,T,z) = 4*V + +@testset "differentials" begin + model = TestModel() + T = 500*rand() + V = 10*rand() + z = rand(3) + f = Clapeyron.eos(model,V,T,z) + v = ∂f∂V_analytical(model,V,T,z) + t = ∂f∂T_analytical(model,V,T,z) + vv = ∂2f∂V2_analytical(model,V,T,z) + vt = ∂2f∂V∂T_analytical(model,V,T,z) + tt = ∂2f∂T2_analytical(model,V,T,z) + vvv = ∂3f∂V3_analytical(model,V,T,z) + + df = [v,t] + d2f = [vv vt;vt tt] + p = -v + pv = -vv + pvv = -vvv + pt = -vt + dp = [pv,pt] + pvt = -∂3f∂V2∂T_analytical(model,V,T,z) + ptt = -∂3f∂V∂T2_analytical(model,V,T,z) + d2p = [pvv pvt;pvt ptt] + + @testset "first order" begin + @test Clapeyron.∂f∂T(model,V,T,z) ≈ t + @test Clapeyron.∂f∂V(model,V,T,z) ≈ v + df1 = Clapeyron.∂f(model,V,T,z) + @test df1[2] ≈ f + @test df1[1][1] ≈ v + @test df1[1][2] ≈ t + end + + @testset "second order" begin + pdp = Clapeyron.p∂p∂V(model,V,T,z) + @test pdp[1] ≈ -v + @test pdp[2] ≈ -vv + h = Clapeyron.f_hess(model,V,T,z) + @test all(h .≈ d2f) + ddf = Clapeyron.∂2f(model,V,T,z) + @test ddf[3] ≈ f + @test all(ddf[2] .≈ df) + @test all(ddf[1] .≈ d2f) + end + + @testset "third order" begin + dddf = Clapeyron.∂²³f(model,V,T,z) + @test dddf[1] ≈ -vv + @test dddf[2] ≈ -vvv + + ddp = Clapeyron.∂2p(model,V,T,z) + @test ddp[3] ≈ p + @test all(ddp[2] .≈ dp) + @test all(ddp[1] .≈ d2p) + end + +end + + diff --git a/test/test_solvers.jl b/test/test_solvers.jl new file mode 100644 index 000000000..f19748278 --- /dev/null +++ b/test/test_solvers.jl @@ -0,0 +1,47 @@ +const SOL = Clapeyron.Solvers + +quadratic(x) = x*x - 4 +rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2 +golden_number_fixpoint(x) = one(x) + one(x)/x +@testset "ad_newton" begin + x0 = 2+rand() + fdf = SOL.f∂f(quadratic,x0) + @test fdf[1] ≈ quadratic(x0) + @test fdf[2] ≈ 2*x0 + @test SOL.ad_newton(quadratic,x0) ≈ 2.0 +end + +@testset "box optimize" begin + #example of FMinBox in Optim.jl + lower = [1.25, -2.1] + upper = [Inf, Inf] + solution = [1.25, 1.5625] + initial_x = [2.0, 2.0] + res = SOL.box_optimize(rosenbrock,initial_x,lower,upper) + @test all(res.info.minimizer .≈ solution) + @test res.info.minimum ≈ 0.0625 +end + +@testset "fixpoint" begin + #example of FMinBox in Optim.jl + x0 = 2+ rand() + ϕ = (1+sqrt(5))/2 + @test SOL.fixpoint(golden_number_fixpoint,x0) ≈ ϕ +end + +function f_diffmcp!(fvec, x) + fvec[1] = (1-x[1])^2-1.01 +end + +# A difficult MCP. +# +# Presented in Miranda and Fackler (2002): "Applied Computational Economics and +# Finance", p. 51 +#does not converge in NLSolve.jl converges here with NLSolvers.jl +@testset "nlsolve" begin + res = SOL.nlsolve(f_diffmcp!,[0.0]) + solution = SOL.x_sol(res) + v = [1.0] + f_diffmcp!(v,solution) + @test v[1] <= sqrt(eps(Float64)) +end \ No newline at end of file