minires

MiniRes is a lean, 2D, two-phase, black-oil, immiscible reservoir simulator using TPFA (two-point flux approximation) written in python.

Governing equations

The simulator solves eqn. (1) and (2) (corresponding to (42) and (43) of the reference paper) :

$$\begin{align} - \nabla \cdot \mathbf{K} \lambda(s) \, \nabla p &= q \,, \tag{1} \cr \phi \frac{\partial s}{\partial t} + \nabla \cdot (f(s)\, \mathbf{v}) &= \frac{q_w}{\rho_w} \,. \tag{2} \end{align}$$

The quantities involved are all 2D-spatial fields, namely

$\phi \in [0, 1]$ is the porosity
$s \in [0, 1]$ is the water saturation
$p$ is the pressure
$v$ is the (volumetric) flow velocity ($\mathbf{v} = \mathbf{v}_o + \mathbf{v}_w$).
$q$ is the sources/sinks
$\rho$ is the density
$\mathbf{K}$ is the (absolute) permeability tensor: the rock's conductivity to flow, here diagonal, $\mathrm{diag}(K_x, K_y)$ per cell.
$\lambda(s)$ is the total mobility (sum of mobilities). Each (relative) mobility is the phase relative permeability divided by the phase viscosity, $\lambda_{\text{phase}} = k_{\text{phase}}/\mu_{\text{phase}}$.
- The relative permeabilities $k_{\text{phase}}(s) \in [0, 1]$ are a constitutive relation, not data: here, Corey (power-law) curves of the saturation rescaled by its residual values, with adjustable exponents and end-points, quadratic by default (ResSim.fluid, a minires.fluids.Fluid). Need not sum to 1.
- $\mu_{\text{phase}}$ is the phase viscosity, here constant.
$f(s) = \lambda_w(s) / \lambda(s) \in [0, 1]$ is the water fractional flow, giving $\mathbf{v}_w = f(s) \, \mathbf{v}$.

The right hand side of (2) is further simplified (relabelled) as $q$, i.e. dropping the $w$ (for "water") subscripts.

Derivation

Single phase

Conservation of mass in a porous ($\phi$) medium is expressed by

$$\frac{∂(\rho \phi)}{∂t} + ∇ \cdot (\rho \mathbf{v}) = q \,. \tag{3}$$

This equation is also called continuity eqn., advection eqn., transport eqn., or even 1st-order wave eqn. (if constant $v$). It says that divergence (or convergence) must be balanced by change in density or porosity, or sinks or sources. If we assume constant porosity, $\phi$, and incompressibility (constant $\rho$), then the time derivative vanishes, yielding the steady-state equation $$\nabla \cdot \mathbf{\mathbf{v}} = \frac{q}{\rho} \,. \tag{4}$$

We now have 1 equation and 2 unknowns (in 2D). Closing the system, Darcy's law provides 2 additional equations and 1 additional unknown, pressure $p$: $$\mathbf{v} = − \frac{\mathbf{K}}{\mu} \nabla u \,, \tag{5}$$ where $u = p - \rho g z \,.$ Analogously to Fourier's heat diffusion and Ohm's conduction law, Darcy's law (5) was initially derived empirically, but can be shown to be a special case of Navier-Stokes' momentum equation. It says that $\mathbf{v}$ is the gradient of the velocity potential, $u$, linearly transformed by the permeability tensor (matrix). Inserting the formula (5) into eqn. (4) yields $$− \nabla \cdot \frac{\mathbf{K}}{\mu} \nabla u = \frac{q}{\rho} \,. \tag{6}$$ which can be solved for $u$. In reservoir engineering, no-flow boundary conditions are most often used, and $u$ is only determined up to a constant (as behoves a potential). Finally, $u$ can be inserted in Darcy's law (5) to yield the (steady-state) velocity.

Two phases

Incompressibility again yields eqn. (4) for the total (volumetric) velocity.
Darcy's law (5) is assumed for each (both) individual phase, with $\mathbf{K}$ replaced by $\mathbf{K} \lambda_{\text{phase}}(s)$.
Neglecting $\nabla z$ (gravity, i.e. hydrostatic pressure), the flow potential, $u$, reduces to the pressure field, $p$.
Summing Darcy's law over the two phases yields $$\mathbf{v} = − \mathbf{K} \lambda (s) \nabla p \,. \tag{7}$$
Repeating the steps right above eqn. (6), one arrives at eqn. (1).
Meanwhile, immiscibility means that conservation of mass (3) must hold for each phase separately, i.e. the density $\rho$ gets replaced by $s_{\text{phase}} \, \rho_{\text{phase}}$, and $\mathbf{v}$ by $\mathbf{v}_{\text{phase}} = f_{\text{phase}}(s)\, \mathbf{v}$, immediately yielding eqn. (2).

How to solve

Equations (1) and (2) are nonlinearly coupled: $s$ and $p$ (yielding $v$ via eqn. (7)) appear in both equations. Trying to solve both equations simultaneously is a nonlinear root-finding problem, requiring Newton iterations and matrix inversions. In this context, it is tempting to use implicit time discretization (like ECLIPSE 100) where $s_{t+1}$ is expressed as a (nonlinear) function of itself, since this would also requires iterations

Here, instead, we apply sequential operator splitting, meaning that the two equations are solved independently, inserting the previous solution of (1) into (2), and vice-versa. Since it yields smaller systems (which can potentially be discretized explicitly) this is faster, but less accurate. When using an explicit (upwind) scheme for the nearly-hyperbolic saturation/transport equation, the strategy is called IMPES (implicit pressure, explicit saturation). The simulator also contains implicit saturation scheme, but it rarely outperforms the explicit one, ref ResSim.saturation_step_implicit.

The spatial discretization is carried out by finite volumes (FV), which is similar to finite differences (FD), but arguably easier to formulate for non-structured (irregular) grids (not our case). For the pressure equation, using only two points two approximate the transmissibility and fluxes at the interfaces is called it is called two-point flux approximation (TPFA); simple, but used widely (nearly default) in oil industry, due to its robustness and efficiency. Consider the equation $$- \nabla \cdot \lambda \nabla u = q \,, \tag{8}$$ where replacing $\lambda \leftarrow \mathbf{K} \lambda(s)$ reproduces eqn. (1), or $\lambda \leftarrow \mathbf{K}/\mu$ and $q \leftarrow q/\rho$ reproduces eqn. (6). FV methods apply the divergence theorem to eqn. (8) to replace point derivatives by integral quantities: interface fluxes and volumetric sources/sinks: $$- \int_{\partial \Omega_i} d x^2 \, \lambda \, (\nabla u) \cdot \mathbf{n} = \int_{\Omega_i} d x^3 \, q \,, \tag{9}$$ where $\Omega_i$ is the domain of cell index $i$, and $\partial \Omega_i$ is its boundary, with normal vector $\mathbf{n}$.

Now, in TPFA we approximate $(\nabla u) \cdot \mathbf{n}$ by a finite difference $$\delta u_{ij} := 2 \frac{u_j - u_i}{\Delta x_i + \Delta x_j}$$ where $u_i, u_i$ are the values of the potential, $u$, at centre of cells $i$ and $j$, which are located either side of the interface $\gamma_{ij}$, which is part of $\partial \Omega_i$. PS: by contrasts, mixed finite-element methods (FEM) do not approximate fluxes over cell edges but considers them unknown. Next, $\lambda$ is approximated by a harmonic average, $\lambda_{ij}$, including weights that account for the distances from the interface to the cell centres. Thus eqn. (9) becomes $$- \sum_j |\gamma_{ij}| \lambda_{ij} \delta u_{ij} = \int_{\Omega_i} d x^3 \, q \,, \tag{10}$$ where the sum is over the indexes $j$ of the interfaces around cell $i$. The left-hand side can be succinctly expressed as $- \sum_j t_{ij} (u_i - u_j)$, where $t_{ij}$ (see above their equation 17) is symmetric. Thus the whole linear system (for all $i$) is symmetric. Moreover, summing over $i$ yields $\sum_{ij} t_{ij} u_i - \sum_{ij} u_j = 0$, meaning that the vector of ones is a null vector for the system (as appropriate for a differential operator), and that $u$ is determined only up to an arbitrary constant (as appropriate for a potential). The constant is fixed, and the system is rendered invertible, by adding to the first element of the diagonal.

The system is thus symmetric positive definite, and is solved by conjugate gradients, preconditioned by a sparse LU factorization that is cached across the time steps (ResSim.cached_precond).

The pressure system need not be factorized afresh each step.

Its matrix changes only through the mobility $λ(s)$, so the factorization of an earlier step is an excellent preconditioner for the current one, at the cost of a back-substitution, and a refactorization only once the iteration stalls. tests/test_precond.py benchmarks.

Units

The units are by default SI (m, s, Pa). But you can switch to metric (m, day, bar, mD) or field-like (ft, day, psi) by changing ResSim.cdarcy.

Compressibility

The above is the default incompressible model, which is what the reference paper treats. Below we derive the so-called slightly compressible approximation, switched on by setting minires.ResSim.ct ($c_t$) $> 0$.

Definition

The compressibility of anything (rock or fluid) is the relative change of its volume by pressure, $c = -\frac{1}{V} \frac{\partial V}{\partial p}$. For rock (pores) it becomes $c_r = \frac{1}{\phi} \frac{\partial \phi}{\partial p}$, while for fluids it is $c_f = \frac{1}{\rho} \frac{\partial \rho}{\partial p}$. With two phases, the fluid in the pores is a mixture, so that the total compressibility is the saturation-weighted sum $c_t = c_r + s_w c_w + s_o c_o$.

Approximation

The model, however, holds it as a single constant, ct, so it is accurate to $O(c_t)$ alone -- the slightly of slightly compressible, which is reasonable for liquids, but not for gas. Thus $\rho \propto e^{c (p - p_0)}$, which the approximation retains only to first order, $\rho \approx \rho_0 [1 + c (p - p_0)]$: a density affine in $p$.

Note that $\rho(p)$ (and thus nonlinearity in $p$) also appears through the source term (wells): a rate fixed at the surface moves a reservoir volume $\propto 1/\rho(p)$. Here it is approximated as a constant (the formation volume factor $B = 1$), or, for a BHP well, as linear in $p$ (Peaceman).

Derivation

Return to the conservation of mass (3), now with $\rho = \rho(p)$ and $\phi = \phi(p)$. By the chain rule and the definitions (constant $c$) above, the accumulation term becomes $$\frac{\partial (\rho \phi)}{\partial t} = \rho \, \phi \, (c_r + c_f) \, \frac{\partial p}{\partial t} \,,$$ where $c_f$ is the compressibility of the fluid filling the pores. Meanwhile, in the flux term, $\nabla \cdot (\rho \mathbf{v}) = \rho \, \nabla \cdot \mathbf{v} + \mathbf{v} \cdot \nabla \rho$, the latter term is $O(c)$ relative to the former (since $\nabla \rho = \rho \, c \, \nabla p$ thanks to constant $c$), and is therefore dropped, an approximation equivalent to the affine one above. Dividing by $\rho$ and inserting Darcy's law (7), we recover eqn. (1) except now with a time derivative: $$\phi \, c_t \frac{\partial p}{\partial t} - \nabla \cdot \mathbf{K} \lambda(s) \, \nabla p = q \,. \tag{11}$$ Eqn. (11) is parabolic: a diffusion equation for pressure, whose coefficient $\eta = \mathbf{K} \lambda / (\phi \, c_t)$ is the (pressure, or hydraulic) diffusivity.

The transport equation (2) needs a corresponding term. The total velocity is no longer divergence-free: by eqn. (11) and reverting Darcy's law (7), $\nabla \cdot \mathbf{v} = q - \phi \, c_t \, \partial p / \partial t$, so the storage must be charged to the phases. This model does so in proportion to their saturation, $$\phi \frac{\partial s}{\partial t} + s \, \phi \, c_t \frac{\partial p}{\partial t} + \nabla \cdot (f(s)\, \mathbf{v}) = q_w \,, \tag{12}$$ which is what makes the water and oil equations sum to eqn. (11), so that e.g. depleting a fully water-saturated reservoir leaves $s = 1$, rather than conjuring oil out of the produced volume. (Deriving each phase equation individually would instead charge the water $s \, (c_r + c_w) \, \phi \, \partial p / \partial t$; the two coincide iff $c_w = c_o$, the difference being within the $O(c_t)$ fidelity anyway.) Ref minires.ResSim.storage_rate. Both new terms vanish for $c_t = 0$, recovering eqns. (1) and (2) exactly.

Consequences

The now parabolic pressure equation (11) is discretized here by backward Euler over the same $\Delta t$ as the saturation step, which adds $\phi \, c_t \, h^2 / \Delta t$ to the diagonal of the system (10), rendering it nonsingular without pinning. Thus the solution method survives. The sequential splitting remains applicable, and the pressure step is still one sparse linear solve, with no Newton iteration on $p$, and no PVT properties ($\rho$, $\mu$, $B$, $\phi$) to update with the pressure.

