From 1743a761275b75b62b5ff07dbe26f90fc43c19ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 3 Oct 2026 21:11:22 +0200 Subject: [PATCH 1/3] Add a static path for implicit function theorem derivatives of flash --- Project.toml | 2 +- src/MultiComponentFlash.jl | 3 + src/implicit_flash_derivatives.jl | 159 ++++++++++++++++++++++++++++++ test/implicit_derivatives.jl | 87 ++++++++++++++++ test/runtests.jl | 1 + 5 files changed, 251 insertions(+), 1 deletion(-) create mode 100644 src/implicit_flash_derivatives.jl create mode 100644 test/implicit_derivatives.jl diff --git a/Project.toml b/Project.toml index 038430d..4c9977f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MultiComponentFlash" uuid = "35e5bd01-9722-4017-9deb-64a5d32478ff" -version = "2.0.0" +version = "2.0.1" authors = ["Olav Møyner "] [deps] diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl index 8794ec2..10a5e6d 100644 --- a/src/MultiComponentFlash.jl +++ b/src/MultiComponentFlash.jl @@ -44,6 +44,7 @@ module MultiComponentFlash export SinglePhaseLiquid, SinglePhaseVapor, TwoPhaseLiquidVapor, FlashedPhase, UnknownPhaseState, FlashedMixture2Phase export phase_saturations, mass_density, mass_densities, molar_volume, phase_is_present export inverse_flash_update! + export implicit_flash_derivatives export lbc_viscosity, lbc_viscosities export set_partials, set_partials_phase_mole_fractions!, set_partials_vapor_fraction @@ -56,6 +57,8 @@ module MultiComponentFlash include("derivatives.jl") include("stability.jl") include("static.jl") + include("implicit_flash_derivatives.jl") + using .ImplicitFlashDerivatives: implicit_flash_derivatives include("tables.jl") include("flow_coupler.jl") diff --git a/src/implicit_flash_derivatives.jl b/src/implicit_flash_derivatives.jl new file mode 100644 index 0000000..e17359e --- /dev/null +++ b/src/implicit_flash_derivatives.jl @@ -0,0 +1,159 @@ +module ImplicitFlashDerivatives + +using LinearAlgebra +using ForwardDiff +using StaticArrays: SVector, SMatrix +using ..MultiComponentFlash: GenericCubicEOS, KValuesEOS, + liquid_mole_fraction, vapor_mole_fraction, objectiveRR, objectiveRR_dV, + static_force_coefficients, static_fugacity_coefficients, prep + +export implicit_flash_derivatives + +""" + implicit_flash_derivatives(eos::GenericCubicEOS, numeric_cond, ad_cond, + V_numeric, K_numeric) + +Differentiate an already converged, interior two-phase cubic flash. Return +`(V_ad, K_ad)`, with derivatives of the vapor fraction and equilibrium K-values +with respect to the AD inputs in `ad_cond`. + +`numeric_cond` and `ad_cond` each contain `p`, `T`, and `z`. Their primal values +must agree: `numeric_cond` is the numeric condition used to solve the flash, +while any of `ad_cond.p`, `ad_cond.T`, and `ad_cond.z` may carry AD derivatives. +Both `z` fields and `K_numeric` must be `SVector`s of the EOS component count. +`V_numeric` and `K_numeric` are the converged numeric vapor fraction and +equilibrium K-values at `numeric_cond`; this function does not solve a flash. + +The function differentiates the fugacity equilibrium and Rachford-Rice +residuals, then solves their static `(N+1) × (N+1)` Jacobian system. It requires +`0 < V_numeric < 1`, positive compositions and K-values, and a nonsingular +Jacobian. It does not differentiate stability decisions or phase boundaries. +For accelerator use, pass an immutable EOS from [`make_eos_immutable`](@ref). +""" +@inline function implicit_flash_derivatives( + eos::GenericCubicEOS{E, R, N}, numeric_cond, ad_cond, + V_numeric, K_numeric::SVector{N}) where {E, R, N} + u = static_flash_unknowns(K_numeric, V_numeric) + + # Seed only the equilibrium unknowns for the primary Jacobian. The + # condition here is numeric, so no nested AD types are needed. + u_seeded = seed_static_flash_unknowns(u) + residual_seeded = static_equilibrium_residual(eos, numeric_cond, u_seeded) + jacobian = static_flash_jacobian(residual_seeded) + residual_numeric = static_flash_residual_values(residual_seeded) + + # The constant residual is subtracted so the returned primals remain the + # converged numeric solution, even when its tolerance is finite. + T = promote_type(typeof(ad_cond.p), typeof(ad_cond.T), eltype(ad_cond.z)) + ad_cond = (p = convert(T, ad_cond.p), T = convert(T, ad_cond.T), + z = SVector{N, T}(ad_cond.z)) + u_ad = SVector{N+1, T}(u) + residual_ad = static_equilibrium_residual(eos, ad_cond, u_ad) + correction = jacobian \ (residual_ad - residual_numeric) + updated = SVector{N+1, T}(u_ad - correction) + return updated[N+1], static_flash_K(updated) +end + +""" + implicit_flash_derivatives(eos::KValuesEOS, numeric_cond, ad_cond, + V_numeric, K_numeric, K_ad) + +Differentiate an already converged, interior two-phase K-value flash. Return +`(V_ad, K_ad)` using the implicit derivative of the Rachford-Rice equation. + +`numeric_cond` and `ad_cond` each contain `p`, `T`, and `z` with matching primal +values. Both `z` fields, `K_numeric`, and `K_ad` must be `SVector`s of the EOS +component count. `V_numeric` is the converged numeric vapor fraction for +`numeric_cond` and `K_numeric`. Evaluate `K_ad` at `ad_cond` before calling this +function; its primal values must match `K_numeric`. Pressure and temperature +derivatives enter through `K_ad`, while composition derivatives enter through +`ad_cond.z`. + +This function does not solve a flash or evaluate K-values. It requires +`0 < V_numeric < 1`, positive K-values, and a nonzero Rachford-Rice slope. It +does not differentiate phase boundaries. +""" +@inline function implicit_flash_derivatives( + eos::KValuesEOS{E, R, N}, numeric_cond, ad_cond, + V_numeric, K_numeric::SVector{N}, K_ad::SVector{N}) where {E, R, N} + residual_numeric = objectiveRR(V_numeric, K_numeric, numeric_cond.z) + T = promote_type(typeof(ad_cond.p), typeof(ad_cond.T), + eltype(ad_cond.z), eltype(K_ad)) + K = SVector{N, T}(K_ad) + z = SVector{N, T}(ad_cond.z) + residual_ad = objectiveRR(V_numeric, K, z) + slope = objectiveRR_dV(V_numeric, K_numeric, numeric_cond.z) + V = convert(T, V_numeric) - convert(T, + (residual_ad - residual_numeric)/slope) + return V, K +end + +@generated function static_flash_unknowns(K::SVector{N, F}, V) where {N, F} + entries = Any[:(K[$i]) for i in 1:N] + push!(entries, :(V)) + return :(SVector{$(N+1)}(($(entries...),))) +end + +@generated function seed_static_flash_unknowns( + u::SVector{M, F}) where {M, F} + entries = [:(D(u[$i], ForwardDiff.single_seed( + ForwardDiff.Partials{$M, F}, Val($i)))) for i in 1:M] + return quote + D = ForwardDiff.Dual{Nothing, F, $M} + SVector{$M, D}(($(entries...),)) + end +end + +@generated function static_flash_jacobian( + residual::SVector{M, <:ForwardDiff.Dual}) where M + entries = [:(residual[$row].partials[$col]) + for col in 1:M for row in 1:M] + return :(SMatrix{$M, $M}(($(entries...),))) +end + +@generated function static_flash_residual_values( + residual::SVector{M, <:ForwardDiff.Dual}) where M + entries = [:(ForwardDiff.value(residual[$i])) for i in 1:M] + return :(SVector{$M}(($(entries...),))) +end + +@generated function static_flash_K(u::SVector{M}) where M + entries = [:(u[$i]) for i in 1:(M-1)] + return :(SVector{$(M-1)}(($(entries...),))) +end + +@generated function static_equilibrium_residual( + eos::GenericCubicEOS{E, R, N}, cond, + unknowns::SVector{M, F}) where {E, R, N, M, F} + M == N + 1 || error("Expected N K-values and one vapor fraction") + k_entries = [:(unknowns[$i]) for i in 1:N] + x_entries = [:(liquid_mole_fraction(z[$i], K[$i], V)) for i in 1:N] + y_entries = [:(vapor_mole_fraction(x[$i], K[$i])) for i in 1:N] + residual_entries = [:(log(K[$i]) + lnphi_v[$i] - lnphi_l[$i]) + for i in 1:N] + rr = :(zero(F)) + for i in 1:N + rr = :($rr + z[$i]*(K[$i] - one(F))/ + (one(F) + V*(K[$i] - one(F)))) + end + push!(residual_entries, rr) + return quote + V = unknowns[$M] + K = SVector{$N, F}(($(k_entries...),)) + z = cond.z + x = SVector{$N, F}(($(x_entries...),)) + 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) + Z_l, scalars_l = prep(eos, liquid, forces) + Z_v, scalars_v = prep(eos, vapor, forces) + lnphi_l = static_fugacity_coefficients( + eos, liquid, Z_l, forces, scalars_l, F) + lnphi_v = static_fugacity_coefficients( + eos, vapor, Z_v, forces, scalars_v, F) + SVector{$M, F}(($(residual_entries...),)) + end +end + +end # module ImplicitFlashDerivatives diff --git a/test/implicit_derivatives.jl b/test/implicit_derivatives.jl new file mode 100644 index 0000000..690289c --- /dev/null +++ b/test/implicit_derivatives.jl @@ -0,0 +1,87 @@ +using ForwardDiff + +primal(x) = x +primal(x::ForwardDiff.Dual) = ForwardDiff.value(x) + +function central_difference(f, u, i, h) + plus = copy(u) + minus = copy(u) + plus[i] += h + minus[i] -= h + return (f(plus) - f(minus))/(2h) +end + +@testset "Implicit static flash derivatives" begin + cubic = make_eos_immutable(get_test_eos()) + function cubic_result(u) + z = SVector(u[3], u[4], 1-u[3]-u[4]) + ad_cond = (p = 1e6*u[1], T = u[2], z = z) + numeric_cond = (p = primal(ad_cond.p), T = primal(ad_cond.T), + z = map(primal, z)) + V, K = flash_2ph_immutable(cubic, numeric_cond; check = false) + return implicit_flash_derivatives(cubic, numeric_cond, ad_cond, V, K) + end + u = [1.0, 300.0, 0.6, 0.1] + @test 0 < cubic_result(u)[1] < 1 + for f in (v -> cubic_result(v)[1], + v -> cubic_result(v)[2][1]) + for (i, h) in ((2, 1e-2), (3, 1e-5)) + derivative = ForwardDiff.derivative(t -> begin + v = [j == i ? t : u[j] for j in eachindex(u)] + f(v) + end, u[i]) + @test derivative ≈ central_difference(f, u, i, h) rtol=1e-4 + end + end + + mixture = MultiComponentMixture(["CarbonDioxide", "Water"]) + kvalue = KValuesEOS(cond -> SVector( + 0.05*(cond.p/1e6)*(cond.T/300)^0.2, + 5.0*(1e6/cond.p)*(cond.T/300)^(-0.1)), mixture) + function kvalue_result(u) + z = SVector(u[3], 1-u[3]) + ad_cond = (p = 1e6*u[1], T = u[2], z = z) + K_ad = SVector(initial_guess_K(kvalue, ad_cond)) + K_numeric = map(primal, K_ad) + numeric_cond = (p = primal(ad_cond.p), T = primal(ad_cond.T), + z = map(primal, z)) + V = solve_rachford_rice(K_numeric, numeric_cond.z) + return implicit_flash_derivatives(kvalue, numeric_cond, ad_cond, + V, K_numeric, K_ad) + end + u_k = [1.0, 300.0, 0.4] + @test 0 < kvalue_result(u_k)[1] < 1 + gradient = ForwardDiff.gradient(v -> kvalue_result(v)[1], u_k) + for i in eachindex(u_k) + f = v -> kvalue_result(v)[1] + @test gradient[i] ≈ central_difference(f, u_k, i, 1e-4) rtol=1e-4 + end + + # A constant K evaluator can leave K_ad numeric even when the condition + # carries AD seeds. The returned vapor fraction and K-values must still + # use the condition's promoted type. + constant_K = @SVector [0.05, 5.0] + numeric_cond = (p = 1e6, T = 300.0, z = @SVector [0.4, 0.6]) + V_numeric = solve_rachford_rice(constant_K, numeric_cond.z) + function constant_k_result(; p = numeric_cond.p, T = numeric_cond.T, + z = numeric_cond.z) + ad_cond = (; p, T, z) + return implicit_flash_derivatives(kvalue, numeric_cond, ad_cond, + V_numeric, constant_K, constant_K) + end + for variable in (:p, :T) + value = getproperty(numeric_cond, variable) + sensitivity = ForwardDiff.derivative(value) do seeded + cond = merge(numeric_cond, NamedTuple{(variable,)}((seeded,))) + constant_k_result(; cond...)[1] + end + @test sensitivity == 0 + end + composition_sensitivity = ForwardDiff.derivative(0.4) do z1 + constant_k_result(z = SVector(z1, 1-z1))[1] + end + composition_finite_difference = (solve_rachford_rice(constant_K, + SVector(0.4+1e-5, 0.6-1e-5)) - solve_rachford_rice(constant_K, + SVector(0.4-1e-5, 0.6+1e-5)))/(2e-5) + @test composition_sensitivity ≈ composition_finite_difference rtol=1e-6 +end diff --git a/test/runtests.jl b/test/runtests.jl index f3ff01d..7ab2dcd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,6 @@ include("test_setup.jl") include("eos_equations.jl") +include("implicit_derivatives.jl") @testset "Rachford-Rice" begin test_rachford_rice() From ea9f0c37c454467ab1fab5049a37b4d20400e82f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 4 Oct 2026 10:03:55 +0200 Subject: [PATCH 2/3] Update utilities.md --- docs/src/api/utilities.md | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/docs/src/api/utilities.md b/docs/src/api/utilities.md index 62c80fa..34efaec 100644 --- a/docs/src/api/utilities.md +++ b/docs/src/api/utilities.md @@ -1,10 +1,12 @@ # Utilities + ```@index Pages = ["utilities.md"] ``` ## Partial derivatives -!!! note "Experimental features" + +!!! note "Experimental features" Functions for obtaining partial derivatives of the flashed results. Please note that this is an experimental feature. Examples of usage are found in the unit tests. ```@autodocs @@ -15,7 +17,8 @@ Private = false ``` ## Coupling utilities -!!! note "Experimental features" + +!!! note "Experimental features" Utilities for coupling flash to simulation codes. Please note that this is an experimental feature. ```@autodocs @@ -25,8 +28,12 @@ Order = [:type, :function] Private = false ``` +```@docs +implicit_flash_derivatives +``` ## Various + ```@autodocs Modules = [MultiComponentFlash] Pages = ["utils.jl"] From 8fe2506b17db7e0e0307f2883067712a50a46ac2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 4 Oct 2026 10:25:01 +0200 Subject: [PATCH 3/3] Update implicit_flash_derivatives.jl --- src/implicit_flash_derivatives.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/implicit_flash_derivatives.jl b/src/implicit_flash_derivatives.jl index e17359e..3f29cc9 100644 --- a/src/implicit_flash_derivatives.jl +++ b/src/implicit_flash_derivatives.jl @@ -27,8 +27,8 @@ equilibrium K-values at `numeric_cond`; this function does not solve a flash. The function differentiates the fugacity equilibrium and Rachford-Rice residuals, then solves their static `(N+1) × (N+1)` Jacobian system. It requires `0 < V_numeric < 1`, positive compositions and K-values, and a nonsingular -Jacobian. It does not differentiate stability decisions or phase boundaries. -For accelerator use, pass an immutable EOS from [`make_eos_immutable`](@ref). +Jacobian. It does not differentiate stability decisions or phase boundaries. For +accelerator use, pass an immutable EOS from `make_eos_immutable`. """ @inline function implicit_flash_derivatives( eos::GenericCubicEOS{E, R, N}, numeric_cond, ad_cond,