Adjoints and gradients

JutulDarcy models are fully differentiable, and can produce gradients with respect to any input parameter through the adjoint method. For more details, we refer you to the examples. It is recommended to go through the Introduction to history matching with gradients example to get an overview of the functionality.

Calculation of gradients

JutulDarcy.reservoir_sensitivities Function

julia

result, sens = reservoir_sensitivities(case::JutulCase, objective::Function; sim_arg = NamedTuple(), kwarg...)

Simulate a case and calculate parameter sensitivities with respect to an objective function on the form:

julia

obj(model, state, dt_n, n, forces_for_step_n)

The objective is summed up for all steps.

julia

reservoir_sensitivities(case::JutulCase, rsr::ReservoirSimResult, objective::Function; kwarg...)

Calculate parameter sensitivities with respect to an objective function on the form for a case and a simulation result from that case. The objective function is on the form:

julia

obj(model, state, dt_n, n, forces_for_step_n)

The objective is summed up for all steps.

julia

reservoir_sensitivities(case, rsr, objective; kwarg...)

Optimization interface

We use the functions from Jutul to free parameters before optimization.

Jutul.DictOptimization.free_optimization_parameter! Function

julia

free_optimization_parameter!(dopt, "parameter_name", rel_min = 0.01, rel_max = 100.0)
free_optimization_parameter!(dopt, ["dict_name", "parameter_name"], abs_min = -8.0, abs_max = 7.0)

Free an existing parameter for optimization in the DictParameters object. This will allow the parameter to be optimized through a call to optimize.

Nesting structures

If your DictParameters has a nesting structure, you can use a vector of strings or symbols to specify the parameter name, e.g. ["dict_name", "parameter_name"] to access the parameter located at ["dict_name"]["parameter_name"].

Setting limits

The limits can be set using the following keyword arguments:

abs_min: Absolute minimum value for the parameter. If not set, no absolute minimum will be applied.
abs_max: Absolute maximum value for the parameter. If not set, no absolute maximum will be applied.
rel_min: Relative minimum value for the parameter. If not set, no relative minimum will be applied.
rel_max: Relative maximum value for the parameter. If not set, no relative maximum will be applied.

For either of these entries it is possible to pass either a scalar, or an array. If an array is passed, it must have the same size as the parameter being set.

Note that if dopt.strict is set to true, at least one of the upper or lower bounds must be set for free parameters. If dopt.strict is set to false, the bounds are optional and the DictParameters object can be used to compute sensitivities, but the built-in optimization routine assumes that finite limits are set for all parameters.

Other keyword arguments

initial: Initial value for the parameter. If not set, the current value in dopt.parameters will be used.
scaler=missing: Optional scaler for the parameter. If not set, no scaling will be applied. Available scalers are:
- :log: Logarithmic scaling. This value uses shifts to avoid issues with zero values.
- :exp: Exponential scaling
- :linear: Linear scaling (scaling to bounds of values, guaranteeing values between between 0 and 1 for initial values.)
- linear_limits: Linear scaling with limits (scaling to bounds of values, guaranteeing values between between 0 and 1 for all values within the limits.)
- reciprocal: Reciprocal scaling
- log10: Base-10 logarithmic scaling
- log: Base-e logarithmic scaling without shifts
- A custom scaler object implementing the DictOptimizationScaler interface.
lumping=missing: Optional lumping array for the parameter. If not set, no lumping will be applied. The lumping array should have the same size as the parameter and contain positive integers. The lumping array defines groups of indices that should be lumped together, i.e. the same value will be used for all indices in the same group. The lumping array should contain all integers from 1 to the maximum value in the array, and all indices in the same group should have the same value in the initial parameter, otherwise an error will be thrown.

Jutul.DictOptimization.add_optimization_multiplier! Function

julia

add_optimization_multiplier!(dprm::DictParameters, name_of_target; abs_min = 0.2, abs_max = 5.0)
add_optimization_multiplier!(dprm, target1, target2, target3; abs_min = 0.2, abs_max = 5.0, initial = 2.0)

