TPFA_ResSim

TPFA-ResSim is a 2D, two-phase, black-oil, immiscible reservoir simulator using TPFA (two-point flux approximation).

One waterflood: permeability, pressure, water front, and the adjoint sensitivity

See examples for more demonstrations.

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, quadratic (Corey) curves of the saturation rescaled by its residual values (ResSim.RelPerm). 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 TPFA_ResSim.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 TPFA_ResSim.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 "Vocabulary"), 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 TPFA_ResSim.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.

Vocabulary of reservoir engineering

Reservoir simulators implement porous media flow on upscaled geophysical parameters typically with grid blocks between 1 - 100 m. They usually parameterize multiphase flow. If only the two phases of oil and water are used it is called black-oil. A common assumption is that the flow is immiscible: not mixing (oil and water). But this does not mean that gas cannot be dissolved in oil.

Fossil fuel hydrocarbons is sedimented, pressurized, organic material (mostly plants?) that used to live on the sub-sea continental shelves On-land organic material turns into coal. ⇒ Saudi-Arabia used to be sub-sea? The energy in oil & gas comes from the sun (photosynthesis), not the compression.

The lightest hydrocarbons (methane, ethane, etc.) usually escapes quickly, while oils moves slowly towards the surface. Sometimes the geology is bends to form caps of non-permeable rock, so that the migrating hydrocarbons are trapped. 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. In the North Sea, these reservoirs lie 1000-3000 meters below the sea bed. Norway is also surrounded by the Norwegian sea, and the Barents sea, towards Murmansk.

Porosity, $\phi$, is the void volume fraction. Depends on pressure, because rock is compressible. Compressibility is the porosity's (relative) gradient wrt. pressure. Usually neglected, so that $\phi$ is a constant, but spatial, field.

Permeability, denoted by tensor $\mathbf{K}$, quantifies transmissibility. Usually SPD, and correlated with $\phi$. 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$). A medium is called isotropic if $\mathbf{K}$ is scalar.

The phases (rock, oil, gas), whose saturations sum to $1$, contains components (e.g. methane, ethane, propane), usually grouped as pseudo-components. Each phase's mass fraction component, $c_{phase,i}$, sums to $1$. Each phase has density, $\rho$ and viscosity, $\mu$, generally functions of the phase pressure, but usually neglected except for gas. The differences in pressure are named capillary pressure because they arise due to interfacial tensions. A phase's compressibility is defined similar as for the rock's. Confusingly, it is also denoted with $c$, but using only a single subscript.

Phases do not really mix. But in macro-scale modelling all phases may be present at the same location. Therefore a phase's permeability should depend on the saturations, to which end we introduce relative permeability, $k_{r,i} = k_{r,i}(s_g, s_o), i = g, o, w$ a nonlinear function, yielding an (effective) permeability $\mathbf{K_i} = \mathbf{K} k_{r,i}$ Relative permeability curves do not extend all over the interval $[0, 1]$. The smallest saturation where a phase is mobile is called the residual saturation. This adsorption effects may vary, and this may have important effects, particularly for simulation of polymer injection. The uncertainty regarding relative permeability is modest compared to the enormous uncertainty of the rock permeability.

Everything depends on thermodynamics, but this is often complex and neglected, except perhaps for the bubble/boiling point pressures, which govern how much of the gas dissolves in oil.

Since compressibility relates volumes to pressure, a volume must be qualified by where it is measured. The formation volume factor, $B$, is the ratio of the volume at reservoir conditions to that of the same mass at the surface ("stock tank"), and is how field rates (measured at the surface) are converted to the reservoir rates that a simulator works in. This model has $B = 1$. Related PVT (pressure-volume-temperature) vocabulary: the bubble point is the pressure below which gas comes out of solution; an oil above it is undersaturated, and the amount of gas it holds is the solution gas-oil ratio, $R_s$.

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), gravity drainage, and compaction drive (which manifests as seabed subsidence). Recovery is staged: primary production runs on the native drive; secondary 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.

