TPFA_ResSim
As mentioned in the the main README this is a 2D, two-phase, black-oil, immiscible, incompressible reservoir simulator, neglecting capillary forces and gravity, using TPFA (two-point flux approximation), equipped with explicit and implicit ode solvers.
Governing equations
The simulator solves eqn. (1) and (2) (corresponding to (42) and (43) of the reference paper) :
$$- \nabla \cdot \mathbf{K} \lambda(s) \, \nabla p = q \,, \tag{1}$$ $$\; \phi \frac{\partial s}{\partial t} + \nabla \cdot (f(s)\, \mathbf{v}) = \frac{q_w}{\rho_w} \,. \tag{2}$$
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
- $\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}}$.
- $f(s) = \lambda_w / \lambda \in [0, 1]$ is the water fractional flow, where both mobilities depend on $s$. It gives $\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.
Relative permeabilities, $k_{\text{phase}} \in [0, 1]$,
are set via a constituent relation that is a function of the (reducible) saturation. They do not generally sum to 1. Their approximation and uncertainty is significant, but usually less important than those of the absolute permeability, $\mathbf{K}$.
Derivation
Single phase
Conservation of mass in a (single-phase) fluid setting is usually derived along with the concept of "material derivative", i.e. $\frac{D}{D t}$. This would be a bit trickier in the case of a porous medium, but the same concepts lead to
$$\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$). Notice that it "simply" says that divergence/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}$$
This still leaves us with one equation for 3 unknowns. The system is closed by Darcy's law: $$\mathbf{v} = − \frac{\mathbf{K}}{\mu} \nabla u \,, \tag{5}$$ where $$u = p - \rho g z $$ is the velocity potential, Darcy's law provides us with 3 additional equations and 1 additional unknown, $p$. It is analogous to Fourier's heat diffusion and Ohm's conduction law, but contains two forces (pressure gradient and gravity).
Darcy's law (5) may be derived from Navier-Stokes momentum equations, but was empirically derived by Darcy. Indeed, it is simply a formula for the velocity, as the gradient of the 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. Still, $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
The reference paper explains how to apply the continuity equation (3) and Darcy's law (5) for each phase in a multiphase (and even multicomponent) flow system. Even the black-oil case involves 27 unknowns and equations. By assuming immiscibility and incompressibility, and constituent relations, and astute combination of the equations, they arrive at eqn. (1) and (2).
Here we shall be more heuristic but brief.
- Incompressibility again yields eqn. (4) for the total (volumetric) velocity.
- Darcy's law (5) is assumed for each (both) individual phase, meaning that $\mathbf{K}$ is 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}$$
- Hence, repeating the derivation for eqn. (6), we obtain eqn. (1).
Meanwhile, conservation of mass (3) for a single, incompressible phase is obtained by replacing the density $\rho$ in eqn. (3) by $s_{\text{phase}} \, \rho_{\text{phase}}$, and $\mathbf{v}$ by $\mathbf{v}_{\text{phase}} = f_{\text{phase}}(s)\, \mathbf{v}$. This immediately yields 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. Given this complication, it is then possible to use implicit time discretization (like ECLIPSE 100) where $s_{t+1}$ is expressed as a (nonlinear) function of itself, which also requires iterative solution.
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. The simulator code contains both an implicit and explicit (upwind) time discretization for the nearly-hyperbolic saturation equation. When using the explicit one, the strategy is called IMPES (implicit pressure, explicit saturation), although I'm not sure why, since the pressure equation itself does not contain a time derivative (with $s$ fixed, equation (1) is a nearly-elliptic boundary value problem for the pressure, $p$).
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. 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.
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.
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, water cut, fissures, fractures.
1""".. include:: README.md""" 2 3from dataclasses import dataclass 4import numpy as np 5from scipy import sparse 6from scipy.sparse.linalg import spsolve 7from struct_tools import DotDict, NicePrint 8from tqdm.auto import tqdm 9 10from TPFA_ResSim.grid import Grid2D 11from TPFA_ResSim.plotting import Plot2D 12 13 14@dataclass 15class ResSim(NicePrint, Grid2D, Plot2D): 16 """Reservoir simulator class. 17 18 Implemented with OOP (instead of passing around dicts) to facilitate 19 bookkeeping of ensemble forecasting 20 (where parameter values of one instance should not influence another) 21 22 Example: 23 >>> model = ResSim(Lx=1, Ly=1, Nx=64, Ny=64) 24 >>> model.inj_xy=[[0, .32]] 25 >>> model.prd_xy=[[1, 1]] 26 >>> model.inj_rates=[[1]] 27 >>> model.prd_rates=[[1]] 28 >>> water_sat0 = np.zeros(model.Nxy) 29 >>> dt = .35 30 >>> nSteps = 2 31 >>> S = model.sim(dt, nSteps, water_sat0, pbar=False) 32 33 This produces the following values (used for automatic testing): 34 >>> S[-1, [100, 1300, 2900]] 35 array([0.9429345 , 0.91358172, 0.71554613]) 36 """ 37 # Dont use dataclass repr 38 __repr__ = NicePrint.__repr__ 39 __str__ = NicePrint.__str__ 40 41 def __post_init__(self): 42 defaults = dict(K=np.ones((2, *self.shape)), 43 por=np.ones(self.shape)) 44 for k, v in defaults.items(): 45 if getattr(self, k) is None: 46 setattr(self, k, v) 47 48 # Prefer __setattr__ approach (over @property get/set-ers) 49 # because @property requires the _private pattern, 50 # which is pretty ugly with dataclasses, 51 # and also because can unify treatment of inj/prd wells. 52 def __setattr__(self, key, val): 53 if val is not None: 54 # Well positions -- collocate at some node 55 if key in ["inj_xy", "prd_xy"]: 56 val = np.array(val, float).reshape((-1, 2)) 57 for i, (x, y) in enumerate(val): 58 val[i] = self.ind2xy(self.xy2ind(x, y)) 59 # Well rates 60 if key in ["inj_rates", "prd_rates"]: 61 nWell = len(getattr(self, key.replace("rates", "xy"))) 62 val = np.array(val, float).reshape((nWell, -1)) 63 # Permeabilities 64 if key == "K": 65 if np.isscalar(val): 66 val = np.full_like(self.shape, val, dtype=float) 67 if val.size == self.size: 68 val = np.stack([val, val]) # both components 69 val = val.reshape((2, *self.shape)) 70 # Set 71 super().__setattr__(key, val) 72 73 name: str = "Unnamed" 74 """Description.""" 75 76 vw: float = 1. 77 """Viscosity for water.""" 78 vo: float = 1. 79 """Viscosity for oil.""" 80 swc: float = 0. 81 """Irreducible saturation, water.""" 82 sor: float = 0. 83 """Irreducible saturation, oil.""" 84 K: np.ndarray = None 85 """Permeabilities (in x and y directions). Array of shape `(2, Nx, Ny)`).""" 86 por: np.ndarray = None 87 """Porosity; Array of shape `(Nx, Ny)`).""" 88 89 nInj = property(lambda self: len(self.inj_xy)) 90 """Num. of injector wells.""" 91 nPrd = property(lambda self: len(self.prd_xy)) 92 """Num. of producer wells.""" 93 94 inj_xy: np.ndarray = None 95 """Array of shape `(nWell, 2)` of x- and y-coords for `nWell` injectors. 96 97 Values should be betwen `0` and `Lx` or `Ly`. 98 99 .. warning:: The wells get co-located with grid nodes, ref `xy2sub`. 100 This is a design choice, not a mathematical necessity. 101 An alternative would be to distribute them over nearby nodes. 102 """ 103 prd_xy: np.ndarray = None 104 """Like `inj_xy`, but for producing wells.""" 105 inj_rates: np.ndarray = None 106 """Array of shape `(nWell, nTime)` -- or `(nWell, 1)` if constant-in-time. 107 108 .. note:: Both `inj_rates` and `prd_rates` are rates should be positive. 109 At each time index, it is asserted that the difference of their sums is 0, 110 otherwise the model would silently input deficit from SW corner. 111 """ 112 prd_rates: np.ndarray = None 113 """Like `prd_rates`, but for producing wells.""" 114 115 def _set_Q(self, S, k): 116 """Populate (for time `k`) the source/sink *field*, `Q`, from well specs.""" 117 Q = np.zeros(self.Nxy) 118 rates = self.dynamic_rate(S, k) 119 for kind in ['inj', 'prd']: 120 # Populate Q 121 xys = getattr(self, f'{kind}_xy') 122 sgn = +1 if kind == "inj" else -1 123 for xy, q in zip(xys, rates[kind]): 124 Q[self.xy2ind(*xy)] += sgn * q # += enables superimposition 125 # Store the computed/dynamic rates 126 if hasattr(self, "actual_rates"): 127 self.actual_rates[kind][:, k] = rates[kind] 128 self._Q = Q 129 130 def _wanted_rates_at(self, k): 131 """Lookup nominal/specified rates. Allows constant-in-time (singleton) spec.""" 132 get_now = lambda arr: np.copy(arr[k] if (len(arr) > 1) else arr[0]) 133 return (get_now(self.inj_rates.T), 134 get_now(self.prd_rates.T)) 135 136 def dynamic_rate(self, S, k): 137 """Compute the `actual_rates` for time index `k`. 138 139 This default implementation simply reads the given well specifications. 140 But you can overwrite (patch/inherit) it, for example to halt production wells 141 if water saturation is too high or simply if the suggested rate is near 0. 142 """ 143 inj, prd = self._wanted_rates_at(k) 144 return dict(inj=inj, prd=prd) 145 146 # Pres() -- listing 5 147 def pressure_step(self, S): 148 """Compute permeabilities then solve Darcy's equation. Returns `[P, V]`.""" 149 # Compute K*λ(S) 150 Mw, Mo = self.RelPerm(S) 151 Mt = Mw + Mo 152 Mt = Mt.reshape(self.shape) 153 KM = Mt * self.K 154 # Compute pressure and extract fluxes 155 [P, V] = self.TPFA(KM) 156 return P, V 157 158 def _spdiags(self, data, diags): 159 return sparse.spdiags(data, diags, self.Nxy, self.Nxy) 160 161 def rescale_sat(self, s): 162 """Account for irreducible saturations. Ref paper, p. 32.""" 163 return (s - self.swc) / (1 - self.swc - self.sor) 164 165 # RelPerm() -- listing 6 166 def RelPerm(self, s): 167 """Rel. permeabilities of oil and water. Return as mobilities (perm/viscocity).""" 168 S = self.rescale_sat(s) 169 Mw = S**2 / self.vw # Water mobility 170 Mo = (1 - S)**2 / self.vo # Oil mobility 171 return Mw, Mo 172 173 def dRelPerm(self, s): 174 """Derivatives of `RelPerm`.""" 175 S = self.rescale_sat(s) 176 dMw = 2 * S / self.vw / (1 - self.swc - self.sor) 177 dMo = -2 * (1 - S) / self.vo / (1 - self.swc - self.sor) 178 return dMw, dMo 179 180 # TPFA() -- Listing 1 181 def TPFA(self, K): 182 """Two-point flux-approximation (TPFA) of Darcy: $ -∇(K ∇u) = q $ 183 184 i.e. steady-state diffusion w/ nonlinear coefficient, $K$. 185 186 After solving for pressure `P`, extract the fluxes `V` 187 by finite differences. 188 """ 189 # Compute transmissibilities by harmonic averaging. 190 L = 1/K 191 TX = np.zeros((self.Nx + 1, self.Ny)) 192 TY = np.zeros((self.Nx, self.Ny + 1)) 193 TX[1:-1, :] = 2 * self.hy / self.hx / (L[0, :-1, :] + L[0, 1:, :]) 194 TY[:, 1:-1] = 2 * self.hx / self.hy / (L[1, :, :-1] + L[1, :, 1:]) 195 196 # Assemble TPFA discretization matrix. 197 x1 = TX[:-1, :].ravel() 198 x2 = TX[1:, :] .ravel() 199 y1 = TY[:, :-1].ravel() 200 y2 = TY[:, 1:] .ravel() 201 202 # Setup linear system 203 DiagVecs = [-x2, -y2, y1 + y2 + x1 + x2, -y1, -x1] 204 DiagIndx = [-self.Ny, -1, 0, 1, self.Ny] 205 DiagVecs[2][0] += np.sum(self.K[:, 0, 0]) # ref article p. 13 206 A = self._spdiags(DiagVecs, DiagIndx) 207 208 # Solve; compute A\q 209 q = self._Q 210 # u = np.linalg.solve(A.A, q) # direct dense solver 211 u = spsolve(A.tocsr(), q) # direct sparse solver 212 # u, _info = cg(A, q) # conjugate gradient 213 # Could also try scipy.linalg.solveh_banded which, according to 214 # https://scicomp.stackexchange.com/a/30074 uses the Thomas algorithm, 215 # as recommended by Aziz and Settari ("Petro. Res. simulation"). 216 # NB: stackexchange also mentions that solve_banded does not work well 217 # when the band offsets large, i.e. higher-dimensional problems. 218 219 # Extract fluxes 220 P = u.reshape(self.shape) 221 V = DotDict( 222 x=np.zeros((self.Nx+1, self.Ny)), 223 y=np.zeros((self.Nx, self.Ny+1)), 224 ) 225 V.x[1:-1, :] = (P[:-1, :] - P[1:, :]) * TX[1:-1, :] 226 V.y[:, 1:-1] = (P[:, :-1] - P[:, 1:]) * TY[:, 1:-1] 227 return P, V 228 229 # GenA() -- listing 7 230 def upwind_diff(self, V): 231 """Upwind finite-volume scheme.""" 232 fp = self._Q.clip(max=0) # production 233 # Flow fluxes, separated into direction (x-y) and sign 234 x1 = V.x.clip(max=0)[:-1, :].ravel() 235 y1 = V.y.clip(max=0)[:, :-1].ravel() 236 x2 = V.x.clip(min=0)[1:, :] .ravel() 237 y2 = V.y.clip(min=0)[:, 1:] .ravel() 238 DiagVecs = [x2, y2, fp + y1 - y2 + x1 - x2, -y1, -x1] 239 DiagIndx = [-self.Ny, -1, 0, 1, self.Ny] 240 A = self._spdiags(DiagVecs, DiagIndx) 241 return A 242 243 # Extracted from Upstream() 244 def estimate_1CFL(self, pv, V, fi): 245 """Estimate 1/CFL for use with `saturation_step_upwind`.""" 246 # In-/Out-flux x-/y- faces 247 XP = V.x.clip(min=0) 248 XN = V.x.clip(max=0) 249 YP = V.y.clip(min=0) 250 YN = V.y.clip(max=0) 251 Vi = XP[:-1, :] + YP[:, :-1] - XN[1:, :] - YN[:, 1:] 252 253 flx = max((Vi.ravel() + fi) / pv) # estimate of influx 254 sat = self.swc + self.sor 255 return 3 / (1 - sat) * flx # NB: 3-->2 since no z-dim ? 256 257 # Upstream() -- listing 8 258 def saturation_step_upwind(self, S, V, dt): 259 """Explicit upwind FV discretisation of conserv. of mass (water sat.).""" 260 A = self.upwind_diff(V) # FV discretized transport operator 261 pv = self.h2 * self.por.ravel() # Pore volume = cell volume * porosity 262 fi = self._Q.clip(min=0) # Well inflow 263 264 # Compute sub/local dt 265 cfl1 = self.estimate_1CFL(pv, V, fi) 266 nT = int(np.ceil(dt * cfl1)) 267 nT = max(1, nT) 268 269 # Scale A 270 dtx = dt / nT / pv # timestep / pore volume 271 B = self._spdiags(dtx, 0) @ A # A * dt/|Omega i| 272 273 for _ in range(nT): 274 Mw, Mo = self.RelPerm(S) # compute mobilities 275 fw = Mw / (Mw + Mo) # compute fractional flow 276 S = S + (B@fw + fi*dtx) # update saturation 277 return S 278 279 # NewtRaph() -- listing 10 280 def saturation_step_implicit(self, S, V, dt, nNewtonMax=10, nTmax_log2=10): 281 """Implicit FV discretisation of conserv. of mass (water sat.).""" 282 A = self.upwind_diff(V) # FV discretized transport operator 283 pv = self.h2 * self.por.ravel() # Pore volume = cell.vol * por 284 fi = self._Q.clip(min=0) # Well inflow 285 286 # For each iter, halve the sub/local dt 287 for nT_log2 in range(0, nTmax_log2): 288 nT = 2**nT_log2 289 290 # Scale A 291 dtx = dt / nT / pv # timestep / pore volume 292 B = self._spdiags(dtx, 0) @ A # A * dt/|Omega i| 293 294 Sn = S 295 for _ in range(nT): 296 Sp = Sn 297 for _ in range(nNewtonMax): 298 Mw, Mo = self.RelPerm(Sn) # mobilities 299 dMw, dMo = self.dRelPerm(Sn) # their derivatives 300 df = dMw/(Mw+Mo) - Mw/(Mw+Mo)**2 * (dMw + dMo) # df w/ds 301 dG = sparse.eye(self.Nxy) - B @ self._spdiags(df, 0) # deriv of G 302 303 fw = Mw / (Mw+Mo) # fract. flow 304 G = Sn - Sp - (B@fw + fi*dtx) # G(s) 305 dS = spsolve(dG, G) # compute dS 306 Sn = Sn - dS # update S 307 308 if np.sqrt(sum(dS**2)) < 1e-3: 309 # If converged: halt Newton iterations 310 break 311 else: 312 # If never converged: increase nT, restart time loop 313 break 314 else: 315 # If completed all time steps, halt 316 break 317 else: 318 # Failed (even with max nT) to complete all time steps 319 print("Warning: did not converge") 320 321 return Sn 322 323 def time_stepper(self, dt, implicit=False): 324 """Get ODE solver (integrator) for model. 325 326 Whatever time step `dt` is given, both schemes will use smaller steps internally. 327 328 - `explicit`: computes sub-`dt` based on CFL esitmate. 329 - `implicit`: reduces sub-`dt` until convergence is achieved. 330 """ 331 def integrate(S, k): 332 self._set_Q(S, k) 333 334 # Catch some common issues before they become mysterious/insidious 335 # (e.g. mass imblance silently inserts deficit in SW corner). 336 assert len(self.inj_rates) == len(self.inj_xy) 337 assert len(self.prd_rates) == len(self.prd_xy) 338 assert np.isclose(self._Q.sum(), 0), "(inj - prd) does not sum to 0" 339 assert np.all(self.inj_rates >= 0) 340 assert np.all(self.prd_rates >= 0) 341 assert np.all((0 <= self.K ) & np.isfinite(self.K)) 342 assert np.all((0 <= self.por) & (self.por <= 1)) 343 344 [P, V] = self.pressure_step(S) 345 if implicit: 346 S = self.saturation_step_implicit(S, V, dt) 347 else: 348 S = self.saturation_step_upwind(S, V, dt) 349 return S 350 return integrate 351 352 def sim(self, dt, nSteps, x0, pbar=True, leave=True, **kwargs): 353 """Recursively (`nSteps` times) apply `time_stepper` with `dt`, from `x0`. 354 355 .. note:: `output[0] == x0`, hence `len(output) = nSteps + 1`. 356 357 .. note:: 358 "Recurse" is a describes a function calling itself. 359 Here we implement it simply by a for loop. 360 "Wecursive" is also an accurate description of causal (Markov) processes, 361 such as nature or its simulators, which build on themselves. 362 """ 363 step = self.time_stepper(dt, **kwargs) 364 365 # pbar 366 kk = np.arange(nSteps) 367 if pbar: 368 kk = tqdm(kk, "Simulation", leave=leave, mininterval=1e-2) 369 370 # Init 371 xx = np.zeros((nSteps+1,)+x0.shape) 372 xx[0] = x0 373 self.actual_rates = dict(inj=np.zeros((self.nInj, nSteps)), 374 prd=np.zeros((self.nPrd, nSteps))) 375 376 # Recurse 377 for k in kk: 378 xx[k+1] = step(xx[k], k) 379 380 return xx
15@dataclass 16class ResSim(NicePrint, Grid2D, Plot2D): 17 """Reservoir simulator class. 18 19 Implemented with OOP (instead of passing around dicts) to facilitate 20 bookkeeping of ensemble forecasting 21 (where parameter values of one instance should not influence another) 22 23 Example: 24 >>> model = ResSim(Lx=1, Ly=1, Nx=64, Ny=64) 25 >>> model.inj_xy=[[0, .32]] 26 >>> model.prd_xy=[[1, 1]] 27 >>> model.inj_rates=[[1]] 28 >>> model.prd_rates=[[1]] 29 >>> water_sat0 = np.zeros(model.Nxy) 30 >>> dt = .35 31 >>> nSteps = 2 32 >>> S = model.sim(dt, nSteps, water_sat0, pbar=False) 33 34 This produces the following values (used for automatic testing): 35 >>> S[-1, [100, 1300, 2900]] 36 array([0.9429345 , 0.91358172, 0.71554613]) 37 """ 38 # Dont use dataclass repr 39 __repr__ = NicePrint.__repr__ 40 __str__ = NicePrint.__str__ 41 42 def __post_init__(self): 43 defaults = dict(K=np.ones((2, *self.shape)), 44 por=np.ones(self.shape)) 45 for k, v in defaults.items(): 46 if getattr(self, k) is None: 47 setattr(self, k, v) 48 49 # Prefer __setattr__ approach (over @property get/set-ers) 50 # because @property requires the _private pattern, 51 # which is pretty ugly with dataclasses, 52 # and also because can unify treatment of inj/prd wells. 53 def __setattr__(self, key, val): 54 if val is not None: 55 # Well positions -- collocate at some node 56 if key in ["inj_xy", "prd_xy"]: 57 val = np.array(val, float).reshape((-1, 2)) 58 for i, (x, y) in enumerate(val): 59 val[i] = self.ind2xy(self.xy2ind(x, y)) 60 # Well rates 61 if key in ["inj_rates", "prd_rates"]: 62 nWell = len(getattr(self, key.replace("rates", "xy"))) 63 val = np.array(val, float).reshape((nWell, -1)) 64 # Permeabilities 65 if key == "K": 66 if np.isscalar(val): 67 val = np.full_like(self.shape, val, dtype=float) 68 if val.size == self.size: 69 val = np.stack([val, val]) # both components 70 val = val.reshape((2, *self.shape)) 71 # Set 72 super().__setattr__(key, val) 73 74 name: str = "Unnamed" 75 """Description.""" 76 77 vw: float = 1. 78 """Viscosity for water.""" 79 vo: float = 1. 80 """Viscosity for oil.""" 81 swc: float = 0. 82 """Irreducible saturation, water.""" 83 sor: float = 0. 84 """Irreducible saturation, oil.""" 85 K: np.ndarray = None 86 """Permeabilities (in x and y directions). Array of shape `(2, Nx, Ny)`).""" 87 por: np.ndarray = None 88 """Porosity; Array of shape `(Nx, Ny)`).""" 89 90 nInj = property(lambda self: len(self.inj_xy)) 91 """Num. of injector wells.""" 92 nPrd = property(lambda self: len(self.prd_xy)) 93 """Num. of producer wells.""" 94 95 inj_xy: np.ndarray = None 96 """Array of shape `(nWell, 2)` of x- and y-coords for `nWell` injectors. 97 98 Values should be betwen `0` and `Lx` or `Ly`. 99 100 .. warning:: The wells get co-located with grid nodes, ref `xy2sub`. 101 This is a design choice, not a mathematical necessity. 102 An alternative would be to distribute them over nearby nodes. 103 """ 104 prd_xy: np.ndarray = None 105 """Like `inj_xy`, but for producing wells.""" 106 inj_rates: np.ndarray = None 107 """Array of shape `(nWell, nTime)` -- or `(nWell, 1)` if constant-in-time. 108 109 .. note:: Both `inj_rates` and `prd_rates` are rates should be positive. 110 At each time index, it is asserted that the difference of their sums is 0, 111 otherwise the model would silently input deficit from SW corner. 112 """ 113 prd_rates: np.ndarray = None 114 """Like `prd_rates`, but for producing wells.""" 115 116 def _set_Q(self, S, k): 117 """Populate (for time `k`) the source/sink *field*, `Q`, from well specs.""" 118 Q = np.zeros(self.Nxy) 119 rates = self.dynamic_rate(S, k) 120 for kind in ['inj', 'prd']: 121 # Populate Q 122 xys = getattr(self, f'{kind}_xy') 123 sgn = +1 if kind == "inj" else -1 124 for xy, q in zip(xys, rates[kind]): 125 Q[self.xy2ind(*xy)] += sgn * q # += enables superimposition 126 # Store the computed/dynamic rates 127 if hasattr(self, "actual_rates"): 128 self.actual_rates[kind][:, k] = rates[kind] 129 self._Q = Q 130 131 def _wanted_rates_at(self, k): 132 """Lookup nominal/specified rates. Allows constant-in-time (singleton) spec.""" 133 get_now = lambda arr: np.copy(arr[k] if (len(arr) > 1) else arr[0]) 134 return (get_now(self.inj_rates.T), 135 get_now(self.prd_rates.T)) 136 137 def dynamic_rate(self, S, k): 138 """Compute the `actual_rates` for time index `k`. 139 140 This default implementation simply reads the given well specifications. 141 But you can overwrite (patch/inherit) it, for example to halt production wells 142 if water saturation is too high or simply if the suggested rate is near 0. 143 """ 144 inj, prd = self._wanted_rates_at(k) 145 return dict(inj=inj, prd=prd) 146 147 # Pres() -- listing 5 148 def pressure_step(self, S): 149 """Compute permeabilities then solve Darcy's equation. Returns `[P, V]`.""" 150 # Compute K*λ(S) 151 Mw, Mo = self.RelPerm(S) 152 Mt = Mw + Mo 153 Mt = Mt.reshape(self.shape) 154 KM = Mt * self.K 155 # Compute pressure and extract fluxes 156 [P, V] = self.TPFA(KM) 157 return P, V 158 159 def _spdiags(self, data, diags): 160 return sparse.spdiags(data, diags, self.Nxy, self.Nxy) 161 162 def rescale_sat(self, s): 163 """Account for irreducible saturations. Ref paper, p. 32.""" 164 return (s - self.swc) / (1 - self.swc - self.sor) 165 166 # RelPerm() -- listing 6 167 def RelPerm(self, s): 168 """Rel. permeabilities of oil and water. Return as mobilities (perm/viscocity).""" 169 S = self.rescale_sat(s) 170 Mw = S**2 / self.vw # Water mobility 171 Mo = (1 - S)**2 / self.vo # Oil mobility 172 return Mw, Mo 173 174 def dRelPerm(self, s): 175 """Derivatives of `RelPerm`.""" 176 S = self.rescale_sat(s) 177 dMw = 2 * S / self.vw / (1 - self.swc - self.sor) 178 dMo = -2 * (1 - S) / self.vo / (1 - self.swc - self.sor) 179 return dMw, dMo 180 181 # TPFA() -- Listing 1 182 def TPFA(self, K): 183 """Two-point flux-approximation (TPFA) of Darcy: $ -∇(K ∇u) = q $ 184 185 i.e. steady-state diffusion w/ nonlinear coefficient, $K$. 186 187 After solving for pressure `P`, extract the fluxes `V` 188 by finite differences. 189 """ 190 # Compute transmissibilities by harmonic averaging. 191 L = 1/K 192 TX = np.zeros((self.Nx + 1, self.Ny)) 193 TY = np.zeros((self.Nx, self.Ny + 1)) 194 TX[1:-1, :] = 2 * self.hy / self.hx / (L[0, :-1, :] + L[0, 1:, :]) 195 TY[:, 1:-1] = 2 * self.hx / self.hy / (L[1, :, :-1] + L[1, :, 1:]) 196 197 # Assemble TPFA discretization matrix. 198 x1 = TX[:-1, :].ravel() 199 x2 = TX[1:, :] .ravel() 200 y1 = TY[:, :-1].ravel() 201 y2 = TY[:, 1:] .ravel() 202 203 # Setup linear system 204 DiagVecs = [-x2, -y2, y1 + y2 + x1 + x2, -y1, -x1] 205 DiagIndx = [-self.Ny, -1, 0, 1, self.Ny] 206 DiagVecs[2][0] += np.sum(self.K[:, 0, 0]) # ref article p. 13 207 A = self._spdiags(DiagVecs, DiagIndx) 208 209 # Solve; compute A\q 210 q = self._Q 211 # u = np.linalg.solve(A.A, q) # direct dense solver 212 u = spsolve(A.tocsr(), q) # direct sparse solver 213 # u, _info = cg(A, q) # conjugate gradient 214 # Could also try scipy.linalg.solveh_banded which, according to 215 # https://scicomp.stackexchange.com/a/30074 uses the Thomas algorithm, 216 # as recommended by Aziz and Settari ("Petro. Res. simulation"). 217 # NB: stackexchange also mentions that solve_banded does not work well 218 # when the band offsets large, i.e. higher-dimensional problems. 219 220 # Extract fluxes 221 P = u.reshape(self.shape) 222 V = DotDict( 223 x=np.zeros((self.Nx+1, self.Ny)), 224 y=np.zeros((self.Nx, self.Ny+1)), 225 ) 226 V.x[1:-1, :] = (P[:-1, :] - P[1:, :]) * TX[1:-1, :] 227 V.y[:, 1:-1] = (P[:, :-1] - P[:, 1:]) * TY[:, 1:-1] 228 return P, V 229 230 # GenA() -- listing 7 231 def upwind_diff(self, V): 232 """Upwind finite-volume scheme.""" 233 fp = self._Q.clip(max=0) # production 234 # Flow fluxes, separated into direction (x-y) and sign 235 x1 = V.x.clip(max=0)[:-1, :].ravel() 236 y1 = V.y.clip(max=0)[:, :-1].ravel() 237 x2 = V.x.clip(min=0)[1:, :] .ravel() 238 y2 = V.y.clip(min=0)[:, 1:] .ravel() 239 DiagVecs = [x2, y2, fp + y1 - y2 + x1 - x2, -y1, -x1] 240 DiagIndx = [-self.Ny, -1, 0, 1, self.Ny] 241 A = self._spdiags(DiagVecs, DiagIndx) 242 return A 243 244 # Extracted from Upstream() 245 def estimate_1CFL(self, pv, V, fi): 246 """Estimate 1/CFL for use with `saturation_step_upwind`.""" 247 # In-/Out-flux x-/y- faces 248 XP = V.x.clip(min=0) 249 XN = V.x.clip(max=0) 250 YP = V.y.clip(min=0) 251 YN = V.y.clip(max=0) 252 Vi = XP[:-1, :] + YP[:, :-1] - XN[1:, :] - YN[:, 1:] 253 254 flx = max((Vi.ravel() + fi) / pv) # estimate of influx 255 sat = self.swc + self.sor 256 return 3 / (1 - sat) * flx # NB: 3-->2 since no z-dim ? 257 258 # Upstream() -- listing 8 259 def saturation_step_upwind(self, S, V, dt): 260 """Explicit upwind FV discretisation of conserv. of mass (water sat.).""" 261 A = self.upwind_diff(V) # FV discretized transport operator 262 pv = self.h2 * self.por.ravel() # Pore volume = cell volume * porosity 263 fi = self._Q.clip(min=0) # Well inflow 264 265 # Compute sub/local dt 266 cfl1 = self.estimate_1CFL(pv, V, fi) 267 nT = int(np.ceil(dt * cfl1)) 268 nT = max(1, nT) 269 270 # Scale A 271 dtx = dt / nT / pv # timestep / pore volume 272 B = self._spdiags(dtx, 0) @ A # A * dt/|Omega i| 273 274 for _ in range(nT): 275 Mw, Mo = self.RelPerm(S) # compute mobilities 276 fw = Mw / (Mw + Mo) # compute fractional flow 277 S = S + (B@fw + fi*dtx) # update saturation 278 return S 279 280 # NewtRaph() -- listing 10 281 def saturation_step_implicit(self, S, V, dt, nNewtonMax=10, nTmax_log2=10): 282 """Implicit FV discretisation of conserv. of mass (water sat.).""" 283 A = self.upwind_diff(V) # FV discretized transport operator 284 pv = self.h2 * self.por.ravel() # Pore volume = cell.vol * por 285 fi = self._Q.clip(min=0) # Well inflow 286 287 # For each iter, halve the sub/local dt 288 for nT_log2 in range(0, nTmax_log2): 289 nT = 2**nT_log2 290 291 # Scale A 292 dtx = dt / nT / pv # timestep / pore volume 293 B = self._spdiags(dtx, 0) @ A # A * dt/|Omega i| 294 295 Sn = S 296 for _ in range(nT): 297 Sp = Sn 298 for _ in range(nNewtonMax): 299 Mw, Mo = self.RelPerm(Sn) # mobilities 300 dMw, dMo = self.dRelPerm(Sn) # their derivatives 301 df = dMw/(Mw+Mo) - Mw/(Mw+Mo)**2 * (dMw + dMo) # df w/ds 302 dG = sparse.eye(self.Nxy) - B @ self._spdiags(df, 0) # deriv of G 303 304 fw = Mw / (Mw+Mo) # fract. flow 305 G = Sn - Sp - (B@fw + fi*dtx) # G(s) 306 dS = spsolve(dG, G) # compute dS 307 Sn = Sn - dS # update S 308 309 if np.sqrt(sum(dS**2)) < 1e-3: 310 # If converged: halt Newton iterations 311 break 312 else: 313 # If never converged: increase nT, restart time loop 314 break 315 else: 316 # If completed all time steps, halt 317 break 318 else: 319 # Failed (even with max nT) to complete all time steps 320 print("Warning: did not converge") 321 322 return Sn 323 324 def time_stepper(self, dt, implicit=False): 325 """Get ODE solver (integrator) for model. 326 327 Whatever time step `dt` is given, both schemes will use smaller steps internally. 328 329 - `explicit`: computes sub-`dt` based on CFL esitmate. 330 - `implicit`: reduces sub-`dt` until convergence is achieved. 331 """ 332 def integrate(S, k): 333 self._set_Q(S, k) 334 335 # Catch some common issues before they become mysterious/insidious 336 # (e.g. mass imblance silently inserts deficit in SW corner). 337 assert len(self.inj_rates) == len(self.inj_xy) 338 assert len(self.prd_rates) == len(self.prd_xy) 339 assert np.isclose(self._Q.sum(), 0), "(inj - prd) does not sum to 0" 340 assert np.all(self.inj_rates >= 0) 341 assert np.all(self.prd_rates >= 0) 342 assert np.all((0 <= self.K ) & np.isfinite(self.K)) 343 assert np.all((0 <= self.por) & (self.por <= 1)) 344 345 [P, V] = self.pressure_step(S) 346 if implicit: 347 S = self.saturation_step_implicit(S, V, dt) 348 else: 349 S = self.saturation_step_upwind(S, V, dt) 350 return S 351 return integrate 352 353 def sim(self, dt, nSteps, x0, pbar=True, leave=True, **kwargs): 354 """Recursively (`nSteps` times) apply `time_stepper` with `dt`, from `x0`. 355 356 .. note:: `output[0] == x0`, hence `len(output) = nSteps + 1`. 357 358 .. note:: 359 "Recurse" is a describes a function calling itself. 360 Here we implement it simply by a for loop. 361 "Wecursive" is also an accurate description of causal (Markov) processes, 362 such as nature or its simulators, which build on themselves. 363 """ 364 step = self.time_stepper(dt, **kwargs) 365 366 # pbar 367 kk = np.arange(nSteps) 368 if pbar: 369 kk = tqdm(kk, "Simulation", leave=leave, mininterval=1e-2) 370 371 # Init 372 xx = np.zeros((nSteps+1,)+x0.shape) 373 xx[0] = x0 374 self.actual_rates = dict(inj=np.zeros((self.nInj, nSteps)), 375 prd=np.zeros((self.nPrd, nSteps))) 376 377 # Recurse 378 for k in kk: 379 xx[k+1] = step(xx[k], k) 380 381 return xx
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)
>>> model.inj_xy=[[0, .32]]
>>> model.prd_xy=[[1, 1]]
>>> model.inj_rates=[[1]]
>>> model.prd_rates=[[1]]
>>> water_sat0 = np.zeros(model.Nxy)
>>> dt = .35
>>> nSteps = 2
>>> S = 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])
137 def dynamic_rate(self, S, k): 138 """Compute the `actual_rates` for time index `k`. 139 140 This default implementation simply reads the given well specifications. 141 But you can overwrite (patch/inherit) it, for example to halt production wells 142 if water saturation is too high or simply if the suggested rate is near 0. 143 """ 144 inj, prd = self._wanted_rates_at(k) 145 return dict(inj=inj, prd=prd)
Compute the actual_rates for time index k.
This default implementation simply reads the given well specifications. But you can overwrite (patch/inherit) it, for example to halt production wells if water saturation is too high or simply if the suggested rate is near 0.
148 def pressure_step(self, S): 149 """Compute permeabilities then solve Darcy's equation. Returns `[P, V]`.""" 150 # Compute K*λ(S) 151 Mw, Mo = self.RelPerm(S) 152 Mt = Mw + Mo 153 Mt = Mt.reshape(self.shape) 154 KM = Mt * self.K 155 # Compute pressure and extract fluxes 156 [P, V] = self.TPFA(KM) 157 return P, V
Compute permeabilities then solve Darcy's equation. Returns [P, V].
162 def rescale_sat(self, s): 163 """Account for irreducible saturations. Ref paper, p. 32.""" 164 return (s - self.swc) / (1 - self.swc - self.sor)
Account for irreducible saturations. Ref paper, p. 32.
167 def RelPerm(self, s): 168 """Rel. permeabilities of oil and water. Return as mobilities (perm/viscocity).""" 169 S = self.rescale_sat(s) 170 Mw = S**2 / self.vw # Water mobility 171 Mo = (1 - S)**2 / self.vo # Oil mobility 172 return Mw, Mo
Rel. permeabilities of oil and water. Return as mobilities (perm/viscocity).
174 def dRelPerm(self, s): 175 """Derivatives of `RelPerm`.""" 176 S = self.rescale_sat(s) 177 dMw = 2 * S / self.vw / (1 - self.swc - self.sor) 178 dMo = -2 * (1 - S) / self.vo / (1 - self.swc - self.sor) 179 return dMw, dMo
Derivatives of RelPerm.
182 def TPFA(self, K): 183 """Two-point flux-approximation (TPFA) of Darcy: $ -∇(K ∇u) = q $ 184 185 i.e. steady-state diffusion w/ nonlinear coefficient, $K$. 186 187 After solving for pressure `P`, extract the fluxes `V` 188 by finite differences. 189 """ 190 # Compute transmissibilities by harmonic averaging. 191 L = 1/K 192 TX = np.zeros((self.Nx + 1, self.Ny)) 193 TY = np.zeros((self.Nx, self.Ny + 1)) 194 TX[1:-1, :] = 2 * self.hy / self.hx / (L[0, :-1, :] + L[0, 1:, :]) 195 TY[:, 1:-1] = 2 * self.hx / self.hy / (L[1, :, :-1] + L[1, :, 1:]) 196 197 # Assemble TPFA discretization matrix. 198 x1 = TX[:-1, :].ravel() 199 x2 = TX[1:, :] .ravel() 200 y1 = TY[:, :-1].ravel() 201 y2 = TY[:, 1:] .ravel() 202 203 # Setup linear system 204 DiagVecs = [-x2, -y2, y1 + y2 + x1 + x2, -y1, -x1] 205 DiagIndx = [-self.Ny, -1, 0, 1, self.Ny] 206 DiagVecs[2][0] += np.sum(self.K[:, 0, 0]) # ref article p. 13 207 A = self._spdiags(DiagVecs, DiagIndx) 208 209 # Solve; compute A\q 210 q = self._Q 211 # u = np.linalg.solve(A.A, q) # direct dense solver 212 u = spsolve(A.tocsr(), q) # direct sparse solver 213 # u, _info = cg(A, q) # conjugate gradient 214 # Could also try scipy.linalg.solveh_banded which, according to 215 # https://scicomp.stackexchange.com/a/30074 uses the Thomas algorithm, 216 # as recommended by Aziz and Settari ("Petro. Res. simulation"). 217 # NB: stackexchange also mentions that solve_banded does not work well 218 # when the band offsets large, i.e. higher-dimensional problems. 219 220 # Extract fluxes 221 P = u.reshape(self.shape) 222 V = DotDict( 223 x=np.zeros((self.Nx+1, self.Ny)), 224 y=np.zeros((self.Nx, self.Ny+1)), 225 ) 226 V.x[1:-1, :] = (P[:-1, :] - P[1:, :]) * TX[1:-1, :] 227 V.y[:, 1:-1] = (P[:, :-1] - P[:, 1:]) * TY[:, 1:-1] 228 return P, V
Two-point flux-approximation (TPFA) of Darcy: $ -∇(K ∇u) = q $
i.e. steady-state diffusion w/ nonlinear coefficient, $K$.
After solving for pressure P, extract the fluxes V
by finite differences.
231 def upwind_diff(self, V): 232 """Upwind finite-volume scheme.""" 233 fp = self._Q.clip(max=0) # production 234 # Flow fluxes, separated into direction (x-y) and sign 235 x1 = V.x.clip(max=0)[:-1, :].ravel() 236 y1 = V.y.clip(max=0)[:, :-1].ravel() 237 x2 = V.x.clip(min=0)[1:, :] .ravel() 238 y2 = V.y.clip(min=0)[:, 1:] .ravel() 239 DiagVecs = [x2, y2, fp + y1 - y2 + x1 - x2, -y1, -x1] 240 DiagIndx = [-self.Ny, -1, 0, 1, self.Ny] 241 A = self._spdiags(DiagVecs, DiagIndx) 242 return A
Upwind finite-volume scheme.
245 def estimate_1CFL(self, pv, V, fi): 246 """Estimate 1/CFL for use with `saturation_step_upwind`.""" 247 # In-/Out-flux x-/y- faces 248 XP = V.x.clip(min=0) 249 XN = V.x.clip(max=0) 250 YP = V.y.clip(min=0) 251 YN = V.y.clip(max=0) 252 Vi = XP[:-1, :] + YP[:, :-1] - XN[1:, :] - YN[:, 1:] 253 254 flx = max((Vi.ravel() + fi) / pv) # estimate of influx 255 sat = self.swc + self.sor 256 return 3 / (1 - sat) * flx # NB: 3-->2 since no z-dim ?
Estimate 1/CFL for use with saturation_step_upwind.
259 def saturation_step_upwind(self, S, V, dt): 260 """Explicit upwind FV discretisation of conserv. of mass (water sat.).""" 261 A = self.upwind_diff(V) # FV discretized transport operator 262 pv = self.h2 * self.por.ravel() # Pore volume = cell volume * porosity 263 fi = self._Q.clip(min=0) # Well inflow 264 265 # Compute sub/local dt 266 cfl1 = self.estimate_1CFL(pv, V, fi) 267 nT = int(np.ceil(dt * cfl1)) 268 nT = max(1, nT) 269 270 # Scale A 271 dtx = dt / nT / pv # timestep / pore volume 272 B = self._spdiags(dtx, 0) @ A # A * dt/|Omega i| 273 274 for _ in range(nT): 275 Mw, Mo = self.RelPerm(S) # compute mobilities 276 fw = Mw / (Mw + Mo) # compute fractional flow 277 S = S + (B@fw + fi*dtx) # update saturation 278 return S
Explicit upwind FV discretisation of conserv. of mass (water sat.).
281 def saturation_step_implicit(self, S, V, dt, nNewtonMax=10, nTmax_log2=10): 282 """Implicit FV discretisation of conserv. of mass (water sat.).""" 283 A = self.upwind_diff(V) # FV discretized transport operator 284 pv = self.h2 * self.por.ravel() # Pore volume = cell.vol * por 285 fi = self._Q.clip(min=0) # Well inflow 286 287 # For each iter, halve the sub/local dt 288 for nT_log2 in range(0, nTmax_log2): 289 nT = 2**nT_log2 290 291 # Scale A 292 dtx = dt / nT / pv # timestep / pore volume 293 B = self._spdiags(dtx, 0) @ A # A * dt/|Omega i| 294 295 Sn = S 296 for _ in range(nT): 297 Sp = Sn 298 for _ in range(nNewtonMax): 299 Mw, Mo = self.RelPerm(Sn) # mobilities 300 dMw, dMo = self.dRelPerm(Sn) # their derivatives 301 df = dMw/(Mw+Mo) - Mw/(Mw+Mo)**2 * (dMw + dMo) # df w/ds 302 dG = sparse.eye(self.Nxy) - B @ self._spdiags(df, 0) # deriv of G 303 304 fw = Mw / (Mw+Mo) # fract. flow 305 G = Sn - Sp - (B@fw + fi*dtx) # G(s) 306 dS = spsolve(dG, G) # compute dS 307 Sn = Sn - dS # update S 308 309 if np.sqrt(sum(dS**2)) < 1e-3: 310 # If converged: halt Newton iterations 311 break 312 else: 313 # If never converged: increase nT, restart time loop 314 break 315 else: 316 # If completed all time steps, halt 317 break 318 else: 319 # Failed (even with max nT) to complete all time steps 320 print("Warning: did not converge") 321 322 return Sn
Implicit FV discretisation of conserv. of mass (water sat.).
324 def time_stepper(self, dt, implicit=False): 325 """Get ODE solver (integrator) for model. 326 327 Whatever time step `dt` is given, both schemes will use smaller steps internally. 328 329 - `explicit`: computes sub-`dt` based on CFL esitmate. 330 - `implicit`: reduces sub-`dt` until convergence is achieved. 331 """ 332 def integrate(S, k): 333 self._set_Q(S, k) 334 335 # Catch some common issues before they become mysterious/insidious 336 # (e.g. mass imblance silently inserts deficit in SW corner). 337 assert len(self.inj_rates) == len(self.inj_xy) 338 assert len(self.prd_rates) == len(self.prd_xy) 339 assert np.isclose(self._Q.sum(), 0), "(inj - prd) does not sum to 0" 340 assert np.all(self.inj_rates >= 0) 341 assert np.all(self.prd_rates >= 0) 342 assert np.all((0 <= self.K ) & np.isfinite(self.K)) 343 assert np.all((0 <= self.por) & (self.por <= 1)) 344 345 [P, V] = self.pressure_step(S) 346 if implicit: 347 S = self.saturation_step_implicit(S, V, dt) 348 else: 349 S = self.saturation_step_upwind(S, V, dt) 350 return S 351 return integrate
Get ODE solver (integrator) for model.
Whatever time step dt is given, both schemes will use smaller steps internally.
explicit: computes sub-dtbased on CFL esitmate.implicit: reduces sub-dtuntil convergence is achieved.
353 def sim(self, dt, nSteps, x0, pbar=True, leave=True, **kwargs): 354 """Recursively (`nSteps` times) apply `time_stepper` with `dt`, from `x0`. 355 356 .. note:: `output[0] == x0`, hence `len(output) = nSteps + 1`. 357 358 .. note:: 359 "Recurse" is a describes a function calling itself. 360 Here we implement it simply by a for loop. 361 "Wecursive" is also an accurate description of causal (Markov) processes, 362 such as nature or its simulators, which build on themselves. 363 """ 364 step = self.time_stepper(dt, **kwargs) 365 366 # pbar 367 kk = np.arange(nSteps) 368 if pbar: 369 kk = tqdm(kk, "Simulation", leave=leave, mininterval=1e-2) 370 371 # Init 372 xx = np.zeros((nSteps+1,)+x0.shape) 373 xx[0] = x0 374 self.actual_rates = dict(inj=np.zeros((self.nInj, nSteps)), 375 prd=np.zeros((self.nPrd, nSteps))) 376 377 # Recurse 378 for k in kk: 379 xx[k+1] = step(xx[k], k) 380 381 return xx
Recursively (nSteps times) apply time_stepper with dt, from x0.
output[0] == x0, hence len(output) = nSteps + 1.
"Recurse" is a describes a function calling itself. Here we implement it simply by a for loop. "Wecursive" is also an accurate description of causal (Markov) processes, such as nature or its simulators, which build on themselves.