Add an optimization multiplier that acts on one or more targets to the DictParameters object. The multiplier will be optimized during the optimization process. All parameters with the same multiplier must have the same dimensions.

JutulDarcy.optimize_reservoir Function

julia

optimize_reservoir(dopt, objective, setup_fn = dopt.setup_function)

Perform optimization for reservoir models using the given DictParameters struct dopt, objective function objective, and setup function setup_fn.

Additional keyword arguments can be provided to customize the simulator setup and optimization process:

simulator_arg: Arguments for the simulator setup (default: (output_substates = true,)).
simulator: Custom simulator instance (default: missing).
config: Configuration for the simulator (default: missing).
deps: Dependencies for the optimization (default: :parameters_and_state0).

Notes

If you are optimizing forces (i.e. well constraints or boundary conditions), you need to set deps = :case to ensure that gradients are correctly computed. The same applies if you change the model itself in the setup function. The defaults of :parameters_and_state0 provide significant performance improvements in most cases where only parameters and initial state are changed.

Jutul.DictOptimization.DictParameters Type

julia

DictParameters(parameters)
DictParameters(parameters::AbstractDict, setup_function = missing;
        strict = true,
        verbose = true,
        active_type = Float64
    )

Set up a DictParameters object for optimization. Optionally, the setup function that takes an instance with the same keys as parameters together with a step_info dictionary can be provided. The setup function should return a JutulCase set up from the parameters in the Dict.

Optional keyword arguments:

strict: If true, the optimization will throw an error if any of the parameters are not set with at least one of the upper or lower bounds.
verbose: If true, the optimization will print information about the optimization process.
active_type: The type of the parameters that are considered active in the optimization. Defaults to Float64. This is used to determine which parameters are active and should be optimized. This means that all entries (and entries in nested dictionaries) of the parameters dictionary must be of this type or an array with this type as element type.

Functionality from Jutul.jl

This optimization interface is a wrapper around the Jutul function optimize with preselected options that are sensible for reservoir models. It can still be useful to look at the inner function to see additional supported arguments.

Jutul.DictOptimization.optimize Function

julia

optimized_dict = optimize(dopt, objective)
optimize(dopt::DictParameters, objective, setup_fn = dopt.setup_function;
    grad_tol = 1e-6,
    obj_change_tol = 1e-6,
    max_it = 25,
    opt_fun = missing,
    maximize = false,
    simulator = missing,
    config = missing,
    solution_history = false,
    backend_arg = (
        use_sparsity = false,
        di_sparse = true,
        single_step_sparsity = false,
        do_prep = true,
    ),
    kwarg...
)

Optimize parameters defined in a DictParameters object using the provided objective function. At least one variable has to be declared to be free using free_optimization_parameter! prior to calling the optimizer.

Arguments

dopt::DictParameters: Container with parameters to optimize
objective: The objective function to minimize (or maximize)
setup_fn: Function to set up the optimization problem. Defaults to dopt.setup_function

Keyword Arguments

