From 7ef9a5484a0758517f1a5ae36bce1531436e3723 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Fri, 11 Sep 2026 10:36:28 +0200
Subject: [PATCH 1/9] Draft of GPU support for flash

---
 src/MultiComponentFlash.jl |   5 +-
 src/eos.jl                 | 107 +++++++++++++++++++++++---------
 src/eos_types.jl           |  28 ++++++++-
 src/flash.jl               | 121 ++++++++++++++++++++++++++++++++++---
 src/kvalues.jl             |  13 +++-
 src/mixture_types.jl       |  25 ++++++--
 src/rachford_rice.jl       |  73 ++++++++++++++++++++++
 src/utils.jl               |  20 ++++--
 test/runtests.jl           |  21 +++++++
 9 files changed, 358 insertions(+), 55 deletions(-)

diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl
index 6514b7e..0d99ff1 100644
--- a/src/MultiComponentFlash.jl
+++ b/src/MultiComponentFlash.jl
@@ -19,14 +19,14 @@ module MultiComponentFlash
     export KValuesEOS
     export number_of_components
     # Flash interfaces
-    export flash_2ph, flash_2ph!, flash_storage
+    export flash_2ph, flash_2ph!, flash_2ph_static, flash_storage, flash_storage_static
     export stability_2ph, stability_2ph!
     # Algorithms for flash
     export SSIFlash, NewtonFlash, SSINewtonFlash
     # Mixtures and their molecular makeup
     export MolecularProperty, MultiComponentMixture
     # K-values
-    export wilson_estimate, wilson_estimate!, initial_guess_K, initial_guess_K!
+    export wilson_estimate, wilson_estimate!, initial_guess_K, initial_guess_K!, initial_guess_K_static
     # Vapor-liquid equilibrium
     export solve_rachford_rice
 
@@ -35,6 +35,7 @@ module MultiComponentFlash
 
     export force_scalars, force_coefficients, force_coefficients!
     export critical_pressure, critical_temperature, critical_volume, acentric_factor, molar_weight
+    export static_eos
     export cubic_benchmark
     export single_phase_label
 
diff --git a/src/eos.jl b/src/eos.jl
index c4ceeed..303949b 100644
--- a/src/eos.jl
+++ b/src/eos.jl
@@ -11,9 +11,12 @@ number_of_components(e::AbstractEOS) = number_of_components(e.mixture)
 forces_per_phase(eos::GenericCubicEOS) = false
 
 function get_phase(cond)
-    return get(cond, :phase, :unknown)::Symbol
+    return phase_symbol(get(cond, :phase, :unknown))
 end
 
+@inline phase_symbol(phase::Symbol) = phase
+@inline phase_symbol(::Val{phase}) where phase = phase
+
 function set_phase(cond, phase::Symbol, throw::Bool = false)
     if throw && haskey(cond, :phase) && cond.phase != :unknown
         throw(ArgumentError("Phase state already set to $(cond.phase), cannot change to $phase."))
@@ -93,34 +96,50 @@ minimum_allowable_root(eos, forces, scalars) = 1e-16
     return roots
 end
 
-function pick_root(eos, roots, cond, forces, scalars)
-    phase = get_phase(cond)
-    r_ϵ = minimum_allowable_root(eos, forces, scalars)
-    max_r = maximum(roots)
-    min_r = minimum((x) -> x > r_ϵ ? x : Inf, roots)
-    if min_r == max_r
-        r = min_r
-    elseif phase == :liquid
-        r = min_r
-    elseif phase == :vapor
-        r = max_r
-    else
-        function Gibbs(Z)
-            E = 0.0
-            z = cond.z
-            @inbounds for i in eachindex(z)
-                ϕ = component_fugacity_coefficient(eos, cond, i, Z, forces, scalars)
-                E += z[i]*ϕ
-            end
-            return E
+@inline function root_bounds(roots, minimum_root)
+    max_root = -Inf
+    min_root = Inf
+    for root in roots
+        max_root = max(max_root, root)
+        if root > minimum_root
+            min_root = min(min_root, root)
         end
-        if Gibbs(min_r) < Gibbs(max_r)
-            r = min_r
-        else
-            r = max_r
+    end
+    return min_root, max_root
+end
+
+@inline function pick_root(eos, roots, cond, forces, scalars)
+    phase = get(cond, :phase, :unknown)
+    return pick_root(eos, roots, cond, forces, scalars, phase)
+end
+
+@inline function pick_root(eos, roots, cond, forces, scalars, ::Val{:liquid})
+    min_r, _ = root_bounds(roots, minimum_allowable_root(eos, forces, scalars))
+    return min_r
+end
+
+@inline function pick_root(eos, roots, cond, forces, scalars, ::Val{:vapor})
+    _, max_r = root_bounds(roots, minimum_allowable_root(eos, forces, scalars))
+    return max_r
+end
+
+function pick_root(eos, roots, cond, forces, scalars, phase::Symbol)
+    min_r, max_r = root_bounds(roots, minimum_allowable_root(eos, forces, scalars))
+    if min_r == max_r || phase == :liquid
+        return min_r
+    elseif phase == :vapor
+        return max_r
+    end
+    function Gibbs(Z)
+        E = 0.0
+        z = cond.z
+        @inbounds for i in eachindex(z)
+            ϕ = component_fugacity_coefficient(eos, cond, i, Z, forces, scalars)
+            E += z[i]*ϕ
         end
+        return E
     end
-    return r
+    return Gibbs(min_r) < Gibbs(max_r) ? min_r : max_r
 end
 
 """
@@ -138,9 +157,9 @@ function force_coefficients(eos::AbstractCubicEOS, cond; static_size = false)
     n = number_of_components(eos)
     eT = Base.promote_eltype(cond.p, cond.T, cond.z[1])
     if static_size
-        A_ij = @MMatrix zeros(eT, n, n)
-        A_i = @MVector zeros(eT, n)
-        B_i = @MVector zeros(eT, n)
+        A_ij = zero(MMatrix{n, n, eT})
+        A_i = zero(MVector{n, eT})
+        B_i = zero(MVector{n, eT})
     else
         A_ij = zeros(eT, n, n)
         A_i = zeros(eT, n)
@@ -151,6 +170,34 @@ function force_coefficients(eos::AbstractCubicEOS, cond; static_size = false)
     return coeff
 end
 
+"""Construct inline force coefficients with the component count from the EOS type."""
+function force_coefficients_static(eos::GenericCubicEOS{E, R, N}, cond) where {E, R, N}
+    T = Base.promote_eltype(cond.p, cond.T, cond.z[1])
+    return force_coefficients_static(eos, cond, T)
+end
+
+function force_coefficients_static(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T}) where {E, R, N, T}
+    coeff = (
+        A_ij = zero(MMatrix{N, N, T}),
+        A_i = zero(MVector{N, T}),
+        B_i = zero(MVector{N, T})
+    )
+    return update_force_coefficients!(coeff, eos, cond)
+end
+
+"""Immutable force coefficients for register-oriented accelerator kernels."""
+@inline function force_coefficients_stack(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T}) where {E, R, N, T}
+    A_i_static = SVector{N, T}(ntuple(i -> A_i(eos, cond, i), Val(N)))
+    B_i_static = SVector{N, T}(ntuple(i -> B_i(eos, cond, i), Val(N)))
+    A_ij_static = SMatrix{N, N, T}(ntuple(Val(N*N)) do index
+        i = mod1(index, N)
+        j = (index - 1) ÷ N + 1
+        sqrt(A_i_static[i]*A_i_static[j]) *
+            (one(T) - binary_interaction(eos, i, j, cond))
+    end)
+    return (A_ij = A_ij_static, A_i = A_i_static, B_i = B_i_static)
+end
+
 function get_force_coefficients(forces, eos, cond)
     if forces_per_phase(eos)
         phase = get_phase(cond)
@@ -166,6 +213,8 @@ function get_force_coefficients(forces, eos, cond)
     end
 end
 
+@inline get_force_coefficients(forces, eos::GenericCubicEOS, cond) = forces
+
 """
     force_coefficients!(coeff, eos, cond)
 
diff --git a/src/eos_types.jl b/src/eos_types.jl
index 748ada2..64517d3 100644
--- a/src/eos_types.jl
+++ b/src/eos_types.jl
@@ -11,9 +11,9 @@ definitions for the terms (they are, after all, all cubic in form). References:
  2. [Simulation of Gas Condensate Reservoir Performance  by K.H. Coats](https://doi.org/10.2118/10512-PA)
 
 """
-struct GenericCubicEOS{T, R, N, V} <: AbstractCubicEOS
+struct GenericCubicEOS{T, R, N, V, M<:MultiComponentMixture{R, N}} <: AbstractCubicEOS
     type::T
-    mixture::MultiComponentMixture{R, N}
+    mixture::M
     m_1::R
     m_2::R
     ω_a::R
@@ -182,6 +182,30 @@ 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
 
+"""
+    static_eos(eos)
+
+Convert a generic cubic EOS to an isbits representation suitable for passing to
+accelerator kernels. The numerical model is unchanged; only descriptive metadata
+and array storage are converted.
+"""
+function static_eos(eos::GenericCubicEOS{T, R, N}) where {T, R, N}
+    mixture = static_mixture(eos.mixture)
+    volume_shift = eos.volume_shift
+    if !isnothing(volume_shift)
+        volume_shift = SVector{N, eltype(volume_shift)}(volume_shift)
+    end
+    return GenericCubicEOS(
+        eos.type,
+        mixture,
+        eos.m_1,
+        eos.m_2,
+        eos.ω_a,
+        eos.ω_b,
+        volume_shift
+    )
+end
+
 struct KValuesEOS{T, R, N, V} <: AbstractEOS
     "Callable on the form `cond -> V` or a set of constants (Tuple/AbstractVector)"
     K_values_evaluator::T
diff --git a/src/flash.jl b/src/flash.jl
index 02fabd6..2a86952 100644
--- a/src/flash.jl
+++ b/src/flash.jl
@@ -159,6 +159,71 @@ function flash_storage(eos, cond = (p = 10e5, T = 273.15, z = zeros(number_of_co
     return NamedTuple(pairs(d))
 end
 
+"""
+    flash_storage_static(eos, cond, [method])
+
+Construct an inline, statically sized SSI workspace without using the `Dict`-based
+host storage builder. This constructor is suitable for use inside accelerator
+kernels when `eos` and `cond` are isbits values.
+"""
+function flash_storage_static(eos::GenericCubicEOS{E, R, N}, cond, method::SSIFlash = SSIFlash()) where {E, R, N}
+    T = Base.promote_eltype(cond.p, cond.T, cond.z[1])
+    forces = force_coefficients_static(eos, cond)
+    return (
+        forces = forces,
+        x = zero(MVector{N, T}),
+        y = zero(MVector{N, T}),
+        buffer1 = zero(MVector{N, T}),
+        buffer2 = zero(MVector{N, T})
+    )
+end
+
+
+"""
+    flash_2ph_static(eos, cond, K, V[, maxiter, tolerance])
+
+Accelerator-compatible SSI equilibrium solve with entirely local static storage.
+`K` must be a static vector and `V` must be a finite initial vapor
+fraction. Stability testing and host-side diagnostics are intentionally excluded.
+Returns `(V, K, iterations, converged)`.
+"""
+@inline flash_2ph_static(eos, cond, K, V) =
+    flash_2ph_static(eos, cond, K, V, 25000, 1e-8)
+
+@inline function flash_2ph_static(
+        eos::GenericCubicEOS{E, R, N}, cond, K::MVector{N, T}, V,
+        maxiter::Int, tolerance
+    ) where {E, R, N, T}
+    forces = force_coefficients_static(eos, cond, T)
+    x = zero(MVector{N, T})
+    y = zero(MVector{N, T})
+    converged = false
+    iteration = 0
+    while iteration < maxiter
+        iteration += 1
+        V, residual = ssi_static!(K, cond.p, cond.T, x, y, cond.z, V, eos, forces)
+        converged = residual <= tolerance
+        converged && break
+    end
+    return (V, K, iteration, converged)
+end
+
+@inline function flash_2ph_static(
+        eos::GenericCubicEOS{E, R, N}, cond, K::SVector{N, T}, V,
+        maxiter::Int, tolerance
+    ) where {E, R, N, T}
+    forces = force_coefficients_stack(eos, cond, T)
+    converged = false
+    iteration = 0
+    while iteration < maxiter
+        iteration += 1
+        V, K, residual = ssi_static(K, cond.p, cond.T, cond.z, V, eos, forces)
+        converged = residual <= tolerance
+        converged && break
+    end
+    return (V, K, iteration, converged)
+end
+
 function flash_storage_internal!(out, eos, cond, method;
         inc_jac = isa(method, AbstractNewtonFlash),
         inc_bypass = false,
@@ -177,7 +242,12 @@ function flash_storage_internal!(out, eos, cond, method;
         out[:forces] = alloc_forces(cond)
     end
     if static_size
-        alloc_vec = () -> @MVector zeros(n)
+        # Do not write `@MVector zeros(n)` here.  That form first creates a
+        # regular Vector and then copies it into the MVector.  Apart from the
+        # unnecessary allocation, this is particularly inconvenient for GPU
+        # callers since the temporary is a host array.  Construct the
+        # tuple-backed static array directly instead.
+        alloc_vec = () -> zero(MVector{n, eltype(cond.z)})
     else
         alloc_vec = () -> zeros(n)
     end
@@ -202,10 +272,10 @@ function flash_storage_internal_newton!(out, eos, cond, method; static_size = fa
     V_ad = primary_ad(np)
     T = typeof(V_ad)
     if static_size
-        x_ad = @MVector zeros(T, n)
-        y_ad = @MVector zeros(T, n)
-        r = @MVector zeros(np)
-        J = @MMatrix zeros(np, np)
+        x_ad = zero(MVector{n, T})
+        y_ad = zero(MVector{n, T})
+        r = zero(MVector{np, Float64})
+        J = zero(MMatrix{np, np, Float64})
     else
         x_ad = zeros(T, n)
         y_ad = zeros(T, n)
@@ -235,8 +305,8 @@ function flash_storage_internal_inverse!(out, eos, cond, method; static_size = f
     T_ad = secondary_ad(2)
     T_cond = typeof(p_ad)
     if static_size
-        z_ad = @MVector zeros(T_cond, n)
-        J_inv = @MMatrix zeros(np, external_partials)
+        z_ad = zero(MVector{n, T_cond})
+        J_inv = zero(MMatrix{np, external_partials, Float64})
     else
         z_ad = zeros(T_cond, n)
         J_inv = zeros(np, external_partials)
@@ -248,7 +318,7 @@ function flash_storage_internal_inverse!(out, eos, cond, method; static_size = f
     cond_ad = (p = p_ad, T = T_ad, z = z_ad, phase = :liquid)
     if !isnothing(npartials)
         if static_size
-            buf = @MVector zeros(npartials)
+            buf = zero(MVector{npartials, Float64})
         else
             buf = zeros(npartials)
         end
@@ -279,12 +349,43 @@ function flash_update!(K, storage, type::SSIFlash, eos, cond, forces, V, iterati
 end
 
 function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces) where {F<:Real}
+    return ssi_with_phases!(K, p, T, x, y, z, V, eos, forces, :liquid, :vapor)
+end
+
+@inline function ssi_static!(K, p::F, T::F, x, y, z, V::F, eos, forces) where {F<:Real}
+    return ssi_with_phases!(K, p, T, x, y, z, V, eos, forces, Val(:liquid), Val(:vapor))
+end
+
+@inline function ssi_static(K::SVector{N, F}, p::F, T::F, z, V::F, eos, forces) where {N, F<:Real}
+    K_input = K
+    V_input = V
+    x = SVector{N, F}(ntuple(i -> liquid_mole_fraction(z[i], K_input[i], V_input), Val(N)))
+    y = SVector{N, F}(ntuple(i -> vapor_mole_fraction(x[i], K_input[i]), Val(N)))
+    liquid = (p = p, T = T, z = x, phase = Val(:liquid))
+    vapor = (p = p, T = T, z = y, phase = Val(:vapor))
+    Z_l, s_l = prep(eos, liquid, forces)
+    Z_v, s_v = prep(eos, vapor, forces)
+    ratios = SVector{N, F}(ntuple(Val(N)) do component
+        f_l = component_fugacity(eos, liquid, component, Z_l, forces, s_l)
+        f_v = component_fugacity(eos, vapor, component, Z_v, forces, s_v)
+        f_l/f_v
+    end)
+    residual = zero(F)
+    @inbounds for i in eachindex(ratios)
+        residual = max(residual, abs(one(F) - ratios[i]))
+    end
+    K_next = SVector{N, F}(ntuple(i -> K_input[i]*ratios[i], Val(N)))
+    V_next = solve_rachford_rice(K_next, z, V_input)
+    return V_next, K_next, residual
+end
+
+function ssi_with_phases!(K, p::F, T::F, x, y, z, V::F, eos, forces, liquid_phase, vapor_phase) where {F<:Real}
     # Initialize conditions for vapor and liquid phases based on K-values
     x = liquid_mole_fraction!(x, z, K, V)
     y = vapor_mole_fraction!(y, x, K)
 
-    liquid = (p = p, T = T, z = x, phase = :liquid)
-    vapor = (p = p, T = T, z = y, phase = :vapor)
+    liquid = (p = p, T = T, z = x, phase = liquid_phase)
+    vapor = (p = p, T = T, z = y, phase = vapor_phase)
 
     Z_l, s_l = prep(eos, liquid, forces)
     Z_v, s_v = prep(eos, vapor, forces)
diff --git a/src/kvalues.jl b/src/kvalues.jl
index 0c4109a..fba6245 100644
--- a/src/kvalues.jl
+++ b/src/kvalues.jl
@@ -60,4 +60,15 @@ end
 
 In-place version of `initial_guess_K`.
 """
-initial_guess_K!(K, eos, cond) = wilson_estimate!(K, eos, cond.p, cond.T)
\ No newline at end of file
+initial_guess_K!(K, eos, cond) = wilson_estimate!(K, eos, cond.p, cond.T)
+
+"""Return inline Wilson K-values for an EOS with a compile-time component count."""
+function initial_guess_K_static(eos::GenericCubicEOS{E, R, N}, cond) where {E, R, N}
+    T = Base.promote_eltype(cond.p, cond.T, cond.z[1])
+    return initial_guess_K_static(eos, cond, T)
+end
+
+@inline function initial_guess_K_static(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T}) where {E, R, N, T}
+    properties = eos.mixture.properties
+    return SVector{N, T}(ntuple(i -> wilson_estimate(properties[i], cond.p, cond.T), Val(N)))
+end
diff --git a/src/mixture_types.jl b/src/mixture_types.jl
index b019226..4b1d751 100644
--- a/src/mixture_types.jl
+++ b/src/mixture_types.jl
@@ -61,11 +61,11 @@ end
 
 Create a multicomponent mixture with an optional binary interaction coefficient matrix `A_ij`.
 """
-struct MultiComponentMixture{R, N}
-    name::String
-    component_names::Vector{String}
+struct MultiComponentMixture{R, N, Name, Names, BIC}
+    name::Name
+    component_names::Names
     properties::NTuple{N, MolecularProperty{R}}
-    binary_interaction::Union{Matrix{R}, Nothing}
+    binary_interaction::BIC
     function MultiComponentMixture(properties; A_ij = nothing, names = ["C$d" for d in 1:length(properties)], name = "UnnamedMixture")
         n = length(properties)
         n > 0 || throw(ArgumentError("At least one property must be present"))
@@ -77,10 +77,25 @@ struct MultiComponentMixture{R, N}
             A_ij = Symmetric(A_ij)
         end
         length(names) == n || throw(ArgumentError("Vector of component names must have same length as mixture."))
-        new{realtype, n}(name, names, properties, A_ij)
+        new{realtype, n, typeof(name), typeof(names), typeof(A_ij)}(name, names, properties, A_ij)
     end
 end
 
+"""
+    static_mixture(mixture)
+
+Return an isbits representation of a mixture for accelerator kernels. Human-readable
+metadata is omitted and binary-interaction coefficients are stored inline.
+"""
+function static_mixture(mixture::MultiComponentMixture{R, N}) where {R, N}
+    names = ntuple(_ -> nothing, Val(N))
+    bic = mixture.binary_interaction
+    if !isnothing(bic)
+        bic = SMatrix{N, N, R}(bic)
+    end
+    return MultiComponentMixture(mixture.properties; A_ij = bic, names = names, name = nothing)
+end
+
 """
     MultiComponentMixture(names::Vector{String}; A_ij = nothing, name = "UnnamedMixture")
 
