minires.wells

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