grad_tol: Gradient tolerance for stopping criterion
obj_change_tol: Objective function change tolerance for stopping criterion
max_it: Maximum number of iterations
optimizer: Symbol defining the optimization algorithm to use. Available options are :lbfgs (default), :lbfgsb_qp and :lbfgsb (requires LBFGSB.jl to be imported)
opt_fun: Optional custom optimization function. If missing, L-BFGS will be used. Takes in a NamedTuple containing fields f, g, x0, min, max. Here, f(x) returns the objective function value at x, g(dFdx, x) fills dFdx with the gradient at x, x0 is the initial guess, and min and max are the lower and upper bounds, respectively. The functions u = F.scale(x) and x = F.descale(u) can be used to convert between scaled and unscaled variables. Nominally, the initial values are scaled to the unit cube and the solution must thus be unscaled before usage. Gradients and internal scaling/descaling is automatically handled.
maximize: Set to true to maximize the objective instead of minimizing
gradient_scaling: If true, internally scales the objective gradient according to the initial 2-norm of the gradient. If a Float64 value is provided, that value is used as a global scaling factor for the gradient. The internal gradient and objective value is divided by the chosen scaling.
simulator: Optional simulator object used in forward simulations
config: Optional configuration for the setup
solution_history: If true, stores all intermediate solutions
deps: One of :case, :parameters, :parameters_and_state0. Defines the dependencies for the adjoint computation. See notes for more details.
backend_arg: Options for the autodiff backend:
- use_sparsity: Enable sparsity detection for the objective function
- di_sparse: Use sparse differentiation
- single_step_sparsity: Enable single step sparsity detection (if sparsity does not change during timesteps). This means that the solver will assume that the sparsity pattern will be determined entirely by the first and last steps of the simulation. Alternatively, this can be set to :unique_forces to use the sparsity pattern determined by all unique force terms in the solve, :firstlast to only use the first and last time steps or :allsteps to use all time steps (the latter is equivalent to setting use_sparsity to true).
- do_prep: Perform preparation step
- output_path: If provided, the optimization results will be stored in the given path as a JLD2 file named final.jld2, with intermediate steps being stored as step_1.jld2, step_2.jld2, etc if solution_history is enabled.
- randomized_start: If true, the initial guess will be randomized within the provided limits instead of using the initial values in dopt.parameters.

Returns

The optimized parameters as a dictionary.

Notes

The function stores the optimization history and optimized parameters in the input dopt object.
If solution_history is true or :x, intermediate solutions are stored in dopt.history.solutions. If it is set to :full, the full states are also copied and stored for each iteration. This can use a lot of memory for large simulations.
The default optimization algorithm is L-BFGS with box constraints.

Type of dependencies in deps

The deps argument is used to set the type of dependency the case setup function has on the active optimization parameters. The default, :case, is fully general and allows dependencies on everything contained within the case instance. This can be slow, however, as the setup function must be called for every time-step. If you know that the model instance and forces are independent of the active parameters, you can use deps = :parameters_and_state0. If there is no dependence on state0, you can set deps = :parameters. This can substantially speed up the optimization process, but as there is no programmatic verification that this assumption is true, it should be used with care.

This interface is dependent on the model supporting use of vectorize_variables! and devectorize_variables! for state0/parameters, which should be the case for most Jutul models.

Jutul.LBFGS.unit_box_bfgs Function

julia

unit_box_bfgs(u0, f; kwargs...)

Iterative line search optimization using BFGS intended for scaled problems where 0 ≤ u ≤ 1 and f ~ O(1).

This is a port of the MRST function unitBoxBFGS, relicensed under MIT.

Arguments

u0: Initial guess vector of length n with 0 ≤ u0 ≤ 1, must be feasible with respect to additional constraints
f: Function handle that returns tuple (v, g) where:
- v: objective value
- g: objective gradient vector of length n

Keywords

maximize::Bool=false: Set to true to maximize objective
step_init::Float64=NaN: Initial step (gradient scaling). If NaN, uses max_initial_update/max(|initial gradient|)
time_limit::Float64=Inf: Time limit for optimizer.
max_initial_update::Float64=0.05: Maximum initial update step
history::Any=nothing: For warm starting based on previous optimization (requires output_hessian=true)

Stopping Criteria

grad_tol::Float64=1e-3: Absolute tolerance of inf-norm of projected gradient
obj_change_tol::Float64=5e-4: Absolute objective update tolerance
max_it::Int=25: Maximum number of iterations

Line Search Options

line_searchmax_it::Int=5: Maximum number of line-search iterations
stepIncreaseTol::Float64=10: Maximum step increase factor between line search iterations
wolfe1::Float64=1e-4: Objective improvement condition
wolfe2::Float64=0.9: Gradient reduction condition

Hessian Approximation

use_bfgs::Bool=true: Use BFGS (false for pure gradient search)
limited_memory::Bool=true: Use L-BFGS instead of full Hessian approximations
lbfgs_num::Int=5: Number of vector-pairs stored for L-BFGS
lbfgs_strategy::Symbol=:dynamic: Strategy for L-BFGS (:static or :dynamic)

