examples.history_match_gradient

  1"""The gradient of a production-history misfit wrt. $\\log K$, and a few descent steps.
  2
  3The history-matching case, in its simplest form. A *truth* permeability field
  4(log-normal, ref `examples.random_fields`) produces the *observations*: the water cut at each of
  5the four producers of a five-spot, at every time step. The *prior* guess is
  6homogeneous, $\\log K = 0$. The objective is the mean squared error of the
  7prior's production history against the observations, and
  8`minires.tlm.adjoint` gives its gradient with respect to every cell's
  9$\\log K$ -- checked against a finite difference in a random direction -- for
 10about the cost of one more simulation, however many cells there are. That is
 11what makes gradient-based history matching feasible at all, and a few steps of
 12steepest descent (each with a coarse line search) are taken to show it works:
 13the misfit falls, and the update goes the way of the truth.
 14
 15The seeds are those of a data misfit
 16(ref the "Seeding" section of `minires.tlm`): for each observed time and
 17producer, the residual, weighted by the derivative of the observation operator
 18-- here $ f_w'(s) $ in the producer's cell, via `minires.fluids.Fluid.fractional_flow`.
 19
 20In the figure:
 21
 22- Top left: the truth $\\log K$, which the observations come from.
 23- Top middle: the gradient of the misfit at the prior. Since the prior is
 24  homogeneous, this is the whole of the first update (with a minus sign): it
 25  is negative (raise $K$) where the truth is permeable *and* the water reaches
 26  a producer through it, positive where the truth is tight -- but only where
 27  the production data are sensitive, i.e. along the flow paths. Hence its
 28  correlation with the truth over the whole field is weak (about 0.2; it is
 29  recorded in the digest): the data are blind to most of the field.
 30- Top right: the estimate after the descent steps (on its own colour scale:
 31  the update is far smaller than the truth's variations). It moves towards
 32  the truth along the flow paths and not elsewhere, the data saying nothing
 33  about the rest: the classic ill-posedness that a prior (regularization)
 34  would address.
 35- Bottom: the production histories -- observations, prior, and the final
 36  estimate -- and the misfit per iteration.
 37
 38.. note:: Steepest descent is used for its simplicity, not its merit.
 39
 40    A Gauss--Newton or quasi-Newton method, and a prior term, would do far
 41    better; the point here is the *gradient*, which any of them needs.
 42"""
 43
 44import matplotlib.pyplot as plt
 45import numpy as np
 46
 47from minires import ResSim
 48from minires.geostat import gaussian_fields
 49from minires.plotting import show
 50from minires.tlm import adjoint
 51
 52rng = np.random.default_rng(1)  # Reproducibility (the values are regression tested)
 53
 54## Model: a five-spot
 55wells = [
 56    dict(xy=[.5, .5], rate=+1  , name="inj"),
 57    dict(xy=[1 , 1 ], rate=-.25, name="NE"),
 58    dict(xy=[0 , 1 ], rate=-.25, name="NW"),
 59    dict(xy=[0 , 0 ], rate=-.25, name="SW"),
 60    dict(xy=[1 , 0 ], rate=-.25, name="SE"),
 61]
 62grid: dict = dict(Lx=1, Ly=1, Nx=32, Ny=32)
 63dt, nSteps = .05, 30
 64
 65
 66def new_model(logK):
 67    model = ResSim(**grid, wells=wells)
 68    model.K = np.exp(logK)  # isotropic: broadcast to both components
 69    return model
 70
 71
 72model = new_model(0)  # the prior: homogeneous
 73S0 = np.zeros(model.Nxy)
 74producers = model.wells.names[1:]
 75prd = model.xy2ind(*model.wells.xy[1:].T)  # their cells
 76
 77
 78def water_cut(model, SS):
 79    """`(nSteps, nPrd)` water cut at each producer, at times `1..nSteps`."""
 80    return np.array([model.fluid.fractional_flow(S)[prd] for S in SS[1:]])
 81
 82
 83## The truth, and the observations it produces
 84logK_true = gaussian_fields(model.mesh, r=.15, rng=rng)[0].reshape(model.shape)
 85# logK_true = 3 * smooth(smooth(rng.standard_normal(model.shape)))  # the old way
 86truth = new_model(logK_true)
 87obs = water_cut(truth, truth.sim(dt, nSteps, S0, pbar=False)[0])
 88
 89
 90## The objective, and its gradient by the adjoint
 91def misfit(logK, gradient=False):
 92    """Mean squared error of the production history of `logK` against `obs`."""
 93    model = new_model(logK)
 94    SS, PP = model.sim(dt, nSteps, S0, pbar=False)
 95    fw = water_cut(model, SS)
 96    residual = fw - obs
 97    J = (residual**2).mean()
 98    if not gradient:
 99        return J
