Merge pull request #47 from sintefmath/dev

Improved performance, better behavior on input cases

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"

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;

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

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
 

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]

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,25 +193,25 @@ 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
         nmix = length(mix)
         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)
 

src/facility/types.jl

@@ -7,16 +7,19 @@ 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
 
 """
     WellGroup(wells::Vector{Symbol}; can_shut_wells = true)
 
 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)

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))

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

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