Skip to content

Commit

Permalink
Merge pull request #21 from sintefmath/dev
Browse files Browse the repository at this point in the history
Better compositional solvers and fixes to simple wells + blackoil diffusion
  • Loading branch information
moyner authored Oct 18, 2023
2 parents 075db60 + a52060b commit a7da60f
Show file tree
Hide file tree
Showing 28 changed files with 499 additions and 187 deletions.
4 changes: 2 additions & 2 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "JutulDarcy"
uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a"
authors = ["Olav Møyner <[email protected]>"]
version = "0.2.9"
version = "0.2.10"

[deps]
AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
Expand Down Expand Up @@ -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"
Expand Down
49 changes: 36 additions & 13 deletions ext/JutulDarcyGLMakieExt/well_plots.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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"
Expand All @@ -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
Expand All @@ -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
Expand Down
18 changes: 9 additions & 9 deletions ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Original file line number Diff line number Diff line change
Expand Up @@ -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
3 changes: 2 additions & 1 deletion src/InputParser/deckinput/keywords/schedule.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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
Expand Down
54 changes: 50 additions & 4 deletions src/blackoil/flux.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
7 changes: 4 additions & 3 deletions src/blackoil/variables/varswitch.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
15 changes: 15 additions & 0 deletions src/ext.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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...)
Expand Down
14 changes: 11 additions & 3 deletions src/facility/controls.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
Expand All @@ -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"
Expand All @@ -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

Expand Down
19 changes: 12 additions & 7 deletions src/facility/cross_terms.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down
Loading

2 comments on commit a7da60f

@moyner
Copy link
Member Author

@moyner moyner commented on a7da60f Oct 18, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/93666

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.2.10 -m "<description of version>" a7da60f069ee02ba60ae9a9f756ce7abd63c1da4
git push origin v0.2.10

Please sign in to comment.