examples.water_cut_gradient

  1"""The gradient of the water cut at one producer, at one time, wrt. the field of $\\log K$ and the BHP schedule.
  2
  3A five-spot: one injector at the centre, at a fixed rate, and four producers in
  4the corners, on BHP control (so their rates are outcomes, and the flow splits
  5among them as the heterogeneous -- log-normal, ref `examples.random_fields` --
  6permeability dictates). The objective is the water cut at the NE producer at a time index
  7just after its breakthrough, while it is rising fast.
  8`minires.tlm.adjoint` returns its gradient with respect to every cell's
  9$\\log K$ (and the initial state) *and* with respect to every producer's BHP
 10at every time step, at the cost of about one more simulation -- here each is
 11checked against a finite difference in a random direction, which agrees to
 12about 1e-6 (relative; the floor is set by the kinks of the discrete map, ref
 13the "What is differentiated" section of `minires.tlm`).
 14
 15The water cut at a producer, $ f_w(s) $ in its cell (the fraction of water in
 16what it produces), is a function of the saturation there alone, so the adjoint
 17is seeded by the single entry $ ∂J/∂s_k[i_\\mathrm{prd}] = f_w'(s) $ -- which
 18`minires.fluids.Fluid.fractional_flow` supplies.
 19
 20In the figure:
 21
 22- Top left: the $\\log K$ field, with the wells.
 23- Top right: the water cut of each producer over time, the objective marked.
 24- Bottom left: the gradient wrt. $\\log K$. It is positive along the flow path
 25  from the injector to the NE producer -- more permeable rock there brings the
 26  water sooner -- and negative along the paths to the *other* producers, and
 27  around the NE path: more permeable rock there diverts the water from the NE
 28  well, or lets it sweep a wider area (later arrival at the well itself).
 29  Outside the drainage area of the injector-NE pair, it vanishes: what happens
 30  there has not yet had time to affect the NE well.
 31- Bottom right: the gradient wrt. the BHP of each producer, per time step
 32  (`tlm.Gradient.bhp`). *Lowering* the NE well's BHP draws more of the
 33  injected water its way, hence an earlier breakthrough and a higher water cut
 34  at the objective's time: its gradient is negative, throughout. The other
 35  three's are positive: drawing more to them starves the NE well. That they
 36  are of a size, the neighbours' (NW, SE) no smaller than the opposite
 37  corner's (SW), is a matter of the heterogeneity: drawing more flow towards a
 38  neighbour also reshapes the stream tube to NE, which may help or hinder (for
 39  another field, these may even change sign over time). All of them vanish
 40  from the objective's time step on -- a control cannot affect what came
 41  before it. Their sum over time is the gradient with respect to a
 42  constant-in-time BHP.
 43
 44.. note:: The gradient is with respect to the *isotropic* $\\log K$.
 45
 46    `model.K` holds both components, so the adjoint returns a gradient for
 47    each; as they are set equal here, the gradient of the single field is
 48    their sum, ref `tlm.Gradient.logK`.
 49
 50.. note:: The well indices are held fixed, as in the adjoint.
 51
 52    They are computed (from the permeability at the wells) once, when the
 53    wells are configured, and stored (`Wells.WI`); the finite differences
 54    below reuse them, as `tlm` assumes (ref its "What is differentiated").
 55"""
 56
 57from dataclasses import replace
 58
 59import matplotlib.pyplot as plt
 60import numpy as np
 61
 62from minires import Grid2D, ResSim
 63from minires.geostat import gaussian_fields
 64from minires.plotting import show
 65from minires.tlm import adjoint
 66
 67rng = np.random.default_rng(1)  # Reproducibility (the values are regression tested)
 68
 69## Model: a five-spot on heterogeneous permeability
 70grid: dict = dict(Lx=1, Ly=1, Nx=32, Ny=32)
 71mesh = Grid2D(**grid).mesh
 72logK = gaussian_fields(mesh, r=.15, rng=rng)[0].reshape(mesh[0].shape)
 73# logK = 3 * smooth(smooth(rng.standard_normal((grid["Nx"], grid["Ny"]))))  # the old way
 74model = ResSim(**grid, K=np.exp(logK), wells=[  # isotropic: K broadcast to both components
 75    dict(xy=[.5, .5], rate=+1, name="inj"),
 76    dict(xy=[1 , 1 ], bhp=0, rw=1e-3, name="NE"),
 77    dict(xy=[0 , 1 ], bhp=0, rw=1e-3, name="NW"),
 78    dict(xy=[0 , 0 ], bhp=0, rw=1e-3, name="SW"),
 79    dict(xy=[1 , 0 ], bhp=0, rw=1e-3, name="SE"),
 80])
 81
 82dt, nSteps = .05, 30
 83S0 = np.zeros(model.Nxy)
 84SS, PP = model.sim(dt, nSteps, S0, pbar=False)
 85
 86## The objective: water cut at the NE producer at time `k`
 87producers = model.wells.names[1:]
 88prd = model.xy2ind(*model.wells.xy[1:].T)  # their cells
 89
 90
 91def water_cut(model, SS):
 92    """`(nSteps+1, nPrd)` water cut at each producer, for each stored time."""
 93    return np.array([model.fluid.fractional_flow(S)[prd] for S in SS])
 94
 95
 96fw = water_cut(model, SS)
 97well, k = 0, 13  # NE, just after breakthrough (fw ≈ .4)
 98J = fw[k, well]
 99
