From dda843777f1bd9cd798fbd8fda395885c226d359 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Tue, 11 Jun 2024 21:18:13 +0200
Subject: [PATCH 01/36] Cap endpoint for Rs/Rv max table

---
 src/deck_types.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/deck_types.jl b/src/deck_types.jl
index 962644ae..82a20bca 100644
--- a/src/deck_types.jl
+++ b/src/deck_types.jl
@@ -462,7 +462,7 @@ function saturated_table(p, r)
         p = vcat([-1.0, 0.0], p)
         r = vcat([0.0, 0.0], r)
     end
-    return get_1d_interpolator(p, r, cap_end = false)
+    return get_1d_interpolator(p, r, cap_end = true)
 end
 
 pvt_table_vectors(pvt::PVTGTable) = (pvt.rv, pvt.pressure, pvt.sat_rv, pvt.pos)

From 5b34bf0fff4e007fb9bd731cdb39d80be2679a83 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Tue, 11 Jun 2024 22:15:21 +0200
Subject: [PATCH 02/36] Relax min pressure - better for history matching

---
 src/input_simulation/data_input.jl | 6 +++---
 src/multiphase.jl                  | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl
index 081f0b58..eb9c3041 100644
--- a/src/input_simulation/data_input.jl
+++ b/src/input_simulation/data_input.jl
@@ -425,12 +425,12 @@ function parse_schedule(domain, runspec, props, schedule, sys; simple_well = tru
             k_is_simple_well = simple_well
         end
         W, wc_base, WI_base, open = parse_well_from_compdat(domain, k, v, wspec, msdata_k, compord; simple_well = k_is_simple_well)
+        wpi_mul = ones(n_wi)
         for (i, c) in enumerate(completions)
             compdat = c[k]
             well_is_shut = controls[i][k] isa DisabledControl
             n_wi = length(WI_base)
             wi_mul = zeros(n_wi)
-            wpi_mul = ones(n_wi)
             if !well_is_shut
                 wc, WI, open = compdat_to_connection_factors(domain, wspec, compdat, sort = false)
                 for (c, wi, is_open) in zip(wc, WI, open)
@@ -1498,9 +1498,9 @@ function producer_control(sys, flag, ctrl, orat, wrat, grat, lrat, bhp; is_hist
             self_symbol = translate_target_to_symbol(t, shortname = true)
             # Put pressure slightly above 1 atm to avoid hard limit.
             if self_symbol == :resv_history
-                lims = (; :bhp => 1.001*si_unit(:atm))
+                lims = (; :bhp => 1.01*si_unit(:atm))
             else
-                lims = (; :bhp => 1.001*si_unit(:atm), self_symbol => self_val)
+                lims = (; :bhp => 1.01*si_unit(:atm), self_symbol => self_val)
             end
         else
             lims = producer_limits(bhp = bhp, orat = orat, wrat = wrat, grat = grat, lrat = lrat)
diff --git a/src/multiphase.jl b/src/multiphase.jl
index fb1bee12..f944709b 100644
--- a/src/multiphase.jl
+++ b/src/multiphase.jl
@@ -57,7 +57,7 @@ phase_name(::VaporPhase) = "Vapor"
 ## Main implementation
 # Primary variable logic
 
-const DEFAULT_MINIMUM_PRESSURE = 101325.0
+const DEFAULT_MINIMUM_PRESSURE = 10000.0 # 0.1 MPa
 
 struct Pressure <: ScalarVariable
     max_abs::Union{Float64, Nothing}

From 8ecddc80bddaf12c95a724bcd4286567ffedb46f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Tue, 11 Jun 2024 22:33:20 +0200
Subject: [PATCH 03/36] Put in place some multiplier logic

---
 src/input_simulation/data_input.jl | 25 ++++++++++++++++---------
 1 file changed, 16 insertions(+), 9 deletions(-)

diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl
index eb9c3041..7f10cada 100644
--- a/src/input_simulation/data_input.jl
+++ b/src/input_simulation/data_input.jl
@@ -182,8 +182,8 @@ function setup_case_from_parsed_data(datafile;
     return JutulCase(model, dt, forces, state0 = state0, parameters = parameters, input_data = datafile)
 end
 
-function parse_well_from_compdat(domain, wname, cdat, wspecs, msdata, compord; simple_well = isnothing(msdata))
-    wc, WI, open = compdat_to_connection_factors(domain, wspecs, cdat, sort = true, order = compord)
+function parse_well_from_compdat(domain, wname, cdat, wspecs, msdata, compord, step; simple_well = isnothing(msdata))
+    wc, WI, open, = compdat_to_connection_factors(domain, wspecs, cdat, step, sort = true, order = compord)
     ref_depth = wspecs.ref_depth
     if isnan(ref_depth)
         ref_depth = nothing
@@ -347,7 +347,7 @@ function map_compdat_to_multisegment_segments(compsegs, branches, tubing_lengths
     return perforation_cells
 end
 
-function compdat_to_connection_factors(domain, wspec, v; sort = true, order = "TRACK")
+function compdat_to_connection_factors(domain, wspec, v, step; sort = true, order = "TRACK")
     G = physical_representation(domain)
     K = domain[:permeability]
 
@@ -361,6 +361,7 @@ function compdat_to_connection_factors(domain, wspec, v; sort = true, order = "T
     WI = getf(:WI)
     skin = getf(:skin)
     dir = getf(:dir)
+    fresh = map(x -> v[x].ctrl == step, ij_ix)
 
     for i in eachindex(WI)
         W_i = WI[i]
@@ -394,8 +395,9 @@ function compdat_to_connection_factors(domain, wspec, v; sort = true, order = "T
         wc = wc[ix]
         WI = WI[ix]
         open = open[ix]
+        fresh = fresh[ix]
     end
-    return (wc, WI, open)
+    return (wc, WI, open, fresh)
 end
 
 function parse_schedule(domain, runspec, props, schedule, sys; simple_well = true)
@@ -424,22 +426,27 @@ function parse_schedule(domain, runspec, props, schedule, sys; simple_well = tru
             msdata_k = nothing
             k_is_simple_well = simple_well
         end
-        W, wc_base, WI_base, open = parse_well_from_compdat(domain, k, v, wspec, msdata_k, compord; simple_well = k_is_simple_well)
+        W, wc_base, WI_base, open = parse_well_from_compdat(domain, k, v, wspec, msdata_k, compord, length(completions); simple_well = k_is_simple_well)
+        n_wi = length(WI_base)
         wpi_mul = ones(n_wi)
         for (i, c) in enumerate(completions)
             compdat = c[k]
             well_is_shut = controls[i][k] isa DisabledControl
-            n_wi = length(WI_base)
             wi_mul = zeros(n_wi)
             if !well_is_shut
-                wc, WI, open = compdat_to_connection_factors(domain, wspec, compdat, sort = false)
-                for (c, wi, is_open) in zip(wc, WI, open)
+                wc, WI, open, fresh = compdat_to_connection_factors(domain, wspec, compdat, i, sort = false)
+                for (c, wi, is_open, is_fresh) in zip(wc, WI, open, fresh)
                     compl_idx = findfirst(isequal(c), wc_base)
                     if isnothing(compl_idx)
                         # This perforation is missing, leave at zero.
                         continue
                     end
-                    wi_mul[compl_idx] = wi*is_open/WI_base[compl_idx]
+                    if is_fresh
+                        # TODO: Make sure that WPIMULT is actually handled, this
+                        # is half-finished and does not have impact.
+                        wpi_mul[compl_idx] = 1.0
+                    end
+                    wi_mul[compl_idx] = wpi_mul[compl_idx]*wi*is_open/WI_base[compl_idx]
                 end
             end
             if all(isapprox(1.0), wi_mul)

From f42cdbadd5379656f8d0c49023cd4da35e19e48a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Wed, 12 Jun 2024 14:24:24 +0200
Subject: [PATCH 04/36] WPIMULT support

---
 src/input_simulation/data_input.jl | 49 +++++++++++++++++++++++++-----
 1 file changed, 41 insertions(+), 8 deletions(-)

diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl
index 7f10cada..5043509e 100644
--- a/src/input_simulation/data_input.jl
+++ b/src/input_simulation/data_input.jl
@@ -361,6 +361,7 @@ function compdat_to_connection_factors(domain, wspec, v, step; sort = true, orde
     WI = getf(:WI)
     skin = getf(:skin)
     dir = getf(:dir)
+    mul = getf(:mul)
     fresh = map(x -> v[x].ctrl == step, ij_ix)
 
     for i in eachindex(WI)
@@ -397,7 +398,7 @@ function compdat_to_connection_factors(domain, wspec, v, step; sort = true, orde
         open = open[ix]
         fresh = fresh[ix]
     end
-    return (wc, WI, open, fresh)
+    return (wc, WI, open, mul, fresh)
 end
 
 function parse_schedule(domain, runspec, props, schedule, sys; simple_well = true)
@@ -434,18 +435,23 @@ function parse_schedule(domain, runspec, props, schedule, sys; simple_well = tru
             well_is_shut = controls[i][k] isa DisabledControl
             wi_mul = zeros(n_wi)
             if !well_is_shut
-                wc, WI, open, fresh = compdat_to_connection_factors(domain, wspec, compdat, i, sort = false)
-                for (c, wi, is_open, is_fresh) in zip(wc, WI, open, fresh)
+                wc, WI, open, multipliers, fresh = compdat_to_connection_factors(domain, wspec, compdat, i, sort = false)
+                for (c, wi, is_open, is_fresh, mul) in zip(wc, WI, open, fresh, multipliers)
                     compl_idx = findfirst(isequal(c), wc_base)
                     if isnothing(compl_idx)
                         # This perforation is missing, leave at zero.
                         continue
                     end
-                    if is_fresh
-                        # TODO: Make sure that WPIMULT is actually handled, this
-                        # is half-finished and does not have impact.
-                        wpi_mul[compl_idx] = 1.0
+                    if !is_fresh
+                        # Perforation is not fresh and could have previously
+                        # gotten a WPIMULT applied. TODO: Some of the logic
+                        # around WPIMULT is not 100% clear. The current
+                        # implementation should be fine up to the stage where
+                        # many different keywords start to interact for the same
+                        # timestep.
+                        mul *= wpi_mul[compl_idx]
                     end
+                    wpi_mul[compl_idx] = mul
                     wi_mul[compl_idx] = wpi_mul[compl_idx]*wi*is_open/WI_base[compl_idx]
                 end
             end
@@ -1256,7 +1262,7 @@ function parse_control_steps(runspec, props, schedule, sys)
         push!(cstep, ctrl_ix)
     end
 
-    skip = ("WELLSTRE", "WINJGAS", "GINJGAS", "GRUPINJE", "WELLINJE", "WEFAC", "WTEMP")
+    skip = ("WELLSTRE", "WINJGAS", "GINJGAS", "GRUPINJE", "WELLINJE", "WEFAC", "WTEMP", "WPIMULT")
     bad_kw = Dict{String, Bool}()
     for (ctrl_ix, step) in enumerate(steps)
         found_time = false
@@ -1316,6 +1322,32 @@ function parse_control_steps(runspec, props, schedule, sys)
                     end
                     entry = compdat[wname]
                     for K in K1:K2
+                        perf_mult = 1.0
+                        current_completion_index = length(keys(entry)) + 1
+                        if haskey(step, "WPIMULT")
+                            for wpimult in step["WPIMULT"]
+                                wname_wpi, mul_i, I_wpi, J_wpi, K_wpi, start_wpi, end_wpi = wpimult
+                                if wname_wpi != wname
+                                    continue
+                                end
+                                if I != I_wpi && I_wpi > 0
+                                    continue
+                                end
+                                if J != J_wpi && J_wpi > 0
+                                    continue
+                                end
+                                if K != K_wpi && K_wpi > 0
+                                    continue
+                                end
+                                if start_wpi > current_completion_index && start_wpi > 0
+                                    continue
+                                end
+                                if end_wpi < current_completion_index && end_wpi > 0
+                                    continue
+                                end
+                                perf_mult *= mul_i
+                            end
+                        end
                         entry[(I, J, K)] = (
                             open = flag == "OPEN",
                             satnum = satnum,
@@ -1324,6 +1356,7 @@ function parse_control_steps(runspec, props, schedule, sys)
                             Kh = Kh,
                             skin = skin,
                             dir = dir,
+                            mul = perf_mult,
                             ctrl = ctrl_ix)
                     end
                 end

From 743f05fc1ff77dee983e2b6454769e89393c5144 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Wed, 12 Jun 2024 14:34:52 +0200
Subject: [PATCH 05/36] Clean up region matching for trans mult

---
 src/input_simulation/data_input.jl | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl
index 5043509e..6b3cbbcf 100644
--- a/src/input_simulation/data_input.jl
+++ b/src/input_simulation/data_input.jl
@@ -830,12 +830,15 @@ function parse_reservoir(data_file; zcorn_depths = true)
                 if do_apply
                     m = regt[3]
                     region_type = only(lowercase(regt[6]))
-                    if region_type == 'm' && pair_matchex(pairt, multnum_pair)
-                        tranmult[fno] *= m
-                    elseif region_type == 'o' && pair_matchex(pairt, opernum_pair)
-                        tranmult[fno] *= m
-                    elseif pair_matchex(pairt, fluxnum_pair)
-                        @assert region_type == 'f'
+                    if region_type == 'm'
+                        pair_to_match = multnum_pair
+                    elseif region_type == 'o'
+                        pair_to_match = opernum_pair
+                    else
+                        @assert region_type == 'f' "Region type was expected to be m, o or f, was $region_type"
+                        pair_to_match = fluxnum_pair
+                    end
+                    if pair_matchex(pairt, pair_to_match)
                         tranmult[fno] *= m
                     end
                 end

From 568ea7a6798d61f30cdbd7e599f80d0453cd0018 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Thu, 13 Jun 2024 13:39:40 +0200
Subject: [PATCH 06/36] Fix right preconditioning, adjust default linear
 tolerances

---
 src/cpr.jl      | 10 ++++++----
 src/linsolve.jl |  4 ++--
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/cpr.jl b/src/cpr.jl
index 69e0cf1e..5cf98c5f 100644
--- a/src/cpr.jl
+++ b/src/cpr.jl
@@ -298,14 +298,16 @@ function operator_nrows(cpr::CPRPreconditioner)
 end
 
 using Krylov
-function apply!(x, cpr::CPRPreconditioner, r, arg...)
+function apply!(x, cpr::CPRPreconditioner, r0, arg...)
     cpr_s = cpr.storage
-    buf = cpr_s.r_ps
-    # x = cpr_s.x_ps
+    # Get buffers and set working values
+    r = cpr_s.r_ps
+    @. r  = r0
+    buf = cpr_s.x_ps
     A_ps = cpr_s.A_ps
     smoother = cpr.system_precond
     bz = cpr_s.block_size
-    # Zero out buffer, just in case
+    # Zero out buffer, just in case (assumed by some solvers)
     @. x = 0.0
     # presmooth
     @tic "cpr smoother" apply_cpr_smoother!(x, r, buf, smoother, A_ps, cpr.npre)
diff --git a/src/linsolve.jl b/src/linsolve.jl
index 75ec5aab..10598078 100644
--- a/src/linsolve.jl
+++ b/src/linsolve.jl
@@ -43,7 +43,7 @@ function reservoir_linsolve(model,
     if !Jutul.is_cell_major(matrix_layout(model.context))
         return nothing
     end
-    default_tol = 0.005
+    default_tol = 0.01
     max_it = 200
     if precond == :cpr
         if isnothing(cpr_type)
@@ -68,7 +68,7 @@ function reservoir_linsolve(model,
             update_interval_partial = update_interval_partial,
             mode = mode
         )
-        default_tol = 1e-3
+        default_tol = 0.005
         max_it = 50
     elseif precond == :ilu0
         prec = ILUZeroPreconditioner()

From 1a570e85fc1190dfeed15ebf9b0d804a2b1613ce Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Thu, 13 Jun 2024 14:40:55 +0200
Subject: [PATCH 07/36] ENDSCALE logic fix

---
 src/input_simulation/mrst_input.jl | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/src/input_simulation/mrst_input.jl b/src/input_simulation/mrst_input.jl
index 3f5e5c79..8b568f4d 100644
--- a/src/input_simulation/mrst_input.jl
+++ b/src/input_simulation/mrst_input.jl
@@ -369,12 +369,17 @@ function deck_relperm(props; oil, water, gas, satnum = nothing)
         # Early return for single-phase.
         return BrooksCoreyRelativePermeabilities(1)
     end
-    if haskey(props, "SCALECRS")
-        scalecrs = props["SCALECRS"]
-        if scalecrs isa String
-            scalecrs = [scalecrs]
+    if haskey(props, "ENDSCALE")
+        if haskey(props, "SCALECRS")
+            scalecrs = props["SCALECRS"]
+            if scalecrs isa String
+                scalecrs = [scalecrs]
+            end
+            two_point_scaling = length(scalecrs) == 0 || lowercase(first(scalecrs)) == "no"
+        else
+            two_point_scaling = true
         end
-        if length(scalecrs) == 0 || lowercase(first(scalecrs)) == "no"
+        if two_point_scaling
             Jutul.jutul_message("Rel. Perm. Scaling", "Two-point scaling active.")
             scaling = TwoPointKrScale
         else

From 1505f3a37fbbbab8c4e561d7aa0b67f00f347be5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Fri, 14 Jun 2024 11:19:49 +0200
Subject: [PATCH 08/36] Don't extrapolate rs/rv vs depth tables

---
 src/init/init.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/init/init.jl b/src/init/init.jl
index b11c1ebe..0d4e59ca 100644
--- a/src/init/init.jl
+++ b/src/init/init.jl
@@ -426,7 +426,7 @@ function parse_state0_equil(model, datafile)
                         pb = vec(pbvd[:, 2])
                         Rs = sys.rs_max[preg].(pb)
                     end
-                    rs = Jutul.LinearInterpolant(z, Rs_scale.*Rs)
+                    rs = get_1d_interpolator(z, Rs_scale.*Rs)
                 else
                     rs = missing
                 end
@@ -438,7 +438,7 @@ function parse_state0_equil(model, datafile)
                         rvvd = sol["RVVD"][ereg]
                         z = rvvd[:, 1]
                         Rv = rvvd[:, 2]
-                        rv = Jutul.LinearInterpolant(z, Rv_scale.*Rv)
+                        rv = get_1d_interpolator(z, Rv_scale.*Rv)
                     else
                         # TODO: Should maybe be the pressure at GOC not datum
                         # depth, but that isn't computed at this stage.

From cc4626cfe0404d72330f52489993f552e893d508 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Fri, 14 Jun 2024 11:41:44 +0200
Subject: [PATCH 09/36] Refactor relperm and reduce code duplication

---
 src/variables/relperm.jl | 41 +++++++++++++++++++++++++++-------------
 1 file changed, 28 insertions(+), 13 deletions(-)

diff --git a/src/variables/relperm.jl b/src/variables/relperm.jl
index 840eceec..1cc7ed9f 100644
--- a/src/variables/relperm.jl
+++ b/src/variables/relperm.jl
@@ -401,22 +401,11 @@ end
 end
 
 @jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeabilities{NoKrScale, :wog}, model, Saturations, ConnateWater, ix)
-    s = Saturations
-    phases = phase_indices(model.system)
-    for c in ix
-        @inbounds update_three_phase_relperm!(kr, relperm, phases, s, c, ConnateWater[c], nothing)
-    end
-    return kr
+    return evaluate_relative_permeability_no_scaling!(kr, relperm, model, Saturations, ix, ConnateWater)
 end
 
 @jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeabilities{NoKrScale, :wo}, model, Saturations, ix)
-    s = Saturations
-    regions = relperm.regions
-    indices = phase_indices(model.system)
-    for c in ix
-        @inbounds two_phase_relperm!(kr, s, regions, relperm.krw, relperm.krow, indices, c)
-    end
-    return kr
+    return evaluate_relative_permeability_no_scaling!(kr, relperm, model, Saturations, ix)
 end
 
 @jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeabilities{NoKrScale, :og}, model, Saturations, ix)
@@ -439,6 +428,32 @@ end
     return kr
 end
 
+function evaluate_relative_permeability_no_scaling!(kr, relperm::ReservoirRelativePermeabilities{<:Any, ph}, model, Saturations, ix, ConnateWater = missing) where ph
+    s = Saturations
+    regions = relperm.regions
+    phases = phase_indices(model.system)
+    if ph == :wog
+        @assert !ismissing(ConnateWater)
+        for c in ix
+            @inbounds update_three_phase_relperm!(kr, relperm, phases, s, c, ConnateWater[c], nothing)
+        end
+    elseif ph == :wo
+        for c in ix
+            @inbounds two_phase_relperm!(kr, s, regions, relperm.krw, relperm.krow, phases, c)
+        end
+    elseif ph == :og
+        for c in ix
+            @inbounds two_phase_relperm!(kr, s, regions, relperm.krg, relperm.krog, reverse(phases), c)
+        end
+    else
+        @assert ph == :wg
+        for c in ix
+            @inbounds two_phase_relperm!(kr, s, regions, relperm.krw, relperm.krg, phases, c)
+        end
+    end
+    return kr
+end
+
 Base.@propagate_inbounds @inline function three_phase_oil_relperm(Krow, Krog, swcon, sg, sw)
     swc = min(swcon, value(sw) - 1e-5)
     d  = (sg + sw - swc)

From 242e5ccd7f1bad1a8c7cbcc97200be718c6cad7b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Fri, 14 Jun 2024 11:44:39 +0200
Subject: [PATCH 10/36] Update initialization to use new rel perm

---
 src/init/init.jl         | 20 +++++++-------------
 src/variables/relperm.jl | 16 ++--------------
 2 files changed, 9 insertions(+), 27 deletions(-)

diff --git a/src/init/init.jl b/src/init/init.jl
index 0d4e59ca..7c85b2fa 100644
--- a/src/init/init.jl
+++ b/src/init/init.jl
@@ -163,7 +163,6 @@ function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pr
                 end
                 s0 = s[:, i]
                 sw_i = sw[i]
-                # sw_i = max(sw[i], s[1, i])
                 s[1, i] = sw_i
                 s_rem = 0.0
                 for ph in 2:nph
@@ -188,21 +187,16 @@ function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pr
         s_eval = zeros(nph, nc_total)
         s_eval[:, cells] .= s
         phases = get_phases(sys)
-        if scaling_type(relperm) == NoKrScale
-            if length(phases) == 3 && AqueousPhase() in phases
-                swcon = zeros(nc_total)
-                if !ismissing(s_min)
-                    swcon[cells] .= s_min[1]
-                end
-                JutulDarcy.update_kr!(kr, relperm, model, s_eval, swcon, cells)
-            else
-                JutulDarcy.update_kr!(kr, relperm, model, s_eval, cells)
+        if length(phases) == 3 && AqueousPhase() in phases
+            swcon = zeros(nc_total)
+            if !ismissing(s_min)
+                swcon[cells] .= s_min[1]
             end
-            kr = kr[:, cells]
+            evaluate_relative_permeability_no_scaling!(kr, relperm, model, s_eval, cells, swcon)
         else
-            # TODO: Improve this.
-            kr = s
+            JutulDarcy.update_kr!(kr, relperm, model, s_eval, cells)
         end
+        kr = kr[:, cells]
         init[:Saturations] = s
         init[:Pressure] = init_reference_pressure(pressures, contacts, kr, pc, 2)
     else
diff --git a/src/variables/relperm.jl b/src/variables/relperm.jl
index 1cc7ed9f..ad1a5e6f 100644
--- a/src/variables/relperm.jl
+++ b/src/variables/relperm.jl
@@ -409,23 +409,11 @@ end
 end
 
 @jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeabilities{NoKrScale, :og}, model, Saturations, ix)
-    s = Saturations
-    regions = relperm.regions
-    indices = phase_indices(model.system)
-    for c in ix
-        @inbounds two_phase_relperm!(kr, s, regions, relperm.krg, relperm.krog, reverse(indices), c)
-    end
-    return kr
+    return evaluate_relative_permeability_no_scaling!(kr, relperm, model, Saturations, ix)
 end
 
 @jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeabilities{NoKrScale, :wg}, model, Saturations, ix)
-    s = Saturations
-    regions = relperm.regions
-    indices = phase_indices(model.system)
-    for c in ix
-        @inbounds two_phase_relperm!(kr, s, regions, relperm.krw, relperm.krg, indices, c)
-    end
-    return kr
+    return evaluate_relative_permeability_no_scaling!(kr, relperm, model, Saturations, ix)
 end
 
 function evaluate_relative_permeability_no_scaling!(kr, relperm::ReservoirRelativePermeabilities{<:Any, ph}, model, Saturations, ix, ConnateWater = missing) where ph

From 33d3d4e6915173c0256e86a947c0a2b24d7e151b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Sun, 16 Jun 2024 18:53:24 +0200
Subject: [PATCH 11/36] A few fixes to compositional injectors gas stream

---
 src/input_simulation/data_input.jl | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl
index 6b3cbbcf..5d2e6538 100644
--- a/src/input_simulation/data_input.jl
+++ b/src/input_simulation/data_input.jl
@@ -1102,7 +1102,7 @@ function parse_physics_types(datafile; pvt_region = missing)
         push!(phases, VaporPhase())
         push!(rhoS, rhoGS)
         if length(props["MW"]) > 1
-            jutul_message("EOSNUM:", "$(length(props["MW"])) regions active. Only one region supported. Taking the first set of values for all EOS properties.", color = :yellow)
+            jutul_message("EOSNUM", "$(length(props["MW"])) regions active. Only one region supported. Taking the first set of values for all EOS properties.", color = :yellow)
         end
 
         cnames = copy(props["CNAMES"])
@@ -1267,9 +1267,10 @@ function parse_control_steps(runspec, props, schedule, sys)
 
     skip = ("WELLSTRE", "WINJGAS", "GINJGAS", "GRUPINJE", "WELLINJE", "WEFAC", "WTEMP", "WPIMULT")
     bad_kw = Dict{String, Bool}()
+    streams = setup_well_streams()
     for (ctrl_ix, step) in enumerate(steps)
         found_time = false
-        streams = parse_well_streams_for_step(step, props)
+        streams = parse_well_streams_for_step!(streams, step, props)
         if haskey(step, "WEFAC")
             for wk in step["WEFAC"]
                 well_factor[wk[1]] = wk[2]
@@ -1403,7 +1404,13 @@ function parse_control_steps(runspec, props, schedule, sys)
     return (dt = tstep, control_step = cstep, controls = all_controls, completions = all_compdat, multisegment = mswell_kw, limits = all_limits)
 end
 
-function parse_well_streams_for_step(step, props)
+function setup_well_streams()
+    return (streams = Dict{String, Any}(), wells = Dict{String, String}())
+end
+
+function parse_well_streams_for_step!(streams, step, props)
+    well_streams = streams.wells
+    streams = streams.streams
     if haskey(props, "STCOND")
         std = props["STCOND"]
         T = convert_to_si(std[1], :Celsius)
@@ -1414,8 +1421,6 @@ function parse_well_streams_for_step(step, props)
         p = convert_to_si(1.0, :atm)
     end
 
-    streams = Dict{String, Any}()
-    well_streams = Dict{String, String}()
     if haskey(step, "WELLSTRE")
         for stream in step["WELLSTRE"]
             mix = Float64.(stream[2:end])
@@ -1672,7 +1677,7 @@ function select_injector_mixture_spec(sys::CompositionalSystem, name, streams, t
             z, props
         )
         z_mass /= sum(z_mass)
-        for i in 1:ncomp
+        for i in eachindex(z_mass)
             mix[i+offset] = z_mass[i]
         end
         @assert sum(mix) ≈ 1.0 "Sum of mixture was $(sum(mix)) != 1 for mole mixture $(z) as mass $z_mass"

From 121ce546e29c1708f376ca7d787282d084759958 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Sun, 16 Jun 2024 19:00:12 +0200
Subject: [PATCH 12/36] Robustness for compositional inputs

---
 src/init/init.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/init/init.jl b/src/init/init.jl
index 7c85b2fa..d3117d19 100644
--- a/src/init/init.jl
+++ b/src/init/init.jl
@@ -291,9 +291,9 @@ function parse_state0_equil(model, datafile)
     for ereg in 1:nequil
         if haskey(sol, "RTEMP") || haskey(props, "RTEMP")
             if haskey(props, "RTEMP")
-                rtmp = only(props["RTEMP"])
+                rtmp = props["RTEMP"][1]
             else
-                rtmp = only(sol["RTEMP"])
+                rtmp = sol["RTEMP"][1]
             end
             Ti = convert_to_si(rtmp, :Celsius)
             T_z = z -> Ti
@@ -450,6 +450,7 @@ function parse_state0_equil(model, datafile)
                     ztab_static = Sz[]
                     for row in size(ztab, 1)
                         mf_i = Sz(ztab[row, 2:end])
+                        mf_i = mf_i./sum(mf_i)
                         z_i = ztab[row, 1]
                         push!(ztab_z, z_i)
                         push!(ztab_static, mf_i)

From c164866f4f8bdf9327379f96293313e14a0996f0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Sun, 16 Jun 2024 20:01:25 +0200
Subject: [PATCH 13/36] Split shutting of wells into two options

One for inj, one for prod
---
 src/facility/types.jl      | 9 ++++++---
 src/facility/wellgroups.jl | 5 ++---
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/facility/types.jl b/src/facility/types.jl
index 42ac7fb0..c5eb6e87 100644
--- a/src/facility/types.jl
+++ b/src/facility/types.jl
@@ -7,7 +7,10 @@ abstract type SurfaceFacilityDomain <: JutulDomain end
 abstract type WellControllerDomain <: SurfaceFacilityDomain end
 mutable struct WellGroup <: WellControllerDomain
     const well_symbols::Vector{Symbol} # Controlled wells
-    can_shut_wells::Bool               # Can temporarily shut wells that try to reach zero rate multiple solves in a row
+    "Can temporarily shut producers that try to reach zero rate multiple solves in a row"
+    can_shut_producers::Bool
+    "Can temporarily shut injectors that try to reach zero rate multiple solves in a row"
+    can_shut_injectors::Bool
 end
 
 """
@@ -15,8 +18,8 @@ end
 
 Create a well group that can control the given set of wells.
 """
-function WellGroup(wells::Vector{Symbol}; can_shut_wells = true)
-    return WellGroup(wells, can_shut_wells)
+function WellGroup(wells::Vector{Symbol}; can_shut_wells = true, can_shut_injectors = can_shut_wells, can_shut_producers = can_shut_wells)
+    return WellGroup(wells, can_shut_producers, can_shut_injectors)
 end
 
 const WellGroupModel = SimulationModel{WellGroup, <:Any, <:Any, <:Any}
diff --git a/src/facility/wellgroups.jl b/src/facility/wellgroups.jl
index ee12cef0..ce1f8622 100644
--- a/src/facility/wellgroups.jl
+++ b/src/facility/wellgroups.jl
@@ -16,7 +16,6 @@ function Jutul.update_primary_variable!(state, massrate::TotalSurfaceMassRate, s
     v = state[state_symbol]
     symbols = model.domain.well_symbols
     cfg = state.WellGroupConfiguration
-    can_shut_wells = model.domain.can_shut_wells
     # Injectors can only have strictly positive injection rates,
     # producers can only have strictly negative and disabled controls give zero rate.
     rel_max = Jutul.absolute_increment_limit(massrate)
@@ -26,8 +25,8 @@ function Jutul.update_primary_variable!(state, massrate::TotalSurfaceMassRate, s
     end
     function do_update!(wcfg, s, v, dx, ctrl::InjectorControl)
         limit_rate = MIN_INITIAL_WELL_RATE
-        if v <= limit_rate && v + dx < limit_rate && can_shut_wells
             jutul_message("$s", "Approaching zero rate, disabling injector for current step", color = :red)
+        if v <= limit_rate && v + dx < limit_rate && model.domain.can_shut_injectors
             wcfg.operating_controls[s] = DisabledControl()
             next = Jutul.update_value(v, -value(v))
         else
@@ -38,7 +37,7 @@ function Jutul.update_primary_variable!(state, massrate::TotalSurfaceMassRate, s
     function do_update!(wcfg, s, v, dx, ctrl::ProducerControl)
         # A significant negative rate is the valid producer control
         limit_rate = -MIN_INITIAL_WELL_RATE
-        if v >= limit_rate && v + dx > limit_rate && can_shut_wells
+        if v >= limit_rate && v + dx > limit_rate && model.domain.can_shut_producers
             jutul_message("$s", "Approaching zero rate, disabling producer for current step", color = :red)
             wcfg.operating_controls[s] = DisabledControl()
             next = Jutul.update_value(v, -value(v))

From 6e33658f7179010e5129817164858923880ed042 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Sun, 16 Jun 2024 21:36:52 +0200
Subject: [PATCH 14/36] Update mrst_input.jl

---
 src/input_simulation/mrst_input.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/input_simulation/mrst_input.jl b/src/input_simulation/mrst_input.jl
index 8b568f4d..3ada0f7a 100644
--- a/src/input_simulation/mrst_input.jl
+++ b/src/input_simulation/mrst_input.jl
@@ -1362,7 +1362,7 @@ function mrst_well_ctrl(model, wdata, is_comp, rhoS)
                 ct = comp_i
             end
             if haskey(wdata, "rhoS") && length(wdata["rhoS"]) > 0
-                rhoSw = tuple(wdata["rhoS"]...)
+                rhoSw = tuple(wdata["rhoS"][1:nph]...)
             else
                 rhoSw = rhoS
             end

From 188c31512bd19104039ffa971f92e9aefd411f91 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Mon, 17 Jun 2024 10:52:55 +0200
Subject: [PATCH 15/36] Printing fix

---
 src/facility/wellgroups.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/facility/wellgroups.jl b/src/facility/wellgroups.jl
index ce1f8622..f9467f0f 100644
--- a/src/facility/wellgroups.jl
+++ b/src/facility/wellgroups.jl
@@ -25,8 +25,8 @@ function Jutul.update_primary_variable!(state, massrate::TotalSurfaceMassRate, s
     end
     function do_update!(wcfg, s, v, dx, ctrl::InjectorControl)
         limit_rate = MIN_INITIAL_WELL_RATE
-            jutul_message("$s", "Approaching zero rate, disabling injector for current step", color = :red)
         if v <= limit_rate && v + dx < limit_rate && model.domain.can_shut_injectors
+            jutul_message("$s", "Approaching zero rate, disabling injector for current step", color = :red)
             wcfg.operating_controls[s] = DisabledControl()
             next = Jutul.update_value(v, -value(v))
         else

From 5a8669754b4e09cd3c839facbab7ceb694cedeae Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Mon, 17 Jun 2024 21:21:34 +0200
Subject: [PATCH 16/36] Fix bug in scaling setup

---
 src/input_simulation/data_input.jl |  2 +-
 src/input_simulation/mrst_input.jl | 16 ++++++++--------
 2 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl
index 5d2e6538..35011ab0 100644
--- a/src/input_simulation/data_input.jl
+++ b/src/input_simulation/data_input.jl
@@ -160,7 +160,7 @@ function setup_case_from_parsed_data(datafile;
                 set_thermal_deck_specialization!(submodel, props, domain[:pvtnum], oil, water, gas)
             end
             if k == :Reservoir
-                set_deck_specialization!(submodel, props, domain[:satnum], oil, water, gas)
+                set_deck_specialization!(submodel, rs, props, domain[:satnum], oil, water, gas)
             end
         end
     end
diff --git a/src/input_simulation/mrst_input.jl b/src/input_simulation/mrst_input.jl
index 3ada0f7a..ed8a244e 100644
--- a/src/input_simulation/mrst_input.jl
+++ b/src/input_simulation/mrst_input.jl
@@ -364,12 +364,12 @@ function deck_pvt_gas(props; scaling = missing)
     return pvt
 end
 
-function deck_relperm(props; oil, water, gas, satnum = nothing)
+function deck_relperm(runspec, props; oil, water, gas, satnum = nothing)
     if (water + oil + gas) == 1
         # Early return for single-phase.
         return BrooksCoreyRelativePermeabilities(1)
     end
-    if haskey(props, "ENDSCALE")
+    if haskey(runspec, "ENDSCALE")
         if haskey(props, "SCALECRS")
             scalecrs = props["SCALECRS"]
             if scalecrs isa String
@@ -652,7 +652,7 @@ function model_from_mat_deck(G, data_domain, mrst_data, res_context)
         if length(unique(T)) == 1
             T = T[1]
         end
-        set_deck_specialization!(model, props, satnum, has_oil, has_wat, has_gas)
+        set_deck_specialization!(model, runspec, props, satnum, has_oil, has_wat, has_gas)
         param = setup_parameters(model, Temperature = T)
     else
         if has_wat
@@ -685,7 +685,7 @@ function model_from_mat_deck(G, data_domain, mrst_data, res_context)
         end
         mu = DeckPhaseViscosities(pvt)
         set_secondary_variables!(model, PhaseViscosities = mu, PhaseMassDensities = rho)
-        set_deck_specialization!(model, props, satnum, has_oil, has_wat, has_gas)
+        set_deck_specialization!(model, runspec, props, satnum, has_oil, has_wat, has_gas)
         param = setup_parameters(model)
     end
 
@@ -711,12 +711,12 @@ function model_from_mat_deck(G, data_domain, mrst_data, res_context)
     return (model, param)
 end
 
-function set_deck_specialization!(model, props, satnum, oil, water, gas)
+function set_deck_specialization!(model, runspec, props, satnum, oil, water, gas)
     sys = model.system
     svar = model.secondary_variables
     param = model.parameters
     if number_of_phases(sys) > 1
-        set_deck_relperm!(svar, param, sys, props; oil = oil, water = water, gas = gas, satnum = satnum)
+        set_deck_relperm!(svar, param, sys, runspec, props; oil = oil, water = water, gas = gas, satnum = satnum)
         set_deck_pc!(svar, param, sys, props; oil = oil, water = water, gas = gas, satnum = satnum)
     end
     set_deck_pvmult!(svar, param, sys, props, model.data_domain)
@@ -778,8 +778,8 @@ function set_deck_pc!(vars, param, sys, props; kwarg...)
     end
 end
 
-function set_deck_relperm!(vars, param, sys, props; kwarg...)
-    kr = deck_relperm(props; kwarg...)
+function set_deck_relperm!(vars, param, sys, runspec,  props; kwarg...)
+    kr = deck_relperm(runspec, props; kwarg...)
     vars[:RelativePermeabilities] = wrap_reservoir_variable(sys, kr, :flow)
     if scaling_type(kr) != NoKrScale
         ph = kr.phases

From 013644a7373d94104353fb92e16eea6ec97b70d3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Mon, 17 Jun 2024 21:21:55 +0200
Subject: [PATCH 17/36] Update init.jl

---
 src/init/init.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/init/init.jl b/src/init/init.jl
index d3117d19..b044dfc9 100644
--- a/src/init/init.jl
+++ b/src/init/init.jl
@@ -194,7 +194,7 @@ function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pr
             end
             evaluate_relative_permeability_no_scaling!(kr, relperm, model, s_eval, cells, swcon)
         else
-            JutulDarcy.update_kr!(kr, relperm, model, s_eval, cells)
+            evaluate_relative_permeability_no_scaling!(kr, relperm, model, s_eval, cells)
         end
         kr = kr[:, cells]
         init[:Saturations] = s

From 12a256d0b80b9ca270155ce61100fe770664c341 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Tue, 18 Jun 2024 11:12:24 +0200
Subject: [PATCH 18/36] compositional: Always compute mass fractions for absent
 phases

Better for output purposes (maybe also derivatives?)
---
 src/multicomponent/variables/others.jl | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/src/multicomponent/variables/others.jl b/src/multicomponent/variables/others.jl
index 29c9298e..54352e60 100644
--- a/src/multicomponent/variables/others.jl
+++ b/src/multicomponent/variables/others.jl
@@ -14,12 +14,9 @@ end
     phase = m.phase
     @inbounds for i in ix
         f = FlashResults[i]
-        if phase_is_present(phase, f.state)
-            # X_i = view(X, :, i)
-            r = phase_data(f, phase)
-            x_i = r.mole_fractions
-            update_mass_fractions!(X, x_i, i, molar_mass, n)
-        end
+        r = phase_data(f, phase)
+        x_i = r.mole_fractions
+        update_mass_fractions!(X, x_i, i, molar_mass, n)
     end
 end
 

From 0f69ffa3b4fc604b770dfe3e6b4b6679b5049c52 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Tue, 18 Jun 2024 13:33:57 +0200
Subject: [PATCH 19/36] Remove molar mass scaling from convergence criterion

I think this is now outdated
---
 src/multicomponent/multicomponent.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/multicomponent/multicomponent.jl b/src/multicomponent/multicomponent.jl
index 8cbe0b80..42d35074 100644
--- a/src/multicomponent/multicomponent.jl
+++ b/src/multicomponent/multicomponent.jl
@@ -178,7 +178,7 @@ function compositional_residual_scale(cell, dt, w, sl, liquid_density, sv, vapor
         denl = liquid_density[cell]
         denv = vapor_density[cell]
         total_density = denl * sl_scaled + denv * (1.0 - sl_scaled)
-        scale = dt * mean(w) * (scale_lv / (vol[cell] * max(total_density, 1e-8)))
+        scale = dt * (scale_lv / (vol[cell] * max(total_density, 1e-8)))
     end
     return scale
 end
@@ -188,7 +188,7 @@ function compositional_criterion(state, dt, active, r, nc, w, sl, liquid_density
     for (ix, i) in enumerate(active)
         scaling = compositional_residual_scale(i, dt, w, sl, liquid_density, sv, vapor_density, sw, water_density, vol)
         for c in 1:(nc-1)
-            val = scaling*abs(r[c, ix])/w[c]
+            val = scaling*abs(r[c, ix])
             if val > e[c]
                 e[c] = val
             end
@@ -206,7 +206,7 @@ function compositional_criterion(state, dt, active, r, nc, w, sl, liquid_density
     for (ix, i) in enumerate(active)
         scaling = compositional_residual_scale(i, dt, w, sl, liquid_density, sv, vapor_density, sw, water_density, vol)
         for c in 1:nc
-            e[c] = max(e[c], scaling*abs(r[c, ix])/w[c])
+            e[c] = max(e[c], scaling*abs(r[c, ix]))
         end
     end
     return e

From 9c4e1cc99997c2705929c8d604c6a1c273d2ef1f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Tue, 18 Jun 2024 13:34:12 +0200
Subject: [PATCH 20/36] Slight cleanup of cpr.jl

---
 src/cpr.jl | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/cpr.jl b/src/cpr.jl
index 5cf98c5f..a6b91fea 100644
--- a/src/cpr.jl
+++ b/src/cpr.jl
@@ -380,7 +380,8 @@ function update_weights!(cpr, cpr_storage::CPRStorage, model, res_storage, J, ps
     r = cpr_storage.w_rhs
     ncomp = size(w, 1)
     scaling = cpr.weight_scaling
-    if cpr.strategy == :true_impes
+    strategy = cpr.strategy
+    if strategy == :true_impes
         eq_s = res_storage.equations[:mass_conservation]
         if eq_s isa ConservationLawTPFAStorage
             acc = eq_s.accumulation.entries
@@ -390,15 +391,15 @@ function update_weights!(cpr, cpr_storage::CPRStorage, model, res_storage, J, ps
             ps = 1.0
         end
         true_impes!(w, acc, r, n, ncomp, ps, scaling)
-    elseif cpr.strategy == :analytical
+    elseif strategy == :analytical
         rstate = res_storage.state
         cpr_weights_no_partials!(w, model, rstate, r, n, ncomp, scaling)
-    elseif cpr.strategy == :quasi_impes
+    elseif strategy == :quasi_impes
         quasi_impes!(w, J, r, n, ncomp, scaling)
-    elseif cpr.strategy == :none
+    elseif strategy == :none
         # Do nothing. Already set to one.
     else
-        error("Unsupported strategy $(cpr.strategy)")
+        error("Unsupported strategy $(strategy)")
     end
     return w
 end

From 54447a93512919bea28c51e66e0abe30c051b9b6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=98ystein=20Klemetsdal?= <oystein.klemetsdal@gmail.com>
Date: Tue, 18 Jun 2024 21:19:15 +0200
Subject: [PATCH 21/36] setup_reservoir_model: Change default value for
 dT_max_abs from nothing to 50.0

---
 src/utils.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/utils.jl b/src/utils.jl
index dc8cb44d..c4c9d906 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -205,7 +205,7 @@ low values can lead to slow convergence.
 - `dT_max_rel=`: Maximum relative change in temperature (JutulDarcy uses Kelvin,
   so comments about changing limits near zero above does not apply to typical
   reservoir temperatures)
-- `dT_max_abs=`: Maximum absolute change in temperature (in °K/°C)
+- `dT_max_abs=50.0`: Maximum absolute change in temperature (in °K/°C)
 - `fast_flash=false`: Shorthand to enable `flash_reuse_guess` and
   `flash_stability_bypass`. These options can together speed up the time spent
   in flash solver for compositional models. Options are based on "Increasing the
@@ -241,7 +241,7 @@ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem;
         p_max = Inf,
         dr_max = Inf,
         dT_max_rel = nothing,
-        dT_max_abs = nothing,
+        dT_max_abs = 50.0,
         fast_flash = false,
         can_shut_wells = true,
         flash_reuse_guess = fast_flash,

From 6c20e7c8ee54c9f77f8a9549ad0088273f31dcab Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Thu, 20 Jun 2024 22:46:46 +0200
Subject: [PATCH 22/36] Refactor well downtime logic

---
 src/facility/controls.jl    |  9 +++++++--
 src/facility/cross_terms.jl | 21 ++++++++-------------
 src/facility/types.jl       |  3 +++
 src/utils.jl                |  6 ++----
 4 files changed, 20 insertions(+), 19 deletions(-)

diff --git a/src/facility/controls.jl b/src/facility/controls.jl
index a38edf90..2075826b 100644
--- a/src/facility/controls.jl
+++ b/src/facility/controls.jl
@@ -198,9 +198,14 @@ function check_limit(current_control, target_limit, target, is_lower::Bool, q_t,
 end
 
 
-function facility_surface_mass_rate_for_well(model::SimulationModel, wsym, fstate)
+function facility_surface_mass_rate_for_well(model::SimulationModel, wsym, fstate; effective::Bool = false)
     pos = get_well_position(model.domain, wsym)
-    return fstate.TotalSurfaceMassRate[pos]
+    q_t = fstate.TotalSurfaceMassRate[pos]
+    if effective
+        control = fstate.WellGroupConfiguration.operating_controls[wsym]
+        q_t = effective_surface_rate(q_t, control)
+    end
+    return q_t
 end
 
 bottom_hole_pressure(ws) = ws.Pressure[1]
diff --git a/src/facility/cross_terms.jl b/src/facility/cross_terms.jl
index b929883d..53409749 100644
--- a/src/facility/cross_terms.jl
+++ b/src/facility/cross_terms.jl
@@ -114,7 +114,7 @@ function Jutul.prepare_cross_term_in_entity!(i,
         limits = current_limits(cfg, well_symbol)
         if !isnothing(limits)
             rhoS, S = surface_density_and_volume_fractions(state_well)
-            q_t = facility_surface_mass_rate_for_well(facility, well_symbol, state_facility)
+            q_t = facility_surface_mass_rate_for_well(facility, well_symbol, state_facility, effective = true)
             apply_well_limit!(cfg, target, well, state_well, well_symbol, rhoS, S, value(q_t), limits)
         end
     end
@@ -137,7 +137,7 @@ function update_cross_term_in_entity!(out, i,
     ctrl = operating_control(cfg, well_symbol)
 
     target = ctrl.target
-    q_t = facility_surface_mass_rate_for_well(facility, well_symbol, state_facility)
+    q_t = facility_surface_mass_rate_for_well(facility, well_symbol, state_facility, effective = true)
     t, t_num = target_actual_pair(target, well, state_well, q_t, ctrl)
     t += 0*bottom_hole_pressure(state_well) + 0*q_t
     scale = target_scaling(target)
@@ -188,17 +188,12 @@ function update_cross_term_in_entity!(out, i,
 
     cfg = state_facility.WellGroupConfiguration
     ctrl = operating_control(cfg, well_symbol)
-    qT = state_facility.TotalSurfaceMassRate[pos] 
+    q_t = facility_surface_mass_rate_for_well(facility, well_symbol, state_facility, effective = true)
     # Hack for sparsity detection
-    qT += 0*bottom_hole_pressure(state_well)
-    if ctrl isa DisabledControl
-        factor = 1.0
-    else
-        factor = ctrl.factor
-    end
-    qT *= factor
+    q_t += 0*bottom_hole_pressure(state_well)
+
     if isa(ctrl, InjectorControl)
-        if value(qT) < 0
+        if value(q_t) < 0
             @warn "Injector $well_symbol is producing?"
         end
         mix = ctrl.injection_mixture
@@ -206,7 +201,7 @@ function update_cross_term_in_entity!(out, i,
         ncomp = number_of_components(flow_system(well.system))
         @assert nmix == ncomp "Injection composition length ($nmix) must match number of components ($ncomp)."
     else
-        if value(qT) > 0 && ctrl isa ProducerControl
+        if value(q_t) > 0 && ctrl isa ProducerControl
             @warn "Producer $well_symbol is injecting?"
         end
         if haskey(state_well, :MassFractions)
@@ -218,7 +213,7 @@ function update_cross_term_in_entity!(out, i,
         end
     end
     for i in eachindex(out)
-        @inbounds out[i] = -mix[i]*qT
+        @inbounds out[i] = -mix[i]*q_t
     end
 end
 
diff --git a/src/facility/types.jl b/src/facility/types.jl
index c5eb6e87..786d953f 100644
--- a/src/facility/types.jl
+++ b/src/facility/types.jl
@@ -366,6 +366,9 @@ function replace_target(f::ProducerControl, target)
     return ProducerControl(target, factor = f.factor)
 end
 
+effective_surface_rate(qts, ::DisabledControl) = qts
+effective_surface_rate(qts, c::Union{InjectorControl, ProducerControl}) = qts*c.factor
+
 mutable struct WellGroupConfiguration{T, O, L}
     const operating_controls::T # Currently operating control
     const requested_controls::O # The requested control (which may be different if limits are hit)
diff --git a/src/utils.jl b/src/utils.jl
index c4c9d906..b05fd171 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -971,7 +971,7 @@ function setup_reservoir_state(rmodel::SimulationModel; kwarg...)
         I = findfirst(isequal(k), pvars)
         if isnothing(I)
             if !(k in svars)
-                @warn "Recieved primary variable $k, but this is not known to reservoir model... Adding anyway."
+                jutul_message("setup_reservoir_state", "Recieved primary variable $k, but this is not known to reservoir model.")
             end
         else
             push!(found, k)
@@ -1158,9 +1158,7 @@ function well_output(model::MultiModel, states, well_symbol, forces, target = Bo
                 gforce = force[Symbol("$(well_symbol)_ctrl")]
             end
             control = gforce.control[well_symbol]
-            if control isa InjectorControl || control isa ProducerControl
-                q_t *= control.factor
-            end
+            q_t = effective_surface_rate(q_t, control)
             if target == :TotalSurfaceMassRate
                 d[i] = q_t
             elseif target isa Int

From 9bc235f469a2d1469fc0eb56d3b767bd1c8567dd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Fri, 21 Jun 2024 11:32:02 +0200
Subject: [PATCH 23/36] Support WELTARG

---
 src/input_simulation/data_input.jl | 65 ++++++++++++++++++++++++++++--
 1 file changed, 61 insertions(+), 4 deletions(-)

diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl
index 35011ab0..ad4e2f6a 100644
--- a/src/input_simulation/data_input.jl
+++ b/src/input_simulation/data_input.jl
@@ -1368,14 +1368,17 @@ function parse_control_steps(runspec, props, schedule, sys)
                 for wk in kword
                     name = wk[1]
                     controls[name], limits[name] = keyword_to_control(sys, streams, wk, key, factor = well_factor[name], temperature = well_temp[name])
-                    if !(controls[name] isa DisabledControl)
-                        active_controls[name] = controls[name]
-                    end
                 end
+                set_active_controls!(active_controls, controls)
             elseif key == "WELOPEN"
                 for wk in kword
                     apply_welopen!(controls, compdat, wk, active_controls)
                 end
+            elseif key == "WELTARG"
+                for wk in kword
+                    controls, limits = apply_weltarg!(controls, limits, wk)
+                end
+                set_active_controls!(active_controls, controls)
             elseif key in skip
                 # Already handled
             elseif key in ("WSEGVALV", "COMPSEGS", "WELSEGS")
@@ -1401,7 +1404,24 @@ function parse_control_steps(runspec, props, schedule, sys)
     for k in keys(bad_kw)
         jutul_message("Unsupported keyword", "Keyword $k was present, but is not supported.", color = :yellow)
     end
-    return (dt = tstep, control_step = cstep, controls = all_controls, completions = all_compdat, multisegment = mswell_kw, limits = all_limits)
+    return (
+        dt = tstep,
+        control_step = cstep,
+        controls = all_controls,
+        completions = all_compdat,
+        multisegment = mswell_kw,
+        limits = all_limits
+    )
+end
+
+function set_active_controls!(active_controls, controls)
+    for (name, ctrl) in pairs(controls)
+        if ctrl isa DisabledControl
+            continue
+        end
+        active_controls[name] = ctrl
+    end
+    return active_controls
 end
 
 function setup_well_streams()
@@ -1852,6 +1872,43 @@ function well_completion_sortperm(domain, wspec, order_t0, wc, dir)
     return sorted
 end
 
+function apply_weltarg!(controls, limits, wk)
+    well, ctype, value = wk
+    ctype = lowercase(ctype)
+    ctrl = controls[well]
+    @assert !(ctrl isa DisabledControl) "Cannot use WELTARG on disabled well."
+    is_injector = ctrl isa InjectorControl
+    limit = limits[well]
+    limit = OrderedDict{Symbol, Float64}(pairs(limit))
+    if ctype == "bhp"
+        new_target = BottomHolePressureTarget(value)
+        limit[Symbol(ctype)] = value
+    else
+        if is_injector
+            rate_target = value
+        else
+            rate_target = -value
+        end
+        if ctype == "orat"
+            new_target = SurfaceOilRateTarget(rate_target)
+        elseif ctype == "grat"
+            new_target = SurfaceGasRateTarget(rate_target)
+        elseif ctype == "wrat"
+            new_target = SurfaceWaterRateTarget(rate_target)
+        elseif ctype == "rate"
+            new_target = TotalRateTarget(rate_target)
+        elseif ctype == "lrat"
+            new_target = SurfaceLiquidRateTarget(rate_target)
+        else
+            error("WELTARG $ctype is not yet supported.")
+        end
+        limit[Symbol(ctype)] = rate_target
+    end
+    limits[well] = convert_to_immutable_storage(limit)
+    controls[well] = replace_target(ctrl, new_target)
+    return (controls, limits)
+end
+
 function apply_welopen!(controls, compdat, wk, controls_if_active)
     name, flag, I, J, K, first_num, last_num = wk
     # TODO: Handle shut in a better way

From b84d6cbd093f4f2d237d901bbff0eff683fa0226 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Fri, 21 Jun 2024 13:08:16 +0200
Subject: [PATCH 24/36] Improve WEFAC behavior

---
 src/facility/cross_terms.jl        | 18 ++++++++++++++----
 src/input_simulation/data_input.jl |  8 +++++++-
 src/utils.jl                       |  1 -
 3 files changed, 21 insertions(+), 6 deletions(-)

diff --git a/src/facility/cross_terms.jl b/src/facility/cross_terms.jl
index 53409749..7dae8aa8 100644
--- a/src/facility/cross_terms.jl
+++ b/src/facility/cross_terms.jl
@@ -114,7 +114,7 @@ function Jutul.prepare_cross_term_in_entity!(i,
         limits = current_limits(cfg, well_symbol)
         if !isnothing(limits)
             rhoS, S = surface_density_and_volume_fractions(state_well)
-            q_t = facility_surface_mass_rate_for_well(facility, well_symbol, state_facility, effective = true)
+            q_t = facility_surface_mass_rate_for_well(facility, well_symbol, state_facility, effective = false)
             apply_well_limit!(cfg, target, well, state_well, well_symbol, rhoS, S, value(q_t), limits)
         end
     end
@@ -137,7 +137,12 @@ function update_cross_term_in_entity!(out, i,
     ctrl = operating_control(cfg, well_symbol)
 
     target = ctrl.target
-    q_t = facility_surface_mass_rate_for_well(facility, well_symbol, state_facility, effective = true)
+    q_t = facility_surface_mass_rate_for_well(
+        facility,
+        well_symbol,
+        state_facility,
+        effective = false
+    )
     t, t_num = target_actual_pair(target, well, state_well, q_t, ctrl)
     t += 0*bottom_hole_pressure(state_well) + 0*q_t
     scale = target_scaling(target)
@@ -188,7 +193,12 @@ function update_cross_term_in_entity!(out, i,
 
     cfg = state_facility.WellGroupConfiguration
     ctrl = operating_control(cfg, well_symbol)
-    q_t = facility_surface_mass_rate_for_well(facility, well_symbol, state_facility, effective = true)
+    q_t = facility_surface_mass_rate_for_well(
+        facility,
+        well_symbol,
+        state_facility,
+        effective = true
+    )
     # Hack for sparsity detection
     q_t += 0*bottom_hole_pressure(state_well)
 
@@ -309,7 +319,7 @@ function update_cross_term_in_entity!(out, i,
 
     cfg = state_facility.WellGroupConfiguration
     ctrl = operating_control(cfg, well_symbol)
-    qT = state_facility.TotalSurfaceMassRate[pos] 
+    qT = state_facility.TotalSurfaceMassRate[pos]
     # Hack for sparsity detection
     qT += 0*bottom_hole_pressure(state_well)
 
diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl
index ad4e2f6a..702f2424 100644
--- a/src/input_simulation/data_input.jl
+++ b/src/input_simulation/data_input.jl
@@ -432,7 +432,8 @@ function parse_schedule(domain, runspec, props, schedule, sys; simple_well = tru
         wpi_mul = ones(n_wi)
         for (i, c) in enumerate(completions)
             compdat = c[k]
-            well_is_shut = controls[i][k] isa DisabledControl
+            well_control = controls[i][k]
+            well_is_shut = well_control isa DisabledControl
             wi_mul = zeros(n_wi)
             if !well_is_shut
                 wc, WI, open, multipliers, fresh = compdat_to_connection_factors(domain, wspec, compdat, i, sort = false)
@@ -455,6 +456,11 @@ function parse_schedule(domain, runspec, props, schedule, sys; simple_well = tru
                     wi_mul[compl_idx] = wpi_mul[compl_idx]*wi*is_open/WI_base[compl_idx]
                 end
             end
+            if !well_is_shut
+                # Multiply by well factor to account for downtime in
+                # pressure drop when well is not always active
+                @. wi_mul *= well_control.factor
+            end
             if all(isapprox(1.0), wi_mul)
                 mask = nothing
             else
diff --git a/src/utils.jl b/src/utils.jl
index b05fd171..e72ad51e 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1158,7 +1158,6 @@ function well_output(model::MultiModel, states, well_symbol, forces, target = Bo
                 gforce = force[Symbol("$(well_symbol)_ctrl")]
             end
             control = gforce.control[well_symbol]
-            q_t = effective_surface_rate(q_t, control)
             if target == :TotalSurfaceMassRate
                 d[i] = q_t
             elseif target isa Int

From 1dd1bc24cb1a6b7af05480aa14bd99c82bbb98c5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Sat, 22 Jun 2024 19:17:28 +0200
Subject: [PATCH 25/36] Input: Handle disconnection of shut/non-existent wells

---
 src/input_simulation/data_input.jl | 19 ++++++++++++++-----
 1 file changed, 14 insertions(+), 5 deletions(-)

diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl
index 702f2424..041b1127 100644
--- a/src/input_simulation/data_input.jl
+++ b/src/input_simulation/data_input.jl
@@ -404,7 +404,7 @@ end
 function parse_schedule(domain, runspec, props, schedule, sys; simple_well = true)
     G = physical_representation(domain)
 
-    dt, cstep, controls, completions, msdata, limits = parse_control_steps(runspec, props, schedule, sys)
+    dt, cstep, controls, status, completions, msdata, limits = parse_control_steps(runspec, props, schedule, sys)
     completions, bad_wells = filter_inactive_completions!(completions, G)
     @assert length(controls) == length(completions)
     handle_wells_without_active_perforations!(bad_wells, completions, controls, limits)
@@ -433,6 +433,7 @@ function parse_schedule(domain, runspec, props, schedule, sys; simple_well = tru
         for (i, c) in enumerate(completions)
             compdat = c[k]
             well_control = controls[i][k]
+            flag = status[i][k]
             well_is_shut = well_control isa DisabledControl
             wi_mul = zeros(n_wi)
             if !well_is_shut
@@ -461,6 +462,9 @@ function parse_schedule(domain, runspec, props, schedule, sys; simple_well = tru
                 # pressure drop when well is not always active
                 @. wi_mul *= well_control.factor
             end
+            if flag == "SHUT"
+                @. wi_mul = 0.0
+            end
             if all(isapprox(1.0), wi_mul)
                 mask = nothing
             else
@@ -1222,7 +1226,8 @@ function parse_control_steps(runspec, props, schedule, sys)
     active_controls = sdict()
     limits = sdict()
     streams = sdict()
-    mswell_kw = sdict()
+    status = sdict()
+    mswell_kw = Dict{String, String}()
     function get_and_create_mswell_kw(k::AbstractString, subkey = missing)
         if !haskey(mswell_kw, k)
             mswell_kw[k] = sdict()
@@ -1247,6 +1252,7 @@ function parse_control_steps(runspec, props, schedule, sys)
         active_controls[k] = DisabledControl()
         limits[k] = nothing
         streams[k] = nothing
+        status[k] = "SHUT"
         well_injection[k] = nothing
         well_factor[k] = 1.0
         # 0 deg C is strange, but the default for .DATA files.
@@ -1255,6 +1261,7 @@ function parse_control_steps(runspec, props, schedule, sys)
     all_compdat = []
     all_controls = []
     all_limits = []
+    all_status = []
 
     if haskey(runspec, "START")
         start_date = runspec["START"]
@@ -1373,7 +1380,7 @@ function parse_control_steps(runspec, props, schedule, sys)
             elseif key in ("WCONINJE", "WCONPROD", "WCONHIST", "WCONINJ", "WCONINJH")
                 for wk in kword
                     name = wk[1]
-                    controls[name], limits[name] = keyword_to_control(sys, streams, wk, key, factor = well_factor[name], temperature = well_temp[name])
+                    controls[name], limits[name], status[name] = keyword_to_control(sys, streams, wk, key, factor = well_factor[name], temperature = well_temp[name])
                 end
                 set_active_controls!(active_controls, controls)
             elseif key == "WELOPEN"
@@ -1403,6 +1410,7 @@ function parse_control_steps(runspec, props, schedule, sys)
             push!(all_compdat, deepcopy(compdat))
             push!(all_controls, deepcopy(controls))
             push!(all_limits, deepcopy(limits))
+            push!(all_status, deepcopy(status))
         else
             @warn "Did not find supported time kw in step $ctrl_ix: Keys were $(keys(step))."
         end
@@ -1414,6 +1422,7 @@ function parse_control_steps(runspec, props, schedule, sys)
         dt = tstep,
         control_step = cstep,
         controls = all_controls,
+        status = all_status,
         completions = all_compdat,
         multisegment = mswell_kw,
         limits = all_limits
@@ -1580,7 +1589,7 @@ function producer_control(sys, flag, ctrl, orat, wrat, grat, lrat, bhp; is_hist
             lims = producer_limits(bhp = bhp, orat = orat, wrat = wrat, grat = grat, lrat = lrat)
         end
     end
-    return (ctrl, lims)
+    return (ctrl, lims, flag)
 end
 
 function injector_limits(; bhp = Inf, surface_rate = Inf, reservoir_rate = Inf)
@@ -1642,7 +1651,7 @@ function injector_control(sys, streams, name, flag, type, ctype, surf_rate, res_
         end
         lims = injector_limits(bhp = bhp_lim, surface_rate = surf_rate, reservoir_rate = res_rate)
     end
-    return (ctrl, lims)
+    return (ctrl, lims, flag)
 end
 
 function select_injector_mixture_spec(sys::Union{ImmiscibleSystem, StandardBlackOilSystem}, name, streams, type)

From bcbc6060ed48ccae770a7bb21f656a671697880f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Mon, 24 Jun 2024 11:09:47 +0200
Subject: [PATCH 26/36] Account for upstream changes

---
 src/facility/controls.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/facility/controls.jl b/src/facility/controls.jl
index 2075826b..7640ea1d 100644
--- a/src/facility/controls.jl
+++ b/src/facility/controls.jl
@@ -3,11 +3,12 @@ function Jutul.initialize_extra_state_fields!(state, domain::WellGroup, model)
     state[:WellGroupConfiguration] = WellGroupConfiguration(domain.well_symbols)
 end
 
-function Jutul.update_before_step_multimodel!(storage_g, model_g::MultiModel, model::WellGroupModel, dt, forces, key;
+function Jutul.update_before_step_multimodel!(storage_g, model_g::MultiModel, model::WellGroupModel, dt, forces_g, key;
         time = NaN,
         recorder = ProgressRecorder(),
         update_explicit = true
     )
+    forces = forces_g[key]
     # Set control to whatever is on the forces
     storage = storage_g[key]
     cfg = storage.state.WellGroupConfiguration

From ae4940de3de91006a23f4cf7f9620e15d30e1f67 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Mon, 24 Jun 2024 11:55:55 +0200
Subject: [PATCH 27/36] Pass mask onto cdp

---
 src/facility/controls.jl       |  3 ++-
 src/facility/wells/stdwells.jl | 24 +++++++++++++++++-------
 2 files changed, 19 insertions(+), 8 deletions(-)

diff --git a/src/facility/controls.jl b/src/facility/controls.jl
index 7640ea1d..6bbb4075 100644
--- a/src/facility/controls.jl
+++ b/src/facility/controls.jl
@@ -52,9 +52,10 @@ function Jutul.update_before_step_multimodel!(storage_g, model_g::MultiModel, mo
     for wname in model.domain.well_symbols
         wmodel = model_g[wname]
         wstate = storage_g[wname].state
+        mask = forces_g[wname].mask
         rmodel = model_g[:Reservoir]
         rstate = storage_g.Reservoir.state
-        update_before_step_well!(wstate, wmodel, rstate, rmodel, op_ctrls[wname], update_explicit = update_explicit)
+        update_before_step_well!(wstate, wmodel, rstate, rmodel, op_ctrls[wname], mask, update_explicit = update_explicit)
     end
 end
 
diff --git a/src/facility/wells/stdwells.jl b/src/facility/wells/stdwells.jl
index 91d26d87..df5ab89b 100644
--- a/src/facility/wells/stdwells.jl
+++ b/src/facility/wells/stdwells.jl
@@ -81,14 +81,14 @@ end
 well_has_explicit_pressure_drop(m::SimpleWellFlowModel) = well_has_explicit_pressure_drop(physical_representation(m.domain))
 well_has_explicit_pressure_drop(w::SimpleWell) = w.explicit_dp
 
-function update_before_step_well!(well_state, well_model::SimpleWellFlowModel, res_state, res_model, ctrl; update_explicit = true)
+function update_before_step_well!(well_state, well_model::SimpleWellFlowModel, res_state, res_model, ctrl, mask; update_explicit = true)
     if well_has_explicit_pressure_drop(well_model) && update_explicit
         dp = well_state.ConnectionPressureDrop
-        update_connection_pressure_drop!(dp, well_state, well_model, res_state, res_model, ctrl)
+        update_connection_pressure_drop!(dp, well_state, well_model, res_state, res_model, ctrl, mask)
     end
 end
 
-function update_connection_pressure_drop!(dp, well_state, well_model, res_state, res_model, ctrl::InjectorControl)
+function update_connection_pressure_drop!(dp, well_state, well_model, res_state, res_model, ctrl::InjectorControl, mask)
     # Traverse down the well, using the phase notion encoded in ctrl and then
     # just accumulate pressure drop as we go assuming no cross flow
     phases = ctrl.phases
@@ -118,7 +118,7 @@ function update_connection_pressure_drop!(dp, well_state, well_model, res_state,
     end
 end
 
-function update_connection_pressure_drop!(dp, well_state, well_model, res_state, res_model, ctrl)
+function update_connection_pressure_drop!(dp, well_state, well_model, res_state, res_model, ctrl, mask)
     # Well is either disabled or producing. Loop over well from the bottom,
     # aggregating mixture density as we go. Then traverse down from the top and
     # accumulate the actual pressure drop due to hydrostatic assumptions.
@@ -129,21 +129,31 @@ function update_connection_pressure_drop!(dp, well_state, well_model, res_state,
     WI = as_value(well_state.WellIndices)
     ρ = as_value(res_state.PhaseMassDensities)
     mob = as_value(res_state.PhaseMobilities)
-
+    p = as_value(res_state.Pressure)
+    bhp = value(well_state.Pressure[1])
     # Integrate up, adding weighted density into well bore and keeping track of
     # current weight
     current_weight = 0.0
     current_density = 0.0
+    first_step = all(isequal(0.0), dp)
     for i in reverse(eachindex(dp))
         rc = res_cells[i]
         wi = WI[i]
-
+        if !isnothing(mask)
+            wi *= mask.values[i]
+        end
+        if first_step
+            pot = 1.0
+        else
+            pot = abs(bhp + dp[i] - p[rc])
+        end
+        q_perf = wi*pot
         # Mixture density along well bore
         local_density = 0
         local_weight = 0
         for ph in axes(ρ, 1)
             λ = mob[ph, rc]
-            weight_ph = wi*λ
+            weight_ph = q_perf*λ
             local_weight += weight_ph
             local_density += weight_ph*ρ[ph, rc]
         end

From 6ebbb7560f201c3416675799a268b88203a18ecd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Mon, 24 Jun 2024 11:56:18 +0200
Subject: [PATCH 28/36] Fix copy bug in compositional names array

---
 src/types.jl | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/src/types.jl b/src/types.jl
index 55f08ccf..456eb206 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -39,7 +39,7 @@ const LVCompositionalModel3Phase = SimulationModel{D, S, F, C} where {D, S<:LVCo
 Set up a compositional system for a given `equation_of_state` from `MultiComponentFlash`.
 """
 function MultiPhaseCompositionalSystemLV(equation_of_state, phases = (LiquidPhase(), VaporPhase()); reference_densities = ones(length(phases)), other_name = "Water")
-    c = MultiComponentFlash.component_names(equation_of_state)
+    c = copy(MultiComponentFlash.component_names(equation_of_state))
     N = length(c)
     phases = tuple(phases...)
     T = typeof(phases)
@@ -60,16 +60,18 @@ function MultiPhaseCompositionalSystemLV(equation_of_state, phases = (LiquidPhas
 end
 
 function Base.show(io::IO, sys::MultiPhaseCompositionalSystemLV)
+    components = copy(sys.components)
     n = number_of_components(sys)
     if has_other_phase(sys)
-        name = "(with water)"
+        name = "(three-phase)"
+        components[end] = "and $(components[end]) as immiscible phase"
         n = n - 1
     else
-        name = "(no water)"
+        name = "(two-phase)"
     end
     eos = sys.equation_of_state
-    cnames = join(MultiComponentFlash.component_names(eos), ", ")
-    print(io, "MultiPhaseCompositionalSystemLV $name with $(MultiComponentFlash.eostype(eos)) EOS with $n components: $cnames")
+    cnames = join(components, ", ")
+    print(io, "MultiPhaseCompositionalSystemLV $name with $(MultiComponentFlash.eostype(eos)) EOS with $n EOS components: $cnames")
 end
 
 struct StandardBlackOilSystem{D, V, W, R, F, T, P, Num} <: BlackOilSystem

From 3b8bda2292738130e6d525459e180448e7d9473a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Mon, 24 Jun 2024 11:57:46 +0200
Subject: [PATCH 29/36] Fix duplicate names in component_names for water LV
 system

---
 src/multicomponent/utils.jl | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/multicomponent/utils.jl b/src/multicomponent/utils.jl
index ba89b5ac..ec299b4f 100644
--- a/src/multicomponent/utils.jl
+++ b/src/multicomponent/utils.jl
@@ -11,11 +11,7 @@ function number_of_components(sys::MultiPhaseCompositionalSystemLV{E, T, O, G, N
 end
 
 function component_names(sys::MultiPhaseCompositionalSystemLV{E, T, O, G, N}) where {E, T, O, G, N}
-    names = copy(sys.components)
-    if has_other_phase(sys)
-        push!(names, phase_name(O()))
-    end
-    return names
+    return copy(sys.components)
 end
 
 phase_index(sys, phase) = only(findfirst(isequal(phase), sys.phases))

From 77153225507f54c60dc1af33bd2b1e23dfbe3161 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Mon, 24 Jun 2024 13:23:00 +0200
Subject: [PATCH 30/36] Update wells.jl

---
 src/facility/wells/wells.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/facility/wells/wells.jl b/src/facility/wells/wells.jl
index 18d3f8b0..ca2ff90a 100644
--- a/src/facility/wells/wells.jl
+++ b/src/facility/wells/wells.jl
@@ -239,7 +239,7 @@ end
 
 include("mswells_equations.jl")
 
-function update_before_step_well!(well_state, well_model, res_state, res_model, ctrl; kwarg...)
+function update_before_step_well!(well_state, well_model, res_state, res_model, ctr, mask; kwarg...)
 
 end
 

From 9103e083afe7d1fd86ac56226a7207d57f8b7d91 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Mon, 24 Jun 2024 14:54:02 +0200
Subject: [PATCH 31/36] Improved initialization code for ZMFVD

---
 src/init/init.jl                   | 20 +++++++++++++++++---
 src/input_simulation/data_input.jl |  9 +++++----
 2 files changed, 22 insertions(+), 7 deletions(-)

diff --git a/src/init/init.jl b/src/init/init.jl
index b044dfc9..69972f07 100644
--- a/src/init/init.jl
+++ b/src/init/init.jl
@@ -216,7 +216,7 @@ function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pr
 end
 
 
-function parse_state0_equil(model, datafile)
+function parse_state0_equil(model, datafile; normalize = :sum)
     has_water = haskey(datafile["RUNSPEC"], "WATER")
     has_oil = haskey(datafile["RUNSPEC"], "OIL")
     has_gas = haskey(datafile["RUNSPEC"], "GAS")
@@ -448,9 +448,20 @@ function parse_state0_equil(model, datafile)
                     Sz = SVector{N, Float64}
                     ztab_z = Float64[]
                     ztab_static = Sz[]
-                    for row in size(ztab, 1)
+                    num_renormalized = 0
+                    for row in 1:size(ztab, 1)
                         mf_i = Sz(ztab[row, 2:end])
-                        mf_i = mf_i./sum(mf_i)
+                        if maximum(mf_i) > 1.0 || minimum(mf) < 0 || !(sum(mf_i) ≈ 1.0)
+                            num_renormalized += 1
+                        end
+                        if normalize == :sum || normalize == true
+                            mf_i = mf_i./sum(mf_i)
+                        elseif normalize == :clampsum
+                            mf_i = clamp.(mf_i, 0.0, 1.0)
+                            mf_i = mf_i./sum(mf_i)
+                        else
+                            @assert normalize == false
+                        end
                         z_i = ztab[row, 1]
                         push!(ztab_z, z_i)
                         push!(ztab_static, mf_i)
@@ -463,6 +474,9 @@ function parse_state0_equil(model, datafile)
                             push!(ztab_static, mf_i)
                         end
                     end
+                    if num_renormalized > 0 && normalize != false
+                        jutul_message("Initialization", "$num_renormalized entries ZMFVD were normalized in equilibriation region $ereg.")
+                    end
                     composition = get_1d_interpolator(ztab_z, ztab_static)
                 else
                     composition = missing
diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl
index 041b1127..14a2a59e 100644
--- a/src/input_simulation/data_input.jl
+++ b/src/input_simulation/data_input.jl
@@ -66,6 +66,7 @@ function setup_case_from_parsed_data(datafile;
         use_ijk_trans = true,
         verbose = false,
         zcorn_depths = true,
+        normalize = true,
         kwarg...
     )
     function msg(s)
@@ -177,7 +178,7 @@ function setup_case_from_parsed_data(datafile;
     msg("Setting up forces.")
     forces = parse_forces(model, wells, controls, limits, cstep, dt, well_forces)
     msg("Setting up initial state.")
-    state0 = parse_state0(model, datafile)
+    state0 = parse_state0(model, datafile, normalize = normalize)
     msg("Setup complete.")
     return JutulCase(model, dt, forces, state0 = state0, parameters = parameters, input_data = datafile)
 end
@@ -558,13 +559,13 @@ function parse_forces(model, wells, controls, limits, cstep, dt, well_forces)
     return forces[cstep]
 end
 
-function parse_state0(model, datafile)
+function parse_state0(model, datafile; normalize = true)
     rmodel = reservoir_model(model)
     init = Dict{Symbol, Any}()
     sol = datafile["SOLUTION"]
 
     if haskey(sol, "EQUIL")
-        init = parse_state0_equil(rmodel, datafile)
+        init = parse_state0_equil(rmodel, datafile; normalize = normalize)
     else
         init = parse_state0_direct_assignment(rmodel, datafile)
     end
@@ -1227,7 +1228,7 @@ function parse_control_steps(runspec, props, schedule, sys)
     limits = sdict()
     streams = sdict()
     status = sdict()
-    mswell_kw = Dict{String, String}()
+    mswell_kw = Dict{String, sdict}()
     function get_and_create_mswell_kw(k::AbstractString, subkey = missing)
         if !haskey(mswell_kw, k)
             mswell_kw[k] = sdict()

From 6c3b73227d55bc651dec2ac4912423cfaef993a8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Mon, 24 Jun 2024 15:13:05 +0200
Subject: [PATCH 32/36] Update controls.jl

---
 src/facility/controls.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/facility/controls.jl b/src/facility/controls.jl
index 6bbb4075..71f065e9 100644
--- a/src/facility/controls.jl
+++ b/src/facility/controls.jl
@@ -52,7 +52,12 @@ function Jutul.update_before_step_multimodel!(storage_g, model_g::MultiModel, mo
     for wname in model.domain.well_symbols
         wmodel = model_g[wname]
         wstate = storage_g[wname].state
-        mask = forces_g[wname].mask
+        forces_w = forces_g[wname]
+        if isnothing(forces_w)
+            mask = nothing
+        else
+            mask = forces_g[wname].mask
+        end
         rmodel = model_g[:Reservoir]
         rstate = storage_g.Reservoir.state
         update_before_step_well!(wstate, wmodel, rstate, rmodel, op_ctrls[wname], mask, update_explicit = update_explicit)

From b39425cec0952f38ee97d799d3d532ca42041417 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Tue, 25 Jun 2024 14:34:21 +0200
Subject: [PATCH 33/36] Move HYPRE into parray extension, bump versions

---
 Project.toml                                  |  9 ++++----
 ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl  | 22 -------------------
 .../JutulDarcyPartitionedArraysExt.jl         |  3 +++
 .../cpr.jl                                    |  0
 4 files changed, 7 insertions(+), 27 deletions(-)
 delete mode 100644 ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl
 rename ext/{JutulDarcyHYPREExt => JutulDarcyPartitionedArraysExt}/cpr.jl (100%)

diff --git a/Project.toml b/Project.toml
index eb7404cd..b8218e03 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "JutulDarcy"
 uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a"
 authors = ["Olav Møyner <olav.moyner@gmail.com>"]
-version = "0.2.24"
+version = "0.2.25"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
@@ -38,7 +38,6 @@ PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 
 [extensions]
 JutulDarcyGLMakieExt = "GLMakie"
-JutulDarcyHYPREExt = "HYPRE"
 JutulDarcyPartitionedArraysExt = ["PartitionedArrays", "MPI", "HYPRE"]
 
 [compat]
@@ -49,14 +48,14 @@ DelimitedFiles = "1.9.1"
 DocStringExtensions = "0.9.3"
 ForwardDiff = "0.10.30"
 GLMakie = "0.9.2"
-GeoEnergyIO = "1.1.2"
+GeoEnergyIO = "1.1.3"
 HYPRE = "1.4.0"
-Jutul = "0.2.33"
+Jutul = "0.2.34"
 Krylov = "0.9.1"
 LinearAlgebra = "1"
 LoopVectorization = "0.12.115"
 MAT = "0.10.3"
-MultiComponentFlash = "1.1.14"
+MultiComponentFlash = "1.1.15"
 OrderedCollections = "1.6.2"
 PartitionedArrays = "0.3.2"
 Polyester = "0.6.11, 0.7.3"
diff --git a/ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl b/ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl
deleted file mode 100644
index f89ab706..00000000
--- a/ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl
+++ /dev/null
@@ -1,22 +0,0 @@
-module JutulDarcyHYPREExt
-    using JutulDarcy
-    using HYPRE
-    using Jutul
-    using SparseArrays
-    using PrecompileTools
-    using TimerOutputs
-
-    timeit_debug_enabled() = Jutul.timeit_debug_enabled()
-
-    include("cpr.jl")
-
-    @compile_workload begin
-    #     targets = [(true, :csc), (true, :csr)]
-    #     # MPI, trivial partition
-    #     JutulDarcy.precompile_darcy_multimodels(targets,
-    #         dims = (4, 1, 1),
-    #         precond = :cpr,
-    #         amg_type = :hypre
-    #     )
-    end
-end
diff --git a/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl b/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl
index 274394a0..2856c1e2 100644
--- a/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl
+++ b/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl
@@ -5,6 +5,7 @@ module JutulDarcyPartitionedArraysExt
     # Specific dependencies
     using PartitionedArrays, MPI, HYPRE
     using LinearAlgebra
+    using SparseArrays
 
     import Jutul: JutulCase, PArrayBackend, JutulConfig, BoomerAMGPreconditioner
     import Jutul: PArraySimulator, MPISimulator, PArrayExecutor
@@ -18,6 +19,8 @@ module JutulDarcyPartitionedArraysExt
 
     timeit_debug_enabled() = Jutul.timeit_debug_enabled()
 
+    include("cpr.jl")
+
     function setup_reservoir_simulator_parray(
             case::JutulCase,
             backend::PArrayBackend;
diff --git a/ext/JutulDarcyHYPREExt/cpr.jl b/ext/JutulDarcyPartitionedArraysExt/cpr.jl
similarity index 100%
rename from ext/JutulDarcyHYPREExt/cpr.jl
rename to ext/JutulDarcyPartitionedArraysExt/cpr.jl

From 5ea6f135f76c8aba84005eb127d2ff5d12401146 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Tue, 25 Jun 2024 15:04:50 +0200
Subject: [PATCH 34/36] Update controls.jl

---
 src/facility/controls.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/facility/controls.jl b/src/facility/controls.jl
index 71f065e9..8f05fdf1 100644
--- a/src/facility/controls.jl
+++ b/src/facility/controls.jl
@@ -53,10 +53,10 @@ function Jutul.update_before_step_multimodel!(storage_g, model_g::MultiModel, mo
         wmodel = model_g[wname]
         wstate = storage_g[wname].state
         forces_w = forces_g[wname]
-        if isnothing(forces_w)
+        if isnothing(forces_w) || !haskey(forces_w, :mask)
             mask = nothing
         else
-            mask = forces_g[wname].mask
+            mask = forces_w.mask
         end
         rmodel = model_g[:Reservoir]
         rstate = storage_g.Reservoir.state

From 6b65c51d06f6f6266d4eb37c65c7694a77f166b9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Tue, 25 Jun 2024 15:20:28 +0200
Subject: [PATCH 35/36] Allow Rs extrapolation

This extrapolation is not entirely consistent. Maybe we want to improve this.
---
 src/deck_types.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/deck_types.jl b/src/deck_types.jl
index 82a20bca..5b611c25 100644
--- a/src/deck_types.jl
+++ b/src/deck_types.jl
@@ -462,7 +462,9 @@ function saturated_table(p, r)
         p = vcat([-1.0, 0.0], p)
         r = vcat([0.0, 0.0], r)
     end
-    return get_1d_interpolator(p, r, cap_end = true)
+    # TODO: This is a bit unclear if it is a good idea, but it is required for
+    # the SPE1 test case.
+    return get_1d_interpolator(p, r, cap_end = false)
 end
 
 pvt_table_vectors(pvt::PVTGTable) = (pvt.rv, pvt.pressure, pvt.sat_rv, pvt.pos)

From 636ed0b4e426b0fda617abe90f126adb75d44dfc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Olav=20M=C3=B8yner?= <olav.moyner@gmail.com>
Date: Tue, 25 Jun 2024 16:16:49 +0200
Subject: [PATCH 36/36] Update test

---
 test/nldd.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/nldd.jl b/test/nldd.jl
index dd161c1f..0a775622 100644
--- a/test/nldd.jl
+++ b/test/nldd.jl
@@ -63,7 +63,7 @@ end
         its_newton = get_its(results_newton)
 
         function test_solver_result(result)
-            @test its_newton > get_its(result)
+            @test its_newton >= get_its(result)
         end
         @testset "NLDD on $name" begin
             for solve_tol_mobility in [nothing, 0.01, 0.05]