Aquifers are beneficial in reservoirs as they act as pressure compensators. Oil production ⇒ pressure decrease ⇒ aquifers expansion ⇒ pressure compensation. Despite consisting of water, the expansion is generally significant because the base volume is so big, or the aquifer might even be connected to the ocean.

Other lingo: water table, facies, channels, fissures, fractures.

View Source

 1""".. include:: README.md"""
 2
 3from TPFA_ResSim.core import ResSim
 4from TPFA_ResSim.grid import Fluxes, Grid2D
 5from TPFA_ResSim.wells import Wells, peaceman_WI, well_path
 6
 7# Also pdoc's table of contents: `ResSim` is documented on the package page (beside the
 8# README), and the listed submodules on their own pages. `core` is deliberately absent,
 9# lest `ResSim` be documented twice; so are the other re-exports, which are documented
10# in their home modules.
11__all__ = ["ResSim", "grid", "wells", "plotting", "tlm"]

@dataclass

class ResSim(TPFA_ResSim._repr.AlignedRepr, TPFA_ResSim.grid.Grid2D, TPFA_ResSim.plotting.Plot2D): View Source

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

vw: float = 1.0

Viscosity for water.

vo: float = 1.0

Viscosity for oil.

swc: float = 0.0

Irreducible saturation, water.

sor: float = 0.0

Irreducible saturation, oil.

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, fidelity ($ c_t \, Δp \ll 1 $, which the voidage sets, not ct) and vocabulary: ref the "Compressibility" section of the docs.

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)).

por: Any = None

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

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 TPFA_ResSim.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.

nComp View Source

183    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 TPFA_ResSim.wells.Wells.nComp.

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

190    def assemble_wells(
191        self, S: np.ndarray | None, P: np.ndarray | None, k: int
192    ) -> None:
193        """Set up (for time `k`) the wells' contributions to the equations.
194
195        The controls are those of `well_controls`, to which `S` and `P` (the
196        state at the *start* of the step) are simply passed on.
197        Rate-controlled wells enter the source/sink *field*, `_Q`, directly.
198        BHP-controlled ones (ref `TPFA_ResSim.wells.Wells.bhp`) cannot: their
199        rate is not yet known. They instead enter the pressure equations in
200        `TPFA`, after which `realize_bhp` folds the resulting rate into `_Q`.
201        """
202        ctrl = self.well_controls(S, P, k)
203        inds = self.xy2ind(*self.wells.xy.T)
204        rates, p_bh = ctrl["rates"], ctrl["bhp"]
205        is_bhp = np.isfinite(p_bh)
206        assert np.isfinite(rates[~is_bhp]).all(), (
207            "A rate-controlled well has a non-finite rate. Give it a number"
208            " (`0` shuts it in), or put it on BHP control; ref `Wells.rates`."
209        )
210
211        # The well model's constant of proportionality, WI * λ_t.
212        # NB: `nan` marks the rate-controlled wells, throughout.
213        WI_lam = np.full(self.nComp, np.nan)
214        if is_bhp.any():
215            WI = self.wells.WI
216            assert WI is not None and np.isfinite(WI[is_bhp]).all(), (
217                "BHP control requires (finite) `Wells.WI`."
218            )
219            assert S is not None, "BHP control requires `S` (for λ_t)."
220            Mw, Mo = self.RelPerm(S)
221            WI_lam[is_bhp] = WI[is_bhp] * (Mw + Mo)[inds[is_bhp]]
222
223        # Translate well conditions for cells.
224        # NB: Dont use `Q[inds] += ...` since `inds` may contain dupes.
225        self._Q, bhp_diag, bhp_rhs = np.zeros((3, self.Nxy))
226        np.add.at(self._Q, inds[~is_bhp], rates[~is_bhp])
227        np.add.at(bhp_diag, inds[is_bhp], WI_lam[is_bhp])
228        np.add.at(bhp_rhs, inds[is_bhp], (WI_lam * p_bh)[is_bhp])
229        rates[is_bhp] = np.nan  # only `realize_bhp` knows these
230        self._wells_now: dict[str, np.ndarray] = dict(
231            inds=inds, rates=rates, p_bh=p_bh,
232            WI_lam=WI_lam, bhp_diag=bhp_diag, bhp_rhs=bhp_rhs,
233        )  # 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 TPFA_ResSim.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

