diff --git a/src/blackoil/variables/varswitch.jl b/src/blackoil/variables/varswitch.jl
index 0d5a8d189..ac340c467 100644
--- a/src/blackoil/variables/varswitch.jl
+++ b/src/blackoil/variables/varswitch.jl
@@ -69,13 +69,18 @@ Base.@propagate_inbounds function varswitch_update_inner!(v, i, dx, dr_max, ds_m
     is_near_bubble = was_near_bubble
     if swi > 1 - 10*ϵ_s
         # Water filled
+        sg_max = 1.0 - swi
         if old_state == OilAndGas
-            next_x = old_x + w*Jutul.choose_increment(value(old_x), dx, ds_max, nothing, 0, 1 - swi)
+            if old_x > sg_max
+                next_x = replace_value(old_x, sg_max)
+            else
+                next_x = old_x + w*Jutul.choose_increment(value(old_x), dx, ds_max, nothing, 0, sg_max)
+            end
         else
             if old_state == OilOnly
-                sg = 0
+                sg = 0.0
             else
-                sg = 1 - swi
+                sg = sg_max
             end
             next_x = replace_value(old_x, sg)
         end
@@ -196,7 +201,7 @@ function blackoil_unknown_init(F_rs::Nothing, F_rv, sw, so, sg, rs, rv, p)
 end
 
 function blackoil_unknown_init(F_rs, F_rv, sw, so, sg, rs, rv, p)
-    if sg > 0 && so > 0
+    if sg > 0 && so > 0 || sw > 1.0 - 10*MINIMUM_COMPOSITIONAL_SATURATION
         rs = F_rs(p)
         rv = F_rv(p)
         x = sg
@@ -252,7 +257,7 @@ end
         sw = ImmiscibleSaturation[i]
         X = BlackOilUnknown[i]
         phases = X.phases_present
-        rem = one(T) - sw + MINIMUM_COMPOSITIONAL_SATURATION
+        rem = one(T) - sw# + MINIMUM_COMPOSITIONAL_SATURATION
         if phases == OilOnly
             sg = zero(T)
             so = rem
diff --git a/src/types.jl b/src/types.jl
index 1e924ac7e..9a61a8aa2 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -127,7 +127,7 @@ function StandardBlackOilSystem(;
         reference_densities = [786.507, 1037.84, 0.969758], 
         saturated_chop = false,
         keep_bubble_flag = true,
-        eps_s = 1e-10,
+        eps_s = MINIMUM_COMPOSITIONAL_SATURATION,
         eps_rs = nothing,
         eps_rv = nothing,
         formulation::Symbol = :varswitch,
@@ -352,7 +352,7 @@ function PhaseRelativePermeability(s, k; label = :w, connate = s[1], epsilon = 1
     crit_ix = findfirst(x -> x > eps(Float64), k) - 1
     crit = s[crit_ix]
     s, k = JutulDarcy.add_missing_endpoints(s, k)
-    JutulDarcy.ensure_endpoints!(s, k, epsilon)
+    # JutulDarcy.ensure_endpoints!(s, k, epsilon)
     kr = get_1d_interpolator(s, k, cap_endpoints = false, constant_dx = false)
     return PhaseRelativePermeability(kr, label, connate, crit, s_max, k_max, s_max_table)
 end
diff --git a/src/utils.jl b/src/utils.jl
index be9874d47..d7957d5e7 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -398,6 +398,7 @@ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem;
         immutable_model = false,
         wells_systems = missing,
         wells_as_cells = false,
+        optimization_level = 0,
         discretization_arg = NamedTuple()
     )
     # Deal with wells, make sure that multisegment wells come last.
@@ -441,7 +442,8 @@ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem;
         context = reservoir_context,
         general_ad = general_ad,
         kgrad = kgrad,
-        upwind = upwind
+        upwind = upwind,
+        optimization_level = optimization_level
     )
     system = rmodel.system
     if thermal
@@ -511,7 +513,10 @@ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem;
                     w_domain[propk] = reservoir[propk][c]
                 end
             end
-            wmodel = SimulationModel(w_domain, system, context = well_context)
+            wmodel = SimulationModel(w_domain, system,
+                context = well_context,
+                optimization_level = optimization_level
+            )
             if thermal
                 wmodel = add_thermal_to_model!(wmodel)
             end
@@ -531,7 +536,11 @@ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem;
             models[wname] = wmodel
             if split_wells
                 wg = WellGroup([wname], can_shut_wells = can_shut_wells)
-                F = SimulationModel(wg, facility_system, context = context, data_domain = DataDomain(wg))
+                F = SimulationModel(wg, facility_system,
+                    context = context,
+                    data_domain = DataDomain(wg),
+                    optimization_level = optimization_level
+                )
                 if thermal
                     add_thermal_to_facility!(F)
                 end
@@ -546,7 +555,11 @@ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem;
             end
         else
             wg = WellGroup(map(x -> physical_representation(x).name, wells), can_shut_wells = can_shut_wells)
-            F = SimulationModel(wg, facility_system, context = context, data_domain = DataDomain(wg))
+            F = SimulationModel(wg, facility_system,
+                context = context,
+                data_domain = DataDomain(wg),
+                optimization_level = optimization_level
+            )
             if thermal
                 add_thermal_to_facility!(F)
             end