Skip to content

Commit

Permalink
Merge pull request #46 from sintefmath/dev
Browse files Browse the repository at this point in the history 
Improved initialization, thermal output and grid processing
  • Loading branch information
@moyner
moyner authored Jun 5, 2024
2 parents 0280f91 + e5b6897 commit 32a11ad
Show file tree
Hide file tree
Showing 42 changed files with 1,311 additions and 208 deletions.
6 changes: 3 additions & 3 deletions Project.toml
View file
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.23"
version = "0.2.24"

[deps]
AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
Expand Down Expand Up @@ -49,9 +49,9 @@ DelimitedFiles = "1.9.1"
DocStringExtensions = "0.9.3"
ForwardDiff = "0.10.30"
GLMakie = "0.9.2"
GeoEnergyIO = "1.1"
GeoEnergyIO = "1.1.2"
HYPRE = "1.4.0"
Jutul = "0.2.31"
Jutul = "0.2.33"
Krylov = "0.9.1"
LinearAlgebra = "1"
LoopVectorization = "0.12.115"
Expand Down
7 changes: 6 additions & 1 deletion docs/make.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,12 @@ function build_jutul_darcy_docs(build_format = nothing; build_examples = true, b
"Intro to compositional flow" => "co2_brine_2d_vertical",
"Compositional with five components" => "compositional_5components",
"Parameter optimization of Buckley-Leverett" => "optimize_simple_bl",
"Validation of reservoir simulator" => "mrst_validation"
"Validation: SPE1" => "validation_spe1",
"Validation: SPE9" => "validation_spe9",
"Validation: Egg" => "validation_egg",
"Validation: OLYMPUS 1" => "validation_olympus_1",
"Validation: Norne" => "validation_norne_nohyst",
"Validation: MRST input files" => "validation_mrst"
]
examples_markdown = []
function update_footer(content, pth)
Expand Down
2 changes: 1 addition & 1 deletion examples/data_input_file.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#
# We begin by loading the SPE9 dataset via the GeoEnergyIO package.
using JutulDarcy
pth = JutulDarcy.GeoEnergyIO.test_input_file_path("spe9", "SPE9_CP.DATA", base = "opm-tests")
pth = JutulDarcy.GeoEnergyIO.test_input_file_path("SPE9", "SPE9.DATA")
# ## Set up and run a simulation
# If we do not need the case, we could also have done:
# ws, states = simulate_data_file(pth)
Expand Down
2 changes: 1 addition & 1 deletion examples/intro_example.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ density = ConstantCompressibilityDensities(
compressibility = c
)
kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 3.0])
replace_variables!(model, PhaseMassDensities = density, PhaseRelativePermeability = kr);
replace_variables!(model, PhaseMassDensities = density, RelativePermeabilities = kr);