diff --git a/src/rachford_rice.jl b/src/rachford_rice.jl
index 9120e70..1384a9e 100644
--- a/src/rachford_rice.jl
+++ b/src/rachford_rice.jl
@@ -20,6 +20,79 @@ julia> solve_rachford_rice([0.5, 1.5], [0.3, 0.7])
 0.8000000000000002
 ```
 """
+@inline function solve_rachford_rice(K::StaticVector{2}, z::StaticVector{2}, V = NaN)
+    z1, z2 = z
+    k1, k2 = K
+    b1, b2 = inv(1 - k1), inv(1 - k2)
+    return (z1*b2 + z2*b1)/(z1 + z2)
+end
+
+@inline function solve_rachford_rice(K::StaticVector{3}, z::StaticVector{3}, V = NaN)
+    z1, z2, z3 = z
+    k1, k2, k3 = K
+    b1, b2, b3 = inv(1-k1), inv(1-k2), inv(1-k3)
+    a2 = z1 + z2 + z3
+    a1 = -b1*(z2 + z3) - b2*(z1 + z3) - b3*(z1 + z2)
+    a0 = b1*b2*z3 + b1*b3*z2 + b2*b3*z1
+    discriminant = a1*a1 - 4*a0*a2
+    if discriminant >= zero(discriminant)
+        inv_2a2 = inv(2*a2)
+        root_offset = sqrt(discriminant)*inv_2a2
+        root_center = -a1*inv_2a2
+        root1 = root_center - root_offset
+        root2 = root_center + root_offset
+        if zero(root1) < root1 < one(root1)
+            return root1
+        elseif zero(root2) < root2 < one(root2)
+            return root2
+        elseif isfinite(root1 + root2)
+            kmin = min(k1, k2, k3)
+            kmax = max(k1, k2, k3)
+            kmin > one(kmin) && return max(root1, root2)
+            kmax < one(kmax) && return min(root1, root2)
+        end
+    end
+    return solve_rachford_rice_static_iterative(K, z, V)
+end
+
+@inline function solve_rachford_rice(K::StaticVector, z::StaticVector, V = NaN)
+    return solve_rachford_rice_static_iterative(K, z, V)
+end
+
+"""GPU-compatible Rachford--Rice fallback for statically sized inputs."""
+@inline function solve_rachford_rice_static_iterative(K, z, V; tol = 1e-12, maxiter = 1000)
+    V_lo = inv(1 - maximum(K))
+    V_hi = inv(1 - minimum(K))
+    if V_hi < V_lo
+        V_lo, V_hi = V_hi, V_lo
+    end
+    if isnan(V)
+        V = (V_lo + V_hi)/2
+    end
+    for _ in 1:maxiter
+        residual = zero(V)
+        denominator = zero(V)
+        @inbounds for i in eachindex(K)
+            delta_K = K[i] - one(K[i])
+            term_denominator = one(V) + V*delta_K
+            residual += z[i]*delta_K/term_denominator
+            denominator += z[i]*delta_K^2/term_denominator^2
+        end
+        abs(residual) < tol && break
+        if residual > zero(residual)
+            V_lo = V
+        else
+            V_hi = V
+        end
+        V_next = V + residual/denominator
+        if !(V_lo < V_next < V_hi) || !isfinite(V_next)
+            V_next = (V_lo + V_hi)/2
+        end
+        V = V_next
+    end
+    return V
+end
+
 function solve_rachford_rice(K, z, V = NaN; tol=1e-12, maxiter=1000, ad = false, analytical = true, verbose = false)
     n = length(z)
     if analytical
diff --git a/src/utils.jl b/src/utils.jl
index 63868eb..c387497 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -107,16 +107,23 @@ end
 function michelsen_critical_point_measure_storage(eos; T = Float64, static_size = true)
     n = number_of_components(eos)
     T∂ = ForwardDiff.Dual{nothing, T, n}
-    mole_numbers_ad = Vector{T∂}(undef, n)
-    for i in 1:n
-        partials = ForwardDiff.Partials{n, T}(Tuple(map(j -> Float64(i == j), 1:n)))
-        mole_numbers_ad[i] = T∂(1.0/n, partials)
+    if static_size
+        mole_numbers_ad = ntuple(n) do i
+            partials = ForwardDiff.Partials{n, T}(ntuple(j -> T(i == j), n))
+            T∂(one(T) / n, partials)
+        end
+    else
+        mole_numbers_ad = Vector{T∂}(undef, n)
+        for i in 1:n
+            partials = ForwardDiff.Partials{n, T}(Tuple(map(j -> T(i == j), 1:n)))
+            mole_numbers_ad[i] = T∂(one(T) / n, partials)
+        end
     end
     z = mole_numbers_ad./sum(mole_numbers_ad)
     if static_size
         z = MVector{n, T∂}(z)
     end
-    B = zeros(T, n, n)
+    B = static_size ? zero(MMatrix{n, n, T}) : zeros(T, n, n)
     c = (p = 101325.0, T = 303.15, z = z)
     forces = force_coefficients(eos, c; static_size = static_size)
     return (B = B, forces = forces, z = z)
@@ -148,7 +155,8 @@ function michelsen_critical_point_measure!(S, eos, p, T, mole_numbers)
         end
     end
     v = Inf
-    for x in eigvals!(B)
+    eigenvalues = B isa MMatrix ? eigvals(Symmetric(SMatrix(B))) : eigvals!(B)
+    for x in eigenvalues
         v = min(v, real(x))
     end
     return v
diff --git a/test/runtests.jl b/test/runtests.jl
index 1984694..1481a0b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -49,6 +49,27 @@ end
     end
 end
 
+@testset "Static accelerator path" begin
+    eos = static_eos(get_test_eos())
+    c = (p = 1e6, T = 300.0, z = @SVector [0.5, 0.3, 0.2])
+    K = zero(MVector{3, Float64})
+    initial_guess_K!(K, eos, c)
+    V, K, iterations, converged = flash_2ph_static(eos, c, K, 0.5)
+
+    @test isbitstype(typeof(eos))
+    @test converged
+    @test iterations == 6
+    @test V ≈ 0.7632068334421974
+    @test K ≈ @SVector [4.553402802323027, 17.73895830809456, 0.0004031451448211194]
+
+    K_static = initial_guess_K_static(eos, c)
+    V_static, K_static, iterations_static, converged_static = flash_2ph_static(eos, c, K_static, 0.5)
+    @test converged_static
+    @test iterations_static == iterations
+    @test V_static ≈ V
+    @test K_static ≈ K
+end
+
 @testset "Partial derivatives" begin
     test_flash_partials()
 end

From 45b853ea40758c1eefb9f16be877ae430ad97f4b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Fri, 11 Sep 2026 10:53:43 +0200
Subject: [PATCH 2/9] stability + unification on GPU

---
 src/MultiComponentFlash.jl |   2 +-
 src/eos.jl                 |  13 ----
 src/flash.jl               | 147 ++++++++++++++++++-------------------
 src/flash_types.jl         |  20 +++++
 src/kvalues.jl             |   2 +-
 src/kvalues_eos.jl         |   9 +++
 src/stability.jl           |  42 ++++++++---
 test/runtests.jl           |  10 +++
 8 files changed, 141 insertions(+), 104 deletions(-)

diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl
index 0d99ff1..84402fc 100644
--- a/src/MultiComponentFlash.jl
+++ b/src/MultiComponentFlash.jl
@@ -22,7 +22,7 @@ module MultiComponentFlash
     export flash_2ph, flash_2ph!, flash_2ph_static, flash_storage, flash_storage_static
     export stability_2ph, stability_2ph!
     # Algorithms for flash
-    export SSIFlash, NewtonFlash, SSINewtonFlash
+    export SSIFlash, NewtonFlash, SSINewtonFlash, FlashConfig
     # Mixtures and their molecular makeup
     export MolecularProperty, MultiComponentMixture
     # K-values
diff --git a/src/eos.jl b/src/eos.jl
index 303949b..15c11c7 100644
--- a/src/eos.jl
+++ b/src/eos.jl
@@ -185,19 +185,6 @@ function force_coefficients_static(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T
     return update_force_coefficients!(coeff, eos, cond)
 end
 
-"""Immutable force coefficients for register-oriented accelerator kernels."""
-@inline function force_coefficients_stack(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T}) where {E, R, N, T}
-    A_i_static = SVector{N, T}(ntuple(i -> A_i(eos, cond, i), Val(N)))
-    B_i_static = SVector{N, T}(ntuple(i -> B_i(eos, cond, i), Val(N)))
-    A_ij_static = SMatrix{N, N, T}(ntuple(Val(N*N)) do index
-        i = mod1(index, N)
-        j = (index - 1) ÷ N + 1
-        sqrt(A_i_static[i]*A_i_static[j]) *
-            (one(T) - binary_interaction(eos, i, j, cond))
-    end)
-    return (A_ij = A_ij_static, A_i = A_i_static, B_i = B_i_static)
-end
-
 function get_force_coefficients(forces, eos, cond)
     if forces_per_phase(eos)
         phase = get_phase(cond)
diff --git a/src/flash.jl b/src/flash.jl
index 2a86952..2dcaffa 100644
--- a/src/flash.jl
+++ b/src/flash.jl
@@ -33,16 +33,22 @@ 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), V = NaN; method = SSIFlash(), kwarg...) where T
+    return flash_2ph(eos, c, K, V, FlashConfig(); method = method, kwarg...)
+end
+
+function flash_2ph(eos, c::T, K, V, config::FlashConfig; method = SSIFlash(), kwarg...) where T
     nc = number_of_components(eos)
-    @assert hasfield(T, :p)
-    @assert hasfield(T, :T)
-    @assert hasfield(T, :z)
-    @assert length(c.z) == nc
-    @assert length(K) == nc
+    if print_output(config)
+        @assert hasfield(T, :p)
+        @assert hasfield(T, :T)
+        @assert hasfield(T, :z)
+        @assert length(c.z) == nc
+        @assert length(K) == nc
+    end
     method::AbstractFlash
 
-    S = flash_storage(eos, c, method = method)
-    flash_2ph!(S, K, eos, c, V, update_forces = false; method = method, kwarg...)
+    S = flash_storage(eos, c, method, config)
+    flash_2ph!(S, K, eos, c, V, config, update_forces = false; method = method, kwarg...)
 end
 
 """
@@ -69,7 +75,11 @@ function flash_2ph!(arg...; extra_out = false, kwarg...)
     end
 end
 
-function flash_2ph_impl!(storage, K, eos, c, V = NaN;
+function flash_2ph_impl!(storage, K, eos, c, V = NaN; kwarg...)
+    return flash_2ph_impl!(storage, K, eos, c, V, FlashConfig(); kwarg...)
+end
+
+function flash_2ph_impl!(storage, K, eos, c, V, config::FlashConfig;
         method = SSIFlash(),
         verbose::Bool = false,
         maxiter::Int = 25000,
@@ -92,7 +102,7 @@ function flash_2ph_impl!(storage, K, eos, c, V = NaN;
     end
     single_phase_init = isnan(V) || V == 1.0 || V == 0.0
     if single_phase_init
-        stable, stability_report = stability_2ph!(storage, K, eos, c;
+        stable, stability_report = stability_2ph!(storage, K, eos, c, config;
             maxiter = maxiter,
             verbose = verbose,
             extra_out = true,
@@ -101,7 +111,9 @@ function flash_2ph_impl!(storage, K, eos, c, V = NaN;
     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"
+        if print_output(config)
+            @assert all(isfinite, K) "K values must be finite: K = $K"
+        end
         stability_report = StabilityReport(stable_liquid = false, stable_vapor = false)
         stable = false
     end
@@ -114,13 +126,13 @@ function flash_2ph_impl!(storage, K, eos, c, V = NaN;
             V = solve_rachford_rice(K, z, V)
         end
         while true
-            V, ϵ = flash_update!(K, storage, method, eos, c, forces, V, i)
+            V, ϵ = flash_update!(K, storage, method, eos, c, forces, V, i, config)
             converged = ϵ ≤ tolerance
             if converged || i == maxiter
-                if verbose
+                if print_output(config) && verbose
                     @info "Flash done in $i iterations." V K converged
                 end
-                if check && !converged
+                if print_output(config) && check && !converged
                     @warn "Flash did not converge in $i iterations. Final ϵ = $ϵ > $tolerance = tolerance" c
                     if !isfinite(V)
                         # Hard error
@@ -152,7 +164,11 @@ Pre-allocate storage for `flash_2ph!`.
 
 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(), kwarg...)
+function flash_storage(eos, cond = (p = 10e5, T = 273.15, z = zeros(number_of_components(eos))); method = SSIFlash(), kwarg...)
+    return flash_storage(eos, cond, method, FlashConfig(); kwarg...)
+end
+
+function flash_storage(eos, cond, method, config::FlashConfig{PrintOutput, true}; kwarg...) where PrintOutput
     out = Dict{Symbol,Any}()
     d = flash_storage_internal!(out, eos, cond, method; kwarg...)
     # Convert to named tuple
@@ -166,9 +182,14 @@ Construct an inline, statically sized SSI workspace without using the `Dict`-bas
 host storage builder. This constructor is suitable for use inside accelerator
 kernels when `eos` and `cond` are isbits values.
 """
-function flash_storage_static(eos::GenericCubicEOS{E, R, N}, cond, method::SSIFlash = SSIFlash()) where {E, R, N}
-    T = Base.promote_eltype(cond.p, cond.T, cond.z[1])
-    forces = force_coefficients_static(eos, cond)
+function flash_storage(eos::GenericCubicEOS{E, R, N}, cond, method::SSIFlash, config::FlashConfig{PrintOutput, false}) where {E, R, N, PrintOutput}
+    T = eltype(cond.z)
+    return flash_storage(eos, cond, method, config, T)
+end
+
+@inline function flash_storage(eos::GenericCubicEOS{E, R, N}, cond, method::SSIFlash,
+        config::FlashConfig{PrintOutput, false}, ::Type{T}) where {E, R, N, PrintOutput, T}
+    forces = force_coefficients_static(eos, cond, T)
     return (
         forces = forces,
         x = zero(MVector{N, T}),
@@ -179,12 +200,15 @@ function flash_storage_static(eos::GenericCubicEOS{E, R, N}, cond, method::SSIFl
 end
 
 
+flash_storage_static(eos, cond, method::SSIFlash = SSIFlash()) =
+    flash_storage(eos, cond, method, FlashConfig{false, false}())
+
+
 """
     flash_2ph_static(eos, cond, K, V[, maxiter, tolerance])
 
-Accelerator-compatible SSI equilibrium solve with entirely local static storage.
-`K` must be a static vector and `V` must be a finite initial vapor
-fraction. Stability testing and host-side diagnostics are intentionally excluded.
+Compatibility wrapper around the config-driven flash implementation using entirely
+local static storage and no host-side diagnostics.
 Returns `(V, K, iterations, converged)`.
 """
 @inline flash_2ph_static(eos, cond, K, V) =
@@ -194,34 +218,12 @@ Returns `(V, K, iterations, converged)`.
         eos::GenericCubicEOS{E, R, N}, cond, K::MVector{N, T}, V,
         maxiter::Int, tolerance
     ) where {E, R, N, T}
-    forces = force_coefficients_static(eos, cond, T)
-    x = zero(MVector{N, T})
-    y = zero(MVector{N, T})
-    converged = false
-    iteration = 0
-    while iteration < maxiter
-        iteration += 1
-        V, residual = ssi_static!(K, cond.p, cond.T, x, y, cond.z, V, eos, forces)
-        converged = residual <= tolerance
-        converged && break
-    end
-    return (V, K, iteration, converged)
-end
-
-@inline function flash_2ph_static(
-        eos::GenericCubicEOS{E, R, N}, cond, K::SVector{N, T}, V,
-        maxiter::Int, tolerance
-    ) where {E, R, N, T}
-    forces = force_coefficients_stack(eos, cond, T)
-    converged = false
-    iteration = 0
-    while iteration < maxiter
-        iteration += 1
-        V, K, residual = ssi_static(K, cond.p, cond.T, cond.z, V, eos, forces)
-        converged = residual <= tolerance
-        converged && break
-    end
-    return (V, K, iteration, converged)
+    config = FlashConfig{false, false}()
+    storage = flash_storage(eos, cond, SSIFlash(), config, T)
+    V, K, report = flash_2ph_impl!(storage, K, eos, cond, V, config;
+        method = SSIFlash(), maxiter = maxiter, tolerance = tolerance,
+        verbose = false, check = false, z_min = nothing)
+    return (V, K, report.its, report.converged)
 end
 
 function flash_storage_internal!(out, eos, cond, method;
@@ -342,41 +344,24 @@ function update_value(v::T, newv::Real) where T<:ForwardDiff.Dual
 end
 
 function flash_update!(K, storage, type::SSIFlash, eos, cond, forces, V, iteration)
+    return flash_update!(K, storage, type, eos, cond, forces, V, iteration, FlashConfig())
+end
+
+function flash_update!(K, storage, type::SSIFlash, eos, cond, forces, V, iteration, config::FlashConfig)
     z = cond.z
     x, y = storage.x, storage.y
     p, T = cond.p, cond.T
-    return ssi!(K, p, T, x, y, z, V, eos, forces)
+    return ssi!(K, p, T, x, y, z, V, eos, forces, config)
 end
 
 function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces) where {F<:Real}
-    return ssi_with_phases!(K, p, T, x, y, z, V, eos, forces, :liquid, :vapor)
+    return ssi!(K, p, T, x, y, z, V, eos, forces, FlashConfig())
 end
 
-@inline function ssi_static!(K, p::F, T::F, x, y, z, V::F, eos, forces) where {F<:Real}
-    return ssi_with_phases!(K, p, T, x, y, z, V, eos, forces, Val(:liquid), Val(:vapor))
-end
-
-@inline function ssi_static(K::SVector{N, F}, p::F, T::F, z, V::F, eos, forces) where {N, F<:Real}
-    K_input = K
-    V_input = V
-    x = SVector{N, F}(ntuple(i -> liquid_mole_fraction(z[i], K_input[i], V_input), Val(N)))
-    y = SVector{N, F}(ntuple(i -> vapor_mole_fraction(x[i], K_input[i]), Val(N)))
-    liquid = (p = p, T = T, z = x, phase = Val(:liquid))
-    vapor = (p = p, T = T, z = y, phase = Val(:vapor))
-    Z_l, s_l = prep(eos, liquid, forces)
-    Z_v, s_v = prep(eos, vapor, forces)
-    ratios = SVector{N, F}(ntuple(Val(N)) do component
-        f_l = component_fugacity(eos, liquid, component, Z_l, forces, s_l)
-        f_v = component_fugacity(eos, vapor, component, Z_v, forces, s_v)
-        f_l/f_v
-    end)
-    residual = zero(F)
-    @inbounds for i in eachindex(ratios)
-        residual = max(residual, abs(one(F) - ratios[i]))
-    end
-    K_next = SVector{N, F}(ntuple(i -> K_input[i]*ratios[i], Val(N)))
-    V_next = solve_rachford_rice(K_next, z, V_input)
-    return V_next, K_next, residual
+@inline function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces, config::FlashConfig) where {F<:Real}
+    liquid_phase = phase_value(config, Val(:liquid))
+    vapor_phase = phase_value(config, Val(:vapor))
+    return ssi_with_phases!(K, p, T, x, y, z, V, eos, forces, liquid_phase, vapor_phase)
 end
 
 function ssi_with_phases!(K, p::F, T::F, x, y, z, V::F, eos, forces, liquid_phase, vapor_phase) where {F<:Real}
@@ -407,6 +392,10 @@ 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)
+    return flash_update!(K, storage, type, eos, cond, forces, V, iteration, FlashConfig())
+end
+
+function flash_update!(K, storage, type::NewtonFlash, eos, cond, forces, V, iteration, config::FlashConfig)
     x, y = storage.x, storage.y
     z = cond.z
     x = liquid_mole_fraction!(x, z, K, V)
@@ -499,10 +488,14 @@ function newton_dampen!(dMax, Δ)
 end
 
 function flash_update!(K, storage, type::SSINewtonFlash, eos, cond, forces, V, iteration)
+    return flash_update!(K, storage, type, eos, cond, forces, V, iteration, FlashConfig())
+end
+
+function flash_update!(K, storage, type::SSINewtonFlash, eos, cond, forces, V, iteration, config::FlashConfig)
     if iteration >= type.swap_iter
-        flash_update!(K, storage, NewtonFlash(type.dMax), eos, cond, forces, V, iteration)
+        flash_update!(K, storage, NewtonFlash(type.dMax), eos, cond, forces, V, iteration, config)
     else
-        flash_update!(K, storage, SSIFlash(), eos, cond, forces, V, iteration)
+        flash_update!(K, storage, SSIFlash(), eos, cond, forces, V, iteration, config)
     end
 end
 
diff --git a/src/flash_types.jl b/src/flash_types.jl
index eba83a5..11dc3f5 100644
--- a/src/flash_types.jl
+++ b/src/flash_types.jl
@@ -3,6 +3,26 @@ abstract type AbstractFlash end
 "Abstract type for all flash types that use Newton in some form"
 abstract type AbstractNewtonFlash <: AbstractFlash end
 
+"""
+    FlashConfig{PrintOutput, UseDictStorage}()
+    FlashConfig(; print_output=true, use_dict_storage=true)
+
+Compile-time configuration for the flash call graph. `PrintOutput=false` removes
+logging, warnings, and checked error paths from accelerator kernels.
+`UseDictStorage=false` selects direct static storage construction instead of the
+host-oriented `Dict` builder.
+"""
+struct FlashConfig{PrintOutput, UseDictStorage} end
+
+FlashConfig(; print_output::Bool = true, use_dict_storage::Bool = true) =
+    FlashConfig{print_output, use_dict_storage}()
+
+@inline print_output(::FlashConfig{PrintOutput}) where PrintOutput = PrintOutput
+@inline use_dict_storage(::FlashConfig{PrintOutput, UseDictStorage}) where {PrintOutput, UseDictStorage} = UseDictStorage
+
+@inline phase_value(::FlashConfig{PrintOutput, true}, ::Val{Phase}) where {PrintOutput, Phase} = Phase
+@inline phase_value(::FlashConfig{PrintOutput, false}, phase::Val) where PrintOutput = phase
+
 """
     ssi = SSIFlash()
 
diff --git a/src/kvalues.jl b/src/kvalues.jl
index fba6245..54e721d 100644
--- a/src/kvalues.jl
+++ b/src/kvalues.jl
@@ -70,5 +70,5 @@ end
 
 @inline function initial_guess_K_static(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T}) where {E, R, N, T}
     properties = eos.mixture.properties
-    return SVector{N, T}(ntuple(i -> wilson_estimate(properties[i], cond.p, cond.T), Val(N)))
+    return MVector{N, T}(ntuple(i -> wilson_estimate(properties[i], cond.p, cond.T), Val(N)))
 end
diff --git a/src/kvalues_eos.jl b/src/kvalues_eos.jl
index f19d3ec..33d3e2b 100644
--- a/src/kvalues_eos.jl
+++ b/src/kvalues_eos.jl
@@ -18,6 +18,15 @@ function flash_storage(eos::KValuesEOS, cond = missing; kwarg...)
     return nothing
 end
 
+flash_storage(eos::KValuesEOS, cond, method, config::FlashConfig{PrintOutput, true}) where PrintOutput = nothing
+flash_storage(eos::KValuesEOS, cond, method, config::FlashConfig{PrintOutput, false}) where PrintOutput = nothing
+
 function flash_2ph!(storage, K, eos::KValuesEOS, cond, V = NaN; kwarg...)
     return solve_rachford_rice(K, cond.z, V)
 end
+
+
+function flash_2ph!(storage, K, eos::KValuesEOS, cond, V,
+        config::FlashConfig{PrintOutput, UseDictStorage}; kwarg...) where {PrintOutput, UseDictStorage}
+    return solve_rachford_rice(K, cond.z, V)
+end
diff --git a/src/stability.jl b/src/stability.jl
index 92cb3f5..4d24bb0 100644
--- a/src/stability.jl
+++ b/src/stability.jl
@@ -8,11 +8,19 @@ This is done using a version of Michelsen's stability test.
 Reference: [The isothermal flash problem. Part I. Stability](https://doi.org/10.1016/0378-3812(82)85001-2)
 """
 function stability_2ph(eos, c, K = initial_guess_K(eos, c); kwarg...)
-    storage = flash_storage(eos, c)
-    stability_2ph!(storage, K, eos, c)
+    return stability_2ph(eos, c, K, FlashConfig(); kwarg...)
 end
 
-function stability_2ph!(storage, K, eos, c;
+function stability_2ph(eos, c, K, config::FlashConfig; kwarg...)
+    storage = flash_storage(eos, c, SSIFlash(), config)
+    stability_2ph!(storage, K, eos, c, config; kwarg...)
+end
+
+function stability_2ph!(storage, K, eos, c; kwarg...)
+    return stability_2ph!(storage, K, eos, c, FlashConfig(); kwarg...)
+end
+
+function stability_2ph!(storage, K, eos, c, config::FlashConfig;
         verbose::Bool = false,
         extra_out::Bool = false,
         check_vapor::Bool = true,
@@ -24,14 +32,16 @@ function stability_2ph!(storage, K, eos, c;
     f_xy = storage.buffer2
     x, y = storage.x, storage.y
     z, p, T = c.z, c.p, c.T
-    liquid = (p = p, T = T, z = x, phase = :liquid)
-    vapor = (p = p, T = T, z = y, phase = :vapor)
-    current_as_liquid = (p = p, T = T, z = z, phase = :liquid)
-    current_as_vapor = (p = p, T = T, z = z, phase = :vapor)
+    liquid_phase = phase_value(config, Val(:liquid))
+    vapor_phase = phase_value(config, Val(:vapor))
+    liquid = (p = p, T = T, z = x, phase = liquid_phase)
+    vapor = (p = p, T = T, z = y, phase = vapor_phase)
+    current_as_liquid = (p = p, T = T, z = z, phase = liquid_phase)
+    current_as_vapor = (p = p, T = T, z = z, phase = vapor_phase)
     mixture_fugacities!(f_z, eos, current_as_vapor, forces)
     if check_vapor
         wilson_estimate!(K, eos, p, T)
-        v = michelsen_test!(vapor, f_z, f_xy, vapor.z, z, K, eos, c, forces, Val(true); kwarg...)
+        v = michelsen_test!(vapor, f_z, f_xy, vapor.z, z, K, eos, c, forces, Val(true), config; kwarg...)
     else
         v = (true, true, 0)
     end
@@ -44,7 +54,7 @@ function stability_2ph!(storage, K, eos, c;
             mixture_fugacities!(f_z, eos, current_as_liquid, forces)
         end
         wilson_estimate!(K, eos, p, T)
-        l = michelsen_test!(liquid, f_z, f_xy, liquid.z, z, K, eos, c, forces, Val(false); kwarg...)
+        l = michelsen_test!(liquid, f_z, f_xy, liquid.z, z, K, eos, c, forces, Val(false), config; kwarg...)
     else
         l = (true, true, 0)
     end
@@ -59,7 +69,7 @@ function stability_2ph!(storage, K, eos, c;
     if !stable
         @. K = y/x
     end
-    if verbose
+    if print_output(config) && verbose
         @info "Stability done. Iterations:\nV: $i_v\nL: $i_l" stable_vapor stable_liquid stable
     end
     if extra_out
@@ -80,7 +90,13 @@ xy_value(z, K, ::Val{false}) = z/K
 
 In-place version of [`stability_2ph`](@ref). `storage` should be allocated by `flash_storage`.
 """
-function michelsen_test!(c_inside, f_z, f_xy, xy, z, K, eos, cond, forces, inside_is_vapor;
+function michelsen_test!(c_inside, f_z, f_xy, xy, z, K, eos, cond, forces, inside_is_vapor; kwarg...)
+    return michelsen_test!(c_inside, f_z, f_xy, xy, z, K, eos, cond, forces,
+        inside_is_vapor, FlashConfig(); kwarg...)
+end
+
+function michelsen_test!(c_inside, f_z, f_xy, xy, z, K, eos, cond, forces,
+        inside_is_vapor, config::FlashConfig;
         tol_equil = 1e-10,
         tol_trivial = tol_equil,
         tol_sat = tol_trivial,
@@ -122,7 +138,9 @@ function michelsen_test!(c_inside, f_z, f_xy, xy, z, K, eos, cond, forces, insid
         done = ok || iter == maxiter
         if done && !ok
             trivial = true
-            @warn "Stability test failed to converge in $maxiter iterations. Assuming stability." cond xy K_norm R_norm K
+            if print_output(config)
+                @warn "Stability test failed to converge in $maxiter iterations. Assuming stability." cond xy K_norm R_norm K
+            end
         end
     end
     stable = trivial || S <= 1.0 + tol_sat
diff --git a/test/runtests.jl b/test/runtests.jl
index 1481a0b..023e6eb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -68,6 +68,16 @@ end
     @test iterations_static == iterations
     @test V_static ≈ V
     @test K_static ≈ K
+
+    config = FlashConfig(print_output=false, use_dict_storage=false)
+    storage = flash_storage(eos, c, SSIFlash(), config)
+    K_config = initial_guess_K_static(eos, c)
+    V_config, K_config, report = flash_2ph!(storage, K_config, eos, c, NaN, config;
+        method=SSIFlash(), extra_out=true, z_min=nothing)
+    @test report.stability isa MultiComponentFlash.StabilityReport
+    @test report.converged
+    @test V_config ≈ V
+    @test K_config ≈ K
 end
 
 @testset "Partial derivatives" begin

From d6f6e66d6ba25bd344edbd53022a9304c66318f0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Fri, 11 Sep 2026 11:49:12 +0200
Subject: [PATCH 3/9] Update flash.jl

---
 src/flash.jl | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/flash.jl b/src/flash.jl
index 2dcaffa..103168f 100644
--- a/src/flash.jl
+++ b/src/flash.jl
@@ -221,8 +221,13 @@ Returns `(V, K, iterations, converged)`.
     config = FlashConfig{false, false}()
     storage = flash_storage(eos, cond, SSIFlash(), config, T)
     V, K, report = flash_2ph_impl!(storage, K, eos, cond, V, config;
-        method = SSIFlash(), maxiter = maxiter, tolerance = tolerance,
-        verbose = false, check = false, z_min = nothing)
+        method = SSIFlash(),
+        maxiter = maxiter,
+        tolerance = tolerance,
+        verbose = false,
+        check = false,
+        z_min = nothing
+    )
     return (V, K, report.its, report.converged)
 end
 

From fa654bb0a8fc4def88d8f33f1ba20d95cbfe7723 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Fri, 11 Sep 2026 13:52:54 +0200
Subject: [PATCH 4/9] add test

---
 Project.toml     |  6 +++++-
 test/runtests.jl | 40 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 45 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index a3d8863..7623787 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,6 +12,8 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 ForwardDiff = "0.10, 1"
+JLArrays = "0.1, 0.2, 0.3"
+KernelAbstractions = "0.9"
 LinearAlgebra = "1"
 Printf = "1"
 Roots = "1, 2, 3"
@@ -20,8 +22,10 @@ julia = "1.6"
 
 [extras]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "ForwardDiff", "StaticArrays"]
+test = ["Test", "ForwardDiff", "JLArrays", "KernelAbstractions", "StaticArrays"]
diff --git a/test/runtests.jl b/test/runtests.jl
index 023e6eb..134ab02 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -146,3 +146,43 @@ end
     @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
 end
+
+using StaticArrays, KernelAbstractions, JLArrays
+@testset "Static flash with KernelAbstractions/JLArrays" begin
+    @kernel function static_flash_kernel!(out, pressure, temperature, z, eos)
+        i = @index(Global)
+        if i <= length(out)
+            @inbounds cond = (p = pressure[i], T = temperature[i], z = z)
+            config = FlashConfig{false, false}()
+            K = initial_guess_K_static(eos, cond)
+            V = flash_2ph(eos, cond, K, NaN, config;
+                method = SSIFlash(), check = false, verbose = false, z_min = nothing)
+            @inbounds out[i] = V
+        end
+    end
+    if isdefined(JLArrays, :JLBackend)
+        host_eos = get_test_eos()
+        eos = static_eos(host_eos)
+        z = @SVector [0.5, 0.3, 0.2]
+        n = 16
+        pressure_host = collect(range(1e5, 4e6, length = n))
+        temperature_host = collect(range(280.0, 320.0, length = n))
+        expected = map(pressure_host, temperature_host) do p, T
+            flash_2ph(host_eos, (p = p, T = T, z = collect(z));
+                method = SSIFlash(), check = false)
+        end
+
+        pressure = JLArray(pressure_host)
+        temperature = JLArray(temperature_host)
+        out = JLArray(zeros(n))
+        backend = JLArrays.JLBackend()
+        kernel! = static_flash_kernel!(backend, 8)
+        kernel!(out, pressure, temperature, z, eos; ndrange = n)
+
+        @test Array(out) ≈ expected rtol = 1e-11
+    else
+        # JLArrays 0.1 supports Julia 1.6 but predates the KernelAbstractions backend.
+        @test_skip false
+    end
+end
+

From 1f7ea61069e122533b6a8243bbbe22d55efb1bf8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Fri, 11 Sep 2026 19:08:27 +0200
Subject: [PATCH 5/9] GPU allocation fix

---
 src/eos.jl         |  13 ++++++
 src/flash.jl       |  73 ++++++++++++++++++++++++++++++++
 src/flash_types.jl |  17 +++++---
 src/kvalues.jl     |   7 +++
 src/stability.jl   | 103 +++++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl   |   2 +-
 6 files changed, 209 insertions(+), 6 deletions(-)

diff --git a/src/eos.jl b/src/eos.jl
index 15c11c7..e974b54 100644
--- a/src/eos.jl
+++ b/src/eos.jl
@@ -185,6 +185,19 @@ function force_coefficients_static(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T
     return update_force_coefficients!(coeff, eos, cond)
 end
 
+"""Immutable force coefficients for heap-free accelerator kernels."""
+@inline function force_coefficients_stack(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T}) where {E, R, N, T}
+    A_i_static = SVector{N, T}(ntuple(i -> A_i(eos, cond, i), Val(N)))
+    B_i_static = SVector{N, T}(ntuple(i -> B_i(eos, cond, i), Val(N)))
+    A_ij_static = SMatrix{N, N, T}(ntuple(Val(N*N)) do index
+        i = mod1(index, N)
+        j = (index - 1) ÷ N + 1
+        sqrt(A_i_static[i]*A_i_static[j]) *
+            (one(T) - binary_interaction(eos, i, j, cond))
+    end)
+    return (A_ij = A_ij_static, A_i = A_i_static, B_i = B_i_static)
+end
+
 function get_force_coefficients(forces, eos, cond)
     if forces_per_phase(eos)
         phase = get_phase(cond)
diff --git a/src/flash.jl b/src/flash.jl
index 103168f..dda135e 100644
--- a/src/flash.jl
+++ b/src/flash.jl
@@ -36,6 +36,53 @@ function flash_2ph(eos, c::T, K = initial_guess_K(eos, c), V = NaN; method = SSI
     return flash_2ph(eos, c, K, V, FlashConfig(); method = method, kwarg...)
 end
 
+@inline function flash_2ph(eos::GenericCubicEOS{E, R, N}, c,
+        K::SVector{N, F}, V, config::FlashConfig{false, false};
+        method::SSIFlash = SSIFlash(),
+        verbose::Bool = false,
+        maxiter::Int = 25000,
+        tolerance::Float64 = 1e-8,
+        extra_out::Bool = false,
+        update_forces::Bool = true,
+        check::Bool = true,
+        z_min = MINIMUM_COMPOSITION,
+        kwarg...
+    ) where {E, R, N, F}
+    z = if isnothing(z_min)
+        c.z
+    else
+        SVector{N, F}(ntuple(i -> max(c.z[i], z_min), Val(N)))
+    end
+    cond = (p = c.p, T = c.T, z = z)
+    forces = force_coefficients_stack(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 = stability_2ph_stack(
+            K, eos, cond, forces; maxiter = maxiter, kwarg...)
+    else
+        stable = false
+        stability_report = StabilityReport(false, false, false, false)
+    end
+    converged = false
+    if stable
+        iteration = 0
+    else
+        iteration = 1
+        if isnan(V)
+            V = solve_rachford_rice(K, z, V)
+        end
+        while true
+            V, K, residual = ssi_stack(K, cond.p, cond.T, z, V, eos, forces)
+            converged = residual <= tolerance
+            (converged || iteration == maxiter) && break
+            iteration += 1
+        end
+    end
+    out = (V, K, (its = iteration, converged = converged, stability = stability_report))
+    return extra_out ? out : out[1]
+end
+
 function flash_2ph(eos, c::T, K, V, config::FlashConfig; method = SSIFlash(), kwarg...) where T
     nc = number_of_components(eos)
     if print_output(config)
@@ -369,6 +416,32 @@ end
     return ssi_with_phases!(K, p, T, x, y, z, V, eos, forces, liquid_phase, vapor_phase)
 end
 
+@inline 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)
+end
+
+@inline function ssi_stack(K::SVector{N, F}, p::F, T::F, z, V::F,
+        eos, forces) where {N, F<:Real}
+    x = SVector{N, F}(ntuple(i -> liquid_mole_fraction(z[i], K[i], V), Val(N)))
+    y = SVector{N, F}(ntuple(i -> vapor_mole_fraction(x[i], K[i]), Val(N)))
+    liquid = (p = p, T = T, z = x, phase = Val(:liquid))
+    vapor = (p = p, T = T, z = y, phase = Val(:vapor))
+    f_l = static_fugacities(eos, liquid, forces, F)
+    f_v = static_fugacities(eos, vapor, forces, F)
+    ratios = SVector{N, F}(ntuple(i -> f_l[i]/f_v[i], Val(N)))
+    residual = zero(F)
+    @inbounds for i in 1:N
+        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)
+    return V_next, K_next, residual
+end
+
 function ssi_with_phases!(K, p::F, T::F, x, y, z, V::F, eos, forces, liquid_phase, vapor_phase) where {F<:Real}
     # Initialize conditions for vapor and liquid phases based on K-values
     x = liquid_mole_fraction!(x, z, K, V)
diff --git a/src/flash_types.jl b/src/flash_types.jl
index 11dc3f5..036346e 100644
--- a/src/flash_types.jl
+++ b/src/flash_types.jl
@@ -66,11 +66,14 @@ end
 struct PhaseStabilityStatus
     stable::Bool
     trivial::Bool
-    function PhaseStabilityStatus(stable = false; trivial = stable)
+    function PhaseStabilityStatus(stable::Bool, trivial::Bool)
         return new(stable, trivial && stable)
     end
 end
 
+PhaseStabilityStatus(stable::Bool = false; trivial::Bool = stable) =
+    PhaseStabilityStatus(stable, trivial)
+
 function Base.show(io::IOContext, sr::PhaseStabilityStatus)
     compact = get(io, :compact, false)
     s = sr.stable ? "single-phase" : "two-phase"
@@ -86,16 +89,20 @@ struct StabilityReport
     stable::Bool
     liquid::PhaseStabilityStatus
     vapor::PhaseStabilityStatus
+    function StabilityReport(stable_liquid::Bool, trivial_liquid::Bool,
+            stable_vapor::Bool, trivial_vapor::Bool)
+        new(stable_liquid && stable_vapor,
+            PhaseStabilityStatus(stable_liquid, trivial_liquid),
+            PhaseStabilityStatus(stable_vapor, trivial_vapor)
+        )
+    end
     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)
-        )
+        StabilityReport(stable_liquid, trivial_liquid, stable_vapor, trivial_vapor)
     end
 end
 
diff --git a/src/kvalues.jl b/src/kvalues.jl
index 54e721d..b42f81e 100644
--- a/src/kvalues.jl
+++ b/src/kvalues.jl
@@ -72,3 +72,10 @@ end
     properties = eos.mixture.properties
     return MVector{N, T}(ntuple(i -> wilson_estimate(properties[i], cond.p, cond.T), Val(N)))
 end
+
+@inline function initial_guess_K_static(eos::GenericCubicEOS{E, R, N}, cond,
+        ::FlashConfig{false, false}) where {E, R, N}
+    T = Base.promote_eltype(cond.p, cond.T, cond.z[1])
+    properties = eos.mixture.properties
+    return SVector{N, T}(ntuple(i -> wilson_estimate(properties[i], cond.p, cond.T), Val(N)))
+end
diff --git a/src/stability.jl b/src/stability.jl
index 4d24bb0..5e7479e 100644
--- a/src/stability.jl
+++ b/src/stability.jl
@@ -85,6 +85,109 @@ f_ratio(f_z, f_xy, ::Val{false}) = f_xy/f_z
 xy_value(z, K, ::Val{true}) = z*K
 xy_value(z, K, ::Val{false}) = z/K
 
+stability_phase(::Val{true}) = Val(:vapor)
+stability_phase(::Val{false}) = Val(:liquid)
+
+@generated function stability_xy(z::SVector{N, F}, K::SVector{N, F}, phase) where {N, F}
+    values = [:(xy_value(z[$i], K[$i], phase)) for i in 1:N]
+    return :(SVector{N, F}(($(values...),)))
+end
+
+@generated function static_scale(v::SVector{N, F}, scale::F) where {N, F}
+    values = [:(v[$i]/scale) for i in 1:N]
+    return :(SVector{N, F}(($(values...),)))
+end
+
+@generated function stability_ratios(f_z::SVector{N, F}, f_xy::SVector{N, F},
+        scale::F, phase) where {N, F}
+    values = [:(f_ratio(f_z[$i], scale*f_xy[$i], phase)) for i in 1:N]
+    return :(SVector{N, F}(($(values...),)))
+end
+
+@generated function static_multiply(a::SVector{N, F}, b::SVector{N, F}) where {N, F}
+    values = [:(a[$i]*b[$i]) for i in 1:N]
+    return :(SVector{N, F}(($(values...),)))
+end
+
+@generated function static_divide(a::SVector{N, F}, b::SVector{N, F}) where {N, F}
+    values = [:(a[$i]/b[$i]) for i in 1:N]
+    return :(SVector{N, F}(($(values...),)))
+end
+
+@inline function michelsen_test_stack(f_z, z::SVector{N, F}, K::SVector{N, F},
+        eos, cond, forces, inside_is_vapor;
+        tol_equil = 1e-10,
+        tol_trivial = tol_equil,
+        tol_sat = tol_trivial,
+        maxiter = 1000
+    ) where {N, F}
+    trivial = false
+    S = one(F)
+    iter = 0
+    xy = zero(SVector{N, F})
+    while true
+        iter += 1
+        unnormalized = stability_xy(z, K, inside_is_vapor)
+        S = zero(F)
+        @inbounds for i in 1:N
+            S += unnormalized[i]
+        end
+        xy = static_scale(unnormalized, S)
+        inside = (p = cond.p, T = cond.T, z = xy,
+            phase = stability_phase(inside_is_vapor))
+        f_xy = static_fugacities(eos, inside, forces, F)
+        ratios = stability_ratios(f_z, f_xy, S, inside_is_vapor)
+        K = static_multiply(K, ratios)
+        R_norm = zero(F)
+        K_norm = zero(F)
+        @inbounds for i in 1:N
+            R_norm += (ratios[i] - one(F))^2
+            K_norm += log(K[i])^2
+        end
+        trivial = K_norm < tol_trivial
+        converged = R_norm < tol_equil
+        if trivial || converged
+            break
+        elseif iter == maxiter
+            # Match the checked host path: failed stability iterations are
+            # conservatively treated as a trivial, stable solution.
+            trivial = true
+            break
+        end
+    end
+    stable = trivial || S <= one(F) + tol_sat
+    return stable, trivial, iter, K, xy
+end
+
+@inline function stability_2ph_stack(K::SVector{N, F}, eos, cond, forces;
+        check_vapor::Bool = true,
+        check_liquid::Bool = true,
+        kwarg...
+    ) where {N, F}
+    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_static(eos, cond, FlashConfig{false, false}())
+    if check_vapor
+        stable_vapor, trivial_vapor, i_v, K_vapor, y = michelsen_test_stack(
+            f_z_vapor, cond.z, K_wilson, eos, cond, forces, Val(true); kwarg...)
+    else
+        stable_vapor, trivial_vapor, i_v, K_vapor, y = true, true, 0, K_wilson, cond.z
+    end
+    if check_liquid
+        liquid_phase = (p = cond.p, T = cond.T, z = cond.z, phase = Val(:liquid))
+        f_z_liquid = forces_per_phase(eos) ?
+            static_fugacities(eos, liquid_phase, forces, F) : f_z_vapor
+        stable_liquid, trivial_liquid, i_l, K_liquid, x = michelsen_test_stack(
+            f_z_liquid, cond.z, K_wilson, eos, cond, forces, Val(false); kwarg...)
+    else
+        stable_liquid, trivial_liquid, i_l, K_liquid, x = true, true, 0, K_wilson, cond.z
+    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
+end
+
 """
     stability_2ph!(storage, eos, c, [K])
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 134ab02..eda6d75 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -154,7 +154,7 @@ using StaticArrays, KernelAbstractions, JLArrays
         if i <= length(out)
             @inbounds cond = (p = pressure[i], T = temperature[i], z = z)
             config = FlashConfig{false, false}()
-            K = initial_guess_K_static(eos, cond)
+            K = initial_guess_K_static(eos, cond, config)
             V = flash_2ph(eos, cond, K, NaN, config;
                 method = SSIFlash(), check = false, verbose = false, z_min = nothing)
             @inbounds out[i] = V

From e39bf4f0d2691bd8cc9a5662f09019e49ae58dbf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Fri, 11 Sep 2026 20:11:53 +0200
Subject: [PATCH 6/9] Clean up and consolidate the static variants into one

---
 docs/src/api/flash.md         |   2 +-
 docs/src/examples/advanced.md |  65 ++-----
 src/MultiComponentFlash.jl    |   7 +-
 src/eos.jl                    |  41 ----
 src/eos_types.jl              |  24 ---
 src/flash.jl                  | 238 ++++-------------------
 src/flash_types.jl            |  23 +--
 src/kvalues.jl                |  18 --
 src/kvalues_eos.jl            |   5 +-
 src/mixture_types.jl          |  15 --
 src/rachford_rice.jl          |  73 -------
 src/stability.jl              | 116 +----------
 src/static.jl                 | 357 ++++++++++++++++++++++++++++++++++
 test/runtests.jl              |  54 +++--
 test/test_setup.jl            |   8 +-
 15 files changed, 456 insertions(+), 590 deletions(-)
 create mode 100644 src/static.jl

diff --git a/docs/src/api/flash.md b/docs/src/api/flash.md
index 0c85aa9..673af48 100644
--- a/docs/src/api/flash.md
+++ b/docs/src/api/flash.md
@@ -13,7 +13,7 @@ Private = false
 ## Two-phase flash
 ```@autodocs
 Modules = [MultiComponentFlash]