100## Its gradient, by the adjoint
101dJ_dSS = np.zeros_like(SS)
102dJ_dSS[k, prd[well]] = model.fluid.dfractional_flow(SS[k])[prd[well]]  # f_w'(s)
103grad = adjoint(model, dt, SS, PP, dJ_dSS)
104G = grad.logK.sum(0)  # isotropic ⇒ sum the components
105G_bhp = grad.bhp      # (nComp, nSteps); zero for the (rate-controlled) injector
106
107
108## Check: finite differences in random directions of log K and of the BHP schedule
109def J_of(logK, bhp):
110    """The objective of a model with `logK`, its producers at the schedule `bhp`."""
111    m = ResSim(**grid, K=np.exp(logK), wells=replace(model.wells, bhp=bhp))  # same `WI`
112    return water_cut(m, m.sim(dt, nSteps, S0, pbar=False)[0])[k, well]
113
114
115eps = 1e-5
116bhp = model.wells.bhp  # (nComp, 1): constant in time
117d_logK = rng.standard_normal(model.shape)
118d_bhp = np.where(np.isfinite(bhp), rng.standard_normal((model.nComp, nSteps)), 0)
119fd_logK = (J_of(logK + eps*d_logK, bhp) - J_of(logK - eps*d_logK, bhp)) / (2*eps)
120fd_bhp = (J_of(logK, bhp + eps*d_bhp) - J_of(logK, bhp - eps*d_bhp)) / (2*eps)
121directional_logK = (G * d_logK).sum()
122directional_bhp = (G_bhp * d_bhp).sum()
123assert abs(fd_logK - directional_logK) < 1e-4 * abs(directional_logK), (fd_logK, directional_logK)
124assert abs(fd_bhp - directional_bhp) < 1e-4 * abs(directional_bhp), (fd_bhp, directional_bhp)
125
126## Plot
127fig, axs = plt.subplots(num="Water-cut gradient", clear=True, ncols=2, nrows=2, figsize=(10, 8),
128                        gridspec_kw={'width_ratios': (1, 1.2)})
129
130ax = axs[0, 0]
131model.plt_field(ax, logK, title="$\\log K$", cmap="viridis", levels=17,
132                wells="color", finalize=False)
133
134tt = dt * np.arange(nSteps + 1)
135ax = axs[0, 1]
136for i, name in enumerate(producers):
137    ax.plot(tt, fw[:, i], label=name, c=f"C{i}")
138ax.plot(tt[k], J, "o", c="k", mfc="none", ms=10, zorder=3,
139        label=f"objective: {producers[well]} @ t={tt[k]:.2f}")
140ax.set(title="Water cut", xlabel="Time", ylabel="$f_w$", ylim=(-.02, 1))
141ax.legend(loc="upper left")
142
143ax = axs[1, 0]
144# The few cells next to the wells dominate; clip the color scale (the cmap's
145# `over`/`under` make `plt_field` extend the colorbar, rather than leave blanks).
146m = np.percentile(abs(G), 98)
147cmap = plt.get_cmap("RdBu_r")
148cmap = cmap.with_extremes(over=cmap(1.0), under=cmap(0.0))
149model.plt_field(ax, G, title="$∂J/∂\\log K$", cmap=cmap,
150                levels=np.linspace(-m, m, 21), cticks=[-m, 0, m],
151                wells="color", finalize=False)
152
153ax = axs[1, 1]
154for i, name in enumerate(producers):
155    ax.step(tt[:-1], G_bhp[1 + i], where="post", label=name, c=f"C{i}")
156ax.axvline(tt[k], c="k", ls=":", lw=1, label="objective's time")
157ax.axhline(0, c="k", lw=.5)
158ax.set(title="$∂J/∂p_\\mathrm{bh}$, per time step", xlabel="Time (of the control)",
159       ylabel="$∂J/∂p_\\mathrm{bh}$")
160ax.legend(loc="upper right")
161
162fig.tight_layout()
163
164# Regression values, checked by `tests/test_examples.py`.
165__digest__ = dict(water_cut    = fw,
166                  gradient     = G,
167                  gradient_bhp = G_bhp,
168                  directional  = [directional_logK, fd_logK, directional_bhp, fd_bhp])
169
170if __name__ == "__main__":
171    show()