Modernize docs + EOS fixes + Clapeyron comparisons

docs/Project.toml

@@ -1,5 +1,7 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Clapeyron = "7c7805af-46cc-48c9-995b-ed0ed2dc909a"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 MultiComponentFlash = "35e5bd01-9722-4017-9deb-64a5d32478ff"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/make.jl

@@ -1,12 +1,15 @@
 using MultiComponentFlash
 using Documenter
-using BenchmarkTools, Plots
+using BenchmarkTools, CairoMakie
+
+# CairoMakie.activate!(type = "svg")
 
 DocMeta.setdocmeta!(MultiComponentFlash, :DocTestSetup, :(using MultiComponentFlash); recursive=true)
 
 makedocs(;
     modules=[MultiComponentFlash],
-    warnonly = true,
+    warnonly = false,
+    checkdocs=:exports,
     authors="Olav Møyner <olav.moyner@gmail.com>",
     repo="https://github.com/moyner/MultiComponentFlash.jl/blob/{commit}{path}#{line}",
     sitename="MultiComponentFlash.jl",
@@ -19,7 +22,8 @@ makedocs(;
         "Home" => "index.md",
         "Examples" => Any[
             "Basic usage" => "examples/basics.md",
-            "Advanced usage" => "examples/advanced.md"
+            "Advanced usage" => "examples/advanced.md",
+            "EOS reference validation" => "examples/eos_validation.md",
         ],
         "API" => Any[
             "Mixtures" => "api/mixtures.md",

docs/src/examples/advanced.md

@@ -86,9 +86,11 @@ between kernel work items.
 
 We create a three-component mixture and flash for a range of pressure and temperature conditions:
 
-```julia
-using MultiComponentFlash, Plots
-ns = 1000
+```@example phase-diagram
+using MultiComponentFlash, CairoMakie
+CairoMakie.activate!(type = "svg")
+
+ns = 100
 ubar = 1e5
 # Pressure range
 p0 = 1*ubar
@@ -110,20 +112,29 @@ cond = (p = p0, T = T0, z = z)
 m = SSIFlash()
 S = flash_storage(eos, cond, method = m)
 K = initial_guess_K(eos, cond)
-data = zeros(ns, ns)
-for ip = 1:ns
-    for iT = 1:ns
-        c = (p = p[ip], T = T[iT], z = z)
-        data[ip, iT] = flash_2ph!(S, K, eos, c, NaN, method = m)
+data = zeros(length(T), length(p))
+for (iT, temperature) in pairs(T)
+    for (ip, pressure) in pairs(p)
+        c = (p = pressure, T = temperature, z = z)
+        V = flash_2ph!(S, K, eos, c, NaN, method = m)
+        data[iT, ip] = V
     end
 end
 
-contour(p./ubar, T .- 273.15, data, levels = 10, fill=(true,cgrad(:hot)))
-ylabel!("Pressure [Bar]")
-xlabel!("T [°Celsius]")
+fig = Figure(size = (760, 480))
+ax = Axis(fig[1, 1];
+    xlabel = "Temperature [°C]",
+    ylabel = "Pressure [bar]",
+    title = "Vapor fraction")
+contours = contourf!(ax, T .- 273.15, p./ubar, data;
+    levels = range(0.0, 1.0, length = 11), colormap = :hot)
+Colorbar(fig[1, 2], contours; label = "Vapor mole fraction")
+fig
 ```
 
-![Phase diagram](../assets/phase_diagram_simple.png)
+The colored region is the two-phase envelope. Unfilled states are stable
+single-phase conditions, for which the two-phase solver reports no intermediate
+vapor fraction.
 
 ### PVT table generation
 
@@ -132,3 +143,11 @@ There is experimental support for generating simulator input blackoil tables (e.
 ```@docs
 generate_pvt_tables
 ```
+
+### Coupling to simulators and other utilities
+
+```@docs
+cubic_benchmark
+FlashedMixture2Phase
+FlashedPhase
+```

docs/src/examples/eos_validation.md

@@ -0,0 +1,175 @@
+# EOS reference validation
+
+This example compares Peng-Robinson 1976 (PR76), Soave-Redlich-Kwong (SRK), and
+Redlich-Kwong (RK) with
+[Clapeyron.jl](https://github.com/ClapeyronThermo/Clapeyron.jl). It checks
+compressibility and fugacity coefficients for a set of cases.
+
+## Comparison matrix
+
+```@example eos-reference
+using Clapeyron, MultiComponentFlash, Test
+
+names = ["carbon dioxide", "methane", "decane"]
+properties = (
+    MolecularProperty(0.0440, 7.38e6, 304.1, 9.412e-5, 0.224),
+    MolecularProperty(0.0160, 4.60e6, 190.6, 9.863e-5, 0.011),
+    MolecularProperty(0.1420, 2.10e6, 617.7, 6.098e-4, 0.488),
+)
+conditions = (
+    (p = 1.0e6, T = 300.0, z = [0.5, 0.3, 0.2]),
+    (p = 5.0e6, T = 350.0, z = [0.2, 0.7, 0.1]),
+    (p = 2.0e7, T = 500.0, z = [0.1, 0.2, 0.7]),
+)
+interactions = (
+    zeros(3, 3),
+    [0.0 0.08 0.03; 0.08 0.0 0.015; 0.03 0.015 0.0],
+)
+
+function compare_cubic_eos(names, properties, conditions, interactions)
+    relative_error(value, reference) =
+        abs(value - reference)/max(abs(reference), eps(Float64))
+    Tc = [property.T_c for property in properties]
+    Pc = [property.p_c for property in properties]
+    Mw = [1000property.mw for property in properties]
+    acentricfactor = [property.ω for property in properties]
+    comparisons = NamedTuple[]
+
+    for k in interactions
+        mixture = MultiComponentMixture(properties; names = names, A_ij = k)
+        parameters = (; Tc, Pc, Mw, acentricfactor, k, l = zeros(size(k)))
+        models = (
+            (label = "PR76", eos = PengRobinson(),
+                reference = Clapeyron.PR(names;
+                    userlocations = parameters, verbose = false)),
+            (label = "SRK", eos = SoaveRedlichKwong(),
+                reference = Clapeyron.SRK(names;
+                    userlocations = parameters, verbose = false)),
+            (label = "RK", eos = RedlichKwong(),
+                reference = Clapeyron.RK(names;
+                    userlocations = parameters, verbose = false)),
+        )
+
+        for model in models, state in conditions
+            eos = GenericCubicEOS(mixture, model.eos)
+            for (phase, reference_phase) in ((:liquid, :l), (:vapor, :v))
+                state_with_phase = (;
+                    p = state.p, T = state.T, z = state.z, phase)
+                forces = force_coefficients(eos, state_with_phase)
+                scalars = force_scalars(eos, state_with_phase, forces)
+                Z = mixture_compressibility_factor(
+                    eos, state_with_phase, forces, scalars)
+                phi = [exp(MultiComponentFlash.component_fugacity_coefficient(
+                    eos, state_with_phase, i, Z, forces, scalars))
+                    for i in eachindex(state.z)]
+
+                reference_volume = Clapeyron.volume(
+                    model.reference, state.p, state.T, state.z;
+                    phase = reference_phase)
+                reference_Z = state.p*reference_volume/
+                    (Clapeyron.Rgas(model.reference)*state.T*sum(state.z))
+                reference_phi = Clapeyron.fugacity_coefficient(
+                    model.reference, state.p, state.T, state.z;
+                    phase = reference_phase)
+
+                push!(comparisons, (;
+                    eos = model.label,
+                    Z,
+                    reference_Z,
+                    Z_relative_error = relative_error(Z, reference_Z),
+                    phi,
+                    reference_phi,
+                    phi_relative_error = relative_error.(phi, reference_phi)))
+            end
+        end
+    end
+    return comparisons
+end
+
+comparisons = compare_cubic_eos(
+    names, properties, conditions, interactions)
+
+@test length(comparisons) == 36
+@test all(row -> isapprox(row.Z, row.reference_Z;
+    rtol = 1.0e-7, atol = 1.0e-10), comparisons)
+@test all(row -> isapprox(row.phi, row.reference_phi;
+    rtol = 1.0e-7, atol = 1.0e-10), comparisons)
+
+length(comparisons)
+```
+
+## Parity plots
+
+The dashed line indicates exact agreement. Fugacity coefficients use logarithmic
+axes because they span a wider range than ``Z``.
+
+```@example eos-reference
+using CairoMakie
+# CairoMakie.activate!(type = "svg")
+
+eos_names = ["PR76", "SRK", "RK"]
+colors = [:dodgerblue3, :darkorange2, :seagreen4]
+
+fig = Figure(size = (920, 410))
+z_axis = Axis(fig[1, 1];
+    xlabel = "Clapeyron Z",
+    ylabel = "MultiComponentFlash Z",
+    title = "Compressibility parity",
+    aspect = DataAspect())
+phi_axis = Axis(fig[1, 2];
+    xlabel = "Clapeyron fugacity coefficient",
+    ylabel = "MultiComponentFlash fugacity coefficient",
+    title = "Fugacity-coefficient parity",
+    xscale = log10,
+    yscale = log10,
+    aspect = DataAspect())
+
+for (eos_name, color) in zip(eos_names, colors)
+    rows = filter(row -> row.eos == eos_name, comparisons)
+    reference_z = [row.reference_Z for row in rows]
+    calculated_z = [row.Z for row in rows]
+    reference_phi = [value for row in rows for value in row.reference_phi]
+    calculated_phi = [value for row in rows for value in row.phi]
+    scatter!(z_axis, reference_z, calculated_z;
+        color = color, markersize = 9, label = eos_name)
+    scatter!(phi_axis, reference_phi, calculated_phi;
+        color = color, markersize = 9, label = eos_name)
+end
+
+z_limits = extrema(vcat(
+    [row.reference_Z for row in comparisons],
+    [row.Z for row in comparisons]))
+phi_limits = extrema(vcat(
+    [value for row in comparisons for value in row.reference_phi],
+    [value for row in comparisons for value in row.phi]))
+z_line = collect(z_limits)
+phi_line = collect(phi_limits)
+lines!(z_axis, z_line, z_line; color = :black, linestyle = :dash)
+lines!(phi_axis, phi_line, phi_line; color = :black, linestyle = :dash)
+axislegend(z_axis; position = :lt)
+fig
+```
+
+The maximum relative error makes the remaining differences visible.
+
+```@example eos-reference
+maximum_z_error = [maximum(row.Z_relative_error
+    for row in comparisons if row.eos == eos_name) for eos_name in eos_names]
+maximum_phi_error = [maximum(maximum(row.phi_relative_error)
+    for row in comparisons if row.eos == eos_name) for eos_name in eos_names]
+
+error_figure = Figure(size = (720, 420))
+error_axis = Axis(error_figure[1, 1];
+    xlabel = "Equation of state",
+    ylabel = "Maximum relative error",
+    title = "Worst case across all states, roots, and BIC matrices",
+    yscale = log10,
+    xticks = (1:length(eos_names), eos_names))
+scatterlines!(error_axis, 1:length(eos_names), maximum_z_error;
+    markersize = 12, label = "Z")
+scatterlines!(error_axis, 1:length(eos_names), maximum_phi_error;
+    markersize = 12,
+    label = "Fugacity coefficients")
+axislegend(error_axis; position = :lt)
+error_figure
+```

docs/src/index.md

@@ -11,9 +11,12 @@ This package implements several equations of state for multicomponent vapor-liqu
  The following equations of state (EOS) are implemented as a class of generic cubics:
 
 * [Peng-Robinson](https://doi.org/10.1021/i160057a011)
+* Corrected Peng-Robinson for large acentric factors
 * [Redlich-Kwong](https://doi.org/10.1021/cr60137a013)
 * [Soave-Redlich-Kwong](https://doi.org/10.1016/0009-2509(72)80096-4)
 * [Zudkevitch-Joffe](https://doi.org/10.1002/aic.690160122)
+* [Soreide-Whitson](https://doi.org/10.1016/0378-3812(92)85105-H) for
+  water and brine mixtures
 
 The code is fully type stable, easy to use and fairly performant, with additional options to avoid allocations if you need to perform many flashes. The main implementation goal is to have a compact, performant and easy to use code suitable for integration in simulators of multiphase flow.
 
@@ -38,15 +41,17 @@ If you want to use the former features, please know that they might be subject t
 * The module currently only supports cubic equations of state. These are limited in accuracy for longer chains of molecules without extensive tuning for a specific mixture.
 * Flash is limited to two-phase pressure-temperature (pT) flash.
 * The flash algorithms are limited to the basics - there are many strategies that could be implemented
-* Peng-Robinson is the only equation of state that has been thoroughly validated.
+* Cubic equations remain approximations whose practical accuracy depends on
+  component data and fitted binary-interaction coefficients. See
+  [EOS reference validation](@ref) for implementation-level comparisons.
 
 # Other packages
 Julia packages:
-* [Clapeyron.jl](https://github.com/ypaul21/Clapeyron.jl) supports more advanced equations of state in addition to the cubics (SAFT-type and empirical EOS) and has a much larger API for thermodynamical properties. At the time of writing, this package does not support a full flash.
+* [Clapeyron.jl](https://github.com/ypaul21/Clapeyron.jl) supports more advanced equations of state in addition to the cubics (SAFT-type and empirical EOS) and has a much larger API for thermodynamical properties. It is also used by this package's independent [EOS reference validation](@ref) suite.
 
 Julia interfaces to other useful packages:
 * [PyThermo.jl](https://github.com/stillyslalom/PyThermo.jl) is an interface to the [Thermo](https://pypi.org/project/thermo/) package.
 * [CoolProp.jl](https://github.com/CoolProp/CoolProp.jl) is an interface to the [CoolProp](http://www.coolprop.org/) library.
 
 # Contact
-You can use the [Github webpage](https://github.com/moyner/MultiComponentFlash.jl) or drop me a line at [Olav Møyner](mailto:olav.moyner@gmail.no).
+You can use the [Github webpage](https://github.com/moyner/MultiComponentFlash.jl) or drop me a line at [Olav Møyner](mailto:olav.moyner@sintef.no).

src/PVTExperiments/PVTExperiments.jl

@@ -21,5 +21,5 @@ module PVTExperiments
     include("tables.jl")
     include("interface.jl")
 
-    export generate_pvt_tables, PVTTableSet
+    export generate_pvt_tables
 end

src/eos.jl

@@ -159,18 +159,17 @@ function force_coefficients(eos::AbstractCubicEOS, cond; static_size = false)
     return coeff
 end
 
-function get_force_coefficients(forces, eos, cond)
-    if forces_per_phase(eos)
-        phase = get_phase(cond)
-        if phase == :liquid
-            return forces.liquid
-        elseif phase == :vapor
-            return forces.vapor
-        else
-            error("Forces per phase are only supported for liquid and vapor phases, not $phase.")
-        end
+@inline get_force_coefficients(forces, eos, cond) = forces
+
+@inline function get_force_coefficients(
+        forces::NamedTuple{(:liquid, :vapor)}, eos, cond)
+    phase = get_phase(cond)
+    if phase == :liquid
+        return forces.liquid
+    elseif phase == :vapor
+        return forces.vapor
     else
-        return forces
+        error("Forces per phase are only supported for liquid and vapor phases, not $phase.")
     end
 end
 
@@ -363,17 +362,24 @@ function solve_cubic_positive_roots(a, b, c)
     M = R^2 - Q^3
     single_root = M > 0
     if single_root
-        # Single real roots
+        # Single real root
         S = -sign(R)*(abs(R) + sqrt(M))^(1/3)
         if S == 0
             T = 0
         else
             T = Q/S
         end
         return S + T - a/3
+    elseif iszero(Q)
+        # Q = R = 0: the polynomial has one triple root. The trigonometric
+        # expression below is otherwise undefined (acos(0/0)).
+        return -a/3
     else
         # Three real roots
-        theta = acos(R/sqrt(Q^3))
+        # Clamp guards against a round-off excursion outside [-1, 1] at a
+        # repeated root, where acos is still well defined analytically.
+        acos_arg = clamp(R/sqrt(Q^3), -one(R), one(R))
+        theta = acos(acos_arg)
         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

src/eos/soreide_whitson.jl

@@ -87,5 +87,5 @@ function soreide_whitson_hc_aqueous_bic(sw::SoreideWhitson, acf_i, T_ri)
     A_0 = a_0 + a_mw_0*sign(acf_i)*abs(acf_i)^(-0.1)
     A_1 = a_1 + a_mw_1*acf_i
     A_2 = a_2 + a_mw_2*acf_i
-    return A_0*(1 + α_0*c_sw) + A_1*T_ri*(1 + α_1*c_sw) + A_2*T_ri*(1 + α_2*c_sw)
+    return A_0*(1 + α_0*c_sw) + A_1*T_ri*(1 + α_1*c_sw) + A_2*T_ri^2*(1 + α_2*c_sw)
 end

src/eos/zudkevitch_joffe.jl

@@ -1,13 +1,13 @@
 # ZudkevitchJoffe
-function weight_ai(eos::GenericCubicEOS{ZudkevitchJoffe}, cond, i)
+function weight_ai(eos::GenericCubicEOS{E}, cond, i) where E<:ZudkevitchJoffe
     zj = eos.type
     mix = eos.mixture
     T = cond.T
     T_r = reduced_temperature(mix, cond, i)
     return eos.ω_a*zj.F_a(T, i)*T_r^(-0.5)
 end
 
-function weight_bi(eos::GenericCubicEOS{ZudkevitchJoffe}, cond, i)
+function weight_bi(eos::GenericCubicEOS{E}, cond, i) where E<:ZudkevitchJoffe
     zj = eos.type
     T = cond.T
     return eos.ω_b*zj.F_b(T, i)