100    # Seed: ∂J/∂s_k[i] = 2/nObs * residual * f_w'(s), at the producers, each time
101    dJ_dSS = np.zeros_like(SS)
102    for k in range(1, nSteps + 1):
103        dJ_dSS[k, prd] = 2 / obs.size * residual[k - 1] * model.fluid.dfractional_flow(SS[k])[prd]
104    G = adjoint(model, dt, SS, PP, dJ_dSS).logK.sum(0)  # isotropic ⇒ sum the components
105    return J, G, fw
106
107
108logK = np.zeros(model.shape)  # the prior
109J0, G0, fw_prior = misfit(logK, gradient=True)
110
111## Check: a finite difference in a random direction of log K
112direction = rng.standard_normal(model.shape)
113eps = 1e-5
114fd = (misfit(logK + eps*direction) - misfit(logK - eps*direction)) / (2*eps)
115directional = (G0 * direction).sum()
116assert abs(fd - directional) < 1e-4 * abs(directional), (fd, directional)
117
118## A few steps of steepest descent, each with a coarse line search
119nIter = 4
120step_sizes = [.1, .2, .4, .8]  # in units of max |Δ log K| per step
121JJ = [J0]
122J, G = J0, G0
123for _ in range(nIter):
124    d = -G / abs(G).max()  # the (normalized) descent direction
125    trials = [misfit(logK + a*d) for a in step_sizes]
126    best = int(np.argmin(trials))
127    if trials[best] >= J:
128        break  # no step size improves: stop
129    logK = logK + step_sizes[best] * d
130    J, G, fw = misfit(logK, gradient=True)
131    JJ.append(J)
132fw_final = water_cut(new_model(logK), new_model(logK).sim(dt, nSteps, S0, pbar=False)[0])
133corr = np.corrcoef(-G0.ravel(), logK_true.ravel())[0, 1]
134
135## Plot
136fig, axs = plt.subplots(num="History-match gradient", clear=True, ncols=3, nrows=2, figsize=(13, 8))
137
138kws = dict(cmap="viridis", wells="color", finalize=False)
139model.plt_field(axs[0, 0], logK_true, title="Truth, $\\log K$",
140                levels=np.linspace(-3, 3, 19), cticks=np.arange(-3, 4), **kws)
141m = abs(logK).max()  # NB: its own scale -- the update is far smaller than the truth
142model.plt_field(axs[0, 2], logK, title=f"Estimate, after {len(JJ) - 1} descent steps",
143                levels=np.linspace(-m, m, 19), cticks=[-m, 0, m], **kws)
144m = abs(G0).max()
145model.plt_field(axs[0, 1], G0, title="$∂J/∂\\log K$ at the prior ($\\log K = 0$)",
146                cmap="RdBu_r", levels=np.linspace(-m, m, 21), cticks=[-m, 0, m],
147                wells="color", finalize=False)
148axs[0, 1].text(.02, .02, f"corr(−∂J/∂log K, truth) = {corr:.2f}",
149               transform=axs[0, 1].transAxes, fontsize=8)
150
151tt = dt * np.arange(1, nSteps + 1)
152for ax, fw, title in [(axs[1, 0], fw_prior, "prior"), (axs[1, 1], fw_final, "estimate")]:
153    for i, name in enumerate(producers):
154        ax.plot(tt, obs[:, i], "*", c=f"C{i}", label=f"{name} obs.")
155        ax.plot(tt, fw[:, i], "-", c=f"C{i}", label=f"{name} {title}")
156    ax.set(title=f"Production history: {title}", xlabel="Time", ylabel="Water cut",
157           ylim=(-.02, 1))
158    ax.legend(loc="upper left", ncol=2, fontsize=8)
159
160ax = axs[1, 2]
161ax.semilogy(JJ, "o-")
162ax.set(title="Misfit (MSE)", xlabel="Iteration")
163ax.xaxis.get_major_locator().set_params(integer=True)
164
165fig.tight_layout()
166
167# Regression values, checked by `tests/test_examples.py`.
168__digest__ = dict(misfit      = JJ,
169                  gradient    = G0,
170                  directional = [directional, fd],
171                  corr        = [corr],
172                  logK_final  = logK)
173
174if __name__ == "__main__":
175    show()