examples.inactive_cells

 1"""An irregular reservoir on the rectangular grid: an outline, and a sealing fault.
 2
 3The grid stays rectangular; the reservoir need not. `active`
 4(ref `minires.ResSim.active`) marks the cells that take part; the rest are
 5inert -- closed off (zero flux across their faces), without equations, and
 6their state carried through unchanged -- as if they were not there.
 7
 8- **An outline**: here an ellipse, cut out of a 48² grid. The flow goes around
 9  the (now curved) boundary as it would around the box.
10- **A sealing fault**: a line of inactive cells. It need only be unbroken --
11  a diagonal staircase suffices, fluxes passing through faces alone -- to block
12  the flow entirely, which must then go around its tip: the front reaches the
13  producer (top right) by the long way round, beneath the fault, and the region
14  behind the fault is swept last.
15  Should it cut the reservoir in two, each half is a reservoir of its own, and
16  (being incompressible) must balance its own rates -- ref `active`.
17
18In the figure: the permeability, painted cell by cell (`cellwise=True`), which
19shows the mask exactly as it is; then the oil saturation at two times, as
20contours -- these interpolate between cell centres, so they also leave blank
21the half-cell margin around the inactive cells (as around the domain).
22"""
23
24from mpl_tools.place import freshfig
25import numpy as np
26from scipy.ndimage import uniform_filter as smooth
27
28from minires import ResSim
29from minires.plotting import show
30
31rng = np.random.default_rng(3)  # Reproducibility (the values are regression tested)
32
33## An irregular reservoir: an ellipse, with a sealing fault
34model = ResSim(Lx=1, Ly=1, Nx=48, Ny=48,
35               wells=[dict(xy=[.2, .5], rate=+1, name="Inj"),
36                      dict(xy=[.8, .65], rate=-1, name="Prd")])
37X, Y = model.mesh
38outline = ((X - .5) / .47)**2 + ((Y - .5) / .4)**2 <= 1
39# The fault: one cell per row (slope < 1 ⇒ an unbroken staircase), from y = .3 up
40x_fault = .55 + .25 * (Y - .5)
41fault = (abs(X - x_fault) <= model.hx / 2 + 1e-9) & (Y > .3)
42model.active = outline & ~fault
43logK = 4 * smooth(smooth(rng.standard_normal(model.shape)))
44model.K = np.exp(logK)
45
46dt = .01
47nSteps = 60
48S0 = np.zeros(model.Nxy)
49SS, PP = model.sim(dt, nSteps, S0, pbar=False)
50
51# The inactive cells are inert: their state is simply carried through.
52inactive = ~model.active.ravel()
53assert (SS[:, inactive] == 0).all() and (PP[:, inactive] == 0).all()
54# The active ones conserve the water: what is injected and not yet produced is in place.
55iprd = model.xy2ind(*model.wells.xy[1])
56water_cut = model.fluid.fractional_flow(SS[:, iprd])
57kBT = int(np.argmax(water_cut > 0))  # breakthrough: during step kBT-1 → kBT
58assert 0 < kBT < nSteps
59pv = model.pore_volume()[~inactive]
60assert np.isclose((pv * SS[kBT - 1][~inactive]).sum(), (kBT - 1) * dt)
61
62## Plot
63fig, axs = freshfig("Inactive cells", ncols=3, figsize=(11, 3.8),
64                    sharex=True, sharey=True)
65kws: dict = dict(finalize=False, wells=dict(size=.5))
66model.plt_field(axs[0], logK, cellwise=True, cmap="viridis", levels=17,
67                title="log-Permeability (cells painted flat)", **kws)
68for i, k in enumerate([nSteps // 3, nSteps]):
69    model.plt_field(axs[1 + i], SS[k], "oil", colorbar=(i == 1), labels=False,
70                    title=f"Oil saturation, t = {k*dt:.2f}", **kws)
71fig.tight_layout()
72
73# Regression values, checked by `tests/test_examples.py`.
74__digest__ = dict(sat_final = SS[-1][~inactive],
75                  p_final   = PP[-1][~inactive],
76                  water_cut = water_cut)
77
78if __name__ == "__main__":
79    show()