-Pages   = ["flash.jl", "flash_types.jl"]
+Pages   = ["flash.jl", "flash_types.jl", "static.jl"]
 Order   = [:type, :function]
 Private = false
 ```
diff --git a/docs/src/examples/advanced.md b/docs/src/examples/advanced.md
index 7b225cc..94149b4 100644
--- a/docs/src/examples/advanced.md
+++ b/docs/src/examples/advanced.md
@@ -38,64 +38,27 @@ S = flash_storage(eos, conditions, method = m)
 16
 ```
 
-See the unit tests for examples where the flash can use `StaticArrays` to avoid allocations entirely.
+## GPU and static-array use
 
-## Performance example
-
-The default interface is designed for ease-of-use with standard Julia types, but the module  also supports further by using `StaticArrays`:
+Use the fully static SSI path for GPU kernels. Construct the isbits EOS and storage
+once outside the kernel, then pass both to it:
 
 ```julia
-using MultiComponentFlash, BenchmarkTools, StaticArrays
-function bench(m, static_size = false)
-    p = 6e6
-    T = 480.0
-    # Take the SPE5 benchmark
-    eos, data = cubic_benchmark("spe5")
-    n = number_of_components(eos)
-    z = repeat([1/n], n)
-    conditions = (p = p, T = T, z = z)
-    S = flash_storage(eos, conditions, method = m, static_size = static_size)
-    K = initial_guess_K(eos, conditions)
-    if static_size
-        N = number_of_components(eos)
-        K = MVector{N}(K)
-    end
-    V, K, status = flash_2ph!(S, K, eos, conditions, NaN, method = m, extra_out = true)
-    println("V = $V (Completed in $(status.its) iterations)")
-    @btime flash_2ph!($S, $K, $eos, $conditions, NaN, method = $m)
-    return nothing
-end
-println("SSI:")
-bench(SSIFlash())
-println("SSI (static arrays):")
-bench(SSIFlash(), true)
-##
-println("Newton:")
-bench(NewtonFlash())
-println("Newton (static arrays):")
-bench(NewtonFlash(), true)
-```
+using MultiComponentFlash, StaticArrays
 