The absolute pressure level is meaningful, so an initial pressure must be given. However, the datum remains arbitrary. Eqn. (11) involves $p$ only through its derivatives, so shifting P0 (and any BHP targets) by a constant shifts the whole pressure trajectory by it, leaving saturations and rates untouched. The level is thus consequential (unlike for $c_t = 0$, it is propagated, not free) but only relative to the initial one. An absolute pressure would enter only through pressure-dependent properties -- precisely what the approximation drops.
Sources and sinks need not balance. The imbalance -- the voidage, production minus injection -- is supplied by expansion, permitting primary depletion by a lone producer. Summing the rows of the system (ref tests/test_compressible.py) yields $c_t \, \Delta \bar{p} = V_{\text{voidage}} / V_{\text{pore}}$, so the fidelity requirement, $c_t \, \Delta p \ll 1$, is a matter of the voidage asked of the fluids, not of choosing ct small. Linearity again: the mean pressure declines in proportion to the cumulative voidage, whatever its distribution in space or time -- the straight line of the material-balance plot (ref minires.ResSim.ct), by which pore volume is estimated.
Pressure is transient rather than instantaneous: $\sqrt{\eta t}$ is the radius of investigation, how far a well has "felt" after time $t$. Flow is called transient while that radius is still growing, and pseudo-steady state (or boundary-dominated) once it has reached the whole of the drainage volume, whereafter the pressure declines uniformly. Well testing is the inverse problem of inferring $\mathbf{K}$ and the skin (ref minires.wells.peaceman_WI) from a measured transient, typically during the build-up after shutting a well in -- as examples.buildup does. By the linear approximations, superposition holds, in space and in time: a shut-in is a flowing well plus an equal and opposite one started at the shut-in, and the pressure anywhere is the sum of the wells' individual transients. This is what makes well testing an analytical inference method -- the line-source solution, the Horner plot, and the semilog-derivative plateau that examples.buildup reads $\mathbf{K}$ off, are all solutions of the linear diffusion equation.

View Source

 1""".. include:: README.md"""
 2
 3from importlib.metadata import version as _version
 4
 5from minires.core import ResSim
 6from minires.grid import Fluxes, Grid2D
 7from minires.fluids import Fluid
 8from minires.wells import Wells, aquifer_WI, peaceman_WI, well_path
 9
10# Also pdoc's table of contents: `ResSim` is documented on the package page (beside the
11# README), and the listed submodules on their own pages. `core` is deliberately absent,
12# lest `ResSim` be documented twice; so are the other re-exports, which are documented
13# in their home modules.
14__all__ = ["ResSim", "grid", "fluids", "wells", "plotting", "tlm", "geostat"]
15
16__version__ = _version("minires")  # single source: `version` in pyproject.toml

@dataclass(repr=False)

