Flash improvements + flash bypass infrastructure

.github/workflows/CI.yml

@@ -13,8 +13,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
           - '1.9'
+          - '1'
         os:
           - ubuntu-latest
         arch:

Project.toml

@@ -1,7 +1,7 @@
 name = "MultiComponentFlash"
 uuid = "35e5bd01-9722-4017-9deb-64a5d32478ff"
 authors = ["Olav Møyner <[email protected]>"]
-version = "1.1.11"
+version = "1.1.12"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

docs/src/examples/basics.md

@@ -116,20 +116,22 @@ we see that the vapor phase is almost entirely made up of the lighter methane at
 
 ## Changing the solution algorithm to Newton
 
-If we examine the third output, we can see output about number of iterations and a verification that the flash converged within the default tolerance:
+If we examine the third output, we can see output about number of iterations, the result of the phase stability test and a verification that the flash converged within the default tolerance:
 
 ```jldoctest
 julia> report
-(its = 8, converged = true)
+(its = 8, converged = true, stability = StabilityReport (two-phase, liquid = stable, vapor = unstable))
 ```
 
-The default method is the `SSIFlash()` method. The name is an abbreviation for successive-substitution, a simple but very robust method.
+Here, we can see that forming a second vapor-like phase led to the solver to conclude that the mixture is two-phase.
+
+The default flash method is the `SSIFlash()` method. The name is an abbreviation for successive-substitution, a simple but very robust method.
 
 We could alternatively switch to `NewtonFlash()` to use Newton's method with AD instead to reduce the number of iterations:
 
 ```jldoctest
 julia> V, K, report = flash_2ph(eos, conditions, extra_out = true, method = NewtonFlash()); report
-(its = 5, converged = true)
+(its = 5, converged = true, stability = StabilityReport (two-phase, liquid = stable, vapor = unstable))
 ```
 
-It is also possible to use `SSINewtonFlash()` that switches from SSI to Newton after a prescribed number of iterations, which is effective around the critical region where SSI has slow convergence.
+Note that Newton's method is not unconditionally stable: It is therefore also possible to use `SSINewtonFlash()` that switches from SSI to Newton after a prescribed number of iterations, which is effective around the critical region where SSI has slow convergence.

src/eos.jl

@@ -25,30 +25,33 @@ end
 
 # PengRobinson specialization
 function static_coefficients(::AbstractPengRobinson)
-    return (0.457235529, 0.077796074, 1 + sqrt(2), 1 - sqrt(2))
+    return (0.457235529, 0.077796074, 1.0 + sqrt(2), 1.0 - sqrt(2))
 end
 
 function weight_ai(eos::GenericCubicEOS{T}, cond, i) where T<:AbstractPengRobinson
     mix = eos.mixture
     m = molecular_property(mix, i)
     a = acentric_factor(m)
     T_r = reduced_temperature(mix, cond, i)
-    return eos.ω_a*(1 + (0.37464 + 1.54226*a - 0.26992*a^2)*(1-T_r^(1/2)))^2;
+    tmp = (1.0 + (0.37464 + 1.54226*a - 0.26992*a*a)*(1-T_r^0.5))
+    return eos.ω_a*(tmp*tmp)
 end
 
 function weight_ai(eos::GenericCubicEOS{PengRobinsonCorrected}, cond, i)
     mix = eos.mixture
     m = molecular_property(mix, i)
     a = acentric_factor(m)
     T_r = reduced_temperature(mix, cond, i)
+    aa = a*a
     if a > 0.49
         # Use alternate expression.
-        D = (0.379642 + 1.48503*a - 0.164423*a^2 + 0.016666*a^3)
+        D = (0.379642 + 1.48503*a - 0.164423*aa + 0.016666*a*aa)
     else
         # Use standard expression.
-        D = (0.37464 + 1.54226*a - 0.26992*a^2)
+        D = (0.37464 + 1.54226*a - 0.26992*aa)
     end
-    return eos.ω_a*(1 + D*(1-T_r^(1/2)))^2;
+    tmp = 1.0 + D*(1.0-T_r^0.5)
+    return eos.ω_a*(tmp*tmp);
 end
 
 # ZudkevitchJoffe
