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
24import matplotlib.pyplot as plt
25import numpy as np
26
27from minires import ResSim
28from minires.geostat import gaussian_fields
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 = gaussian_fields(model.mesh, r=.15, rng=rng)[0].reshape(model.shape)
44# logK = 4 * smooth(smooth(rng.standard_normal(model.shape)))  # the old way
45model.K = np.exp(logK)
46
47dt = .01
48nSteps = 60
49S0 = np.zeros(model.Nxy)
50SS, PP = model.sim(dt, nSteps, S0, pbar=False)
51
52# The inactive cells are inert: their state is simply carried through.
53inactive = ~model.active.ravel()
54assert (SS[:, inactive] == 0).all() and (PP[:, inactive] == 0).all()
55# The active ones conserve the water: what is injected and not yet produced is in place.
56iprd = model.xy2ind(*model.wells.xy[1])
57water_cut = model.fluid.fractional_flow(SS[:, iprd])
58kBT = int(np.argmax(water_cut > 0))  # breakthrough: during step kBT-1 → kBT
59assert 0 < kBT < nSteps
60pv = model.pore_volume()[~inactive]
61assert np.isclose((pv * SS[kBT - 1][~inactive]).sum(), (kBT - 1) * dt)
62
63## Plot
64fig, axs = plt.subplots(num="Inactive cells", clear=True, ncols=3, figsize=(11, 3.8),
65                        sharex=True, sharey=True)
66kws: dict = dict(finalize=False, wells=dict(size=.5))
67model.plt_field(axs[0], logK, cellwise=True, cmap="viridis", levels=17,
68                title="log-Permeability (cells painted flat)", **kws)
69for i, k in enumerate([nSteps // 3, nSteps]):
70    model.plt_field(axs[1 + i], SS[k], "oil", colorbar=(i == 1), labels=False,
71                    title=f"Oil saturation, t = {k*dt:.2f}", **kws)
72fig.tight_layout()
73
74# Regression values, checked by `tests/test_examples.py`.
75__digest__ = dict(sat_final = SS[-1][~inactive],
76                  p_final   = PP[-1][~inactive],
77                  water_cut = water_cut)
78
79if __name__ == "__main__":
80    show()