diff --git a/src/eos.jl b/src/eos.jl index 78b9e71..1456a7e 100644 --- a/src/eos.jl +++ b/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) diff --git a/src/eos_types.jl b/src/eos_types.jl index 84b095c..fc4500f 100644 --- a/src/eos_types.jl +++ b/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}}; diff --git a/src/flow_coupler.jl b/src/flow_coupler.jl index e31affb..11efb28 100644 --- a/src/flow_coupler.jl +++ b/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 diff --git a/src/implicit_flash_derivatives.jl b/src/implicit_flash_derivatives.jl index 3f29cc9..e8fcba2 100644 --- a/src/implicit_flash_derivatives.jl +++ b/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( diff --git a/src/static.jl b/src/static.jl index 2472cbf..dbc3181 100644 --- a/src/static.jl +++ b/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) diff --git a/test/eos_equations.jl b/test/eos_equations.jl index e2e2f47..b258238 100644 --- a/test/eos_equations.jl +++ b/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