-The output will be a bit different on other CPUs, but this flash generally takes around 20 microseconds to complete, including both stability test and flash.
+eos_gpu = static_eos(eos)
+conditions_gpu = (p = p, T = T, z = SVector{length(z)}(z))
+storage = flash_storage(eos_gpu, conditions_gpu; method = SSIFlash(), static = true)
 
-```julia
-SSI:
-V = 0.03279769425318795 (Completed in 14 iterations)
-  18.500 μs (0 allocations: 0 bytes)
-SSI (static arrays):
-V = 0.03279769425318795 (Completed in 14 iterations)
-  16.500 μs (0 allocations: 0 bytes)
-
-Newton:
-V = 0.032797694260046494 (Completed in 4 iterations)
-  20.100 μs (0 allocations: 0 bytes)
-Newton (static arrays):
-V = 0.032797694260046494 (Completed in 4 iterations)
-  19.900 μs (0 allocations: 0 bytes)
+# Inside each kernel work item:
+K = initial_guess_K(eos_gpu, conditions_gpu, storage)
+V = flash_2ph!(storage, K, eos_gpu, conditions_gpu, NaN)
 ```
 
-!!! note "Use of `StaticArrays`"
-    Switching to statically sized arrays can improve the speed, at the cost of longer compilation times. Please note that for `StaticArrays` there will be compilation that is dependent on the number of components in your mixture. For example, switching from a five to six component mixture will trigger a full recompilation of your chosen flash.
+`static=false` returns the ordinary reusable mutable storage. `static=true` returns
+a `StaticConfig` and keeps all per-flash arrays immutable and local to the work item.
+The static path currently supports `GenericCubicEOS` with `SSIFlash`; compilation is
+specialized on the component count.
 
 ## Generate and plot a phase diagram
 
diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl
index 84402fc..218cb14 100644
--- a/src/MultiComponentFlash.jl
+++ b/src/MultiComponentFlash.jl
@@ -19,14 +19,14 @@ module MultiComponentFlash
     export KValuesEOS
     export number_of_components
     # Flash interfaces
-    export flash_2ph, flash_2ph!, flash_2ph_static, flash_storage, flash_storage_static
+    export flash_2ph, flash_2ph!, flash_storage
     export stability_2ph, stability_2ph!
     # Algorithms for flash
-    export SSIFlash, NewtonFlash, SSINewtonFlash, FlashConfig
+    export SSIFlash, NewtonFlash, SSINewtonFlash, FlashConfig, StaticConfig
     # Mixtures and their molecular makeup
     export MolecularProperty, MultiComponentMixture
     # K-values
-    export wilson_estimate, wilson_estimate!, initial_guess_K, initial_guess_K!, initial_guess_K_static
+    export wilson_estimate, wilson_estimate!, initial_guess_K, initial_guess_K!
     # Vapor-liquid equilibrium
     export solve_rachford_rice
 
@@ -54,6 +54,7 @@ module MultiComponentFlash
     include("flash.jl")
     include("derivatives.jl")
     include("stability.jl")
+    include("static.jl")
     include("tables.jl")
 
     include("flow_coupler.jl")
diff --git a/src/eos.jl b/src/eos.jl
index e974b54..d8264ef 100644
--- a/src/eos.jl
+++ b/src/eos.jl
@@ -15,7 +15,6 @@ function get_phase(cond)
 end
 
 @inline phase_symbol(phase::Symbol) = phase
-@inline phase_symbol(::Val{phase}) where phase = phase
 
 function set_phase(cond, phase::Symbol, throw::Bool = false)
     if throw && haskey(cond, :phase) && cond.phase != :unknown
@@ -113,16 +112,6 @@ end
     return pick_root(eos, roots, cond, forces, scalars, phase)
 end
 
-@inline function pick_root(eos, roots, cond, forces, scalars, ::Val{:liquid})
-    min_r, _ = root_bounds(roots, minimum_allowable_root(eos, forces, scalars))
-    return min_r
-end
-
-@inline function pick_root(eos, roots, cond, forces, scalars, ::Val{:vapor})
-    _, max_r = root_bounds(roots, minimum_allowable_root(eos, forces, scalars))
-    return max_r
-end
-
 function pick_root(eos, roots, cond, forces, scalars, phase::Symbol)
     min_r, max_r = root_bounds(roots, minimum_allowable_root(eos, forces, scalars))
     if min_r == max_r || phase == :liquid
@@ -170,34 +159,6 @@ function force_coefficients(eos::AbstractCubicEOS, cond; static_size = false)
     return coeff
 end
 
-"""Construct inline force coefficients with the component count from the EOS type."""
-function force_coefficients_static(eos::GenericCubicEOS{E, R, N}, cond) where {E, R, N}
-    T = Base.promote_eltype(cond.p, cond.T, cond.z[1])
-    return force_coefficients_static(eos, cond, T)
-end
-
-function force_coefficients_static(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T}) where {E, R, N, T}
-    coeff = (
-        A_ij = zero(MMatrix{N, N, T}),
-        A_i = zero(MVector{N, T}),
-        B_i = zero(MVector{N, T})
-    )
-    return update_force_coefficients!(coeff, eos, cond)
-end
-
-"""Immutable force coefficients for heap-free accelerator kernels."""
-@inline function force_coefficients_stack(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T}) where {E, R, N, T}
-    A_i_static = SVector{N, T}(ntuple(i -> A_i(eos, cond, i), Val(N)))
-    B_i_static = SVector{N, T}(ntuple(i -> B_i(eos, cond, i), Val(N)))
-    A_ij_static = SMatrix{N, N, T}(ntuple(Val(N*N)) do index
-        i = mod1(index, N)
-        j = (index - 1) ÷ N + 1
-        sqrt(A_i_static[i]*A_i_static[j]) *
-            (one(T) - binary_interaction(eos, i, j, cond))
-    end)
-    return (A_ij = A_ij_static, A_i = A_i_static, B_i = B_i_static)
-end
-
 function get_force_coefficients(forces, eos, cond)
     if forces_per_phase(eos)
         phase = get_phase(cond)
@@ -213,8 +174,6 @@ function get_force_coefficients(forces, eos, cond)
     end
 end
 
-@inline get_force_coefficients(forces, eos::GenericCubicEOS, cond) = forces
-
 """
     force_coefficients!(coeff, eos, cond)
 
diff --git a/src/eos_types.jl b/src/eos_types.jl
index 64517d3..55f57ef 100644
--- a/src/eos_types.jl
+++ b/src/eos_types.jl
@@ -182,30 +182,6 @@ 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
 
-"""
-    static_eos(eos)
-
-Convert a generic cubic EOS to an isbits representation suitable for passing to
-accelerator kernels. The numerical model is unchanged; only descriptive metadata
-and array storage are converted.
-"""
-function static_eos(eos::GenericCubicEOS{T, R, N}) where {T, R, N}
-    mixture = static_mixture(eos.mixture)
-    volume_shift = eos.volume_shift
-    if !isnothing(volume_shift)
-        volume_shift = SVector{N, eltype(volume_shift)}(volume_shift)
-    end
-    return GenericCubicEOS(
-        eos.type,
-        mixture,
-        eos.m_1,
-        eos.m_2,
-        eos.ω_a,
-        eos.ω_b,
-        volume_shift
-    )
-end
-
 struct KValuesEOS{T, R, N, V} <: AbstractEOS
     "Callable on the form `cond -> V` or a set of constants (Tuple/AbstractVector)"
     K_values_evaluator::T
diff --git a/src/flash.jl b/src/flash.jl
index dda135e..412adce 100644
--- a/src/flash.jl
+++ b/src/flash.jl
@@ -36,54 +36,10 @@ function flash_2ph(eos, c::T, K = initial_guess_K(eos, c), V = NaN; method = SSI
     return flash_2ph(eos, c, K, V, FlashConfig(); method = method, kwarg...)
 end
 
-@inline function flash_2ph(eos::GenericCubicEOS{E, R, N}, c,
-        K::SVector{N, F}, V, config::FlashConfig{false, false};
-        method::SSIFlash = SSIFlash(),
-        verbose::Bool = false,
-        maxiter::Int = 25000,
-        tolerance::Float64 = 1e-8,
-        extra_out::Bool = false,
-        update_forces::Bool = true,
-        check::Bool = true,
-        z_min = MINIMUM_COMPOSITION,
-        kwarg...
-    ) where {E, R, N, F}
-    z = if isnothing(z_min)
-        c.z
-    else
-        SVector{N, F}(ntuple(i -> max(c.z[i], z_min), Val(N)))
-    end
-    cond = (p = c.p, T = c.T, z = z)
-    forces = force_coefficients_stack(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 = stability_2ph_stack(
-            K, eos, cond, forces; maxiter = maxiter, kwarg...)
-    else
-        stable = false
-        stability_report = StabilityReport(false, false, false, false)
-    end
-    converged = false
-    if stable
-        iteration = 0
-    else
-        iteration = 1
-        if isnan(V)
-            V = solve_rachford_rice(K, z, V)
-        end
-        while true
-            V, K, residual = ssi_stack(K, cond.p, cond.T, z, V, eos, forces)
-            converged = residual <= tolerance
-            (converged || iteration == maxiter) && break
-            iteration += 1
-        end
-    end
-    out = (V, K, (its = iteration, converged = converged, stability = stability_report))
-    return extra_out ? out : out[1]
-end
-
 function flash_2ph(eos, c::T, K, V, config::FlashConfig; method = SSIFlash(), kwarg...) where T
+    if !use_dict_storage(config)
+        return flash_2ph(eos, c, K, V, StaticConfig(); method = method, kwarg...)
+    end
     nc = number_of_components(eos)
     if print_output(config)
         @assert hasfield(T, :p)
@@ -173,7 +129,7 @@ function flash_2ph_impl!(storage, K, eos, c, V, config::FlashConfig;
             V = solve_rachford_rice(K, z, V)
         end
         while true
-            V, ϵ = flash_update!(K, storage, method, eos, c, forces, V, i, config)
+            V, ϵ = flash_update!(K, storage, method, eos, c, forces, V, i)
             converged = ϵ ≤ tolerance
             if converged || i == maxiter
                 if print_output(config) && verbose
@@ -205,87 +161,38 @@ Pre-allocate storage for `flash_2ph!`.
 
 # Keyword arguments
 - `method = SSIFlash()`: Flash method to use. Can be `SSIFlash()`, `NewtonFlash()` or `SSINewtonFlash()`.
-- `static_size = false`: Use `SArrays` and `MArrays` for fast flash, but slower compile times.
+- `static = false`: Return fully static, GPU-compatible storage when `true`. The
+  static path currently supports `SSIFlash`.
 - `inc_jac`: Allocate storage for Newton/Jacobian. Required for Newton (and defaults to `true` for that method) or for `diff_externals`.
 - `diff_externals = false`: Allocate storage for matrix inversion required to produce partial derivatives of flash using `set_partials`.
 
 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(), kwarg...)
-    return flash_storage(eos, cond, method, FlashConfig(); kwarg...)
+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."))
+    config = static ? StaticConfig() : FlashConfig()
+    return flash_storage(eos, cond, method, config; kwarg...)
 end
 
-function flash_storage(eos, cond, method, config::FlashConfig{PrintOutput, true}; kwarg...) where PrintOutput
+function flash_storage(eos, cond, method, config::FlashConfig; kwarg...)
+    if !use_dict_storage(config)
+        return flash_storage(eos, cond, method, StaticConfig(); kwarg...)
+    end
     out = Dict{Symbol,Any}()
     d = flash_storage_internal!(out, eos, cond, method; kwarg...)
     # Convert to named tuple
     return NamedTuple(pairs(d))
 end
 
-"""
-    flash_storage_static(eos, cond, [method])
-
-Construct an inline, statically sized SSI workspace without using the `Dict`-based
-host storage builder. This constructor is suitable for use inside accelerator
-kernels when `eos` and `cond` are isbits values.
-"""
-function flash_storage(eos::GenericCubicEOS{E, R, N}, cond, method::SSIFlash, config::FlashConfig{PrintOutput, false}) where {E, R, N, PrintOutput}
-    T = eltype(cond.z)
-    return flash_storage(eos, cond, method, config, T)
-end
-
-@inline function flash_storage(eos::GenericCubicEOS{E, R, N}, cond, method::SSIFlash,
-        config::FlashConfig{PrintOutput, false}, ::Type{T}) where {E, R, N, PrintOutput, T}
-    forces = force_coefficients_static(eos, cond, T)
-    return (
-        forces = forces,
-        x = zero(MVector{N, T}),
-        y = zero(MVector{N, T}),
-        buffer1 = zero(MVector{N, T}),
-        buffer2 = zero(MVector{N, T})
-    )
-end
-
-
-flash_storage_static(eos, cond, method::SSIFlash = SSIFlash()) =
-    flash_storage(eos, cond, method, FlashConfig{false, false}())
-
-
-"""
-    flash_2ph_static(eos, cond, K, V[, maxiter, tolerance])
-
-Compatibility wrapper around the config-driven flash implementation using entirely
-local static storage and no host-side diagnostics.
-Returns `(V, K, iterations, converged)`.
-"""
-@inline flash_2ph_static(eos, cond, K, V) =
-    flash_2ph_static(eos, cond, K, V, 25000, 1e-8)
-
-@inline function flash_2ph_static(
-        eos::GenericCubicEOS{E, R, N}, cond, K::MVector{N, T}, V,
-        maxiter::Int, tolerance
-    ) where {E, R, N, T}
-    config = FlashConfig{false, false}()
-    storage = flash_storage(eos, cond, SSIFlash(), config, T)
-    V, K, report = flash_2ph_impl!(storage, K, eos, cond, V, config;
-        method = SSIFlash(),
-        maxiter = maxiter,
-        tolerance = tolerance,
-        verbose = false,
-        check = false,
-        z_min = nothing
-    )
-    return (V, K, report.its, report.converged)
-end
-
 function flash_storage_internal!(out, eos, cond, method;
         inc_jac = isa(method, AbstractNewtonFlash),
         inc_bypass = false,
-        static_size = false,
         kwarg...
     )
     n = number_of_components(eos)
-    alloc_forces(c) = force_coefficients(eos, c, static_size = static_size)
+    alloc_forces(c) = force_coefficients(eos, c)
     if forces_per_phase(eos)
         cond = set_phase(cond, :liquid)
         lforces = alloc_forces(cond)
@@ -295,47 +202,31 @@ function flash_storage_internal!(out, eos, cond, method;
     else
         out[:forces] = alloc_forces(cond)
     end
-    if static_size
-        # Do not write `@MVector zeros(n)` here.  That form first creates a
-        # regular Vector and then copies it into the MVector.  Apart from the
-        # unnecessary allocation, this is particularly inconvenient for GPU
-        # callers since the temporary is a host array.  Construct the
-        # tuple-backed static array directly instead.
-        alloc_vec = () -> zero(MVector{n, eltype(cond.z)})
-    else
-        alloc_vec = () -> zeros(n)
-    end
+    alloc_vec = () -> zeros(n)
     out[:x] = alloc_vec()
     out[:y] = alloc_vec()
 
     out[:buffer1] = alloc_vec()
     out[:buffer2] = alloc_vec()
     if inc_jac
-        flash_storage_internal_newton!(out, eos, cond, method, static_size = static_size; kwarg...)
+        flash_storage_internal_newton!(out, eos, cond, method; kwarg...)
     end
     if inc_bypass
-        out[:bypass] = michelsen_critical_point_measure_storage(eos, static_size = static_size)
+        out[:bypass] = michelsen_critical_point_measure_storage(eos, static_size = false)
     end
     return out
 end
 
-function flash_storage_internal_newton!(out, eos, cond, method; static_size = false, diff_externals = false, kwarg...)
+function flash_storage_internal_newton!(out, eos, cond, method; diff_externals = false, kwarg...)
     n = number_of_components(eos)
     np = 2*n + 1
     primary_ad(ix) = get_ad(0.0, np, typeof(ForwardDiff.Tag(Val(:Flash),Nothing)), ix)
     V_ad = primary_ad(np)
     T = typeof(V_ad)
-    if static_size
-        x_ad = zero(MVector{n, T})
-        y_ad = zero(MVector{n, T})
-        r = zero(MVector{np, Float64})
-        J = zero(MMatrix{np, np, Float64})
-    else
-        x_ad = zeros(T, n)
-        y_ad = zeros(T, n)
-        r = zeros(np)
-        J = zeros(np, np)
-    end
+    x_ad = zeros(T, n)
+    y_ad = zeros(T, n)
+    r = zeros(np)
+    J = zeros(np, np)
     out[:r] = r
     out[:J] = J
 
@@ -345,12 +236,12 @@ function flash_storage_internal_newton!(out, eos, cond, method; static_size = fa
     end
     out[:AD] = (x = x_ad, y = y_ad, V = V_ad)
     if diff_externals
-        flash_storage_internal_inverse!(out, eos, cond, method, static_size = static_size; kwarg...)
+        flash_storage_internal_inverse!(out, eos, cond, method; kwarg...)
     end
     return out
 end
 
-function flash_storage_internal_inverse!(out, eos, cond, method; static_size = false, npartials = nothing)
+function flash_storage_internal_inverse!(out, eos, cond, method; npartials = nothing)
     n = number_of_components(eos)
     np = length(out[:r])
     external_partials = n + 2 # p, T, z_1, ... z_n
@@ -358,28 +249,19 @@ function flash_storage_internal_inverse!(out, eos, cond, method; static_size = f
     p_ad = secondary_ad(1)
     T_ad = secondary_ad(2)
     T_cond = typeof(p_ad)
-    if static_size
-        z_ad = zero(MVector{n, T_cond})
-        J_inv = zero(MMatrix{np, external_partials, Float64})
-    else
-        z_ad = zeros(T_cond, n)
-        J_inv = zeros(np, external_partials)
-    end
+    z_ad = zeros(T_cond, n)
+    J_inv = zeros(np, external_partials)
     out[:J_inv] = J_inv
     for i = 1:n
         z_ad[i] = secondary_ad(i+2)
     end
     cond_ad = (p = p_ad, T = T_ad, z = z_ad, phase = :liquid)
     if !isnothing(npartials)
-        if static_size
-            buf = zero(MVector{npartials, Float64})
-        else
-            buf = zeros(npartials)
-        end
+        buf = zeros(npartials)
         out[:buf_inv] = buf
     end
     out[:AD_cond] = cond_ad
-    out[:forces_secondary] = force_coefficients(eos, cond_ad, static_size = static_size)
+    out[:forces_secondary] = force_coefficients(eos, cond_ad)
 end
 
 function get_ad(v::T, npartials, tag, diag_pos = nothing) where {T<:Real}
@@ -396,59 +278,19 @@ function update_value(v::T, newv::Real) where T<:ForwardDiff.Dual
 end
 
 function flash_update!(K, storage, type::SSIFlash, eos, cond, forces, V, iteration)
-    return flash_update!(K, storage, type, eos, cond, forces, V, iteration, FlashConfig())
-end
-
-function flash_update!(K, storage, type::SSIFlash, eos, cond, forces, V, iteration, config::FlashConfig)
     z = cond.z
     x, y = storage.x, storage.y
     p, T = cond.p, cond.T
-    return ssi!(K, p, T, x, y, z, V, eos, forces, config)
+    return ssi!(K, p, T, x, y, z, V, eos, forces)
 end
 
 function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces) where {F<:Real}
-    return ssi!(K, p, T, x, y, z, V, eos, forces, FlashConfig())
-end
-
-@inline function ssi!(K, p::F, T::F, x, y, z, V::F, eos, forces, config::FlashConfig) where {F<:Real}
-    liquid_phase = phase_value(config, Val(:liquid))
-    vapor_phase = phase_value(config, Val(:vapor))
-    return ssi_with_phases!(K, p, T, x, y, z, V, eos, forces, liquid_phase, vapor_phase)
-end
-
-@inline 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)
-end
-
-@inline function ssi_stack(K::SVector{N, F}, p::F, T::F, z, V::F,
-        eos, forces) where {N, F<:Real}
-    x = SVector{N, F}(ntuple(i -> liquid_mole_fraction(z[i], K[i], V), Val(N)))
-    y = SVector{N, F}(ntuple(i -> vapor_mole_fraction(x[i], K[i]), Val(N)))
-    liquid = (p = p, T = T, z = x, phase = Val(:liquid))
-    vapor = (p = p, T = T, z = y, phase = Val(:vapor))
-    f_l = static_fugacities(eos, liquid, forces, F)
-    f_v = static_fugacities(eos, vapor, forces, F)
-    ratios = SVector{N, F}(ntuple(i -> f_l[i]/f_v[i], Val(N)))
-    residual = zero(F)
-    @inbounds for i in 1:N
-        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)
-    return V_next, K_next, residual
-end
-
-function ssi_with_phases!(K, p::F, T::F, x, y, z, V::F, eos, forces, liquid_phase, vapor_phase) where {F<:Real}
     # Initialize conditions for vapor and liquid phases based on K-values
     x = liquid_mole_fraction!(x, z, K, V)
     y = vapor_mole_fraction!(y, x, K)
 
-    liquid = (p = p, T = T, z = x, phase = liquid_phase)
-    vapor = (p = p, T = T, z = y, phase = vapor_phase)
+    liquid = (p = p, T = T, z = x, phase = :liquid)
+    vapor = (p = p, T = T, z = y, phase = :vapor)
 
     Z_l, s_l = prep(eos, liquid, forces)
     Z_v, s_v = prep(eos, vapor, forces)
@@ -470,10 +312,6 @@ 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)
-    return flash_update!(K, storage, type, eos, cond, forces, V, iteration, FlashConfig())
-end
-
-function flash_update!(K, storage, type::NewtonFlash, eos, cond, forces, V, iteration, config::FlashConfig)
     x, y = storage.x, storage.y
     z = cond.z
     x = liquid_mole_fraction!(x, z, K, V)
@@ -566,14 +404,10 @@ function newton_dampen!(dMax, Δ)
 end
 
 function flash_update!(K, storage, type::SSINewtonFlash, eos, cond, forces, V, iteration)
-    return flash_update!(K, storage, type, eos, cond, forces, V, iteration, FlashConfig())
-end
-
-function flash_update!(K, storage, type::SSINewtonFlash, eos, cond, forces, V, iteration, config::FlashConfig)
     if iteration >= type.swap_iter
-        flash_update!(K, storage, NewtonFlash(type.dMax), eos, cond, forces, V, iteration, config)
+        flash_update!(K, storage, NewtonFlash(type.dMax), eos, cond, forces, V, iteration)
     else
-        flash_update!(K, storage, SSIFlash(), eos, cond, forces, V, iteration, config)
+        flash_update!(K, storage, SSIFlash(), eos, cond, forces, V, iteration)
     end
 end
 
diff --git a/src/flash_types.jl b/src/flash_types.jl
index 036346e..23fd6e5 100644
--- a/src/flash_types.jl
+++ b/src/flash_types.jl
@@ -4,24 +4,19 @@ abstract type AbstractFlash end
 abstract type AbstractNewtonFlash <: AbstractFlash end
 
 """
