Add safeguarded Newton flashes for the static path

Project.toml

@@ -1,6 +1,6 @@
 name = "MultiComponentFlash"
 uuid = "35e5bd01-9722-4017-9deb-64a5d32478ff"
-version = "2.0.2"
+version = "2.1.0"
 authors = ["Olav Møyner <olav.moyner@gmail.com>"]
 
 [deps]

src/MultiComponentFlash.jl

@@ -59,6 +59,7 @@ module MultiComponentFlash
     include("static.jl")
     include("implicit_flash_derivatives.jl")
     using .ImplicitFlashDerivatives: implicit_flash_derivatives
+    include("static_newton.jl")
     include("tables.jl")
 
     include("flow_coupler.jl")

src/eos.jl

@@ -97,8 +97,8 @@ minimum_allowable_root(eos, forces, scalars) = 1e-16
 end
 
 @inline function root_bounds(roots, minimum_root)
-    max_root = -Inf
-    min_root = Inf
+    max_root = oftype(minimum_root, -Inf)
+    min_root = oftype(minimum_root, Inf)
     for root in roots
         max_root = max(max_root, root)
         if root > minimum_root
@@ -364,9 +364,9 @@ function solve_cubic_positive_roots(a, b, c)
     single_root = M > 0
     if single_root
         # Single real root
-        S = -sign(R)*(abs(R) + sqrt(M))^(1/3)
+        S = -sign(R)*cbrt(abs(R) + sqrt(M))
         if S == 0
-            T = 0
+            T = zero(S)
         else
             T = Q/S
         end
@@ -381,9 +381,10 @@ function solve_cubic_positive_roots(a, b, c)
         # repeated root, where acos is still well defined analytically.
         acos_arg = clamp(R/sqrt(Q^3), -one(R), one(R))
         theta = acos(acos_arg)
+        two_pi = oftype(theta, 2*pi)
         r1 = -(2*sqrt(Q)*cos(theta/3)) - a/3
-        r2 = -(2*sqrt(Q)*cos((theta + 2*pi)/3)) - a/3
-        r3 = -(2*sqrt(Q)*cos((theta - 2*pi)/3)) - a/3
+        r2 = -(2*sqrt(Q)*cos((theta + two_pi)/3)) - a/3
+        r3 = -(2*sqrt(Q)*cos((theta - two_pi)/3)) - a/3
         return (r1, r2, r3)
     end
 end

src/flash.jl

@@ -183,7 +183,7 @@ Pre-allocate storage for `flash_2ph!`.
 # Keyword arguments
 - `method = SSIFlash()`: Flash method to use. Can be `SSIFlash()`, `NewtonFlash()` or `SSINewtonFlash()`.
 - `static = false`: Return fully static, GPU-compatible storage when `true`. The
-  static path currently supports `SSIFlash`.
+  static path supports `SSIFlash` and `NewtonFlash`.
 - `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`.
 

src/implicit_flash_derivatives.jl

@@ -34,13 +34,8 @@ accelerator use, pass an immutable EOS from `make_eos_immutable`.
         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)
+    residual_numeric, jacobian = static_flash_system(eos, numeric_cond,
+        K_numeric, V_numeric)
 
     # The constant residual is subtracted so the returned primals remain the
     # converged numeric solution, even when its tolerance is finite.
@@ -54,6 +49,15 @@ accelerator use, pass an immutable EOS from `make_eos_immutable`.
     return updated[N+1], static_flash_K(updated)
 end
 