class ResSim(struct_tools.AlignedRepr, minires.grid.Grid2D, minires.plotting.Plot2D): View Source

 20@dataclass(repr=False)
 21class ResSim(AlignedRepr, Grid2D, Plot2D):
 22    """Reservoir simulator class.
 23
 24    Implemented with OOP (instead of passing around dicts) to facilitate
 25    bookkeeping of ensemble forecasting
 26    (where parameter values of one instance should not influence another)
 27
 28    Example:
 29    >>> model = ResSim(Lx=1, Ly=1, Nx=64, Ny=64, wells=[
 30    ...     dict(xy=[0, .32], rate=+1),   # injector
 31    ...     dict(xy=[1, 1], rate=-1),     # producer
 32    ... ])
 33    >>> water_sat0 = np.zeros(model.Nxy)
 34    >>> dt = .35
 35    >>> nSteps = 2
 36    >>> S, P = model.sim(dt, nSteps, water_sat0, pbar=False)
 37
 38    This produces the following values (used for automatic testing):
 39    >>> S[-1, [100, 1300, 2900]]
 40    array([0.9429345 , 0.91358172, 0.71554613])
 41    """
 42
 43    # Prefer __setattr__ approach (over @property get/set-ers)
 44    # because @property requires the _private pattern,
 45    # which is pretty ugly with dataclasses.
 46    def __setattr__(self, key: str, val: Any) -> None:
 47        # Defaults that the dataclass cannot express, depending as they do on the grid
 48        if val is None:
 49            if key == "K":
 50                val = np.ones((2, *self.shape))
 51            elif key == "por":
 52                val = np.ones(self.shape)
 53            elif key == "active":
 54                val = np.ones(self.shape, bool)
 55        # Permeabilities
 56        if key == "K" and val is not None:
 57            if np.isscalar(val):
 58                val = np.full(self.shape, val, dtype=float)
 59            if val.size == self.size:
 60                val = np.stack([val, val])  # both components
 61            val = val.reshape((2, *self.shape))
 62        # Active cells: the first is where `TPFA` pins the incompressible pressure
 63        if key == "active" and val is not None:
 64            val = np.asarray(val, bool).reshape(self.shape)
 65            assert val.any(), "No active cells."
 66            self._pin = int(np.argmax(val))
 67        # Fluid: a dict (or `None`) builds the `Fluid` it parameterizes
 68        if key == "fluid":
 69            if val is None:
 70                val = Fluid()
 71            elif isinstance(val, dict):
 72                val = Fluid(**val)
 73        # Wells -- records (or `None`) get assembled into a `Wells`, which then
 74        # gets bound, whereupon it snaps its completions onto this grid.
 75        # NB: the wells' own normalization is `minires.wells.Wells.__setattr__`.
 76        if key == "wells":
 77            if not isinstance(val, Wells):
 78                val = Wells.from_records(self, val)
 79            val._bind(self)
 80        # Set
 81        super().__setattr__(key, val)
 82
 83    def __getstate__(self) -> dict:
 84        # The cached factorization (ref `cached_precond`) is a `SuperLU`, which
 85        # cannot be pickled -- and is a mere cache, so `deepcopy` and
 86        # multiprocessing (as HistoryMatching does) simply leave it behind.
 87        state = self.__dict__.copy()
 88        state.pop("_pLU", None)
 89        return state
 90
 91    name: str = "Unnamed"
 92    """Description."""
 93
 94    cdarcy: float = 1.0
 95    """Unit conversion factor for Darcy's law, $C$ -- ECLIPSE's `CDARCY`.
 96
 97    If you want to change unit system you not only need to manually convert
 98    the dimensional input quantities to the new units, but also change $C$ according to
 99    $$ C = \\frac{u_k \\, u_p \\, u_t}{u_μ \\, u_L^2} \\,, $$
100    (with $u_k$ the SI magnitude of the unit chosen for $k$).
101    Any *coherent* system gives `1`: choose base units for length, time and mass,
102    derive $u_p = M/(L T^2)$, $u_μ = M/(L T)$ and $u_k = L^2$ from them.
103
104    | System | $u_L$ | $u_t$ | $u_p$ | $u_k$ | $u_μ$ | rate | $C$ |
105    |---|---|---|---|---|---|---|---|
106    | SI | m | s | Pa | m² | Pa·s | m²/s | `1` |
107    | CGS | cm | s | barye | cm² | poise | cm²/s | `1` |
108    | MTS | m | s | pièze | m² | pz·s | m²/s | `1` |
109    | mm-ms-g | mm | ms | MPa | mm² | kPa·s | mm²/ms | `1` |
110    | Darcy's own | cm | s | atm | darcy | cP | cm²/s | `1` |
111    | metric | m | day | bar | mD | cP | m²/day | `0.008527` |
112    | field-like | ft | day | psi | mD | cP | ft²/day | `0.006328` |
113    | lab | cm | hour | atm | mD | cP | cm²/hour | `3.6` |
114
115    .. note:: The rate unit is forced to $u_L^2/u_t$ -- an areal rate.
116
117        A well rate of `20` for a 25 m thick reservoir means 500 m³/day.
118
119    .. note:: $C$ enters at exactly 2 sites, both of them Darcy's law.
120
121        The transmissibilities of `TPFA` and the well index of
122        `minires.wells.peaceman_WI`. Everything else is derivative, and
123        already consistent.
124    """
125
126    fluid: Any = None
127    """The two-phase fluid: a `minires.fluids.Fluid`, holding the viscosities
128    and the Corey relative permeability parameters, and computing the mobilities
129    and fractional flow from them.
130
131    Assigning a `dict` builds one (`ResSim(fluid=dict(vo=5, swc=.2))`), and `None`
132    the default (unit viscosities, the reference paper's quadratic curves). Its
133    fields stay writable (`model.fluid.vo = 5`). A `Fluid` subclass (e.g. with
134    tabulated curves) may be assigned instead.
135    """
136    ct: float = 0.0
137    """Total (rock + fluids) compressibility, $c_t$, as a single constant.
138
139    The default, `0`, yields the incompressible model, whose pressure eqn. is
140    elliptic: pressure is defined only up to an additive constant, and the
141    sources/sinks must balance. Setting `ct > 0` yields the *slightly
142    compressible* model: the pressure eqn. gains the accumulation term
143    $ φ \\, c_t \\, ∂p/∂t $ (discretized by backward Euler over the same `dt` as
144    the saturation step, which is what makes `P0` of `sim` consequential), and
145    the transport eqn. the matching storage term, charged to the phases in
146    proportion to their saturation (ref `storage_rate`). Injection and
147    production then need not balance, enabling e.g. primary depletion.
148
149    Derivation and fidelity ($ c_t \\, Δp \\ll 1 $, which the voidage sets, not
150    `ct`): ref the "Compressibility" section of the docs.
151
152    The **drive mechanism** is whatever supplies the energy that pushes the
153    hydrocarbons to the well. *Fluid and rock expansion* (a.k.a. **depletion
154    drive**), which is what $ c_t > 0 $ enables here in the absence of injection,
155    is the weakest, recovering only a few percent, because $c_t$ is so small.
156    Stronger ones are *solution gas drive*, *gas cap drive*, *water drive*
157    (aquifers, ref `minires.wells.aquifer_WI`), *gravity drainage*, and
158    *compaction drive* (which manifests as seabed subsidence).
159    Recovery is staged: **primary** production runs on the native drive.
160    Upon drilling, unless valves are in place, the pressure of the initial
161    *equilibrium* will cause a *blow out*; a new equilibrium is usually attained
162    when 20% of the hydrocarbons have been produced, which marks the end of the
163    primary production.
164    **Secondary** recovery adds *pressure support* by injecting water or gas
165    (**waterflooding** being the case simulated here);
166    **tertiary**, or **EOR** (enhanced oil recovery), alters the flow physics itself,
167    e.g. by polymer, surfactant, or CO₂ injection.
168    The **voidage replacement ratio** is the injected reservoir volume divided by the
169    produced one; $\\mathrm{VRR} = 1$ is exactly the balance, $\\sum q = 0$,
170    that the incompressible model is obliged to impose.
171    The zero-dimensional (single tank) accounting of all of the above,
172    used to estimate reserves without a grid, is called **material balance**.
173    """
174    cached_precond: bool = True
175    """Solve the pressure system iteratively, preconditioned by a cached factorization.
176
177    The alternative (`False`) is a fresh sparse direct factorization each
178    time step. But the system changes slowly -- only through the mobility
179    $λ(s)$, i.e. where the front has moved -- so the factorization of an
180    *earlier* step remains an excellent preconditioner: with it, conjugate
181    gradients converges in 1--2 iterations if the saturation is (nearly)
182    static, and in ~12 behind a moving front (each the cost of one
183    back-substitution, i.e. 1/30 of a factorization). The factorization is
184    refreshed only when convergence fails (as it does when the mobility has
185    drifted too far, e.g. after many steps at a strong viscosity contrast),
186    so the result is exact to the solver tolerance (`1e-10`) either way. It
187    requires the system to be SPD, which the TPFA system is (ref the "How to
188    solve" section of the docs).
189
190    The cache is per instance (`_pLU`) and is dropped on pickling and
191    `deepcopy`, a `SuperLU` not being picklable; a copy simply refactorizes
192    on its next step. The measurements (2--6x on the well-test and depletion
193    examples, 15--40% on the waterfloods), and what was tried besides, are
194    recorded in `tests/test_precond.py`.
195    """
196
197    # NB: the array attributes are typed `Any` since `__setattr__` normalizes
198    # whatever array-like (nested lists, scalars) is assigned to them.
199    K: Any = None
200    """Permeabilities (in x and y directions). Array of shape `(2, Nx, Ny)`).
201
202    **Permeability**, denoted by tensor $\\mathbf{K}$, quantifies transmissibility.
203    Usually SPD -- here diagonal, $\\mathrm{diag}(K_x, K_y)$ per cell, which is
204    what the two components of the first axis are -- and correlated with `por`.
205    Among the reservoir rocks,
206    *sandstone* usually have large, well-connected pores,
207    and high permeability, *shale* is nearly impermeable,
208    like cap rock and bed rock.
209    Permeability is measured in Darcy ($≈ 10^{-12} m^2$), ref `cdarcy`.
210    A medium is called *isotropic* if $\\mathbf{K}$ is scalar.
211    Its spatial structure goes by yet more lingo: *facies* (bodies of rock of
212    like character, e.g. the *channels* of a fluvial deposit), and the *fissures*
213    and *fractures* whose conductivity dwarfs that of the rock around them.
214    Being the most uncertain of the reservoir's properties (ref `minires.fluids`),
215    it is the parameter that history matching estimates -- and the one that
216    `minires.tlm` differentiates (as $\\log \\mathbf{K}$) with respect to.
217    """
218    por: Any = None
219    """Porosity; Array of shape `(Nx, Ny)`).
220
221    **Porosity**, $φ$, is the *void volume fraction*.
222    Depends on pressure, because rock is compressible.
223    *Compressibility* is the porosity's (relative) gradient wrt. pressure.
224    Usually neglected, as it is here in the pore volume (ref `pore_volume`),
225    so that $φ$ is a constant, but spatial, field; the rock's compressibility
226    enters instead through `ct`, which lumps it with the fluids'.
227    """
228    active: Any = None
229    """Mask of the active cells, `(Nx, Ny)`, boolean. Default: all `True`.
230
231    Setting some cells inactive carves an irregular reservoir out of the
232    rectangular grid -- an outline, holes, or a (sealing) fault: a line of
233    inactive cells, which blocks flow so long as it is unbroken (a diagonal
234    staircase suffices, fluxes passing only through faces). Inactive cells
235    take no part in the physics:
236
237    - The faces to them carry zero transmissibility, hence zero flux, so the
238      reservoir is closed along their perimeter as it is along the boundary.
239    - They have no equations. The pressure system carries an identity row for
240      each (keeping it well conditioned, unlike a tiny permeability would), and
241      their pore volume is taken as infinite (ref `pore_volume`; so they never
242      bind the CFL, as a tiny porosity would). So their state is simply carried
243      through `sim` unchanged from `S0` and `P0`. Plots mask them out.
244    - A well may not be completed in one (`_validate` checks). Should the
245      active cells form several disconnected regions, each is a reservoir of
246      its own, and only the first is pinned (ref `TPFA`): if `ct == 0`, the
247      others must balance their own rates, or hold a BHP well -- which is not
248      checked as such, but the singular system it would otherwise make is
249      caught by the residual check of `_solve_pressure`.
250
251    The arrays (`K`, `por`, `S0`, ...) keep the full grid shape, and the flat
252    index (`xy2ind`) runs over all cells: the mask selects the active ones.
253    """
254
255    wells: Any = None
256    """The wells: a `Wells`, holding the flat, per-*completion* arrays --
257    positions, rates, pressures, well indices -- that the model runs on.
258
259    Assigning a list (or `dict`) of records -- one per well -- assembles one,
260    which is the convenient way to configure them; the record format is
261    documented in `minires.wells.Wells.from_records`. Assigning `None`
262    empties it. A `Wells` may also be given directly, in which case it is
263    *bound* to this model, whereupon its completions snap onto the grid.
264
265    The arrays remain writable throughout (`model.wells.rates = ...`), as an
266    ensemble or optimisation loop requires.
267    """
268
269    def pore_volume(self) -> np.ndarray:
270        """Pore volume (per unit thickness) of each cell, `h2 * por`. Flat.
271
272        `inf` for inactive cells (ref `active`): dividing by it, the transport
273        schemes leave their saturation alone, and the CFL estimate ignores them.
274        """
275        pv = self.h2 * self.por.ravel()
276        return np.where(self.active.ravel(), pv, np.inf)
277
278    nComp = property(lambda self: self.wells.nComp)
279    """Num. of well *completions*, i.e. the rows of every array the model
280    indexes by them -- which is what it actually solves for, the equations
281    being assembled per completion.
282    Forwarded from `minires.wells.Wells.nComp`.
283    """
284
285    def assemble_wells(
286        self, S: np.ndarray | None, P: np.ndarray | None, k: int
287    ) -> None:
288        """Set up (for time `k`) the wells' contributions to the equations.
289
290        The controls are those of `well_controls`, to which `S` and `P` (the
291        state at the *start* of the step) are simply passed on.
292        Rate-controlled wells enter the source/sink *field*, `_Q`, directly.
293        BHP-controlled ones (ref `minires.wells.Wells.bhp`) cannot: their
294        rate is not yet known. They instead enter the pressure equations in
295        `TPFA`, after which `realize_bhp` folds the resulting rate into `_Q`.
296        """
297        ctrl = self.well_controls(S, P, k)
298        inds = self.xy2ind(*self.wells.xy.T)
299        rates, p_bh = ctrl["rates"], ctrl["bhp"]
300        is_bhp = np.isfinite(p_bh)
301        assert np.isfinite(rates[~is_bhp]).all(), (
302            "A rate-controlled well has a non-finite rate. Give it a number"
303            " (`0` shuts it in), or put it on BHP control; ref `Wells.rates`."
304        )
305
306        # The well model's constant of proportionality, WI * λ_t.
307        # NB: `nan` marks the rate-controlled wells, throughout.
308        WI_lam = np.full(self.nComp, np.nan)
309        if is_bhp.any():
310            WI = self.wells.WI
311            assert WI is not None and np.isfinite(WI[is_bhp]).all(), (
312                "BHP control requires (finite) `Wells.WI`."
313            )
314            assert S is not None, "BHP control requires `S` (for λ_t)."
315            Mw, Mo = self.fluid.RelPerm(S)
316            WI_lam[is_bhp] = WI[is_bhp] * (Mw + Mo)[inds[is_bhp]]
317
318        # Translate well conditions for cells.
319        # NB: Dont use `Q[inds] += ...` since `inds` may contain dupes.
320        self._Q, bhp_diag, bhp_rhs = np.zeros((3, self.Nxy))
321        np.add.at(self._Q, inds[~is_bhp], rates[~is_bhp])
322        np.add.at(bhp_diag, inds[is_bhp], WI_lam[is_bhp])
323        np.add.at(bhp_rhs, inds[is_bhp], (WI_lam * p_bh)[is_bhp])
324        rates[is_bhp] = np.nan  # only `realize_bhp` knows these
325        self._wells_now: dict[str, np.ndarray] = dict(
326            inds=inds, rates=rates, p_bh=p_bh,
327            WI_lam=WI_lam, bhp_diag=bhp_diag, bhp_rhs=bhp_rhs,
328        )  # fmt: skip
329
330    def realize_bhp(self, P: np.ndarray) -> None:
331        """Compute rates for BHP wells. Enter into `_Q` and `_wells_now["rates"]`.
332
333        The rate, $ WI λ_t (p_\\mathrm{bh} - p_\\mathrm{cell}) $, is signed by
334        nature (ref the `minires.wells.Wells.bhp` warning).
335
336        By construction of the linear system of `TPFA`, this leaves `_Q` equal
337        to the *total* well flux, which is what keeps `storage_rate` -- and
338        hence the transport step -- consistent with the pressure solution.
339        """
340        wls = self._wells_now
341        WI_lam = wls["WI_lam"]  # `nan` marks the rate-controlled wells
342        # Insert in cell source/sink field
343        self._Q = self._Q + wls["bhp_rhs"] - wls["bhp_diag"] * P
344        # Insert in per-well rates
345        is_bhp = np.isfinite(WI_lam)
346        wls["rates"][is_bhp] = (WI_lam * (wls["p_bh"] - P[wls["inds"]]))[is_bhp]
347
348    def _record_actual_well_operation(
349        self, S: np.ndarray, P: np.ndarray, k: int
350    ) -> None:
351        """Record `actual_rates`/`actual_bhp`. Warn about flow direction flip."""
352        wls = self._wells_now
353        if k:
354            is_bhp = np.isfinite(wls["WI_lam"])
355            flipped = is_bhp & (wls["rates"] * self.wells.actual_rates[:, k - 1] < 0)
356            if flipped.any():
357                warnings.warn(
358                    f"BHP-controlled well(s) {np.flatnonzero(flipped).tolist()}"
359                    f" reversed flow direction at step {k}"
360                    " (an inflow injects water); ref `Wells.bhp`.",
361                    stacklevel=2,
362                )
363        self.wells.actual_rates[:, k] = wls["rates"]
364        self.wells.actual_bhp[:, k] = self.bhp(S, P, wls["rates"])
365
366    def well_controls(self, S: np.ndarray | None, P: np.ndarray | None, k: int) -> dict:
367        """Compute the wells' controls for time `k`: `dict(rates=..., bhp=...)`.
368
369        Each is a `(nComp,)` array, read off the specifications --
370        `minires.wells.Wells.rates`, `minires.wells.Wells.bhp` -- which
371        are *open-loop*: fixed before the simulation begins. Overriding
372        (patching/subclassing) this method is therefore how to do *feedback*
373        control, the controls being free to depend on the state at the *start*
374        of the step: the saturation `S` and the pressure `P`.
375        The returned arrays are copies, so they may be modified in place.
376
377        Most feedback concerns the rates alone -- e.g. shutting the wells upon
378        water breakthrough at the producer:
379
380        >>> class Shutter(ResSim):
381        ...     def well_controls(self, S, P, k):
382        ...         ctrl = super().well_controls(S, P, k)
383        ...         if S is not None and S[self.xy2ind(1, 1)] > .5:
384        ...             ctrl["rates"][:] = 0    # NB: all of them! See warning
385        ...         return ctrl
386        >>> model = Shutter(Lx=1, Ly=1, Nx=16, Ny=16,
387        ...                 wells=Wells(xy=[[0, 0], [1, 1]], rates=[[1], [-1]]))
388        >>> SS, PP = model.sim(.05, 20, model.fluid.swc*np.ones(model.Nxy), pbar=False)
389        >>> int((model.wells.actual_rates[1] == 0).argmax())  # step of breakthrough
390        16
391
392        But the `bhp` is here too, and with it each well's *control mode*
393        (`nan` => rate-controlled, ref `minires.wells.Wells.bhp`) -- which
394        is what an approximate mode *switch* requires. For example, rate
395        control with a BHP limit -- the industrial default -- wherein a
396        producer holds its rate only for as long as that does not draw it
397        below some `p_min`:
398
399        >>> class Limited(ResSim):
400        ...     p_min = .5
401        ...     def well_controls(self, S, P, k):
402        ...         ctrl = super().well_controls(S, P, k)
403        ...         if P is None:
404        ...             return ctrl                    # nothing to switch on
405        ...         p_bh = self.bhp(S, P, ctrl["rates"])
406        ...         switch = p_bh < self.p_min         # the rate is unsustainable
407        ...         ctrl["bhp"] = np.where(switch, self.p_min, np.nan)
408        ...         return ctrl
409        >>> model = Limited(Lx=1, Ly=1, Nx=16, Ny=16, ct=.1,
410        ...                 wells=Wells(xy=[[.5, .5]], rates=[[-.25]]))
411        >>> from minires import peaceman_WI
412        >>> model.wells.WI = peaceman_WI(model, model.wells.xy, rw=1e-3)
413        >>> SS, PP = model.sim(.02, 25, np.zeros(model.Nxy),
414        ...                    P0=np.ones(model.Nxy), pbar=False)
415
416        The well delivers its target rate until the limit binds, and declines
417        thereafter -- at constant $ p_\\mathrm{bh} $, exponentially so
418        (ref `examples/well_control.py`, which plots all three modes):
419
420        >>> (-model.wells.actual_rates[0, [0, 5, 6, -1]]).round(3)
421        array([0.25 , 0.25 , 0.182, 0.005])
422
423        .. warning:: With `ct == 0` the rates must still sum to 0 at every step.
424
425            Ref `minires.wells.Wells.rates`. So shutting one well requires
426            matching it on the other side -- as above.
427
428        .. note:: The mode switch lags the solve by one step.
429
430            It is decided from the previous step's pressure, whereas the well
431            model itself is solved *simultaneously* with the new one. So the
432            limit is breached for the one step in which it comes to bind.
433            Shorten `dt` to refine.
434
435        .. note:: `S` and `P` may be `None`, so an override should tolerate that.
436
437            They are `None` if the caller has none to offer -- as when
438            `assemble_wells` is used merely to set up a plot.
439
440        .. note:: Setting both controls for a well is not an error, just pointless.
441
442            `assemble_wells` discards the rate of a BHP-controlled well -- it
443            is `realize_bhp` that fills it in.
444        """
445        return dict(
446            rates=self.wells.at_time("rates", 0.0, k),
447            bhp=self.wells.at_time("bhp", np.nan, k),
448        )
449
450    def bhp(self, S: np.ndarray, P: np.ndarray, rates: np.ndarray) -> np.ndarray:
451        """Bottom-hole pressures implied by the (signed) `rates`, via the well indices.
452
453        I.e. the well model of `minires.wells.Wells.WI`, solved for
454        $ p_\\mathrm{bh} $:
455        the rate's sign puts an injector above, a producer below, its cell
456        pressure. `nan` wherever the well index is unset.
457
458        `S` and `P` (both flat) should be the saturation and the pressure of the
459        *same* `pressure_step`, i.e. `SS[k]` and `PP[k+1]` of `sim` -- which is
460        what `actual_bhp` records, so prefer reading that.
461
462        .. warning:: $ λ_t $ is that of the well's *cell*.
463
464            So an injector's injectivity is governed by the mobility of
465            whatever the cell currently holds, rather than by that of the
466            injectant.
467        """
468        if self.wells.WI is None:
469            return np.full(self.nComp, np.nan)
470        Mw, Mo = self.fluid.RelPerm(S)
471        ii = self.xy2ind(*self.wells.xy.T)
472        return P[ii] + rates / (self.wells.WI * (Mw + Mo)[ii])
473
474    # Pres() -- listing 5
475    def pressure_step(
476        self,
477        S: np.ndarray,
478        P: np.ndarray | None = None,
479        dt: float | None = None,
480    ) -> tuple[np.ndarray, Fluxes]:
481        """Compute permeabilities then solve Darcy's equation. Returns `[P, V]`.
482
483        `P` (flat, like `S`) is the *previous* step's pressure: used (and
484        required) only if `ct > 0`, along with `dt`. The new one replaces it.
485        """
486        # Compute K*λ(S)
487        Mw, Mo = self.fluid.RelPerm(S)
488        Mt = Mw + Mo
489        Mt = Mt.reshape(self.shape)
490        KM = Mt * self.K
491        # Compute pressure and extract fluxes
492        [P, V] = self.TPFA(KM, P, dt)
493        return P, V
494
495    def _spdiags(self, data: Any, diags: Any) -> sparse.dia_matrix:
496        """`sparse.spdiags` of the `(Nxy, Nxy)` matrix -- tolerating coincident offsets.
497
498        `TPFA` and `upwind_diff` place the x- and y-neighbours at offsets
499        $ ±N_y $ and $ ±1 $, which coincide when $ N_y = 1 $ (a 1D row of
500        cells). `scipy` rejects duplicate offsets, so they are summed here --
501        the y-diagonals then holding only zeros, there being no y-faces.
502        """
503        diags = np.atleast_1d(diags)
504        if len(diags) > len(set(diags)):
505            uniq, inv = np.unique(diags, return_inverse=True)
506            summed = np.zeros((len(uniq), self.Nxy))
507            np.add.at(summed, inv, np.atleast_2d(data))
508            data, diags = summed, uniq
509        return sparse.spdiags(data, diags, self.Nxy, self.Nxy)
510
511    # TPFA() -- Listing 1
512    def TPFA(
513        self,
514        K: np.ndarray,
515        P: np.ndarray | None = None,
516        dt: float | None = None,
517    ) -> tuple[np.ndarray, Fluxes]:
518        """Two-point flux-approximation (TPFA) of Darcy: $ -∇(K ∇u) = q $
519
520        i.e. steady-state diffusion w/ nonlinear coefficient, $K$,
521        if `ct == 0`. Otherwise (slightly compressible model) solve
522        the backward-Euler step of $ φ c_t ∂u/∂t - ∇(K ∇u) = q $,
523        which requires the previous pressure, `P`, and `dt`.
524
525        After solving for pressure `P`, extract the fluxes `V`
526        by finite differences.
527        """
528        # Compute transmissibilities by harmonic averaging.
529        C = self.cdarcy
530        L = 1 / K
531        TX = np.zeros((self.Nx + 1, self.Ny))
532        TY = np.zeros((self.Nx, self.Ny + 1))
533        TX[1:-1, :] = C * 2 * self.hy / self.hx / (L[0, :-1, :] + L[0, 1:, :])
534        TY[:, 1:-1] = C * 2 * self.hx / self.hy / (L[1, :, :-1] + L[1, :, 1:])
535        # No flow across the faces of inactive cells (ref `active`)
536        act = self.active
537        TX[1:-1, :] *= act[:-1, :] & act[1:, :]
538        TY[:, 1:-1] *= act[:, :-1] & act[:, 1:]
539
540        # Assemble TPFA discretization matrix.
541        x1 = TX[:-1, :].ravel()
542        x2 = TX[1:, :].ravel()
543        y1 = TY[:, :-1].ravel()
544        y2 = TY[:, 1:].ravel()
545
546        # Setup linear system
547        DiagVecs = [-x2, -y2, y1 + y2 + x1 + x2, -y1, -x1]
548        DiagIndx = [-self.Ny, -1, 0, 1, self.Ny]
549        q = self._Q
550        if self.ct > 0:
551            # Accumulation term (φ ct h²/dt) of backward Euler.
552            # Renders the system nonsingular (unlike the pure-Neumann problem).
553            accum = self.por.ravel() * self.ct * self.h2 / dt
554            DiagVecs[2] = DiagVecs[2] + accum
555            q = q + accum * P
556        elif not self._wells_now["bhp_diag"].any():
557            # Pin the (o/w pure-Neumann & singular) problem, at the 1st active cell
558            i = self._pin
559            DiagVecs[2][i] += np.sum(self.K.reshape(2, -1)[:, i])  # ref article p. 13
560        # Well model of the BHP-controlled wells
561        DiagVecs[2] = DiagVecs[2] + self._wells_now["bhp_diag"]
562        q = q + self._wells_now["bhp_rhs"]
563        # Inactive cells have no equation. Identity rows (their faces being
564        # already closed) keep them at their previous pressure, and the system SPD.
565        act = act.ravel()
566        DiagVecs[2] = np.where(act, DiagVecs[2], 1.0)
567        q = np.where(act, q, 0.0 if P is None else P)
568
569        # Solve; compute A\q to update P
570        A = self._spdiags(DiagVecs, DiagIndx)
571        P = self._solve_pressure(A.tocsr(), q, P)
572        # P = np.linalg.solve(A.A, q) # direct dense solver
573        # Could also try scipy.linalg.solveh_banded which, according to
574        # https://scicomp.stackexchange.com/a/30074 uses the Thomas algorithm,
575        # as recommended by Aziz and Settari ("Petro. Res. simulation").
576        # NB: stackexchange also mentions that solve_banded does not work well
577        # when the band offsets large, i.e. higher-dimensional problems.
578
579        # Extract fluxes, via a grid-shaped view of the (flat) pressure.
580        P2d = P.reshape(self.shape)
581        V = Fluxes(
582            x=np.zeros((self.Nx + 1, self.Ny)),
583            y=np.zeros((self.Nx, self.Ny + 1)),
584        )
585        V.x[1:-1, :] = (P2d[:-1, :] - P2d[1:, :]) * TX[1:-1, :]
586        V.y[:, 1:-1] = (P2d[:, :-1] - P2d[:, 1:]) * TY[:, 1:-1]
587        return P, V
588
589    def _solve_pressure(
590        self, A: sparse.csr_matrix, q: np.ndarray, P0: np.ndarray | None
591    ) -> np.ndarray:
592        """Solve the (SPD) pressure system `A P = q`, ref `cached_precond`.
593
594        `P0` is the initial guess of the iterative solver (the previous pressure).
595        The cached factorization is `_pLU`; being a mere preconditioner, it is
596        safe to hold stale (across `sim` calls, or a change of `K`), since the
597        iteration converges to the solution of the *current* `A` regardless,
598        and refactorizes when it does not converge.
599
600        Either way, the factorization orders the columns by MMD on `A + A'`
601        (the ordering for a symmetric matrix), which halves the fill, hence
602        the cost, of the COLAMD that `spsolve` would use: 1.5x on the direct
603        solve, for free.
604        """
605        if self.cached_precond:
606            LU = getattr(self, "_pLU", None)
607            if LU is not None and LU.shape == A.shape:
608                LinOp: Any = LinearOperator  # (ty cannot see its factory `__new__`)
609                M = LinOp(A.shape, matvec=LU.solve, dtype=A.dtype)
610                P, info = cg(A, q, x0=P0, M=M, rtol=1e-10, maxiter=30)
611                if info == 0:
612                    return P
613        # Direct solve: first call, non-convergence, or `not cached_precond`.
614        LU = splu(A.tocsc(), permc_spec="MMD_AT_PLUS_A")
615        if self.cached_precond:
616            self._pLU = LU
617        P = LU.solve(q)
618        # A singular system (an incompressible region of `active` cells whose
619        # rates do not balance) does not make `splu` raise, but leaves a
620        # residual of O(1) that no `P` can remove, whereas O(1e-14) is normal.
621        # NB: explicit connectivity checks of `active` (graph search; a cursory
622        # pressure step) were tried and found laborious, while also refusing
623        # the valid disconnected configs (each region balanced, or well-less).
624        # The tolerance is loose because the achievable residual degrades with
625        # the conditioning, i.e. with the permeability contrast: an ensemble of
626        # log-normal fields (HistoryMatching's, contrast up to 1e10) reaches
627        # 6e-8, so 1e-8 cried wolf on it -- while the singular case it is
628        # actually after stays 5 orders of magnitude above 1e-6.
629        assert np.linalg.norm(A @ P - q) <= 1e-6 * np.linalg.norm(q), (
630            "The pressure solve failed. Is a disconnected region of `active`"
631            " cells left without balanced rates (ref `active`)?"
632        )
633        return P
634
635    # GenA() -- listing 7
636    def upwind_diff(self, V: Fluxes) -> sparse.dia_matrix:
637        """Upwind finite-volume scheme."""
638        fp = self._Q.clip(max=0)  # production
639        # Flow fluxes, separated into direction (x-y) and sign
640        x1 = V.x.clip(max=0)[:-1, :].ravel()
641        y1 = V.y.clip(max=0)[:, :-1].ravel()
642        x2 = V.x.clip(min=0)[1:, :].ravel()
643        y2 = V.y.clip(min=0)[:, 1:].ravel()
644        DiagVecs = [x2, y2, fp + y1 - y2 + x1 - x2, -y1, -x1]
645        DiagIndx = [-self.Ny, -1, 0, 1, self.Ny]
646        A = self._spdiags(DiagVecs, DiagIndx)
647        return A
648
649    def storage_rate(self, V: Fluxes) -> np.ndarray:
650        """The volume rate, per cell, that goes into storage: $ q - ∇ ⋅ V $.
651
652        For the incompressible model this is `0`: the fluxes balance the wells
653        exactly, cell by cell. With `ct > 0` it is (by construction of the
654        linear system of `TPFA`) the accumulation term of the backward-Euler
655        step, $ φ \\, c_t \\, h^2 \\, (p^{n+1} - p^n) / Δt $, which the saturation
656        steps charge to the phases (ref `ct`).
657
658        Computing it from `V` (rather than from $p^{n+1} - p^n$) means it is
659        *exactly* the imbalance seen by the transport scheme, whose `upwind_diff`
660        is assembled from the same fluxes.
661        """
662        if self.ct == 0:
663            return np.zeros(self.Nxy)
664        divV = (V.x[1:, :] - V.x[:-1, :]) + (V.y[:, 1:] - V.y[:, :-1])
665        return self._Q - divV.ravel()
666
667    # Extracted from Upstream()
668    def estimate_1CFL(self, pv: np.ndarray, V: Fluxes, fi: np.ndarray) -> float:
669        """Estimate 1/CFL for use with `saturation_step_upwind`."""
670        # In-/Out-flux x-/y- faces
671        XP = V.x.clip(min=0)
672        XN = V.x.clip(max=0)
673        YP = V.y.clip(min=0)
674        YN = V.y.clip(max=0)
675        Vi = XP[:-1, :] + YP[:, :-1] - XN[1:, :] - YN[:, 1:]
676
677        flx = max((Vi.ravel() + fi) / pv)  # estimate of influx
678        # NB: `storage_rate` is not counted here. In practice it is a small
679        # fraction of the fluxes that are (under 20% even at `ct = 10`),
680        # so the safety factor below covers it.
681        # The characteristic speed is f_w'(s) times the velocity: bound it by the
682        # maximal slope of f_w, sampled at 1001 points of the mobile range (which
683        # undershoots the true maximum by O(1e-6), nothing against the safety
684        # factor of 1.5). For the default curves at equal viscosities this is the
685        # paper's 3 / (1 - swc - sor) (Listing 8), which steeper curves (unequal
686        # viscosities, higher Corey exponents) exceed.
687        ss = np.linspace(self.fluid.swc, 1 - self.fluid.sor, 1001)
688        dfw_max = self.fluid.dfractional_flow(ss).max()
689        cfl = 1.5 * dfw_max * flx
690        # NB: the ceiling is nudged down by a relative epsilon, so that a `dt`
691        # sitting *on* an integer multiple of the CFL limit (as the examples'
692        # round numbers tend to) does not gain a whole extra sub-step from the
693        # last bits of the linear solve -- which differ across platforms and
694        # library versions, and would make the results irreproducible.
695        # The safety factor covers the shaving, the sampling's undershoot, and
696        # the storage rate.
697        return cfl * (1 - 1e-9)
698
699    # Upstream() -- listing 8
700    def saturation_step_upwind(self, S: np.ndarray, V: Fluxes, dt: float) -> np.ndarray:
701        """Explicit upwind FV discretisation of conserv. of mass (water sat.)."""
702        # fmt: off
703        A  = self.upwind_diff(V)                 # FV discretized transport operator
704        pv = self.pore_volume()                  # Pore volume (per thickness)
705        fi = self._Q.clip(min=0)                 # Well inflow
706        st = self.storage_rate(V)                # Storage (0 if incompressible)
707
708        # Compute sub/local dt
709        cfl1 = self.estimate_1CFL(pv, V, fi)
710        nT = int(np.ceil(dt * cfl1))
711        nT = max(1, nT)
712
713        # Scale A
714        dtx = dt / nT / pv                       # timestep / pore volume
715        B   = self._spdiags(dtx, 0) @ A          # A * dt/|Omega i|
716
717        for _ in range(nT):
718            fw = self.fluid.fractional_flow(S)      # fractional flow
719            S = S + (B@fw + (fi - S*st)*dtx)     # update saturation
720        # fmt: on
721        return S
722
723    # NewtRaph() -- listing 10
724    def saturation_step_implicit(
725        self,
726        S: np.ndarray,
727        V: Fluxes,
728        dt: float,
729        nNewtonMax: int = 10,
730        nTmax_log2: int = 10,
731    ) -> np.ndarray:
732        """Implicit FV discretisation of conserv. of mass (water sat.).
733
734        .. warning:: The Newton iteration can converge to a spurious root.
735
736            Far outside the $ c_t \\, Δp \\ll 1 $ regime (ref `ct`), it may
737            converge -- silently -- to a root of the residual outside $[0, 1]$:
738            the polynomial `minires.fluids.Fluid.RelPerm` extends smoothly
739            beyond the unit interval, and the sub-`dt` halving only triggers on
740            *non*-convergence.
741            The explicit scheme
742            (`saturation_step_upwind`), being monotone, stays within $[0, 1]$
743            even for extreme `ct`.
744
745        .. note:: This scheme rarely earns its keep.
746
747            It is usually both slower & less accurate than `saturation_step_upwind`.
748            Both schemes sub-divide `dt` internally, so it is not `dt` that decides
749            their cost but the stiffness of the grid -- and the well cells, being
750            normally the stiffest, hold the two requirements within a factor 3 of each
751            other (that being the safety margin of `ResSim.estimate_1CFL`), while an
752            implicit sub-step -- a sparse solve, or several -- 10 or 100x more.
753            It does pay off where the stiffest cell is *not* a well cell -- a tight
754            streak, a fracture, a locally refined region -- running 7 times faster at a
755            1000x pore-volume contrast. The branch `implicit-transport-scheme` says more.
756        """
757        # fmt: off
758        A  = self.upwind_diff(V)                 # FV discretized transport operator
759        pv = self.pore_volume()                  # Pore volume (per thickness)
760        fi = self._Q.clip(min=0)                 # Well inflow
761        st = self.storage_rate(V)                # Storage (0 if incompressible)
762
763        # For each iter, halve the sub/local dt
764        for nT_log2 in range(0, nTmax_log2):
765            nT = 2**nT_log2
766
767            # Scale A
768            dtx = dt / nT / pv                   # timestep / pore volume
769            B   = self._spdiags(dtx, 0) @ A      # A * dt/|Omega i|
770            C   = self._spdiags(dtx*st, 0)       # storage, likewise scaled
771
772            Sn = S
773            for _ in range(nT):
774                Sp = Sn
775                for _ in range(nNewtonMax):
776                    fw = self.fluid.fractional_flow(Sn)     # fract. flow
777                    df = self.fluid.dfractional_flow(Sn)    # its derivative
778                    dG = (sparse.eye(self.Nxy) + C                        # deriv of G
779                          - B @ self._spdiags(df, 0))
780                    G  = Sn - Sp - (B@fw + (fi - Sn*st)*dtx)  # G(s)
781                    dS = spsolve(dG, G)             # compute dS
782                    Sn = Sn - dS                    # update S
783
784                    if np.sqrt(sum(dS**2)) < 1e-3:
785                        # If converged: halt Newton iterations
786                        break
787                else:
788                    # If never converged: increase nT, restart time loop
789                    break
790            else:
791                # If completed all time steps, halt
792                break
793        else:
794            # Failed (even with max nT) to complete all time steps
795            print("Warning: did not converge")
796        # fmt: on
797
798        return Sn
799
800    def _validate(self):
801        # Catch some common issues before they become mysterious/insidious
802        act = self.active.ravel()
803        assert act[self._wells_now["inds"]].all(), (
804            "A well is completed in an inactive cell (ref `active`)."
805        )
806        if self.ct == 0 and not self._wells_now["bhp_diag"].any():
807            # No storage, no anchor ⇒ src/sinks must balance (ref `Wells.rates`)
808            SA = np.abs(self._Q).sum()
809            AS = abs(self._Q.sum())
810            assert AS <= 1e-10 * SA, "well rates do not sum to 0"
811        assert np.all((0 <= self.K) & np.isfinite(self.K))
812        assert np.all((0 <= self.por) & (self.por <= 1))
813
814    def time_stepper(self, dt: float, implicit: bool = False) -> Callable:
815        """Get ODE solver (integrator) for model.
816
817        Whatever time step `dt` is given, both schemes will use smaller steps internally.
818
819        - `explicit`: computes sub-`dt` based on CFL esitmate.
820        - `implicit`: reduces sub-`dt` until convergence is achieved.
821        """
822
823        def integrate(S, P, k):
824            self.assemble_wells(S, P, k)
825            self._validate()
826            [P, V] = self.pressure_step(S, P, dt)
827            self.realize_bhp(P)
828            self._record_actual_well_operation(S, P, k)
829            if implicit:
830                S = self.saturation_step_implicit(S, V, dt)
831            else:
832                S = self.saturation_step_upwind(S, V, dt)
833            return S, P
834
835        return integrate
836
837    def sim(
838        self,
839        dt: float,
840        nSteps: int,
841        S0: np.ndarray,
842        P0: np.ndarray | None = None,
843        pbar: bool = True,
844        leave: bool = True,
845        **kwargs,
846    ) -> tuple:
847        """Recursively (`nSteps` times) apply `time_stepper` with `dt`, from `S0`.
848
849        Returns the saturation and pressure trajectories, `(SS, PP)`.
850
851        .. note:: `SS[0] == S0` and `PP[0] == P0`, hence both have `len = nSteps + 1`.
852
853            `P0` defaults to zeros. It is only consequential if `ct > 0`.
854        """
855        step = self.time_stepper(dt, **kwargs)
856
857        # pbar
858        kk = np.arange(nSteps)
859        if pbar:
860            kk = tqdm(kk, "Simulation", leave=leave, mininterval=1e-2)
861
862        # Allocate
863        SS = np.zeros((nSteps + 1,) + S0.shape)
864        PP = np.zeros((nSteps + 1, self.Nxy))
865        self.wells.actual_rates = np.zeros((self.nComp, nSteps))
866        self.wells.actual_bhp = np.full((self.nComp, nSteps), np.nan)
867
868        # Init
869        SS[0] = S0
870        if P0 is not None:
871            PP[0] = P0
872
873        # Recurse
874        for k in kk:
875            SS[k + 1], PP[k + 1] = step(SS[k], PP[k], k)
876
877        return SS, PP