-    FlashConfig{PrintOutput, UseDictStorage}()
     FlashConfig(; print_output=true, use_dict_storage=true)
 
-Compile-time configuration for the flash call graph. `PrintOutput=false` removes
-logging, warnings, and checked error paths from accelerator kernels.
-`UseDictStorage=false` selects direct static storage construction instead of the
-host-oriented `Dict` builder.
+Runtime options for the ordinary flash implementation. Set `print_output=false`
+to suppress diagnostics. Setting `use_dict_storage=false` selects the fully
+static implementation; new code should prefer `flash_storage(...; static=true)`.
 """
-struct FlashConfig{PrintOutput, UseDictStorage} end
-
-FlashConfig(; print_output::Bool = true, use_dict_storage::Bool = true) =
-    FlashConfig{print_output, use_dict_storage}()
-
-@inline print_output(::FlashConfig{PrintOutput}) where PrintOutput = PrintOutput
-@inline use_dict_storage(::FlashConfig{PrintOutput, UseDictStorage}) where {PrintOutput, UseDictStorage} = UseDictStorage
+Base.@kwdef struct FlashConfig
+    print_output::Bool = true
+    use_dict_storage::Bool = true
+end
 
-@inline phase_value(::FlashConfig{PrintOutput, true}, ::Val{Phase}) where {PrintOutput, Phase} = Phase
-@inline phase_value(::FlashConfig{PrintOutput, false}, phase::Val) where PrintOutput = phase
+@inline print_output(config::FlashConfig) = config.print_output
+@inline use_dict_storage(config::FlashConfig) = config.use_dict_storage
 
 """
     ssi = SSIFlash()
diff --git a/src/kvalues.jl b/src/kvalues.jl
index b42f81e..8d164cd 100644
--- a/src/kvalues.jl
+++ b/src/kvalues.jl
@@ -61,21 +61,3 @@ end
 In-place version of `initial_guess_K`.
 """
 initial_guess_K!(K, eos, cond) = wilson_estimate!(K, eos, cond.p, cond.T)
-
-"""Return inline Wilson K-values for an EOS with a compile-time component count."""
-function initial_guess_K_static(eos::GenericCubicEOS{E, R, N}, cond) where {E, R, N}
-    T = Base.promote_eltype(cond.p, cond.T, cond.z[1])
-    return initial_guess_K_static(eos, cond, T)
-end
-
-@inline function initial_guess_K_static(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T}) where {E, R, N, T}
-    properties = eos.mixture.properties
-    return MVector{N, T}(ntuple(i -> wilson_estimate(properties[i], cond.p, cond.T), Val(N)))
-end
-
-@inline function initial_guess_K_static(eos::GenericCubicEOS{E, R, N}, cond,
-        ::FlashConfig{false, false}) where {E, R, N}
-    T = Base.promote_eltype(cond.p, cond.T, cond.z[1])
-    properties = eos.mixture.properties
-    return SVector{N, T}(ntuple(i -> wilson_estimate(properties[i], cond.p, cond.T), Val(N)))
-end
diff --git a/src/kvalues_eos.jl b/src/kvalues_eos.jl
index 33d3e2b..62307f3 100644
--- a/src/kvalues_eos.jl
+++ b/src/kvalues_eos.jl
@@ -18,8 +18,7 @@ function flash_storage(eos::KValuesEOS, cond = missing; kwarg...)
     return nothing
 end
 
-flash_storage(eos::KValuesEOS, cond, method, config::FlashConfig{PrintOutput, true}) where PrintOutput = nothing
-flash_storage(eos::KValuesEOS, cond, method, config::FlashConfig{PrintOutput, false}) where PrintOutput = nothing
+flash_storage(eos::KValuesEOS, cond, method, config::FlashConfig) = nothing
 
 function flash_2ph!(storage, K, eos::KValuesEOS, cond, V = NaN; kwarg...)
     return solve_rachford_rice(K, cond.z, V)
@@ -27,6 +26,6 @@ end
 
 
 function flash_2ph!(storage, K, eos::KValuesEOS, cond, V,
-        config::FlashConfig{PrintOutput, UseDictStorage}; kwarg...) where {PrintOutput, UseDictStorage}
+        config::FlashConfig; kwarg...)
     return solve_rachford_rice(K, cond.z, V)
 end
diff --git a/src/mixture_types.jl b/src/mixture_types.jl
index 4b1d751..ad1a667 100644
--- a/src/mixture_types.jl
+++ b/src/mixture_types.jl
@@ -81,21 +81,6 @@ struct MultiComponentMixture{R, N, Name, Names, BIC}
     end
 end
 
-"""
-    static_mixture(mixture)
-
-Return an isbits representation of a mixture for accelerator kernels. Human-readable
-metadata is omitted and binary-interaction coefficients are stored inline.
-"""
-function static_mixture(mixture::MultiComponentMixture{R, N}) where {R, N}
-    names = ntuple(_ -> nothing, Val(N))
-    bic = mixture.binary_interaction
-    if !isnothing(bic)
-        bic = SMatrix{N, N, R}(bic)
-    end
-    return MultiComponentMixture(mixture.properties; A_ij = bic, names = names, name = nothing)
-end
-
 """
     MultiComponentMixture(names::Vector{String}; A_ij = nothing, name = "UnnamedMixture")
 
