examples.logo

 1"""A logo: the saturation field of a reservoir shaped like a smiley.
 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
 9The face is a disc, with the smile *cut out* of it (inactive). The eyes are
10the injectors -- one well, two completions -- whose water pools around them,
11and the nose the producer, drawing it down. Along the bottom lies an
12**aquifer** (ref `examples.aquifer`), whose water rises *around* the smile to
13meet the producer -- which is set to take more than is injected, the aquifer
14supplying the rest. Set `aquifer = False` to remove it: the rates are then
15balanced, so the picture can be compared with and without. The aquifer's
16contact is stroked in blue (`minires.plotting.Plot2D.plt_faces`).
17
18The figure is the plain `plt_field` (oil saturation: water in teal, oil in
19coral), stripped of axes, colorbar, and labels; the well markers remain (the
20eyes and nose are the markers) -- each coloured individually, by name
21(and by completion, for the two eyes), with the little centre dot turned off
22(`wells=dict(color=..., dot=False)`, ref `minires.plotting.Plot2D.well_scatter`). The stroke of the aquifer contact runs
23along the boundary *faces* of the contact cells, so it has no width but its `lw`,
24and the contours stop half a cell short of it; `cellwise=True` would instead paint
25the cells flat, up to those faces, with a pixelated outline.
26"""
27
28from typing import Any
29
30import matplotlib.pyplot as plt
31import numpy as np
32from minires import ResSim
33from minires.plotting import show
34from minires.wells import boundary_faces
35
36aquifer = True  # toggle: water beyond the bottom boundary, at pressure 1
37p_aq = 1.
38q_aq = .5 if aquifer else 0  # the deficit of injection, which the aquifer covers
39y_aq = .18  # the aquifer lies below this
40
41## The outline: a 32² unit square cut to a disc, less the smile
42model = ResSim(Lx=1, Ly=1, Nx=32, Ny=32)
43X, Y = model.mesh
44face = (X - .5)**2 + (Y - .5)**2 <= .46**2
45r, angle = np.hypot(X - .5, Y - .52), np.arctan2(Y - .52, X - .5)  # polar about the nose
46smile = (abs(r - .27) <= .03) & (-.75 * np.pi < angle) & (angle < -.25 * np.pi)
47model.active = face & ~smile
48
49## The wells (`Any`, so that the type checker sees `model.wells` as the `Wells` that
50## `__setattr__` makes of these records, rather than as the records)
51wells: Any = [dict(name="Nose", xy=[.5, .5], rate=-(1 + q_aq)),
52              dict(name="Eyes", xy=[[.25, .62], [.75, .62]], rate=+1)]  # 1 well, 2 compl.
53if aquifer:
54    # The contact: the boundary cells (those with a face to an inactive cell
55    # or off-grid) of the lower part of the outline.
56    xy = np.column_stack([X[model.active], Y[model.active]])
57    contact = xy[boundary_faces(model, xy).any(-1) & (xy[:, 1] < y_aq)]
58    wells = wells + [dict(name="Aq", xy=contact, aquifer=True, bhp=p_aq)]
59model.wells = wells
60
61S0 = np.zeros(model.Nxy)
62SS, PP = model.sim(.004, 55, S0, pbar=False)
63
64# The aquifer supplies exactly the deficit (incompressible), and the inactive cells stay inert.
65assert np.allclose(model.wells.rates_by_well.sum(0), 0)
66if aquifer:
67    assert np.allclose(model.wells.rates_by_well[-1], q_aq)
68assert (SS[:, ~model.active.ravel()] == 0).all()
69
70## Plot
71# A colour per marker: keyed by well name, and, the eyes being one well of two
72# completions, a colour per *completion* for them (left eye first, ref
73# `minires.wells.Wells.xy`).
74fig, ax = plt.subplots(num="Logo", clear=True, figsize=(5, 5))
75model.plt_field(ax, SS[-1], "oil", colorbar=False, labels=False, finalize=False, title="",
76                wells=dict(exclude=["Aq"], size=1.5, text=False, dot=False))
77if "Aq" in model.wells.names:
78    contact = model.wells.xy[model.wells.group == model.wells.nWell - 1]
79    model.plt_faces(ax, contact, lw=12)
80ax.axis("off")
81fig.tight_layout()
82
83# Regression values, checked by `tests/test_examples.py`.
84__digest__ = dict(smiley = SS[-1][model.active.ravel()])
85
86if __name__ == "__main__":
87    show()