Reservoir simulator class.

Implemented with OOP (instead of passing around dicts) to facilitate bookkeeping of ensemble forecasting (where parameter values of one instance should not influence another)

Example:

>>> model = ResSim(Lx=1, Ly=1, Nx=64, Ny=64, wells=[
...     dict(xy=[0, .32], rate=+1),   # injector
...     dict(xy=[1, 1], rate=-1),     # producer
... ])
>>> water_sat0 = np.zeros(model.Nxy)
>>> dt = .35
>>> nSteps = 2
>>> S, P = model.sim(dt, nSteps, water_sat0, pbar=False)

This produces the following values (used for automatic testing):

>>> S[-1, [100, 1300, 2900]]
array([0.9429345 , 0.91358172, 0.71554613])

ResSim( Lx: float = 1.0, Ly: float = 1.0, Nx: int = 32, Ny: int = 32, name: str = 'Unnamed', cdarcy: float = 1.0, fluid: Any = None, ct: float = 0.0, cached_precond: bool = True, K: Any = None, por: Any = None, active: Any = None, wells: Any = None)

name: str = 'Unnamed'

Description.

cdarcy: float = 1.0

Unit conversion factor for Darcy's law, $C$ -- ECLIPSE's CDARCY.

If you want to change unit system you not only need to manually convert the dimensional input quantities to the new units, but also change $C$ according to $$ C = \frac{u_k \, u_p \, u_t}{u_μ \, u_L^2} \,, $$ (with $u_k$ the SI magnitude of the unit chosen for $k$). Any coherent system gives 1: choose base units for length, time and mass, derive $u_p = M/(L T^2)$, $u_μ = M/(L T)$ and $u_k = L^2$ from them.

