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 -- smoothed, log-normal -- permeability
  6dictates). 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 opposite
 35  corner's (SW) is positive throughout: drawing more to it starves the NE
 36  well. The neighbours' (NW, SE) are smaller and change sign over time:
 37  drawing more flow towards a neighbour reshapes the stream tube to NE as much
 38  as it starves it, and which effect wins depends on the heterogeneity and on
 39  where the front is. All of them vanish from the objective's time step on --
 40  a control cannot affect what came before it. Their sum over time is the
 41  gradient with respect to a constant-in-time BHP.
 42
 43.. note:: The gradient is with respect to the *isotropic* $\\log K$.
 44
 45    `model.K` holds both components, so the adjoint returns a gradient for
 46    each; as they are set equal here, the gradient of the single field is
 47    their sum, ref `tlm.Gradient.logK`.
 48
 49.. note:: The well indices are held fixed, as in the adjoint.
 50
 51    They are computed (from the permeability at the wells) once, when the
 52    wells are configured, and stored (`Wells.WI`); the finite differences
 53    below reuse them, as `tlm` assumes (ref its "What is differentiated").
 54"""
 55
 56from dataclasses import replace
 57
 58from mpl_tools.place import freshfig
 59import matplotlib.pyplot as plt
 60import numpy as np
 61from scipy.ndimage import uniform_filter as smooth
 62
 63from minires import ResSim
 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)
 71logK = 3 * smooth(smooth(rng.standard_normal((grid["Nx"], grid["Ny"]))))
 72model = ResSim(**grid, K=np.exp(logK), wells=[  # isotropic: K broadcast to both components
 73    dict(xy=[.5, .5], rate=+1, name="inj"),
 74    dict(xy=[1 , 1 ], bhp=0, rw=1e-3, name="NE"),
 75    dict(xy=[0 , 1 ], bhp=0, rw=1e-3, name="NW"),
 76    dict(xy=[0 , 0 ], bhp=0, rw=1e-3, name="SW"),
 77    dict(xy=[1 , 0 ], bhp=0, rw=1e-3, name="SE"),
 78])
 79
 80dt, nSteps = .05, 30
 81S0 = np.zeros(model.Nxy)
 82SS, PP = model.sim(dt, nSteps, S0, pbar=False)
 83
 84## The objective: water cut at the NE producer at time `k`
 85producers = model.wells.names[1:]
 86prd = model.xy2ind(*model.wells.xy[1:].T)  # their cells
 87
 88
 89def water_cut(model, SS):
 90    """`(nSteps+1, nPrd)` water cut at each producer, for each stored time."""
 91    return np.array([model.fluid.fractional_flow(S)[prd] for S in SS])
 92
 93
 94fw = water_cut(model, SS)
 95well, k = 0, 12  # NE, just after breakthrough (fw ≈ .5)
 96J = fw[k, well]
 97
 98## Its gradient, by the adjoint
 99dJ_dSS = np.zeros_like(SS)
100dJ_dSS[k, prd[well]] = model.fluid.dfractional_flow(SS[k])[prd[well]]  # f_w'(s)
101grad = adjoint(model, dt, SS, PP, dJ_dSS)
102G = grad.logK.sum(0)  # isotropic ⇒ sum the components
103G_bhp = grad.bhp      # (nComp, nSteps); zero for the (rate-controlled) injector
104
105
106## Check: finite differences in random directions of log K and of the BHP schedule
107def J_of(logK, bhp):
108    """The objective of a model with `logK`, its producers at the schedule `bhp`."""
109    m = ResSim(**grid, K=np.exp(logK), wells=replace(model.wells, bhp=bhp))  # same `WI`
110    return water_cut(m, m.sim(dt, nSteps, S0, pbar=False)[0])[k, well]
111
112
113eps = 1e-5
114bhp = model.wells.bhp  # (nComp, 1): constant in time
115d_logK = rng.standard_normal(model.shape)
116d_bhp = np.where(np.isfinite(bhp), rng.standard_normal((model.nComp, nSteps)), 0)
117fd_logK = (J_of(logK + eps*d_logK, bhp) - J_of(logK - eps*d_logK, bhp)) / (2*eps)
118fd_bhp = (J_of(logK, bhp + eps*d_bhp) - J_of(logK, bhp - eps*d_bhp)) / (2*eps)
119directional_logK = (G * d_logK).sum()
120directional_bhp = (G_bhp * d_bhp).sum()
121assert abs(fd_logK - directional_logK) < 1e-4 * abs(directional_logK), (fd_logK, directional_logK)
122assert abs(fd_bhp - directional_bhp) < 1e-4 * abs(directional_bhp), (fd_bhp, directional_bhp)
123
124## Plot
125fig, axs = freshfig("Water-cut gradient", ncols=2, nrows=2, figsize=(10, 8),
126                    gridspec_kw={'width_ratios': (1, 1.2)})
127
128ax = axs[0, 0]
129model.plt_field(ax, logK, title="$\\log K$", cmap="viridis", levels=17,
130                wells="color", finalize=False)
131
132tt = dt * np.arange(nSteps + 1)
133ax = axs[0, 1]
134for i, name in enumerate(producers):
135    ax.plot(tt, fw[:, i], label=name, c=f"C{i}")
136ax.plot(tt[k], J, "o", c="k", mfc="none", ms=10, zorder=3,
137        label=f"objective: {producers[well]} @ t={tt[k]:.2f}")
138ax.set(title="Water cut", xlabel="Time", ylabel="$f_w$", ylim=(-.02, 1))
139ax.legend(loc="upper left")
140
141ax = axs[1, 0]
142# The few cells next to the wells dominate; clip the color scale (the cmap's
143# `over`/`under` make `plt_field` extend the colorbar, rather than leave blanks).
144m = np.percentile(abs(G), 98)
145cmap = plt.get_cmap("RdBu_r")
146cmap = cmap.with_extremes(over=cmap(1.0), under=cmap(0.0))
147model.plt_field(ax, G, title="$∂J/∂\\log K$", cmap=cmap,
148                levels=np.linspace(-m, m, 21), cticks=[-m, 0, m],
149                wells="color", finalize=False)
150
151ax = axs[1, 1]
152for i, name in enumerate(producers):
153    ax.step(tt[:-1], G_bhp[1 + i], where="post", label=name, c=f"C{i}")
154ax.axvline(tt[k], c="k", ls=":", lw=1, label="objective's time")
155ax.axhline(0, c="k", lw=.5)
156ax.set(title="$∂J/∂p_\\mathrm{bh}$, per time step", xlabel="Time (of the control)",
157       ylabel="$∂J/∂p_\\mathrm{bh}$")
158ax.legend(loc="upper right")
159
160fig.tight_layout()
161
162# Regression values, checked by `tests/test_examples.py`.
163__digest__ = dict(water_cut    = fw,
164                  gradient     = G,
165                  gradient_bhp = G_bhp,
166                  directional  = [directional_logK, fd_logK, directional_bhp, fd_bhp])
167
168if __name__ == "__main__":
169    show()