Flash bypass and improved performance for many components

Project.toml

@@ -56,7 +56,7 @@ Krylov = "0.9.1"
 LinearAlgebra = "1"
 LoopVectorization = "0.12.115"
 MAT = "0.10.3"
-MultiComponentFlash = "1.1.10"
+MultiComponentFlash = "1.1.13"
 OrderedCollections = "1.6.2"
 PartitionedArrays = "0.3.2"
 Polyester = "0.6.11, 0.7.3"

src/CO2Properties/setup.jl

@@ -1,10 +1,10 @@
-function setup_reservoir_model_co2_brine(reservoir::DataDomain,
-    T = missing;
-    thermal = false,
-    composite = thermal,
-    kwarg...
+function setup_reservoir_model_co2_brine(reservoir::DataDomain;
+        temperature = missing,
+        thermal = false,
+        composite = thermal,
+        kwarg...
     )
-    tables = co2_brine_property_tables(T)
+    tables = co2_brine_property_tables(temperature)
     rho = JutulDarcy.BrineCO2MixingDensities(tables[:density])
     mu = JutulDarcy.PTViscosities(tables[:viscosity])
     if thermal

src/NLDD/simulator.jl

@@ -335,7 +335,11 @@ function get_solve_status(rep, ok)
             if length(rep) > 1
                 status = local_solved_in_multiple_steps
             else
-                if length(only(rep)[:steps]) == 0
+                rep = only(rep)
+                nstep = length(rep[:steps])
+                if nstep == 0
+                    status = local_already_converged
+                elseif nstep == 1 && !haskey(rep[:steps][1], :linear_solver)
                     status = local_already_converged
                 else
                     status = local_solved_in_single_step
@@ -673,6 +677,10 @@ function solve_subdomain(sim, i, dt, forces, cfg)
         ok, report = Jutul.solve_timestep!(sim, dt, f, max_iter, cfg,
             finalize = true, prepare = true, update_explicit = false)
     catch e
+        if e isa InterruptException
+            # Don't leave the user trapped...
+            rethrow(e)
+        end
         @warn "Subdomain $i threw exception: $e"
     end
     return (ok, report)

src/NLDD/utils.jl

@@ -5,6 +5,7 @@ import JutulDarcy: set_default_cnv_mb!
 
 function simulator_config(sim::NLDDSimulator;
         method = :nldd,
+        subdomain_info_level = -1,
         inner_tol_mul = 1.0,
         inner_tol_final = 10.0,
         inner_max_timestep_cuts = 2,
@@ -19,6 +20,9 @@ function simulator_config(sim::NLDDSimulator;
         inc_tol_dz = Inf,
         kwarg...
     )
+    if subdomain_info_level isa Real
+        subdomain_info_level = i -> i == subdomain_info_level
+    end
     check_before_solve = isinf(inc_tol_dp_abs) && isinf(inc_tol_dp_rel) && isinf(inc_tol_dz)
     cfg = simulator_config(sim.simulator)
     set_default_cnv_mb!(cfg, sim.simulator.model,
@@ -68,7 +72,7 @@ function simulator_config(sim::NLDDSimulator;
             min_nonlinear_iterations = inner_min_nonlinear_iterations,
             tol_factor_final_iteration = inner_tol_final,
             check_before_solve = check_before_solve,
-            info_level = -1
+            info_level = Int(subdomain_info_level(i))
         )
         # Small hack for reservoir stuff
         set_default_cnv_mb!(subconfigs[i], subsims[i].model, tol_cnv = inner_tol_mul*tol_cnv, 

src/cpr.jl

@@ -367,7 +367,7 @@ function true_impes!(w, acc, r, n, bz, arg...)
     elseif bz == 8
         true_impes_8!(w, acc, r, n, bz, arg...)
     else
-        true_impes_gen!(w, acc, r, n, bz, arg...)
+        true_impes_gen!(w, acc, r, n, Val(bz), arg...)
     end
 end
 
@@ -447,9 +447,10 @@ function true_impes_8!(w, acc, r, n, bz, s, scaling)
     end
 end
 
-function true_impes_gen!(w, acc, r, n, bz, p_scale, scaling)
+function true_impes_gen!(w, acc, r, n, ::Val{bz}, p_scale, scaling) where bz
     r_p = SVector{bz}(r)
-    A = MMatrix{bz, bz, eltype(r)}(zeros(bz, bz))
+    r_T = eltype(r)
+    A = MMatrix{bz, bz, r_T}(undef)
     for cell in 1:n
         @inbounds for i = 1:bz
             v = acc[i, cell]
@@ -458,7 +459,7 @@ function true_impes_gen!(w, acc, r, n, bz, p_scale, scaling)
                 A[j, i] = v.partials[j]
             end
         end
-        invert_w!(w, A, r_p, cell, bz, scaling)
+        invert_w!(w, SMatrix{bz, bz, r_T}(A), r_p, cell, bz, scaling)
     end
 end
 

src/facility/well_presolve.jl

@@ -38,7 +38,7 @@ end
 
 function Jutul.simulator_config!(cfg, sim, ::PrepareStepWellSolver, storage)
     add_option!(cfg, :well_iterations, 25, "Well iterations to be performed before step.", types = Int, values = 0:10000)
-    add_option!(cfg, :well_outer_iterations_limit, Inf, "Number of outer iterations for each solve where wells are to be solved.", types = Union{Int, Float64})
+    add_option!(cfg, :well_outer_iterations_limit, 1, "Number of outer iterations for each solve where wells are to be solved.", types = Union{Int, Float64})
     add_option!(cfg, :well_acceptance_factor, 10.0, "Accept well pre-solve results at this relaxed factor.", types = Float64)
     add_option!(cfg, :well_info_level, -1, "Info level for well solver.", types = Int)
 end

src/flux.jl

@@ -89,7 +89,8 @@ function face_average_density(model::Union{CompositionalModel, ThermalCompositio
         elseif s_r_tiny
             ρ_avg = ρ_l
         else
-            ρ_avg = (s_l*ρ_r + s_r*ρ_l)/(s_l + s_r)
+            # alt def: (s_l*ρ_r + s_r*ρ_l)/(s_l + s_r)
+            ρ_avg = 0.5*(ρ_r + ρ_l)
         end
     end
     return ρ_avg

src/init/init.jl

@@ -5,6 +5,7 @@ function equilibriate_state(model, contacts,
         cells = missing,
         rs = missing,
         rv = missing,
+        composition = missing,
         kwarg...
     )
     model = reservoir_model(model)
@@ -25,7 +26,7 @@ function equilibriate_state(model, contacts,
 
     init = Dict{Symbol, Any}()
     init = equilibriate_state!(init, pts, model, sys, contacts, datum_depth, datum_pressure;
-        cells = cells, rv = rv, rs = rs, kwarg...
+        cells = cells, rv = rv, rs = rs, composition = composition, kwarg...
     )
 
     is_blackoil = sys isa StandardBlackOilSystem
@@ -51,13 +52,29 @@ function equilibriate_state(model, contacts,
             init[:Rv] = Rv
         end
     end
+    is_compositional = sys isa CompositionalSystem
+    if is_compositional
+        @assert !ismissing(composition)
+        has_wat = has_other_phase(sys)
+        nc = length(init[:Pressure])
+        ncomp = number_of_components(sys) - has_wat
+        zmf = zeros(ncomp, nc)
+        for i in 1:nc
+            zmf[:, i] = composition(pts[i])
+        end
+        init[:OverallMoleFractions] = zmf
+        if has_wat
+            init[:ImmiscibleSaturation] = init[:Saturations][1, :]
+        end
+    end
     return init
 end
 
 function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pressure;
         cells = 1:length(depths),
         rs = missing,
         rv = missing,
+        composition = missing,
         T_z = missing,
         s_min = missing,
         contacts_pc = missing,
@@ -95,24 +112,43 @@ function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pr
         relperm = last(relperm)
     end
 
-    pvt = rho.pvt
     reg = Int[pvtnum]
     function density_f(p, z, ph)
-        pvt_i = pvt[ph]
-        if phases[ph] == LiquidPhase() && disgas
-            rs_max = table_by_region(sys.rs_max, pvtnum)
-            Rs = min(rs(z), rs_max(p))
-            b = JutulDarcy.shrinkage(pvt_i, reg, p, Rs, 1)
-            rho = b*(rhoOS + Rs*rhoGS)
-        elseif phases[ph] == VaporPhase() && vapoil
-            rv_max = table_by_region(sys.rv_max, pvtnum)
-            Rv = min(rv(z), rv_max(p))
-            b = JutulDarcy.shrinkage(pvt_i, reg, p, Rv, 1)
-            rho = b*(rhoGS + Rv*rhoOS)
+        if rho isa ThreePhaseCompositionalDensitiesLV || rho isa TwoPhaseCompositionalDensities
+            if phases[ph] == AqueousPhase()
+                phase_density = rho_s[ph]*JutulDarcy.shrinkage(rho.immiscible_pvt, p)
+            else
+                @assert !ismissing(composition) "Composition must be present for equilibrium."
+                @assert !ismissing(T_z) "Temperature must be present for equilibrium."
+                z_i = Vector{Float64}(composition(z))
+                T = T_z(z)
+                eos = model.system.equation_of_state
+                f = MultiComponentFlash.flashed_mixture_2ph(eos, (p = p, T = T, z = z_i))
+                rho_l, rho_v = mass_densities(eos, p, T, f)
+                if phases[ph] == VaporPhase() || rho_l == 0
+                    phase_density = rho_v
+                else
+                    phase_density = rho_l
+                end
+            end
         else
-            rho = rho_s[ph]*JutulDarcy.shrinkage(pvt_i, reg, p, 1)
+            pvt = rho.pvt
+            pvt_i = pvt[ph]
+            if phases[ph] == LiquidPhase() && disgas
+                rs_max = table_by_region(sys.rs_max, pvtnum)
+                Rs = min(rs(z), rs_max(p))
+                b = JutulDarcy.shrinkage(pvt_i, reg, p, Rs, 1)
+                phase_density = b*(rhoOS + Rs*rhoGS)
+            elseif phases[ph] == VaporPhase() && vapoil
+                rv_max = table_by_region(sys.rv_max, pvtnum)
+                Rv = min(rv(z), rv_max(p))
+                b = JutulDarcy.shrinkage(pvt_i, reg, p, Rv, 1)
+                phase_density = b*(rhoGS + Rv*rhoOS)
+            else
+                phase_density = rho_s[ph]*JutulDarcy.shrinkage(pvt_i, reg, p, 1)
+            end
         end
-        return rho
+        return phase_density
     end
     pressures = determine_hydrostatic_pressures(depths, depth, zmin, zmax, contacts, datum_pressure, density_f, contacts_pc)
     if nph > 1
@@ -207,6 +243,8 @@ function parse_state0_equil(model, datafile)
 
     init = Dict{Symbol, Any}()
     sol = datafile["SOLUTION"]
+    props = datafile["PROPS"]
+
     G = physical_representation(model.data_domain)
     nc = number_of_cells(G)
     nph = number_of_phases(model.system)
@@ -220,12 +258,17 @@ function parse_state0_equil(model, datafile)
     npvt = GeoEnergyIO.InputParser.number_of_tables(datafile, :pvtnum)
     nsat = GeoEnergyIO.InputParser.number_of_tables(datafile, :satnum)
 
-    if haskey(sol, "RTEMP")
-        Ti = convert_to_si(only(sol["RTEMP"]), :Celsius)
+    if haskey(sol, "RTEMP") || haskey(props, "RTEMP")
+        if haskey(props, "RTEMP")
+            rtmp = only(props["RTEMP"])
+        else
+            rtmp = only(sol["RTEMP"])
+        end
+        Ti = convert_to_si(rtmp, :Celsius)
         T_z = z -> Ti
-    elseif haskey(sol, "TEMPVD")
-        z = vec(sol["TEMPVD"][:, 1])
-        Tvd = vec(sol["TEMPVD"][:, 2] + 273.15)
+    elseif haskey(props, "TEMPVD")
+        z = vec(props["TEMPVD"][:, 1])
+        Tvd = vec(props["TEMPVD"][:, 2] + 273.15)
         T_z = get_1d_interpolator(z, Tvd)
     else
         T_z = missing
@@ -359,9 +402,34 @@ function parse_state0_equil(model, datafile)
                 else
                     rv = missing
                 end
+                if haskey(props, "ZMFVD")
+                    ztab = props["ZMFVD"][ereg]
+                    N = size(ztab, 2) - 1
+                    Sz = SVector{N, Float64}
+                    ztab_z = Float64[]
+                    ztab_static = Sz[]
+                    for row in size(ztab, 1)
+                        mf_i = Sz(ztab[row, 2:end])
+                        z_i = ztab[row, 1]
+                        push!(ztab_z, z_i)
+                        push!(ztab_static, mf_i)
+                        if row == 1
+                            push!(ztab_z, z_i-1.0)
+                            push!(ztab_static, mf_i)
+                        end
+                        if row == size(ztab, 1)
+                            push!(ztab_z, z_i+1.0)
+                            push!(ztab_static, mf_i)
+                        end
+                    end
+                    composition = get_1d_interpolator(ztab_z, ztab_static)
+                else
+                    composition = missing
+                end
 
                 subinit = equilibriate_state(
                         model, contacts, datum_depth, datum_pressure,
+                        composition = composition,
                         cells = cells,
                         pvtnum = preg,
                         contacts_pc = contacts_pc,

src/input_simulation/data_input.jl

@@ -426,8 +426,8 @@ function filter_inactive_completions!(completions_vector, g)
                 if !haskey(tuple_active, k)
                     ix = cell_index(g, tupl, throw = false)
                     if isnothing(ix)
-                        jutul_message("Removing well",
-                            "Well $well completion $tupl was declared using COMPDAT, but that cell is not active in model. Skipped.",
+                        jutul_message("$well completion",
+                            "Removed COMPDAT as $tupl is not active in processed mesh.",
                             color = :yellow
                         )
                         tuple_active[k] = false
@@ -824,7 +824,7 @@ function parse_physics_types(datafile; pvt_region = missing)
         push!(phases, VaporPhase())
         push!(rhoS, rhoGS)
 
-        cnames = props["CNAMES"]
+        cnames = copy(props["CNAMES"])
         acf = props["ACF"]
         mw = props["MW"]
         p_c = props["PCRIT"]
@@ -837,6 +837,7 @@ function parse_physics_types(datafile; pvt_region = missing)
             A_ij = nothing
         end
         mp = MolecularProperty.(mw, p_c, T_c, V_c, acf)
+        @assert length(cnames) == length(mp)
         mixture = MultiComponentMixture(mp, A_ij = A_ij, names = cnames)
         if haskey(props, "EOS")
             eos_str = uppercase(props["EOS"])
@@ -853,7 +854,12 @@ function parse_physics_types(datafile; pvt_region = missing)
         else
             eos_type = PengRobinson()
         end
-        eos = GenericCubicEOS(mixture, eos_type)
+        if haskey(props, "SSHIFT")
+            vshift = Tuple(props["SSHIFT"])
+        else
+            vshift = nothing
+        end
+        eos = GenericCubicEOS(mixture, eos_type, volume_shift = vshift)
         sys = MultiPhaseCompositionalSystemLV(eos, phases, reference_densities = rhoS)
     else
         if has_oil