state0 = setup_reservoir_state(model,
Pressure = 120bar,
Expand Down
99 changes: 99 additions & 0 deletions examples/validation_egg.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,99 @@
# # Validation: The Egg model (oil-water compressible)
# A two-phase model that is taken from the first member of the EGG ensemble. The
# model is a synthetic case with channelized permeability and water injection
# with fixed controls. For more details, see the paper where the ensemble is
# introduced:
#
# [Jansen, Jan-Dirk, et al. "The egg model–a geological ensemble for reservoir
# simulation." Geoscience Data Journal 1.2 (2014): 192-195.](https://doi.org/10.1002/gdj3.21)

using Jutul, JutulDarcy, GLMakie, DelimitedFiles, HYPRE
egg_dir = JutulDarcy.GeoEnergyIO.test_input_file_path("EGG")
case = setup_case_from_data_file(joinpath(egg_dir, "EGG.DATA"))
ws, states = simulate_reservoir(case, output_substates = true)
# ## Plot the reservoir solution
plot_reservoir(case.model, states, step = 135, key = :Saturations)

# ## Load reference solution (OPM Flow)
# We load a CSV file with the reference solution and set up plotting
csv_path = joinpath(egg_dir, "REFERENCE.CSV")
data, header = readdlm(csv_path, ',', header = true)
time_ref = data[:, 1]
time_jutul = deepcopy(ws.time)
wells = deepcopy(ws.wells)
wnames = collect(keys(wells))
nw = length(wnames)
day = si_unit(:day)
cmap = :tableau_hue_circle

inj = Symbol[]
prod = Symbol[]
for (wellname, well) in pairs(wells)
qts = well[:wrat] + well[:orat]
if sum(qts) > 0
push!(inj, wellname)
else
push!(prod, wellname)
end
end

function plot_well_comparison(response, well_names, reponse_name = "$response")
fig = Figure(size = (1000, 400))
if response == :bhp
ys = 1/si_unit(:bar)
yl = "Bottom hole pressure / Bar"
elseif response == :wrat
ys = si_unit(:day)
yl = "Surface water rate / m³/day"
elseif response == :orat
ys = si_unit(:day)/(1000*si_unit(:stb))
yl = "Surface oil rate / 10³ stb/day"
else
error("$response not ready.")
end
welltypes = []
ax = Axis(fig[1:4, 1], xlabel = "Time / days", ylabel = yl)
i = 1
linehandles = []
linelabels = []
for well_name in well_names
well = wells[well_name]
label_in_csv = "$well_name:$response"
ref_pos = findfirst(x -> x == label_in_csv, vec(header))
qoi = copy(well[response]).*ys
qoi_ref = data[:, ref_pos].*ys
tot_rate = copy(well[:rate])
@. qoi[tot_rate == 0] = NaN
orat_ref = data[:, findfirst(x -> x == "$well_name:orat", vec(header))]
wrat_ref = data[:, findfirst(x -> x == "$well_name:wrat", vec(header))]
tot_rate_ref = orat_ref + wrat_ref
@. qoi_ref[tot_rate_ref == 0] = NaN
crange = (1, max(length(well_names), 2))
lh = lines!(ax, time_jutul./day, abs.(qoi),
color = i,
colorrange = crange,
label = "$well_name", colormap = cmap
)
push!(linehandles, lh)
push!(linelabels, "$well_name")
lines!(ax, time_ref./day, abs.(qoi_ref),
color = i,
colorrange = crange,
linestyle = :dash,
colormap = cmap
)
i += 1
end
l1 = LineElement(color = :black, linestyle = nothing)
l2 = LineElement(color = :black, linestyle = :dash)

Legend(fig[1:3, 2], linehandles, linelabels, nbanks = 3)
Legend(fig[4, 2], [l1, l2], ["JutulDarcy.jl", "E100"])
fig
end
# ## Well responses and comparison
# ### Water production rates
plot_well_comparison(:wrat, prod, "Water surface rate")
# ### Oil production rates
plot_well_comparison(:orat, prod, "Oil surface rate")

0 examples/mrst_validation.jl → examples/validation_mrst.jl
View file
File renamed without changes.
114 changes: 114 additions & 0 deletions examples/validation_norne_nohyst.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,114 @@
# # Norne: Real field black-oil model
# The Norne model is a real field model. The model has been adapted so that the
# input file only contains features present in JutulDarcy, with the most notable
# omissions being removal of hysteresis and threshold pressures between
# equilibriation reqgions. For more details, see the [OPM data
# webpage](https://opm-project.org/?page_id=559)
using Jutul, JutulDarcy, GLMakie, DelimitedFiles, HYPRE
norne_dir = JutulDarcy.GeoEnergyIO.test_input_file_path("NORNE_NOHYST")
case = setup_case_from_data_file(joinpath(norne_dir, "NORNE_NOHYST.DATA"))
ws, states = simulate_reservoir(case, output_substates = true)
# ## Plot the reservoir solution
fig = plot_reservoir(case.model, states, step = 247, key = :Saturations)
ax = fig.current_axis[]
ax.azimuth[] = -3.0
ax.elevation[] = 0.5
fig
# ## Plot the reservoir static properties
fig = plot_reservoir(case.model, key = :porosity)
ax = fig.current_axis[]
ax.azimuth[] = -3.0
ax.elevation[] = 0.5
fig
# ## Load reference and set up plotting
csv_path = joinpath(norne_dir, "REFERENCE.CSV")
data, header = readdlm(csv_path, ',', header = true)
time_ref = data[:, 1]
time_jutul = deepcopy(ws.time)
wells = deepcopy(ws.wells)
wnames = collect(keys(wells))
nw = length(wnames)
day = si_unit(:day)
cmap = :tableau_hue_circle

inj = Symbol[]
prod = Symbol[]
for (wellname, well) in pairs(wells)
qts = well[:wrat] + well[:orat] + well[:grat]
if sum(qts) > 0
push!(inj, wellname)
else
push!(prod, wellname)
end
end

function plot_well_comparison(response, well_names, reponse_name = "$response")
fig = Figure(size = (1000, 400))
if response == :bhp
ys = 1/si_unit(:bar)
yl = "Bottom hole pressure / Bar"
elseif response == :wrat
ys = si_unit(:day)
yl = "Surface water rate / m³/day"
elseif response == :grat
ys = si_unit(:day)/1e6
yl = "Surface gas rate / 10⁶ m³/day"
elseif response == :orat
ys = si_unit(:day)/(1000*si_unit(:stb))
yl = "Surface oil rate / 10³ stb/day"
else
error("$response not ready.")
end
welltypes = []
ax = Axis(fig[1:4, 1], xlabel = "Time / days", ylabel = yl)
i = 1
linehandles = []
linelabels = []
for well_name in well_names
well = wells[well_name]
label_in_csv = "$well_name:$response"
ref_pos = findfirst(x -> x == label_in_csv, vec(header))
qoi = copy(well[response]).*ys
qoi_ref = data[:, ref_pos].*ys
tot_rate = copy(well[:rate])
@. qoi[tot_rate == 0] = NaN
grat_ref = data[:, findfirst(x -> x == "$well_name:grat", vec(header))]
orat_ref = data[:, findfirst(x -> x == "$well_name:orat", vec(header))]
wrat_ref = data[:, findfirst(x -> x == "$well_name:wrat", vec(header))]
tot_rate_ref = grat_ref + orat_ref + wrat_ref
@. qoi_ref[tot_rate_ref == 0] = NaN
crange = (1, max(length(well_names), 2))
lh = lines!(ax, time_jutul./day, abs.(qoi),
color = i,
colorrange = crange,
label = "$well_name", colormap = cmap
)
push!(linehandles, lh)
push!(linelabels, "$well_name")
lines!(ax, time_ref./day, abs.(qoi_ref),
color = i,
colorrange = crange,
linestyle = :dash,
colormap = cmap
)
i += 1
end
l1 = LineElement(color = :black, linestyle = nothing)
l2 = LineElement(color = :black, linestyle = :dash)

Legend(fig[1:3, 2], linehandles, linelabels, nbanks = 3)
Legend(fig[4, 2], [l1, l2], ["JutulDarcy.jl", "OPM Flow"])
fig
end
# ## Injector bhp
plot_well_comparison(:bhp, inj, "Bottom hole pressure")
# ## Gas injection rates
plot_well_comparison(:grat, inj, "Gas surface injection rate")
# ## Water injection rates
plot_well_comparison(:wrat, inj, "Water surface injection rate")
# ## Producer bhp
plot_well_comparison(:bhp, prod, "Bottom hole pressure")
# ## Oil production rates
plot_well_comparison(:orat, prod, "Oil surface production rate")
# ## Gas production rates
plot_well_comparison(:grat, prod, "Gas surface production rate")
95 changes: 95 additions & 0 deletions examples/validation_olympus_1.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,95 @@
# # OLYMPUS_1 model
# Model from the [ISAPP Optimization challenge](https://www.isapp2.com/optimization-challenge/reservoir-model-description.html)
using Jutul, JutulDarcy, GLMakie, DelimitedFiles, HYPRE
olympus_dir = JutulDarcy.GeoEnergyIO.test_input_file_path("OLYMPUS_1")
case = setup_case_from_data_file(joinpath(olympus_dir, "OLYMPUS_1.DATA"), backend = :csr)
ws, states = simulate_reservoir(case, output_substates = true)
# ## Plot the reservoir
plot_reservoir(case.model, key = :porosity)
# ## Plot the saturations
plot_reservoir(case.model, states, step = 10, key = :Saturations)
# ## Load reference and set up plotting
csv_path = joinpath(olympus_dir, "REFERENCE.CSV")
data, header = readdlm(csv_path, ',', header = true)
time_ref = data[:, 1]
time_jutul = deepcopy(ws.time)
wells = deepcopy(ws.wells)
wnames = collect(keys(wells))
nw = length(wnames)
day = si_unit(:day)
cmap = :tableau_hue_circle

inj = Symbol[]
prod = Symbol[]
for (wellname, well) in pairs(wells)
qts = well[:wrat] + well[:orat]
if sum(qts) > 0
push!(inj, wellname)
else
push!(prod, wellname)
end
end

function plot_well_comparison(response, well_names, reponse_name = "$response"; ylims = missing)
fig = Figure(size = (1000, 400))
if response == :bhp
ys = 1/si_unit(:bar)
yl = "Bottom hole pressure / Bar"
elseif response == :wrat
ys = si_unit(:day)
yl = "Surface water rate / m³/day"
elseif response == :orat
ys = si_unit(:day)/(1000*si_unit(:stb))
yl = "Surface oil rate / 10³ stb/day"
else
error("$response not ready.")
end
welltypes = []
ax = Axis(fig[1:4, 1], xlabel = "Time / days", ylabel = yl)
i = 1
linehandles = []
linelabels = []
for well_name in well_names
well = wells[well_name]
label_in_csv = "$well_name:$response"
ref_pos = findfirst(x -> x == label_in_csv, vec(header))
qoi = copy(well[response]).*ys
qoi_ref = data[:, ref_pos].*ys
tot_rate = copy(well[:rate])
@. qoi[tot_rate == 0] = NaN
orat_ref = data[:, findfirst(x -> x == "$well_name:orat", vec(header))]
wrat_ref = data[:, findfirst(x -> x == "$well_name:wrat", vec(header))]
tot_rate_ref = orat_ref + wrat_ref
@. qoi_ref[tot_rate_ref == 0] = NaN
crange = (1, max(length(well_names), 2))
lh = lines!(ax, time_jutul./day, abs.(qoi),
color = i,
colorrange = crange,
label = "$well_name", colormap = cmap
)
push!(linehandles, lh)
push!(linelabels, "$well_name")
lines!(ax, time_ref./day, abs.(qoi_ref),
color = i,
colorrange = crange,
linestyle = :dash,
colormap = cmap
)
i += 1
if !ismissing(ylims)
ylims!(ax, ylims)
end
end
l1 = LineElement(color = :black, linestyle = nothing)
l2 = LineElement(color = :black, linestyle = :dash)

Legend(fig[1:3, 2], linehandles, linelabels, nbanks = 3)
Legend(fig[4, 2], [l1, l2], ["JutulDarcy.jl", "OPM Flow"])
fig
end
# ## Plot water production rates
plot_well_comparison(:wrat, prod, "Water surface production rate")
# ## Plot oil production rates
plot_well_comparison(:orat, prod, "Oil surface production rate", ylims = (0, 5))
# ## Plot water injection rates
plot_well_comparison(:wrat, inj, "Water injection rate")
Loading

4 comments on commit 32a11ad

@moyner
Copy link
Member Author

@moyner moyner commented on 32a11ad Jun 5, 2024

Choose a reason for hiding this comment

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

@JuliaRegistrator register

@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/108311

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

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.24 -m "<description of version>" 32a11ad7f315a1dff26fb73551366456c8760d2e
git push origin v0.2.24

@moyner
Copy link
Member Author

@moyner moyner commented on 32a11ad Jun 5, 2024

Choose a reason for hiding this comment

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

@JuliaRegistrator register

Release notes:

Updates

  • Improvements to equilibriation (maximum water saturation, non-unique data points for pressure, ...)
  • Better plotting of reservoirs and wells
  • More accurate well indices (effects from NTG + drainage radius)
  • Precompilation workflow now includes grid processing
  • Support for AMGCL in CPR
  • Better nonlinear performance for blackoil models
  • Added examples for Eclipse-type input: Norne, SPE1, SPE9, Egg and Olympus models with validation
  • Improved handling of many types of multipliers
  • New unexported routine reservoir_sensitivities computes sensitivities with respect to static properties (perm/poro/geometry/multipliers)
  • Output of thermal well results
  • Better NLDD defaults

@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 updated: JuliaRegistries/General/108311

Tagging

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.24 -m "<description of version>" 32a11ad7f315a1dff26fb73551366456c8760d2e
git push origin v0.2.24

Please sign in to comment.