235    def realize_bhp(self, P: np.ndarray) -> None:
236        """Compute rates for BHP wells. Enter into `_Q` and `_wells_now["rates"]`.
237
238        The rate, $ WI λ_t (p_\\mathrm{bh} - p_\\mathrm{cell}) $, is signed by
239        nature (ref the `TPFA_ResSim.wells.Wells.bhp` warning).
240
241        By construction of the linear system of `TPFA`, this leaves `_Q` equal
242        to the *total* well flux, which is what keeps `storage_rate` -- and
243        hence the transport step -- consistent with the pressure solution.
244        """
245        wls = self._wells_now
246        WI_lam = wls["WI_lam"]  # `nan` marks the rate-controlled wells
247        # Insert in cell source/sink field
248        self._Q = self._Q + wls["bhp_rhs"] - wls["bhp_diag"] * P
249        # Insert in per-well rates
250        is_bhp = np.isfinite(WI_lam)
251        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 TPFA_ResSim.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

271    def well_controls(self, S: np.ndarray | None, P: np.ndarray | None, k: int) -> dict:
272        """Compute the wells' controls for time `k`: `dict(rates=..., bhp=...)`.
273
274        Each is a `(nComp,)` array, read off the specifications --
275        `TPFA_ResSim.wells.Wells.rates`, `TPFA_ResSim.wells.Wells.bhp` -- which
276        are *open-loop*: fixed before the simulation begins. Overriding
277        (patching/subclassing) this method is therefore how to do *feedback*
278        control, the controls being free to depend on the state at the *start*
279        of the step: the saturation `S` and the pressure `P`.
280        The returned arrays are copies, so they may be modified in place.
281
282        Most feedback concerns the rates alone -- e.g. shutting the wells upon
283        water breakthrough at the producer:
284
285        >>> class Shutter(ResSim):
286        ...     def well_controls(self, S, P, k):
287        ...         ctrl = super().well_controls(S, P, k)
288        ...         if S is not None and S[self.xy2ind(1, 1)] > .5:
289        ...             ctrl["rates"][:] = 0    # NB: all of them! See warning
290        ...         return ctrl
291        >>> model = Shutter(Lx=1, Ly=1, Nx=16, Ny=16,
292        ...                 wells=Wells(xy=[[0, 0], [1, 1]], rates=[[1], [-1]]))
293        >>> SS, PP = model.sim(.05, 20, model.swc*np.ones(model.Nxy), pbar=False)
294        >>> int((model.wells.actual_rates[1] == 0).argmax())  # step of breakthrough
295        16
296
297        But the `bhp` is here too, and with it each well's *control mode*
298        (`nan` => rate-controlled, ref `TPFA_ResSim.wells.Wells.bhp`) -- which
299        is what an approximate mode *switch* requires. For example, rate
300        control with a BHP limit -- the industrial default -- wherein a
301        producer holds its rate only for as long as that does not draw it
302        below some `p_min`:
303
304        >>> class Limited(ResSim):
305        ...     p_min = .5
306        ...     def well_controls(self, S, P, k):
307        ...         ctrl = super().well_controls(S, P, k)
308        ...         if P is None:
309        ...             return ctrl                    # nothing to switch on
310        ...         p_bh = self.bhp(S, P, ctrl["rates"])
311        ...         switch = p_bh < self.p_min         # the rate is unsustainable
312        ...         ctrl["bhp"] = np.where(switch, self.p_min, np.nan)
313        ...         return ctrl
314        >>> model = Limited(Lx=1, Ly=1, Nx=16, Ny=16, ct=.1,
315        ...                 wells=Wells(xy=[[.5, .5]], rates=[[-.25]]))
316        >>> from TPFA_ResSim import peaceman_WI
317        >>> model.wells.WI = peaceman_WI(model, model.wells.xy, rw=1e-3)
318        >>> SS, PP = model.sim(.02, 25, np.zeros(model.Nxy),
319        ...                    P0=np.ones(model.Nxy), pbar=False)
320
321        The well delivers its target rate until the limit binds, and declines
322        thereafter -- at constant $ p_\\mathrm{bh} $, exponentially so
323        (ref `examples/well_control.py`, which plots all three modes):
324
325        >>> (-model.wells.actual_rates[0, [0, 5, 6, -1]]).round(3)
326        array([0.25 , 0.25 , 0.182, 0.005])
327
328        .. warning:: With `ct == 0` the rates must still sum to 0 at every step.
329
330            Ref `TPFA_ResSim.wells.Wells.rates`. So shutting one well requires
331            matching it on the other side -- as above.
332
333        .. note:: The mode switch lags the solve by one step.
334
335            It is decided from the previous step's pressure, whereas the well
336            model itself is solved *simultaneously* with the new one. So the
337            limit is breached for the one step in which it comes to bind.
338            Shorten `dt` to refine.
339
340        .. note:: `S` and `P` may be `None`, so an override should tolerate that.
341
342            They are `None` if the caller has none to offer -- as when
343            `assemble_wells` is used merely to set up a plot.
344
345        .. note:: Setting both controls for a well is not an error, just pointless.
346
347            `assemble_wells` discards the rate of a BHP-controlled well -- it
348            is `realize_bhp` that fills it in.
349        """
350        return dict(
351            rates=self.wells.at_time("rates", 0.0, k),
352            bhp=self.wells.at_time("bhp", np.nan, k),
353        )

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