System	$u_L$	$u_t$	$u_p$	$u_k$	$u_μ$	rate	$C$
SI	m	s	Pa	m²	Pa·s	m²/s	`1`
CGS	cm	s	barye	cm²	poise	cm²/s	`1`
MTS	m	s	pièze	m²	pz·s	m²/s	`1`
mm-ms-g	mm	ms	MPa	mm²	kPa·s	mm²/ms	`1`
Darcy's own	cm	s	atm	darcy	cP	cm²/s	`1`
metric	m	day	bar	mD	cP	m²/day	`0.008527`
field-like	ft	day	psi	mD	cP	ft²/day	`0.006328`
lab	cm	hour	atm	mD	cP	cm²/hour	`3.6`

The rate unit is forced to $u_L^2/u_t$ -- an areal rate.

A well rate of 20 for a 25 m thick reservoir means 500 m³/day.

$C$ enters at exactly 2 sites, both of them Darcy's law.

The transmissibilities of TPFA and the well index of minires.wells.peaceman_WI. Everything else is derivative, and already consistent.

fluid: Any = None

The two-phase fluid: a minires.fluids.Fluid, holding the viscosities and the Corey relative permeability parameters, and computing the mobilities and fractional flow from them.

Assigning a dict builds one (ResSim(fluid=dict(vo=5, swc=.2))), and None the default (unit viscosities, the reference paper's quadratic curves). Its fields stay writable (model.fluid.vo = 5). A Fluid subclass (e.g. with tabulated curves) may be assigned instead.

ct: float = 0.0

Total (rock + fluids) compressibility, $c_t$, as a single constant.

The default, 0, yields the incompressible model, whose pressure eqn. is elliptic: pressure is defined only up to an additive constant, and the sources/sinks must balance. Setting ct > 0 yields the slightly compressible model: the pressure eqn. gains the accumulation term $ φ \, c_t \, ∂p/∂t $ (discretized by backward Euler over the same dt as the saturation step, which is what makes P0 of sim consequential), and the transport eqn. the matching storage term, charged to the phases in proportion to their saturation (ref storage_rate). Injection and production then need not balance, enabling e.g. primary depletion.

Derivation and fidelity ($ c_t \, Δp \ll 1 $, which the voidage sets, not ct): ref the "Compressibility" section of the docs.

The drive mechanism is whatever supplies the energy that pushes the hydrocarbons to the well. Fluid and rock expansion (a.k.a. depletion drive), which is what $ c_t > 0 $ enables here in the absence of injection, is the weakest, recovering only a few percent, because $c_t$ is so small. Stronger ones are solution gas drive, gas cap drive, water drive (aquifers, ref minires.wells.aquifer_WI), gravity drainage, and compaction drive (which manifests as seabed subsidence). Recovery is staged: primary production runs on the native drive. Upon drilling, unless valves are in place, the pressure of the initial equilibrium will cause a blow out; a new equilibrium is usually attained when 20% of the hydrocarbons have been produced, which marks the end of the primary production. Secondary recovery adds pressure support by injecting water or gas (waterflooding being the case simulated here); tertiary, or EOR (enhanced oil recovery), alters the flow physics itself, e.g. by polymer, surfactant, or CO₂ injection. The voidage replacement ratio is the injected reservoir volume divided by the produced one; $\mathrm{VRR} = 1$ is exactly the balance, $\sum q = 0$, that the incompressible model is obliged to impose. The zero-dimensional (single tank) accounting of all of the above, used to estimate reserves without a grid, is called material balance.

cached_precond: bool = True

Solve the pressure system iteratively, preconditioned by a cached factorization.

The alternative (False) is a fresh sparse direct factorization each time step. But the system changes slowly -- only through the mobility $λ(s)$, i.e. where the front has moved -- so the factorization of an earlier step remains an excellent preconditioner: with it, conjugate gradients converges in 1--2 iterations if the saturation is (nearly) static, and in ~12 behind a moving front (each the cost of one back-substitution, i.e. 1/30 of a factorization). The factorization is refreshed only when convergence fails (as it does when the mobility has drifted too far, e.g. after many steps at a strong viscosity contrast), so the result is exact to the solver tolerance (1e-10) either way. It requires the system to be SPD, which the TPFA system is (ref the "How to solve" section of the docs).

The cache is per instance (_pLU) and is dropped on pickling and deepcopy, a SuperLU not being picklable; a copy simply refactorizes on its next step. The measurements (2--6x on the well-test and depletion examples, 15--40% on the waterfloods), and what was tried besides, are recorded in tests/test_precond.py.

K: Any = None

Permeabilities (in x and y directions). Array of shape (2, Nx, Ny)).

