Make sure that phases are passed onto static flash and add static path for Søreide-Whitson

src/eos.jl

@@ -52,7 +52,8 @@ Base.@propagate_inbounds function binary_interaction(eos::AbstractEOS, i::Int, j
 end
 
 Base.@propagate_inbounds function binary_interaction(mixture::MultiComponentMixture{R}, i, j) where {R}
-    return binary_interaction(mixture.binary_interaction, i, j)::R
+    bic = mixture.binary_interaction
+    return isnothing(bic) ? zero(R) : binary_interaction(bic, i, j)::R
 end
 
 function binary_interaction(::Nothing, i, j)

src/eos_types.jl

@@ -95,14 +95,21 @@ Note that this EOS requires water to be present in the mixture by name (it is
 after all the whole point of the Søreide-Whitson EOS). Special treatment is
 used for H2S, CO2, and N2 when given by name.
 """
-struct SoreideWhitson{T} <: AbstractPengRobinson
+struct SoreideWhitson{T, C<:AbstractVector{COMPONENT_ENUM}} <: AbstractPengRobinson
     A::NTuple{3, T}
     A_mw::NTuple{3, T}
     alphas::NTuple{3, T}
     water_coefficients::NTuple{3, T}
     molality::T
     T_co2::T
-    component_types::Vector{COMPONENT_ENUM}
+    component_types::C
+end
+
+# Retain the explicitly parameterized constructor used by existing callers.
+function SoreideWhitson{T}(A, A_mw, alphas, water_coefficients,
+        molality, T_co2, component_types::C) where {T, C<:AbstractVector{COMPONENT_ENUM}}
+    return SoreideWhitson{T, C}(A, A_mw, alphas, water_coefficients,
+        molality, T_co2, component_types)
 end
 
 function SoreideWhitson(mixture_or_cnames::Union{MultiComponentMixture, Vector{String}};

src/flow_coupler.jl

@@ -31,22 +31,24 @@ function flashed_mixture_2ph!(storage, eos, conditions, K; kwarg...)
     z = conditions.z
     forces = storage.forces
     if isnan(V)
-        V = single_phase_label(eos.mixture, conditions)
+        V = single_phase_label(eos, conditions)
         if V == 0
             state = single_phase_l
+            phase = :liquid
         else
             state = single_phase_v
+            phase = :vapor
         end
-        liquid = vapor = cond
-        Z_L = Z_V = mixture_compressibility_factor(eos, conditions, forces)
+        phase_cond = set_phase(conditions, phase)
+        Z_L = Z_V = mixture_compressibility_factor(eos, phase_cond, forces)
         x = copy(z)
         y = copy(z)
     else
         x = @. liquid_mole_fraction(z, K, V)
         y = @. vapor_mole_fraction(x, K)
         state = two_phase_lv
-        liquid = (p = p, T = T, z = x)
-        vapor = (p = p, T = T, z = y)
+        liquid = set_phase((p = p, T = T, z = x), :liquid)
+        vapor = set_phase((p = p, T = T, z = y), :vapor)
         Z_L = mixture_compressibility_factor(eos, liquid, forces)
         Z_V = mixture_compressibility_factor(eos, vapor, forces)
     end

src/implicit_flash_derivatives.jl

@@ -145,7 +145,11 @@ end
         y = SVector{$N, F}(($(y_entries...),))
         liquid = (p = cond.p, T = cond.T, z = x, phase = Val(:liquid))
         vapor = (p = cond.p, T = cond.T, z = y, phase = Val(:vapor))
-        forces = static_force_coefficients(eos, cond, F)
+        # Coefficients depend on pressure and temperature, not the seeded
+        # K-values or vapor fraction. Keep them numeric in the unknown Jacobian
+        # instead of filling both phase matrices with zero-derivative duals.
+        force_type = Base.promote_eltype(cond.p, cond.T)
+        forces = static_force_coefficients(eos, cond, force_type)
         Z_l, scalars_l = prep(eos, liquid, forces)
         Z_v, scalars_v = prep(eos, vapor, forces)
         lnphi_l = static_fugacity_coefficients(

src/static.jl

@@ -107,7 +107,9 @@ end
 Convert a generic cubic EOS to an isbits representation for accelerator kernels.
 """
 function make_eos_immutable(eos::GenericCubicEOS{T, R, N}; float_type = missing) where {T, R, N}
-    mixture = static_mixture(eos.mixture)
+    mixture = static_mixture(eos.mixture; float_type = float_type)
+    type = static_cubic_type(eos.type; float_type = float_type)
+    C = ismissing(float_type) ? identity : float_type
     volume_shift = eos.volume_shift
     if !isnothing(volume_shift)
         if ismissing(float_type)
@@ -118,16 +120,29 @@ function make_eos_immutable(eos::GenericCubicEOS{T, R, N}; float_type = missing)
         volume_shift = SVector{N, v_type}(volume_shift)
     end
     return GenericCubicEOS(
-        eos.type,
+        type,
         mixture,
-        eos.m_1,
-        eos.m_2,
-        eos.ω_a,
-        eos.ω_b,
+        C(eos.m_1),
+        C(eos.m_2),
+        C(eos.ω_a),
+        C(eos.ω_b),
         volume_shift
     )
 end
 
+# Most cubic types carry no mutable data. Søreide-Whitson additionally stores
+# component classifications, which must be embedded in the kernel argument.
+static_cubic_type(type; float_type = missing) = type
+
+function static_cubic_type(sw::SoreideWhitson{T}; float_type = missing) where T
+    R = ismissing(float_type) ? T : float_type
+    N = length(sw.component_types)
+    return SoreideWhitson{R}(
+        map(R, sw.A), map(R, sw.A_mw), map(R, sw.alphas),
+        map(R, sw.water_coefficients), R(sw.molality), R(sw.T_co2),
+        SVector{N, COMPONENT_ENUM}(sw.component_types))
+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}
@@ -139,6 +154,10 @@ end
 # Val phase tags keep Symbol construction and dynamic dispatch out of kernels.
 @inline phase_symbol(::Val{phase}) where phase = phase
 