Each is a (nComp,) array, read off the specifications -- TPFA_ResSim.wells.Wells.rates, TPFA_ResSim.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.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 TPFA_ResSim.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 TPFA_ResSim 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 TPFA_ResSim.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

355    def bhp(self, S: np.ndarray, P: np.ndarray, rates: np.ndarray) -> np.ndarray:
356        """Bottom-hole pressures implied by the (signed) `rates`, via the well indices.
357
358        I.e. the well model of `TPFA_ResSim.wells.Wells.WI`, solved for
359        $ p_\\mathrm{bh} $:
360        the rate's sign puts an injector above, a producer below, its cell
361        pressure. `nan` wherever the well index is unset.
362
363        `S` and `P` (both flat) should be the saturation and the pressure of the
364        *same* `pressure_step`, i.e. `SS[k]` and `PP[k+1]` of `sim` -- which is
365        what `actual_bhp` records, so prefer reading that.
366
367        .. warning:: $ λ_t $ is that of the well's *cell*.
368
369            So an injector's injectivity is governed by the mobility of
370            whatever the cell currently holds, rather than by that of the
371            injectant.
372        """
373        if self.wells.WI is None:
374            return np.full(self.nComp, np.nan)
375        Mw, Mo = self.RelPerm(S)
376        ii = self.xy2ind(*self.wells.xy.T)
377        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 TPFA_ResSim.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, TPFA_ResSim.grid.Fluxes]: View Source

380    def pressure_step(
381        self,
382        S: np.ndarray,
383        P: np.ndarray | None = None,
384        dt: float | None = None,
385    ) -> tuple[np.ndarray, Fluxes]:
386        """Compute permeabilities then solve Darcy's equation. Returns `[P, V]`.
387
388        `P` (flat, like `S`) is the *previous* step's pressure: used (and
389        required) only if `ct > 0`, along with `dt`. The new one replaces it.
390        """
391        # Compute K*λ(S)
392        Mw, Mo = self.RelPerm(S)
393        Mt = Mw + Mo
394        Mt = Mt.reshape(self.shape)
395        KM = Mt * self.K
396        # Compute pressure and extract fluxes
397        [P, V] = self.TPFA(KM, P, dt)
398        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 rescale_sat(self, s: numpy.ndarray) -> numpy.ndarray: View Source

