GitHub

History matching of simulation models

We demonstrate history matching of a Direct Column Breakthrough (DCB) simulation model in Mocca. We leverage the powerful and flexible optimization functionality of Jutul to set up and perform the history matching. For more details about the DCB modelling, see the Haghpanah DCB example.

Import necessary modules

import Jutul
import Mocca

We create a function for setting up new simulation cases from the value of the parameter we wish to tune

function setup_case(prm, step_info=missing)

    param_dict_symb = Dict(Symbol(k) => v for (k, v) in prm);
    RealT = valtype(param_dict_symb);
    constants, info = Mocca.parse_input(Mocca.haghpanah_DCB_input(); typeT=RealT)

    for (k, v) in param_dict_symb
        print(k)
        print(v)
        setproperty!(constants, Symbol(k), v)
    end
    case,  = Mocca.setup_mocca_case(constants, info)

    return case
end;

Create synthetic reference data

constants_ref, = Mocca.parse_input(Mocca.haghpanah_DCB_input(); typeT=Float64)

prm_ref = Dict("v_feed" => constants_ref.v_feed);
case_ref = setup_case(prm_ref);
v_feed0.37

Configure simulator which will be used in the history matching

timestep_selector_cfg = (y=0.01, Temperature=10.0, Pressure=10.0)
sim, cfg = Mocca.setup_process_simulator(case_ref.model, case_ref.state0, case_ref.parameters;
    timestep_selector_cfg = timestep_selector_cfg,
    initial_dt = 1.0,
    output_substates = true,
    info_level = -1
);

Run reference simulation to generate and generate "ground truth" data from the result

states, timesteps_out = Mocca.simulate_process(case_ref;
    simulator = sim,
    config = cfg
);

times_ref = cumsum(timesteps_out)
total_time = times_ref[end]
last_cell_idx = Jutul.number_of_cells(case_ref.model.domain)
qCO2_ref = map(s -> getindex(s[:AdsorbedConcentration], 1, last_cell_idx), states)
qCO2_ref_by_time = Jutul.get_1d_interpolator(times_ref, qCO2_ref);

Setting up and solving the optimization problem

We define a suitable objective function to quantify the match between our simulations and the reference solution. Here we choose deviation of adsorbed CO2 in the last grid cell.

function objective_function(model, state, dt, step_info, forces)
    current_time = step_info[:time]
    q_co2 = getindex(state[:AdsorbedConcentration], 1, last_cell_idx)
    q_co2_ref = qCO2_ref_by_time(current_time)
    v = dt/total_time*(q_co2 - q_co2_ref)^2
    return v
end;

Perturb the known parameter $v_{feed}$ to form our initial guess for the optimization

prm_guess = Dict("v_feed" => constants_ref.v_feed+0.2);

Activate $v_{feed}$ as a free parameter

dprm = Jutul.DictOptimization.DictParameters(prm_guess)
Jutul.DictOptimization.free_optimization_parameter!(dprm, "v_feed"; rel_min = 0.001, rel_max = 100.0)
DictParameters with 1 parameters (1 active), and 0 multipliers:
Active optimization parameters
┌────────┬───────────────┬───────┬─────────┬──────┐
│   Name │ Initial value │ Count │     Min │  Max │
├────────┼───────────────┼───────┼─────────┼──────┤
│ v_feed │ 0.57          │     1 │ 0.00057 │ 57.0 │
└────────┴───────────────┴───────┴─────────┴──────┘
No inactive optimization parameters.
No multipliers set.

Run the optimization

prm_opt = Jutul.DictOptimization.optimize(dprm, objective_function, setup_case;
    config = cfg,
    max_it = 10,
    obj_change_tol = 1e-6,
    solution_history = true
);
Optimization: Starting calibration of 1 parameters.
Optimization: Setting up adjoint storage.
Optimization: Finished setup in 43.098538081 seconds.
Optimization: Adjoint solve took 29.350653884 seconds.

Optimization: Objective #1: 7.61286e+05, gradient 2-norm: 4.71893e+06
It.  | Objective  | Proj. grad | Linesearch-its
-----------------------------------------------
   0 | 1.6133e-01 | 5.6999e+01 | -
Optimization: Adjoint solve took 0.515855883 seconds.

Optimization: Objective #2: 6.36994e+06 (f/f0=8.367e+00), gradient 2-norm: 1.21894e-04
Optimization: Adjoint solve took 7.358947922 seconds.

