Merge pull request #51 from sintefmath/hyst

Hysteresis & relperm refactoring

Project.toml

@@ -1,7 +1,7 @@
 name = "JutulDarcy"
 uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a"
 authors = ["Olav Møyner <[email protected]>"]
-version = "0.2.26"
+version = "0.2.27"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
@@ -47,10 +47,10 @@ Dates = "1"
 DelimitedFiles = "1.9.1"
 DocStringExtensions = "0.9.3"
 ForwardDiff = "0.10.30"
-GLMakie = "0.9.2"
+GLMakie = "0.10.0"
 GeoEnergyIO = "1.1.3"
 HYPRE = "1.4.0"
-Jutul = "0.2.34"
+Jutul = "0.2.35"
 Krylov = "0.9.1"
 LinearAlgebra = "1"
 LoopVectorization = "0.12.115"
@@ -64,7 +64,7 @@ ProgressMeter = "1.7.2"
 SparseArrays = "1"
 StaticArrays = "1.4.4"
 Statistics = "1"
-TestItems = "0.1.0, 1"
+TestItems = "1.0.0"
 TimerOutputs = "0.5.19"
 Tullio = "0.3.4"
 julia = "1.7"

docs/src/man/basics/parameters.md

@@ -3,7 +3,7 @@
 ## General
 
 ```@docs
-FluidVolume
+JutulDarcy.FluidVolume
 ```
 
 ## Reservoir parameters
@@ -15,20 +15,20 @@ JutulDarcy.Transmissibilities
 ```
 
 ```@docs
-reservoir_transmissibility
+JutulDarcy.reservoir_transmissibility
 ```
 
 ### Other
 
 ```@docs
-TwoPointGravityDifference
-ConnateWater
-EndPointScalingCoefficients
+JutulDarcy.TwoPointGravityDifference
+JutulDarcy.ConnateWater
+JutulDarcy.EndPointScalingCoefficients
 ```
 
 ## Well parameters
 
 ```@docs
-WellIndices
-PerforationGravityDifference
+JutulDarcy.WellIndices
+JutulDarcy.PerforationGravityDifference
 ```

src/facility/controls.jl

@@ -179,8 +179,11 @@ function translate_limit(control::InjectorControl, name, val)
         # Lower limit, total volumetric surface rate
         target_limit = TotalRateTarget(val)
         is_lower = true
+    elseif name == :resv_rate
+        # Upper limit, total volumetric reservoir rate
+        target_limit = TotalReservoirRateTarget(val)
     else
-        error("$name limit not supported for well acting as producer.")
+        error("$name limit not supported for well acting as injector.")
     end
     return (target_limit, is_lower)
 end
@@ -246,6 +249,20 @@ function well_target(control::InjectorControl, target::SurfaceVolumeTarget, well
     return t/control.mixture_density
 end
 
+
+"""
+Well target contribution from well itself (reservoir volume, injector)
+"""
+function well_target(control::InjectorControl, target::TotalReservoirRateTarget, well_model, well_state, surface_densities, surface_volume_fractions)
+    w_phases = get_phases(flow_system(well_model.system))
+    t = 0.0
+    rho = well_state.PhaseMassDensities
+    for (ix, mix) in control.phases
+        t += mix/rho[ix, 1]
+    end
+    return t
+end
+
 """
 Well target contribution from well itself (surface volume, injector)
 """

src/facility/types.jl

@@ -215,6 +215,23 @@ struct TotalRateTarget{T} <: SurfaceVolumeTarget where T<:AbstractFloat
 end
 Base.show(io::IO, t::TotalRateTarget) = print(io, "TotalRateTarget with value $(t.value) [m^3/s]")
 
+"""
+    TotalReservoirRateTarget(q)
+
+Well target of specified total rate (sum of all phases) with value `q` at reservoir
+conditions.
+
+Often used for both [`InjectorControl`](@ref) [`ProducerControl`](@ref).
+"""
+struct TotalReservoirRateTarget{T} <: WellTarget where T<:AbstractFloat
+    value::T
+    function TotalReservoirRateTarget(v::T) where T
+        isfinite(v) || throw(ArgumentError("Rate must be finite, was $v"))'
+        return new{T}(v)
+    end
+end
+Base.show(io::IO, t::TotalReservoirRateTarget) = print(io, "TotalReservoirRateTarget with value $(t.value) [m^3/s]")
+
 """
     HistoricalReservoirVoidageTarget(q, weights)
 
@@ -565,6 +582,17 @@ function well_target_information(t::Union{TotalRateTarget, Val{:rate}})
     )
 end
 
+
+function well_target_information(t::Union{TotalReservoirRateTarget, Val{:resv_rate}})
+    return well_target_information(
+        symbol = :resv_rate,
+        description = "Surface total rate",
+        explanation = "Total volumetric rate at reservoir conditions. This is the sum of all phases. For most models, it is the sum of the mass rates divided by the prescribed surface densities.",
+        unit_type = :reservoir_volume_per_time,
+        unit_label = "m³/s"
+    )
+end
+
 function well_target_information(t::Union{SurfaceWaterRateTarget, Val{:wrat}})
     return well_target_information(
         symbol = :wrat,

src/init/init.jl

@@ -187,14 +187,28 @@ 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 length(phases) == 3 && AqueousPhase() in phases
+        phase_ind = phase_indices(model.system)
+        if length(phases) == 3
             swcon = zeros(nc_total)
-            if !ismissing(s_min)
+            if AqueousPhase() in phases && !ismissing(s_min)
                 swcon[cells] .= s_min[1]
             end
-            evaluate_relative_permeability_no_scaling!(kr, relperm, model, s_eval, cells, swcon)
+            for c in cells
+                update_three_phase_relperm!(kr, relperm, phase_ind, s_eval, nothing, c, swcon[c], nothing, nothing)
+            end
         else
-            evaluate_relative_permeability_no_scaling!(kr, relperm, model, s_eval, cells)
+            ph = relperm.phases
+            if ph == :wo
+                kr1, kr2 = relperm.krw, relperm.krow
+            elseif ph == :og
+                kr1, kr2 = relperm.krog, relperm.krg
+            else
+                @assert ph == :wg
+                kr1, kr2 = relperm.krw, relperm.krg
+            end
+            for c in cells
+                update_two_phase_relperm!(kr, relperm, kr1, kr2, nothing, nothing, phase_ind, s_eval, nothing, c, nothing, nothing)
+            end
         end
         kr = kr[:, cells]
         init[:Saturations] = s

src/input_simulation/data_input.jl

@@ -1601,8 +1601,8 @@ function injector_limits(; bhp = Inf, surface_rate = Inf, reservoir_rate = Inf)
     if isfinite(surface_rate)
         lims[:rate] = surface_rate
     end
-    if !isinf(reservoir_rate)
-        @warn "Non-defaulted reservoir rate limit not supported: $reservoir_rate"
+    if isfinite(reservoir_rate)
+        lims[:resv_rate] = reservoir_rate
     end
     return NamedTuple(pairs(lims))
 end
@@ -1629,6 +1629,13 @@ function injector_control(sys, streams, name, flag, type, ctype, surf_rate, res_
             else
                 well_is_disabled = true
             end
+        elseif ctype == "RESV"
+            is_rate = false
+            if isfinite(res_rate) && res_rate > MIN_ACTIVE_WELL_RATE
+                t = TotalReservoirRateTarget(res_rate)
+            else
+                well_is_disabled = true
+            end
         else
             # RESV, GRUP, THP
             error("$ctype control not supported")