416    def rescale_sat(self, s: np.ndarray) -> np.ndarray:
417        """Account for irreducible saturations. Ref paper, p. 32."""
418        return (s - self.swc) / (1 - self.swc - self.sor)

Account for irreducible saturations. Ref paper, p. 32.

def RelPerm(self, s: numpy.ndarray) -> tuple: View Source

421    def RelPerm(self, s: np.ndarray) -> tuple:
422        """Rel. permeabilities of oil and water. Return as mobilities (perm/viscocity)."""
423        S = self.rescale_sat(s)
424        Mw = S**2 / self.vw  # Water mobility
425        Mo = (1 - S) ** 2 / self.vo  # Oil mobility
426        return Mw, Mo

Rel. permeabilities of oil and water. Return as mobilities (perm/viscocity).

def dRelPerm(self, s: numpy.ndarray) -> tuple: View Source

428    def dRelPerm(self, s: np.ndarray) -> tuple:
429        """Derivatives of `RelPerm`."""
430        S = self.rescale_sat(s)
431        dMw = 2 * S / self.vw / (1 - self.swc - self.sor)
432        dMo = -2 * (1 - S) / self.vo / (1 - self.swc - self.sor)
433        return dMw, dMo

Derivatives of RelPerm.

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

436    def TPFA(
437        self,
438        K: np.ndarray,
439        P: np.ndarray | None = None,
440        dt: float | None = None,
441    ) -> tuple[np.ndarray, Fluxes]:
442        """Two-point flux-approximation (TPFA) of Darcy: $ -∇(K ∇u) = q $
443
444        i.e. steady-state diffusion w/ nonlinear coefficient, $K$,
445        if `ct == 0`. Otherwise (slightly compressible model) solve
446        the backward-Euler step of $ φ c_t ∂u/∂t - ∇(K ∇u) = q $,
447        which requires the previous pressure, `P`, and `dt`.
448
449        After solving for pressure `P`, extract the fluxes `V`
450        by finite differences.
451        """
452        # Compute transmissibilities by harmonic averaging.
453        C = self.cdarcy
454        L = 1 / K
455        TX = np.zeros((self.Nx + 1, self.Ny))
456        TY = np.zeros((self.Nx, self.Ny + 1))
457        TX[1:-1, :] = C * 2 * self.hy / self.hx / (L[0, :-1, :] + L[0, 1:, :])
458        TY[:, 1:-1] = C * 2 * self.hx / self.hy / (L[1, :, :-1] + L[1, :, 1:])
459
460        # Assemble TPFA discretization matrix.
461        x1 = TX[:-1, :].ravel()
462        x2 = TX[1:, :].ravel()
463        y1 = TY[:, :-1].ravel()
464        y2 = TY[:, 1:].ravel()
465
466        # Setup linear system
467        DiagVecs = [-x2, -y2, y1 + y2 + x1 + x2, -y1, -x1]
468        DiagIndx = [-self.Ny, -1, 0, 1, self.Ny]
469        q = self._Q
470        if self.ct > 0:
471            # Accumulation term (φ ct h²/dt) of backward Euler.
472            # Renders the system nonsingular (unlike the pure-Neumann problem).
473            assert P is not None and dt is not None, (
474                "Compressible model (ct > 0) requires the previous P, and dt."
475            )
476            accum = self.por.ravel() * self.ct * self.h2 / dt
477            DiagVecs[2] = DiagVecs[2] + accum
478            q = q + accum * P
479        elif not self._wells_now["bhp_diag"].any():
480            # Pin the (o/w pure-Neumann & singular) problem.
481            DiagVecs[2][0] += np.sum(self.K[:, 0, 0])  # ref article p. 13
482        # Well model of the BHP-controlled wells
483        DiagVecs[2] = DiagVecs[2] + self._wells_now["bhp_diag"]
484        q = q + self._wells_now["bhp_rhs"]
485
486        # Solve; compute A\q to update P
487        A = self._spdiags(DiagVecs, DiagIndx)
488        P = self._solve_pressure(A.tocsr(), q, P)
489        # P = np.linalg.solve(A.A, q) # direct dense solver
490        # Could also try scipy.linalg.solveh_banded which, according to
491        # https://scicomp.stackexchange.com/a/30074 uses the Thomas algorithm,
492        # as recommended by Aziz and Settari ("Petro. Res. simulation").
493        # NB: stackexchange also mentions that solve_banded does not work well
494        # when the band offsets large, i.e. higher-dimensional problems.
495
496        # Extract fluxes, via a grid-shaped view of the (flat) pressure.
497        P2d = P.reshape(self.shape)
498        V = Fluxes(
499            x=np.zeros((self.Nx + 1, self.Ny)),
500            y=np.zeros((self.Nx, self.Ny + 1)),
501        )
502        V.x[1:-1, :] = (P2d[:-1, :] - P2d[1:, :]) * TX[1:-1, :]
503        V.y[:, 1:-1] = (P2d[:, :-1] - P2d[:, 1:]) * TY[:, 1:-1]
504        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: TPFA_ResSim.grid.Fluxes) -> scipy.sparse._dia.dia_matrix: View Source