Permeability, denoted by tensor $\mathbf{K}$, quantifies transmissibility. Usually SPD -- here diagonal, $\mathrm{diag}(K_x, K_y)$ per cell, which is what the two components of the first axis are -- and correlated with por. Among the reservoir rocks, sandstone usually have large, well-connected pores, and high permeability, shale is nearly impermeable, like cap rock and bed rock. Permeability is measured in Darcy ($≈ 10^{-12} m^2$), ref cdarcy. A medium is called isotropic if $\mathbf{K}$ is scalar. Its spatial structure goes by yet more lingo: facies (bodies of rock of like character, e.g. the channels of a fluvial deposit), and the fissures and fractures whose conductivity dwarfs that of the rock around them. Being the most uncertain of the reservoir's properties (ref minires.fluids), it is the parameter that history matching estimates -- and the one that minires.tlm differentiates (as $\log \mathbf{K}$) with respect to.

por: Any = None

Porosity; Array of shape (Nx, Ny)).

Porosity, $φ$, is the void volume fraction. Depends on pressure, because rock is compressible. Compressibility is the porosity's (relative) gradient wrt. pressure. Usually neglected, as it is here in the pore volume (ref pore_volume), so that $φ$ is a constant, but spatial, field; the rock's compressibility enters instead through ct, which lumps it with the fluids'.

active: Any = None

Mask of the active cells, (Nx, Ny), boolean. Default: all True.

Setting some cells inactive carves an irregular reservoir out of the rectangular grid -- an outline, holes, or a (sealing) fault: a line of inactive cells, which blocks flow so long as it is unbroken (a diagonal staircase suffices, fluxes passing only through faces). Inactive cells take no part in the physics:

The faces to them carry zero transmissibility, hence zero flux, so the reservoir is closed along their perimeter as it is along the boundary.
They have no equations. The pressure system carries an identity row for each (keeping it well conditioned, unlike a tiny permeability would), and their pore volume is taken as infinite (ref pore_volume; so they never bind the CFL, as a tiny porosity would). So their state is simply carried through sim unchanged from S0 and P0. Plots mask them out.
A well may not be completed in one (_validate checks). Should the active cells form several disconnected regions, each is a reservoir of its own, and only the first is pinned (ref TPFA): if ct == 0, the others must balance their own rates, or hold a BHP well -- which is not checked as such, but the singular system it would otherwise make is caught by the residual check of _solve_pressure.

The arrays (K, por, S0, ...) keep the full grid shape, and the flat index (xy2ind) runs over all cells: the mask selects the active ones.

wells: Any = None

The wells: a Wells, holding the flat, per-completion arrays -- positions, rates, pressures, well indices -- that the model runs on.

Assigning a list (or dict) of records -- one per well -- assembles one, which is the convenient way to configure them; the record format is documented in minires.wells.Wells.from_records. Assigning None empties it. A Wells may also be given directly, in which case it is bound to this model, whereupon its completions snap onto the grid.

The arrays remain writable throughout (model.wells.rates = ...), as an ensemble or optimisation loop requires.

def pore_volume(self) -> numpy.ndarray: View Source

269    def pore_volume(self) -> np.ndarray:
270        """Pore volume (per unit thickness) of each cell, `h2 * por`. Flat.
271
272        `inf` for inactive cells (ref `active`): dividing by it, the transport
273        schemes leave their saturation alone, and the CFL estimate ignores them.
274        """
275        pv = self.h2 * self.por.ravel()
276        return np.where(self.active.ravel(), pv, np.inf)

Pore volume (per unit thickness) of each cell, h2 * por. Flat.

inf for inactive cells (ref active): dividing by it, the transport schemes leave their saturation alone, and the CFL estimate ignores them.

nComp View Source

278    nComp = property(lambda self: self.wells.nComp)

Num. of well completions, i.e. the rows of every array the model indexes by them -- which is what it actually solves for, the equations being assembled per completion. Forwarded from minires.wells.Wells.nComp.

def assemble_wells(self, S: numpy.ndarray | None, P: numpy.ndarray | None, k: int) -> None: View Source

285    def assemble_wells(
286        self, S: np.ndarray | None, P: np.ndarray | None, k: int
287    ) -> None:
288        """Set up (for time `k`) the wells' contributions to the equations.
289
290        The controls are those of `well_controls`, to which `S` and `P` (the
291        state at the *start* of the step) are simply passed on.
292        Rate-controlled wells enter the source/sink *field*, `_Q`, directly.
293        BHP-controlled ones (ref `minires.wells.Wells.bhp`) cannot: their
294        rate is not yet known. They instead enter the pressure equations in
295        `TPFA`, after which `realize_bhp` folds the resulting rate into `_Q`.
296        """
297        ctrl = self.well_controls(S, P, k)
298        inds = self.xy2ind(*self.wells.xy.T)
299        rates, p_bh = ctrl["rates"], ctrl["bhp"]
300        is_bhp = np.isfinite(p_bh)
301        assert np.isfinite(rates[~is_bhp]).all(), (
302            "A rate-controlled well has a non-finite rate. Give it a number"
303            " (`0` shuts it in), or put it on BHP control; ref `Wells.rates`."
304        )
305
306        # The well model's constant of proportionality, WI * λ_t.
307        # NB: `nan` marks the rate-controlled wells, throughout.
308        WI_lam = np.full(self.nComp, np.nan)
309        if is_bhp.any():
310            WI = self.wells.WI
311            assert WI is not None and np.isfinite(WI[is_bhp]).all(), (
312                "BHP control requires (finite) `Wells.WI`."
313            )
314            assert S is not None, "BHP control requires `S` (for λ_t)."
315            Mw, Mo = self.fluid.RelPerm(S)
316            WI_lam[is_bhp] = WI[is_bhp] * (Mw + Mo)[inds[is_bhp]]
317
318        # Translate well conditions for cells.
319        # NB: Dont use `Q[inds] += ...` since `inds` may contain dupes.
320        self._Q, bhp_diag, bhp_rhs = np.zeros((3, self.Nxy))
321        np.add.at(self._Q, inds[~is_bhp], rates[~is_bhp])
322        np.add.at(bhp_diag, inds[is_bhp], WI_lam[is_bhp])
323        np.add.at(bhp_rhs, inds[is_bhp], (WI_lam * p_bh)[is_bhp])
324        rates[is_bhp] = np.nan  # only `realize_bhp` knows these
325        self._wells_now: dict[str, np.ndarray] = dict(
326            inds=inds, rates=rates, p_bh=p_bh,
327            WI_lam=WI_lam, bhp_diag=bhp_diag, bhp_rhs=bhp_rhs,
328        )  # fmt: skip

Set up (for time k) the wells' contributions to the equations.

The controls are those of well_controls, to which S and P (the state at the start of the step) are simply passed on. Rate-controlled wells enter the source/sink field, _Q, directly. BHP-controlled ones (ref minires.wells.Wells.bhp) cannot: their rate is not yet known. They instead enter the pressure equations in TPFA, after which realize_bhp folds the resulting rate into _Q.

def realize_bhp(self, P: numpy.ndarray) -> None: View Source

330    def realize_bhp(self, P: np.ndarray) -> None:
331        """Compute rates for BHP wells. Enter into `_Q` and `_wells_now["rates"]`.
332
333        The rate, $ WI λ_t (p_\\mathrm{bh} - p_\\mathrm{cell}) $, is signed by
334        nature (ref the `minires.wells.Wells.bhp` warning).
335
336        By construction of the linear system of `TPFA`, this leaves `_Q` equal
337        to the *total* well flux, which is what keeps `storage_rate` -- and
338        hence the transport step -- consistent with the pressure solution.
339        """
340        wls = self._wells_now
341        WI_lam = wls["WI_lam"]  # `nan` marks the rate-controlled wells
342        # Insert in cell source/sink field
343        self._Q = self._Q + wls["bhp_rhs"] - wls["bhp_diag"] * P
344        # Insert in per-well rates
345        is_bhp = np.isfinite(WI_lam)
346        wls["rates"][is_bhp] = (WI_lam * (wls["p_bh"] - P[wls["inds"]]))[is_bhp]

Compute rates for BHP wells. Enter into _Q and _wells_now["rates"].

The rate, $ WI λ_t (p_\mathrm{bh} - p_\mathrm{cell}) $, is signed by nature (ref the minires.wells.Wells.bhp warning).

By construction of the linear system of TPFA, this leaves _Q equal to the total well flux, which is what keeps storage_rate -- and hence the transport step -- consistent with the pressure solution.

def well_controls(self, S: numpy.ndarray | None, P: numpy.ndarray | None, k: int) -> dict: View Source

366    def well_controls(self, S: np.ndarray | None, P: np.ndarray | None, k: int) -> dict:
367        """Compute the wells' controls for time `k`: `dict(rates=..., bhp=...)`.
368
369        Each is a `(nComp,)` array, read off the specifications --
370        `minires.wells.Wells.rates`, `minires.wells.Wells.bhp` -- which
371        are *open-loop*: fixed before the simulation begins. Overriding
372        (patching/subclassing) this method is therefore how to do *feedback*
373        control, the controls being free to depend on the state at the *start*
374        of the step: the saturation `S` and the pressure `P`.
375        The returned arrays are copies, so they may be modified in place.
376
377        Most feedback concerns the rates alone -- e.g. shutting the wells upon
378        water breakthrough at the producer:
379
380        >>> class Shutter(ResSim):
381        ...     def well_controls(self, S, P, k):
382        ...         ctrl = super().well_controls(S, P, k)
383        ...         if S is not None and S[self.xy2ind(1, 1)] > .5:
384        ...             ctrl["rates"][:] = 0    # NB: all of them! See warning
385        ...         return ctrl
386        >>> model = Shutter(Lx=1, Ly=1, Nx=16, Ny=16,
387        ...                 wells=Wells(xy=[[0, 0], [1, 1]], rates=[[1], [-1]]))
388        >>> SS, PP = model.sim(.05, 20, model.fluid.swc*np.ones(model.Nxy), pbar=False)
389        >>> int((model.wells.actual_rates[1] == 0).argmax())  # step of breakthrough
390        16
391
392        But the `bhp` is here too, and with it each well's *control mode*
393        (`nan` => rate-controlled, ref `minires.wells.Wells.bhp`) -- which
394        is what an approximate mode *switch* requires. For example, rate
395        control with a BHP limit -- the industrial default -- wherein a
396        producer holds its rate only for as long as that does not draw it
397        below some `p_min`:
398
399        >>> class Limited(ResSim):
400        ...     p_min = .5
401        ...     def well_controls(self, S, P, k):
402        ...         ctrl = super().well_controls(S, P, k)
403        ...         if P is None:
404        ...             return ctrl                    # nothing to switch on
405        ...         p_bh = self.bhp(S, P, ctrl["rates"])
406        ...         switch = p_bh < self.p_min         # the rate is unsustainable
407        ...         ctrl["bhp"] = np.where(switch, self.p_min, np.nan)
408        ...         return ctrl
409        >>> model = Limited(Lx=1, Ly=1, Nx=16, Ny=16, ct=.1,
410        ...                 wells=Wells(xy=[[.5, .5]], rates=[[-.25]]))
411        >>> from minires import peaceman_WI
412        >>> model.wells.WI = peaceman_WI(model, model.wells.xy, rw=1e-3)
413        >>> SS, PP = model.sim(.02, 25, np.zeros(model.Nxy),
414        ...                    P0=np.ones(model.Nxy), pbar=False)
415
416        The well delivers its target rate until the limit binds, and declines
417        thereafter -- at constant $ p_\\mathrm{bh} $, exponentially so
418        (ref `examples/well_control.py`, which plots all three modes):
419
420        >>> (-model.wells.actual_rates[0, [0, 5, 6, -1]]).round(3)
421        array([0.25 , 0.25 , 0.182, 0.005])
422
423        .. warning:: With `ct == 0` the rates must still sum to 0 at every step.
424
425            Ref `minires.wells.Wells.rates`. So shutting one well requires
426            matching it on the other side -- as above.
427
428        .. note:: The mode switch lags the solve by one step.
429
430            It is decided from the previous step's pressure, whereas the well
431            model itself is solved *simultaneously* with the new one. So the
432            limit is breached for the one step in which it comes to bind.
433            Shorten `dt` to refine.
434
435        .. note:: `S` and `P` may be `None`, so an override should tolerate that.
436
437            They are `None` if the caller has none to offer -- as when
438            `assemble_wells` is used merely to set up a plot.
439
440        .. note:: Setting both controls for a well is not an error, just pointless.
441
442            `assemble_wells` discards the rate of a BHP-controlled well -- it
443            is `realize_bhp` that fills it in.
444        """
445        return dict(
446            rates=self.wells.at_time("rates", 0.0, k),
447            bhp=self.wells.at_time("bhp", np.nan, k),
448        )

