Add a static path for implicit function theorem derivatives of flash

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 <olav.moyner@gmail.com>"]
 
 [deps]

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"]

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")

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`.
+"""
+@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

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

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()