TPFA_ResSim.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(`TPFA_ResSim.grid.Grid2D`), and -- for the well index -- on the permeability
  5`K`. What couples the wells to the *fluids* (the mobilities of
  6`TPFA_ResSim.ResSim.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 `TPFA_ResSim.wells.Wells.rates`), and
 48under BHP control not even by that, the direction being left to the pressures
 49(ref `TPFA_ResSim.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 `TPFA_ResSim.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 `TPFA_ResSim.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 `TPFA_ResSim.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 `TPFA_ResSim.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 TPFA_ResSim._repr import AlignedRepr
130
131if TYPE_CHECKING:
132    from TPFA_ResSim 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 `TPFA_ResSim.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 TPFA_ResSim 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. `TPFA_ResSim.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: `TPFA_ResSim.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 TPFA_ResSim 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 `TPFA_ResSim.ResSim.TPFA` (Listing 1) is not set up for. Use
222        `TPFA_ResSim.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
258@dataclass
259class Wells(AlignedRepr):
260    """The wells of a `TPFA_ResSim.ResSim`: the flat, per-completion arrays.
261
262    These arrays *are* the configuration (ref `from_records`), and they are
263    meant to be written to, as an ensemble or optimisation loop does:
264
265    >>> from TPFA_ResSim import ResSim
266    >>> model = ResSim(Lx=1, Ly=1, Nx=16, Ny=16,
267    ...                wells=Wells(xy=[[0, 0], [1, 1]], rates=[[1], [-1]]))
268    >>> model.wells.rates = [[2], [-2]]
269    >>> model.wells.nComp
270    2
271
272    Assigning one normalizes it (`__setattr__`): the positions get snapped onto
273    the grid nodes, the schedules reshaped to `(nComp, nTime)`, and so on.
274    A `Wells` may also be built before it has a grid to snap to -- as above,
275    where it is the assignment to the model that binds (and snaps) it.
276
277    Ref `from_records` for the convenient, well-shaped way to configure them,
278    and `TPFA_ResSim.ResSim.wells` for the attribute that holds them.
279    """
280
281    # Dont use dataclass repr
282    __repr__ = AlignedRepr.__repr__
283
284    # NB: the array attributes are typed `Any` since `__setattr__` normalizes
285    # whatever array-like (nested lists, scalars) is assigned to them.
286    xy: Any = None
287    """Array of shape `(nComp, 2)` of x- and y-coords for the completions.
288
289    Values should be betwen `0` and `Lx` or `Ly`. Empty, shape `(0, 2)`, if
290    there are no wells.
291
292    .. warning:: The wells get co-located with grid nodes, ref `xy2sub`.
293
294        This is a design choice, not a mathematical necessity.
295        An alternative would be to distribute them over nearby nodes.
296    """
297    rates: Any = None
298    """Array of shape `(nComp, nTime)` -- or `(nComp, 1)` if constant-in-time.
299    Ref `from_records` for the convenient way to set this (and the other specs).
300
301    **Signed**: a positive rate *injects* (water), a negative one *produces*
302    (at the well cell's fractional flow). There is no other distinction
303    between injectors and producers, anywhere in the model.
304    An **areal** rate, i.e. volumetric *per unit thickness* -- ref
305    `TPFA_ResSim.ResSim.cdarcy`.
306
307    .. note:: With `ct == 0` the rates must sum to 0 at each time index.
308
309        This is asserted (unless a well is BHP-controlled), otherwise the model
310        would silently input the deficit from the SW corner.
311
312    .. note:: Prefer `0` as the (ignored) fill value where control is by `bhp`.
313
314        The array being shared by all the wells, an entry must be given there
315        regardless.
316    """
317    bhp: Any = None
318    """Bottom-hole pressures for the wells. `None`, or an array shaped like
319    `rates`, i.e. `(nComp, nTime)` -- or `(nComp, 1)` if constant-in-time.
320
321    A well whose entry is finite is **BHP-controlled** at that time. This is
322    solved *simultaneously* with the pressure field (not lagged by a time step):
323    `TPFA_ResSim.ResSim.TPFA` puts $ WI λ_t $ on its diagonal and
324    $ WI λ_t \\, p_\\mathrm{bh} $ on its right-hand side.
325    Only once $ p $ is known is the resulting rate folded into the source field,
326    `_Q`, for the transport step -- ref `TPFA_ResSim.ResSim.assemble_wells`
327    and `TPFA_ResSim.ResSim.realize_bhp`.
328
329    Entries left as `nan` -- which is the default, for every well -- keep the
330    well rate-controlled, per `rates`. So the two mechanisms can be mixed
331    freely, both across wells and in time.
332
333    Requires `WI` (finite, for the BHP-controlled wells). The realized
334    rates are recorded in `actual_rates`, and the corresponding entries of
335    `rates` are ignored (it may be left as `None`).
336
337    .. note:: BHP control also *anchors the pressure*.
338
339        With `ct == 0` the pressure equation is otherwise a pure-Neumann
340        problem: solvable only up to a constant (which
341        `TPFA_ResSim.ResSim.TPFA` pins arbitrarily, ref article p. 13), and
342        only if injection balances production. A single BHP
343        well lifts both restrictions -- the level is then set by
344        $ p_\\mathrm{bh} $, and the voidage by the well model.
345
346    .. warning:: A BHP well's flow direction is *emergent*, not declared.
347
348        It flows whichever way $ p_\\mathrm{bh} $ vs. its cell pressure dictates
349        (ref `TPFA_ResSim.ResSim.realize_bhp`) -- and, like any inflow, an
350        inflow through a BHP well injects *water*. Since a reversal mid-`sim`
351        may nonetheless be a surprise, a `UserWarning` is emitted when a BHP
352        well's realized rate flips sign between steps of `sim`.
353        Nor is there native switching of control modes (e.g. rate control with
354        a BHP limit), but `TPFA_ResSim.ResSim.well_controls` can approximate it.
355    """
356    WI: Any = None
357    """Well indices: `None`, or an array of shape `(nComp,)`, `nan` allowed.
358
359    Compute $ WI $ with `peaceman_WI`, or set it directly
360    (it need not come from any particular formula).
361    A well whose entry is `nan` (or all of them, if `None`) has no well model:
362    its `actual_bhp` is `nan`, and BHP control (`bhp`) unavailable.
363
364    The well index is the *sub-grid* well model (ref the "Theory" section of
365    `TPFA_ResSim.wells`), relating a well's (signed) flow rate to its drawdown,
366    $$ q = WI \\, λ_t \\, (p_\\mathrm{bh} - p_\\mathrm{cell}) \\,,$$
367    with $ λ_t $ the total mobility (ref `TPFA_ResSim.ResSim.RelPerm`) of the
368    well's cell. Re-arranging,
369    $$ p_\\mathrm{bh} = p_\\mathrm{cell} + q / (WI \\, λ_t) \\,, $$
370    which is what `TPFA_ResSim.ResSim.bhp` computes.
371
372    .. warning:: The drawdown is *not* a fixed offset.
373
374        It is not to be calibrated away once and for all: being
375        $ q / (WI \\, λ_t) $, it tracks the mobility --
376        which, in a waterflood, dips as the front arrives (by half, for equal
377        viscosities). So the gap doubles at breakthrough: precisely when the
378        well is most interesting.
379    """
380    group: Any = None
381    """Which well each completion belongs to: `None`, or an int array of shape
382    `(nComp,)` whose values index the wells, i.e. `names`.
383
384    The model itself is indifferent to it: the equations are assembled per
385    *completion* (ref `TPFA_ResSim.ResSim.assemble_wells`), and the arrays --
386    `xy`, `rates`, `actual_rates`, ... -- are all indexed likewise. The grouping
387    is what lets the *reporting* speak of wells nonetheless: ref
388    `rates_by_well`, and the labels of `TPFA_ResSim.plotting.Plot2D.plt_field`.
389    """
390    names: Any = None
391    """Names of wells (*not* completions): `None`, or a list of `nWell` strings."""
392
393    actual_rates: Any = None
394    """The *realized* well rates: array of shape `(nComp, nSteps)`. Signed.
395
396    Mostly used as a diagnostic in case of `bhp`. But even for
397    rate-control it only coincides with `rates` up to broadcasting
398    and assuming `TPFA_ResSim.ResSim.well_controls` did not override it.
399    """
400    actual_bhp: Any = None
401    """Like `actual_rates`, but the bottom-hole pressures.
402
403    `nan` wherever the well index (`WI`) is unset.
404    """
405
406    _grid: Any = None
407    """The model these wells are in -- set when assigned to it. Used only for
408    the grid geometry (which is why it is typed as such, ref `_bind`)."""
409
410    def __setattr__(self, key: str, val: Any) -> None:
411        # NB: the single normalization layer for the wells -- `from_records`
412        # routes its assignments through it rather than writing past it.
413        if key == "xy":
414            # Completion positions -- collocate at some node
415            val = (
416                np.zeros((0, 2))
417                if val is None
418                else np.array(val, float).reshape((-1, 2))
419            )
420            if self._grid is not None:
421                for i, (x, y) in enumerate(val):
422                    val[i] = self._grid.ind2xy(self._grid.xy2ind(x, y))
423        elif val is not None:
424            # Rates and/or pressures
425            if key in ["rates", "bhp"]:
426                val = np.array(val, float).reshape((self.nComp, -1))
427            # Well indices
428            elif key == "WI":
429                val = np.broadcast_to(np.asarray(val, float).ravel(), self.nComp).copy()
430            # Completion-to-well map
431            elif key == "group":
432                val = np.asarray(val, int).reshape(self.nComp)
433        super().__setattr__(key, val)
434
435    def _bind(self, grid: Any) -> None:
436        """Attach to `grid` (a `TPFA_ResSim.grid.Grid2D`, i.e. the model).
437
438        Whereupon the completions snap onto its nodes -- which an unbound
439        `Wells`, having no grid to snap to, could not do.
440        """
441        self._grid = grid
442        self.xy = self.xy  # re-normalize, now that there is a grid
443
444    nComp = property(lambda self: len(self.xy))
445    """Num. of *completions*, i.e. rows of `xy`, which is what the model
446    actually solves for. Several completions may compose a single well
447    (ref `group`, `well_path`)."""
448
449    nWell = property(
450        lambda self: self.nComp if self.group is None else 1 + int(self.group.max())
451    )
452    """Num. of *wells*, i.e. groups of completions (ref `group`)."""
453
454    @property
455    def rates_by_well(self) -> np.ndarray:
456        """`actual_rates`, summed over each well's completions: `(nWell, nSteps)`."""
457        group = np.arange(self.nComp) if self.group is None else self.group
458        out = np.zeros((self.nWell, self.actual_rates.shape[1]))
459        np.add.at(out, group, self.actual_rates)  # NB: `+=` would skip the dupes
460        return out
461
462    @property
463    def signs(self) -> np.ndarray:
464        """The sign (`+1` inject, `-1` produce, `0` unknown) of each well's rate.
465
466        Read off the *spec*, `rates`, summed over time (`nan` entries --
467        which a BHP-controlled well may well have -- being skipped). Wells left
468        undecided by it, i.e. those with no spec or a vanishing one (as under
469        pure BHP control), fall back on the `actual_rates` of the latest `sim`,
470        if there has been one. Only the truly undecided are then `0`.
471        """
472        sgn = np.zeros(self.nComp, int)
473        for rates in [self.rates, self.actual_rates]:
474            if rates is not None:
475                q = np.nansum(rates, axis=1)
476                sgn = np.where(sgn, sgn, (q > 0).astype(int) - (q < 0))
477        return sgn
478
479    def at_time(self, spec: str, absent: float, k: int) -> np.ndarray:
480        """Lookup the `spec` (`"rates"`/`"bhp"`) at time `k`.
481
482        Allows a constant-in-time (singleton) spec, and an unset (`None`) one
483        (for which `absent`, i.e. `0`/`nan` for rate/bhp-controlled wells
484        respectively, is returned).
485        Avoids broadcast (and potentially stale copies) to `(nComp, nSteps)`,
486        which requires `nSteps`, i.e. `sim()`.
487        """
488        arr = getattr(self, spec)
489        if arr is None:
490            return np.full(self.nComp, absent)
491        assert len(arr) == self.nComp, (
492            f"`wells.{spec}` has {len(arr)} rows, but there are"
493            f" {self.nComp} completions (ref `Wells.xy`)."
494        )
495        # Copy, lest `well_controls` write into the spec itself
496        return np.copy(arr[:, k if arr.shape[1] > 1 else 0])
497
498    @classmethod
499    def from_records(cls, model: "ResSim", wells: Any) -> "Wells":
500        """Assemble the flat, per-completion arrays from one record (`dict`) per well.
501
502        This is the convenient way to configure the wells, and assigning the
503        records to `TPFA_ResSim.ResSim.wells` is what applies it. Each record
504        may specify
505
506        - `xy`: the well's position, `[x, y]` -- or positions,
507          `[[x, y], ...]`, for a multi-completion well.
508        - `path`: alternatively, a polyline, `[[x, y], ...]`, to be discretized
509          into one completion per cell it traverses, ref
510          `well_path`. Needs `rw`.
511        - `rate`: the well's (signed, ref `Wells.rates`) rate: a scalar, or a
512          schedule (an array over time). Apportioned among its completions in
513          proportion to their well indices (uniformly, absent those).
514        - `bhp`: alternatively (or, in time, additionally) the bottom-hole
515          pressure, ref `Wells.bhp`. Scalar or schedule. Shared -- as a wellbore
516          does -- by all of the well's completions.
517        - `rw`, `skin`: the wellbore radius and skin, whence the well index, via
518          `peaceman_WI`. Without them (or `WI`) the well has no well model.
519        - `WI`: alternatively, the well index itself, given directly.
520        - `name`: for the reporting. Defaults to the well's index.
521
522        The concise cases stay concise -- a position and a rate is a well. A
523        `dict` of records names them by its keys (as `name` does otherwise), and
524        the constructor takes the same thing:
525
526        >>> from TPFA_ResSim import ResSim
527        >>> model = ResSim(Lx=1, Ly=1, Nx=16, Ny=16, wells={
528        ...     "I1": dict(xy=[0, 0], rate=+1),
529        ...     "P1": dict(xy=[1, 1], rate=-1, rw=1e-3),
530        ... })
531        >>> model.wells.names
532        ['I1', 'P1']
533        >>> model.wells.rates
534        array([[ 1.],
535               [-1.]])
536        >>> model.wells.WI.round(3)  # `P1` alone asked for a well model
537        array([  nan, 2.498])
538
539        A `path` becomes several completions of a single well, whose rate it
540        shares out (ref `well_path`) and whose name they share:
541
542        >>> model = ResSim(Lx=1, Ly=1, Nx=10, Ny=10, wells=[
543        ...     dict(name="I1", path=[[.05, .05], [.45, .05]], rate=+1, rw=1e-2),
544        ...     dict(name="P1", xy=[.95, .95], rate=-1),
545        ... ])
546        >>> model.wells.nComp, model.wells.nWell
547        (6, 2)
548        >>> model.wells.group
549        array([0, 0, 0, 0, 0, 1])
550        >>> model.wells.rates.ravel().round(3)
551        array([ 0.125,  0.25 ,  0.25 ,  0.25 ,  0.125, -1.   ])
552
553        Schedules and control modes may be mixed freely across the wells: a
554        constant is broadcast to the length of the longest schedule, while a
555        spec that no well varies in time stays a singleton, which
556        `TPFA_ResSim.ResSim.well_controls` reads at any `k`. A BHP-controlled
557        well leaves the (ignored) `0` in `Wells.rates`, and `nan` marks the
558        rate-controlled ones in `Wells.bhp` -- the conventions that the specs,
559        being shared arrays, call for (ref `Wells.rates`).
560
561        .. note:: A well must be given a control (`rate` and/or `bhp`).
562
563            `rate=0` shuts it in. This is deliberate: an uncontrolled well would
564            silently be a shut one.
565
566        .. note:: The records are *not* retained.
567
568            The arrays they produce are the whole of the configuration, so
569            there is nothing to fall out of step with a subsequent edit of them
570            (which is what the `repr` therefore reports).
571        """
572        keys = ("name", "xy", "path", "rate", "bhp", "rw", "skin", "WI")
573        if isinstance(wells, dict):
574            wells = [dict(spec, name=name) for name, spec in wells.items()]
575        names, xy, WI, group, rates, bhp = [], [], [], [], [], []
576        specified: set = set()
577        for i, well in enumerate(wells or []):
578            spec = dict(well)
579            if unknown := set(spec) - set(keys):
580                raise TypeError(
581                    f"Unknown key(s) in the spec of well {i}: {sorted(unknown)}."
582                    f" Valid ones: {list(keys)}."
583                )
584            specified |= set(spec)
585            name = str(spec.pop("name", i))
586            rw, skin = spec.pop("rw", None), spec.pop("skin", 0.0)
587
588            # Completions: their positions, and their well indices
589            if (path := spec.pop("path", None)) is not None:
590                assert "xy" not in spec, (
591                    f"Well '{name}': give it `xy` or `path`, not both."
592                )
593                assert rw is not None, f"Well '{name}': a `path` requires `rw`."
594                _xy, _WI, _ = well_path(model, path, rw, skin)
595            else:
596                assert "xy" in spec, f"Well '{name}': give it an `xy` (or a `path`)."
597                _xy = np.array(spec.pop("xy"), float).reshape((-1, 2))
598                _WI = (
599                    np.full(len(_xy), np.nan)
600                    if rw is None
601                    else peaceman_WI(model, _xy, rw, skin)
602                )
603            if (given := spec.pop("WI", None)) is not None:
604                _WI = np.broadcast_to(np.asarray(given, float).ravel(), len(_xy)).copy()
605            nc = len(_xy)
606
607            # Apportion the rate by well index -- the standard, ref `well_path`
608            alloc = np.full(nc, 1 / nc)
609            if nc > 1 and np.isfinite(_WI).all() and _WI.sum() > 0:
610                alloc = _WI / _WI.sum()
611
612            # Controls. NB: the BHP is shared by the completions, the rate split
613            rate, p_bh = spec.pop("rate", None), spec.pop("bhp", None)
614            assert rate is not None or p_bh is not None, (
615                f"Well '{name}' has no control: give it a `rate`"
616                " (`0` shuts it in), or a `bhp`."
617            )
618            rate = 0.0 if rate is None else rate
619            p_bh = np.nan if p_bh is None else p_bh
620            rates.append(np.outer(alloc, np.ravel(rate)))
621            bhp.append(np.outer(np.ones(nc), np.ravel(p_bh)))
622
623            names.append(name)
624            xy.append(_xy)
625            WI.append(_WI)
626            group.append(np.full(nc, i))
627
628        if not names:
629            return cls()
630
631        def stack(specs):
632            """Stack the wells' `(nComp_i, nTime_i)` specs, widening the constants.
633
634            NB: `at_time` broadcasts a *wholly* singleton spec, but the array is
635            shared, so a well held constant beside a scheduled one is widened here.
636            """
637            nTime = max(spec.shape[1] for spec in specs)
638            assert all(spec.shape[1] in [1, nTime] for spec in specs), (
639                "The wells' schedules must be of equal length (or constant):"
640                f" got {sorted({spec.shape[1] for spec in specs})}."
641            )
642            return np.vstack([np.broadcast_to(s, (len(s), nTime)) for s in specs])
643
644        WI = np.concatenate(WI)
645        return cls(
646            # NB: `xy` first -- it is what defines `nComp`, by which the
647            # `__setattr__` normalization shapes the others.
648            xy     = np.vstack(xy),
649            # Leave a spec unset (`None`) if no well made use of it
650            rates  = stack(rates) if "rate" in specified else None,
651            bhp    = stack(bhp) if "bhp" in specified else None,
652            WI     = WI if np.isfinite(WI).any() else None,
653            group  = np.concatenate(group),
654            names  = names,
655        )  # fmt: off