From a03d805b378d76000fe5a447ac6fbbb3260112c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 12 Sep 2026 16:32:12 +0200 Subject: [PATCH 1/5] Stability stuff --- docs/src/examples/advanced.md | 21 ++++ src/MultiComponentFlash.jl | 3 +- src/flash.jl | 4 +- src/flash_types.jl | 32 ++++++ src/static.jl | 189 ++++++++++++++++++++++++++++++---- src/utils.jl | 6 ++ test/runtests.jl | 135 +++++++++++++++++++++++- test/test_setup.jl | 1 + 8 files changed, 368 insertions(+), 23 deletions(-) diff --git a/docs/src/examples/advanced.md b/docs/src/examples/advanced.md index a72d4d5..505736f 100644 --- a/docs/src/examples/advanced.md +++ b/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. diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl index 0e3b251..eb6a398 100644 --- a/src/MultiComponentFlash.jl +++ b/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 diff --git a/src/flash.jl b/src/flash.jl index 412adce..7952083 100644 --- a/src/flash.jl +++ b/src/flash.jl @@ -303,12 +303,12 @@ function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces) where {F<:Real} K[c] *= r ϵ = max(ϵ, abs(1-r)) end - V = solve_rachford_rice(K, z, V) + V = cap_unit(solve_rachford_rice(K, z, 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_unit(v) = min(max(v, zero(v)), one(v)) cap_VL(v) = min(max(v, MINIMUM_COMPOSITION), 1 - MINIMUM_COMPOSITION) function flash_update!(K, storage, type::NewtonFlash, eos, cond, forces, V, iteration) diff --git a/src/flash_types.jl b/src/flash_types.jl index 23fd6e5..622ea39 100644 --- a/src/flash_types.jl +++ b/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 diff --git a/src/static.jl b/src/static.jl index 7be38de..f04810e 100644 --- a/src/static.jl +++ b/src/static.jl @@ -22,24 +22,50 @@ Run the immutable, accelerator-friendly two-phase flash implementation. scalar vapor fraction. When `storage` is omitted, a static storage marker is created automatically. +Set `return_stability=true` to additionally return a +[`StaticStabilityResult`](@ref). Its immutable `storage` can be supplied as +`stability_storage` on the next call to enable the Michelsen bypass. The +default `bypass_tolerance=10` uses the same conservative pressure, temperature +and composition bounds as the mutable simulator integration. + The immutable path currently supports `SSIFlash` and generic cubic EOS values converted with [`make_eos_immutable`](@ref). """ -@inline function flash_2ph_immutable(eos, c; method = SSIFlash(), kwarg...) +@inline function flash_2ph_immutable(eos, c; + method = SSIFlash(), + stability_storage = nothing, + stability_bypass::Bool = !isnothing(stability_storage), + return_stability::Bool = false, + kwarg...) return flash_2ph_immutable(eos, c, flash_storage(eos, c; method = method, static = true); - method = method, kwarg...) + method = method, + stability_storage = stability_storage, + stability_bypass = stability_bypass, + return_stability = return_stability, + kwarg...) end @inline function flash_2ph_immutable(eos, c, storage::StaticConfig; - method = SSIFlash(), kwarg...) + method = SSIFlash(), + stability_storage = nothing, + stability_bypass::Bool = !isnothing(stability_storage), + return_stability::Bool = false, + kwarg...) c.z isa SVector || throw(ArgumentError( "flash_2ph_immutable requires c.z to be an SVector")) - V, K, _ = flash_2ph!(storage, initial_guess_K(eos, c, storage), eos, c, - NaN; method = method, extra_out = true, kwarg...) - return V, K + V, K, report = flash_2ph!(storage, initial_guess_K(eos, c, storage), eos, c, + NaN; method = method, extra_out = true, + stability_storage = stability_storage, + stability_bypass = stability_bypass || return_stability, + kwarg...) + return immutable_flash_output(V, K, report.stability_result, + Val(return_stability)) end +@inline immutable_flash_output(V, K, stability, ::Val{false}) = (V, K) +@inline immutable_flash_output(V, K, stability, ::Val{true}) = (V, K, stability) + """Return an isbits representation of a mixture for accelerator kernels.""" function static_mixture(mixture::MultiComponentMixture{R, N}) where {R, N} names = ntuple(_ -> nothing, Val(N)) @@ -109,6 +135,69 @@ end return (A_ij = A_ij_static, A_i = A_i_static, B_i = B_i_static) end +@inline function static_condition(c, ::Type{F}, ::Val{N}; + z_min = nothing) where {F, N} + z = SVector{N, F}(ntuple(Val(N)) do i + isnothing(z_min) ? c.z[i] : max(c.z[i], z_min) + end) + return (p = convert(F, c.p), T = convert(F, c.T), z = z) +end + +@inline invalid_stability_storage(cond) = + StaticStabilityStorage(cond, convert(typeof(cond.p), NaN)) + +@inline function stability_bypass_available(storage::StaticStabilityStorage, + cond; tolerance::Real = 10.0) + tolerance > zero(tolerance) || throw(ArgumentError( + "bypass_tolerance must be positive")) + b = storage.critical_distance + if !(isfinite(b) && b > zero(b)) + return false + end + reference = storage.reference + return maximum(abs, reference.z - cond.z) < b/tolerance && + abs(reference.p - cond.p) < b*abs(cond.p)/tolerance && + abs(reference.T - cond.T) < b*tolerance +end + +@inline static_minimum_eigenvalue(B::SMatrix) = + minimum(eigvals(Symmetric(B))) + +"""Immutable Michelsen critical-point distance used by the stability bypass.""" +@generated function static_fugacity_coefficients( + eos::GenericCubicEOS{E, R, N}, cond, Z, forces, scalars, + ::Type{D}) where {E, R, N, D} + values = [:(component_fugacity_coefficient( + eos, cond, $i, Z, forces, scalars)) for i in 1:N] + return :(SVector{N, D}(($(values...),))) +end + +@inline function static_michelsen_critical_point_measure( + eos::GenericCubicEOS{E, R, N}, p, temperature, + mole_numbers::SVector{N, F}) where {E, R, N, F} + D = ForwardDiff.Dual{Nothing, F, N} + mole_numbers_ad = SVector{N, D}(ntuple(Val(N)) do i + partials = ForwardDiff.single_seed(ForwardDiff.Partials{N, F}, Val(i)) + D(mole_numbers[i], partials) + end) + z = mole_numbers_ad/sum(mole_numbers_ad) + cond = (p = convert(F, p), T = convert(F, temperature), z = z) + forces = static_force_coefficients(eos, cond, D) + scalars = force_scalars(eos, cond, forces) + Z = mixture_compressibility_factor(eos, cond, forces, scalars) + coefficients = static_fugacity_coefficients( + eos, cond, Z, forces, scalars, D) + B = SMatrix{N, N, F}(ntuple(Val(N*N)) do index + i = mod1(index, N) + j = (index - 1) ÷ N + 1 + F(i == j) + sqrt(mole_numbers[i]*mole_numbers[j])*coefficients[i].partials[j] + end) + # Roundoff in the AD construction can make the theoretically symmetric + # matrix differ by a few ulps. Symmetrize before finding its eigenvalues. + B = (B + transpose(B))/2 + return static_minimum_eigenvalue(B) +end + @inline function solve_rachford_rice(K::StaticVector{2}, z::StaticVector{2}, V = NaN) z1, z2 = z k1, k2 = K @@ -203,7 +292,7 @@ end residual = max(residual, abs(one(F) - ratios[i])) end K_next = SVector{N, F}(ntuple(i -> K[i]*ratios[i], Val(N))) - V_next = solve_rachford_rice(K_next, z, V) + V_next = cap_unit(solve_rachford_rice(K_next, z, V)) return V_next, K_next, residual end @@ -230,23 +319,33 @@ end check::Bool = true, update_forces::Bool = true, z_min = MINIMUM_COMPOSITION, + stability_storage = nothing, + stability_bypass::Bool = !isnothing(stability_storage), + bypass_tolerance::Real = 10.0, kwarg... ) where {E, R, N} F = Base.promote_eltype(c.p, c.T, c.z[1], K[1]) - z = SVector{N, F}(ntuple(Val(N)) do i - isnothing(z_min) ? c.z[i] : max(c.z[i], z_min) - end) K = SVector{N, F}(K) - cond = (p = convert(F, c.p), T = convert(F, c.T), z = z) + cond = static_condition(c, F, Val(N); z_min = z_min) + z = cond.z forces = static_force_coefficients(eos, cond, F) V = convert(F, V) single_phase_init = isnan(V) || V == one(F) || V == zero(F) if single_phase_init - stable, stability_report, K = static_stability_2ph( - K, eos, cond, forces; maxiter = maxiter, kwarg...) + stability_result = static_stability_2ph(K, eos, cond, forces; + storage = stability_storage_value(stability_storage), + update_bypass = stability_bypass, + bypass_tolerance = bypass_tolerance, + maxiter = maxiter, + kwarg...) + stable = stability_result.stable + stability_report = stability_result.report + K = stability_result.K else stable = false stability_report = StabilityReport(false, false, false, false) + stability_result = StaticStabilityResult(stable, stability_report, K, + invalid_stability_storage(cond), false) end converged = false if stable @@ -263,7 +362,8 @@ end iteration += 1 end end - report = (its = iteration, converged = converged, stability = stability_report) + report = (its = iteration, converged = converged, + stability = stability_report, stability_result = stability_result) return V, K, report end @@ -342,8 +442,17 @@ end @inline function static_stability_2ph(K::SVector{N, F}, eos, cond, forces; check_vapor::Bool = true, check_liquid::Bool = true, + storage = nothing, + update_bypass::Bool = false, + bypass_tolerance::Real = 10.0, kwarg... ) where {N, F} + if update_bypass && !isnothing(storage) && + stability_bypass_available(storage, cond; + tolerance = bypass_tolerance) + report = StabilityReport(true, true, true, true) + return StaticStabilityResult(true, report, K, storage, true) + end vapor_phase = (p = cond.p, T = cond.T, z = cond.z, phase = Val(:vapor)) f_z_vapor = static_fugacities(eos, vapor_phase, forces, F) K_wilson = initial_guess_K(eos, cond, StaticConfig()) @@ -365,19 +474,61 @@ end report = StabilityReport(stable_liquid, trivial_liquid, stable_vapor, trivial_vapor) K_out = report.stable ? K_liquid : static_divide(y, x) - return report.stable, report, K_out + if update_bypass && report.stable && report.liquid.trivial && + report.vapor.trivial + critical_distance = static_michelsen_critical_point_measure( + eos, cond.p, cond.T, cond.z) + next_storage = StaticStabilityStorage(cond, critical_distance) + else + next_storage = invalid_stability_storage(cond) + end + return StaticStabilityResult(report.stable, report, K_out, + next_storage, false) end @inline function stability_2ph(eos::GenericCubicEOS{E, R, N}, c, K, config::StaticConfig; extra_out::Bool = false, kwarg...) where {E, R, N} F = Base.promote_eltype(c.p, c.T, c.z[1], K[1]) - cond = (p = convert(F, c.p), T = convert(F, c.T), - z = SVector{N, F}(c.z)) + cond = static_condition(c, F, Val(N)) K = SVector{N, F}(K) forces = static_force_coefficients(eos, cond, F) - stable, report, _ = static_stability_2ph(K, eos, cond, forces; kwarg...) - return extra_out ? (stable, report) : stable + result = static_stability_2ph(K, eos, cond, forces; kwarg...) + return extra_out ? (result.stable, result.report) : result.stable end @inline stability_2ph!(::StaticConfig, K, eos::GenericCubicEOS, c; kwarg...) = stability_2ph(eos, c, K, StaticConfig(); kwarg...) + +""" + result = stability_2ph_immutable(eos, c[, storage]; ) + +Run the immutable stability test independently of a flash. The returned +[`StaticStabilityResult`](@ref) contains the stability report, K-values and +updated [`StaticStabilityStorage`](@ref). Pass `result.storage` to a later call +to enable Michelsen's stability bypass for nearby conditions. + +The bypass is only armed after both trial phases converge to trivial stable +solutions. Calls inside the shadow region retain the full stability test. +""" +@inline function stability_2ph_immutable(eos::GenericCubicEOS{E, R, N}, c, + storage = nothing; + K = initial_guess_K(eos, c, StaticConfig()), + bypass_tolerance::Real = 10.0, + z_min = MINIMUM_COMPOSITION, + kwarg...) where {E, R, N} + c.z isa SVector || throw(ArgumentError( + "stability_2ph_immutable requires c.z to be an SVector")) + F = Base.promote_eltype(c.p, c.T, c.z[1], K[1]) + cond = static_condition(c, F, Val(N); z_min = z_min) + K = SVector{N, F}(K) + forces = static_force_coefficients(eos, cond, F) + return static_stability_2ph(K, eos, cond, forces; + storage = stability_storage_value(storage), + update_bypass = true, + bypass_tolerance = bypass_tolerance, + kwarg...) +end + +@inline stability_storage_value(::Nothing) = nothing +@inline stability_storage_value(storage::StaticStabilityStorage) = storage +@inline stability_storage_value(result::StaticStabilityResult) = result.storage diff --git a/src/utils.jl b/src/utils.jl index c387497..7d22e9f 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -134,6 +134,12 @@ function michelsen_critical_point_measure(eos, p, T, mole_numbers; kwarg...) michelsen_critical_point_measure!(S, eos, p, T, mole_numbers) end +@inline function michelsen_critical_point_measure( + eos::GenericCubicEOS{E, R, N}, p, T, + mole_numbers::SVector{N, F}; kwarg...) where {E, R, N, F} + return static_michelsen_critical_point_measure(eos, p, T, mole_numbers) +end + function michelsen_critical_point_measure!(S, eos, p, T, mole_numbers) # From Michelsen (1982): The Isothermal Flash Problem. Part I: Stability # Estimate the distance to the critical point through smallest eigenvalue. diff --git a/test/runtests.jl b/test/runtests.jl index 1e17c26..b68797b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,7 +42,8 @@ end end @testset "Static accelerator path" begin - eos = make_eos_immutable(get_test_eos()) + host_eos = get_test_eos() + eos = make_eos_immutable(host_eos) c = (p = 1e6, T = 300.0, z = @SVector [0.5, 0.3, 0.2]) storage = flash_storage(eos, c; method = SSIFlash(), static = true) K = initial_guess_K(eos, c, storage) @@ -79,6 +80,102 @@ end @test flash_2ph_immutable(eos, c, storage) == (V_immutable, K_immutable) @test_throws ArgumentError flash_2ph_immutable(eos, (p = c.p, T = c.T, z = collect(c.z))) + + @testset "Standalone stability and bypass" begin + stable_cond = (p = 1e5, T = 800.0, z = c.z) + stability = stability_2ph_immutable(eos, stable_cond) + @test stability isa MultiComponentFlash.StaticStabilityResult + @test stability.storage isa MultiComponentFlash.StaticStabilityStorage + @test isbitstype(typeof(stability)) + @test stability.stable + @test stability.report.liquid.trivial + @test stability.report.vapor.trivial + @test isfinite(stability.storage.critical_distance) + @test !stability.bypassed + @test test_static_stability_allocs(eos, stable_cond) == 0 + + nearby = (p = 1.001e5, T = 800.01, + z = @SVector [0.50001, 0.29999, 0.2]) + bypassed = stability_2ph_immutable(eos, nearby, stability) + @test bypassed.stable + @test bypassed.bypassed + @test bypassed.storage == stability.storage + + far_away = (p = 1e6, T = 800.0, z = c.z) + retested = stability_2ph_immutable(eos, far_away, stability.storage) + @test !retested.bypassed + @test retested.storage.reference == far_away + + # A stable state inside the shadow region is deliberately not armed. + shadow = stability_2ph_immutable(eos, + (p = 1e5, T = 500.0, z = c.z)) + @test shadow.stable + @test !shadow.report.liquid.trivial + @test isnan(shadow.storage.critical_distance) + + V_stable, K_stable, flash_stability = flash_2ph_immutable( + eos, stable_cond; return_stability = true) + @test isnan(V_stable) + @test all(isfinite, K_stable) + @test flash_stability.stable + V_nearby, K_nearby, nearby_stability = flash_2ph_immutable( + eos, nearby; + stability_storage = flash_stability, + return_stability = true) + @test isnan(V_nearby) + @test all(isfinite, K_nearby) + @test nearby_stability.bypassed + @test_throws ArgumentError stability_2ph_immutable(eos, nearby, + stability.storage; bypass_tolerance = 0.0) + end + + @testset "Nearly absent components" begin + tiny = 1e-16 + compositions = ( + SVector(tiny, 0.3, 0.7 - tiny), + SVector(0.5, tiny, 0.5 - tiny), + SVector(0.3, 0.7 - tiny, tiny), + SVector(tiny, tiny, 1.0 - 2tiny) + ) + for z_tiny in compositions + cond_tiny = (p = 1e5, T = 250.0, z = z_tiny) + result = stability_2ph_immutable(eos, cond_tiny) + @test isbitstype(typeof(result)) + @test all(isfinite, result.K) + + V_tiny, K_tiny, report_tiny = flash_2ph!(storage, + initial_guess_K(eos, cond_tiny, storage), eos, cond_tiny, NaN; + extra_out = true) + @test report_tiny.converged || result.stable + @test isnan(V_tiny) || 0.0 <= V_tiny <= 1.0 + @test all(isfinite, K_tiny) + + host_cond_tiny = (p = cond_tiny.p, T = cond_tiny.T, + z = collect(z_tiny)) + V_host, K_host, host_report = flash_2ph(host_eos, + host_cond_tiny; extra_out = true) + @test host_report.converged || host_report.stability.stable + @test isnan(V_host) || 0.0 <= V_host <= 1.0 + @test all(isfinite, K_host) + @test isapprox(V_tiny, V_host; nans = true, atol = 1e-12) + + # Exercise the AD/eigenvalue part of the bypass without flooring + # away the trace component. + distance = MultiComponentFlash.michelsen_critical_point_measure( + eos, cond_tiny.p, cond_tiny.T, z_tiny) + @test isfinite(distance) + + stable_tiny = (p = 5e7, T = 800.0, z = z_tiny) + tiny_stability = stability_2ph_immutable(eos, stable_tiny; + z_min = nothing) + @test tiny_stability.stable + @test isfinite(tiny_stability.storage.critical_distance) + tiny_nearby = (p = 5.001e7, T = 800.01, z = z_tiny) + tiny_bypass = stability_2ph_immutable(eos, tiny_nearby, + tiny_stability; z_min = nothing) + @test tiny_bypass.bypassed + end + end end @testset "Partial derivatives" begin @@ -146,6 +243,21 @@ end z = [0.4, 0.6] @test MultiComponentFlash.michelsen_critical_point_measure(equation_of_state, 5e6, 303.15, z) ≈ 0.776435 atol = 1e-4 @test MultiComponentFlash.michelsen_critical_point_measure(equation_of_state, 5e6, 303.15, z, static_size = false) ≈ 0.776435 atol = 1e-4 + equation_of_state_static = make_eos_immutable(equation_of_state) + @test MultiComponentFlash.michelsen_critical_point_measure( + equation_of_state_static, 5e6, 303.15, SVector{2}(z)) ≈ 0.776435 atol = 1e-4 + + ethane = MolecularProperty(0.03007, 4.872e6, 305.32, 1.455e-4, 0.099) + carbon_dioxide = MolecularProperty(0.0440, 7.38e6, 304.1, 9.412e-5, 0.224) + mixture_4 = MultiComponentMixture((methane, ethane, carbon_dioxide, decane)) + eos_4 = GenericCubicEOS(mixture_4, PengRobinson()) + eos_4_static = make_eos_immutable(eos_4) + z_4 = @SVector [1e-16, 0.2, 0.3, 0.5 - 1e-16] + dynamic_distance = MultiComponentFlash.michelsen_critical_point_measure( + eos_4, 5e7, 700.0, collect(z_4); static_size = false) + static_distance = MultiComponentFlash.michelsen_critical_point_measure( + eos_4_static, 5e7, 700.0, z_4) + @test static_distance ≈ dynamic_distance rtol = 1e-12 end using StaticArrays, KernelAbstractions, JLArrays @@ -160,6 +272,19 @@ using StaticArrays, KernelAbstractions, JLArrays @inbounds out[i] = V end end + @kernel function static_bypass_kernel!(distance, bypassed, eos, z) + i = @index(Global) + if i <= length(distance) + initial = (p = 5e7, T = 800.0, z = z) + stability = stability_2ph_immutable(eos, initial; + z_min = nothing) + nearby = (p = 5.001e7, T = 800.01, z = z) + next_stability = stability_2ph_immutable(eos, nearby, + stability; z_min = nothing) + @inbounds distance[i] = stability.storage.critical_distance + @inbounds bypassed[i] = next_stability.bypassed + end + end if isdefined(JLArrays, :JLBackend) host_eos = get_test_eos() eos = make_eos_immutable(host_eos) @@ -182,6 +307,14 @@ using StaticArrays, KernelAbstractions, JLArrays kernel!(out, pressure, temperature, z, eos, storage; ndrange = n) @test Array(out) ≈ expected rtol = 1e-11 + + trace_z = @SVector [1e-16, 0.3, 0.7 - 1e-16] + distance = JLArray(zeros(1)) + bypassed = JLArray(falses(1)) + bypass_kernel! = static_bypass_kernel!(backend, 1) + bypass_kernel!(distance, bypassed, eos, trace_z; ndrange = 1) + @test isfinite(only(Array(distance))) + @test only(Array(bypassed)) else # JLArrays 0.1 supports Julia 1.6 but predates the KernelAbstractions backend. @test_skip false diff --git a/test/test_setup.jl b/test/test_setup.jl index 0ec3b42..4fbb2f4 100644 --- a/test/test_setup.jl +++ b/test/test_setup.jl @@ -12,6 +12,7 @@ end test_conditions() = (p = 10e5, T = 300.0, z = [0.5, 0.3, 0.2]) test_allocs(S, K, eos, c, m) = @allocated flash_2ph!(S, K, eos, c, NaN, method = m) +test_static_stability_allocs(eos, c) = @allocated stability_2ph_immutable(eos, c) function test_flash_inplace(m, do_test = true) eos = get_test_eos() From 75c887ff54875d705289369bdc1c40f18e82207e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 12 Sep 2026 21:27:41 +0200 Subject: [PATCH 2/5] Compositional KA --- src/eos_types.jl | 4 ++-- src/flow_coupler.jl | 3 +++ src/flow_coupler_types.jl | 4 +++- src/kvalues.jl | 5 +++-- src/static.jl | 25 ++++++++++++++++++++----- test/runtests.jl | 17 +++++++++++++++++ 6 files changed, 48 insertions(+), 10 deletions(-) diff --git a/src/eos_types.jl b/src/eos_types.jl index 55f57ef..84b095c 100644 --- a/src/eos_types.jl +++ b/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 diff --git a/src/flow_coupler.jl b/src/flow_coupler.jl index 68e8f12..9b7e1bb 100644 --- a/src/flow_coupler.jl +++ b/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) diff --git a/src/flow_coupler_types.jl b/src/flow_coupler_types.jl index 1658f21..01e2f8a 100644 --- a/src/flow_coupler_types.jl +++ b/src/flow_coupler_types.jl @@ -82,7 +82,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(x), critical_distance = b, cond = cond, + stability_report = stability) end function FlashedMixture2Phase(eos::AbstractEOS, T = Float64, T_num = Float64, b = NaN, cond = missing, stability = StabilityReport()) diff --git a/src/kvalues.jl b/src/kvalues.jl index 8d164cd..3c3e650 100644 --- a/src/kvalues.jl +++ b/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 """ diff --git a/src/static.jl b/src/static.jl index f04810e..86054ed 100644 --- a/src/static.jl +++ b/src/static.jl @@ -135,6 +135,19 @@ end return (A_ij = A_ij_static, A_i = A_i_static, B_i = B_i_static) end +function make_eos_immutable(eos::KValuesEOS{T, R, N}) where {T, R, N} + mixture = static_mixture(eos.mixture) + evaluator = eos.K_values_evaluator + if evaluator isa AbstractVector + evaluator = SVector{N, eltype(evaluator)}(evaluator) + end + volume_shift = eos.volume_shift + if !isnothing(volume_shift) + volume_shift = SVector{N, eltype(volume_shift)}(volume_shift) + end + return KValuesEOS(evaluator, mixture; volume_shift = volume_shift) +end + @inline function static_condition(c, ::Type{F}, ::Val{N}; z_min = nothing) where {F, N} z = SVector{N, F}(ntuple(Val(N)) do i @@ -270,12 +283,14 @@ end return V end -@inline function static_fugacities(eos::GenericCubicEOS{E, R, N}, cond, forces, +@generated function static_fugacities(eos::GenericCubicEOS{E, R, N}, cond, forces, ::Type{F}) where {E, R, N, F} - Z, scalars = prep(eos, cond, forces) - return SVector{N, F}(ntuple(Val(N)) do component - component_fugacity(eos, cond, component, Z, forces, scalars) - end) + values = [:(component_fugacity( + eos, cond, $i, Z, forces, scalars)) for i in 1:N] + return quote + Z, scalars = prep(eos, cond, forces) + SVector{N, F}(($(values...),)) + end end @inline function static_ssi(K::SVector{N, F}, p::F, T::F, z, V::F, diff --git a/test/runtests.jl b/test/runtests.jl index b68797b..5bf0b5d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -196,6 +196,23 @@ end @test number_of_components(eos) == 2 eos2 = KValuesEOS(cond -> [0.01, 2.0], mixture) @test round(flash_2ph(eos2, cond), digits = 4) ≈ 0.8091 + + static_eos = make_eos_immutable(eos) + static_cond = (p = cond.p, T = cond.T, z = @SVector [0.1, 0.9]) + @test isbitstype(typeof(static_eos)) + @test static_eos.K_values_evaluator isa SVector{2, Float64} + @test round(flash_2ph(static_eos, static_cond), digits = 4) ≈ 0.8091 +end + +@testset "Static flashed mixture storage" begin + x = @SVector [0.8, 0.2] + y = @SVector [0.1, 0.9] + flashed = FlashedMixture2Phase( + MultiComponentFlash.two_phase_lv, SVector(0.125, 4.5), + 0.4, x, y, 0.9, 1.1) + @test isbitstype(typeof(flashed)) + @test phase_data(flashed, Val(:liquid)).mole_fractions === x + @test phase_data(flashed, Val(:vapor)).mole_fractions === y end using ForwardDiff From f570a590ee51607dedcf087db776bb96a6d51a37 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 12 Sep 2026 22:49:33 +0200 Subject: [PATCH 3/5] Fix promotion issue --- src/flow_coupler_types.jl | 3 ++- test/runtests.jl | 15 +++++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/src/flow_coupler_types.jl b/src/flow_coupler_types.jl index 01e2f8a..aeb66f1 100644 --- a/src/flow_coupler_types.jl +++ b/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 @@ -83,7 +84,7 @@ function FlashedMixture2Phase(state, K, V, x, y, Z_L, Z_V, b = NaN, cond = missi liquid = FlashedPhase(x, Z_L) vapor = FlashedPhase(y, Z_V) return FlashedMixture2Phase(state, K, V, liquid, vapor, - vec_type = typeof(x), critical_distance = b, cond = cond, + vec_type = typeof(liquid.mole_fractions), critical_distance = b, cond = cond, stability_report = stability) end diff --git a/test/runtests.jl b/test/runtests.jl index 5bf0b5d..cb13454 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -216,6 +216,21 @@ end end using ForwardDiff +@testset "Static flashed mixture promotion" begin + x = @SVector [0.8, 0.2] + y = @SVector [0.1, 0.9] + dZ = ForwardDiff.derivative(0.4) do V + flashed = FlashedMixture2Phase( + MultiComponentFlash.two_phase_lv, SVector(0.125, 4.5), + V, x, y, V + 0.5, V + 0.7) + @test eltype(flashed.liquid.mole_fractions) === typeof(V) + @test eltype(flashed.vapor.mole_fractions) === typeof(V) + @test isbitstype(typeof(flashed)) + return flashed.liquid.Z + end + @test dZ == 1.0 +end + @testset "Rachford-Rice derivatives" begin N = 25 for z_light in range(0.0, 1.0, length = N) From 85dadf72d5cd3d5ab97652e89ee0227cf1707d9f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 13 Sep 2026 11:39:58 +0200 Subject: [PATCH 4/5] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7623787..37e4d3a 100644 --- a/Project.toml +++ b/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 "] [deps] From ec9e483a5ba8679aa7e439a4822830d417d3d79b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 13 Sep 2026 12:05:57 +0200 Subject: [PATCH 5/5] Deprecate --- src/flash.jl | 11 +++++++---- src/static.jl | 3 ++- test/runtests.jl | 5 ++++- 3 files changed, 13 insertions(+), 6 deletions(-) diff --git a/src/flash.jl b/src/flash.jl index 7952083..f9c98c9 100644 --- a/src/flash.jl +++ b/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 @@ -303,12 +306,12 @@ function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces) where {F<:Real} K[c] *= r ϵ = max(ϵ, abs(1-r)) end - V = cap_unit(solve_rachford_rice(K, z, V)) + 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(v)), one(v)) cap_VL(v) = min(max(v, MINIMUM_COMPOSITION), 1 - MINIMUM_COMPOSITION) function flash_update!(K, storage, type::NewtonFlash, eos, cond, forces, V, iteration) diff --git a/src/static.jl b/src/static.jl index 728540c..f3055da 100644 --- a/src/static.jl +++ b/src/static.jl @@ -305,7 +305,8 @@ end residual = max(residual, abs(one(F) - ratios[i])) end K_next = SVector{N, F}(ntuple(i -> K[i]*ratios[i], Val(N))) - V_next = cap_unit(solve_rachford_rice(K_next, z, V)) + V_next = solve_rachford_rice(K_next, z, V) + V_next = clamp(V_next, zero(V_next), one(V_next)) return V_next, K_next, residual end diff --git a/test/runtests.jl b/test/runtests.jl index be2b480..0c5c80e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -61,7 +61,10 @@ end normal_config = MultiComponentFlash.FlashConfig(print_output=false) @test typeof(normal_config) == MultiComponentFlash.FlashConfig @test flash_storage(eos, c, SSIFlash(), normal_config).x isa Vector - @test_throws ArgumentError flash_storage(eos, c; static_size = true) + deprecated_storage = @test_deprecated flash_storage(eos, c; static_size = true) + @test deprecated_storage isa MultiComponentFlash.StaticConfig + dynamic_storage = @test_deprecated flash_storage(eos, c; static_size = false) + @test dynamic_storage.x isa Vector config = MultiComponentFlash.FlashConfig(print_output=false, use_dict_storage=false) @test typeof(config) == MultiComponentFlash.FlashConfig