From d6fb9fcb71a42eede94a8151eaee130837ff5366 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Mon, 25 Sep 2023 14:45:22 +0200 Subject: [PATCH 01/50] Add option for using true well lengths for friction --- src/input_simulation/mrst_input.jl | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/input_simulation/mrst_input.jl b/src/input_simulation/mrst_input.jl index c5880035..12448247 100644 --- a/src/input_simulation/mrst_input.jl +++ b/src/input_simulation/mrst_input.jl @@ -94,6 +94,7 @@ function get_well_from_mrst_data( mrst_data, system, ix; volume = 1e-3, extraout = false, + use_lengths = false, well_type = :ms, W_data = mrst_data["W"], kwarg... @@ -195,12 +196,17 @@ function get_well_from_mrst_data( else pvol, accumulator_volume, perf_cells, well_topo, z, dz, reservoir_cells = simple_ms_setup(n, volume, well_cell_volume, rc, ref_depth, z_res) end + if use_lengths + L_i = nothing + else + L_i = 1.0 + end W = MultiSegmentWell(rc, pvol, centers, WI = WI, reference_depth = ref_depth, dz = dz, N = well_topo, name = Symbol(nm), segment_models = segment_models, - segment_length = 1.0, + segment_length = L_i, perforation_cells = perf_cells, accumulator_volume = accumulator_volume, surface_conditions = cond) @@ -783,6 +789,7 @@ function setup_case_from_mrst(casename; wells = :ms, backend = :csc, block_backend = true, split_wells = false, + use_well_lengths = false, facility_grouping = :onegroup, minbatch = 1000, steps = :full, @@ -847,7 +854,7 @@ function setup_case_from_mrst(casename; wells = :ms, for i = 1:num_wells sym = well_symbols[i] wi, wdata, res_cells = get_well_from_mrst_data(mrst_data, sys, i, W_data = first_well_set, - extraout = true, well_type = wells, context = w_context) + extraout = true, well_type = wells, context = w_context, use_lengths = use_well_lengths) param_w = setup_parameters(wi) sv = wi.secondary_variables From d4d03be99d629f365930e78191dfcc003b7dff66 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Mon, 25 Sep 2023 14:45:29 +0200 Subject: [PATCH 02/50] Add plotting function --- src/ext.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/ext.jl b/src/ext.jl index 0204c8f5..01ba5a7d 100644 --- a/src/ext.jl +++ b/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...) From dfa794d9871848ace28c2ca9ebc28fe83e3b7eaa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Mon, 25 Sep 2023 17:00:42 +0200 Subject: [PATCH 03/50] Update compositional increments to account for inactive cells --- src/multicomponent/variables/primary.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/multicomponent/variables/primary.jl b/src/multicomponent/variables/primary.jl index 964d34c0..d3546c9b 100644 --- a/src/multicomponent/variables/primary.jl +++ b/src/multicomponent/variables/primary.jl @@ -40,17 +40,19 @@ function Jutul.increment_norm(dX, state, model, X, pvar::OverallMoleFractions) else sw = missing end + M = global_map(model.domain) + active = active_entities(model.domain, M, Cells()) get_scaling(::Missing, i) = 1.0 get_scaling(s, i) = 1.0 - value(s[i]) T = eltype(dX) scale = @something Jutul.variable_scale(pvar) one(T) max_v = sum_v = max_v_scaled = sum_v_scaled = zero(T) N = degrees_of_freedom_per_entity(model, pvar) - dx_mat = reshape(dX, N, length(dX) ÷ N) - for i in axes(dx_mat, 2) - for inner in axes(dx_mat, 1) - dx = dx_mat[inner, i] - s = get_scaling(sw, i) + for i in axes(dX, 1) + for j in axes(dX, 2) + cell = active[i] + dx = dX[i, j] + s = get_scaling(sw, cell) dx_abs = abs(dx) # Scale by 1-water saturation From 4534f38528415d10cd14994bc35ce2eac8b34008 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 28 Sep 2023 10:50:41 +0200 Subject: [PATCH 04/50] Rewrite compositional criterion Fix issue of no scaling with immiscible weight, cleanup --- src/multicomponent/multicomponent.jl | 32 +++++++++++++++++----------- 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/src/multicomponent/multicomponent.jl b/src/multicomponent/multicomponent.jl index b9a9998d..ca472593 100644 --- a/src/multicomponent/multicomponent.jl +++ b/src/multicomponent/multicomponent.jl @@ -129,20 +129,28 @@ function immiscible_increment(model, state, update_report) end function compositional_residual_scale(cell, dt, w, sl, liquid_density, sv, vapor_density, sw, water_density, vol) - function sat_scale(::Nothing) - return 1.0 - end - function sat_scale(sw) + if isnothing(sw) + sw_i = 0.0 + else sw_i = sw[cell] - if sw_i > 1.0 - 1e-4 - scale = 0.0 - else - scale = 1.0 - sw_i - end - return scale end - total_density = liquid_density[cell] * sl[cell] + vapor_density[cell] * sv[cell] - return dt * mean(w) * (sat_scale(sw) / (vol[cell] * max(total_density, 1e-3))) + + if sw_i > 1.0 - 1e-4 + scale = 0.0 + else + # The convergence criterion is taken to be dimensionless in a similar + # manner to the standard CNV type criterion. We scale everything by the + # immiscible phase so that the convergence criterion is more relaxed in + # cells that are close to or completely filled with the immiscible + # phase. + scale_lv = 1.0 - sw_i + sl_scaled = sl[cell]/scale_lv + 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-3))) + end + return scale end function compositional_criterion(state, dt, active, r, nc, w, sl, liquid_density, sv, vapor_density, sw, water_density, vol) From 694488a0e0e316057e3eedeefdc31def5adc3a4f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 28 Sep 2023 20:06:06 +0200 Subject: [PATCH 05/50] Fix simple wells --- src/facility/wells/stdwells.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/facility/wells/stdwells.jl b/src/facility/wells/stdwells.jl index f4534f98..c98b37a1 100644 --- a/src/facility/wells/stdwells.jl +++ b/src/facility/wells/stdwells.jl @@ -69,6 +69,13 @@ function Jutul.initialize_extra_state_fields!(state, d::DiscretizedDomain, m::Si end end +function Jutul.select_minimum_output_variables!(vars, domain::DiscretizedDomain, model::SimpleWellFlowModel) + push!(vars, :PhaseMassDensities) + push!(vars, :Saturations) + push!(vars, :SurfaceWellConditions) + return vars +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 From 364cbf980c64d3733a8b3e2d0ea3819388bc81d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 28 Sep 2023 20:33:41 +0200 Subject: [PATCH 06/50] dp drop: Handle inactive wells --- src/facility/wells/stdwells.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/facility/wells/stdwells.jl b/src/facility/wells/stdwells.jl index c98b37a1..7602ae88 100644 --- a/src/facility/wells/stdwells.jl +++ b/src/facility/wells/stdwells.jl @@ -150,7 +150,12 @@ function update_connection_pressure_drop!(dp, well_state, well_model, res_state, end current_weight += local_weight current_density += local_density - dp[i] = current_density/current_weight + if abs(current_weight) > 0.0 + next_dp = current_density/current_weight + else + next_dp = 0.0 + end + dp[i] = next_dp end # Integrate down, using the mixture densities (temporarily stored in dp) to # calculate the pressure drop from the top. From eeb1d7b6e3e7e98dde5fe1244287c2b7e2ecb98f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 28 Sep 2023 20:40:40 +0200 Subject: [PATCH 07/50] Fix bug in connection pressure drop for simple well --- src/facility/wells/stdwells.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/facility/wells/stdwells.jl b/src/facility/wells/stdwells.jl index 7602ae88..606e7d5e 100644 --- a/src/facility/wells/stdwells.jl +++ b/src/facility/wells/stdwells.jl @@ -167,7 +167,7 @@ function update_connection_pressure_drop!(dp, well_state, well_model, res_state, Δgdz = gdz_next - gdz_current dp_current += local_density*Δgdz - gdz_current += gdz_next + gdz_current = gdz_next dp[i] = dp_current end From 25f2a3cef381a8b2639ea333abcaa404bd3d4673 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 29 Sep 2023 15:04:18 +0200 Subject: [PATCH 08/50] Simple well dp is enumerated by perf, not well cell --- src/facility/cross_terms.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/facility/cross_terms.jl b/src/facility/cross_terms.jl index c7cead4d..30fa61bb 100644 --- a/src/facility/cross_terms.jl +++ b/src/facility/cross_terms.jl @@ -36,6 +36,7 @@ function update_cross_term_in_entity!(out, i, conn = (dp = p_well[well_cell] - p_res[reservoir_cell], 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) @@ -49,7 +50,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) From 1251be136cb05149bb7476cd1c7a0c2cd869d684 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 29 Sep 2023 15:25:08 +0200 Subject: [PATCH 09/50] Use WI from state for cdp --- src/facility/wells/stdwells.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/facility/wells/stdwells.jl b/src/facility/wells/stdwells.jl index 606e7d5e..1deba1c7 100644 --- a/src/facility/wells/stdwells.jl +++ b/src/facility/wells/stdwells.jl @@ -122,14 +122,11 @@ function update_connection_pressure_drop!(dp, well_state, well_model, res_state, # accumulate the actual pressure drop due to hydrostatic assumptions. perf = physical_representation(well_model).perforations res_cells = perf.reservoir - WI = perf.WI - gdz = perf.gdz - + # Explicit update, take value. + WI = as_value(well_state.WellIndices) ρ = as_value(res_state.PhaseMassDensities) mob = as_value(res_state.PhaseMobilities) - # kr = as_value(res_state.RelativePermeabilities) - # mu = as_value(res_state.PhaseViscosities) # Integrate up, adding weighted density into well bore and keeping track of # current weight From 42c12941dc4ce88cc4896870443fba6f960ba12c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 3 Oct 2023 12:29:54 +0200 Subject: [PATCH 10/50] Alternate form of saturations for water compositional --- src/multicomponent/variables/saturations.jl | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/multicomponent/variables/saturations.jl b/src/multicomponent/variables/saturations.jl index cb3b5879..9e6bcfb3 100644 --- a/src/multicomponent/variables/saturations.jl +++ b/src/multicomponent/variables/saturations.jl @@ -15,10 +15,20 @@ end eos = model.system.equation_of_state @inbounds for i in ix S_other = ImmiscibleSaturation[i] + fr = FlashResults[i] rem = one(T) - S_other - S_l, S_v = phase_saturations(eos, Pressure[i], Temperature[i], FlashResults[i]) - Sat[l, i] = rem*S_l - Sat[v, i] = rem*S_v + if fr.state == MultiComponentFlash.two_phase_lv + S_v_pure = MultiComponentFlash.two_phase_vapor_saturation(eos, Pressure[i], Temperature[i], fr) + S_v = rem*S_v_pure + elseif fr.state == MultiComponentFlash.single_phase_v + S_v = rem + else + S_v = zero(T) + end + # S_l, S_v = phase_saturations(eos, Pressure[i], Temperature[i], fr) + S_l = one(T) - S_v - S_other + Sat[l, i] = S_l + Sat[v, i] = S_v Sat[a, i] = S_other end end From 40d867eb3a9eb86b80ee5b3883bbf30f3a2f7f63 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 3 Oct 2023 12:30:20 +0200 Subject: [PATCH 11/50] Increase minimum saturation --- src/multicomponent/multicomponent.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/multicomponent/multicomponent.jl b/src/multicomponent/multicomponent.jl index ca472593..30ab982a 100644 --- a/src/multicomponent/multicomponent.jl +++ b/src/multicomponent/multicomponent.jl @@ -2,7 +2,7 @@ using MultiComponentFlash export MultiPhaseCompositionalSystemLV export StandardVolumeSource, VolumeSource, MassSource -const MINIMUM_COMPOSITIONAL_SATURATION = 1e-10 +const MINIMUM_COMPOSITIONAL_SATURATION = 1e-5 include("variables/variables.jl") include("utils.jl") From 7ff9b8773042e98ef3250e24da2fb44ecafbce4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 3 Oct 2023 12:30:31 +0200 Subject: [PATCH 12/50] Disable parray precompilation --- .../JutulDarcyPartitionedArraysExt.jl | 48 +++++++++---------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl b/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl index 25396cd4..2ae5bb75 100644 --- a/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl +++ b/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 From 95f5ff723fa3f6cf552cf82075600fbcd662a702 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 3 Oct 2023 15:31:26 +0200 Subject: [PATCH 13/50] Fix simple well constructor --- src/facility/wells/wells.jl | 2 +- src/utils.jl | 7 +------ 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/src/facility/wells/wells.jl b/src/facility/wells/wells.jl index 3462893d..2bb7851c 100644 --- a/src/facility/wells/wells.jl +++ b/src/facility/wells/wells.jl @@ -125,7 +125,7 @@ function setup_well(g, K, reservoir_cells::AbstractVector; volumes[i] = h*π*r_i^2 end if simple_well - W = SimpleWell(reservoir_cells, WI = WI_computed, dz = dz, reference_depth = reference_depth, kwarg...) + W = SimpleWell(reservoir_cells; WI = WI_computed, dz = dz, kwarg...) else # Depth differences are taken care of via centers. dz *= 0.0 diff --git a/src/utils.jl b/src/utils.jl index 551a5df3..1bad3782 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -63,14 +63,9 @@ function setup_reservoir_model(reservoir::DataDomain, system; mode = PredictionMode() if length(wells) > 0 for w in wells - if w isa SimpleWell - well_context = reservoir_context - else - well_context = context - end w_domain = DataDomain(w) wname = w.name - models[wname] = SimulationModel(w_domain, system, context = well_context) + models[wname] = SimulationModel(w_domain, system, context = context) if split_wells wg = WellGroup([wname]) F = SimulationModel(wg, mode, context = context, data_domain = DataDomain(wg)) From 6ac10bfa27763886a98badbd3d65a8082ebd9256 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 3 Oct 2023 15:31:41 +0200 Subject: [PATCH 14/50] Clean up perforation cross term --- src/facility/cross_terms.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/facility/cross_terms.jl b/src/facility/cross_terms.jl index 30fa61bb..2bac7d8c 100644 --- a/src/facility/cross_terms.jl +++ b/src/facility/cross_terms.jl @@ -27,17 +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, - perforation = i, - 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) From 0c64adecfa9759852e5b744809238644d4fa8da1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 4 Oct 2023 09:33:37 +0200 Subject: [PATCH 15/50] Work on skipping flash for pure water cells --- src/input_simulation/mrst_input.jl | 3 ++- src/multicomponent/multicomponent.jl | 4 ++++ src/multicomponent/variables/flash.jl | 27 ++++++++++++++++++++----- src/multicomponent/variables/primary.jl | 2 +- src/multicomponent/wells.jl | 2 +- 5 files changed, 30 insertions(+), 8 deletions(-) diff --git a/src/input_simulation/mrst_input.jl b/src/input_simulation/mrst_input.jl index 12448247..18d655d2 100644 --- a/src/input_simulation/mrst_input.jl +++ b/src/input_simulation/mrst_input.jl @@ -728,7 +728,8 @@ function init_from_mat(mrst_data, model, param) s = copy(state0["s"]) if size(s, 2) == 3 sw = vec(s[:, 1]) - sw = min.(sw, 1 - MINIMUM_COMPOSITIONAL_SATURATION) + # sw = min.(sw, 1 - MINIMUM_COMPOSITIONAL_SATURATION) + sw = min.(sw, 1.0) init[:ImmiscibleSaturation] = sw else @assert size(s, 2) == 2 diff --git a/src/multicomponent/multicomponent.jl b/src/multicomponent/multicomponent.jl index 30ab982a..f168a8c7 100644 --- a/src/multicomponent/multicomponent.jl +++ b/src/multicomponent/multicomponent.jl @@ -4,6 +4,10 @@ export StandardVolumeSource, VolumeSource, MassSource const MINIMUM_COMPOSITIONAL_SATURATION = 1e-5 +@inline function is_pure_single_phase(s_immiscible) + return s_immiscible > 1.0 - MINIMUM_COMPOSITIONAL_SATURATION +end + include("variables/variables.jl") include("utils.jl") include("flux.jl") diff --git a/src/multicomponent/variables/flash.jl b/src/multicomponent/variables/flash.jl index 7273d678..f292fee6 100644 --- a/src/multicomponent/variables/flash.jl +++ b/src/multicomponent/variables/flash.jl @@ -75,6 +75,10 @@ function initialize_variable_ad!(state, model, pvar::FlashResults, symb, npartia end @jutul_secondary function update_flash!(flash_results, fr::FlashResults, model, Pressure, Temperature, OverallMoleFractions, ix) + flash_entity_loop!(flash_results, fr, model, Pressure, Temperature, OverallMoleFractions, nothing, ix) +end + +function flash_entity_loop!(flash_results, fr, model, Pressure, Temperature, OverallMoleFractions, sw, ix) storage, m, buffers = fr.storage, fr.method, fr.update_buffer eos = model.system.equation_of_state @@ -83,7 +87,7 @@ end buf_z = buf.z buf_forces = buf.forces @inbounds for i in ix - flash_results[i] = internal_flash!(flash_results[i], S, m, eos, buf_z, buf_forces, Pressure, Temperature, OverallMoleFractions, i) + flash_results[i] = internal_flash!(flash_results[i], S, m, eos, buf_z, buf_forces, Pressure, Temperature, OverallMoleFractions, sw, i) end end @@ -128,27 +132,40 @@ function update_flash_buffer!(buf, eos, Pressure, Temperature, OverallMoleFracti end end -function internal_flash!(f, S, m, eos, buf_z, buf_forces, Pressure, Temperature, OverallMoleFractions, i) +function internal_flash!(f, S, m, eos, buf_z, buf_forces, Pressure, Temperature, OverallMoleFractions, sw, i) + @inline function immiscible_sat(::Nothing, i) + return 0.0 + end + + @inline function immiscible_sat(s, i) + return @inbounds s[i] + end + @inbounds begin P = Pressure[i] T = Temperature[i] Z = @view OverallMoleFractions[:, i] + Sw = immiscible_sat(sw, i) K = f.K x = f.liquid.mole_fractions y = f.vapor.mole_fractions - return update_flash_result(S, m, eos, K, x, y, buf_z, buf_forces, P, T, Z) + return update_flash_result(S, m, eos, K, x, y, buf_z, buf_forces, P, T, Z, Sw) end end -function update_flash_result(S, m, eos, K, x, y, z, forces, P, T, Z) +function update_flash_result(S, m, eos, K, x, y, z, forces, P, T, Z, Sw = 0.0) @. z = max(value(Z), 1e-8) # Conditions c = (p = value(P), T = value(T), z = z) # Perform flash - vapor_frac = flash_2ph!(S, K, eos, c, NaN, method = m, extra_out = false, z_min = nothing) + if is_pure_single_phase(Sw) + vapor_frac = NaN + else + vapor_frac = flash_2ph!(S, K, eos, c, NaN, method = m, extra_out = false, z_min = nothing) + end force_coefficients!(forces, eos, (p = P, T = T, z = Z)) if isnan(vapor_frac) # Single phase condition. Life is easy. diff --git a/src/multicomponent/variables/primary.jl b/src/multicomponent/variables/primary.jl index d3546c9b..53d74701 100644 --- a/src/multicomponent/variables/primary.jl +++ b/src/multicomponent/variables/primary.jl @@ -75,6 +75,6 @@ Base.@kwdef struct ImmiscibleSaturation <: ScalarVariable ds_max::Float64 = 0.2 end -maximum_value(::ImmiscibleSaturation) = 1.0 - MINIMUM_COMPOSITIONAL_SATURATION +maximum_value(::ImmiscibleSaturation) = 1.0# - MINIMUM_COMPOSITIONAL_SATURATION minimum_value(::ImmiscibleSaturation) = 0.0 absolute_increment_limit(s::ImmiscibleSaturation) = s.ds_max diff --git a/src/multicomponent/wells.jl b/src/multicomponent/wells.jl index e7024ac5..88fd0a38 100644 --- a/src/multicomponent/wells.jl +++ b/src/multicomponent/wells.jl @@ -154,7 +154,7 @@ function separator_flash!(flash, eos, cond, z) Temperature = cond.T update_flash_buffer!(buf, eos, Pressure, Temperature, z) forces = buf.forces - return update_flash_result(fstorage, SSIFlash(), eos, f.K, x, y, buf.z, forces, Pressure, Temperature, z) + return update_flash_result(fstorage, SSIFlash(), eos, f.K, x, y, buf.z, forces, Pressure, Temperature, z, 0.0) end function flash_stream!(moles::SVector{N, T}, flash, eos, cond) where {N, T} From fcf2ec1b884b5730e3636f5bea9643cccbd82407 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 4 Oct 2023 09:51:48 +0200 Subject: [PATCH 16/50] Split up a function --- src/multicomponent/variables/others.jl | 39 +++++++++++++++++--------- 1 file changed, 26 insertions(+), 13 deletions(-) diff --git a/src/multicomponent/variables/others.jl b/src/multicomponent/variables/others.jl index d848bdac..436a1821 100644 --- a/src/multicomponent/variables/others.jl +++ b/src/multicomponent/variables/others.jl @@ -27,19 +27,31 @@ end end # Total masses -@jutul_secondary function update_total_masses!(totmass, tv::TotalMasses, model::SimulationModel{G,S}, - FlashResults, - PhaseMassDensities, - Saturations, - VaporMassFractions, - LiquidMassFractions, - FluidVolume, ix) where {G,S<:CompositionalSystem} - pv = FluidVolume - ρ = PhaseMassDensities - X = LiquidMassFractions - Y = VaporMassFractions - Sat = Saturations - F = FlashResults +@jutul_secondary function update_total_masses!(totmass, + tv::TotalMasses, + model::SimulationModel{G,S}, + FlashResults, + PhaseMassDensities, + Saturations, + VaporMassFractions, + LiquidMassFractions, + FluidVolume, + ix + ) where {G,S<:CompositionalSystem} + compositiona_mass_update_loop!( + totmass, + model, + FlashResults, + PhaseMassDensities, + Saturations, + LiquidMassFractions, + VaporMassFractions, + FluidVolume, + ix + ) +end + +function compositiona_mass_update_loop!(totmass, model, F, ρ, Sat, X, Y, pv, ix) sys = model.system phase_ix = phase_indices(sys) has_other = Val(has_other_phase(sys)) @@ -47,6 +59,7 @@ end for cell in ix @inbounds two_phase_compositional_mass!(totmass, F[cell].state, pv, ρ, X, Y, Sat, cell, N, has_other, phase_ix) end + return totmass end function degrees_of_freedom_per_entity(model::SimulationModel{G,S}, v::TotalMasses) where {G<:Any,S<:MultiComponentSystem} From e5b71da806061f19c01286c80cd50a1ab436e422 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 4 Oct 2023 11:00:21 +0200 Subject: [PATCH 17/50] Add some helpful consts for compositional --- src/multiphase.jl | 3 --- src/types.jl | 9 ++++++++- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/multiphase.jl b/src/multiphase.jl index db05ac6b..0abfb277 100644 --- a/src/multiphase.jl +++ b/src/multiphase.jl @@ -40,9 +40,6 @@ number_of_phases(::SinglePhaseSystem) = 1 number_of_phases(sys::CompositeSystem) = number_of_phases(sys.systems.flow) ## Phases -# Abstract phase -abstract type AbstractPhase end - function get_short_name(phase::AbstractPhase) return get_name(phase)[] end diff --git a/src/types.jl b/src/types.jl index e6d782fe..d018b38c 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,3 +1,4 @@ +abstract type AbstractPhase end abstract type MultiPhaseSystem <: JutulSystem end abstract type MultiComponentSystem <: MultiPhaseSystem end @@ -17,7 +18,13 @@ struct MultiPhaseCompositionalSystemLV{E, T, O, R} <: CompositionalSystem where equation_of_state::E rho_ref::R end -const LVCompositionalModel = SimulationModel{D, S, F, C} where {D, S<:MultiPhaseCompositionalSystemLV{<:Any, <:Any, <:Any}, F, C} + +const LVCompositional2PhaseSystem = MultiPhaseCompositionalSystemLV{<:Any, <:Any, Nothing, <:Any} +const LVCompositional3PhaseSystem = MultiPhaseCompositionalSystemLV{<:Any, <:Any, <:AbstractPhase, <:Any} + +const LVCompositionalModel = SimulationModel{D, S, F, C} where {D, S<:MultiPhaseCompositionalSystemLV{<:Any, <:Any, <:Any, <:Any}, F, C} +const LVCompositionalModel2Phase = SimulationModel{D, S, F, C} where {D, S<:LVCompositional2PhaseSystem, F, C} +const LVCompositionalModel3Phase = SimulationModel{D, S, F, C} where {D, S<:LVCompositional3PhaseSystem, F, C} """ MultiPhaseCompositionalSystemLV(equation_of_state, phases = (LiquidPhase(), VaporPhase()); reference_densities = ones(length(phases)), other_name = "Water") From abc847f834571daf1ab0c5eb5d0424e800712e85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 5 Oct 2023 21:13:54 +0200 Subject: [PATCH 18/50] Experimental water changes --- src/multicomponent/variables/flash.jl | 8 ++- src/multicomponent/variables/others.jl | 78 +++++++++++++++------ src/multicomponent/variables/saturations.jl | 2 +- 3 files changed, 65 insertions(+), 23 deletions(-) diff --git a/src/multicomponent/variables/flash.jl b/src/multicomponent/variables/flash.jl index f292fee6..a1f1456f 100644 --- a/src/multicomponent/variables/flash.jl +++ b/src/multicomponent/variables/flash.jl @@ -78,6 +78,10 @@ end flash_entity_loop!(flash_results, fr, model, Pressure, Temperature, OverallMoleFractions, nothing, ix) end +@jutul_secondary function update_flash!(flash_results, fr::FlashResults, model::LVCompositionalModel3Phase, Pressure, Temperature, OverallMoleFractions, ImmiscibleSaturation, ix) + flash_entity_loop!(flash_results, fr, model, Pressure, Temperature, OverallMoleFractions, ImmiscibleSaturation, ix) +end + function flash_entity_loop!(flash_results, fr, model, Pressure, Temperature, OverallMoleFractions, sw, ix) storage, m, buffers = fr.storage, fr.method, fr.update_buffer eos = model.system.equation_of_state @@ -161,7 +165,9 @@ function update_flash_result(S, m, eos, K, x, y, z, forces, P, T, Z, Sw = 0.0) # Conditions c = (p = value(P), T = value(T), z = z) # Perform flash - if is_pure_single_phase(Sw) + @info "Cell water = $(value(Sw))" + if is_pure_single_phase(Sw) && false + @info "Setting to single phase: Sw = $Sw" vapor_frac = NaN else vapor_frac = flash_2ph!(S, K, eos, c, NaN, method = m, extra_out = false, z_min = nothing) diff --git a/src/multicomponent/variables/others.jl b/src/multicomponent/variables/others.jl index 436a1821..4a35a493 100644 --- a/src/multicomponent/variables/others.jl +++ b/src/multicomponent/variables/others.jl @@ -27,18 +27,19 @@ end end # Total masses -@jutul_secondary function update_total_masses!(totmass, - tv::TotalMasses, - model::SimulationModel{G,S}, - FlashResults, - PhaseMassDensities, - Saturations, - VaporMassFractions, - LiquidMassFractions, - FluidVolume, - ix - ) where {G,S<:CompositionalSystem} - compositiona_mass_update_loop!( +@jutul_secondary function update_total_masses!( + totmass, + tmvar::TotalMasses, + model::LVCompositionalModel2Phase, + FlashResults, + PhaseMassDensities, + Saturations, + VaporMassFractions, + LiquidMassFractions, + FluidVolume, + ix + ) + compositional_mass_update_loop!( totmass, model, FlashResults, @@ -46,18 +47,46 @@ end Saturations, LiquidMassFractions, VaporMassFractions, + nothing, FluidVolume, ix ) end -function compositiona_mass_update_loop!(totmass, model, F, ρ, Sat, X, Y, pv, ix) +@jutul_secondary function update_total_masses!( + totmass, + tmvar::TotalMasses, + model::LVCompositionalModel3Phase, + FlashResults, + PhaseMassDensities, + Saturations, + VaporMassFractions, + LiquidMassFractions, + ImmiscibleSaturation, + FluidVolume, + ix + ) + compositional_mass_update_loop!( + totmass, + model, + FlashResults, + PhaseMassDensities, + Saturations, + LiquidMassFractions, + VaporMassFractions, + ImmiscibleSaturation, + FluidVolume, + ix + ) +end + + +function compositional_mass_update_loop!(totmass, model, F, ρ, Sat, X, Y, sw, pv, ix) sys = model.system phase_ix = phase_indices(sys) - has_other = Val(has_other_phase(sys)) N = size(totmass, 1) for cell in ix - @inbounds two_phase_compositional_mass!(totmass, F[cell].state, pv, ρ, X, Y, Sat, cell, N, has_other, phase_ix) + @inbounds two_phase_compositional_mass!(totmass, F[cell].state, sw, pv, ρ, X, Y, Sat, cell, N, phase_ix) end return totmass end @@ -69,20 +98,20 @@ end """ Update total masses for two-phase compositional """ -function two_phase_compositional_mass!(M, state, Φ, ρ, X, Y, S, cell, N, aqua::Val{false}, phase_ix) - update_mass_two_phase_compositional!(M, state, Φ, ρ, X, Y, S, cell, phase_ix, N) +function two_phase_compositional_mass!(M, state, sw::Nothing, Φ, ρ, X, Y, S, cell, N, phase_ix) + update_mass_two_phase_compositional!(M, state, sw, Φ, ρ, X, Y, S, cell, phase_ix, N) end """ Update total masses for two-phase compositional where another immiscible phase is present """ -function two_phase_compositional_mass!(M, state, Φ, ρ, X, Y, S, cell, N, aqua::Val{true}, phase_ix) - update_mass_two_phase_compositional!(M, state, Φ, ρ, X, Y, S, cell, phase_ix[2:end], N - 1) +function two_phase_compositional_mass!(M, state, sw, Φ, ρ, X, Y, S, cell, N, phase_ix) + update_mass_two_phase_compositional!(M, state, sw, Φ, ρ, X, Y, S, cell, phase_ix[2:end], N - 1) a, = phase_ix @inbounds M[end, cell] = ρ[a, cell]*S[a, cell]*Φ[cell] end -function update_mass_two_phase_compositional!(M, state, Φ, ρ, X, Y, S, cell, phase_ix, N) +function update_mass_two_phase_compositional!(M, state, sw, Φ, ρ, X, Y, S, cell, phase_ix, N) l, v = phase_ix has_liquid = liquid_phase_present(state) has_vapor = vapor_phase_present(state) @@ -96,7 +125,14 @@ function update_mass_two_phase_compositional!(M, state, Φ, ρ, X, Y, S, cell, p end function single_phase_mass!(M, ρ, S, mass_fractions, Φ, cell, N, phase) - @inbounds M_l = ρ[phase, cell] * S[phase, cell] + S_eos = S[phase, cell] + if S_eos < MINIMUM_COMPOSITIONAL_SATURATION + S_eos = replace_value(S_eos, MINIMUM_COMPOSITIONAL_SATURATION) + end + # S_eos = max(S[phase, cell], MINIMUM_COMPOSITIONAL_SATURATION) + @info "?? $phase" S[phase, cell] MINIMUM_COMPOSITIONAL_SATURATION S_eos + + @inbounds M_l = ρ[phase, cell] * S_eos for c in 1:N @inbounds M[c, cell] = M_l*mass_fractions[c, cell]*Φ[cell] end diff --git a/src/multicomponent/variables/saturations.jl b/src/multicomponent/variables/saturations.jl index 9e6bcfb3..2411c093 100644 --- a/src/multicomponent/variables/saturations.jl +++ b/src/multicomponent/variables/saturations.jl @@ -16,7 +16,7 @@ end @inbounds for i in ix S_other = ImmiscibleSaturation[i] fr = FlashResults[i] - rem = one(T) - S_other + rem = one(T) - S_other - MINIMUM_COMPOSITIONAL_SATURATION if fr.state == MultiComponentFlash.two_phase_lv S_v_pure = MultiComponentFlash.two_phase_vapor_saturation(eos, Pressure[i], Temperature[i], fr) S_v = rem*S_v_pure From 0a7dedbbd85913fe133b31f13719998d1962fb13 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 6 Oct 2023 12:07:38 +0200 Subject: [PATCH 19/50] Restore old behavior for test --- src/multicomponent/variables/flash.jl | 2 -- src/multicomponent/variables/others.jl | 4 ++-- src/multicomponent/variables/primary.jl | 2 +- src/multicomponent/variables/saturations.jl | 5 ++--- 4 files changed, 5 insertions(+), 8 deletions(-) diff --git a/src/multicomponent/variables/flash.jl b/src/multicomponent/variables/flash.jl index a1f1456f..abc7d989 100644 --- a/src/multicomponent/variables/flash.jl +++ b/src/multicomponent/variables/flash.jl @@ -165,9 +165,7 @@ function update_flash_result(S, m, eos, K, x, y, z, forces, P, T, Z, Sw = 0.0) # Conditions c = (p = value(P), T = value(T), z = z) # Perform flash - @info "Cell water = $(value(Sw))" if is_pure_single_phase(Sw) && false - @info "Setting to single phase: Sw = $Sw" vapor_frac = NaN else vapor_frac = flash_2ph!(S, K, eos, c, NaN, method = m, extra_out = false, z_min = nothing) diff --git a/src/multicomponent/variables/others.jl b/src/multicomponent/variables/others.jl index 4a35a493..d59f58f4 100644 --- a/src/multicomponent/variables/others.jl +++ b/src/multicomponent/variables/others.jl @@ -127,10 +127,10 @@ end function single_phase_mass!(M, ρ, S, mass_fractions, Φ, cell, N, phase) S_eos = S[phase, cell] if S_eos < MINIMUM_COMPOSITIONAL_SATURATION - S_eos = replace_value(S_eos, MINIMUM_COMPOSITIONAL_SATURATION) + # S_eos = replace_value(S_eos, MINIMUM_COMPOSITIONAL_SATURATION) end # S_eos = max(S[phase, cell], MINIMUM_COMPOSITIONAL_SATURATION) - @info "?? $phase" S[phase, cell] MINIMUM_COMPOSITIONAL_SATURATION S_eos + # @info "?? $phase" S[phase, cell] MINIMUM_COMPOSITIONAL_SATURATION S_eos @inbounds M_l = ρ[phase, cell] * S_eos for c in 1:N diff --git a/src/multicomponent/variables/primary.jl b/src/multicomponent/variables/primary.jl index 53d74701..d3546c9b 100644 --- a/src/multicomponent/variables/primary.jl +++ b/src/multicomponent/variables/primary.jl @@ -75,6 +75,6 @@ Base.@kwdef struct ImmiscibleSaturation <: ScalarVariable ds_max::Float64 = 0.2 end -maximum_value(::ImmiscibleSaturation) = 1.0# - MINIMUM_COMPOSITIONAL_SATURATION +maximum_value(::ImmiscibleSaturation) = 1.0 - MINIMUM_COMPOSITIONAL_SATURATION minimum_value(::ImmiscibleSaturation) = 0.0 absolute_increment_limit(s::ImmiscibleSaturation) = s.ds_max diff --git a/src/multicomponent/variables/saturations.jl b/src/multicomponent/variables/saturations.jl index 2411c093..d84d35eb 100644 --- a/src/multicomponent/variables/saturations.jl +++ b/src/multicomponent/variables/saturations.jl @@ -16,7 +16,7 @@ end @inbounds for i in ix S_other = ImmiscibleSaturation[i] fr = FlashResults[i] - rem = one(T) - S_other - MINIMUM_COMPOSITIONAL_SATURATION + rem = one(T) - S_other + 0*MINIMUM_COMPOSITIONAL_SATURATION if fr.state == MultiComponentFlash.two_phase_lv S_v_pure = MultiComponentFlash.two_phase_vapor_saturation(eos, Pressure[i], Temperature[i], fr) S_v = rem*S_v_pure @@ -25,8 +25,7 @@ end else S_v = zero(T) end - # S_l, S_v = phase_saturations(eos, Pressure[i], Temperature[i], fr) - S_l = one(T) - S_v - S_other + S_l = rem - S_v Sat[l, i] = S_l Sat[v, i] = S_v Sat[a, i] = S_other From cddeafe4d7003a915e25a534adc46ef16c18be29 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 6 Oct 2023 12:27:16 +0200 Subject: [PATCH 20/50] Allow saturation to go to 0 --- src/multicomponent/variables/flash.jl | 2 +- src/multicomponent/variables/others.jl | 2 +- src/multicomponent/variables/primary.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/multicomponent/variables/flash.jl b/src/multicomponent/variables/flash.jl index abc7d989..4d8e132e 100644 --- a/src/multicomponent/variables/flash.jl +++ b/src/multicomponent/variables/flash.jl @@ -165,7 +165,7 @@ function update_flash_result(S, m, eos, K, x, y, z, forces, P, T, Z, Sw = 0.0) # Conditions c = (p = value(P), T = value(T), z = z) # Perform flash - if is_pure_single_phase(Sw) && false + if is_pure_single_phase(Sw) vapor_frac = NaN else vapor_frac = flash_2ph!(S, K, eos, c, NaN, method = m, extra_out = false, z_min = nothing) diff --git a/src/multicomponent/variables/others.jl b/src/multicomponent/variables/others.jl index d59f58f4..aae26959 100644 --- a/src/multicomponent/variables/others.jl +++ b/src/multicomponent/variables/others.jl @@ -127,7 +127,7 @@ end function single_phase_mass!(M, ρ, S, mass_fractions, Φ, cell, N, phase) S_eos = S[phase, cell] if S_eos < MINIMUM_COMPOSITIONAL_SATURATION - # S_eos = replace_value(S_eos, MINIMUM_COMPOSITIONAL_SATURATION) + S_eos = replace_value(S_eos, MINIMUM_COMPOSITIONAL_SATURATION) end # S_eos = max(S[phase, cell], MINIMUM_COMPOSITIONAL_SATURATION) # @info "?? $phase" S[phase, cell] MINIMUM_COMPOSITIONAL_SATURATION S_eos diff --git a/src/multicomponent/variables/primary.jl b/src/multicomponent/variables/primary.jl index d3546c9b..53d74701 100644 --- a/src/multicomponent/variables/primary.jl +++ b/src/multicomponent/variables/primary.jl @@ -75,6 +75,6 @@ Base.@kwdef struct ImmiscibleSaturation <: ScalarVariable ds_max::Float64 = 0.2 end -maximum_value(::ImmiscibleSaturation) = 1.0 - MINIMUM_COMPOSITIONAL_SATURATION +maximum_value(::ImmiscibleSaturation) = 1.0# - MINIMUM_COMPOSITIONAL_SATURATION minimum_value(::ImmiscibleSaturation) = 0.0 absolute_increment_limit(s::ImmiscibleSaturation) = s.ds_max From 9082461182bbab4e39b59a5ca553b8b5211471e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 6 Oct 2023 13:42:21 +0200 Subject: [PATCH 21/50] Experimental well disabling when unable to converge --- src/facility/wellgroups.jl | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/src/facility/wellgroups.jl b/src/facility/wellgroups.jl index a656e2c1..1ff0657b 100644 --- a/src/facility/wellgroups.jl +++ b/src/facility/wellgroups.jl @@ -20,22 +20,36 @@ function Jutul.update_primary_variable!(state, massrate::TotalSurfaceMassRate, s # producers can only have strictly negative and disabled controls give zero rate. rel_max = Jutul.absolute_increment_limit(massrate) abs_max = Jutul.relative_increment_limit(massrate) - function do_update(v, dx, ctrl) + function do_update!(wcfg, s, v, dx, ctrl) return Jutul.update_value(v, dx) end - function do_update(v, dx, ctrl::InjectorControl) - return Jutul.update_value(v, dx, abs_max, rel_max, MIN_ACTIVE_WELL_RATE, nothing) + function do_update!(wcfg, s, v, dx, ctrl::InjectorControl) + if v == MIN_ACTIVE_WELL_RATE && v + dx < MIN_ACTIVE_WELL_RATE + @info "Disabling injector $s" + wcfg.operating_controls[s] = DisabledControl() + next = Jutul.update_value(v, -value(v)) + else + next = Jutul.update_value(v, dx, abs_max, rel_max, MIN_ACTIVE_WELL_RATE, nothing) + end + return next end - function do_update(v, dx, ctrl::ProducerControl) - return Jutul.update_value(v, dx, abs_max, rel_max, nothing, -MIN_ACTIVE_WELL_RATE) + function do_update!(wcfg, s, v, dx, ctrl::ProducerControl) + if v == MIN_ACTIVE_WELL_RATE && v + dx > MIN_ACTIVE_WELL_RATE + @info "Disabling producer $s" + wcfg.operating_controls[s] = DisabledControl() + next = Jutul.update_value(v, -value(v)) + else + next = Jutul.update_value(v, dx, abs_max, rel_max, nothing, -MIN_ACTIVE_WELL_RATE) + end + return next end - function do_update(v, dx, ctrl::DisabledControl) + function do_update!(wcfg, s, v, dx, ctrl::DisabledControl) # Set value to zero since we know it is correct. return Jutul.update_value(v, -value(v)) end @inbounds for i in eachindex(v) s = symbols[i] - v[i] = do_update(v[i], w*dx[i], operating_control(cfg, s)) + v[i] = do_update!(cfg, s, v[i], w*dx[i], operating_control(cfg, s)) end end From e508566e196fb6c0974a38b2437db259d608d125 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 8 Oct 2023 21:08:54 +0200 Subject: [PATCH 22/50] Revive wells if in new step --- src/facility/controls.jl | 8 ++++++-- src/facility/cross_terms.jl | 2 +- src/facility/types.jl | 13 +++++++------ 3 files changed, 14 insertions(+), 9 deletions(-) diff --git a/src/facility/controls.jl b/src/facility/controls.jl index e2689b66..6c984716 100644 --- a/src/facility/controls.jl +++ b/src/facility/controls.jl @@ -3,7 +3,7 @@ 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; recorder, kwarg...) # Set control to whatever is on the forces storage = storage_g[key] cfg = storage.state.WellGroupConfiguration @@ -14,6 +14,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 +23,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,6 +38,7 @@ 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 diff --git a/src/facility/cross_terms.jl b/src/facility/cross_terms.jl index 2bac7d8c..c405089a 100644 --- a/src/facility/cross_terms.jl +++ b/src/facility/cross_terms.jl @@ -28,7 +28,7 @@ function update_cross_term_in_entity!(out, i, WI = state_s.WellIndices[i] gdz = state_s.PerforationGravityDifference[i] p_well = state_s.Pressure - p_res = state_t.Pressure + p_res = state_t.Pressure dp = p_well[well_cell] - p_res[reservoir_cell] end rhoS = reference_densities(sys) diff --git a/src/facility/types.jl b/src/facility/types.jl index 5a8dfda1..1380fdb4 100644 --- a/src/facility/types.jl +++ b/src/facility/types.jl @@ -265,11 +265,12 @@ function replace_target(f::ProducerControl, target) return ProducerControl(target) end -struct WellGroupConfiguration - operating_controls # Currently operating control - requested_controls # The requested control (which may be different if limits are hit) - limits # Operating limits for the wells - function WellGroupConfiguration(well_symbols, control = nothing, limits = nothing) +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) + const limits::L # Operating limits for the wells + step_index::Int + function WellGroupConfiguration(well_symbols, control = nothing, limits = nothing, step = 0) if isnothing(control) control = Dict{Symbol, WellControlForce}() for s in well_symbols @@ -283,7 +284,7 @@ struct WellGroupConfiguration limits[s] = nothing end end - new(control, requested, limits) + new{typeof(control), typeof(requested), typeof(limits)}(control, requested, limits, step) end end From 2d0cc99f12758344936ce0ff807bb6b688ab1c24 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 8 Oct 2023 21:09:07 +0200 Subject: [PATCH 23/50] Reset back minim sat --- src/multicomponent/multicomponent.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/multicomponent/multicomponent.jl b/src/multicomponent/multicomponent.jl index f168a8c7..88cc99a6 100644 --- a/src/multicomponent/multicomponent.jl +++ b/src/multicomponent/multicomponent.jl @@ -2,7 +2,7 @@ using MultiComponentFlash export MultiPhaseCompositionalSystemLV export StandardVolumeSource, VolumeSource, MassSource -const MINIMUM_COMPOSITIONAL_SATURATION = 1e-5 +const MINIMUM_COMPOSITIONAL_SATURATION = 1e-10 @inline function is_pure_single_phase(s_immiscible) return s_immiscible > 1.0 - MINIMUM_COMPOSITIONAL_SATURATION From 7527cc209c8538e3dbf6c19726bb84d9400892fb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 10 Oct 2023 11:07:10 +0200 Subject: [PATCH 24/50] make disabled wells a configurable option --- src/facility/types.jl | 11 ++++++++--- src/facility/wellgroups.jl | 9 +++++---- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/facility/types.jl b/src/facility/types.jl index 1380fdb4..311689ea 100644 --- a/src/facility/types.jl +++ b/src/facility/types.jl @@ -6,8 +6,13 @@ struct HistoryMode <: FacilitySystem end abstract type SurfaceFacilityDomain <: JutulDomain end abstract type WellControllerDomain <: SurfaceFacilityDomain end -struct WellGroup <: WellControllerDomain - well_symbols::Vector{Symbol} +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 +end + +function WellGroup(wells::Vector{Symbol}; can_shut_wells = true) + return WellGroup(wells, can_shut_wells) end const WellGroupModel = SimulationModel{WellGroup, <:Any, <:Any, <:Any} @@ -269,7 +274,7 @@ 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) const limits::L # Operating limits for the wells - step_index::Int + step_index::Int # Internal book-keeping of what step we are at function WellGroupConfiguration(well_symbols, control = nothing, limits = nothing, step = 0) if isnothing(control) control = Dict{Symbol, WellControlForce}() diff --git a/src/facility/wellgroups.jl b/src/facility/wellgroups.jl index 1ff0657b..dcccce75 100644 --- a/src/facility/wellgroups.jl +++ b/src/facility/wellgroups.jl @@ -16,6 +16,7 @@ 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) @@ -24,8 +25,8 @@ function Jutul.update_primary_variable!(state, massrate::TotalSurfaceMassRate, s return Jutul.update_value(v, dx) end function do_update!(wcfg, s, v, dx, ctrl::InjectorControl) - if v == MIN_ACTIVE_WELL_RATE && v + dx < MIN_ACTIVE_WELL_RATE - @info "Disabling injector $s" + if v == MIN_ACTIVE_WELL_RATE && v + dx < MIN_ACTIVE_WELL_RATE && can_shut_wells + 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)) else @@ -34,8 +35,8 @@ function Jutul.update_primary_variable!(state, massrate::TotalSurfaceMassRate, s return next end function do_update!(wcfg, s, v, dx, ctrl::ProducerControl) - if v == MIN_ACTIVE_WELL_RATE && v + dx > MIN_ACTIVE_WELL_RATE - @info "Disabling producer $s" + if v == MIN_ACTIVE_WELL_RATE && v + dx > MIN_ACTIVE_WELL_RATE && can_shut_wells + 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)) else From e8245346b105ecce6380914d022a44fa108bc953 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 10 Oct 2023 13:15:12 +0200 Subject: [PATCH 25/50] Update varswitch.jl --- src/blackoil/variables/varswitch.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/blackoil/variables/varswitch.jl b/src/blackoil/variables/varswitch.jl index a9b9b07b..f774a166 100644 --- a/src/blackoil/variables/varswitch.jl +++ b/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 From 2560e8391d4b826e3f44e56a2ebd44c24fa87039 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 10 Oct 2023 13:15:24 +0200 Subject: [PATCH 26/50] Update multicomponent.jl --- src/multicomponent/multicomponent.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/multicomponent/multicomponent.jl b/src/multicomponent/multicomponent.jl index 88cc99a6..17eca08d 100644 --- a/src/multicomponent/multicomponent.jl +++ b/src/multicomponent/multicomponent.jl @@ -139,7 +139,7 @@ function compositional_residual_scale(cell, dt, w, sl, liquid_density, sv, vapor sw_i = sw[cell] end - if sw_i > 1.0 - 1e-4 + if sw_i > 1.0 - MINIMUM_COMPOSITIONAL_SATURATION scale = 0.0 else # The convergence criterion is taken to be dimensionless in a similar @@ -147,12 +147,12 @@ function compositional_residual_scale(cell, dt, w, sl, liquid_density, sv, vapor # immiscible phase so that the convergence criterion is more relaxed in # cells that are close to or completely filled with the immiscible # phase. - scale_lv = 1.0 - sw_i + scale_lv = 1.0 - sw_i + MINIMUM_COMPOSITIONAL_SATURATION sl_scaled = sl[cell]/scale_lv 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-3))) + scale = dt * mean(w) * (scale_lv / (vol[cell] * max(total_density, 1e-8))) end return scale end From 23b8c12cee4aa5439ef0fbc98c1415143c14e624 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 10 Oct 2023 15:00:31 +0200 Subject: [PATCH 27/50] Improve types of partitioner input --- src/utils.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 1bad3782..7ad1e10c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -696,12 +696,13 @@ function partitioner_input(model, parameters; conn = :trans) @assert conn == :unit T = ones(Int, length(trans)) end - groups = [] + groups = Vector{Vector{Int}}() if model isa MultiModel for (k, m) in pairs(model.models) wg = physical_representation(m.domain) if wg isa WellDomain - push!(groups, copy(wg.perforations.reservoir)) + rcells = vec(Int.(wg.perforations.reservoir)) + push!(groups, rcells) end end end From 1d42d30a795839ed67ea22e375204e03c00dc3a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 10 Oct 2023 15:00:46 +0200 Subject: [PATCH 28/50] Handle lower case months in input file --- src/InputParser/deckinput/keywords/schedule.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/InputParser/deckinput/keywords/schedule.jl b/src/InputParser/deckinput/keywords/schedule.jl index 5beb0d8e..9aef77dc 100644 --- a/src/InputParser/deckinput/keywords/schedule.jl +++ b/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 From 4def7423b33c9a2716f854eff610264a388f1008 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 11 Oct 2023 15:10:26 +0200 Subject: [PATCH 29/50] Handle missing saturations for avg density --- src/multicomponent/flux.jl | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/multicomponent/flux.jl b/src/multicomponent/flux.jl index e8cb21aa..4ce3dcd1 100644 --- a/src/multicomponent/flux.jl +++ b/src/multicomponent/flux.jl @@ -48,9 +48,22 @@ function face_average_density(model::CompositionalModel, state, tpfa, phase) s = state.Saturations l = tpfa.left r = tpfa.right + ϵ = MINIMUM_COMPOSITIONAL_SATURATION @inbounds s_l = s[phase, l] @inbounds s_r = s[phase, r] @inbounds ρ_l = ρ[phase, l] @inbounds ρ_r = ρ[phase, r] - return (s_l*ρ_r + s_r*ρ_l)/max(s_l + s_r, 1e-8) + + s_l_tiny = s_l <= ϵ + s_r_tiny = s_r <= ϵ + if s_l_tiny && s_r_tiny + ρ_avg = zero(s_l) + elseif s_l_tiny + ρ_avg = ρ_r + elseif s_r_tiny + ρ_avg = ρ_l + else + ρ_avg = (s_l*ρ_r + s_r*ρ_l)/(s_l + s_r) + end + return ρ_avg end From 240999dbb0d0e2e194ecc0a9530fde110ea74895 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 13 Oct 2023 15:11:00 +0200 Subject: [PATCH 30/50] Handle step_index in state0 --- src/facility/types.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/facility/types.jl b/src/facility/types.jl index 311689ea..39022193 100644 --- a/src/facility/types.jl +++ b/src/facility/types.jl @@ -298,7 +298,7 @@ function Jutul.numerical_type(tc::WellGroupConfiguration) end function Jutul.update_values!(old::WellGroupConfiguration, new::WellGroupConfiguration) - return WellGroupConfiguration(copy(new.operating_controls), copy(new.requested_controls), copy(new.limits)) + return WellGroupConfiguration(copy(new.operating_controls), copy(new.requested_controls), copy(new.limits), new.step_index) end operating_control(cfg::WellGroupConfiguration, well::Symbol) = cfg.operating_controls[well] From 1984814daae46711421574b9e3700488f7d09708 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 13 Oct 2023 15:11:07 +0200 Subject: [PATCH 31/50] Remove unused function --- src/flux.jl | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/flux.jl b/src/flux.jl index 6e3beb94..bd4ed32e 100644 --- a/src/flux.jl +++ b/src/flux.jl @@ -115,16 +115,6 @@ end end end -@inline function saturation_averaged_density(ρ, ph, sat, c1, c2) - @inbounds ρ_1 = ρ[ph, c1] - @inbounds ρ_2 = ρ[ph, c2] - @inbounds S_1 = sat[ph, c1] - @inbounds S_2 = sat[ph, c2] - - avg = (ρ_1*S_1 + ρ_2*S_2)/max(S_1 + S_2, 1e-12) - return avg -end - export component_mass_fluxes!, update_total_masses! """ component_mass_fluxes!(q, face, state, model, kgrad, upw) From b3cd74ac0a6cb086876fd2feaae81299672d9a22 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 14 Oct 2023 20:44:10 +0200 Subject: [PATCH 32/50] Robustness fixes to init --- src/init/init.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/init/init.jl b/src/init/init.jl index 72a32028..c427c880 100644 --- a/src/init/init.jl +++ b/src/init/init.jl @@ -3,6 +3,7 @@ function equilibriate_state(model, contacts, datum_depth = missing, datum_pressu cells = missing, rs = missing, kwarg...) + model = reservoir_model(model) D = model.data_domain G = physical_representation(D) cc = D[:cell_centroids][3, :] @@ -200,7 +201,7 @@ function init_reference_pressure(pressures, contacts, kr, pc, ref_ix = 2) p[i] = pressures[ref_ix, i] if kr[ref_ix, i] <= 0 for ph in 1:nph - if kr[ph, i] > 0 + if kr[ph, i] > 1e-12 p[i] = pressures[ph, i]# - pc[ph, i] end end @@ -274,11 +275,11 @@ end function determine_saturations(depths, contacts, pressures; s_min = missing, s_max = missing, pc = nothing) nc = length(depths) nph = length(contacts) + 1 - if ismissing(s_max) - s_max = zeros(nph) + if ismissing(s_min) + s_min = zeros(nph) end if ismissing(s_max) - s_max = zeros(nph) + s_max = ones(nph) end sat = zeros(nph, nc) sat_pc = similar(sat) From 14914d21f8e68f703f16ad21a266fe2b70531c4b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 14 Oct 2023 20:50:59 +0200 Subject: [PATCH 33/50] Fix calling order for reference pressure with immobile phase --- src/init/init.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/init/init.jl b/src/init/init.jl index c427c880..24202518 100644 --- a/src/init/init.jl +++ b/src/init/init.jl @@ -84,7 +84,7 @@ function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pr JutulDarcy.update_kr!(kr, relperm, model, s, cells) init[:Saturations] = s - init[:Pressure] = init_reference_pressure(pressures, contacts, pc, kr, 2) + init[:Pressure] = init_reference_pressure(pressures, contacts, kr, pc, 2) return init end @@ -197,12 +197,15 @@ end function init_reference_pressure(pressures, contacts, kr, pc, ref_ix = 2) nph, nc = size(kr) p = zeros(nc) + ϵ = 1e-12 for i in eachindex(p) p[i] = pressures[ref_ix, i] - if kr[ref_ix, i] <= 0 + kr_ref = kr[ref_ix, i] + @assert kr_ref >= -ϵ + if kr[ref_ix, i] <= ϵ for ph in 1:nph - if kr[ph, i] > 1e-12 - p[i] = pressures[ph, i]# - pc[ph, i] + if kr[ph, i] > ϵ + p[i] = pressures[ph, i] end end end From be43c9f6b5c25545b83022baf71bdc0d5997c353 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 14 Oct 2023 21:22:44 +0200 Subject: [PATCH 34/50] Put in place basics of diffusion --- src/multiphase.jl | 39 +++++++++++++++++++++++++++++++++++++++ src/utils.jl | 13 ++++++++++--- 2 files changed, 49 insertions(+), 3 deletions(-) diff --git a/src/multiphase.jl b/src/multiphase.jl index 0abfb277..958b1f28 100644 --- a/src/multiphase.jl +++ b/src/multiphase.jl @@ -158,6 +158,45 @@ function Jutul.default_parameter_values(data_domain, model, param::Transmissibil return T end +struct Diffusivities <: VectorVariables end +Jutul.variable_scale(::Diffusivities) = 1e-10 +Jutul.minimum_value(::Diffusivities) = 0.0 + +Jutul.associated_entity(::Diffusivities) = Faces() +Jutul.values_per_entity(model, ::Diffusivities) = number_of_phases(model.system) + +function Jutul.default_parameter_values(data_domain, model, param::Diffusivities, symb) + nf = number_of_faces(model.domain) + nph = number_of_phases(model.system) + if haskey(data_domain, :diffusivities, Faces()) + # This takes precedence + T = data_domain[:diffusivities] + elseif haskey(data_domain, :diffusion, Cells()) + T = zeros(nph, nf) + U = data_domain[:diffusion] + g = physical_representation(data_domain) + if U isa AbstractVector + T_i = compute_face_trans(g, U) + for i in 1:nph + T[i, :] .= T_i + end + else + for i in 1:nph + T_i = compute_face_trans(g, U[i, :]) + T[i, :] .= T_i + end + end + if any(x -> x < 0, T) + c = count(x -> x < 0, T) + @warn "$c negative diffusivities detected." + end + else + error(":diffusion or :diffusivities symbol must be present in DataDomain to initialize parameter $symb, had keys: $(keys(data_domain))") + end + @assert size(T) == (nph, nf) + return T +end + struct TwoPointGravityDifference <: ScalarVariable end Jutul.associated_entity(::TwoPointGravityDifference) = Faces() diff --git a/src/utils.jl b/src/utils.jl index 7ad1e10c..ff4bfae9 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -8,7 +8,7 @@ reservoir_storage(model::MultiModel, storage) = storage.Reservoir export reservoir_domain """ - reservoir_domain(g; permeability = 9.869232667160130e-14, porosity = 0.1, kwarg...) + reservoir_domain(g; permeability = convert_to_si(0.1, :darcy), porosity = 0.1, kwarg...) Set up a `DataDomain` instance for given mesh or other representation `g`. `permeability` and `porosity` are then added to the domain. If scalars are @@ -17,7 +17,10 @@ cells. Permeability is either one value per cell (diagonal scalar), one value per dimension given in each row (for a diagonal tensor) or a vector that represents a compact full tensor representation (6 elements in 3D, 3 in 2D). """ -function reservoir_domain(g; permeability = 9.869232667160130e-14, porosity = 0.1, kwarg...) +function reservoir_domain(g; permeability = convert_to_si(0.1, :darcy), porosity = 0.1, diffusion = missing, kwarg...) + if !ismissing(diffusion) + kwarg = (diffusion = diffusion, kwarg...) + end return DataDomain(g; permeability = permeability, porosity = porosity, kwarg...) end @@ -53,12 +56,16 @@ function setup_reservoir_model(reservoir::DataDomain, system; kwarg... ) # We first set up the reservoir - models[:Reservoir] = SimulationModel( + rmodel = SimulationModel( reservoir, system, context = reservoir_context, general_ad = general_ad ) + if haskey(reservoir, :diffusion) || haskey(reservoir, :diffusivity) + rmodel.parameters[:Diffusivities] = Diffusivities() + end + models[:Reservoir] = rmodel # Then we set up all the wells mode = PredictionMode() if length(wells) > 0 From c74adfbfb079e241dd0cb6c50acf5b289e1719a6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 14 Oct 2023 21:58:01 +0200 Subject: [PATCH 35/50] Add diffusion impl for BO --- src/blackoil/flux.jl | 32 ++++++++++++++++++++++++++++++++ src/flux.jl | 7 ++++++- src/multiphase.jl | 4 +++- 3 files changed, 41 insertions(+), 2 deletions(-) diff --git a/src/blackoil/flux.jl b/src/blackoil/flux.jl index 44d5ea1c..71b00a72 100644 --- a/src/blackoil/flux.jl +++ b/src/blackoil/flux.jl @@ -48,6 +48,38 @@ else q_v = rhoVS*λb_v*∇ψ_v end + + if haskey(state, :Diffusivities) + S = state.Saturations + density = state.PhaseMassDensities + D = state.Diffusivities + @inbounds D_l = D[l, face] + @inbounds D_v = D[v, face] + + if has_disgas(sys) + X_o = cell -> @inbounds rhoLS/(Rs[cell]*rhoVS) + X_g = cell -> @inbounds (Rs[cell]*rhoLS)/(Rs[cell]*rhoVS) + + ΔX_o = -gradient(X_o, kgrad) + ΔX_g = -gradient(X_g, kgrad) + + mass_l = cell -> density[l, cell]*S[l, cell] + q_l += D_l*upwind(upw, mass_l, ΔX_o) + q_v += D_l*upwind(upw, mass_l, ΔX_g) + end + + if has_vapoil(sys) + Y_o = cell -> @inbounds rhoVS/(Rv[cell]*rhoLS) + Y_g = cell -> @inbounds (Rv[cell]*rhoVS)/(Rv[cell]*rhoLS) + + ΔY_o = -gradient(Y_o, kgrad) + ΔY_g = -gradient(Y_g, kgrad) + + mass_v = cell -> density[v, cell]*S[v, cell] + q_l += D_v*upwind(upw, mass_v, ΔY_o) + q_v += D_v*upwind(upw, mass_v, ΔX_g) + end + end q = setindex(q, q_l, l) q = setindex(q, q_v, v) return q diff --git a/src/flux.jl b/src/flux.jl index bd4ed32e..5c3f2147 100644 --- a/src/flux.jl +++ b/src/flux.jl @@ -63,7 +63,7 @@ function face_average_density(model, state, tpfa, phase) return 0.5*(ρ_i + ρ_c) end -@inline function gradient(X, tpfa::TPFA) +@inline function gradient(X::AbstractVector, tpfa::TPFA) return @inbounds X[tpfa.right] - X[tpfa.left] end @@ -71,6 +71,11 @@ end return @inbounds X[i, tpfa.right] - X[i, tpfa.left] end +@inline function gradient(F, tpfa::TPFA) + # Function handle version + return F(tpfa.right) - F(tpfa.left) +end + pressure_gradient(state, disc) = gradient(state.Pressure, disc) @inline function upwind(upw::SPU, F, q) diff --git a/src/multiphase.jl b/src/multiphase.jl index 958b1f28..44c9ca4b 100644 --- a/src/multiphase.jl +++ b/src/multiphase.jl @@ -173,7 +173,9 @@ function Jutul.default_parameter_values(data_domain, model, param::Diffusivities T = data_domain[:diffusivities] elseif haskey(data_domain, :diffusion, Cells()) T = zeros(nph, nf) - U = data_domain[:diffusion] + ϕ = data_domain[:porosity] + D = data_domain[:diffusion] + U = ϕ'.*D g = physical_representation(data_domain) if U isa AbstractVector T_i = compute_face_trans(g, U) From 731ce579be501a033803b5b07d235142c2936e79 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 14 Oct 2023 22:16:43 +0200 Subject: [PATCH 36/50] Correct typo in diffusion --- src/blackoil/flux.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blackoil/flux.jl b/src/blackoil/flux.jl index 71b00a72..1f4f1aba 100644 --- a/src/blackoil/flux.jl +++ b/src/blackoil/flux.jl @@ -77,7 +77,7 @@ mass_v = cell -> density[v, cell]*S[v, cell] q_l += D_v*upwind(upw, mass_v, ΔY_o) - q_v += D_v*upwind(upw, mass_v, ΔX_g) + q_v += D_v*upwind(upw, mass_v, ΔY_g) end end q = setindex(q, q_l, l) From bb20d400852ad7a1cacfe542fb59d50f157f488a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 14 Oct 2023 22:41:03 +0200 Subject: [PATCH 37/50] Work on diffusion --- src/blackoil/flux.jl | 33 ++++++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 9 deletions(-) diff --git a/src/blackoil/flux.jl b/src/blackoil/flux.jl index 1f4f1aba..baed805a 100644 --- a/src/blackoil/flux.jl +++ b/src/blackoil/flux.jl @@ -55,29 +55,33 @@ D = state.Diffusivities @inbounds D_l = D[l, face] @inbounds D_v = D[v, face] - + ϵ = 1e-10 if has_disgas(sys) - X_o = cell -> @inbounds rhoLS/(Rs[cell]*rhoVS) - X_g = cell -> @inbounds (Rs[cell]*rhoLS)/(Rs[cell]*rhoVS) + X_o = cell -> black_oil_phase_mass_fraction(rhoLS, rhoVS, Rs, cell) + X_g = cell -> 1.0 - black_oil_phase_mass_fraction(rhoLS, rhoVS, Rs, cell) ΔX_o = -gradient(X_o, kgrad) ΔX_g = -gradient(X_g, kgrad) mass_l = cell -> density[l, cell]*S[l, cell] - q_l += D_l*upwind(upw, mass_l, ΔX_o) - q_v += D_l*upwind(upw, mass_l, ΔX_g) + q_l += D_l*upwind(upw, mass_l, ΔX_o)*ΔX_o + q_v += D_l*upwind(upw, mass_l, ΔX_g)*ΔX_g + # q_l += D_v*face_average(mass_l, kgrad)*ΔX_o + # q_v += D_v*face_average(mass_l, kgrad)*ΔX_g end if has_vapoil(sys) - Y_o = cell -> @inbounds rhoVS/(Rv[cell]*rhoLS) - Y_g = cell -> @inbounds (Rv[cell]*rhoVS)/(Rv[cell]*rhoLS) + Y_o = cell -> 1.0 - black_oil_phase_mass_fraction(rhoVS, rhoLS, Rv, cell) + Y_g = cell -> black_oil_phase_mass_fraction(rhoVS, rhoLS, Rv, cell) ΔY_o = -gradient(Y_o, kgrad) ΔY_g = -gradient(Y_g, kgrad) mass_v = cell -> density[v, cell]*S[v, cell] - q_l += D_v*upwind(upw, mass_v, ΔY_o) - q_v += D_v*upwind(upw, mass_v, ΔY_g) + q_l += D_v*upwind(upw, mass_v, ΔY_o)*ΔY_o + q_v += D_v*upwind(upw, mass_v, ΔY_g)*ΔY_g + # q_l += D_v*face_average(mass_v, kgrad)*ΔY_o + # q_v += D_v*face_average(mass_v, kgrad)*ΔY_g end end q = setindex(q, q_l, l) @@ -85,6 +89,17 @@ return q 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 = zero(rs) + else + v = rhoLS/(rs*rhoVS) + end + return v +end + function apply_flow_bc!(acc, q, bc, model::StandardBlackOilModel, state, time) mu = state.PhaseViscosities b = state.ShrinkageFactors From d780f93c916aaaac2e6ae634a90aef5fe3fd1d02 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 14 Oct 2023 22:43:03 +0200 Subject: [PATCH 38/50] Add missing density to mass fraction --- src/blackoil/flux.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/blackoil/flux.jl b/src/blackoil/flux.jl index baed805a..ca8f9094 100644 --- a/src/blackoil/flux.jl +++ b/src/blackoil/flux.jl @@ -55,7 +55,6 @@ D = state.Diffusivities @inbounds D_l = D[l, face] @inbounds D_v = D[v, face] - ϵ = 1e-10 if has_disgas(sys) X_o = cell -> black_oil_phase_mass_fraction(rhoLS, rhoVS, Rs, cell) X_g = cell -> 1.0 - black_oil_phase_mass_fraction(rhoLS, rhoVS, Rs, cell) @@ -95,7 +94,7 @@ end if rs < 1e-10 v = zero(rs) else - v = rhoLS/(rs*rhoVS) + v = rhoLS/(rhoLS + rs*rhoVS) end return v end From 95d459a4b92c963ee3afc84b545ef09b7f7afc26 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 14 Oct 2023 22:57:34 +0200 Subject: [PATCH 39/50] Update flux.jl --- src/blackoil/flux.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/blackoil/flux.jl b/src/blackoil/flux.jl index ca8f9094..9d7a007e 100644 --- a/src/blackoil/flux.jl +++ b/src/blackoil/flux.jl @@ -63,10 +63,12 @@ ΔX_g = -gradient(X_g, kgrad) mass_l = cell -> density[l, cell]*S[l, cell] - q_l += D_l*upwind(upw, mass_l, ΔX_o)*ΔX_o - q_v += D_l*upwind(upw, mass_l, ΔX_g)*ΔX_g - # q_l += D_v*face_average(mass_l, kgrad)*ΔX_o - # q_v += D_v*face_average(mass_l, kgrad)*ΔX_g + # 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 + q_l += D_v*face_average(mass_l, kgrad)*ΔX_o + q_v += D_v*face_average(mass_l, kgrad)*ΔX_g end if has_vapoil(sys) @@ -92,7 +94,7 @@ end # TODO: Should have molar weights here maybe, but not part of standard input @inbounds rs = Rs[cell] if rs < 1e-10 - v = zero(rs) + v = one(rs) else v = rhoLS/(rhoLS + rs*rhoVS) end From f6993e4535258262fd9270c34d7866b9a5f2c6b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 14 Oct 2023 22:57:38 +0200 Subject: [PATCH 40/50] Update flux.jl --- src/flux.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/flux.jl b/src/flux.jl index 5c3f2147..3e46ec7c 100644 --- a/src/flux.jl +++ b/src/flux.jl @@ -76,6 +76,15 @@ end return F(tpfa.right) - F(tpfa.left) end +function face_average(F, tpfa) + return 0.5*(F(tpfa.right) + F(tpfa.left)) +end + +function phase_face_average(phase_property, tpfa, cell) + F(cell) = @inbounds phase_property[phase, cell] + return face_average(F, tpfa) +end + pressure_gradient(state, disc) = gradient(state.Pressure, disc) @inline function upwind(upw::SPU, F, q) From 1bd619d764ca17852a5ccd7bc4a6c6c1d323b89b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 15 Oct 2023 21:26:55 +0200 Subject: [PATCH 41/50] Simplify BO diffusive flux --- src/blackoil/flux.jl | 52 ++++++++++++++++++++++---------------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/src/blackoil/flux.jl b/src/blackoil/flux.jl index 9d7a007e..40f6a406 100644 --- a/src/blackoil/flux.jl +++ b/src/blackoil/flux.jl @@ -53,36 +53,16 @@ S = state.Saturations density = state.PhaseMassDensities D = state.Diffusivities - @inbounds D_l = D[l, face] - @inbounds D_v = D[v, face] if has_disgas(sys) - X_o = cell -> black_oil_phase_mass_fraction(rhoLS, rhoVS, Rs, cell) - X_g = cell -> 1.0 - black_oil_phase_mass_fraction(rhoLS, rhoVS, Rs, cell) - - ΔX_o = -gradient(X_o, kgrad) - ΔX_g = -gradient(X_g, kgrad) - - mass_l = cell -> density[l, cell]*S[l, 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 - q_l += D_v*face_average(mass_l, kgrad)*ΔX_o - q_v += D_v*face_average(mass_l, kgrad)*ΔX_g + 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) - Y_o = cell -> 1.0 - black_oil_phase_mass_fraction(rhoVS, rhoLS, Rv, cell) - Y_g = cell -> black_oil_phase_mass_fraction(rhoVS, rhoLS, Rv, cell) - - ΔY_o = -gradient(Y_o, kgrad) - ΔY_g = -gradient(Y_g, kgrad) - - mass_v = cell -> density[v, cell]*S[v, cell] - q_l += D_v*upwind(upw, mass_v, ΔY_o)*ΔY_o - q_v += D_v*upwind(upw, mass_v, ΔY_g)*ΔY_g - # q_l += D_v*face_average(mass_v, kgrad)*ΔY_o - # q_v += D_v*face_average(mass_v, kgrad)*ΔY_g + 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) @@ -90,6 +70,26 @@ 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) +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] From 2370dc6d2c1b9e0c73a31e814b259628a538a1e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 15 Oct 2023 21:47:45 +0200 Subject: [PATCH 42/50] Add array version of upwind --- src/blackoil/flux.jl | 8 +++----- src/flux.jl | 13 ++++++++++++- 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/src/blackoil/flux.jl b/src/blackoil/flux.jl index 40f6a406..803022f8 100644 --- a/src/blackoil/flux.jl +++ b/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,8 +40,7 @@ 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 @@ -87,7 +85,7 @@ function blackoil_diffusion(R, S, density, rhoS_self, rhoS_dissolved, face, D, 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) + return (diff_self::T, diff_dissolved::T)::Tuple{T, T} end @inline function black_oil_phase_mass_fraction(rhoLS, rhoVS, Rs, cell) diff --git a/src/flux.jl b/src/flux.jl index 3e46ec7c..d8e473ef 100644 --- a/src/flux.jl +++ b/src/flux.jl @@ -88,7 +88,7 @@ end pressure_gradient(state, disc) = gradient(state.Pressure, disc) @inline function upwind(upw::SPU, F, q) - flag = q < 0 + flag = q < zero(q) if flag up = upw.right else @@ -97,6 +97,17 @@ pressure_gradient(state, disc) = gradient(state.Pressure, disc) return F(up) end + +@inline function upwind(upw::SPU, X::AbstractArray, q) + flag = q < zero(q) + if flag + up = upw.right + else + up = upw.left + end + return @inbounds X[up] +end + @inline function phase_upwind(upw, m::AbstractMatrix, phase::Integer, q) F(cell) = @inbounds m[phase, cell] return upwind(upw, F, q) From a72b3c9be2137c1e0c596357199063c5e0bdc5b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 15 Oct 2023 21:47:55 +0200 Subject: [PATCH 43/50] Improve TopConditions constructor allocations --- src/types.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/types.jl b/src/types.jl index d018b38c..b58decef 100644 --- a/src/types.jl +++ b/src/types.jl @@ -444,19 +444,19 @@ struct TopConditions{N, R} density::SVector{N, R} volume_fractions::SVector{N, R} function TopConditions(n::Int, R::DataType = Float64; density = missing, volume_fractions = missing) - function internal_convert(x::Missing) + function internal_convert(x::Missing, N) x = @SVector ones(R, n) - return x./n + return x./N end - function internal_convert(x) + function internal_convert(x, N) x0 = x - @assert length(x0) == n - x = @MVector zeros(R, n) + @assert length(x0) == N + x = @MVector zeros(R, N) @. x = x0 return SVector(x) end - density = internal_convert(density) - volume_fractions = internal_convert(volume_fractions) + density = internal_convert(density, n) + volume_fractions = internal_convert(volume_fractions, n) return new{n, R}(density, volume_fractions) end end From 9de3f45aaae191757136ef160bc0e63be26d3cb1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Mon, 16 Oct 2023 13:27:03 +0200 Subject: [PATCH 44/50] Add unit label to well plot --- ext/JutulDarcyGLMakieExt/well_plots.jl | 49 +++++++++++++++++++------- 1 file changed, 36 insertions(+), 13 deletions(-) diff --git a/ext/JutulDarcyGLMakieExt/well_plots.jl b/ext/JutulDarcyGLMakieExt/well_plots.jl index 650ad22a..2641079a 100644 --- a/ext/JutulDarcyGLMakieExt/well_plots.jl +++ b/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,9 +129,6 @@ 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) @@ -114,6 +136,7 @@ function JutulDarcy.plot_well_results(well_data::Vector, time = nothing; start_d on(type_menu.selection) do s val = findfirst(isequal(s), respstr) + y_l[] = response_label_to_unit(s) response_ix[] = val autolimits!(ax) end From 5ce825496a9f1b36e624117339d8c74461572630 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 17 Oct 2023 12:58:50 +0200 Subject: [PATCH 45/50] Avoid floating point equality for well disable --- src/facility/wellgroups.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/facility/wellgroups.jl b/src/facility/wellgroups.jl index dcccce75..1b8d3736 100644 --- a/src/facility/wellgroups.jl +++ b/src/facility/wellgroups.jl @@ -25,7 +25,7 @@ function Jutul.update_primary_variable!(state, massrate::TotalSurfaceMassRate, s return Jutul.update_value(v, dx) end function do_update!(wcfg, s, v, dx, ctrl::InjectorControl) - if v == MIN_ACTIVE_WELL_RATE && v + dx < MIN_ACTIVE_WELL_RATE && can_shut_wells + if v <= MIN_ACTIVE_WELL_RATE && v + dx < MIN_ACTIVE_WELL_RATE && can_shut_wells 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)) @@ -35,7 +35,8 @@ function Jutul.update_primary_variable!(state, massrate::TotalSurfaceMassRate, s return next end function do_update!(wcfg, s, v, dx, ctrl::ProducerControl) - if v == MIN_ACTIVE_WELL_RATE && v + dx > MIN_ACTIVE_WELL_RATE && can_shut_wells + # A significant negative rate is the valid producer control + if v >= MIN_ACTIVE_WELL_RATE && v + dx > MIN_ACTIVE_WELL_RATE && can_shut_wells 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)) From e8797a6ce84aee2a115ade35b39dd808ddd28edf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 17 Oct 2023 14:22:12 +0200 Subject: [PATCH 46/50] Make it possible to turn off explicit dp updates --- src/facility/controls.jl | 4 ++-- src/facility/wells/stdwells.jl | 4 ++-- src/facility/wells/wells.jl | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/facility/controls.jl b/src/facility/controls.jl index 6c984716..c762b7c0 100644 --- a/src/facility/controls.jl +++ b/src/facility/controls.jl @@ -3,7 +3,7 @@ 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; recorder, kwarg...) +function Jutul.update_before_step_multimodel!(storage_g, model_g::MultiModel, model::WellGroupModel, dt, forces, key; time = NaN, recorder, update_explicit = true) # Set control to whatever is on the forces storage = storage_g[key] cfg = storage.state.WellGroupConfiguration @@ -44,7 +44,7 @@ function Jutul.update_before_step_multimodel!(storage_g, model_g::MultiModel, mo 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 diff --git a/src/facility/wells/stdwells.jl b/src/facility/wells/stdwells.jl index 1deba1c7..98b5cd29 100644 --- a/src/facility/wells/stdwells.jl +++ b/src/facility/wells/stdwells.jl @@ -79,8 +79,8 @@ 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) - if well_has_explicit_pressure_drop(well_model) +function update_before_step_well!(well_state, well_model::SimpleWellFlowModel, res_state, res_model, ctrl; 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) end diff --git a/src/facility/wells/wells.jl b/src/facility/wells/wells.jl index 2bb7851c..8f416fcc 100644 --- a/src/facility/wells/wells.jl +++ b/src/facility/wells/wells.jl @@ -197,7 +197,7 @@ end include("mswells_equations.jl") -function update_before_step_well!(well_state, well_model, res_state, res_model, ctrl) +function update_before_step_well!(well_state, well_model, res_state, res_model, ctrl; kwarg...) end From af4b72891350eb14bdd629f5b2b49b7557764b37 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 18 Oct 2023 15:38:44 +0200 Subject: [PATCH 47/50] Make progress recorder optional for well update --- src/facility/controls.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/facility/controls.jl b/src/facility/controls.jl index c762b7c0..2e81b80f 100644 --- a/src/facility/controls.jl +++ b/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; time = NaN, recorder, update_explicit = true) +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 From 7989b0f735da8d46215fa76d9682c371eb4a066e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 18 Oct 2023 15:38:52 +0200 Subject: [PATCH 48/50] Comment on test --- test/multiphase.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/multiphase.jl b/test/multiphase.jl index 86f8c272..60172f7e 100644 --- a/test/multiphase.jl +++ b/test/multiphase.jl @@ -1,6 +1,7 @@ using JutulDarcy, Jutul using Test +# Note: This test can get failures if HYPRE is loaded on some of the simple cases. function test_multiphase(grid = CartesianMesh((2, 2), (2.0, 2.0)); setup = "two_phase_simple", debug_level = 1, linear_solver = nothing, kwarg...) state0, model, prm, f, t = get_test_setup(grid, case_name = setup; kwarg...) sim = Simulator(model, state0 = state0, parameters = prm) From a42ff9c8804a47fbeb183381378b0cafff849a72 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 18 Oct 2023 15:39:03 +0200 Subject: [PATCH 49/50] Update JutulDarcyHYPREExt.jl Disable precompilation for extension --- ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl b/ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl index 1725c2d9..707efe1e 100644 --- a/ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl +++ b/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 From a52060b8d47167f12759f1de88faeb16d6a3060d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 18 Oct 2023 16:01:27 +0200 Subject: [PATCH 50/50] Bump version + compat --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index d1562821..937930ba 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "JutulDarcy" uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a" authors = ["Olav Møyner "] -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"