+@inline function set_phase(cond, phase::Val)
+    return (p = cond.p, T = cond.T, z = cond.z, phase = phase)
+end
+
 @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
@@ -150,7 +169,17 @@ end
 end
 
 """Immutable force coefficients for accelerator kernels."""
-@inline function static_force_coefficients(eos::GenericCubicEOS{E, R, N}, cond,
+@inline function static_force_coefficients(eos::GenericCubicEOS, cond, ::Type{T}) where T
+    if forces_per_phase(eos)
+        liquid = static_phase_force_coefficients(eos, set_phase(cond, Val(:liquid)), T)
+        vapor = static_phase_force_coefficients(eos, set_phase(cond, Val(:vapor)), T)
+        return (liquid = liquid, vapor = vapor)
+    else
+        return static_phase_force_coefficients(eos, cond, T)
+    end
+end
+
+@inline function static_phase_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)))
@@ -520,7 +549,7 @@ end
     K_out = report.stable ? K_liquid : static_divide(y, x)
     # Only a complete, converged two-sided test can establish a new reference.
     # Skipped trial phases are reported as stable above, but prove no such thing.
-    if update_bypass && check_liquid && check_vapor && report.stable &&
+    if update_bypass && !forces_per_phase(eos) && check_liquid && check_vapor && report.stable &&
             report.liquid.trivial && report.vapor.trivial
         critical_distance = static_michelsen_critical_point_measure(
             eos, cond.p, cond.T, cond.z)

test/eos_equations.jl

@@ -120,6 +120,66 @@ end
     @test report.converged
     @test isfinite(vapor_fraction)
     @test all(isfinite, K)
