Better compositional solvers and fixes to simple wells + blackoil diffusion

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.9"
+version = "0.2.10"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
@@ -43,7 +43,7 @@ DataStructures = "0.18.13"
 DelimitedFiles = "1.9.1"
 ForwardDiff = "0.10.30"
 HYPRE = "1.4.0"
-Jutul = "0.2.14"
+Jutul = "0.2.15"
 Krylov = "0.9.1"
 LoopVectorization = "0.12.115"
 MAT = "0.10.3"

ext/JutulDarcyGLMakieExt/well_plots.jl

@@ -58,15 +58,17 @@ function JutulDarcy.plot_well_results(well_data::Dict, arg...; name = "Data", kw
     JutulDarcy.plot_well_results([well_data], arg...; names = [name], kwarg...)
 end
 
-function JutulDarcy.plot_well_results(well_data::Vector, time = nothing; start_date = nothing,
-                                                              names =["Dataset $i" for i in 1:length(well_data)], 
-                                                              linewidth = 3,
-                                                              cmap = nothing, 
-                                                              dashwidth = 1,
-                                                              new_window = false,
-                                                              styles = [:solid, :scatter, :dash, :dashdot, :dot, :dashdotdot],
-                                                              resolution = (1600, 900),
-                                                              kwarg...)
+function JutulDarcy.plot_well_results(well_data::Vector, time = nothing;
+        start_date = nothing,
+        names =["Dataset $i" for i in 1:length(well_data)], 
+        linewidth = 3,
+        cmap = nothing, 
+        dashwidth = 1,
+        new_window = false,
+        styles = [:solid, :scatter, :dash, :dashdot, :dot, :dashdotdot],
+        resolution = (1600, 900),
+        kwarg...
+    )
     # Figure part
     names = Vector{String}(names)
     ndata = length(well_data)
@@ -77,6 +79,9 @@ function JutulDarcy.plot_well_results(well_data::Vector, time = nothing; start_d
     if nw == 0
         return nothing
     end
+    # Type of plot (bhp, rate...)
+    responses = collect(keys(wd[first(wells)]))
+    respstr = [String(x) for x in responses]
 
     is_inj = is_injectors(wd)
     @assert ndata <= length(styles) "Can't plot more datasets than styles provided"
@@ -90,7 +95,27 @@ function JutulDarcy.plot_well_results(well_data::Vector, time = nothing; start_d
     else
         t_l = "Date"
     end
-    ax = Axis(fig[1, 1], xlabel = t_l)
+    function response_label_to_unit(s)
+        s = "$s"
+        rate_labels = [
+            "Surface total rate", "rate",
+            "Surface water rate", "wrat",
+            "Surface liquid rate (water + oil)", "lrat",
+            "Surface oil rate", "orat",
+            "Surface gas rate", "grat",
+            "Reservoir voidage rate", "resv",
+            "Historical reservoir voidage rate", "resv"
+        ]
+        if s in rate_labels
+            return "m^3/s"
+        elseif s in ["bhp", "Bottom hole pressure"]
+            return "Pa"
+        else
+            return ""
+        end
+    end
+    y_l = Observable(response_label_to_unit(first(responses)))
+    ax = Axis(fig[1, 1], xlabel = t_l, ylabel = y_l)
 
     if isnothing(cmap)
         if nw > 20
@@ -104,16 +129,14 @@ function JutulDarcy.plot_well_results(well_data::Vector, time = nothing; start_d
     end
     wellstr = [String(x) for x in wells]
 
-    # Type of plot (bhp, rate...)
-    responses = collect(keys(wd[first(wells)]))
-    respstr = [String(x) for x in responses]
     response_ix = Observable(1)
     is_accum = Observable(false)
     is_abs = Observable(false)
     type_menu = Menu(fig, options = respstr, prompt = respstr[1])
 
     on(type_menu.selection) do s
         val = findfirst(isequal(s), respstr)
+        y_l[] = response_label_to_unit(s)
         response_ix[] = val
         autolimits!(ax)
     end

ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl

@@ -6,13 +6,13 @@ module JutulDarcyHYPREExt
     using PrecompileTools
     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
+    # @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

ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl

@@ -131,28 +131,28 @@ module JutulDarcyPartitionedArraysExt
         return (cpr, preconditioners)
     end
 
-    @compile_workload begin
-        targets = [(true, :csc), (true, :csr)]
-        # MPI, trivial partition
-        JutulDarcy.precompile_darcy_multimodels(targets,
-            dims = (4, 1, 1),
-            default_linsolve = false,
-            setuparg = (
-                mode = :mpi,
-                precond = :ilu0
-                ),
-            split_wells = true
-        )
-        # Native PArray, non-trivial partition
-        JutulDarcy.precompile_darcy_multimodels(targets,
-            dims = (4, 1, 1),
-            default_linsolve = false,
-            setuparg = (
-                mode = :parray,
-                parray_arg = (np = 2, ),
-                precond = :ilu0
-                ),
-            split_wells = true
-        )
-    end
+    # @compile_workload begin
+    #     targets = [(true, :csc), (true, :csr)]
+    #     # MPI, trivial partition
+    #     JutulDarcy.precompile_darcy_multimodels(targets,
+    #         dims = (4, 1, 1),
+    #         default_linsolve = false,
+    #         setuparg = (
+    #             mode = :mpi,
+    #             precond = :ilu0
+    #             ),
+    #         split_wells = true
+    #     )
+    #     # Native PArray, non-trivial partition
+    #     JutulDarcy.precompile_darcy_multimodels(targets,
+    #         dims = (4, 1, 1),
+    #         default_linsolve = false,
+    #         setuparg = (
+    #             mode = :parray,
+    #             parray_arg = (np = 2, ),
+    #             precond = :ilu0
+    #             ),
+    #         split_wells = true
+    #     )
+    # end
 end

src/InputParser/deckinput/keywords/schedule.jl

@@ -160,6 +160,7 @@ end
 function convert_date_kw(t)
     @assert length(t) == 4
     function get_month(s)
+        s = uppercase(s)
         if s == "JAN"
             return 1
         elseif s == "FEB"
@@ -183,7 +184,7 @@ function convert_date_kw(t)
         elseif s == "NOV"
             return 11
         else
-            @assert s == "DEC"
+            @assert s == "DEC" "Did not understand month format $s"
             return 12
         end
     end

src/blackoil/flux.jl

@@ -30,8 +30,7 @@
     if has_vapoil(sys)
         # Rv (vaporized oil) upwinded by vapor potential
         Rv = state.Rv
-        f_rv = cell -> @inbounds Rv[cell]
-        rv = upwind(upw, f_rv, ∇ψ_v)
+        rv = upwind(upw, Rv, ∇ψ_v)
         # Final flux = oil phase flux + oil-in-gas flux
         q_l = rhoLS*(λb_l*∇ψ_l + rv*λb_v*∇ψ_v)
     else
@@ -41,18 +40,65 @@
     if has_disgas(sys)
         # Rs (solute gas) upwinded by liquid potential
         Rs = state.Rs
-        f_rs = cell -> @inbounds Rs[cell]
-        rs = upwind(upw, f_rs, ∇ψ_l)
+        rs = upwind(upw, Rs, ∇ψ_l)
         # Final flux = gas phase flux + gas-in-oil flux
         q_v = rhoVS*(λb_v*∇ψ_v + rs*λb_l*∇ψ_l)
     else
         q_v = rhoVS*λb_v*∇ψ_v
     end
+
+    if haskey(state, :Diffusivities)
+        S = state.Saturations
+        density = state.PhaseMassDensities
+        D = state.Diffusivities
+        if has_disgas(sys)
+            qo_diffusive_l, qo_diffusive_v = blackoil_diffusion(Rs, S, density, rhoLS, rhoVS, face, D, l, kgrad, upw)
+            q_l += qo_diffusive_l
+            q_v += qo_diffusive_v
+        end
+
+        if has_vapoil(sys)
+            qg_diffusive_v, qg_diffusive_l = blackoil_diffusion(Rv, S, density, rhoVS, rhoLS, face, D, v, kgrad, upw)
+            q_l += qg_diffusive_l
+            q_v += qg_diffusive_v
+        end
+    end
     q = setindex(q, q_l, l)
     q = setindex(q, q_v, v)
     return q
 end
 
+function blackoil_diffusion(R, S, density, rhoS_self, rhoS_dissolved, face, D, α, kgrad, upw)
+    @inbounds D_α = D[α, face]
+    X_self = cell -> black_oil_phase_mass_fraction(rhoS_self, rhoS_dissolved, R, cell)
+    # Two components: 1 - X_l - (1 - X_r) = - X_l + X_r = -(X_l - X_r) = ΔX
+    ΔX_self = -gradient(X_self, kgrad)
+    ΔX_other = -ΔX_self
+
+    T = typeof(ΔX_self)
+    mass_l = cell -> density[α, cell]*S[α, cell]
+    # TODO: Upwind or average here? Maybe doesn't matter, should be in
+    # parabolic limit for diffusion
+    # q_l += D_l*upwind(upw, mass_l, ΔX_o)*ΔX_o
+    # q_v += D_l*upwind(upw, mass_l, ΔX_g)*ΔX_g
+
+    diffused_mass = D_α*face_average(mass_l, kgrad)
+    diff_self = convert(T, diffused_mass*ΔX_self)
+    diff_dissolved = convert(T, diffused_mass*ΔX_other)
+    return (diff_self::T, diff_dissolved::T)::Tuple{T, T}
+end
+
+@inline function black_oil_phase_mass_fraction(rhoLS, rhoVS, Rs, cell)
+    # TODO: Should have molar weights here maybe, but not part of standard input
+    @inbounds rs = Rs[cell]
+    if rs < 1e-10
+        v = one(rs)
+    else
+        v = rhoLS/(rhoLS + rs*rhoVS)
+    end
+    return v
+end
+
 function apply_flow_bc!(acc, q, bc, model::StandardBlackOilModel, state, time)
     mu = state.PhaseViscosities
     b = state.ShrinkageFactors

src/blackoil/variables/varswitch.jl

@@ -222,15 +222,16 @@ end
         sw = ImmiscibleSaturation[i]
         X = BlackOilUnknown[i]
         phases = X.phases_present
+        rem = one(T) - sw + MINIMUM_COMPOSITIONAL_SATURATION
         if phases == OilOnly
             sg = zero(T)
-            so = one(T) - sw
+            so = rem
         elseif phases == GasOnly
-            sg = one(T) - sw
+            sg = rem
             so = zero(T)
         else
             sg = X.val
-            so = one(T) - sw - sg
+            so = rem - sg
         end
         s[a, i] = sw
         s[l, i] = so

src/ext.jl

@@ -40,6 +40,21 @@ function plot_reservoir_simulation_result(model::MultiModel, res::ReservoirSimRe
     return fig
 end
 
+export plot_reservoir
+function plot_reservoir(model, arg...; kwarg...)
+    rmodel = reservoir_model(model)
+    fig = plot_interactive(rmodel, arg...; kwarg...)
+    g = physical_representation(rmodel.data_domain)
+    ax = fig.current_axis[]
+    for (k, m) in pairs(model.models)
+        w = physical_representation(m.data_domain)
+        if w isa WellDomain
+            plot_well!(ax, g, w)
+        end
+    end
+    return fig
+end
+
 export simulate_reservoir_parray
 function simulate_reservoir_parray(case, mode = :mpi; kwarg...)
     sim, cfg = setup_reservoir_simulator(case; mode = mode, kwarg...)

src/facility/controls.jl

@@ -3,7 +3,11 @@ 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, key;
+        time = NaN,
+        recorder = ProgressRecorder(),
+        update_explicit = true
+    )
     # Set control to whatever is on the forces
     storage = storage_g[key]
     cfg = storage.state.WellGroupConfiguration
@@ -14,6 +18,7 @@ function Jutul.update_before_step_multimodel!(storage_g, model_g::MultiModel, mo
     for key in keys(forces.limits)
         cfg.limits[key] = forces.limits[key]
     end
+    current_step = recorder.recorder.step
     # Set operational controls
     for key in keys(forces.control)
         # If the requested control in forces differ from the one we are presently using, we need to switch.
@@ -22,7 +27,9 @@ function Jutul.update_before_step_multimodel!(storage_g, model_g::MultiModel, mo
         rstate = storage_g.Reservoir.state
         newctrl, changed = realize_control_for_reservoir(rstate, forces.control[key], rmodel, dt)
         oldctrl = req_ctrls[key]
-        if newctrl != oldctrl
+        is_new_step = cfg.step_index != current_step
+        well_was_disabled = op_ctrls[key] == DisabledControl()
+        if is_new_step && (newctrl != oldctrl || well_was_disabled)
             # We have a new control. Any previous control change is invalid.
             # Set both operating and requested control to the new one.
             @debug "Well $key switching from $oldctrl to $newctrl"
@@ -35,12 +42,13 @@ function Jutul.update_before_step_multimodel!(storage_g, model_g::MultiModel, mo
             cfg.limits[key] = merge(cfg.limits[key], as_limit(newctrl.target))
         end
     end
+    cfg.step_index = current_step
     for wname in model.domain.well_symbols
         wmodel = model_g[wname]
         wstate = storage_g[wname].state
         rmodel = model_g[:Reservoir]
         rstate = storage_g.Reservoir.state
-        update_before_step_well!(wstate, wmodel, rstate, rmodel, op_ctrls[wname])
+        update_before_step_well!(wstate, wmodel, rstate, rmodel, op_ctrls[wname], update_explicit = update_explicit)
     end
 end
 

src/facility/cross_terms.jl

@@ -27,16 +27,21 @@ function update_cross_term_in_entity!(out, i,
         well_cell = ct.well_cells[i]
         WI = state_s.WellIndices[i]
         gdz = state_s.PerforationGravityDifference[i]
+        p_well = state_s.Pressure
+        p_res = state_t.Pressure
+        dp = p_well[well_cell] - p_res[reservoir_cell]
     end
     rhoS = reference_densities(sys)
 
-    p_well = state_s.Pressure
-    p_res = state_t.Pressure
     # Wrap the key connection data in tuple for easy extension later
-    conn = (dp = p_well[well_cell] - p_res[reservoir_cell],
-           WI = WI, gdz = gdz,
-           well = well_cell,
-           reservoir = reservoir_cell)
+    conn = (
+        dp = dp,
+        WI = WI,
+        gdz = gdz,
+        well = well_cell,
+        perforation = i,
+        reservoir = reservoir_cell
+    )
     # Call smaller interface that is easy to specialize
     if haskey(state_s, :MassFractions)
         @inbounds simple_well_perforation_flux!(out, sys, state_t, state_s, rhoS, conn)
@@ -49,7 +54,7 @@ function perforation_phase_potential_difference(conn, state_res, state_well, ix)
     dp = conn.dp
     WI = conn.WI
     if haskey(state_well, :ConnectionPressureDrop)
-        dp += state_well.ConnectionPressureDrop[conn.well]
+        dp += state_well.ConnectionPressureDrop[conn.perforation]
     elseif conn.gdz != 0
         ρ_r = state_res.PhaseMassDensities[ix, conn.reservoir]
         if haskey(state_well, :PhaseMassDensities)