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/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..a72d4d5 100644
--- a/docs/src/examples/advanced.md
+++ b/docs/src/examples/advanced.md
@@ -38,64 +38,49 @@ S = flash_storage(eos, conditions, method = m)
 16
 ```
 
-See the unit tests for examples where the flash can use `StaticArrays` to avoid allocations entirely.
+## Immutable performance and GPU use
 
-## Performance example
+[`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 [`make_eos_immutable`](@ref), and provide the
+overall composition as an `SVector`:
 
-The default interface is designed for ease-of-use with standard Julia types, but the module  also supports further by using `StaticArrays`:
+```julia
+using BenchmarkTools, MultiComponentFlash, StaticArrays
+
+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)
+@btime flash_2ph_immutable($eos_static, $conditions_static)
+```
+
+`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
-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)
+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)
 ```
 
-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.
+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
-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)
+conditions_i = (p = pressure[i], T = temperature[i], z = z_static)
+V, K = flash_2ph_immutable(eos_static, conditions_i, storage)
 ```
 
-!!! 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.
+Do not share ordinary mutable storage from `flash_storage(...; static = false)`
+between kernel work items.
 
 ## Generate and plot a phase diagram
 
diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl
index 6514b7e..0e3b251 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
@@ -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 make_eos_immutable
     export cubic_benchmark
     export single_phase_label
 
@@ -53,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 c4ceeed..d8264ef 100644
--- a/src/eos.jl
+++ b/src/eos.jl
@@ -11,9 +11,11 @@ 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
+
 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 +95,40 @@ 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
+
+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 +146,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)
diff --git a/src/eos_types.jl b/src/eos_types.jl
index 748ada2..55f57ef 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
diff --git a/src/flash.jl b/src/flash.jl
index 02fabd6..412adce 100644
--- a/src/flash.jl
+++ b/src/flash.jl
@@ -33,16 +33,25 @@ 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
+    if !use_dict_storage(config)
+        return flash_2ph(eos, c, K, V, StaticConfig(); method = method, kwarg...)
+    end
     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 +78,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 +105,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 +114,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
@@ -117,10 +132,10 @@ function flash_2ph_impl!(storage, K, eos, c, V = NaN;
             V, ϵ = flash_update!(K, storage, method, eos, c, forces, V, i)
             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
@@ -146,13 +161,25 @@ 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...)
+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; 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
@@ -162,11 +189,10 @@ 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)
@@ -176,42 +202,31 @@ function flash_storage_internal!(out, eos, cond, method;
     else
         out[:forces] = alloc_forces(cond)
     end
-    if static_size
-        alloc_vec = () -> @MVector zeros(n)
-    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 = @MVector zeros(T, n)
-        y_ad = @MVector zeros(T, n)
-        r = @MVector zeros(np)
-        J = @MMatrix zeros(np, np)
-    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
 
@@ -221,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
@@ -234,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 = @MVector zeros(T_cond, n)
-        J_inv = @MMatrix zeros(np, external_partials)
-    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 = @MVector zeros(npartials)
-        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}
diff --git a/src/flash_types.jl b/src/flash_types.jl
index eba83a5..23fd6e5 100644
--- a/src/flash_types.jl
+++ b/src/flash_types.jl
@@ -3,6 +3,21 @@ abstract type AbstractFlash end
 "Abstract type for all flash types that use Newton in some form"
 abstract type AbstractNewtonFlash <: AbstractFlash end
 
+"""
+    FlashConfig(; print_output=true, use_dict_storage=true)
+
+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)`.
+"""
+Base.@kwdef struct FlashConfig
+    print_output::Bool = true
+    use_dict_storage::Bool = true
+end
+
+@inline print_output(config::FlashConfig) = config.print_output
+@inline use_dict_storage(config::FlashConfig) = config.use_dict_storage
+
 """
     ssi = SSIFlash()
 
@@ -46,11 +61,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"
@@ -66,16 +84,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 0c4109a..8d164cd 100644
--- a/src/kvalues.jl
+++ b/src/kvalues.jl
@@ -60,4 +60,4 @@ 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)
diff --git a/src/kvalues_eos.jl b/src/kvalues_eos.jl
index f19d3ec..62307f3 100644
--- a/src/kvalues_eos.jl
+++ b/src/kvalues_eos.jl
@@ -18,6 +18,14 @@ function flash_storage(eos::KValuesEOS, cond = missing; kwarg...)
     return nothing
 end
 
+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)
 end
+
+
+function flash_2ph!(storage, K, eos::KValuesEOS, cond, V,
+        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 b019226..ad1a667 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,7 +77,7 @@ 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
 
diff --git a/src/stability.jl b/src/stability.jl
index 92cb3f5..3433a22 100644
--- a/src/stability.jl
+++ b/src/stability.jl
@@ -8,11 +8,22 @@ 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...)
+    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
+
+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,
@@ -31,7 +42,7 @@ function stability_2ph!(storage, K, eos, c;
     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 +55,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 +70,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 +91,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 +139,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/src/static.jl b/src/static.jl
new file mode 100644
index 0000000..7be38de
--- /dev/null
+++ b/src/static.jl
@@ -0,0 +1,383 @@
+"""
+    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
+
+"""
+    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 [`make_eos_immutable`](@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))
+    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
+
+"""
+    make_eos_immutable(eos)
+
+Convert a generic cubic EOS to an isbits representation for accelerator kernels.
+"""
+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)
+        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/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..1e17c26 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -38,15 +38,47 @@ 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 = 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)
+    V, K, report = flash_2ph!(storage, K, eos, c, 0.5; extra_out = true)
+
+    @test isbitstype(typeof(eos))
+    @test isbitstype(typeof(storage))
+    @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 = 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 = 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)
+    @test config_report.stability isa MultiComponentFlash.StabilityReport
+    @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
@@ -115,3 +147,44 @@ 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, storage)
+        i = @index(Global)
+        if i <= length(out)
+            @inbounds cond = (p = pressure[i], T = temperature[i], z = z)
+            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
+    end
+    if isdefined(JLArrays, :JLBackend)
+        host_eos = get_test_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)
+        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, storage; 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
+
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