diff --git a/src/rachford_rice.jl b/src/rachford_rice.jl
index 1384a9e..9120e70 100644
--- a/src/rachford_rice.jl
+++ b/src/rachford_rice.jl
@@ -20,79 +20,6 @@ julia> solve_rachford_rice([0.5, 1.5], [0.3, 0.7])
 0.8000000000000002
 ```
 """
-@inline function solve_rachford_rice(K::StaticVector{2}, z::StaticVector{2}, V = NaN)
-    z1, z2 = z
-    k1, k2 = K
-    b1, b2 = inv(1 - k1), inv(1 - k2)
-    return (z1*b2 + z2*b1)/(z1 + z2)
-end
-
-@inline function solve_rachford_rice(K::StaticVector{3}, z::StaticVector{3}, V = NaN)
-    z1, z2, z3 = z
-    k1, k2, k3 = K
-    b1, b2, b3 = inv(1-k1), inv(1-k2), inv(1-k3)
-    a2 = z1 + z2 + z3
-    a1 = -b1*(z2 + z3) - b2*(z1 + z3) - b3*(z1 + z2)
-    a0 = b1*b2*z3 + b1*b3*z2 + b2*b3*z1
-    discriminant = a1*a1 - 4*a0*a2
-    if discriminant >= zero(discriminant)
-        inv_2a2 = inv(2*a2)
-        root_offset = sqrt(discriminant)*inv_2a2
-        root_center = -a1*inv_2a2
-        root1 = root_center - root_offset
-        root2 = root_center + root_offset
-        if zero(root1) < root1 < one(root1)
-            return root1
-        elseif zero(root2) < root2 < one(root2)
-            return root2
-        elseif isfinite(root1 + root2)
-            kmin = min(k1, k2, k3)
-            kmax = max(k1, k2, k3)
-            kmin > one(kmin) && return max(root1, root2)
-            kmax < one(kmax) && return min(root1, root2)
-        end
-    end
-    return solve_rachford_rice_static_iterative(K, z, V)
-end
-
-@inline function solve_rachford_rice(K::StaticVector, z::StaticVector, V = NaN)
-    return solve_rachford_rice_static_iterative(K, z, V)
-end
-
-"""GPU-compatible Rachford--Rice fallback for statically sized inputs."""
-@inline function solve_rachford_rice_static_iterative(K, z, V; tol = 1e-12, maxiter = 1000)
-    V_lo = inv(1 - maximum(K))
-    V_hi = inv(1 - minimum(K))
-    if V_hi < V_lo
-        V_lo, V_hi = V_hi, V_lo
-    end
-    if isnan(V)
-        V = (V_lo + V_hi)/2
-    end
-    for _ in 1:maxiter
-        residual = zero(V)
-        denominator = zero(V)
-        @inbounds for i in eachindex(K)
-            delta_K = K[i] - one(K[i])
-            term_denominator = one(V) + V*delta_K
-            residual += z[i]*delta_K/term_denominator
-            denominator += z[i]*delta_K^2/term_denominator^2
-        end
-        abs(residual) < tol && break
-        if residual > zero(residual)
-            V_lo = V
-        else
-            V_hi = V
-        end
-        V_next = V + residual/denominator
-        if !(V_lo < V_next < V_hi) || !isfinite(V_next)
-            V_next = (V_lo + V_hi)/2
-        end
-        V = V_next
-    end
-    return V
-end
-
 function solve_rachford_rice(K, z, V = NaN; tol=1e-12, maxiter=1000, ad = false, analytical = true, verbose = false)
     n = length(z)
     if analytical
diff --git a/src/stability.jl b/src/stability.jl
index 5e7479e..3433a22 100644
--- a/src/stability.jl
+++ b/src/stability.jl
@@ -12,6 +12,9 @@ function stability_2ph(eos, c, K = initial_guess_K(eos, c); kwarg...)
 end
 
 function stability_2ph(eos, c, K, config::FlashConfig; kwarg...)
+    if !use_dict_storage(config)
+        return stability_2ph(eos, c, K, StaticConfig(); kwarg...)
+    end
     storage = flash_storage(eos, c, SSIFlash(), config)
     stability_2ph!(storage, K, eos, c, config; kwarg...)
 end
@@ -32,12 +35,10 @@ function stability_2ph!(storage, K, eos, c, config::FlashConfig;
     f_xy = storage.buffer2
     x, y = storage.x, storage.y
     z, p, T = c.z, c.p, c.T
-    liquid_phase = phase_value(config, Val(:liquid))
-    vapor_phase = phase_value(config, Val(:vapor))
-    liquid = (p = p, T = T, z = x, phase = liquid_phase)
-    vapor = (p = p, T = T, z = y, phase = vapor_phase)
-    current_as_liquid = (p = p, T = T, z = z, phase = liquid_phase)
-    current_as_vapor = (p = p, T = T, z = z, phase = vapor_phase)
+    liquid = (p = p, T = T, z = x, phase = :liquid)
+    vapor = (p = p, T = T, z = y, phase = :vapor)
+    current_as_liquid = (p = p, T = T, z = z, phase = :liquid)
+    current_as_vapor = (p = p, T = T, z = z, phase = :vapor)
     mixture_fugacities!(f_z, eos, current_as_vapor, forces)
     if check_vapor
         wilson_estimate!(K, eos, p, T)
@@ -85,109 +86,6 @@ f_ratio(f_z, f_xy, ::Val{false}) = f_xy/f_z
 xy_value(z, K, ::Val{true}) = z*K
 xy_value(z, K, ::Val{false}) = z/K
 
-stability_phase(::Val{true}) = Val(:vapor)
-stability_phase(::Val{false}) = Val(:liquid)
-
-@generated function stability_xy(z::SVector{N, F}, K::SVector{N, F}, phase) where {N, F}
-    values = [:(xy_value(z[$i], K[$i], phase)) for i in 1:N]
-    return :(SVector{N, F}(($(values...),)))
-end
-
-@generated function static_scale(v::SVector{N, F}, scale::F) where {N, F}
-    values = [:(v[$i]/scale) for i in 1:N]
-    return :(SVector{N, F}(($(values...),)))
-end
-
-@generated function stability_ratios(f_z::SVector{N, F}, f_xy::SVector{N, F},
-        scale::F, phase) where {N, F}
-    values = [:(f_ratio(f_z[$i], scale*f_xy[$i], phase)) for i in 1:N]
-    return :(SVector{N, F}(($(values...),)))
-end
-
-@generated function static_multiply(a::SVector{N, F}, b::SVector{N, F}) where {N, F}
-    values = [:(a[$i]*b[$i]) for i in 1:N]
-    return :(SVector{N, F}(($(values...),)))
-end
-
-@generated function static_divide(a::SVector{N, F}, b::SVector{N, F}) where {N, F}
-    values = [:(a[$i]/b[$i]) for i in 1:N]
-    return :(SVector{N, F}(($(values...),)))
-end
-
-@inline function michelsen_test_stack(f_z, z::SVector{N, F}, K::SVector{N, F},
-        eos, cond, forces, inside_is_vapor;
-        tol_equil = 1e-10,
-        tol_trivial = tol_equil,
-        tol_sat = tol_trivial,
-        maxiter = 1000
-    ) where {N, F}
-    trivial = false
-    S = one(F)
-    iter = 0
-    xy = zero(SVector{N, F})
-    while true
-        iter += 1
-        unnormalized = stability_xy(z, K, inside_is_vapor)
-        S = zero(F)
-        @inbounds for i in 1:N
-            S += unnormalized[i]
-        end
-        xy = static_scale(unnormalized, S)
-        inside = (p = cond.p, T = cond.T, z = xy,
-            phase = stability_phase(inside_is_vapor))
-        f_xy = static_fugacities(eos, inside, forces, F)
-        ratios = stability_ratios(f_z, f_xy, S, inside_is_vapor)
-        K = static_multiply(K, ratios)
-        R_norm = zero(F)
-        K_norm = zero(F)
-        @inbounds for i in 1:N
-            R_norm += (ratios[i] - one(F))^2
-            K_norm += log(K[i])^2
-        end
-        trivial = K_norm < tol_trivial
-        converged = R_norm < tol_equil
-        if trivial || converged
-            break
-        elseif iter == maxiter
-            # Match the checked host path: failed stability iterations are
-            # conservatively treated as a trivial, stable solution.
-            trivial = true
-            break
-        end
-    end
-    stable = trivial || S <= one(F) + tol_sat
-    return stable, trivial, iter, K, xy
-end
-
-@inline function stability_2ph_stack(K::SVector{N, F}, eos, cond, forces;
-        check_vapor::Bool = true,
-        check_liquid::Bool = true,
-        kwarg...
-    ) where {N, F}
-    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_static(eos, cond, FlashConfig{false, false}())
-    if check_vapor
-        stable_vapor, trivial_vapor, i_v, K_vapor, y = michelsen_test_stack(
-            f_z_vapor, cond.z, K_wilson, eos, cond, forces, Val(true); kwarg...)
-    else
-        stable_vapor, trivial_vapor, i_v, K_vapor, y = true, true, 0, K_wilson, cond.z
-    end
-    if check_liquid
-        liquid_phase = (p = cond.p, T = cond.T, z = cond.z, phase = Val(:liquid))
-        f_z_liquid = forces_per_phase(eos) ?
-            static_fugacities(eos, liquid_phase, forces, F) : f_z_vapor
-        stable_liquid, trivial_liquid, i_l, K_liquid, x = michelsen_test_stack(
-            f_z_liquid, cond.z, K_wilson, eos, cond, forces, Val(false); kwarg...)
-    else
-        stable_liquid, trivial_liquid, i_l, K_liquid, x = true, true, 0, K_wilson, cond.z
-    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
-end
-
 """
     stability_2ph!(storage, eos, c, [K])
 
diff --git a/src/static.jl b/src/static.jl
new file mode 100644
index 0000000..904297d
--- /dev/null
+++ b/src/static.jl
@@ -0,0 +1,357 @@
+"""
+    StaticConfig()
+
+Marker returned by `flash_storage(...; static=true)`. It selects the immutable,
+stack-oriented SSI implementation used in accelerator kernels.
+"""
+struct StaticConfig end
+
+@inline print_output(::StaticConfig) = false
+@inline use_dict_storage(::StaticConfig) = false
+
+function flash_storage(eos::GenericCubicEOS, cond, method, config::StaticConfig; kwarg...)
+    method isa SSIFlash || throw(ArgumentError("The static flash currently supports SSIFlash only."))
+    return config
+end
+
+"""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))
+    bic = mixture.binary_interaction
+    if !isnothing(bic)
+        bic = SMatrix{N, N, R}(bic)
+    end
+    return MultiComponentMixture(mixture.properties; A_ij = bic, names = names, name = nothing)
+end
+
+"""
+    static_eos(eos)
+
+Convert a generic cubic EOS to an isbits representation for accelerator kernels.
+"""
+function static_eos(eos::GenericCubicEOS{T, R, N}) where {T, R, N}
+    mixture = static_mixture(eos.mixture)
+    volume_shift = eos.volume_shift
+    if !isnothing(volume_shift)
+        volume_shift = SVector{N, eltype(volume_shift)}(volume_shift)
+    end
+    return GenericCubicEOS(
+        eos.type,
+        mixture,
+        eos.m_1,
+        eos.m_2,
+        eos.ω_a,
+        eos.ω_b,
+        volume_shift
+    )
+end
+
+"""Return immutable Wilson K-values for static storage."""
+@inline function initial_guess_K(eos::GenericCubicEOS{E, R, N}, cond,
+        ::StaticConfig) where {E, R, N}
+    T = Base.promote_eltype(cond.p, cond.T, cond.z[1])
+    properties = eos.mixture.properties
+    return SVector{N, T}(ntuple(i -> wilson_estimate(properties[i], cond.p, cond.T), Val(N)))
+end
+
+# Val phase tags keep Symbol construction and dynamic dispatch out of kernels.
+@inline phase_symbol(::Val{phase}) where phase = phase
+
+@inline function pick_root(eos, roots, cond, forces, scalars, ::Val{:liquid})
+    min_root, _ = root_bounds(roots, minimum_allowable_root(eos, forces, scalars))
+    return min_root
+end
+
+@inline function pick_root(eos, roots, cond, forces, scalars, ::Val{:vapor})
+    _, max_root = root_bounds(roots, minimum_allowable_root(eos, forces, scalars))
+    return max_root
+end
+
+@inline get_force_coefficients(forces, eos::GenericCubicEOS, cond) = forces
+
+"""Immutable force coefficients for accelerator kernels."""
+@inline function static_force_coefficients(eos::GenericCubicEOS{E, R, N}, cond,
+        ::Type{T}) where {E, R, N, T}
+    A_i_static = SVector{N, T}(ntuple(i -> A_i(eos, cond, i), Val(N)))
+    B_i_static = SVector{N, T}(ntuple(i -> B_i(eos, cond, i), Val(N)))
+    A_ij_static = SMatrix{N, N, T}(ntuple(Val(N*N)) do index
+        i = mod1(index, N)
+        j = (index - 1) ÷ N + 1
+        sqrt(A_i_static[i]*A_i_static[j]) *
+            (one(T) - binary_interaction(eos, i, j, cond))
+    end)
+    return (A_ij = A_ij_static, A_i = A_i_static, B_i = B_i_static)
+end
+
+@inline function solve_rachford_rice(K::StaticVector{2}, z::StaticVector{2}, V = NaN)
+    z1, z2 = z
+    k1, k2 = K
+    b1, b2 = inv(1 - k1), inv(1 - k2)
+    return (z1*b2 + z2*b1)/(z1 + z2)
+end
+
+@inline function solve_rachford_rice(K::StaticVector{3}, z::StaticVector{3}, V = NaN)
+    z1, z2, z3 = z
+    k1, k2, k3 = K
+    b1, b2, b3 = inv(1-k1), inv(1-k2), inv(1-k3)
+    a2 = z1 + z2 + z3
+    a1 = -b1*(z2 + z3) - b2*(z1 + z3) - b3*(z1 + z2)
+    a0 = b1*b2*z3 + b1*b3*z2 + b2*b3*z1
+    discriminant = a1*a1 - 4*a0*a2
+    if discriminant >= zero(discriminant)
+        inv_2a2 = inv(2*a2)
+        root_offset = sqrt(discriminant)*inv_2a2
+        root_center = -a1*inv_2a2
+        root1 = root_center - root_offset
+        root2 = root_center + root_offset
+        if zero(root1) < root1 < one(root1)
+            return root1
+        elseif zero(root2) < root2 < one(root2)
+            return root2
+        elseif isfinite(root1 + root2)
+            kmin = min(k1, k2, k3)
+            kmax = max(k1, k2, k3)
+            kmin > one(kmin) && return max(root1, root2)
+            kmax < one(kmax) && return min(root1, root2)
+        end
+    end
+    return solve_rachford_rice_static_iterative(K, z, V)
+end
+
+@inline solve_rachford_rice(K::StaticVector, z::StaticVector, V = NaN) =
+    solve_rachford_rice_static_iterative(K, z, V)
+
+@inline function solve_rachford_rice_static_iterative(K, z, V;
+        tol = 1e-12, maxiter = 1000)
+    V_lo = inv(1 - maximum(K))
+    V_hi = inv(1 - minimum(K))
+    if V_hi < V_lo
+        V_lo, V_hi = V_hi, V_lo
+    end
+    if isnan(V)
+        V = (V_lo + V_hi)/2
+    end
+    for _ in 1:maxiter
+        residual = zero(V)
+        denominator = zero(V)
+        @inbounds for i in eachindex(K)
+            delta_K = K[i] - one(K[i])
+            term_denominator = one(V) + V*delta_K
+            residual += z[i]*delta_K/term_denominator
+            denominator += z[i]*delta_K^2/term_denominator^2
+        end
+        abs(residual) < tol && break
+        if residual > zero(residual)
+            V_lo = V
+        else
+            V_hi = V
+        end
+        V_next = V + residual/denominator
+        if !(V_lo < V_next < V_hi) || !isfinite(V_next)
+            V_next = (V_lo + V_hi)/2
+        end
+        V = V_next
+    end
+    return V
+end
+
+@inline 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)
+end
+
+@inline function static_ssi(K::SVector{N, F}, p::F, T::F, z, V::F,
+        eos, forces) where {N, F<:Real}
+    x = SVector{N, F}(ntuple(i -> liquid_mole_fraction(z[i], K[i], V), Val(N)))
+    y = SVector{N, F}(ntuple(i -> vapor_mole_fraction(x[i], K[i]), Val(N)))
+    liquid = (p = p, T = T, z = x, phase = Val(:liquid))
+    vapor = (p = p, T = T, z = y, phase = Val(:vapor))
+    f_l = static_fugacities(eos, liquid, forces, F)
+    f_v = static_fugacities(eos, vapor, forces, F)
+    ratios = SVector{N, F}(ntuple(i -> f_l[i]/f_v[i], Val(N)))
+    residual = zero(F)
+    @inbounds for i in 1:N
+        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)
+    return V_next, K_next, residual
+end
+
+@inline function flash_2ph(eos::GenericCubicEOS, c, K, V,
+        config::StaticConfig; kwarg...)
+    return flash_2ph!(config, K, eos, c, V; kwarg...)
+end
+
+@inline function flash_2ph!(config::StaticConfig, K, eos::GenericCubicEOS, c,
+        V = NaN; extra_out::Bool = false, kwarg...)
+    out = flash_2ph_impl!(config, K, eos, c, V; kwarg...)
+    return static_flash_output(out, Val(extra_out))
+end
+
+@inline static_flash_output(out, ::Val{true}) = out
+@inline static_flash_output(out, ::Val{false}) = out[1]
+
+@inline function flash_2ph_impl!(::StaticConfig, K,
+        eos::GenericCubicEOS{E, R, N}, c, V;
+        method::SSIFlash = SSIFlash(),
+        maxiter::Int = 25000,
+        tolerance::Float64 = 1e-8,
+        verbose::Bool = false,
+        check::Bool = true,
+        update_forces::Bool = true,
+        z_min = MINIMUM_COMPOSITION,
+        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)
+    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...)
+    else
+        stable = false
+        stability_report = StabilityReport(false, false, false, false)
+    end
+    converged = false
+    if stable
+        iteration = 0
+    else
+        iteration = 1
+        if isnan(V)
+            V = solve_rachford_rice(K, z, V)
+        end
+        while true
+            V, K, residual = static_ssi(K, cond.p, cond.T, z, V, eos, forces)
+            converged = residual <= tolerance
+            (converged || iteration == maxiter) && break
+            iteration += 1
+        end
+    end
+    report = (its = iteration, converged = converged, stability = stability_report)
+    return V, K, report
+end
+
+stability_phase(::Val{true}) = Val(:vapor)
+stability_phase(::Val{false}) = Val(:liquid)
+
+@generated function stability_xy(z::SVector{N, F}, K::SVector{N, F}, phase) where {N, F}
+    values = [:(xy_value(z[$i], K[$i], phase)) for i in 1:N]
+    return :(SVector{N, F}(($(values...),)))
+end
+
+@generated function static_scale(v::SVector{N, F}, scale::F) where {N, F}
+    values = [:(v[$i]/scale) for i in 1:N]
+    return :(SVector{N, F}(($(values...),)))
+end
+
+@generated function stability_ratios(f_z::SVector{N, F}, f_xy::SVector{N, F},
+        scale::F, phase) where {N, F}
+    values = [:(f_ratio(f_z[$i], scale*f_xy[$i], phase)) for i in 1:N]
+    return :(SVector{N, F}(($(values...),)))
+end
+
+@generated function static_multiply(a::SVector{N, F}, b::SVector{N, F}) where {N, F}
+    values = [:(a[$i]*b[$i]) for i in 1:N]
+    return :(SVector{N, F}(($(values...),)))
+end
+
+@generated function static_divide(a::SVector{N, F}, b::SVector{N, F}) where {N, F}
+    values = [:(a[$i]/b[$i]) for i in 1:N]
+    return :(SVector{N, F}(($(values...),)))
+end
+
+@inline function static_michelsen_test(f_z, z::SVector{N, F}, K::SVector{N, F},
+        eos, cond, forces, inside_is_vapor;
+        tol_equil = 1e-10,
+        tol_trivial = tol_equil,
+        tol_sat = tol_trivial,
+        maxiter = 1000
+    ) where {N, F}
+    trivial = false
+    S = one(F)
+    iter = 0
+    xy = zero(SVector{N, F})
+    while true
+        iter += 1
+        unnormalized = stability_xy(z, K, inside_is_vapor)
+        S = zero(F)
+        @inbounds for i in 1:N
+            S += unnormalized[i]
+        end
+        xy = static_scale(unnormalized, S)
+        inside = (p = cond.p, T = cond.T, z = xy,
+            phase = stability_phase(inside_is_vapor))
+        f_xy = static_fugacities(eos, inside, forces, F)
+        ratios = stability_ratios(f_z, f_xy, S, inside_is_vapor)
+        K = static_multiply(K, ratios)
+        R_norm = zero(F)
+        K_norm = zero(F)
+        @inbounds for i in 1:N
+            R_norm += (ratios[i] - one(F))^2
+            K_norm += log(K[i])^2
+        end
+        trivial = K_norm < tol_trivial
+        converged = R_norm < tol_equil
+        if trivial || converged
+            break
+        elseif iter == maxiter
+            trivial = true
+            break
+        end
+    end
+    stable = trivial || S <= one(F) + tol_sat
+    return stable, trivial, iter, K, xy
+end
+
+@inline function static_stability_2ph(K::SVector{N, F}, eos, cond, forces;
+        check_vapor::Bool = true,
+        check_liquid::Bool = true,
+        kwarg...
+    ) where {N, F}
+    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())
+    if check_vapor
+        stable_vapor, trivial_vapor, i_v, K_vapor, y = static_michelsen_test(
+            f_z_vapor, cond.z, K_wilson, eos, cond, forces, Val(true); kwarg...)
+    else
+        stable_vapor, trivial_vapor, i_v, K_vapor, y = true, true, 0, K_wilson, cond.z
+    end
+    if check_liquid
+        liquid_phase = (p = cond.p, T = cond.T, z = cond.z, phase = Val(:liquid))
+        f_z_liquid = forces_per_phase(eos) ?
+            static_fugacities(eos, liquid_phase, forces, F) : f_z_vapor
+        stable_liquid, trivial_liquid, i_l, K_liquid, x = static_michelsen_test(
+            f_z_liquid, cond.z, K_wilson, eos, cond, forces, Val(false); kwarg...)
+    else
+        stable_liquid, trivial_liquid, i_l, K_liquid, x = true, true, 0, K_wilson, cond.z
+    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
+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))
+    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
+end
+
+@inline stability_2ph!(::StaticConfig, K, eos::GenericCubicEOS, c; kwarg...) =
+    stability_2ph(eos, c, K, StaticConfig(); kwarg...)
diff --git a/test/runtests.jl b/test/runtests.jl
index eda6d75..2fc81bb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -38,44 +38,37 @@ end
     end
 end
 @testset "Zero allocating flash" begin
-    for m in flash_methods
-        name = typeof(m)
-        @testset "$name - Arrays" begin
-            test_flash_inplace(m, static_size = true)
-        end
-        @testset "$name - StaticArrays" begin
-            test_flash_inplace(m, static_size = true)
-        end
-    end
+    test_flash_inplace(SSIFlash())
 end
 
 @testset "Static accelerator path" begin
     eos = static_eos(get_test_eos())
     c = (p = 1e6, T = 300.0, z = @SVector [0.5, 0.3, 0.2])
-    K = zero(MVector{3, Float64})
-    initial_guess_K!(K, eos, c)
-    V, K, iterations, converged = flash_2ph_static(eos, c, K, 0.5)
+    storage = flash_storage(eos, c; method = SSIFlash(), static = true)
+    K = initial_guess_K(eos, c, storage)
+    V, K, report = flash_2ph!(storage, K, eos, c, 0.5; extra_out = true)
 
     @test isbitstype(typeof(eos))
-    @test converged
-    @test iterations == 6
+    @test isbitstype(typeof(storage))
+    @test storage isa StaticConfig
+    @test report.converged
+    @test report.its == 6
     @test V ≈ 0.7632068334421974
     @test K ≈ @SVector [4.553402802323027, 17.73895830809456, 0.0004031451448211194]
 
-    K_static = initial_guess_K_static(eos, c)
-    V_static, K_static, iterations_static, converged_static = flash_2ph_static(eos, c, K_static, 0.5)
-    @test converged_static
-    @test iterations_static == iterations
-    @test V_static ≈ V
-    @test K_static ≈ K
+    normal_config = FlashConfig(print_output=false)
+    @test typeof(normal_config) == FlashConfig
+    @test flash_storage(eos, c, SSIFlash(), normal_config).x isa Vector
+    @test_throws ArgumentError flash_storage(eos, c; static_size = true)
 
     config = FlashConfig(print_output=false, use_dict_storage=false)
-    storage = flash_storage(eos, c, SSIFlash(), config)
-    K_config = initial_guess_K_static(eos, c)
-    V_config, K_config, report = flash_2ph!(storage, K_config, eos, c, NaN, config;
+    @test typeof(config) == FlashConfig
+    @test flash_storage(eos, c, SSIFlash(), config) isa StaticConfig
+    K_config = initial_guess_K(eos, c, storage)
+    V_config, K_config, config_report = flash_2ph(eos, c, K_config, NaN, config;
         method=SSIFlash(), extra_out=true, z_min=nothing)
-    @test report.stability isa MultiComponentFlash.StabilityReport
-    @test report.converged
+    @test config_report.stability isa MultiComponentFlash.StabilityReport
+    @test config_report.converged
     @test V_config ≈ V
     @test K_config ≈ K
 end
@@ -149,13 +142,12 @@ end
 
 using StaticArrays, KernelAbstractions, JLArrays
 @testset "Static flash with KernelAbstractions/JLArrays" begin
-    @kernel function static_flash_kernel!(out, pressure, temperature, z, eos)
+    @kernel function static_flash_kernel!(out, pressure, temperature, z, eos, storage)
         i = @index(Global)
         if i <= length(out)
             @inbounds cond = (p = pressure[i], T = temperature[i], z = z)
-            config = FlashConfig{false, false}()
-            K = initial_guess_K_static(eos, cond, config)
-            V = flash_2ph(eos, cond, K, NaN, config;
+            K = initial_guess_K(eos, cond, storage)
+            V = flash_2ph!(storage, K, eos, cond, NaN;
                 method = SSIFlash(), check = false, verbose = false, z_min = nothing)
             @inbounds out[i] = V
         end
@@ -164,6 +156,8 @@ using StaticArrays, KernelAbstractions, JLArrays
         host_eos = get_test_eos()
         eos = static_eos(host_eos)
         z = @SVector [0.5, 0.3, 0.2]
+        storage = flash_storage(eos, (p = 1e5, T = 300.0, z = z);
+            method = SSIFlash(), static = true)
         n = 16
         pressure_host = collect(range(1e5, 4e6, length = n))
         temperature_host = collect(range(280.0, 320.0, length = n))
@@ -177,7 +171,7 @@ using StaticArrays, KernelAbstractions, JLArrays
         out = JLArray(zeros(n))
         backend = JLArrays.JLBackend()
         kernel! = static_flash_kernel!(backend, 8)
-        kernel!(out, pressure, temperature, z, eos; ndrange = n)
+        kernel!(out, pressure, temperature, z, eos, storage; ndrange = n)
 
         @test Array(out) ≈ expected rtol = 1e-11
     else
diff --git a/test/test_setup.jl b/test/test_setup.jl
index 6fb36c9..0ec3b42 100644
--- a/test/test_setup.jl
+++ b/test/test_setup.jl
@@ -13,15 +13,11 @@ 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)
 
-function test_flash_inplace(m, do_test = true; static_size = false)
+function test_flash_inplace(m, do_test = true)
     eos = get_test_eos()
     c = test_conditions()
-    S = flash_storage(eos, c, method = m, static_size = static_size)
+    S = flash_storage(eos, c, method = m)
     K = initial_guess_K(eos, c)
-    if static_size
-        n = number_of_components(eos)
-        K = MVector{n}(K)
-    end
     # Just in case of compilation
     test_allocs(S, K, eos, c, m)
     # Then evaluation

From e45df03377396087f9645876ec2d389225d12928 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Sat, 12 Sep 2026 12:27:50 +0200
Subject: [PATCH 7/9] Add immutable flash variant

---
 src/MultiComponentFlash.jl |  2 +-
 src/static.jl              | 26 ++++++++++++++++++++++++++
 test/runtests.jl           |  8 ++++++++
 3 files changed, 35 insertions(+), 1 deletion(-)

diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl
index 218cb14..15d02fc 100644
--- a/src/MultiComponentFlash.jl
+++ b/src/MultiComponentFlash.jl
@@ -19,7 +19,7 @@ module MultiComponentFlash
     export KValuesEOS
     export number_of_components
     # Flash interfaces
-    export flash_2ph, flash_2ph!, flash_storage
+    export flash_2ph, flash_2ph!, flash_2ph_immutable, flash_storage
     export stability_2ph, stability_2ph!
     # Algorithms for flash
     export SSIFlash, NewtonFlash, SSINewtonFlash, FlashConfig, StaticConfig
diff --git a/src/static.jl b/src/static.jl
index 904297d..ae8ceb6 100644
--- a/src/static.jl
+++ b/src/static.jl
@@ -14,6 +14,32 @@ function flash_storage(eos::GenericCubicEOS, cond, method, config::StaticConfig;
     return config
 end
 
+"""
+    V, K = flash_2ph_immutable(eos, c[, storage]; <keyword arguments>)
+
+Run the immutable, accelerator-friendly two-phase flash implementation.
+`c.z` must be an `SVector`; `K` is returned as an `SVector` and `V` is the
+scalar vapor fraction. When `storage` is omitted, a static storage marker is
+created automatically.
+
+The immutable path currently supports `SSIFlash` and generic cubic EOS values
+converted with [`static_eos`](@ref).
+"""
+@inline function flash_2ph_immutable(eos, c; method = SSIFlash(), kwarg...)
+    return flash_2ph_immutable(eos, c,
+        flash_storage(eos, c; method = method, static = true);
+        method = method, kwarg...)
+end
+
+@inline function flash_2ph_immutable(eos, c, storage::StaticConfig;
+        method = SSIFlash(), 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
+end
+
 """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))
diff --git a/test/runtests.jl b/test/runtests.jl
index 2fc81bb..6519621 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -71,6 +71,14 @@ end
     @test config_report.converged
     @test V_config ≈ V
     @test K_config ≈ K
+
+    V_immutable, K_immutable = flash_2ph_immutable(eos, c)
+    @test V_immutable ≈ V
+    @test K_immutable ≈ K
+    @test K_immutable isa SVector{3, Float64}
+    @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)))
 end
 
 @testset "Partial derivatives" begin

From b0d965c125185f15a14eabd5a29a76fbedabc2ae Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Sat, 12 Sep 2026 13:01:27 +0200
Subject: [PATCH 8/9] Update advanced.md

---
 docs/src/examples/advanced.md | 50 +++++++++++++++++++++++++----------
 1 file changed, 36 insertions(+), 14 deletions(-)

diff --git a/docs/src/examples/advanced.md b/docs/src/examples/advanced.md
index 94149b4..56a4f40 100644
--- a/docs/src/examples/advanced.md
+++ b/docs/src/examples/advanced.md
@@ -38,27 +38,49 @@ S = flash_storage(eos, conditions, method = m)
 16
 ```
 