Linear Constraints

lin_eq::NamedTuple{(:A,:b)}: Linear equality constraints A*u=b
lin_ineq::NamedTuple{(:A,:b)}: Additional linear inequality constraints A*u≤b
enforce_feasible::Bool=false: Attempt to repair constraint violations by projection

Display Options

plotEvolution::Bool=true: Plot optimization progress
logPlot::Bool=false: Use logarithmic y-axis for objective plotting
output_hessian::Bool=false: Include Hessian approximation in history output

Returns

v: Optimal or best objective value
u: Control/parameter vector corresponding to v
history: Named tuple containing iteration history with fields:
- val: objective values
- u: control/parameter vectors
- pg: projected gradient norms
- alpha: line-search step lengths
- lsit: number of line-search iterations
- lsfl: line-search flags
- hess: Hessian inverse approximations (if requested)

History matching module

The history matching module makes it easy to define objective functions from typical reservoir observations (rates, pressures and fractions per well as time series).

JutulDarcy.HistoryMatching.history_match_objective Function

julia

history_match_objective(case::JutulCase)
history_match_objective(case::JutulCase, res::ReservoirSimResult)
history_match_objective(case::JutulCase, states, summary_to_match_against)

Set up a history match objective. The history match objective computes the mismatch between simulation results and observed data for a given case. The objective can be either global or sum, depending on the is_global keyword argument. Sum is generally faster to compute, but may be more sensitive to sharp gradients. Cumulative matches are only available for the global option.

The objective is written as (for each well that is being matched):

julia

`f(x) = sum_i w_i * (sim_i(x) - obs_i)^n` (n defaults to 2)

where w_i is a weight for the i-th well match, sim_i(x) is the simulated value for the i-th step well match, and obs_i is the observed value for the i-th well match. The sum is taken over all well matches in the history match object. You must first instantiate a history match object with history_match_objective, and then add well matches with match_injectors! and match_producers! where you can also set weights per step.

JutulDarcy.HistoryMatching.match_well! Function

julia

match_well!(hm_obj, well_name, quantity; is_injector = true)
match_well!(hm_obj, "WellName", "WBHP"; weight = 3.0, is_injector = true)

Add a well match to the history match object hm_obj for the well with name well_name and quantity quantity. Additional keyword arguments:

weight: Weighting factor for the well match. Can be a scalar or a vector with length equal to the number of report steps in the simulation case. Default is 1.0. If you are matching multiple wells and quantities (especially rates and pressures), you may want to adjust the weights to balance the contributions to the overall objective.
is_injector: Set to true if the well is an injector, false if it is a producer. Mandatory. Use match_injectors! or match_producers! for convenience.
data: Optionally provide observation data as a vector with values per report step or a function (t -> value_at_t). If missing, data is taken from the case summary embedded in the history matching object.
scale: Optionally provide a scaling factor for the well match. If missing, a default scale is used based on the quantity and phase to make the value roughly dimensionless.
exponent: Exponent for the mismatch calculation. Default is 2.0 (squared difference).

JutulDarcy.HistoryMatching.match_injectors! Function

julia

match_injectors!(obj, "WBHP", weight = 3.0)

Match a quantity for all injectors in the history match object obj. See match_well! for details on the keyword arguments.

Possible quantities for injectors:

:bhp or "WBHP": Bottom hole pressure
:rate or "RATE": Total volumetric injection rate at standard conditions
:orat or "WOIR": Oil volumetric injection rate at standard conditions
:wrat or "WWIR": Water volumetric injection rate at standard conditions
:grat or "WGIR": Gas volumetric injection rate at standard conditions

JutulDarcy.HistoryMatching.match_producers! Function

julia

match_producers!(obj, "ORAT", weight = 3.0)

Match a quantity for all producers in the history match object obj. See match_well! for details on the keyword arguments.

Possible quantities for producers:

