examples.quarter_five_spot

  1"""Reproduce Fig. 6 of the reference paper, i.e. listing 9 -- and then vary it.
  2
  3This runs the familiar 5-spot well pattern on a homogeneous and isotropic permeability
  4which, thanks to symmetry, only requires computing one of the 4 quadrants,
  5giving it the quarter-five spot problem.
  6
  7There are some minor discrepancies compared with their Fig. 6.
  8
  9- They claim to plot the initial pressure, but it rather seems like the final one to me.
 10- They panels portray `t` values are not available from the chosen time steps.
 11- Their water front has a corner that is more protruding (not due to the previous issue).
 12
 13However, since we generate the very same output as the matlab code, we believe
 14the issue lies with the description in the paper, not with any error in the code.
 15
 16The same flood is then run twice more: with the implicit transport scheme (whose
 17output is likewise verified against matlab), and with **scheduled** rates. In the
 18latter, a second injector (NW corner) shares the load: for the first 10 steps all of
 19the water enters from the SW corner, thereafter the two injectors split it equally,
 20and the front is correspondingly lopsided. A rate may be specified as a *schedule*
 21(an array over time), as done there -- whereupon the wells held constant get
 22broadcast along with it -- or (for feedback control, e.g. shutting wells on water
 23breakthrough) by overriding `ResSim.well_controls`.
 24
 25In the figures:
 26
 27- "Fig. 6": the water, injected in the SW corner, advances on the producer in
 28  the NE, in the shape that gives the quarter five-spot its name -- essentially
 29  that of the isobars of the first panel.
 30- "final saturation": the implicit scheme is the more diffusive of the two: its
 31  transition band ($0.2 < S < 0.8$) is some 15% wider than the explicit one's.
 32  Under the schedule, the front is no longer symmetric about the diagonal: the
 33  NW injector, switched on late, has flooded a band along the north edge, pinching
 34  the SW injector's front, which fills the rest. The quadrant *between* the two
 35  injectors is thus swept hardest (mean saturation 0.80 in the NW, against 0.57 in
 36  the SE; 0.66 in each for the base case).
 37- "schedule": the schedule itself (left), and the oil saturation in the producer
 38  (right). Both floods break through only in the last two of the 28 steps, the
 39  scheduled one a step ahead: its water arrives along the north edge.
 40"""
 41
 42from mpl_tools.place import freshfig
 43import numpy as np
 44
 45from minires import ResSim
 46from minires.plotting import show
 47
 48## Setup
 49grid: dict = dict(Lx=1, Ly=1, Nx=64, Ny=64)
 50# Fluid properties are left at their defaults: vw = vo = 1, swc = sor = 0.
 51
 52# The base case: an injector (SW) and a producer (NE), at constant, opposite rates.
 53# A well is a record of its position and its (signed) rate; ref `ResSim.wells`.
 54model = ResSim(**grid, wells=[dict(name="SW", xy=[0, 0], rate=+1),
 55                              dict(name="NE", xy=[1, 1], rate=-1)])
 56
 57water_sat0 = model.fluid.swc * np.ones(model.Nxy)
 58nSteps = 28
 59dt = 0.7/nSteps
 60
 61# Scheduled: the SW injector carries everything until step 10, then they share.
 62rate_sw = .5*np.ones(nSteps)
 63rate_nw = .5*np.ones(nSteps)
 64rate_sw[:10] = 1
 65rate_nw[:10] = 0
 66model_sch = ResSim(**grid, wells=[
 67    dict(name="SW", xy=[0, 0], rate=+rate_sw),
 68    dict(name="NW", xy=[0, 1], rate=+rate_nw),
 69    dict(name="NE", xy=[1, 1], rate=-1),   # constant: broadcast to the schedule
 70])
 71
 72## Simulate
 73SS_exp, PP_exp = model.sim(dt, nSteps, water_sat0, pbar=False)
 74SS_imp, PP_imp = model.sim(dt, nSteps, water_sat0, implicit=True, pbar=False)
 75SS_sch, PP_sch = model_sch.sim(dt, nSteps, water_sat0, pbar=False)
 76
 77## Plot: the paper's Fig. 6 (from the explicit scheme)
 78kws = dict(levels=17, cmap="jet", origin=None, extent=(0, model.Lx, 0, model.Ly))
 79
 80fig, axs = freshfig("Fig. 6", nrows=2, ncols=3, sharex=True, sharey=True,
 81                    subplot_kw={'aspect': 'equal'})
 82
 83for ax, t in zip(axs.ravel(), [None, .14, .28, .42, .56, .70]):
 84    if ax.get_subplotspec().is_last_row() : ax.set_xlabel("x")  # noqa
 85    if ax.get_subplotspec().is_first_col(): ax.set_ylabel("y")  # noqa
 86
 87    if t is None:
 88        ax.set_title("Pressure")
 89        [P, V] = model.pressure_step(SS_exp[-1])  # Final+1 pressure
 90        ax.contourf(P.reshape(model.shape).T, **kws)
 91
 92    else:
 93        k = int(t/dt)
 94        ax.set_title("t = {:.2f}".format(k * dt))
 95        Z = SS_exp[k].reshape(model.shape).T  # transpose/flip for plot orientation
 96
 97        # Puts the values in gridcell centers (agrees w/ finite-vol. interpretation)
 98        # ax.imshow(Z[::-1], **kws)
 99
