add a bunch of tests for flash (#320)

src/methods/PH.jl

@@ -9,7 +9,6 @@ function PH_property(model,p,h,z,f::F,phase,T0,threaded) where F
 
     if !is_unknown(phase)
         T,calc_phase = _Tproperty(model,p,h,z,T0 = T0,phase = phase,threaded = threaded)
-        #TODO: handle equilibria conditions
         if calc_phase != :eq && calc_phase != :failure
             return f(model,p,T,z;phase = calc_phase)
         elseif calc_phase == :eq && !supports_lever_rule(f)
@@ -22,14 +21,10 @@ function PH_property(model,p,h,z,f::F,phase,T0,threaded) where F
         end
     else
         res = ph_flash(model,p,h,z,T0 = T0)
-        np = numphases(res)
         if f == temperature
             return temperature(res)
-        end
-        if np > 1
-            return f(model,res)
         else
-            return f(model,p,T,z)
+            return f(model,res)
         end
     end
 end

src/methods/PS.jl

@@ -6,17 +6,26 @@ function PS_property(model,p,s,z,f::F,phase,T0,threaded) where F
     if f == pressure
         return p
     end
-    T,calc_phase = _Tproperty(model,p,s,z,entropy,T0 = T0,phase = phase,threaded = threaded)
-        #TODO: handle equilibria conditions
-    if calc_phase != :eq || calc_phase != :failure
-        return f(model,p,T,z;phase = calc_phase)
-    elseif calc_phase == :eq && !supports_lever_rule(f)
-        thow(invalid_property_multiphase_error(f))
-    elseif calc_phase == :eq && supports_lever_rule(f) && length(model) == 1 #TODO: lift single component condition
-        result = ps_flash(model,p,s,z,T0)
-        return f(model,result)
+
+    if !is_unknown(phase)
+        T,calc_phase = _Tproperty(model,p,h,z,entropy,T0 = T0,phase = phase,threaded = threaded)
+        if calc_phase != :eq && calc_phase != :failure
+            return f(model,p,T,z;phase = calc_phase)
+        elseif calc_phase == :eq && !supports_lever_rule(f)
+            thow(invalid_property_multiphase_error(f))
+        elseif calc_phase == :eq && supports_lever_rule(f)
+            result = ps_flash(model,p,s,z,T0)
+            return f(model,result)
+        else
+            return f(model,p,T,z;phase = phase)
+        end
     else
-        return f(model,p,T,z;phase = phase)
+        res = ps_flash(model,p,s,z,T0 = T0)
+        if f == temperature
+            return temperature(res)
+        else
+            return f(model,res)
+        end
     end
 end
 

src/methods/property_solvers/multicomponent/flash.jl

@@ -274,6 +274,12 @@ is_vle(method::FlashMethod) = is_vle(method.equilibrium)
 is_lle(method::FlashMethod) = is_lle(method.equilibrium)
 is_unknown(method::FlashMethod) = is_unknown(method.equilibrium)
 
+@noinline function incorrect_np_flash_error(method,result)
+    np = numphases(result)
+    s = np == 1 ? "" : "s"
+    throw(ArgumentError("$method does not support an input with $np phase$s as an initial point. Got the following input: \n\n $result"))
+end
+
 include("flash/general_flash.jl")
 include("flash/PT.jl")
 include("flash/PH.jl")

src/methods/property_solvers/multicomponent/flash/QT.jl

@@ -10,7 +10,7 @@ function qt_flash_x0(model,β,T,z,method::FlashMethod)
             x = z ./ sum(z)
             p,vl,vv,y = __x0_bubble_pressure(model,T,x)
             y ./= sum(y)
-            
+
             return FlashResult(p,T,SA[x,y],SA[1.0-β,1.0*β],SA[vl,vv],sort = false)
         elseif 0.99 <= β <= 1.0
             y = z ./ sum(z)
@@ -20,12 +20,23 @@ function qt_flash_x0(model,β,T,z,method::FlashMethod)
         else
             pures = split_model(model)
             sat = extended_saturation_pressure.(pures,T)
-            ps = first.(sat)    
+            ps = first.(sat)
             K = similar(ps)
             pmin,pmax = extrema(ps)
-            p0 = β*pmin + (1-β)*pmax
             fp(p) = qt_f0_p!(K,z,p,ps,β)
-            prob = Roots.ZeroProblem(fp,p0)
+            pm = β*pmin + (1-β)*pmax
+            pr1 = range(pmin,pm,5*length(model))
+            pr2 = range(pm,pmax,5*length(model))
+            δβ1 = abs.(fp.(pr1))
+            δβ2 = abs.(fp.(pr2))
+            δβ1_min,i1 = findmin(δβ1)
+            δβ2_min,i2 = findmin(δβ2)
+            if δβ1_min < δβ2_min
+                p00 = pr1[i1]
+            else
+                p00 = pr2[i2]
+            end
+            prob = Roots.ZeroProblem(fp,p00)
             p = Roots.solve(prob)
         end
     else
@@ -40,7 +51,7 @@ function qt_flash(model::EoSModel,β,T,z;kwargs...)
     return qt_flash(model,β,T,z,method)
 end
 
-function init_preferred_method(method::typeof(qt_flash),model::EoSModel,kwargs) 
+function init_preferred_method(method::typeof(qt_flash),model::EoSModel,kwargs)
     GeneralizedXYFlash(;kwargs...)
 end
 
@@ -58,7 +69,7 @@ function qt_flash(model,β,T,z,method::FlashMethod)
         result1 = βflash_pure(model,temperature,T,βv,z)
         return index_expansion(result1,idx_r)
     end
-    
+
     result = qt_flash_impl(model_r,β,T,z_r,method_r)
     if !issorted(result.volumes)
         #this is in case we catch a bad result.

src/methods/property_solvers/multicomponent/flash/general_flash.jl

@@ -387,7 +387,7 @@ function xy_flash_neq(output,model,zbulk,np,input,state::F,μconfig) where F
         i_spec = spec1.k
         output[end-1] = β[i_spec] - val1
     else #general caloric term, valid for enthalpy , entropy, internal energy, gibbs, helmholtz
-        output[end-1] = (val1 - ∑a - p_end*∑v - T*∑s)*RTinv
+        output[end-1] = (val1 - ∑a - p_end*∑v - T*∑s)/val1
     end
 
     if spec2 == temperature
@@ -400,7 +400,7 @@ function xy_flash_neq(output,model,zbulk,np,input,state::F,μconfig) where F
         i_spec = spec2.k
         output[end-1] = β[i_spec] - val2
     else
-        output[end] = (val2 - ∑a - p_end*∑v - T*∑s)*RTinv
+        output[end] = (val2 - ∑a - p_end*∑v - T*∑s)/val1
     end
     return output
 end
@@ -431,7 +431,7 @@ function xy_flash(model::EoSModel,spec::FlashSpecifications,z,flash::FlashResult
     return xy_flash(model,spec,z,comps0,β0,volumes0,T0;rtol,atol,max_iters)
 end
 
-function xy_flash(model::EoSModel,spec::FlashSpecifications,z,flash::FlashResult;rtol = 1e-12,atol = 1e-10,max_iters = 50)
+function xy_flash(model::EoSModel,spec::FlashSpecifications,z,flash::FlashResult;rtol = 1e-14,atol = 1e-12,max_iters = 50)
     if numphases(flash) == 1
         throw(ArgumentError("xy_flash cannot use single phase initial points as starting points."))
     end
@@ -542,15 +542,13 @@ function xy_flash(model::EoSModel,spec::FlashSpecifications,z,comps0,β0,volumes
         x_old .= x
         #bound positivity
         α0 = Solvers.positive_linesearch(x,s)
-
         #backtrack linesearch, so the next result is strictly better than the last
         α,Θx = Solvers.backtracking_linesearch!(Θ,F,x_old,s,Θx,x,α0)
         Fnorm = sqrt(2*Θx)
         spec_norm = norm(viewlast(F,2),Inf)
         snorm_old = snorm
         snorm = α*norm(s,2)
         δs = max(abs(snorm-snorm_old),norm((F[end-1],F[end]),Inf))
-
         Fnorm = norm(F,Inf)
         δs < rtol && (converged = true)
         xnorm = Solvers.dnorm(x,x_old,Inf)
@@ -561,11 +559,10 @@ function xy_flash(model::EoSModel,spec::FlashSpecifications,z,comps0,β0,volumes
         isnan(δs) && (nan_converged = true)
         i == max_iters && (max_iters_reached = true)
     end
-
+    max_iters_reached,converged,nan_converged
     if !converged || nan_converged || max_iters_reached
         x .= NaN
     end
-
     if spec_norm > sqrt(rtol)
         x .= NaN
     end
@@ -639,14 +636,15 @@ function GeneralizedXYFlash(;equilibrium = :unknown,
                         x0 = nothing,
                         y0 = nothing,
                         v0 = nothing,
-                        rtol = 1e-12,
-                        atol = 1e-10,
+                        rtol = 1e-14,
+                        atol = 1e-12,
                         max_iters = 100,
                         flash_result = nothing)
     !(is_vle(equilibrium) | is_lle(equilibrium) | is_unknown(equilibrium))  && throw(error("invalid equilibrium specification for GeneralizedXYFlash"))
     if flash_result isa FlashResult
         comps,β,volumes = flash_result.compositions,flash_result.fractions,flash_result.volumes
-        @assert length(comps) == 2
+        np = numphases(flash_result)
+        np != 2 && incorrect_np_flash_error(GeneralizedXYFlash,flash_result)
         w1,w2 = comps[1],comps[2]
         v = (volumes[1],volumes[2])
         P00 = flash_result.data.p

src/methods/property_solvers/multicomponent/tp_flash/Michelsentp_flash.jl

@@ -61,7 +61,8 @@ function MichelsenTPFlash(;equilibrium = :unknown,
 
     if flash_result isa FlashResult
         comps,β,volumes = flash_result.compositions,flash_result.fractions,flash_result.volumes
-        @assert length(comps) == 2
+        np = numphases(flash_result)
+        np != 2 && incorrect_np_flash_error(MichelsenTPFlash,flash_result)
         w1,w2 = comps[1],comps[2]
         v = (volumes[1],volumes[2])
         return Michelsentp_flash(;equilibrium,x0 = w1,y0 = w2,vol0 = v,K_tol,ss_iters,nacc,second_order,noncondensables,nonvolatiles)

src/methods/property_solvers/multicomponent/tp_flash/RachfordRicetp_flash.jl

@@ -55,6 +55,10 @@ function RRTPFlash(;equilibrium = :unknown,
     nonvolatiles = nothing,
     flash_result = nothing)
     ss_iters = max_iters
+    if flash_result isa FlashResult
+        np = numphases(flash_result)
+        np != 2 && incorrect_np_flash_error(RRTPFlash,flash_result)
+    end
     #we call Michelsen to check if the arguments are correct.
     m = MichelsenTPFlash(;equilibrium,K0,x0,y0,v0,K_tol,ss_iters,nacc,noncondensables,nonvolatiles,flash_result)
     return RRTPFlash{eltype(m)}(m.equilibrium,m.K0,m.x0,m.y0,m.v0,m.K_tol,max_iters,m.nacc,m.noncondensables,m.nonvolatiles)

src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl

@@ -58,9 +58,11 @@ function MultiPhaseTPFlash(;
 
     if flash_result !== nothing
         comps = flash_result.compositions
+        np = numphases(flash_result)
+        np < 2 && incorrect_np_flash_error(MultiPhaseTPFlash,flash_result)
         ∑n = sum(flash_result.fractions)
         n00 = Vector{eltype(comps)}[]
-        for i in 1:numphases(flash_result)
+        for i in 1:np
             push!(n00,collect(comps[i]))
         end
         return MultiPhaseTPFlash(;K0 = nothing,x0 = nothing, y0 = nothing,n0 = n00,K_tol,ss_iters,nacc,second_order,max_phases,phase_iters)

test/test_methods_api.jl

@@ -363,6 +363,41 @@ end
     end
 end
 
+@testset "XY flash" begin
+    #1 phase (#320)
+    model = cPR(["ethane","methane"],idealmodel = ReidIdeal)
+    p = 101325.0
+    z = [1.2,1.2]
+    T = 350.0
+    h = enthalpy(model,p,T,z)
+    res0 = ph_flash(model,p,h,z)
+    @test Clapeyron.temperature(res0) ≈ T rtol = 1e-6
+    @test enthalpy(model,res1) ≈ h rtol = 1e-6
+
+    #2 phases
+    h = -13831.0
+    res1 = ph_flash(model,p,h,z)
+    @test enthalpy(model,res1) ≈ h rtol = 1e-6
+
+
+    #examples for qt, qp flash (#314)
+    model = cPR(["ethane","propane"],idealmodel=ReidIdeal)
+    res2 = qt_flash(model,0.5,208.0,[0.5,0.5])
+    @test Clapeyron.pressure(res2) ≈ 101634.82435966855 rtol = 1e-6
+    res3 = qp_flash(model,0.5,120000.0,[0.5,0.5])
+    @test Clapeyron.temperature(res3) ≈ 211.4972567716822 rtol = 1e-6
+
+    #1 phase input should error
+    model = PR(["IsoButane", "n-Butane", "n-Pentane", "n-Hexane"])
+    z = [0.25, 0.25, 0.25, 0.25]
+    p = 1e5
+    h = enthalpy(model, 1e5, 303.15, z)
+    r = Clapeyron.ph_flash(model, p, h, z)
+    @test_throws ArgumentError qt_flash(model,0.5,308,z,flash_result = r)
+    res4 = qp_flash(model,0.7,60000.0,z)
+    @test Clapeyron.temperature(res4)  ≈ 289.6991395244328 rtol = 1e-6
+end
+
 @testset "Saturation Methods" begin
     model = PR(["water"])
     vdw = vdW(["water"])