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
diff --git a/src/cpr.jl b/src/cpr.jl
index 69e0cf1e..a6b91fea 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)
@@ -378,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
@@ -388,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
diff --git a/src/deck_types.jl b/src/deck_types.jl
index 962644ae..5b611c25 100644
--- a/src/deck_types.jl
+++ b/src/deck_types.jl
@@ -462,6 +462,8 @@ function saturated_table(p, r)
         p = vcat([-1.0, 0.0], p)
         r = vcat([0.0, 0.0], r)
     end
+    # 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
 
diff --git a/src/facility/controls.jl b/src/facility/controls.jl
index a38edf90..8f05fdf1 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
@@ -51,9 +52,15 @@ 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
+        forces_w = forces_g[wname]
+        if isnothing(forces_w) || !haskey(forces_w, :mask)
+            mask = nothing
+        else
+            mask = forces_w.mask
+        end
         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
 
@@ -198,9 +205,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..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)
+            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)
+    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,17 +193,17 @@ 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 +211,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 +223,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
 
@@ -314,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/facility/types.jl b/src/facility/types.jl
index 42ac7fb0..786d953f 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}
@@ -363,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/facility/wellgroups.jl b/src/facility/wellgroups.jl
index ee12cef0..f9467f0f 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,7 +25,7 @@ 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
+        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))
@@ -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))
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
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
 
diff --git a/src/init/init.jl b/src/init/init.jl
index b11c1ebe..69972f07 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
+            evaluate_relative_permeability_no_scaling!(kr, relperm, model, s_eval, cells)
         end
+        kr = kr[:, cells]
         init[:Saturations] = s
         init[:Pressure] = init_reference_pressure(pressures, contacts, kr, pc, 2)
     else
@@ -222,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")
@@ -297,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
@@ -426,7 +420,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 +432,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.
@@ -454,8 +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])
+                        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)
@@ -468,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 081f0b58..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)
@@ -160,7 +161,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
@@ -177,13 +178,13 @@ 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
 
-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 +348,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 +362,8 @@ function compdat_to_connection_factors(domain, wspec, v; sort = true, order = "T
     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)
         W_i = WI[i]
@@ -394,14 +397,15 @@ 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, mul, fresh)
 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)
@@ -424,24 +428,44 @@ 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)
+            well_control = controls[i][k]
+            flag = status[i][k]
+            well_is_shut = well_control isa DisabledControl
             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)
+                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
-                    wi_mul[compl_idx] = wi*is_open/WI_base[compl_idx]
+                    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
+            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 flag == "SHUT"
+                @. wi_mul = 0.0
+            end
             if all(isapprox(1.0), wi_mul)
                 mask = nothing
             else
@@ -535,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
@@ -817,12 +841,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
@@ -1086,7 +1113,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"])
@@ -1200,7 +1227,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, sdict}()
     function get_and_create_mswell_kw(k::AbstractString, subkey = missing)
         if !haskey(mswell_kw, k)
             mswell_kw[k] = sdict()
@@ -1225,6 +1253,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.
@@ -1233,6 +1262,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"]
@@ -1249,11 +1279,12 @@ 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}()
+    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]
@@ -1309,6 +1340,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,
@@ -1317,21 +1374,25 @@ function parse_control_steps(runspec, props, schedule, sys)
                             Kh = Kh,
                             skin = skin,
                             dir = dir,
+                            mul = perf_mult,
                             ctrl = ctrl_ix)
                     end
                 end
             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])
-                    if !(controls[name] isa DisabledControl)
-                        active_controls[name] = controls[name]
-                    end
+                    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"
                 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")
@@ -1350,6 +1411,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
@@ -1357,10 +1419,34 @@ 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,
+        status = all_status,
+        completions = all_compdat,
+        multisegment = mswell_kw,
+        limits = all_limits
+    )
 end
 
-function parse_well_streams_for_step(step, props)
+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()
+    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)
@@ -1371,8 +1457,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])
@@ -1498,15 +1582,15 @@ 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)
         end
     end
-    return (ctrl, lims)
+    return (ctrl, lims, flag)
 end
 
 function injector_limits(; bhp = Inf, surface_rate = Inf, reservoir_rate = Inf)
@@ -1568,7 +1652,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)
@@ -1629,7 +1713,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"
@@ -1804,6 +1888,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
diff --git a/src/input_simulation/mrst_input.jl b/src/input_simulation/mrst_input.jl
index 3f5e5c79..ed8a244e 100644
--- a/src/input_simulation/mrst_input.jl
+++ b/src/input_simulation/mrst_input.jl
@@ -364,17 +364,22 @@ 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, "SCALECRS")
-        scalecrs = props["SCALECRS"]
-        if scalecrs isa String
-            scalecrs = [scalecrs]
+    if haskey(runspec, "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
@@ -647,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
@@ -680,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
 
@@ -706,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)
@@ -773,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
@@ -1357,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
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()
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
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))
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
 
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}
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
diff --git a/src/utils.jl b/src/utils.jl
index dc8cb44d..e72ad51e 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,
@@ -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,6 @@ 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
             if target == :TotalSurfaceMassRate
                 d[i] = q_t
             elseif target isa Int
diff --git a/src/variables/relperm.jl b/src/variables/relperm.jl
index 840eceec..ad1a5e6f 100644
--- a/src/variables/relperm.jl
+++ b/src/variables/relperm.jl
@@ -401,40 +401,43 @@ 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)
-    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)
+    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
     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)
+    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
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]