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 42.714645379 seconds.
Optimization: Adjoint solve took 28.423541359 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.521956028 seconds.

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

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

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

Optimization: Objective #5: 1.60276e+05 (f/f0=2.105e-01), gradient 2-norm: 4.25403e+06
Optimization: Adjoint solve took 0.513817548 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.151634595 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.507065942 seconds.

Optimization: Objective #8: 6.36994e+06 (f/f0=8.367e+00), gradient 2-norm: 1.21894e-04
Optimization: Adjoint solve took 6.154278612 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.300011735 seconds.

Optimization: Objective #10: 1.56515e+05 (f/f0=2.056e-01), gradient 2-norm: 6.37615e+06
Optimization: Adjoint solve took 6.203241298 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 6.216263184 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.17139298 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.18030875 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 177.664140645 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.