Compute the wells' controls for time k: dict(rates=..., bhp=...).

Each is a (nComp,) array, read off the specifications -- minires.wells.Wells.rates, minires.wells.Wells.bhp -- which are open-loop: fixed before the simulation begins. Overriding (patching/subclassing) this method is therefore how to do feedback control, the controls being free to depend on the state at the start of the step: the saturation S and the pressure P. The returned arrays are copies, so they may be modified in place.

Most feedback concerns the rates alone -- e.g. shutting the wells upon water breakthrough at the producer:

>>> class Shutter(ResSim):
...     def well_controls(self, S, P, k):
...         ctrl = super().well_controls(S, P, k)
...         if S is not None and S[self.xy2ind(1, 1)] > .5:
...             ctrl["rates"][:] = 0    # NB: all of them! See warning
...         return ctrl
>>> model = Shutter(Lx=1, Ly=1, Nx=16, Ny=16,
...                 wells=Wells(xy=[[0, 0], [1, 1]], rates=[[1], [-1]]))
>>> SS, PP = model.sim(.05, 20, model.fluid.swc*np.ones(model.Nxy), pbar=False)
>>> int((model.wells.actual_rates[1] == 0).argmax())  # step of breakthrough
16

But the bhp is here too, and with it each well's control mode (nan => rate-controlled, ref minires.wells.Wells.bhp) -- which is what an approximate mode switch requires. For example, rate control with a BHP limit -- the industrial default -- wherein a producer holds its rate only for as long as that does not draw it below some p_min:

>>> class Limited(ResSim):
...     p_min = .5
...     def well_controls(self, S, P, k):
...         ctrl = super().well_controls(S, P, k)
...         if P is None:
...             return ctrl                    # nothing to switch on
...         p_bh = self.bhp(S, P, ctrl["rates"])
...         switch = p_bh < self.p_min         # the rate is unsustainable
...         ctrl["bhp"] = np.where(switch, self.p_min, np.nan)
...         return ctrl
>>> model = Limited(Lx=1, Ly=1, Nx=16, Ny=16, ct=.1,
...                 wells=Wells(xy=[[.5, .5]], rates=[[-.25]]))
>>> from minires import peaceman_WI
>>> model.wells.WI = peaceman_WI(model, model.wells.xy, rw=1e-3)
>>> SS, PP = model.sim(.02, 25, np.zeros(model.Nxy),
...                    P0=np.ones(model.Nxy), pbar=False)

The well delivers its target rate until the limit binds, and declines thereafter -- at constant $ p_\mathrm{bh} $, exponentially so (ref examples/well_control.py, which plots all three modes):

>>> (-model.wells.actual_rates[0, [0, 5, 6, -1]]).round(3)
array([0.25 , 0.25 , 0.182, 0.005])

With `ct == 0` the rates must still sum to 0 at every step.

Ref minires.wells.Wells.rates. So shutting one well requires matching it on the other side -- as above.

The mode switch lags the solve by one step.

It is decided from the previous step's pressure, whereas the well model itself is solved simultaneously with the new one. So the limit is breached for the one step in which it comes to bind. Shorten dt to refine.

`S` and `P` may be `None`, so an override should tolerate that.

They are None if the caller has none to offer -- as when assemble_wells is used merely to set up a plot.

Setting both controls for a well is not an error, just pointless.

assemble_wells discards the rate of a BHP-controlled well -- it is realize_bhp that fills it in.

def bhp( self, S: numpy.ndarray, P: numpy.ndarray, rates: numpy.ndarray) -> numpy.ndarray: View Source

450    def bhp(self, S: np.ndarray, P: np.ndarray, rates: np.ndarray) -> np.ndarray:
451        """Bottom-hole pressures implied by the (signed) `rates`, via the well indices.
452
453        I.e. the well model of `minires.wells.Wells.WI`, solved for
454        $ p_\\mathrm{bh} $:
455        the rate's sign puts an injector above, a producer below, its cell
456        pressure. `nan` wherever the well index is unset.
457
458        `S` and `P` (both flat) should be the saturation and the pressure of the
459        *same* `pressure_step`, i.e. `SS[k]` and `PP[k+1]` of `sim` -- which is
460        what `actual_bhp` records, so prefer reading that.
461
462        .. warning:: $ λ_t $ is that of the well's *cell*.
463
464            So an injector's injectivity is governed by the mobility of
465            whatever the cell currently holds, rather than by that of the
466            injectant.
467        """
468        if self.wells.WI is None:
469            return np.full(self.nComp, np.nan)
470        Mw, Mo = self.fluid.RelPerm(S)
471        ii = self.xy2ind(*self.wells.xy.T)
472        return P[ii] + rates / (self.wells.WI * (Mw + Mo)[ii])

Bottom-hole pressures implied by the (signed) rates, via the well indices.

I.e. the well model of minires.wells.Wells.WI, solved for $ p_\mathrm{bh} $: the rate's sign puts an injector above, a producer below, its cell pressure. nan wherever the well index is unset.

S and P (both flat) should be the saturation and the pressure of the same pressure_step, i.e. SS[k] and PP[k+1] of sim -- which is what actual_bhp records, so prefer reading that.

$ λ_t $ is that of the well's cell.

So an injector's injectivity is governed by the mobility of whatever the cell currently holds, rather than by that of the injectant.

def pressure_step( self, S: numpy.ndarray, P: numpy.ndarray | None = None, dt: float | None = None) -> tuple[numpy.ndarray, minires.grid.Fluxes]: View Source

475    def pressure_step(
476        self,
477        S: np.ndarray,
478        P: np.ndarray | None = None,
479        dt: float | None = None,
480    ) -> tuple[np.ndarray, Fluxes]:
481        """Compute permeabilities then solve Darcy's equation. Returns `[P, V]`.
482
483        `P` (flat, like `S`) is the *previous* step's pressure: used (and
484        required) only if `ct > 0`, along with `dt`. The new one replaces it.
485        """
486        # Compute K*λ(S)
487        Mw, Mo = self.fluid.RelPerm(S)
488        Mt = Mw + Mo
489        Mt = Mt.reshape(self.shape)
490        KM = Mt * self.K
491        # Compute pressure and extract fluxes
492        [P, V] = self.TPFA(KM, P, dt)
493        return P, V

Compute permeabilities then solve Darcy's equation. Returns [P, V].

P (flat, like S) is the previous step's pressure: used (and required) only if ct > 0, along with dt. The new one replaces it.

def TPFA( self, K: numpy.ndarray, P: numpy.ndarray | None = None, dt: float | None = None) -> tuple[numpy.ndarray, minires.grid.Fluxes]: View Source

512    def TPFA(
513        self,
514        K: np.ndarray,
515        P: np.ndarray | None = None,
516        dt: float | None = None,
517    ) -> tuple[np.ndarray, Fluxes]:
518        """Two-point flux-approximation (TPFA) of Darcy: $ -∇(K ∇u) = q $
519
520        i.e. steady-state diffusion w/ nonlinear coefficient, $K$,
521        if `ct == 0`. Otherwise (slightly compressible model) solve
522        the backward-Euler step of $ φ c_t ∂u/∂t - ∇(K ∇u) = q $,
523        which requires the previous pressure, `P`, and `dt`.
524
525        After solving for pressure `P`, extract the fluxes `V`
526        by finite differences.
527        """
528        # Compute transmissibilities by harmonic averaging.
529        C = self.cdarcy
530        L = 1 / K
531        TX = np.zeros((self.Nx + 1, self.Ny))
532        TY = np.zeros((self.Nx, self.Ny + 1))
533        TX[1:-1, :] = C * 2 * self.hy / self.hx / (L[0, :-1, :] + L[0, 1:, :])
534        TY[:, 1:-1] = C * 2 * self.hx / self.hy / (L[1, :, :-1] + L[1, :, 1:])
535        # No flow across the faces of inactive cells (ref `active`)
536        act = self.active
537        TX[1:-1, :] *= act[:-1, :] & act[1:, :]
538        TY[:, 1:-1] *= act[:, :-1] & act[:, 1:]
539
540        # Assemble TPFA discretization matrix.
541        x1 = TX[:-1, :].ravel()
542        x2 = TX[1:, :].ravel()
543        y1 = TY[:, :-1].ravel()
544        y2 = TY[:, 1:].ravel()
545
546        # Setup linear system
547        DiagVecs = [-x2, -y2, y1 + y2 + x1 + x2, -y1, -x1]
548        DiagIndx = [-self.Ny, -1, 0, 1, self.Ny]
549        q = self._Q
550        if self.ct > 0:
551            # Accumulation term (φ ct h²/dt) of backward Euler.
552            # Renders the system nonsingular (unlike the pure-Neumann problem).
553            accum = self.por.ravel() * self.ct * self.h2 / dt
554            DiagVecs[2] = DiagVecs[2] + accum
555            q = q + accum * P
556        elif not self._wells_now["bhp_diag"].any():
557            # Pin the (o/w pure-Neumann & singular) problem, at the 1st active cell
558            i = self._pin
559            DiagVecs[2][i] += np.sum(self.K.reshape(2, -1)[:, i])  # ref article p. 13
560        # Well model of the BHP-controlled wells
561        DiagVecs[2] = DiagVecs[2] + self._wells_now["bhp_diag"]
562        q = q + self._wells_now["bhp_rhs"]
563        # Inactive cells have no equation. Identity rows (their faces being
564        # already closed) keep them at their previous pressure, and the system SPD.
565        act = act.ravel()
566        DiagVecs[2] = np.where(act, DiagVecs[2], 1.0)
567        q = np.where(act, q, 0.0 if P is None else P)
568
569        # Solve; compute A\q to update P
570        A = self._spdiags(DiagVecs, DiagIndx)
571        P = self._solve_pressure(A.tocsr(), q, P)
572        # P = np.linalg.solve(A.A, q) # direct dense solver
573        # Could also try scipy.linalg.solveh_banded which, according to
574        # https://scicomp.stackexchange.com/a/30074 uses the Thomas algorithm,
575        # as recommended by Aziz and Settari ("Petro. Res. simulation").
576        # NB: stackexchange also mentions that solve_banded does not work well
577        # when the band offsets large, i.e. higher-dimensional problems.
578
579        # Extract fluxes, via a grid-shaped view of the (flat) pressure.
580        P2d = P.reshape(self.shape)
581        V = Fluxes(
582            x=np.zeros((self.Nx + 1, self.Ny)),
583            y=np.zeros((self.Nx, self.Ny + 1)),
584        )
585        V.x[1:-1, :] = (P2d[:-1, :] - P2d[1:, :]) * TX[1:-1, :]
586        V.y[:, 1:-1] = (P2d[:, :-1] - P2d[:, 1:]) * TY[:, 1:-1]
587        return P, V

Two-point flux-approximation (TPFA) of Darcy: $ -∇(K ∇u) = q $

i.e. steady-state diffusion w/ nonlinear coefficient, $K$, if ct == 0. Otherwise (slightly compressible model) solve the backward-Euler step of $ φ c_t ∂u/∂t - ∇(K ∇u) = q $, which requires the previous pressure, P, and dt.

After solving for pressure P, extract the fluxes V by finite differences.

def upwind_diff(self, V: minires.grid.Fluxes) -> scipy.sparse._dia.dia_matrix: View Source

636    def upwind_diff(self, V: Fluxes) -> sparse.dia_matrix:
637        """Upwind finite-volume scheme."""
638        fp = self._Q.clip(max=0)  # production
639        # Flow fluxes, separated into direction (x-y) and sign
640        x1 = V.x.clip(max=0)[:-1, :].ravel()
641        y1 = V.y.clip(max=0)[:, :-1].ravel()
642        x2 = V.x.clip(min=0)[1:, :].ravel()
643        y2 = V.y.clip(min=0)[:, 1:].ravel()
644        DiagVecs = [x2, y2, fp + y1 - y2 + x1 - x2, -y1, -x1]
645        DiagIndx = [-self.Ny, -1, 0, 1, self.Ny]
646        A = self._spdiags(DiagVecs, DiagIndx)
647        return A

Upwind finite-volume scheme.

def storage_rate(self, V: minires.grid.Fluxes) -> numpy.ndarray: View Source

649    def storage_rate(self, V: Fluxes) -> np.ndarray:
650        """The volume rate, per cell, that goes into storage: $ q - ∇ ⋅ V $.
651
652        For the incompressible model this is `0`: the fluxes balance the wells
653        exactly, cell by cell. With `ct > 0` it is (by construction of the
654        linear system of `TPFA`) the accumulation term of the backward-Euler
655        step, $ φ \\, c_t \\, h^2 \\, (p^{n+1} - p^n) / Δt $, which the saturation
656        steps charge to the phases (ref `ct`).
657
658        Computing it from `V` (rather than from $p^{n+1} - p^n$) means it is
659        *exactly* the imbalance seen by the transport scheme, whose `upwind_diff`
660        is assembled from the same fluxes.
661        """
662        if self.ct == 0:
663            return np.zeros(self.Nxy)
664        divV = (V.x[1:, :] - V.x[:-1, :]) + (V.y[:, 1:] - V.y[:, :-1])
665        return self._Q - divV.ravel()

