examples.logo

  1"""A logo: the saturation field of a reservoir shaped like a smiley, or a yin-yang.
  2
  3No physics is illustrated here that the other examples do not; the goal is a
  4picture. But it does show how little a shape costs on the rectangular grid:
  5an outline is a boolean expression in the mesh coordinates, assigned to
  6`minires.ResSim.active` (ref `examples.inactive_cells`), holes included,
  7and the wells and the aquifer are records like any other's.
  8
  9- **The smiley**: a disc, the smile *cut out* of it (inactive). The eyes are
 10  the injectors -- one well, two completions -- whose water pools around them,
 11  and the nose the producer, drawing it down. Along the bottom lies an
 12  **aquifer** (ref `examples.aquifer`), whose water rises *around* the smile to
 13  meet the producer -- which is set to take more than is injected, the aquifer
 14  supplying the rest. Set `aquifer = False` to remove it: the rates are then
 15  balanced, so the picture can be compared with and without. The aquifer's
 16  contact is stroked in blue (`minires.plotting.Plot2D.plt_faces`).
 17- **The yin-yang**: a disc, the S-curve through it -- two semicircles of half
 18  its radius -- a *barrier* of inactive cells, like the fault of
 19  `examples.inactive_cells`, except that it stops short of the rim at the
 20  bottom. The two dots are the wells: the injector in the upper one, the
 21  producer in the lower. So the injector's half floods, and the water can reach
 22  the other half only by way of the gap at the bottom -- where the front stands
 23  at the time shown. The water is made ten times as viscous as the oil, so that
 24  the displacement is piston-like (the Buckley-Leverett shock is then at nearly
 25  full water saturation, ref `examples.buckley_leverett`): the flooded half is
 26  uniformly teal, and the front is sharp. With equal viscosities the front is a
 27  rarefaction -- a spread of pale contours -- and the two halves contrast poorly.
 28
 29The figures are the plain `plt_field` (oil saturation: water in teal, oil in
 30coral), stripped of axes, colorbar, and labels; the well markers remain (the
 31eyes and nose, the dots, are the markers). The stroke of the aquifer contact runs
 32along the boundary *faces* of the contact cells, so it has no width but its `lw`,
 33and the contours stop half a cell short of it; `cellwise=True` would instead paint
 34the cells flat, up to those faces, with a pixelated outline.
 35"""
 36
 37from mpl_tools.place import freshfig
 38import numpy as np
 39
 40from minires import ResSim
 41from minires.plotting import show
 42from minires.wells import boundary_faces
 43
 44aquifer = True  # toggle: water beyond the bottom boundary, at pressure 1
 45p_aq = 1.
 46
 47def make(outline, wells, y_aq=None, **kws):
 48    """A 32² unit square cut to `outline(X, Y)`, with `wells`, and the aquifer below `y_aq`."""
 49    model = ResSim(Lx=1, Ly=1, Nx=32, Ny=32, **kws)
 50    X, Y = model.mesh
 51    model.active = outline(X, Y)
 52    if aquifer and y_aq is not None:
 53        # The contact: the boundary cells (those with a face to an inactive cell
 54        # or off-grid) of the lower part of the outline.
 55        xy = np.column_stack([X[model.active], Y[model.active]])
 56        contact = xy[boundary_faces(model, xy).any(-1) & (xy[:, 1] < y_aq)]
 57        wells = wells + [dict(name="Aq", xy=contact, aquifer=True, bhp=p_aq)]
 58    model.wells = wells
 59    return model
 60
 61q_aq = .5 if aquifer else 0  # the deficit of injection, which the aquifer covers
 62
 63## The smiley
 64def smiley(X, Y):
 65    face = (X - .5)**2 + (Y - .5)**2 <= .46**2
 66    r, angle = np.hypot(X - .5, Y - .52), np.arctan2(Y - .52, X - .5)  # polar about the nose
 67    smile = (abs(r - .27) <= .03) & (-.75 * np.pi < angle) & (angle < -.25 * np.pi)
 68    return face & ~smile
 69
 70model = make(smiley, y_aq=.18, wells=[
 71    dict(name="Nose", xy=[.5, .5], rate=-(1 + q_aq)),
 72    dict(name="Eyes", xy=[[.25, .62], [.75, .62]], rate=+1),  # one well, 2 completions
 73])
 74S0 = np.zeros(model.Nxy)
 75SS, PP = model.sim(.004, 55, S0, pbar=False)
 76
 77# The aquifer supplies exactly the deficit (incompressible), and the inactive cells stay inert.
 78assert np.allclose(model.wells.rates_by_well.sum(0), 0)
 79if aquifer:
 80    assert np.allclose(model.wells.rates_by_well[-1], q_aq)
 81assert (SS[:, ~model.active.ravel()] == 0).all()
 82
 83## The yin-yang
 84R, w, gap = .46, .03, .07  # radius; the barrier's half-width; the opening at the bottom
 85
 86def yinyang(X, Y):
 87    disc = np.hypot(X - .5, Y - .5) <= R
 88    upper = (abs(np.hypot(X - .5, Y - .5 - R/2) - R/2) <= w) & (X >= .5) & (Y >= .5)
 89    lower = (abs(np.hypot(X - .5, Y - .5 + R/2) - R/2) <= w) & (X <= .5) & (Y <= .5)
 90    return disc & ~((upper | lower) & (Y > .5 - R + gap))
 91
 92yy_model = make(yinyang, fluid=dict(vw=10), wells=[  # viscous water ⇒ a piston-like front
 93    dict(name="Inj", xy=[.5, .75], rate=+1),
 94    dict(name="Prd", xy=[.5, .25], rate=-1),
 95])
 96SS_yy, PP_yy = yy_model.sim(.005, 62, S0, pbar=False)
 97
 98# The barrier holds: the water enters the producer's half from the bottom, so the
 99# upper part of that half (right of the upper semicircle) is still dry.
100X, Y = yy_model.mesh
101far = yy_model.active & (X > .5) & (Y > .5) & (np.hypot(X - .5, Y - .5 - R/2) > R/2)
102assert SS_yy[-1].reshape(yy_model.shape)[far].max() < 1e-6
103
104## Plot
105for name, m, S in [("smiley", model, SS[-1]), ("yin-yang", yy_model, SS_yy[-1])]:
106    fig, ax = freshfig(f"Logo ({name})", figsize=(5, 5))
107    m.plt_field(ax, S, "oil", colorbar=False, labels=False, finalize=False,
108                title="", wells=dict(exclude=["Aq"], size=1.8, text=False))
109    if "Aq" in m.wells.names:
110        m.plt_faces(ax, m.wells.xy[m.wells.group == m.wells.nWell - 1], color="darkblue", lw=8)
111    ax.axis("off")
112    fig.tight_layout()
113
114# Regression values, checked by `tests/test_examples.py`.
115__digest__ = dict(smiley  = SS[-1][model.active.ravel()],
116                  yinyang = SS_yy[-1][yy_model.active.ravel()])
117
118if __name__ == "__main__":
119    show()