diff --git a/Project.toml b/Project.toml index 4e80756..4dd4e8d 100644 --- a/Project.toml +++ b/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 "] [deps] diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl index 10a5e6d..63fff4e 100644 --- a/src/MultiComponentFlash.jl +++ b/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") diff --git a/src/eos.jl b/src/eos.jl index 1456a7e..0c841cc 100644 --- a/src/eos.jl +++ b/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 diff --git a/src/flash.jl b/src/flash.jl index 0f85dc2..a021888 100644 --- a/src/flash.jl +++ b/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`. diff --git a/src/implicit_flash_derivatives.jl b/src/implicit_flash_derivatives.jl index e8fcba2..a49d923 100644 --- a/src/implicit_flash_derivatives.jl +++ b/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) diff --git a/src/static.jl b/src/static.jl index dbc3181..924e6b2 100644 --- a/src/static.jl +++ b/src/static.jl @@ -2,7 +2,7 @@ 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 @@ -10,7 +10,8 @@ struct StaticConfig end @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 diff --git a/src/static_newton.jl b/src/static_newton.jl new file mode 100644 index 0000000..1c6463d --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 7ab2dcd..211c0c4 100644 --- a/test/runtests.jl +++ b/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() diff --git a/test/static_newton.jl b/test/static_newton.jl new file mode 100644 index 0000000..55de34c --- /dev/null +++ b/test/static_newton.jl @@ -0,0 +1,120 @@ +using Test, MultiComponentFlash, StaticArrays, ForwardDiff, LinearAlgebra +using KernelAbstractions, JLArrays + +@kernel function static_newton_kernel!(out, converged, eos, cond, K0, V0, tolerance) + i = @index(Global) + V, _, report = flash_2ph!( + MultiComponentFlash.StaticConfig(), K0, eos, cond, V0; + method = NewtonFlash(dMax = 1.0), maxiter = 3, extra_out = true, + tolerance = tolerance, z_min = nothing, check = false + ) + @inbounds out[i] = V + @inbounds converged[i] = report.converged +end + +static_newton_allocations(eos, cond, K, V) = @allocated flash_2ph!( + MultiComponentFlash.StaticConfig(), K, eos, cond, V; + method = NewtonFlash(dMax = 1.0), maxiter = 3, extra_out = true, + z_min = nothing, check = false +) + +@testset "Static Gibbs Newton flash" begin + for F in (Float64, Float32) + eos = make_eos_immutable(get_test_eos(); float_type = F) + cond = (p = F(1.0e6), T = F(300), z = SVector{3, F}(0.5, 0.3, 0.2)) + reference_eos = make_eos_immutable(get_test_eos()) + V_reference, K_reference = flash_2ph_immutable( + reference_eos, + (p = Float64(cond.p), T = Float64(cond.T), z = SVector{3, Float64}(cond.z)) + ) + V0, K0 = F(V_reference), SVector{3, F}(K_reference) + nearby = ( + p = F(1.0001e6), T = F(300.01), + z = SVector{3, F}(0.50001, 0.29999, 0.2), + ) + tolerance = F == Float32 ? 1.0e-4 : 1.0e-8 + V, K, report = flash_2ph!( + MultiComponentFlash.StaticConfig(), + K0, eos, nearby, V0; method = NewtonFlash(dMax = 1.0), + maxiter = 3, extra_out = true, tolerance = tolerance, + z_min = nothing, check = false + ) + expected_V, expected_K = flash_2ph_immutable( + reference_eos, + (p = Float64(nearby.p), T = Float64(nearby.T), z = SVector{3, Float64}(nearby.z)) + ) + @test report.converged + @test report.its <= 3 + @test V isa F + @test K isa SVector{3, F} + @test V ≈ expected_V atol = 2tolerance + @test K ≈ expected_K rtol = 10tolerance + x = nearby.z ./ (1 .+ V .* (K .- 1)) + y = K .* x + @test sum(x) ≈ one(F) atol = 10eps(F) + @test sum(y) ≈ one(F) atol = 10eps(F) + @test (1 - V) * x + V * y ≈ nearby.z + if isdefined(JLArrays, :JLBackend) + out = JLArray(zeros(F, 8)) + converged = JLArray(falses(8)) + backend = JLArrays.JLBackend() + static_newton_kernel!(backend, 8)( + out, converged, eos, nearby, + K0, V0, tolerance; ndrange = 8 + ) + @test all(Array(converged)) + @test all(isapprox.(Array(out), V; atol = tolerance, rtol = 0)) + end + if F == Float64 + static_newton_allocations(eos, nearby, K0, V0) + @test static_newton_allocations(eos, nearby, K0, V0) == 0 + end + + # A finite but stale V is normalized before evaluating fugacities. + other_V, other_K, other_report = flash_2ph!( + MultiComponentFlash.StaticConfig(), + K0, eos, nearby, F(0.1); method = NewtonFlash(dMax = 1.0), + maxiter = 3, extra_out = true, tolerance = tolerance, + z_min = nothing, check = false + ) + @test other_report.converged + @test other_V ≈ V + @test other_K ≈ K + + _, _, unfinished = flash_2ph!( + MultiComponentFlash.StaticConfig(), + K0, eos, nearby, V0; method = NewtonFlash(dMax = 1.0), + maxiter = 0, extra_out = true, tolerance = 1.0e-12, + z_min = nothing, check = false + ) + @test !unfinished.converged + _, _, trivial = flash_2ph!( + MultiComponentFlash.StaticConfig(), + one.(K0), eos, nearby, V0; method = NewtonFlash(), + maxiter = 3, extra_out = true, z_min = nothing, check = false + ) + @test !trivial.converged + end +end + +@testset "Gibbs Hessian from the two-phase derivative Jacobian" begin + MCF = MultiComponentFlash + eos = make_eos_immutable(get_test_eos()) + cond = (p = 1.0e6, T = 300.0, z = @SVector [0.5, 0.3, 0.2]) + V, K = flash_2ph_immutable(eos, cond) + # Check away from the stationary solution so the residual and its + # derivatives are both nonzero. The component amounts still sum to z. + vapor = V .* K .* (cond.z ./ (1 .+ V .* (K .- 1))) + vapor += @SVector [0.001, -0.0001, 0.0002] + liquid = cond.z - vapor + gradient, hessian, scale = MCF.static_newton_system(eos, cond, liquid, vapor) + forces = MCF.static_force_coefficients(eos, cond, Float64) + potentials(v) = MCF.static_newton_phase_potentials(eos, cond, v, Val(:vapor), forces) - + MCF.static_newton_phase_potentials(eos, cond, cond.z - v, Val(:liquid), forces) + independent_hessian = ForwardDiff.jacobian(potentials, vapor) + @test gradient ≈ potentials(vapor) rtol = 1.0e-10 atol = 1.0e-12 + @test hessian ≈ Diagonal(scale) * independent_hessian * Diagonal(scale) rtol = 1.0e-10 + @test isposdef(Symmetric(hessian)) + @test scale .* (hessian \ (-scale .* gradient)) ≈ + independent_hessian \ (-gradient) rtol = 1.0e-10 +end