+
+    # Both simulator flash paths must retain the phase-specific interactions.
+    static_cond = (p = cond.p, T = cond.T, z = SVector{2}(cond.z))
+    static_eos = make_eos_immutable(eos)
+    @test isbitstype(typeof(static_eos))
+    @test static_eos.type.component_types == eos.type.component_types
+    @test static_eos.type.molality == molality
+    @test static_eos.type.water_coefficients == eos.type.water_coefficients
+    @test make_eos_immutable(static_eos) == static_eos
+    static_forces = MCF.static_force_coefficients(static_eos, static_cond, Float64)
+    for phase in (:liquid, :vapor)
+        phase_cond = MCF.set_phase(static_cond, phase)
+        host_forces = MCF.force_coefficients(eos, phase_cond)
+        immutable_forces = MCF.get_force_coefficients(static_forces, static_eos, phase_cond)
+        @test immutable_forces.A_ij ≈ host_forces.A_ij
+        @test immutable_forces.B_i ≈ host_forces.B_i
+    end
+    V_static, K_static = MCF.flash_2ph_immutable(static_eos, static_cond)
+    @test V_static ≈ vapor_fraction rtol = 1e-7
+    @test K_static ≈ K rtol = 1e-7
+    mixture_flash = MCF.flashed_mixture_2ph(eos, cond)
+    @test mixture_flash.V ≈ vapor_fraction
+    @test isfinite(mixture_flash.liquid.Z)
+    @test isfinite(mixture_flash.vapor.Z)
+
+    ad_cond = (p = MCF.ForwardDiff.Dual{Nothing}(cond.p, 1.0),
+        T = cond.T, z = static_cond.z)
+    V_ad, _ = MCF.implicit_flash_derivatives(static_eos, static_cond, ad_cond, V_static, K_static)
+    h = 100.0
+    V_plus, _ = MCF.flash_2ph_immutable(static_eos, (; static_cond..., p = cond.p + h))
+    V_minus, _ = MCF.flash_2ph_immutable(static_eos, (; static_cond..., p = cond.p - h))
+    @test MCF.ForwardDiff.partials(V_ad)[1] ≈ (V_plus - V_minus)/(2h) rtol = 1e-4
+end
+
+@testset "Immutable Søreide–Whitson component classifications" begin
+    mixture = MultiComponentMixture(["Water", "CarbonDioxide", "Nitrogen",
+        "HydrogenSulfide", "Methane", "Hydrogen"])
+    sw = SoreideWhitson(mixture; molality = 1.5, T_co2 = 301.0,
+        A = (1.2, 1.3, -0.2), water_coefficients = (0.45, 0.98, 0.004))
+    eos = GenericCubicEOS(mixture, sw; volume_shift = fill(0.1, 6))
+    for float_type in (Float64, Float32)
+        immutable = make_eos_immutable(eos; float_type = float_type)
+        @test isbitstype(typeof(immutable))
+        @test eltype(immutable.type.A) == float_type
+        @test typeof(immutable.type.molality) == float_type
+        @test eltype(immutable.volume_shift) == float_type
+        @test typeof(immutable.mixture.properties[1].mw) == float_type
+        @test immutable.type.component_types == sw.component_types
+        @test immutable.type.T_co2 == float_type(sw.T_co2)
+        @test immutable.type.A == map(float_type, sw.A)
+        @test immutable.type.A_mw == map(float_type, sw.A_mw)
+        @test immutable.type.alphas == map(float_type, sw.alphas)
+        @test immutable.type.water_coefficients == map(float_type, sw.water_coefficients)
+        cond = (p = 1e6, T = 350.0, z = SVector{6}(fill(1/6, 6)))
+        for phase in (:liquid, :vapor), i in 1:6, j in 1:6
+            phase_cond = MCF.set_phase(cond, phase)
+            @test MCF.binary_interaction(immutable, i, j, phase_cond) ≈
+                MCF.binary_interaction(eos, i, j, phase_cond) rtol = 1e-6
+        end
+    end
 end
 
 @testset "Cubic root degeneracies" begin