Optimization: Objective #3: 6.05491e+05 (f/f0=7.954e-01), gradient 2-norm: 5.00075e+06
Optimization: Adjoint solve took 7.337086057 seconds.

Optimization: Objective #4: 1.24529e+06 (f/f0=1.636e+00), gradient 2-norm: 6.30573e+06
Optimization: Adjoint solve took 6.446183328 seconds.

Optimization: Objective #5: 1.60276e+05 (f/f0=2.105e-01), gradient 2-norm: 4.25403e+06
Optimization: Adjoint solve took 0.54129832 seconds.

Optimization: Objective #6: 6.36994e+06 (f/f0=8.367e+00), gradient 2-norm: 1.21894e-04
LBFGS: Line search unable to succeed in 5 iterations ...
Optimization: Adjoint solve took 6.393882964 seconds.

Optimization: Objective #7: 1.60276e+05 (f/f0=2.105e-01), gradient 2-norm: 4.25403e+06
LBFGS: Resetting 'm' to number of parameters: m = 1
   1 | 3.3964e-02 | 5.6999e+01 | 5
Optimization: Adjoint solve took 0.56305323 seconds.

Optimization: Objective #8: 6.36994e+06 (f/f0=8.367e+00), gradient 2-norm: 1.21894e-04
Optimization: Adjoint solve took 6.360798124 seconds.

Optimization: Objective #9: 9.17200e+04 (f/f0=1.205e-01), gradient 2-norm: 3.55824e+06
   2 | 1.9437e-02 | 5.1384e+01 | 2
Optimization: Adjoint solve took 6.496560096 seconds.

Optimization: Objective #10: 1.56515e+05 (f/f0=2.056e-01), gradient 2-norm: 6.37615e+06
Optimization: Adjoint solve took 6.51306456 seconds.

Optimization: Objective #11: 5.98917e+01 (f/f0=7.867e-05), gradient 2-norm: 1.18220e+05
   3 | 1.2692e-05 | 4.2980e+01 | 2
Optimization: Adjoint solve took 7.307420457 seconds.

Optimization: Objective #12: 1.30250e+01 (f/f0=1.711e-05), gradient 2-norm: 4.91360e+04
   4 | 2.7602e-06 | 1.4280e+00 | 1
Optimization: Adjoint solve took 6.442963115 seconds.

Optimization: Objective #13: 4.89064e-01 (f/f0=6.424e-07), gradient 2-norm: 1.02865e+04
   5 | 1.0364e-07 | 5.9351e-01 | 1
Optimization: Adjoint solve took 6.397948268 seconds.

Optimization: Objective #14: 1.31567e-02 (f/f0=1.728e-08), gradient 2-norm: 1.63327e+03
   6 | 2.7881e-09 | 1.2425e-01 | 1
Optimization: Finished in 183.736794389 seconds.

We can see a clear reduction of the objective function value throughout the optimization iterations, indicating a close match between the reference solution and our simulation.

f = Mocca.plot_optimization_history(dprm)
Example block output

We can look at the optimization result:

dprm
DictParameters with 1 parameters (1 active), and 0 multipliers:
Active optimization parameters
┌────────┬───────────────┬───────┬─────────┬──────┬─────────────────┬────────┐
│   Name │ Initial value │ Count │     Min │  Max │ Optimized value │ Change │
├────────┼───────────────┼───────┼─────────┼──────┼─────────────────┼────────┤
│ v_feed │ 0.57          │     1 │ 0.00057 │ 57.0 │ 0.37            │ -35.0% │
└────────┴───────────────┴───────┴─────────┴──────┴─────────────────┴────────┘
No inactive optimization parameters.
No multipliers set.

and see that the value matches the reference parameter value

constants_ref.v_feed
0.37

Check size of objective mismatch

dprm.history.objectives
14-element Vector{Float64}:
 761285.8833491696
      6.369944548281407e6
 605490.5988767222
      1.2452850490274746e6
 160276.16543350072
      6.369944548281407e6
 160276.16543350072
      6.369944548281407e6
  91720.04717463219
 156514.61420140337
     59.891651093548234
     13.024970004656637
      0.4890638735177114
      0.013156736822509599

Example on GitHub

If you would like to run this example yourself, it can be downloaded from the Mocca.jl GitHub repository.

This page was generated using Literate.jl.