Support phase stability bypass in immutable, bump version

Project.toml

@@ -1,6 +1,6 @@
 name = "MultiComponentFlash"
 uuid = "35e5bd01-9722-4017-9deb-64a5d32478ff"
-version = "1.1.20"
+version = "1.2.0"
 authors = ["Olav Møyner <olav.moyner@gmail.com>"]
 
 [deps]

docs/src/examples/advanced.md

@@ -79,6 +79,27 @@ conditions_i = (p = pressure[i], T = temperature[i], z = z_static)
 V, K = flash_2ph_immutable(eos_static, conditions_i, storage)
 ```
 
+For a sequence of nearby states, the immutable Michelsen stability bypass can
+reuse the last fully tested single-phase condition:
+
+```julia
+V, K, stability = flash_2ph_immutable(eos_static, conditions_static;
+    return_stability = true)
+
+next_conditions = (p = 1.001p, T = T + 0.01,
+    z = SVector{length(z)}(z))
+V, K, stability = flash_2ph_immutable(eos_static, next_conditions;
+    stability_storage = stability,
+    return_stability = true)
+```
+
+The stability calculation is also available on its own with
+`stability_2ph_immutable`. Its result and nested storage are isbits values and
+can be passed through accelerator kernels. `stability.bypassed` indicates
+whether the full Michelsen test was skipped. Storage is only armed when both
+trial phases converged to trivial stable solutions; the shadow region is
+always retested.
+
 Do not share ordinary mutable storage from `flash_storage(...; static = false)`
 between kernel work items.
 

src/MultiComponentFlash.jl

@@ -20,7 +20,8 @@ module MultiComponentFlash
     export number_of_components
     # Flash interfaces
     export flash_2ph, flash_2ph!, flash_2ph_immutable, flash_storage
-    export stability_2ph, stability_2ph!
+    export stability_2ph, stability_2ph!, stability_2ph_immutable
+    export StaticStabilityStorage, StaticStabilityResult
     # Algorithms for flash
     export SSIFlash, NewtonFlash, SSINewtonFlash
     # Mixtures and their molecular makeup

src/eos_types.jl

@@ -182,10 +182,10 @@ function GenericCubicEOS(setup::NamedTuple, mixture; volume_shift = nothing)
     return GenericCubicEOS(setup.type, mixture, setup.m_1, setup.m_2, setup.ω_a, setup.ω_b, volume_shift)
 end
 
-struct KValuesEOS{T, R, N, V} <: AbstractEOS
+struct KValuesEOS{T, R, N, V, M<:MultiComponentMixture{R, N}} <: AbstractEOS
     "Callable on the form `cond -> V` or a set of constants (Tuple/AbstractVector)"
     K_values_evaluator::T
-    mixture::MultiComponentMixture{R, N}
+    mixture::M
     volume_shift::V
 end
 

src/flash.jl

@@ -170,8 +170,11 @@ See also: [`flash_2ph!`](@ref) [`set_partials`](@ref)
 """
 function flash_storage(eos, cond = (p = 10e5, T = 273.15, z = zeros(number_of_components(eos)));
         method = SSIFlash(), static::Bool = false, static_size = nothing, kwarg...)
-    isnothing(static_size) || throw(ArgumentError(
-        "`static_size` has been replaced by `static`; use `static=true` for the fully static path."))
+    if !isnothing(static_size)
+        Base.depwarn("`static_size` is deprecated; use `static=$(static_size)` instead.",
+            :flash_storage)
+        static = static_size
+    end
     config = static ? StaticConfig() : FlashConfig()
     return flash_storage(eos, cond, method, config; kwarg...)
 end
@@ -304,11 +307,11 @@ function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces) where {F<:Real}
         ϵ = max(ϵ, abs(1-r))
     end
     V = solve_rachford_rice(K, z, V)
+    V = clamp(V, zero(V), one(V))
     return (V, ϵ)::Tuple{F, F}
 end
 
 cap_z(z) = min(max(z, MINIMUM_COMPOSITION), one(z))
-cap_unit(v) = min(max(v, zero(z)), one(z))
 cap_VL(v) = min(max(v, MINIMUM_COMPOSITION), 1 - MINIMUM_COMPOSITION)
 
 function flash_update!(K, storage, type::NewtonFlash, eos, cond, forces, V, iteration)

src/flash_types.jl

@@ -112,3 +112,35 @@ function Base.show(io::IOContext, sr::StabilityReport)
         print(io, "StabilityReport ($s, liquid = $ls, vapor = $vs)")
     end
 end
+
+"""
+    StaticStabilityStorage(reference, critical_distance)
+
+Immutable state used by the Michelsen stability bypass. `reference` is the
+last condition at which a complete stability test was performed outside the
+shadow region, and `critical_distance` is Michelsen's smallest-eigenvalue
+measure at that condition.
+
+The storage is returned after each immutable stability or flash call. Keeping
+the reference condition unchanged across bypassed calls prevents a sequence of
+small updates from drifting across a phase boundary without a new test.
+"""
+struct StaticStabilityStorage{C, T}
+    reference::C
+    critical_distance::T
+end
+
+"""
+    StaticStabilityResult
+
+Result of [`stability_2ph_immutable`](@ref). In addition to the ordinary
+stability report and K-values, it contains the updated immutable bypass
+`storage` and records whether the full stability test was `bypassed`.
+"""
+struct StaticStabilityResult{K, S}
+    stable::Bool
+    report::StabilityReport
+    K::K
+    storage::S
+    bypassed::Bool
+end

src/flow_coupler.jl

@@ -62,6 +62,9 @@ function phase_data(mix::FlashedMixture2Phase{T, A, E}, phase) where {T, A, E}
     return out::FlashedPhase{T, A}
 end
 
+@inline phase_data(mix::FlashedMixture2Phase, ::Val{:liquid}) = mix.liquid
+@inline phase_data(mix::FlashedMixture2Phase, ::Val{:vapor}) = mix.vapor
+
 """
     phase_saturations(eos, p, T, flashed_mixture)
 

src/flow_coupler_types.jl

@@ -8,6 +8,7 @@ struct FlashedPhase{T, A<:AbstractVector{T}}
     Z::T
     function FlashedPhase(mole_fractions::AbstractVector, Z::Tz) where Tz
         T = Base.promote_type(Tz, eltype(mole_fractions))
+        mole_fractions = map(x -> convert(T, x), mole_fractions)
         Z = convert(T, Z)
         new{T, typeof(mole_fractions)}(mole_fractions, Z)
     end
@@ -82,7 +83,9 @@ end
 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, cond = cond, stability_report = stability)
+    return FlashedMixture2Phase(state, K, V, liquid, vapor,
+        vec_type = typeof(liquid.mole_fractions), critical_distance = b, cond = cond,
+        stability_report = stability)
 end
 
 function FlashedMixture2Phase(eos::AbstractEOS, T = Float64, T_num = Float64, b = NaN, cond = missing, stability = StabilityReport())

src/kvalues.jl

@@ -9,9 +9,10 @@ Estimate K-values for a given acentric factor ω and pressure and temperature at
 
 Reference [Vapor-Liquid Equilibrium. XI. A New Expression for the Excess Free Energy of Mixing by GM Wilson](https://doi.org/10.1021/ja01056a002)
 """
-function wilson_estimate(p::R, T::R, ω::R, p_c::R, T_c::R) where R<:Real
+function wilson_estimate(p, T, ω, p_c, T_c)
+    R = Base.promote_typeof(p, T, ω, p_c, T_c)
     K = exp(5.37*(1.0 + ω)*(1.0 - T_c/T))*(p_c/p)
-    return K::R
+    return convert(R, K)
 end
 
 """