537    def upwind_diff(self, V: Fluxes) -> sparse.dia_matrix:
538        """Upwind finite-volume scheme."""
539        fp = self._Q.clip(max=0)  # production
540        # Flow fluxes, separated into direction (x-y) and sign
541        x1 = V.x.clip(max=0)[:-1, :].ravel()
542        y1 = V.y.clip(max=0)[:, :-1].ravel()
543        x2 = V.x.clip(min=0)[1:, :].ravel()
544        y2 = V.y.clip(min=0)[:, 1:].ravel()
545        DiagVecs = [x2, y2, fp + y1 - y2 + x1 - x2, -y1, -x1]
546        DiagIndx = [-self.Ny, -1, 0, 1, self.Ny]
547        A = self._spdiags(DiagVecs, DiagIndx)
548        return A

Upwind finite-volume scheme.

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

550    def storage_rate(self, V: Fluxes) -> np.ndarray:
551        """The volume rate, per cell, that goes into storage: $ q - ∇ ⋅ V $.
552
553        For the incompressible model this is `0`: the fluxes balance the wells
554        exactly, cell by cell. With `ct > 0` it is (by construction of the
555        linear system of `TPFA`) the accumulation term of the backward-Euler
556        step, $ φ \\, c_t \\, h^2 \\, (p^{n+1} - p^n) / Δt $, which the saturation
557        steps charge to the phases (ref `ct`).
558
559        Computing it from `V` (rather than from $p^{n+1} - p^n$) means it is
560        *exactly* the imbalance seen by the transport scheme, whose `upwind_diff`
561        is assembled from the same fluxes.
562        """
563        if self.ct == 0:
564            return np.zeros(self.Nxy)
565        divV = (V.x[1:, :] - V.x[:-1, :]) + (V.y[:, 1:] - V.y[:, :-1])
566        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: TPFA_ResSim.grid.Fluxes, fi: numpy.ndarray) -> float: View Source