100        # Also colocates with gridcell centers, but does not extend to edges.
101        # ax.contourf(Z, levels=17, cmap="jet", origin="lower")
102
103        # Artificially stretches the field
104        ax.contourf(Z, **kws)
105
106fig.tight_layout()
107
108## Plot: scheme and schedule comparison (at the final time)
109fig, axs = freshfig("Quarter five-spot -- final saturation", ncols=3,
110                    sharex=True, sharey=True, subplot_kw={'aspect': 'equal'})
111
112for ax, (S, title) in zip(axs, [(SS_exp[-1], "Explicit (upwind)"),
113                                (SS_imp[-1], "Implicit (Newton)"),
114                                (SS_sch[-1], "Explicit, scheduled rates")]):
115    ax.set_title(title)
116    ax.set_xlabel("x")
117    ax.contourf(S.reshape(model.shape).T, **kws)
118axs[0].set_ylabel("y")
119# The implicit scheme is (as usual) the more diffusive: its front is smeared.
120
121fig.tight_layout()
122
123## Plot: the schedule itself, and the resulting production
124fig, axs = freshfig("Quarter five-spot -- schedule", ncols=2, figsize=(9, 3.5))
125
126tt = dt*(1 + np.arange(nSteps))
127for i, rate in enumerate(model_sch.wells.actual_rates[:2]):
128    x, y = model_sch.wells.xy[i]
129    name = model_sch.wells.names[i]
130    axs[0].step(tt, rate, where="post", label=f"Injector {name} @ ({x:.2f}, {y:.2f})")
131axs[0].set(title="Injection rates", xlabel="Time", ylabel="Rate", ylim=(-.05, 1.05))
132axs[0].legend()
133
134prd = model.xy2ind(*model.wells.xy[1])  # the NE producer's cell (in both models)
135model.plt_production(axs[1], np.column_stack([SS_exp[1:, prd], SS_sch[1:, prd]]),
136                     finalize=False, labels=["NE (base case)", "NE (scheduled)"])
137fig.tight_layout()
138
139## Animation
140animation = model.anim(SS_exp, SS_exp[1:, [prd]])
141
142# Regression values, checked by `tests/test_examples.py`.
143# The sub-sampling `[::600]` matches that of the matlab reference values.
144__digest__ = dict(explicit  = SS_exp[-1, ::600],
145                  implicit  = SS_imp[-1, ::600],
146                  scheduled = SS_sch[-1, ::600])
147
148if __name__ == "__main__":
149    show()