@@ -57,7 +60,7 @@ function weight_ai(eos::GenericCubicEOS{ZudkevitchJoffe}, cond, i)
     mix = eos.mixture
     T = cond.T
     T_r = reduced_temperature(mix, cond, i)
-    return eos.ω_a*zj.F_a(T, i)*T_r^(-1/2)
+    return eos.ω_a*zj.F_a(T, i)*T_r^(-0.5)
 end
 
 function weight_bi(eos::GenericCubicEOS{ZudkevitchJoffe}, cond, i)
@@ -272,14 +275,14 @@ end
 """
 Fugacity
 """
-function component_fugacity_coefficient(eos::AbstractCubicEOS, cond, i, Z, forces, scalars)
+Base.@propagate_inbounds function component_fugacity_coefficient(eos::AbstractCubicEOS, cond, i, Z, forces, scalars)
     # NOTE: This returns ln(ψ), not ψ!
     m1 = eos.m_1
     m2 = eos.m_2
     return component_fugacity_coefficient_cubic(m1, m2, cond.z, Z, scalars.A, scalars.B, forces.A_ij, forces.B_i, i)
 end
 
-function component_fugacity_coefficient_cubic(m1, m2, x, Z, A, B, A_mat, B_i, i)
+Base.@propagate_inbounds function component_fugacity_coefficient_cubic(m1, m2, x, Z, A, B, A_mat, B_i, i)
     Δm = m1 - m2
     T = promote_type(eltype(x), eltype(A_mat))
     A_s = zero(T)

src/flash.jl

@@ -31,7 +31,7 @@ Two outcomes are possible:
 
 See also: [`flash_2ph!`](@ref), [`single_phase_label`](@ref)
 """
-function flash_2ph(eos, c::T, K = initial_guess_K(eos, c); method = SSIFlash(), kwarg...) where T
+function flash_2ph(eos, c::T, K = initial_guess_K(eos, c), V = NaN; method = SSIFlash(), kwarg...) where T
     nc = number_of_components(eos)
     @assert hasfield(T, :p)
     @assert hasfield(T, :T)
@@ -41,7 +41,7 @@ function flash_2ph(eos, c::T, K = initial_guess_K(eos, c); method = SSIFlash(),
     method::AbstractFlash
 
     S = flash_storage(eos, c, method = method)
-    flash_2ph!(S, K, eos, c, update_forces = false; method = method, kwarg...)
+    flash_2ph!(S, K, eos, c, V, update_forces = false; method = method, kwarg...)
 end
 
 """
@@ -59,25 +59,49 @@ Remaining arguments documented in [`flash_2ph`](@ref).
 - `update_forces = true`: Update the `p, T` dependent forces in storage initially.
 
 """
-function flash_2ph!(storage, K, eos, c, V = NaN;
-                                method = SSIFlash(), verbose = false, maxiter = 20000,
-                                tolerance = 1e-8, extra_out = false, update_forces = true,
-                                check = true, z_min = MINIMUM_COMPOSITION, kwarg...)
+function flash_2ph!(arg...; extra_out = false, kwarg...)
+    out = flash_2ph_impl!(arg...; kwarg...)
+    if extra_out
+        return out
+    else
+        return out[1]
+    end
+end
+
+function flash_2ph_impl!(storage, K, eos, c, V = NaN;
+        method = SSIFlash(),
+        verbose::Bool = false,
+        maxiter::Int = 20000,
+        tolerance::Float64 = 1e-8,
+        extra_out::Bool = false,
+        update_forces::Bool = true,
+        check::Bool = true,
+        z_min = MINIMUM_COMPOSITION,
+        kwarg...
+    )
     z = c.z
     if !isnothing(z_min)
         for i in eachindex(z)
             @inbounds z[i] = max(z[i], z_min)
         end
     end
-    @assert all(isfinite, K)
     forces = storage.forces
     if update_forces
         force_coefficients!(forces, eos, c)
     end
     single_phase_init = isnan(V) || V == 1.0 || V == 0.0
     if single_phase_init
-        stable = stability_2ph!(storage, K, eos, c, maxiter = maxiter, verbose = verbose; kwarg...)
+        stable, stability_report = stability_2ph!(storage, K, eos, c;
+            maxiter = maxiter,
+            verbose = verbose,
+            extra_out = true,
+            kwarg...
+        )
     else
+        # We check this here - if the stability test is performed, K is
+        # already overwritten with a reasonable finite initial guess.
+        @assert all(isfinite, K) "K values must be finite: K = $K"
+        stability_report = StabilityReport(stable_liquid = false, stable_vapor = false)
         stable = false
     end
     converged = false
@@ -88,29 +112,26 @@ function flash_2ph!(storage, K, eos, c, V = NaN;
         if isnan(V)
             V = solve_rachford_rice(K, z, V)
         end
-        while !converged
+        while true
             V, ϵ = flash_update!(K, storage, method, eos, c, forces, V, i)
             converged = ϵ ≤ tolerance
-            if converged
+            if converged || i == maxiter
                 if verbose
-                    @info "Flash done in $i iterations." V K
+                    @info "Flash done in $i iterations." V K converged
                 end
-                break
-            end
-            if check && i == maxiter
-                @warn "Flash did not converge in $i iterations. Final ϵ = $ϵ > $tolerance = tolerance" c
-                if !isfinite(V)
-                    error("Bad vapor fraction from flash")
+                if check && !converged
+                    @warn "Flash did not converge in $i iterations. Final ϵ = $ϵ > $tolerance = tolerance" c
+                    if !isfinite(V)
+                        # Hard error
+                        error("Bad vapor fraction from flash")
+                    end
                 end
+                break
             end
             i += 1
         end
     end
-    if extra_out
-        return (V, K, (its = i, converged = converged))
-    else
-        return V
-    end
+    return (V, K, (its = i, converged = converged, stability = stability_report))
 end
 
 """