+@inline function static_flash_system(eos, numeric_cond, K, V)
+    # Seed only the equilibrium unknowns. The condition is numeric, so no
+    # nested AD types are needed. Shared by implicit derivatives and Newton.
+    u = static_flash_unknowns(K, V)
+    residual = static_equilibrium_residual(eos, numeric_cond,
+        seed_static_flash_unknowns(u))
+    return static_flash_residual_values(residual), static_flash_jacobian(residual)
+end
+
 """
     implicit_flash_derivatives(eos::KValuesEOS, numeric_cond, ad_cond,
                                V_numeric, K_numeric, K_ad)

src/static.jl

@@ -2,15 +2,16 @@
     StaticConfig()
 
 Marker returned by `flash_storage(...; static=true)`. It selects the immutable,
-stack-oriented SSI implementation used in accelerator kernels.
+stack-oriented flash 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."))
+    method isa Union{SSIFlash, NewtonFlash} || throw(ArgumentError(
+        "The static flash supports SSIFlash and NewtonFlash."))
     return config
 end
 
@@ -22,8 +23,8 @@ Run the immutable, accelerator-friendly two-phase flash implementation.
 scalar vapor fraction. When `storage` is omitted, a static storage marker is
 created automatically.
 
-Pass `V0=Inf` to skip stability testing and perform a negative flash, allowing
-the vapor fraction to converge outside `[0, 1]`.
+With `SSIFlash`, pass `V0=Inf` to skip stability testing and perform a negative
+flash, allowing the vapor fraction to converge outside `[0, 1]`.
 
 Set `return_stability=true` to additionally return a
 [`StaticStabilityResult`](@ref). Its immutable `storage` can be supplied as
@@ -33,8 +34,14 @@ and composition bounds as the mutable simulator integration.
 Set `stability_bypass=false` to force a fresh stability test even when storage
 is supplied; `return_stability=true` still records a new reference if eligible.
 
-The immutable path currently supports `SSIFlash` and generic cubic EOS values
-converted with [`make_eos_immutable`](@ref).
+The immutable path supports `SSIFlash` and `NewtonFlash` with generic cubic
+EOS values converted with [`make_eos_immutable`](@ref). The static Newton method
+minimizes Gibbs energy in component phase amounts, rejecting a nonpositive
+Hessian, nonpositive phase amounts or an increase in Gibbs energy. It returns
+`report.converged=false` on a rejected step or iteration limit. Use a short
+Newton attempt from a previous split and fall back to SSI on failure.
+Newton requires positive component amounts in both phases and uses static
+factorization for up to 24 components; larger systems return nonconvergence.
 """
 @inline function flash_2ph_immutable(eos, c, V0 = NaN;
         method = SSIFlash(),
@@ -367,7 +374,7 @@ end
 
 @inline function flash_2ph_impl!(::StaticConfig, K,
         eos::GenericCubicEOS{E, R, N}, c, V;
-        method::SSIFlash = SSIFlash(),
+        method::Union{SSIFlash, NewtonFlash} = SSIFlash(),
         maxiter::Int = 25000,
         tolerance::Float64 = 1e-8,
         verbose::Bool = false,
@@ -408,6 +415,9 @@ end
     converged = false
     if stable
         iteration = 0
+    elseif method isa NewtonFlash
+        V, K, converged, iteration = static_newton_flash(
+            eos, cond, K, V, tolerance, maxiter, method.dMax)
     else
         iteration = 1
         if isnan(V) || negative_flash

src/static_newton.jl

@@ -0,0 +1,126 @@
+# Phase-amount Gibbs minimization for the short reuse attempt described in
+# Belkadi et al. (2011), SPE 142132, Section 2.1. The caller handles fallback.
+@inline function static_newton_flash(eos, cond, K, V, tolerance, maxiter, dMax)
+    if isinf(V) || !(
+            all(k -> isfinite(k) && k > 0, K) &&
+                all(z -> isfinite(z) && z > 0, cond.z)
+        )
+        return V, K, false, 0
+    end
+    # Always restore the material-balance root before the first EOS evaluation.
+    # A caller's finite V may belong to a previous composition or set of K values.
+    initial_V = solve_rachford_rice(
+        K, cond.z, V;
+        tol = max(1.0e-12, eps(eltype(K)))
+    )
+    if !(isfinite(initial_V) && zero(initial_V) < initial_V < one(initial_V))
+        return initial_V, K, false, 0
+    end
+    liquid = map(cond.z, K) do zi, ki
+        (one(initial_V) - initial_V) * liquid_mole_fraction(zi, ki, initial_V)
+    end
+    vapor = map(cond.z, K) do zi, ki
+        initial_V * vapor_mole_fraction(liquid_mole_fraction(zi, ki, initial_V), ki)
+    end
+    liquid, vapor, converged, iterations = static_gibbs_newton(
+        eos, cond, liquid, vapor, tolerance; maxiter = maxiter, dMax = dMax
+    )
+    total_liquid, total_vapor = sum(liquid), sum(vapor)
+    next_K = (vapor / total_vapor) ./ (liquid / total_liquid)
+    return total_vapor / (total_liquid + total_vapor), next_K, converged, iterations
+end
+
+@inline function static_newton_phase_potentials(eos, cond, amounts, phase, forces)
+    composition = amounts / sum(amounts)
+    phase_cond = (p = cond.p, T = cond.T, z = composition, phase = phase)
+    Z, scalars = prep(eos, phase_cond, forces)
+    lnphi = static_fugacity_coefficients(
+        eos, phase_cond, Z, forces, scalars, eltype(amounts)
+    )
+    return log.(composition) + lnphi
+end
+
+@inline function static_newton_system(
+        eos, cond, liquid::SVector{N, F},
+        vapor::SVector{N, F}
+    ) where {N, F}
+    L, V = sum(liquid), sum(vapor)
+    K = (vapor / V) ./ (liquid / L)
+    residual, jacobian = ImplicitFlashDerivatives.static_flash_system(
+        eos, cond, K, V / (L + V)
+    )
+    gradient = SVector{N, F}(ntuple(i -> residual[i], Val(N)))
+    scale = sqrt.(liquid .* (vapor ./ cond.z))
+    # Reuse the (K, beta) Jacobian of the two-phase implicit derivatives.
+    # For v_i = z_i - l_i, d(log K_i)/dv_j = delta_ij*(1/v_i + 1/l_i)
+    # - 1/V - 1/L, and d(beta)/dv_j = 1/(L + V). Apply this chain rule
+    # directly to the scaled Hessian to avoid large trace-component reciprocals.
+    row_correction = SVector{N, F}(
+        ntuple(Val(N)) do i
+            total = sum(j -> jacobian[i, j] * K[j], 1:N)
+            jacobian[i, N + 1] / (L + V) - (inv(L) + inv(V)) * total
+        end
+    )
+    hessian = SMatrix{N, N, F}(
+        ntuple(Val(N * N)) do index
+            row = mod1(index, N)
+            col = (index - 1) ÷ N + 1
+            (scale[row] / scale[col]) * jacobian[row, col] * K[col] +
+                scale[row] * scale[col] * row_correction[row]
+        end
+    )
+    return gradient, (hessian + transpose(hessian)) / 2, scale
+end
+
+@inline static_positive_phase_amounts(amounts) = all(x -> isfinite(x) && x > 0, amounts)
+
+@inline function static_gibbs_newton(
+        eos, cond, liquid::SVector{N, F},
+        vapor::SVector{N, F}, tolerance; maxiter = 3, dMax = 1.0
+    ) where {N, F}
+    forces = static_force_coefficients(eos, cond, F)
+    # StaticArrays' small Cholesky is allocation-free and accelerator compatible.
+    # Larger systems retain the ordinary flash instead of a dynamic factorization.
+    N <= 24 || return liquid, vapor, false, 0
+    for iteration in 0:maxiter
+        if !(static_positive_phase_amounts(liquid) && static_positive_phase_amounts(vapor))
+            return liquid, vapor, false, iteration
+        end
+        gradient, H, scale = static_newton_system(eos, cond, liquid, vapor)
+        all(isfinite, gradient) && all(isfinite, H) ||
+            return liquid, vapor, false, iteration
+        factor = cholesky(Symmetric(H); check = false)
+        if !issuccess(factor) || !all(
+                i -> isfinite(factor.U[i, i]) &&
+                    factor.U[i, i] > 0, 1:N
+            )
+            return liquid, vapor, false, iteration
+        end
+        if maximum(abs, gradient) <= tolerance
+            return liquid, vapor, true, iteration
+        elseif iteration == maxiter
+            return liquid, vapor, false, iteration
+        end
+        increment = scale .* (factor \ (-scale .* gradient))
+        max_increment = maximum(abs, increment)
+        if max_increment > F(dMax)
+            increment *= F(dMax) / max_increment
+        end
+        next_liquid = liquid - increment
+        next_vapor = vapor + increment
+        if !(static_positive_phase_amounts(next_liquid) && static_positive_phase_amounts(next_vapor))
+            return liquid, vapor, false, iteration + 1
+        end
+        next_mu_l = static_newton_phase_potentials(eos, cond, next_liquid, Val(:liquid), forces)
+        next_mu_v = static_newton_phase_potentials(eos, cond, next_vapor, Val(:vapor), forces)
+        mu_l = static_newton_phase_potentials(eos, cond, liquid, Val(:liquid), forces)
+        mu_v = static_newton_phase_potentials(eos, cond, vapor, Val(:vapor), forces)
+        gibbs = sum(liquid .* mu_l) + sum(vapor .* mu_v)
+        next_gibbs = sum(next_liquid .* next_mu_l) + sum(next_vapor .* next_mu_v)
+        if !(isfinite(next_gibbs) && next_gibbs <= gibbs + 8eps(F) * max(1, abs(gibbs)))
+            return liquid, vapor, false, iteration + 1
+        end
+        liquid, vapor = next_liquid, next_vapor
+    end
+    return liquid, vapor, false, maxiter
+end

test/runtests.jl

@@ -1,6 +1,7 @@
 include("test_setup.jl")
 include("eos_equations.jl")
 include("implicit_derivatives.jl")
+include("static_newton.jl")
 
 @testset "Rachford-Rice" begin
     test_rachford_rice()