569    def estimate_1CFL(self, pv: np.ndarray, V: Fluxes, fi: np.ndarray) -> float:
570        """Estimate 1/CFL for use with `saturation_step_upwind`."""
571        # In-/Out-flux x-/y- faces
572        XP = V.x.clip(min=0)
573        XN = V.x.clip(max=0)
574        YP = V.y.clip(min=0)
575        YN = V.y.clip(max=0)
576        Vi = XP[:-1, :] + YP[:, :-1] - XN[1:, :] - YN[:, 1:]
577
578        flx = max((Vi.ravel() + fi) / pv)  # estimate of influx
579        # NB: `storage_rate` is not counted here. In practice it is a small
580        # fraction of the fluxes that are (under 20% even at `ct = 10`),
581        # so the safety factor below covers it.
582        sat = self.swc + self.sor
583        cfl = 3 / (1 - sat) * flx  # NB: 3-->2 since no z-dim ?
584        # NB: the ceiling is nudged down by a relative epsilon, so that a `dt`
585        # sitting *on* an integer multiple of the CFL limit (as the examples'
586        # round numbers tend to) does not gain a whole extra sub-step from the
587        # last bits of the linear solve -- which differ across platforms and
588        # library versions, and would make the results irreproducible.
589        # The CFL estimate carries a safety factor of 3, so shaving 1e-9 off it
590        # cannot cost stability.
591        return cfl * (1 - 1e-9)

Estimate 1/CFL for use with saturation_step_upwind.

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

594    def saturation_step_upwind(self, S: np.ndarray, V: Fluxes, dt: float) -> np.ndarray:
595        """Explicit upwind FV discretisation of conserv. of mass (water sat.)."""
596        # fmt: off
597        A  = self.upwind_diff(V)                 # FV discretized transport operator
598        pv = self.h2 * self.por.ravel()          # Pore volume (per thickness)
599        fi = self._Q.clip(min=0)                 # Well inflow
600        st = self.storage_rate(V)                # Storage (0 if incompressible)
601
602        # Compute sub/local dt
603        cfl1 = self.estimate_1CFL(pv, V, fi)
604        nT = int(np.ceil(dt * cfl1))
605        nT = max(1, nT)
606
607        # Scale A
608        dtx = dt / nT / pv                       # timestep / pore volume
609        B   = self._spdiags(dtx, 0) @ A          # A * dt/|Omega i|
610
611        for _ in range(nT):
612            Mw, Mo = self.RelPerm(S)             # compute mobilities
613            fw = Mw / (Mw + Mo)                  # compute fractional flow
614            S = S + (B@fw + (fi - S*st)*dtx)     # update saturation
615        # fmt: on
616        return S

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

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

619    def saturation_step_implicit(
620        self,
621        S: np.ndarray,
622        V: Fluxes,
623        dt: float,
624        nNewtonMax: int = 10,
625        nTmax_log2: int = 10,
626    ) -> np.ndarray:
627        """Implicit FV discretisation of conserv. of mass (water sat.).
628
629        .. warning:: The Newton iteration can converge to a spurious root.
630
631            Far outside the $ c_t \\, Δp \\ll 1 $ regime (ref `ct`), it may
632            converge -- silently -- to a root of the residual outside $[0, 1]$:
633            the polynomial `RelPerm` extends smoothly beyond the unit interval,
634            and the sub-`dt` halving only triggers on *non*-convergence.
635            The explicit scheme
636            (`saturation_step_upwind`), being monotone, stays within $[0, 1]$
637            even for extreme `ct`.
638
639        .. note:: This scheme rarely earns its keep.
640
641            It is usually both slower & less accurate than `saturation_step_upwind`.
642            Both schemes sub-divide `dt` internally, so it is not `dt` that decides
643            their cost but the stiffness of the grid -- and the well cells, being
644            normally the stiffest, hold the two requirements within a factor 3 of each
645            other (that being the safety margin of `ResSim.estimate_1CFL`), while an
646            implicit sub-step -- a sparse solve, or several -- 10 or 100x more.
647            It does pay off where the stiffest cell is *not* a well cell -- a tight
648            streak, a fracture, a locally refined region -- running 7 times faster at a
649            1000x pore-volume contrast. The branch `implicit-transport-scheme` says more.
650        """
651        # fmt: off
652        A  = self.upwind_diff(V)                 # FV discretized transport operator
653        pv = self.h2 * self.por.ravel()          # Pore volume (per thickness)
654        fi = self._Q.clip(min=0)                 # Well inflow
655        st = self.storage_rate(V)                # Storage (0 if incompressible)
656
657        # For each iter, halve the sub/local dt
658        for nT_log2 in range(0, nTmax_log2):
659            nT = 2**nT_log2
660
661            # Scale A
662            dtx = dt / nT / pv                   # timestep / pore volume
663            B   = self._spdiags(dtx, 0) @ A      # A * dt/|Omega i|
664            C   = self._spdiags(dtx*st, 0)       # storage, likewise scaled
665
666            Sn = S
667            for _ in range(nT):
668                Sp = Sn
669                for _ in range(nNewtonMax):
670                    Mw, Mo   = self.RelPerm(Sn)    # mobilities
671                    dMw, dMo = self.dRelPerm(Sn)   # their derivatives
672                    df = dMw/(Mw+Mo) - Mw/(Mw+Mo)**2 * (dMw + dMo)        # df w/ds
673                    dG = (sparse.eye(self.Nxy) + C                        # deriv of G
674                          - B @ self._spdiags(df, 0))
675
676                    fw = Mw / (Mw+Mo)               # fract. flow
677                    G  = Sn - Sp - (B@fw + (fi - Sn*st)*dtx)  # G(s)
678                    dS = spsolve(dG, G)             # compute dS
679                    Sn = Sn - dS                    # update S
680
681                    if np.sqrt(sum(dS**2)) < 1e-3:
682                        # If converged: halt Newton iterations
683                        break
684                else:
685                    # If never converged: increase nT, restart time loop
686                    break
687            else:
688                # If completed all time steps, halt
689                break
690        else:
691            # Failed (even with max nT) to complete all time steps
692            print("Warning: did not converge")
693        # fmt: on
694
695        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 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