:bhp or "WBHP": Bottom hole pressure
:rate or "RATE": Total volumetric production rate at standard conditions
:orat or "WOPR": Oil volumetric production rate at standard conditions
:wrat or "WWPR": Water volumetric production rate at standard conditions
:grat or "WGPR": Gas volumetric production rate at standard conditions
:lrat or "WLPR": Liquid volumetric production rate at standard conditions
:wcut or "WWCT": Water cut (ratio of water to total liquid production)
:gor or "WGOR": Gas-oil ratio (ratio of gas to oil production)
:wgr or "WWGR": Water-gas ratio (ratio of water to gas production)
:glr or "WGLR": Gas-liquid ratio (ratio of gas to liquid production)
:cumulative_oil or "WOPT": Cumulative oil production at standard conditions (requires global objective)
:cumulative_gas or "WGPT": Cumulative gas production at standard conditions (requires global objective)
:cumulative_water or "WWPT": Cumulative water production at standard conditions (requires global objective)
:cumulative_liquid or "WLPT": Cumulative liquid production at standard conditions (requires global objective)

Utilities

JutulDarcy.setup_reservoir_dict_optimization Function

julia

setup_reservoir_dict_optimization(case::JutulCase)
setup_reservoir_dict_optimization(case::JutulCase;
    use_trans = false,
    use_pore_volume = false,
    use_multipliers = false,
    strict = false,
    verbose = true,
    do_copy = true,
    parameters = Symbol[],
    kwarg...
)

Set up a DictParameters struct for reservoir model optimization with a setup function. The function extracts relevant model and parameter data from the provided JutulCase instance, puts them in a dict for use in optimization, and defines a setup function to reconstruct the JutulCase from modified parameters.

Options:

use_trans: If true, include transmissibilities as optimization parameters instead of permeability.
use_pore_volume: If true, include pore volumes as optimization parameters instead of porosity.
use_multipliers: If true, use multipliers for permeability/transmissibility, well indices and pore volume instead of absolute values.

JutulDarcy.well_mismatch Function

julia

well_mismatch(qoi, wells, model_f, states_f, model_c, state_c, dt, step_info, forces; <keyword arguments>)

Compute well mismatch for a set of qoi's (well targets) and a set of well symbols.

JutulDarcy.compute_well_qoi Function

julia

compute_well_qoi(model::MultiModel, state, forces, well::Symbol, target::Union{WellTarget, Type})

Compute the quantity of interest (QoI) for a specified well in a reservoir simulation.

Arguments

model::MultiModel: The simulation model (from setup_reservoir_model).
state: The current state of the simulation.
forces: The forces applied in the simulation.
well::Symbol: The symbol representing the well for which the QoI is computed.
target::Union{WellTarget, Type, Symbol}: The target type (as type or symbol) or specific well target for the QoI computation.

Possible targets (as symbols):

:gor (Gas-to-Oil Ratio)
:wcut (Water Cut)
:temperature (Well Temperature)
:mass_rate (Total Surface Mass Rate)
:bhp (Bottom Hole Pressure)
:lrat (Surface Liquid Rate)
:wrat (Surface Water Rate)
:orat (Surface Oil Rate)
:grat (Surface Gas Rate)
:rate (Total volumetric rate at surface conditions)

Returns

The computed QoI for the specified well.

JutulDarcy.npv_objective Function

julia

npv_objective(model, state, dt, step_info, forces;
    timesteps,
    injectors,
    producers,
    oil_price = 60.0,
    gas_price = 10.0,
    water_price = -3.0,
    water_cost = 5.0,
    oil_cost = oil_price,
    gas_cost = gas_price,
    liquid_unit = si_unit(:stb),
    gas_unit = si_unit(:kilo)*si_unit(:feet)^3,
    discount_rate = 0.025,
    discount_unit = si_unit(:year),
    scale = 1.0
)

Evaluate the contribution to net-present-value for a given step, with the given costs and prices for oil, gas, and water and discount rate. Costs are assumed to be the cost of injecting a given fluid and prices are the revenue from producing the corresponding fluids. Prices and costs can be set to negative (e.g. to account for water being a cost when produced).

Setting maximize to false will make the objective function negative-valued, with larger negative values corresponding to better results. This is useful for some optimizers.