JutulDarcy.jl

Appearance

Adjoint gradients for the SPE1 model ​

This is a brief example of how to set up a case and compute adjoint gradients for the SPE1 model in a few different ways. These techniques are in general applicable to DATA-type models.

Load the model ​

julia
using Jutul, JutulDarcy, GeoEnergyIO, GLMakie
data_pth = joinpath(GeoEnergyIO.test_input_file_path("SPE1"), "SPE1.DATA")
data = parse_data_file(data_pth);
case = setup_case_from_data_file(data);

Set up a function to set up the case with custom porosity ​

We create a setup function that takes in a parameter dictionary prm and returns a case with the porosity set to the value in prm["poro"]. This is a convenient way to customize a case prior to setup, allowing direct interaction with keywords like PERMX, PORO, etc.

julia
function F(prm, step_info = missing)
    data_c = deepcopy(data)
    data_c["GRID"]["PORO"] = fill(prm["poro"], size(data_c["GRID"]["PORO"]))
    case = setup_case_from_data_file(data_c)
    return case
end
F (generic function with 2 methods)

Define and simulate truth case (poro from input) ​

julia
x_truth = only(unique(data["GRID"]["PORO"]))
prm_truth = Dict("poro" => x_truth)
case_truth = F(prm_truth)
ws, states = simulate_reservoir(case_truth)
ReservoirSimResult with 120 entries:

  wells (2 present):
    :INJ
    :PROD
    Results per well:
       :wrat => Vector{Float64} of size (120,)
       :Aqueous_mass_rate => Vector{Float64} of size (120,)
       :orat => Vector{Float64} of size (120,)
       :bhp => Vector{Float64} of size (120,)
       :gor => Vector{Float64} of size (120,)
       :lrat => Vector{Float64} of size (120,)
       :mass_rate => Vector{Float64} of size (120,)
       :rate => Vector{Float64} of size (120,)
       :Vapor_mass_rate => Vector{Float64} of size (120,)
       :control => Vector{Symbol} of size (120,)
       :Liquid_mass_rate => Vector{Float64} of size (120,)
       :wcut => Vector{Float64} of size (120,)
       :grat => Vector{Float64} of size (120,)

  states (Vector with 120 entries, reservoir variables for each state)
    :Pressure => Vector{Float64} of size (300,)
    :ImmiscibleSaturation => Vector{Float64} of size (300,)
    :BlackOilUnknown => Vector{BlackOilX{Float64}} of size (300,)
    :TotalMasses => Matrix{Float64} of size (3, 300)
    :Rs => Vector{Float64} of size (300,)
    :Saturations => Matrix{Float64} of size (3, 300)

  time (report time for each state)
     Vector{Float64} of length 120

  result (extended states, reports)
     SimResult with 120 entries

  extra
     Dict{Any, Any} with keys :simulator, :config

  Completed at Jun. 18 2025 14:06 after 1 second, 757 milliseconds, 430.4 microseconds.

Define pressure difference as objective function ​

julia
function pdiff(p, p0)
    v = 0.0
    for i in eachindex(p)
        v += (p[i] - p0[i])^2
    end
    return v
end

step_times = cumsum(case.dt)
function mismatch_objective(m, s, dt, step_info, forces)
    t = step_info[:time] + dt
    step = findmin(x -> abs(x - t), step_times)[2]
    p = s[:Reservoir][:Pressure]
    v = pdiff(p, states[step][:Pressure])
    return dt*(v/(si_unit(:bar)*100)^2)
end
mismatch_objective (generic function with 1 method)

Create a perturbed initial guess and optimize ​

We create a perturbed initial guess for the porosity and optimize the case using the mismatch objective function defined above. The optimization will adjust the porosity to minimize the mismatch between the simulated pressure and the truth case, and recover the original porosity value of 0.3.

julia
prm = Dict("poro" => x_truth .+ 0.25)
dprm = setup_reservoir_dict_optimization(prm, F)
free_optimization_parameter!(dprm, "poro", abs_max = 1.0, abs_min = 0.1)
prm_opt = optimize_reservoir(dprm, mismatch_objective);
dprm
DictParameters with 1 active parameters and 1 inactive:
Active optimization parameters
┌──────┬───────────────┬───────┬─────┬─────┬─────────────────┬────────┐
│ Name │ Initial value │ Count │ Min │ Max │ Optimized value │ Change │
├──────┼───────────────┼───────┼─────┼─────┼─────────────────┼────────┤
│ poro │ 0.55          │     1 │ 0.1 │ 1.0 │ 0.3             │ -45.0% │
└──────┴───────────────┴───────┴─────┴─────┴─────────────────┴────────┘
No inactive optimization parameters.

Plot the optimization history ​

julia
using GLMakie
fig = Figure()
ax = Axis(fig[1, 1], xlabel = "LBFGS iteration", ylabel = "Objective function", yscale = log10)
scatter!(ax, dprm.history.val)
fig

Compute sensitivities outside the optimization ​

We can also compute the sensitivities outside the optimization process. As our previous setup funciton only has a single parameter (the porosity), we instead switch to the setup_reservoir_dict_optimization function that can set up a set of "typical" tunable parameters for any reservoir model. This saves us the hassle of writing this function ourselves when we want to optimize e.g. permeability, porosity and well indices.

julia
dprm_case = setup_reservoir_dict_optimization(case)
free_optimization_parameters!(dprm_case)
dprm_grad = parameters_gradient_reservoir(dprm_case, mismatch_objective);
Optimization: Setting up adjoint storage.
Optimization: Finished setup in 10.294152822 seconds.

Plot the gradient of the mismatch objective with respect to the porosity ​

We see, as expected, that the gradient is largest in magnitude around the wells and near the front of the displacement.

julia
m = physical_representation(reservoir_domain(case.model))
plot_cell_data(m, dprm_grad[:model][:porosity])
(Scene (1600px, 900px):
  0 Plots
  2 Child Scenes:
    ├ Scene (1600px, 900px)
    └ Scene (1600px, 900px), Makie.Axis3(), MakieCore.Mesh{Tuple{GeometryBasics.Mesh{3, Float64, GeometryBasics.NgonFace{3, GeometryBasics.OffsetInteger{-1, UInt32}}, (:position, :normal), Tuple{Vector{GeometryBasics.Point{3, Float64}}, Vector{GeometryBasics.Vec{3, Float32}}}, Vector{GeometryBasics.NgonFace{3, GeometryBasics.OffsetInteger{-1, UInt32}}}}}})

Plot the sensitivities in the interactive viewer ​

If you are running the example yourself, you can now explore the sensitivities in the interactive viewer. This is useful for understanding how the model responds to changes in the parameters.

julia
plot_reservoir(case.model, dprm_grad[:model])

Example on GitHub ​

If you would like to run this example yourself, it can be downloaded from the JutulDarcy.jl GitHub repository as a script, or as a Jupyter Notebook

This example took 144.563134004 seconds to complete.

This page was generated using Literate.jl.