@@ -391,14 +412,11 @@ function update_flash_jacobian!(J, r, eos, p, T, z, x, y, V, forces)
     Z_l, s_l = prep(eos, liquid, forces,:liquid)
     Z_v, s_v = prep(eos, vapor, forces,:vapor)
 
-    if isa(V, ForwardDiff.Dual)
-        np = length(V.partials)
-    else
-        np = length(p.partials)
-    end
+    p, V = Base.promote(p, V)
+    np = length(V.partials)
     # Isofugacity constraint
     # f_li - f_vi ∀ i
-    @inbounds for c in 1:n
+    @inbounds for c in eachindex(z)
         f_l = component_fugacity(eos, liquid, c, Z_l, forces, s_l)
         f_v = component_fugacity(eos, vapor, c, Z_v, forces, s_v)
         Δf = f_l - f_v
@@ -412,9 +430,9 @@ function update_flash_jacobian!(J, r, eos, p, T, z, x, y, V, forces)
     # x_i*(1-V) - V*y_i - z_i = 0 ∀ i
     # Σ_i x_i - y_i = 0
     T = Base.promote_type(eltype(x), eltype(y))
-    L = 1 - V
+    L = one(T) - V
     Σxy = zero(T)
-    @inbounds for c in 1:n
+    @inbounds for c in eachindex(z)
         xc, yc = x[c], y[c]
         Σxy += (xc - yc)
         M = L*xc + V*yc - z[c]

src/flash_types.jl

@@ -40,3 +40,51 @@ See also: [`flash_2ph!`](@ref), [`SSIFlash`](@ref) [`NewtonFlash`](@ref)
     dMax::Float64 = 0.2
 end
 
+
+struct PhaseStabilityStatus
+    stable::Bool
+    trivial::Bool
+    function PhaseStabilityStatus(stable = false; trivial = stable)
+        return new(stable, trivial && stable)
+    end
+end
+
+function Base.show(io::IOContext, sr::PhaseStabilityStatus)
+    compact = get(io, :compact, false)
+    s = sr.stable ? "single-phase" : "two-phase"
+    e = sr.trivial ? ", trivial" : ""
+    if compact
+        print(io, s)
+    else
+        print(io, "PhaseStabilityStatus ($s$e)")
+    end
+end
+
+struct StabilityReport
+    stable::Bool
+    liquid::PhaseStabilityStatus
+    vapor::PhaseStabilityStatus
+    function StabilityReport(;
+            stable_liquid::Bool = false,
+            trivial_liquid::Bool = stable_liquid,
+            stable_vapor::Bool = false,
+            trivial_vapor::Bool = stable_vapor,
+        )
+        new(stable_liquid && stable_vapor,
+            PhaseStabilityStatus(stable_liquid, trivial = trivial_liquid),
+            PhaseStabilityStatus(stable_vapor, trivial = trivial_vapor)
+        )
+    end
+end
+
+function Base.show(io::IOContext, sr::StabilityReport)
+    compact = get(io, :compact, false)
+    s = sr.stable ? "single-phase" : "two-phase"
+    ls = sr.liquid.stable ? "stable" : "unstable"
+    vs = sr.vapor.stable ? "stable" : "unstable"
+    if compact
+        print(io, s)
+    else
+        print(io, "StabilityReport ($s, liquid = $ls, vapor = $vs)")
+    end
+end

