minires.wells

  1"""The wells: their configuration, their bookkeeping, and the well model.
  2
  3Everything here depends on the wells alone, on the *geometry* they sit in
  4(`minires.grid.Grid2D`), and -- for the well index -- on the permeability
  5`K`. What couples the wells to the *fluids* (the mobilities of
  6`minires.fluids.Fluid.RelPerm`) or to the linear system (the source field, the
  7BHP contributions) stays with the simulator, which reads the arrays assembled
  8here.
  9
 10The unit of account is the **completion**, not the well: a well may span
 11several cells (ref `well_path`), and it is the completions that the model
 12solves for. Hence `Wells.nComp` counts the rows of every array here, whereas
 13`Wells.nWell` counts the wells that `Wells.group` says they compose --
 14a distinction that serves the *reporting* alone (`Wells.rates_by_well`, the
 15plot labels); the physics never groups.
 16
 17
 18.. note:: Spreading a well over its neighbouring cells was implemented, then omitted.
 19
 20    Snapped to a cell centre, a well has its every effect a staircase function
 21    of its coordinates -- constant within a cell -- which leaves an
 22    optimisation over well *positions* with no gradient to work with until the
 23    well crosses into the next cell (as in the EnOpt of
 24    [HistoryMatching](https://github.com/patnr/HistoryMatching)). Distributing
 25    it over the 4 surrounding cells by bilinear weights -- a *mollified* rather
 26    than a rounded delta -- fixes that, and is not merely cosmetic: the
 27    position dependence so obtained tracks that of a 3-times-refined grid to
 28    within that grid's own discretization spread, i.e. some 5 times closer than
 29    rounding manages. It does require the well index to be corrected for the
 30    cells dividing the load, Peaceman's equivalent radius being derived for the
 31    whole source in one cell: a well divided 4 ways under-reports its drawdown
 32    by 23%, non-convergently, unless the weighted geometric mean of the
 33    intercell distances is substituted for $r_e$ (which recovers 0.3%).
 34
 35    It was nonetheless judged not to earn its complexity -- a stencil threaded
 36    through the well assembly, the well model and the plotting, for a
 37    convenience of the optimiser rather than a fidelity of the simulator. The
 38    work is preserved on the branch `well-spread` (`ResSim.spread_wells`,
 39    `Grid2D.xy2stencil`, `wells._share_WI`, and `tests/test_spread.py`, which
 40    pins each of the measurements quoted above).
 41
 42## Theory
 43
 44A well is either an **injector** or a **producer** -- the terminal states of the
 45sources and sinks, $q$, of the governing equations.
 46The two are not distinct objects here: a well is an injector or a producer
 47merely by the *sign* of its rate (ref `minires.wells.Wells.rates`), and
 48under BHP control not even by that, the direction being left to the pressures
 49(ref `minires.wells.Wells.bhp`).
 50Its **completion** is the equipment that connects the **wellbore** to the rock,
 51whose interface is the **sandface**; it may be *open hole*, or cased and
 52**perforated**. A well may have several completions, e.g. one per layer,
 53or (here) one per grid cell traversed by the well path -- which the model
 54assembles individually, grouping them back into wells only for the reporting
 55(ref `minires.ResSim.wells`).
 56
 57Because a wellbore (radius $r_w \\sim 0.1$ m) is orders of magnitude smaller than a
 58grid block, its pressure is not resolved by the grid: the radial solution
 59$p \\sim \\ln r$ spends most of its variation within the well's own cell.
 60The **bottom-hole pressure** (BHP), $p_\\mathrm{bh}$, is the pressure at the
 61sandface, which is thus a *sub-grid* quantity;
 62measured at the surface instead, it is the *tubing head pressure* (THP), the
 63difference being hydrostatic and friction losses in the tubing -- neither of which
 64a 2D areal model has, which is why BHP is where this model stops.
 65The **drawdown** is the pressure difference driving the flow,
 66$p_\\mathrm{cell} - p_\\mathrm{bh}$: positive for a producer (whence its name),
 67while for an injector it is negative, being an *overpressure*.
 68
 69The **productivity index** (PI) is the resulting constant of proportionality,
 70$q = \\mathrm{PI} \\cdot \\Delta p$ (the *injectivity index* for injectors),
 71familiar from well testing. Its counterpart in a simulator is the
 72**well index** (WI), which is the same relation with the fluid factored out,
 73$q = \\mathrm{WI} \\, \\lambda_t \\, \\Delta p$,
 74so that $\\mathrm{WI}$ depends only on geometry and rock,
 75and the mobility $\\lambda_t$ carries the (time-varying) fluid dependence.
 76It is the well's analogue of the transmissibility $t_{ij}$ of eqn. (10).
 77Two things enter it, beyond $r_w$ and $\\mathbf{K}$:
 78
 79- The **skin**, $S$, is a dimensionless lumping of all *near-wellbore* effects that
 80  the model does not resolve. It is *positive* for damage (drilling mud invasion,
 81  fines migration, scale) and *negative* for stimulation (acidizing, or hydraulic
 82  **fracturing**), and enters additively to a logarithm, so a skin of $5$ is a lot.
 83- The **equivalent radius**, $r_e$, is the purely *numerical* ingredient: the radius
 84  at which the analytic radial pressure equals the numerical *cell* pressure,
 85  $r_e \\approx 0.2 h$ for the 5-point stencil of TPFA. Beware that the same symbol
 86  is used, in well testing, for the (physical) *drainage radius*.
 87
 88Ref `minires.wells.peaceman_WI` for the formula combining these.
 89
 90A well is **controlled** either by prescribing its rate, or its BHP,
 91the other then being an outcome (ref `minires.wells.Wells.bhp`).
 92Reality is closer to the latter -- one sets a pump speed or a **choke** opening,
 93and the reservoir decides the rate -- but the *rate* is what is usually planned for.
 94Field practice is therefore rate control subject to a BHP *constraint*
 95(from the fracture pressure of an injector, or the lift capacity, or the bubble
 96point, of a producer), **switching** control mode whenever the constraint binds.
 97A well producing at a rate too low to be worthwhile is **shut in**;
 98if it cannot flow unaided it needs **artificial lift** (gas lift, or a downhole pump).
 99
100The **water cut** of a producer is the water fraction of what it produces,
101i.e. the $f(s)$ of its cell; **breakthrough** is when the injected water first
102arrives, after which the water cut climbs and the well eventually becomes uneconomic.
103How much of the oil the flood has contacted by then is the **sweep efficiency**,
104which is governed by the **mobility ratio**, $M = \\lambda_w / \\lambda_o$
105evaluated behind and ahead of the front. $M > 1$ is *unfavourable*:
106the (less viscous) water outruns the oil in **viscous fingering**,
107and even more so along the high-permeability *channels* -- the motivation for
108polymer injection, which fixes $M$ by thickening the water.
109The corresponding vertical phenomenon (which a 2D areal model cannot see)
110is **coning**, of water up, or gas down, into the completion.
111
112Wells are drilled in repeated **patterns**, of which the *five-spot* (a producer
113at the centre of four injectors, or vice-versa if *inverted*) is the classic;
114by symmetry it suffices to simulate the *quarter five-spot*, as in the examples here.
115They need not be vertical: *deviated*, *horizontal* and *multilateral* wells
116contact more rock per well, at the price of an **allocation** problem,
117namely how the total rate distributes itself among the completions
118(ref `minires.wells.well_path`).
119Later interventions to restore or improve a well are **workovers**,
120and drilling extra wells between the existing ones is **infill drilling**.
121
122"""
123
124from dataclasses import dataclass
125from typing import TYPE_CHECKING, Any
126
127import numpy as np
128
129from minires._repr import AlignedRepr
130
131if TYPE_CHECKING:
132    from minires import ResSim
133
134
135def peaceman_WI(model: "ResSim", xy: Any, rw: float, skin: float = 0.0) -> np.ndarray:
136    """Peaceman's well index for wells at `xy`, of radius `rw`, in `model`.
137
138    Applied for you to a well of `Wells.from_records` given an `rw`,
139    which is the convenient way to use it.
140
141    $$ WI = \\frac{2 π \\sqrt{k_x k_y}}{\\ln(r_e / r_w) + \\mathrm{skin}} $$
142
143    where the *equivalent radius*, $ r_e $, is the distance from the well at
144    which the (analytic, radial) pressure equals the (numerical) pressure of
145    the well's cell:
146    $$ r_e = 0.28 \\,
147       \\frac{\\sqrt{\\sqrt{k_y/k_x} \\, h_x^2
148                      + \\sqrt{k_x/k_y} \\, h_y^2}}
149              {(k_y/k_x)^{1/4} + (k_x/k_y)^{1/4}} \\,, $$
150    which reduces to the familiar $ r_e = 0.198 \\, h $ on an isotropic,
151    square grid. That constant is not a fudge factor: it is a property of
152    the 5-point stencil that `minires.ResSim.TPFA` assembles, and this model
153    reproduces it (`tests/test_wells.py` recovers $ r_e / h → 0.198 $ from the
154    simulated drawdown, and thereby the analytic, radial well pressure to
155    within 0.2%, on grids from 16² to 64²).
156
157    >>> from minires import ResSim
158    >>> model = ResSim(Lx=1, Ly=1, Nx=32, Ny=32)
159    >>> peaceman_WI(model, [[.5, .5]], rw=1e-3).round(4)
160    array([3.4476])
161
162    .. note:: `model` is used for its grid and its `K` alone.
163
164        Which is why this is a free function, not a method: the well index is a
165        property of a *location*, not of the well configuration, and it is
166        evaluated once, when asked for -- so a later edit of `K` does not
167        retroactively change a `Wells.WI` computed from it.
168
169    .. note:: $ WI \\, λ_t \\, Δp $ comes out as a rate *per unit thickness*.
170
171        Ref. `minires.ResSim.cdarcy`.
172
173    .. note:: `rw` must be given in the same length unit as `Lx`.
174    """
175    xy = np.asarray(xy, float).reshape((-1, 2))
176    ix, iy = model.xy2sub(*xy.T)
177    kx, ky = model.K[0][ix, iy], model.K[1][ix, iy]
178    # fmt: off
179    a, b = np.sqrt(ky/kx), np.sqrt(kx/ky)
180    r_e  = .28 * np.sqrt(a*model.hx**2 + b*model.hy**2) / (a**.5 + b**.5)
181    return model.cdarcy * 2*np.pi*np.sqrt(kx*ky) / (np.log(r_e/rw) + skin)
182    # fmt: on
183
184
185def well_path(model: "ResSim", vertices: Any, rw: float, skin: float = 0.0) -> tuple:
186    """Discretize a well *path* (a polyline): 1 weighted completion per cell.
187
188    Applied for you to a well of `Wells.from_records` given a `path`, which is
189    the convenient way to use it: the three returned arrays then need not be
190    assembled (with those of the other wells) by hand.
191
192    Returns `(xy, WI, alloc)`:
193
194    - `xy`: centres of the cells that the path traverses -- i.e. a value for
195      `Wells.xy`. Several completions act as a single well simply by
196      being several wells: `minires.ResSim.assemble_wells` superimposes them.
197    - `WI`: their well indices, i.e. a value for `Wells.WI`. Each is
198      `peaceman_WI` for its cell, scaled by the fraction of
199      that cell which the path actually traverses (so a cell merely clipped
200      by the path contributes proportionally less).
201    - `alloc`: `WI / WI.sum()`, for apportioning the rate among its completions:
202      `rates = rate * alloc[:, None]` (the rate signed as usual).
203      This is the standard (static) allocation -- proportional to the well index,
204      hence to both the contacted length and the local permeability.
205
206    >>> from minires import ResSim
207    >>> model = ResSim(Lx=1, Ly=1, Nx=10, Ny=10)
208    >>> xy, WI, alloc = well_path(model, [[.05, .05], [.45, .05]], rw=1e-2)
209    >>> xy.T[0]  # the traversed cells, in x
210    array([0.05, 0.15, 0.25, 0.35, 0.45])
211    >>> alloc  # the end cells, entered mid-way, get half the rate
212    array([0.125, 0.25 , 0.25 , 0.25 , 0.125])
213
214    .. note:: `alloc` is exact under BHP control, approximate under rate control.
215
216        Under BHP control (assuming 0 gravity and friction) the completions
217        simply share a `p_bh`. Under *rate* control, the allocation holds only
218        if the cell pressures are equal. Solving for it would make
219        $ p_\\mathrm{bh} $ an extra
220        unknown, i.e. a bordered linear system -- which the 5-diagonal
221        assembly of `minires.ResSim.TPFA` (Listing 1) is not set up for. Use
222        `minires.ResSim.well_controls` to reallocate per step, if it matters.
223
224    .. warning:: The completions are treated as independent *vertical* wells.
225
226        Which seems reasonable in a 2D areal model. Thus they count towards
227        `Wells.nComp`, not `Wells.nWell` -- ref `Wells.group`, which
228        `Wells.from_records` sets for you.
229    """
230    V = np.asarray(vertices, float).reshape((-1, 2))
231    assert len(V) >= 2, "A well path needs at least 2 vertices."
232    # Walk the polyline, accumulating traversed length per cell
233    lengths: dict = {}
234    for p0, p1 in zip(V[:-1], V[1:]):
235        d = p1 - p0
236        L = float(np.hypot(*d))
237        if L == 0:
238            continue
239        ts = model._crossings(p0, d)
240        mids = p0 + np.outer((ts[:-1] + ts[1:]) / 2, d)
241        for mid, dt in zip(mids, np.diff(ts)):
242            sub = tuple(int(i) for i in model.xy2sub(*mid))
243            lengths[sub] = lengths.get(sub, 0.0) + L * dt
244    # Discard the slivers left by corner crossings
245    total = sum(lengths.values())
246    lengths = {k: v for k, v in lengths.items() if v > 1e-9 * total}
247
248    subs = np.array(list(lengths))
249    xy = model.sub2xy(*subs.T).T
250    # Scale each WI by how much of its cell the path traverses, relative
251    # to the cell size -- so an axis-aligned full crossing scores exactly 1
252    # (and a diagonal one √2, it contacting that much more rock).
253    frac = np.array(list(lengths.values())) / np.sqrt(model.h2)
254    WI = frac * peaceman_WI(model, xy, rw, skin)
255    return xy, WI, WI / WI.sum()
256
257
258def aquifer_WI(model: "ResSim", xy: Any, faces: str = "WESN") -> np.ndarray:
259    """The "well index" of an aquifer contact: the transmissibility of the
260    boundary face(s) of each cell at `xy` -- from its centre out to the face.
261
262    An aquifer -- water-bearing rock beyond the reservoir's boundary, at a
263    pressure $ p_\\mathrm{aq} $ of its own -- feeds each reservoir cell that
264    touches it at the rate $ T \\, (p_\\mathrm{aq} - p) $, which is the law of a
265    BHP-controlled well, $ WI \\, λ_t \\, (p_\\mathrm{bh} - p) $. So an aquifer *is*
266    a BHP-controlled well, completed in every cell it touches, with
267    `bhp = p_aq` and this for `WI` -- and needs nothing else of the model:
268    the influx enters `minires.ResSim.assemble_wells` like any well's,
269    anchors the pressure (so that, if incompressible, a lone producer is fine:
270    the aquifer supplies it), is reported in `Wells.actual_rates`, and is
271    handled by the adjoint, `minires.tlm`, as any BHP well is. `Wells.from_records` applies it to a
272    well given `aquifer=True` (or `aquifer=faces`):
273
274    >>> from minires import ResSim
275    >>> model = ResSim(Lx=1, Ly=1, Nx=4, Ny=4, wells=[
276    ...     dict(name="Aq", xy=[[0, .3], [0, .5], [0, .7]], aquifer=True, bhp=2),
277    ...     dict(name="P1", xy=[1, 1], rate=-1),
278    ... ])
279    >>> model.wells.WI
280    array([ 2.,  2.,  2., nan])
281    >>> SS, PP = model.sim(.1, 3, np.zeros(model.Nxy), pbar=False)
282    >>> model.wells.rates_by_well.round(12)  # the aquifer supplies the producer
283    array([[ 1.,  1.,  1.],
284           [-1., -1., -1.]])
285
286    A cell's *boundary faces* are those across which the neighbour is inactive
287    (ref `minires.ResSim.active`) or outside the grid; it must have at
288    least one, and be active itself. Each contributes the half-cell
289    transmissibility, $ C \\, k \\, h_⊥ / (h_∥ / 2) $ (compare the whole-cell
290    one of `minires.ResSim.TPFA`), so that the aquifer pressure is
291    imposed *at the face* -- a Dirichlet condition, its flux discretized as
292    the interior ones are. A corner cell, with two boundary faces, gets both --
293    unless `faces` (a string of compass directions) leaves one out, as it must
294    when that edge is sealed, or on a 1D strip, whose every cell is a boundary
295    cell to the north and south. On a curved outline, keep them all.
296
297    .. note:: The mobility is the cell's total one, $ λ_t(S) $, as for any well.
298
299        Whereas a boundary face of the TPFA scheme would upwind it from the
300        aquifer side (water, $ S = 1 $). The difference is the injector's
301        well model, no more; the influx is water either way.
302        For an aquifer of a given *strength* -- a productivity index, $ J $,
303        as in the Fetkovich model -- give `WI` directly instead, or scale this.
304    """
305    xy = np.asarray(xy, float).reshape((-1, 2))
306    ix, iy = model.xy2sub(*xy.T)
307    # The number of boundary faces (among those selected) in each direction
308    n = boundary_faces(model, xy, faces).astype(int)  # NB: `bool + bool` is an `or`
309    nx, ny = n[:, :2].sum(1), n[:, 2:].sum(1)
310    assert (nx + ny > 0).all(), (
311        "An aquifer cell must lie on the boundary: have a face to an inactive"
312        " cell, or to outside the grid (ref `aquifer_WI`)."
313    )
314    kx, ky = model.K[0][ix, iy], model.K[1][ix, iy]
315    return model.cdarcy * 2 * (nx * kx * model.hy / model.hx + ny * ky * model.hx / model.hy)
316
317
318def boundary_faces(model: "ResSim", xy: Any, faces: str = "WESN") -> np.ndarray:
319    """Which faces of the cells at `xy` are *boundary* faces: to an inactive
320    cell (ref `minires.ResSim.active`), or to outside the grid.
321
322    Boolean, `(nCells, 4)`, the columns being the directions W, E, S, N --
323    of which `faces` (a string of them) selects the ones considered at all.
324    Serves `aquifer_WI`, and `minires.plotting.Plot2D.plt_faces`.
325
326    >>> from minires import ResSim
327    >>> model = ResSim(Lx=1, Ly=1, Nx=4, Ny=4)
328    >>> boundary_faces(model, [[0, 0], [0, .5], [.5, .5]]).astype(int)
329    array([[1, 0, 1, 0],
330           [1, 0, 0, 0],
331           [0, 0, 0, 0]])
332    """
333    xy = np.asarray(xy, float).reshape((-1, 2))
334    ix, iy = model.xy2sub(*xy.T)
335    act = np.pad(model.active, 1, constant_values=False)  # off-grid ⇒ inactive
336    assert act[ix + 1, iy + 1].all(), "The cells must be active (ref `active`)."
337    ix, iy = ix + 1, iy + 1  # (in the padded mask)
338    nbrs = [act[ix - 1, iy], act[ix + 1, iy], act[ix, iy - 1], act[ix, iy + 1]]
339    selected = np.array([d in faces for d in "WESN"])
340    return ~np.stack(nbrs, -1) & selected
341
342
343@dataclass
344class Wells(AlignedRepr):
345    """The wells of a `minires.ResSim`: the flat, per-completion arrays.
346
347    These arrays *are* the configuration (ref `from_records`), and they are
348    meant to be written to, as an ensemble or optimisation loop does:
349
350    >>> from minires import ResSim
351    >>> model = ResSim(Lx=1, Ly=1, Nx=16, Ny=16,
352    ...                wells=Wells(xy=[[0, 0], [1, 1]], rates=[[1], [-1]]))
353    >>> model.wells.rates = [[2], [-2]]
354    >>> model.wells.nComp
355    2
356
357    Assigning one normalizes it (`__setattr__`): the positions get snapped onto
358    the grid nodes, the schedules reshaped to `(nComp, nTime)`, and so on.
359    A `Wells` may also be built before it has a grid to snap to -- as above,
360    where it is the assignment to the model that binds (and snaps) it.
361
362    Ref `from_records` for the convenient, well-shaped way to configure them,
363    and `minires.ResSim.wells` for the attribute that holds them.
364    """
365
366    # Dont use dataclass repr
367    __repr__ = AlignedRepr.__repr__
368
369    # NB: the array attributes are typed `Any` since `__setattr__` normalizes
370    # whatever array-like (nested lists, scalars) is assigned to them.
371    xy: Any = None
372    """Array of shape `(nComp, 2)` of x- and y-coords for the completions.
373
374    Values should be betwen `0` and `Lx` or `Ly`. Empty, shape `(0, 2)`, if
375    there are no wells.
376
377    .. warning:: The wells get co-located with grid nodes, ref `xy2sub`.
378
379        This is a design choice, not a mathematical necessity.
380        An alternative would be to distribute them over nearby nodes.
381    """
382    rates: Any = None
383    """Array of shape `(nComp, nTime)` -- or `(nComp, 1)` if constant-in-time.
384    Ref `from_records` for the convenient way to set this (and the other specs).
385
386    **Signed**: a positive rate *injects* (water), a negative one *produces*
387    (at the well cell's fractional flow). There is no other distinction
388    between injectors and producers, anywhere in the model.
389    An **areal** rate, i.e. volumetric *per unit thickness* -- ref
390    `minires.ResSim.cdarcy`.
391
392    .. note:: With `ct == 0` the rates must sum to 0 at each time index.
393
394        This is asserted (unless a well is BHP-controlled), otherwise the model
395        would silently input the deficit from the SW corner.
396
397    .. note:: Prefer `0` as the (ignored) fill value where control is by `bhp`.
398
399        The array being shared by all the wells, an entry must be given there
400        regardless.
401    """
402    bhp: Any = None
403    """Bottom-hole pressures for the wells. `None`, or an array shaped like
404    `rates`, i.e. `(nComp, nTime)` -- or `(nComp, 1)` if constant-in-time.
405
406    A well whose entry is finite is **BHP-controlled** at that time. This is
407    solved *simultaneously* with the pressure field (not lagged by a time step):
408    `minires.ResSim.TPFA` puts $ WI λ_t $ on its diagonal and
409    $ WI λ_t \\, p_\\mathrm{bh} $ on its right-hand side.
410    Only once $ p $ is known is the resulting rate folded into the source field,
411    `_Q`, for the transport step -- ref `minires.ResSim.assemble_wells`
412    and `minires.ResSim.realize_bhp`.
413
414    Entries left as `nan` -- which is the default, for every well -- keep the
415    well rate-controlled, per `rates`. So the two mechanisms can be mixed
416    freely, both across wells and in time.
417
418    Requires `WI` (finite, for the BHP-controlled wells). The realized
419    rates are recorded in `actual_rates`, and the corresponding entries of
420    `rates` are ignored (it may be left as `None`).
421
422    .. note:: BHP control also *anchors the pressure*.
423
424        With `ct == 0` the pressure equation is otherwise a pure-Neumann
425        problem: solvable only up to a constant (which
426        `minires.ResSim.TPFA` pins arbitrarily, ref article p. 13), and
427        only if injection balances production. A single BHP
428        well lifts both restrictions -- the level is then set by
429        $ p_\\mathrm{bh} $, and the voidage by the well model.
430
431    .. warning:: A BHP well's flow direction is *emergent*, not declared.
432
433        It flows whichever way $ p_\\mathrm{bh} $ vs. its cell pressure dictates
434        (ref `minires.ResSim.realize_bhp`) -- and, like any inflow, an
435        inflow through a BHP well injects *water*. Since a reversal mid-`sim`
436        may nonetheless be a surprise, a `UserWarning` is emitted when a BHP
437        well's realized rate flips sign between steps of `sim`.
438        Nor is there native switching of control modes (e.g. rate control with
439        a BHP limit), but `minires.ResSim.well_controls` can approximate it.
440    """
441    WI: Any = None
442    """Well indices: `None`, or an array of shape `(nComp,)`, `nan` allowed.
443
444    Compute $ WI $ with `peaceman_WI`, or set it directly
445    (it need not come from any particular formula).
446    A well whose entry is `nan` (or all of them, if `None`) has no well model:
447    its `actual_bhp` is `nan`, and BHP control (`bhp`) unavailable.
448
449    The well index is the *sub-grid* well model (ref the "Theory" section of
450    `minires.wells`), relating a well's (signed) flow rate to its drawdown,
451    $$ q = WI \\, λ_t \\, (p_\\mathrm{bh} - p_\\mathrm{cell}) \\,,$$
452    with $ λ_t $ the total mobility (ref `minires.fluids.Fluid.RelPerm`) of the
453    well's cell. Re-arranging,
454    $$ p_\\mathrm{bh} = p_\\mathrm{cell} + q / (WI \\, λ_t) \\,, $$
455    which is what `minires.ResSim.bhp` computes.
456
457    .. warning:: The drawdown is *not* a fixed offset.
458
459        It is not to be calibrated away once and for all: being
460        $ q / (WI \\, λ_t) $, it tracks the mobility --
461        which, in a waterflood, dips as the front arrives (by half, for equal
462        viscosities). So the gap doubles at breakthrough: precisely when the
463        well is most interesting.
464    """
465    group: Any = None
466    """Which well each completion belongs to: `None`, or an int array of shape
467    `(nComp,)` whose values index the wells, i.e. `names`.
468
469    The model itself is indifferent to it: the equations are assembled per
470    *completion* (ref `minires.ResSim.assemble_wells`), and the arrays --
471    `xy`, `rates`, `actual_rates`, ... -- are all indexed likewise. The grouping
472    is what lets the *reporting* speak of wells nonetheless: ref
473    `rates_by_well`, and the labels of `minires.plotting.Plot2D.plt_field`.
474    """
475    names: Any = None
476    """Names of wells (*not* completions): `None`, or a list of `nWell` strings."""
477
478    actual_rates: Any = None
479    """The *realized* well rates: array of shape `(nComp, nSteps)`. Signed.
480
481    Mostly used as a diagnostic in case of `bhp`. But even for
482    rate-control it only coincides with `rates` up to broadcasting
483    and assuming `minires.ResSim.well_controls` did not override it.
484    """
485    actual_bhp: Any = None
486    """Like `actual_rates`, but the bottom-hole pressures.
487
488    `nan` wherever the well index (`WI`) is unset.
489    """
490
491    _grid: Any = None
492    """The model these wells are in -- set when assigned to it. Used only for
493    the grid geometry (which is why it is typed as such, ref `_bind`)."""
494
495    def __setattr__(self, key: str, val: Any) -> None:
496        # NB: the single normalization layer for the wells -- `from_records`
497        # routes its assignments through it rather than writing past it.
498        if key == "xy":
499            # Completion positions -- collocate at some node
500            val = (
501                np.zeros((0, 2))
502                if val is None
503                else np.array(val, float).reshape((-1, 2))
504            )
505            if self._grid is not None:
506                for i, (x, y) in enumerate(val):
507                    val[i] = self._grid.ind2xy(self._grid.xy2ind(x, y))
508        elif val is not None:
509            # Rates and/or pressures
510            if key in ["rates", "bhp"]:
511                val = np.array(val, float).reshape((self.nComp, -1))
512            # Well indices
513            elif key == "WI":
514                val = np.broadcast_to(np.asarray(val, float).ravel(), self.nComp).copy()
515            # Completion-to-well map
516            elif key == "group":
517                val = np.asarray(val, int).reshape(self.nComp)
518        super().__setattr__(key, val)
519
520    def _bind(self, grid: Any) -> None:
521        """Attach to `grid` (a `minires.grid.Grid2D`, i.e. the model).
522
523        Whereupon the completions snap onto its nodes -- which an unbound
524        `Wells`, having no grid to snap to, could not do.
525        """
526        self._grid = grid
527        self.xy = self.xy  # re-normalize, now that there is a grid
528
529    nComp = property(lambda self: len(self.xy))
530    """Num. of *completions*, i.e. rows of `xy`, which is what the model
531    actually solves for. Several completions may compose a single well
532    (ref `group`, `well_path`)."""
533
534    nWell = property(
535        lambda self: self.nComp if self.group is None else 1 + int(self.group.max())
536    )
537    """Num. of *wells*, i.e. groups of completions (ref `group`)."""
538
539    @property
540    def rates_by_well(self) -> np.ndarray:
541        """`actual_rates`, summed over each well's completions: `(nWell, nSteps)`."""
542        group = np.arange(self.nComp) if self.group is None else self.group
543        out = np.zeros((self.nWell, self.actual_rates.shape[1]))
544        np.add.at(out, group, self.actual_rates)  # NB: `+=` would skip the dupes
545        return out
546
547    @property
548    def signs(self) -> np.ndarray:
549        """The sign (`+1` inject, `-1` produce, `0` unknown) of each well's rate.
550
551        Read off the *spec*, `rates`, summed over time (`nan` entries --
552        which a BHP-controlled well may well have -- being skipped). Wells left
553        undecided by it, i.e. those with no spec or a vanishing one (as under
554        pure BHP control), fall back on the `actual_rates` of the latest `sim`,
555        if there has been one. Only the truly undecided are then `0`.
556        """
557        sgn = np.zeros(self.nComp, int)
558        for rates in [self.rates, self.actual_rates]:
559            if rates is not None:
560                q = np.nansum(rates, axis=1)
561                sgn = np.where(sgn, sgn, (q > 0).astype(int) - (q < 0))
562        return sgn
563
564    def at_time(self, spec: str, absent: float, k: int) -> np.ndarray:
565        """Lookup the `spec` (`"rates"`/`"bhp"`) at time `k`.
566
567        Allows a constant-in-time (singleton) spec, and an unset (`None`) one
568        (for which `absent`, i.e. `0`/`nan` for rate/bhp-controlled wells
569        respectively, is returned).
570        Avoids broadcast (and potentially stale copies) to `(nComp, nSteps)`,
571        which requires `nSteps`, i.e. `sim()`.
572        """
573        arr = getattr(self, spec)
574        if arr is None:
575            return np.full(self.nComp, absent)
576        assert len(arr) == self.nComp, (
577            f"`wells.{spec}` has {len(arr)} rows, but there are"
578            f" {self.nComp} completions (ref `Wells.xy`)."
579        )
580        # Copy, lest `well_controls` write into the spec itself
581        return np.copy(arr[:, k if arr.shape[1] > 1 else 0])
582
583    @classmethod
584    def from_records(cls, model: "ResSim", wells: Any) -> "Wells":
585        """Assemble the flat, per-completion arrays from one record (`dict`) per well.
586
587        This is the convenient way to configure the wells, and assigning the
588        records to `minires.ResSim.wells` is what applies it. Each record
589        may specify
590
591        - `xy`: the well's position, `[x, y]` -- or positions,
592          `[[x, y], ...]`, for a multi-completion well.
593        - `path`: alternatively, a polyline, `[[x, y], ...]`, to be discretized
594          into one completion per cell it traverses, ref
595          `well_path`. Needs `rw`.
596        - `rate`: the well's (signed, ref `Wells.rates`) rate: a scalar, or a
597          schedule (an array over time). Apportioned among its completions in
598          proportion to their well indices (uniformly, absent those).
599        - `bhp`: alternatively (or, in time, additionally) the bottom-hole
600          pressure, ref `Wells.bhp`. Scalar or schedule. Shared -- as a wellbore
601          does -- by all of the well's completions.
602        - `rw`, `skin`: the wellbore radius and skin, whence the well index, via
603          `peaceman_WI`. Without them (or `WI`) the well has no well model.
604        - `WI`: alternatively, the well index itself, given directly.
605        - `aquifer`: `True` makes the well an aquifer contact: its `WI` is then
606          `aquifer_WI` of its cells (which must lie on the boundary), and its
607          `bhp` the aquifer pressure. A string of compass directions instead
608          (`"W"`, `"NE"`, ...) selects the boundary faces that count. Ref `aquifer_WI`.
609        - `name`: for the reporting. Defaults to the well's index.
610
611        The concise cases stay concise -- a position and a rate is a well. A
612        `dict` of records names them by its keys (as `name` does otherwise), and
613        the constructor takes the same thing:
614
615        >>> from minires import ResSim
616        >>> model = ResSim(Lx=1, Ly=1, Nx=16, Ny=16, wells={
617        ...     "I1": dict(xy=[0, 0], rate=+1),
618        ...     "P1": dict(xy=[1, 1], rate=-1, rw=1e-3),
619        ... })
620        >>> model.wells.names
621        ['I1', 'P1']
622        >>> model.wells.rates
623        array([[ 1.],
624               [-1.]])
625        >>> model.wells.WI.round(3)  # `P1` alone asked for a well model
626        array([  nan, 2.498])
627
628        A `path` becomes several completions of a single well, whose rate it
629        shares out (ref `well_path`) and whose name they share:
630
631        >>> model = ResSim(Lx=1, Ly=1, Nx=10, Ny=10, wells=[
632        ...     dict(name="I1", path=[[.05, .05], [.45, .05]], rate=+1, rw=1e-2),
633        ...     dict(name="P1", xy=[.95, .95], rate=-1),
634        ... ])
635        >>> model.wells.nComp, model.wells.nWell
636        (6, 2)
637        >>> model.wells.group
638        array([0, 0, 0, 0, 0, 1])
639        >>> model.wells.rates.ravel().round(3)
640        array([ 0.125,  0.25 ,  0.25 ,  0.25 ,  0.125, -1.   ])
641
642        Schedules and control modes may be mixed freely across the wells: a
643        constant is broadcast to the length of the longest schedule, while a
644        spec that no well varies in time stays a singleton, which
645        `minires.ResSim.well_controls` reads at any `k`. A BHP-controlled
646        well leaves the (ignored) `0` in `Wells.rates`, and `nan` marks the
647        rate-controlled ones in `Wells.bhp` -- the conventions that the specs,
648        being shared arrays, call for (ref `Wells.rates`).
649
650        .. note:: A well must be given a control (`rate` and/or `bhp`).
651
652            `rate=0` shuts it in. This is deliberate: an uncontrolled well would
653            silently be a shut one.
654
655        .. note:: The records are *not* retained.
656
657            The arrays they produce are the whole of the configuration, so
658            there is nothing to fall out of step with a subsequent edit of them
659            (which is what the `repr` therefore reports).
660        """
661        keys = ("name", "xy", "path", "rate", "bhp", "rw", "skin", "WI", "aquifer")
662        if isinstance(wells, dict):
663            wells = [dict(spec, name=name) for name, spec in wells.items()]
664        names, xy, WI, group, rates, bhp = [], [], [], [], [], []
665        specified: set = set()
666        for i, well in enumerate(wells or []):
667            spec = dict(well)
668            if unknown := set(spec) - set(keys):
669                raise TypeError(
670                    f"Unknown key(s) in the spec of well {i}: {sorted(unknown)}."
671                    f" Valid ones: {list(keys)}."
672                )
673            specified |= set(spec)
674            name = str(spec.pop("name", i))
675            rw, skin = spec.pop("rw", None), spec.pop("skin", 0.0)
676            aquifer = spec.pop("aquifer", False)
677
678            # Completions: their positions, and their well indices
679            if (path := spec.pop("path", None)) is not None:
680                assert "xy" not in spec, (
681                    f"Well '{name}': give it `xy` or `path`, not both."
682                )
683                assert rw is not None, f"Well '{name}': a `path` requires `rw`."
684                assert not aquifer, f"Well '{name}': an aquifer is given by `xy`."
685                _xy, _WI, _ = well_path(model, path, rw, skin)
686            else:
687                assert "xy" in spec, f"Well '{name}': give it an `xy` (or a `path`)."
688                _xy = np.array(spec.pop("xy"), float).reshape((-1, 2))
689                if aquifer:
690                    assert rw is None, f"Well '{name}': an aquifer has no `rw`."
691                    _WI = aquifer_WI(model, _xy, "WESN" if aquifer is True else aquifer)
692                elif rw is not None:
693                    _WI = peaceman_WI(model, _xy, rw, skin)
694                else:
695                    _WI = np.full(len(_xy), np.nan)
696            if (given := spec.pop("WI", None)) is not None:
697                _WI = np.broadcast_to(np.asarray(given, float).ravel(), len(_xy)).copy()
698            nc = len(_xy)
699
700            # Apportion the rate by well index -- the standard, ref `well_path`
701            alloc = np.full(nc, 1 / nc)
702            if nc > 1 and np.isfinite(_WI).all() and _WI.sum() > 0:
703                alloc = _WI / _WI.sum()
704
705            # Controls. NB: the BHP is shared by the completions, the rate split
706            rate, p_bh = spec.pop("rate", None), spec.pop("bhp", None)
707            assert rate is not None or p_bh is not None, (
708                f"Well '{name}' has no control: give it a `rate`"
709                " (`0` shuts it in), or a `bhp`."
710            )
711            rate = 0.0 if rate is None else rate
712            p_bh = np.nan if p_bh is None else p_bh
713            rates.append(np.outer(alloc, np.ravel(rate)))
714            bhp.append(np.outer(np.ones(nc), np.ravel(p_bh)))
715
716            names.append(name)
717            xy.append(_xy)
718            WI.append(_WI)
719            group.append(np.full(nc, i))
720
721        if not names:
722            return cls()
723
724        def stack(specs):
725            """Stack the wells' `(nComp_i, nTime_i)` specs, widening the constants.
726
727            NB: `at_time` broadcasts a *wholly* singleton spec, but the array is
728            shared, so a well held constant beside a scheduled one is widened here.
729            """
730            nTime = max(spec.shape[1] for spec in specs)
731            assert all(spec.shape[1] in [1, nTime] for spec in specs), (
732                "The wells' schedules must be of equal length (or constant):"
733                f" got {sorted({spec.shape[1] for spec in specs})}."
734            )
735            return np.vstack([np.broadcast_to(s, (len(s), nTime)) for s in specs])
736
737        WI = np.concatenate(WI)
738        return cls(
739            # NB: `xy` first -- it is what defines `nComp`, by which the
740            # `__setattr__` normalization shapes the others.
741            xy     = np.vstack(xy),
742            # Leave a spec unset (`None`) if no well made use of it
743            rates  = stack(rates) if "rate" in specified else None,
744            bhp    = stack(bhp) if "bhp" in specified else None,
745            WI     = WI if np.isfinite(WI).any() else None,
746            group  = np.concatenate(group),
747            names  = names,
748        )  # fmt: off