The volume rate, per cell, that goes into storage: $ q - ∇ ⋅ V $.

For the incompressible model this is 0: the fluxes balance the wells exactly, cell by cell. With ct > 0 it is (by construction of the linear system of TPFA) the accumulation term of the backward-Euler step, $ φ \, c_t \, h^2 \, (p^{n+1} - p^n) / Δt $, which the saturation steps charge to the phases (ref ct).

Computing it from V (rather than from $p^{n+1} - p^n$) means it is exactly the imbalance seen by the transport scheme, whose upwind_diff is assembled from the same fluxes.

def estimate_1CFL( self, pv: numpy.ndarray, V: minires.grid.Fluxes, fi: numpy.ndarray) -> float: View Source

668    def estimate_1CFL(self, pv: np.ndarray, V: Fluxes, fi: np.ndarray) -> float:
669        """Estimate 1/CFL for use with `saturation_step_upwind`."""
670        # In-/Out-flux x-/y- faces
671        XP = V.x.clip(min=0)
672        XN = V.x.clip(max=0)
673        YP = V.y.clip(min=0)
674        YN = V.y.clip(max=0)
675        Vi = XP[:-1, :] + YP[:, :-1] - XN[1:, :] - YN[:, 1:]
676
677        flx = max((Vi.ravel() + fi) / pv)  # estimate of influx
678        # NB: `storage_rate` is not counted here. In practice it is a small
679        # fraction of the fluxes that are (under 20% even at `ct = 10`),
680        # so the safety factor below covers it.
681        # The characteristic speed is f_w'(s) times the velocity: bound it by the
682        # maximal slope of f_w, sampled at 1001 points of the mobile range (which
683        # undershoots the true maximum by O(1e-6), nothing against the safety
684        # factor of 1.5). For the default curves at equal viscosities this is the
685        # paper's 3 / (1 - swc - sor) (Listing 8), which steeper curves (unequal
686        # viscosities, higher Corey exponents) exceed.
687        ss = np.linspace(self.fluid.swc, 1 - self.fluid.sor, 1001)
688        dfw_max = self.fluid.dfractional_flow(ss).max()
689        cfl = 1.5 * dfw_max * flx
690        # NB: the ceiling is nudged down by a relative epsilon, so that a `dt`
691        # sitting *on* an integer multiple of the CFL limit (as the examples'
692        # round numbers tend to) does not gain a whole extra sub-step from the
693        # last bits of the linear solve -- which differ across platforms and
694        # library versions, and would make the results irreproducible.
695        # The safety factor covers the shaving, the sampling's undershoot, and
696        # the storage rate.
697        return cfl * (1 - 1e-9)

Estimate 1/CFL for use with saturation_step_upwind.

def saturation_step_upwind( self, S: numpy.ndarray, V: minires.grid.Fluxes, dt: float) -> numpy.ndarray: View Source

700    def saturation_step_upwind(self, S: np.ndarray, V: Fluxes, dt: float) -> np.ndarray:
701        """Explicit upwind FV discretisation of conserv. of mass (water sat.)."""
702        # fmt: off
703        A  = self.upwind_diff(V)                 # FV discretized transport operator
704        pv = self.pore_volume()                  # Pore volume (per thickness)
705        fi = self._Q.clip(min=0)                 # Well inflow
706        st = self.storage_rate(V)                # Storage (0 if incompressible)
707
708        # Compute sub/local dt
709        cfl1 = self.estimate_1CFL(pv, V, fi)
710        nT = int(np.ceil(dt * cfl1))
711        nT = max(1, nT)
712
713        # Scale A
714        dtx = dt / nT / pv                       # timestep / pore volume
715        B   = self._spdiags(dtx, 0) @ A          # A * dt/|Omega i|
716
717        for _ in range(nT):
718            fw = self.fluid.fractional_flow(S)      # fractional flow
719            S = S + (B@fw + (fi - S*st)*dtx)     # update saturation
720        # fmt: on
721        return S

Explicit upwind FV discretisation of conserv. of mass (water sat.).

def saturation_step_implicit( self, S: numpy.ndarray, V: minires.grid.Fluxes, dt: float, nNewtonMax: int = 10, nTmax_log2: int = 10) -> numpy.ndarray: View Source

724    def saturation_step_implicit(
725        self,
726        S: np.ndarray,
727        V: Fluxes,
728        dt: float,
729        nNewtonMax: int = 10,
730        nTmax_log2: int = 10,
731    ) -> np.ndarray:
732        """Implicit FV discretisation of conserv. of mass (water sat.).
733
734        .. warning:: The Newton iteration can converge to a spurious root.
735
736            Far outside the $ c_t \\, Δp \\ll 1 $ regime (ref `ct`), it may
737            converge -- silently -- to a root of the residual outside $[0, 1]$:
738            the polynomial `minires.fluids.Fluid.RelPerm` extends smoothly
739            beyond the unit interval, and the sub-`dt` halving only triggers on
740            *non*-convergence.
741            The explicit scheme
742            (`saturation_step_upwind`), being monotone, stays within $[0, 1]$
743            even for extreme `ct`.
744
745        .. note:: This scheme rarely earns its keep.
746
747            It is usually both slower & less accurate than `saturation_step_upwind`.
748            Both schemes sub-divide `dt` internally, so it is not `dt` that decides
749            their cost but the stiffness of the grid -- and the well cells, being
750            normally the stiffest, hold the two requirements within a factor 3 of each
751            other (that being the safety margin of `ResSim.estimate_1CFL`), while an
752            implicit sub-step -- a sparse solve, or several -- 10 or 100x more.
753            It does pay off where the stiffest cell is *not* a well cell -- a tight
754            streak, a fracture, a locally refined region -- running 7 times faster at a
755            1000x pore-volume contrast. The branch `implicit-transport-scheme` says more.
756        """
757        # fmt: off
758        A  = self.upwind_diff(V)                 # FV discretized transport operator
759        pv = self.pore_volume()                  # Pore volume (per thickness)
760        fi = self._Q.clip(min=0)                 # Well inflow
761        st = self.storage_rate(V)                # Storage (0 if incompressible)
762
763        # For each iter, halve the sub/local dt
764        for nT_log2 in range(0, nTmax_log2):
765            nT = 2**nT_log2
766
767            # Scale A
768            dtx = dt / nT / pv                   # timestep / pore volume
769            B   = self._spdiags(dtx, 0) @ A      # A * dt/|Omega i|
770            C   = self._spdiags(dtx*st, 0)       # storage, likewise scaled
771
772            Sn = S
773            for _ in range(nT):
774                Sp = Sn
775                for _ in range(nNewtonMax):
776                    fw = self.fluid.fractional_flow(Sn)     # fract. flow
777                    df = self.fluid.dfractional_flow(Sn)    # its derivative
778                    dG = (sparse.eye(self.Nxy) + C                        # deriv of G
779                          - B @ self._spdiags(df, 0))
780                    G  = Sn - Sp - (B@fw + (fi - Sn*st)*dtx)  # G(s)
781                    dS = spsolve(dG, G)             # compute dS
782                    Sn = Sn - dS                    # update S
783
784                    if np.sqrt(sum(dS**2)) < 1e-3:
785                        # If converged: halt Newton iterations
786                        break
787                else:
788                    # If never converged: increase nT, restart time loop
789                    break
790            else:
791                # If completed all time steps, halt
792                break
793        else:
794            # Failed (even with max nT) to complete all time steps
795            print("Warning: did not converge")
796        # fmt: on
797
798        return Sn

Implicit FV discretisation of conserv. of mass (water sat.).

The Newton iteration can converge to a spurious root.

Far outside the $ c_t \, Δp \ll 1 $ regime (ref ct), it may converge -- silently -- to a root of the residual outside $[0, 1]$: the polynomial minires.fluids.Fluid.RelPerm extends smoothly beyond the unit interval, and the sub-dt halving only triggers on non-convergence. The explicit scheme (saturation_step_upwind), being monotone, stays within $[0, 1]$ even for extreme ct.

This scheme rarely earns its keep.

It is usually both slower & less accurate than saturation_step_upwind. Both schemes sub-divide dt internally, so it is not dt that decides their cost but the stiffness of the grid -- and the well cells, being normally the stiffest, hold the two requirements within a factor 3 of each other (that being the safety margin of ResSim.estimate_1CFL), while an implicit sub-step -- a sparse solve, or several -- 10 or 100x more. It does pay off where the stiffest cell is not a well cell -- a tight streak, a fracture, a locally refined region -- running 7 times faster at a 1000x pore-volume contrast. The branch implicit-transport-scheme says more.

def time_stepper(self, dt: float, implicit: bool = False) -> Callable: View Source

814    def time_stepper(self, dt: float, implicit: bool = False) -> Callable:
815        """Get ODE solver (integrator) for model.
816
817        Whatever time step `dt` is given, both schemes will use smaller steps internally.
818
819        - `explicit`: computes sub-`dt` based on CFL esitmate.
820        - `implicit`: reduces sub-`dt` until convergence is achieved.
821        """
822
823        def integrate(S, P, k):
824            self.assemble_wells(S, P, k)
825            self._validate()
826            [P, V] = self.pressure_step(S, P, dt)
827            self.realize_bhp(P)
828            self._record_actual_well_operation(S, P, k)
829            if implicit:
830                S = self.saturation_step_implicit(S, V, dt)
831            else:
832                S = self.saturation_step_upwind(S, V, dt)
833            return S, P
834
835        return integrate

Get ODE solver (integrator) for model.

Whatever time step dt is given, both schemes will use smaller steps internally.

explicit: computes sub-dt based on CFL esitmate.
implicit: reduces sub-dt until convergence is achieved.

def sim( self, dt: float, nSteps: int, S0: numpy.ndarray, P0: numpy.ndarray | None = None, pbar: bool = True, leave: bool = True, **kwargs) -> tuple: View Source

837    def sim(
838        self,
839        dt: float,
840        nSteps: int,
841        S0: np.ndarray,
842        P0: np.ndarray | None = None,
843        pbar: bool = True,
844        leave: bool = True,
845        **kwargs,
846    ) -> tuple:
847        """Recursively (`nSteps` times) apply `time_stepper` with `dt`, from `S0`.
848
849        Returns the saturation and pressure trajectories, `(SS, PP)`.
850
851        .. note:: `SS[0] == S0` and `PP[0] == P0`, hence both have `len = nSteps + 1`.
852
853            `P0` defaults to zeros. It is only consequential if `ct > 0`.
854        """
855        step = self.time_stepper(dt, **kwargs)
856
857        # pbar
858        kk = np.arange(nSteps)
859        if pbar:
860            kk = tqdm(kk, "Simulation", leave=leave, mininterval=1e-2)
861
862        # Allocate
863        SS = np.zeros((nSteps + 1,) + S0.shape)
864        PP = np.zeros((nSteps + 1, self.Nxy))
865        self.wells.actual_rates = np.zeros((self.nComp, nSteps))
866        self.wells.actual_bhp = np.full((self.nComp, nSteps), np.nan)
867
868        # Init
869        SS[0] = S0
870        if P0 is not None:
871            PP[0] = P0
872
873        # Recurse
874        for k in kk:
875            SS[k + 1], PP[k + 1] = step(SS[k], PP[k], k)
876
877        return SS, PP

Recursively (nSteps times) apply time_stepper with dt, from S0.

Returns the saturation and pressure trajectories, (SS, PP).

`SS[0] == S0` and `PP[0] == P0`, hence both have `len = nSteps + 1`.

P0 defaults to zeros. It is only consequential if ct > 0.

Inherited Members

minires.grid.Grid2D: Lx; Ly; Nx; Ny; shape; size; domain; Nxy; hx; hy; h2; mesh; sub2ind; ind2sub; xy2sub; xy2ind; sub2xy; ind2xy
minires.plotting.Plot2D: plt_field; plt_faces; well_scatter; plt_production; anim

minires

Governing equations

Derivation

Single phase

Two phases

How to solve

The pressure system need not be factorized afresh each step.

Units

Compressibility

Definition

Approximation

Derivation

Consequences

The rate unit is forced to $u_L^2/u_t$ -- an areal rate.

$C$ enters at exactly 2 sites, both of them Darcy's law.

With ct == 0 the rates must still sum to 0 at every step.

The mode switch lags the solve by one step.

S and P may be None, so an override should tolerate that.

Setting both controls for a well is not an error, just pointless.

$ λ_t $ is that of the well's cell.

The Newton iteration can converge to a spurious root.

This scheme rarely earns its keep.

SS[0] == S0 and PP[0] == P0, hence both have len = nSteps + 1.

Inherited Members

With `ct == 0` the rates must still sum to 0 at every step.

`S` and `P` may be `None`, so an override should tolerate that.

`SS[0] == S0` and `PP[0] == P0`, hence both have `len = nSteps + 1`.