707    def time_stepper(self, dt: float, implicit: bool = False) -> Callable:
708        """Get ODE solver (integrator) for model.
709
710        Whatever time step `dt` is given, both schemes will use smaller steps internally.
711
712        - `explicit`: computes sub-`dt` based on CFL esitmate.
713        - `implicit`: reduces sub-`dt` until convergence is achieved.
714        """
715
716        def integrate(S, P, k):
717            self.assemble_wells(S, P, k)
718            self._validate()
719            [P, V] = self.pressure_step(S, P, dt)
720            self.realize_bhp(P)
721            self._record_actual_well_operation(S, P, k)
722            if implicit:
723                S = self.saturation_step_implicit(S, V, dt)
724            else:
725                S = self.saturation_step_upwind(S, V, dt)
726            return S, P
727
728        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

730    def sim(
731        self,
732        dt: float,
733        nSteps: int,
734        S0: np.ndarray,
735        P0: np.ndarray | None = None,
736        pbar: bool = True,
737        leave: bool = True,
738        **kwargs,
739    ) -> tuple:
740        """Recursively (`nSteps` times) apply `time_stepper` with `dt`, from `S0`.
741
742        Returns the saturation and pressure trajectories, `(SS, PP)`.
743
744        .. note:: `SS[0] == S0` and `PP[0] == P0`, hence both have `len = nSteps + 1`.
745
746            `P0` defaults to zeros. It is only consequential if `ct > 0`.
747        """
748        step = self.time_stepper(dt, **kwargs)
749
750        # pbar
751        kk = np.arange(nSteps)
752        if pbar:
753            kk = tqdm(kk, "Simulation", leave=leave, mininterval=1e-2)
754
755        # Allocate
756        SS = np.zeros((nSteps + 1,) + S0.shape)
757        PP = np.zeros((nSteps + 1, self.Nxy))
758        self.wells.actual_rates = np.zeros((self.nComp, nSteps))
759        self.wells.actual_bhp = np.full((self.nComp, nSteps), np.nan)
760
761        # Init
762        SS[0] = S0
763        if P0 is not None:
764            PP[0] = P0
765
766        # Recurse
767        for k in kk:
768            SS[k + 1], PP[k + 1] = step(SS[k], PP[k], k)
769
770        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

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

TPFA_ResSim

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

Vocabulary of reservoir engineering

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`.