src/flow_coupler_types.jl

@@ -41,10 +41,19 @@ struct FlashedMixture2Phase{T, A<:AbstractVector{T}, E}
     liquid::FlashedPhase{T, A}
     vapor::FlashedPhase{T, A}
     critical_distance::Float64
-    p::Float64
-    T::Float64
-    function FlashedMixture2Phase(state::PhaseState2Phase, K::K_t, V::V_t, liquid, vapor; vec_type = Vector{V_t}, critical_distance = NaN, p = NaN, T = NaN) where {V_t, K_t}
-        new{V_t, vec_type, K_t}(state, K, V, liquid, vapor, critical_distance, p, T)
+    flash_cond::@NamedTuple{p::Float64, T::Float64, z::E}
+    flash_stability::StabilityReport
+    function FlashedMixture2Phase(state::PhaseState2Phase, K::K_t, V::V_t, liquid, vapor;
+            vec_type = Vector{V_t},
+            critical_distance = NaN,
+            cond = missing,
+            stability_report = StabilityReport()
+        ) where {V_t, K_t}
+        if ismissing(cond)
+            z0 = convert(K_t, fill(NaN, length(K)))
+            cond = (p = NaN, T = NaN, z = z0)
+        end
+        new{V_t, vec_type, K_t}(state, K, V, liquid, vapor, critical_distance, cond, stability_report)
     end
 end
 
@@ -57,22 +66,32 @@ function Base.convert(::Type{FlashedMixture2Phase{T, Vector{T}, F}}, mixture::Fl
     V = to_T(mixture.V)
     # Kv = map(to_T, mixture.K)
     Kv = mixture.K
-    converted_mixture = FlashedMixture2Phase(mixture.state, Kv, V, liquid, vapor; vec_type = Vector{T}, critical_distance = mixture.critical_distance)
+    converted_mixture = FlashedMixture2Phase(
+        mixture.state,
+        Kv,
+        V,
+        liquid,
+        vapor;
+        vec_type = Vector{T},
+        critical_distance = mixture.critical_distance,
+        cond = mixture.cond,
+        stability_report = mixture.flash_stability
+    )
     return converted_mixture
 end
 
-function FlashedMixture2Phase(state, K, V, x, y, Z_L, Z_V, b = NaN, p = NaN, Temp = NaN)
+function FlashedMixture2Phase(state, K, V, x, y, Z_L, Z_V, b = NaN, cond = missing, stability = StabilityReport())
     liquid = FlashedPhase(x, Z_L)
     vapor = FlashedPhase(y, Z_V)
-    return FlashedMixture2Phase(state, K, V, liquid, vapor, critical_distance = b, p = p, T = Temp)
+    return FlashedMixture2Phase(state, K, V, liquid, vapor, critical_distance = b, cond = cond, stability_report = stability)
 end
 
-function FlashedMixture2Phase(eos::AbstractEOS, T = Float64, T_num = Float64, b = NaN, p = NaN, Temp = NaN)
+function FlashedMixture2Phase(eos::AbstractEOS, T = Float64, T_num = Float64, b = NaN, cond = missing, stability = StabilityReport())
     n = number_of_components(eos)
     V = zero(T)
     # K values are always doubles
     K = zeros(T_num, n)
     liquid = FlashedPhase(n, T)
     vapor = FlashedPhase(n, T)
-    return FlashedMixture2Phase(unknown_phase_state_lv, K, V, liquid, vapor, critical_distance = b, p = p, T = Temp)
+    return FlashedMixture2Phase(unknown_phase_state_lv, K, V, liquid, vapor, critical_distance = b, cond = cond, stability_report = stability)
 end