-## GPU and static-array use
+## Immutable performance and GPU use
 
-Use the fully static SSI path for GPU kernels. Construct the isbits EOS and storage
-once outside the kernel, then pass both to it:
+[`flash_2ph_immutable`](@ref) is the public interface to the fully static SSI path.
+It is useful when the component count is small and fixed, particularly inside CPU
+or GPU kernels. Convert the EOS once with [`static_eos`](@ref), and provide the
+overall composition as an `SVector`:
 
 ```julia
-using MultiComponentFlash, StaticArrays
+using BenchmarkTools, MultiComponentFlash, StaticArrays
 
-eos_gpu = static_eos(eos)
-conditions_gpu = (p = p, T = T, z = SVector{length(z)}(z))
-storage = flash_storage(eos_gpu, conditions_gpu; method = SSIFlash(), static = true)
+eos_static = static_eos(eos)
+conditions_static = (p = p, T = T, z = SVector{length(z)}(z))
 
-# Inside each kernel work item:
-K = initial_guess_K(eos_gpu, conditions_gpu, storage)
-V = flash_2ph!(storage, K, eos_gpu, conditions_gpu, NaN)
+V, K = flash_2ph_immutable(eos_static, conditions_static)
+@btime flash_2ph_immutable($eos_static, $conditions_static)
 ```
 
-`static=false` returns the ordinary reusable mutable storage. `static=true` returns
-a `StaticConfig` and keeps all per-flash arrays immutable and local to the work item.
-The static path currently supports `GenericCubicEOS` with `SSIFlash`; compilation is
-specialized on the component count.
+`V` is the scalar vapor fraction and `K` is an `SVector`. All working vectors are
+immutable values local to the call; the input `conditions_static.z` must also be an
+`SVector`. The implementation currently supports `GenericCubicEOS` with
+`SSIFlash`, and compilation is specialized on the number of components.
+
+The two-argument form creates the static storage marker automatically. It can also
+be constructed once and passed as the final positional argument:
+
+```julia
+storage = flash_storage(eos_static, conditions_static; static = true)
+V, K = flash_2ph_immutable(eos_static, conditions_static, storage)
+@btime flash_2ph_immutable($eos_static, $conditions_static, $storage)
+```
+
+Static storage is a zero-size immutable marker rather than a mutable work buffer,
+so constructing it inline normally compiles away and does not allocate. Passing it
+explicitly can still be convenient when setting up a kernel. For example, each
+kernel work item can construct its conditions and call:
+
+```julia
+conditions_i = (p = pressure[i], T = temperature[i], z = z_static)
+V, K = flash_2ph_immutable(eos_static, conditions_i, storage)
+```
+
+Do not share ordinary mutable storage from `flash_storage(...; static = false)`
+between kernel work items.
 
 ## Generate and plot a phase diagram
 

From b928a5f244fff298acb3db1022d37e2080f061e7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Sat, 12 Sep 2026 13:34:25 +0200
Subject: [PATCH 9/9] Renames

---
 docs/src/examples/advanced.md |  4 ++--
 src/MultiComponentFlash.jl    |  4 ++--
 src/static.jl                 |  6 +++---
 test/runtests.jl              | 16 ++++++++--------
 4 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/docs/src/examples/advanced.md b/docs/src/examples/advanced.md
index 56a4f40..a72d4d5 100644
--- a/docs/src/examples/advanced.md
+++ b/docs/src/examples/advanced.md
@@ -42,13 +42,13 @@ S = flash_storage(eos, conditions, method = m)
 
 [`flash_2ph_immutable`](@ref) is the public interface to the fully static SSI path.
 It is useful when the component count is small and fixed, particularly inside CPU
-or GPU kernels. Convert the EOS once with [`static_eos`](@ref), and provide the
+or GPU kernels. Convert the EOS once with [`make_eos_immutable`](@ref), and provide the
 overall composition as an `SVector`:
 
 ```julia
 using BenchmarkTools, MultiComponentFlash, StaticArrays
 
-eos_static = static_eos(eos)
+eos_static = make_eos_immutable(eos)
 conditions_static = (p = p, T = T, z = SVector{length(z)}(z))
 
 V, K = flash_2ph_immutable(eos_static, conditions_static)
diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl
index 15d02fc..0e3b251 100644
--- a/src/MultiComponentFlash.jl
+++ b/src/MultiComponentFlash.jl
@@ -22,7 +22,7 @@ module MultiComponentFlash
     export flash_2ph, flash_2ph!, flash_2ph_immutable, flash_storage
     export stability_2ph, stability_2ph!
     # Algorithms for flash
-    export SSIFlash, NewtonFlash, SSINewtonFlash, FlashConfig, StaticConfig
+    export SSIFlash, NewtonFlash, SSINewtonFlash
     # Mixtures and their molecular makeup
     export MolecularProperty, MultiComponentMixture
     # K-values
@@ -35,7 +35,7 @@ module MultiComponentFlash
 
     export force_scalars, force_coefficients, force_coefficients!
     export critical_pressure, critical_temperature, critical_volume, acentric_factor, molar_weight
-    export static_eos
+    export make_eos_immutable
     export cubic_benchmark
     export single_phase_label
 
diff --git a/src/static.jl b/src/static.jl
index ae8ceb6..7be38de 100644
--- a/src/static.jl
+++ b/src/static.jl
@@ -23,7 +23,7 @@ scalar vapor fraction. When `storage` is omitted, a static storage marker is
 created automatically.
 
 The immutable path currently supports `SSIFlash` and generic cubic EOS values
-converted with [`static_eos`](@ref).
+converted with [`make_eos_immutable`](@ref).
 """
 @inline function flash_2ph_immutable(eos, c; method = SSIFlash(), kwarg...)
     return flash_2ph_immutable(eos, c,
@@ -51,11 +51,11 @@ function static_mixture(mixture::MultiComponentMixture{R, N}) where {R, N}
 end
 
 """
-    static_eos(eos)
+    make_eos_immutable(eos)
 
 Convert a generic cubic EOS to an isbits representation for accelerator kernels.
 """
-function static_eos(eos::GenericCubicEOS{T, R, N}) where {T, R, N}
+function make_eos_immutable(eos::GenericCubicEOS{T, R, N}) where {T, R, N}
     mixture = static_mixture(eos.mixture)
     volume_shift = eos.volume_shift
     if !isnothing(volume_shift)
diff --git a/test/runtests.jl b/test/runtests.jl
index 6519621..1e17c26 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -42,7 +42,7 @@ end
 end
 
 @testset "Static accelerator path" begin
-    eos = static_eos(get_test_eos())
+    eos = make_eos_immutable(get_test_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)
@@ -50,20 +50,20 @@ end
 
     @test isbitstype(typeof(eos))
     @test isbitstype(typeof(storage))
-    @test storage isa StaticConfig
+    @test storage isa MultiComponentFlash.StaticConfig
     @test report.converged
     @test report.its == 6
     @test V ≈ 0.7632068334421974
     @test K ≈ @SVector [4.553402802323027, 17.73895830809456, 0.0004031451448211194]
 
-    normal_config = FlashConfig(print_output=false)
-    @test typeof(normal_config) == FlashConfig
+    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)
 
-    config = FlashConfig(print_output=false, use_dict_storage=false)
-    @test typeof(config) == FlashConfig
-    @test flash_storage(eos, c, SSIFlash(), config) isa StaticConfig
+    config = MultiComponentFlash.FlashConfig(print_output=false, use_dict_storage=false)
+    @test typeof(config) == MultiComponentFlash.FlashConfig
+    @test flash_storage(eos, c, SSIFlash(), config) isa MultiComponentFlash.StaticConfig
     K_config = initial_guess_K(eos, c, storage)
     V_config, K_config, config_report = flash_2ph(eos, c, K_config, NaN, config;
         method=SSIFlash(), extra_out=true, z_min=nothing)
@@ -162,7 +162,7 @@ using StaticArrays, KernelAbstractions, JLArrays
     end
     if isdefined(JLArrays, :JLBackend)
         host_eos = get_test_eos()
-        eos = static_eos(host_eos)
+        eos = make_eos_immutable(host_eos)
         z = @SVector [0.5, 0.3, 0.2]
         storage = flash_storage(eos, (p = 1e5, T = 300.0, z = z);
             method = SSIFlash(), static = true)