Exercise 11.10
Capstone - A 2D Quarter-Five-Spot Pressure Solver
Extend the 1D simulator to 2D by adding a second spatial dimension. Use a 20×20 grid. Place an injector in one corner and a producer in the opposite corner. Visualize the pressure field as a 2D contour map at several time steps. This is a significant coding challenge, so take it step by step.
---
Time to put the whole chapter together. Every exercise so far has lived on a single line of cells; a real waterflood is areal. On OML 58 the team is screening a quarter-five-spot pilot: one water injector in a corner of a square sand, one producer in the diagonally opposite corner, and wants to see the steady-state pressure field before committing to a pattern. You will build the minimal 2D solver that produces it.
The physics is the same single-phase, slightly-compressible Darcy flow you already know, just on a 2D grid with a 5-point stencil: every interior cell connects to its four neighbours (left, right, up, down). We solve one implicit (steady-state) step (no time-marching) so the linear system
is assembled once with scipy.sparse and solved once with spsolve. The transmissibility and Peaceman well-index helpers are embedded; use them, do not re-derive them.
The grid is nx by ny square cells, indexed so that cell (i, j) (column i along x, row j along y) maps to the flat unknown idx(i, j) = j nx + i. The injector sits at cell (0, 0) held at a high BHP p_inj_psi; the producer sits at the opposite corner (nx-1, ny-1) held at a low BHP p_prod_psi. Both are BHP-controlled wells, so each adds a Peaceman well index wi to its own diagonal and wi p_bhp to the right-hand side.
Your task: write solve_2d(nx, ny, ...) (signature in the starter):
- Build the
N x Nsystem withN = nx * ny, using a sparse builder
(scipy.sparse.lil_matrix is easiest to assemble; convert to CSR/CSC before solving). Loop over every cell (i, j):
- For each in-grid neighbour, compute the transmissibility with the embedded
tx(a, b) (x-direction) or ty(a, b) (y-direction) helper, put -t in the off-diagonal A[p, neighbour], and add t to the cell's diagonal.
- With no source term and no accumulation, an interior cell's row sums to
zero. The diagonal is just the sum of its neighbour transmissibilities.
- Add the two BHP wells. For a BHP well at cell
c:A[c, c] += wi(c)and
b[c] += wi(c) * p_bhp, exactly like the 1D Peaceman well you have seen.
- Solve
A p = bwithspsolveandreturn p.reshape(ny, nx), a 2D
array indexed [j, i] (row = y, column = x).
> Think about it: with the injector pinned high and the producer pinned low, > pressure must be highest at (0,0) and lowest at (nx-1, ny-1). Because > the grid is square and homogeneous, the quarter-five-spot is symmetric across > the main diagonal: swapping x and y maps the injector onto itself and the > producer onto itself, so P[i, j] == P[j, i]. If your field is not symmetric, > your x- and y-transmissibilities (or your idx) are inconsistent. A correct > steady field has no NaNs and varies smoothly from corner to corner.
Stuck? Reveal hints one at a time — they progress from nudge to near-solution.
visibilityReveal reference solutionexpand_more
Try solving it yourself first — the hints walk you through it. The solution below is one valid approach; yours may differ and still be correct.
import numpy as np
from scipy.sparse import lil_matrix, csr_matrix
from scipy.sparse.linalg import spsolve
# ── OML-58 quarter-five-spot pilot constants (do not edit) ───────────────
NX = 20 # grid cells in x
NY = 20 # grid cells in y
LENGTH_FT = 2000.0 # reservoir extent in x
WIDTH_FT = 2000.0 # reservoir extent in y
HEIGHT_FT = 50.0 # net pay thickness
K_MD = 120.0 # homogeneous permeability, md
MU_CP = 1.0 # fluid viscosity, cp
BO = 1.0 # formation volume factor, rb/stb
P_INJ_PSI = 4000.0 # injector bottomhole pressure (corner 0,0)
P_PROD_PSI = 1000.0 # producer bottomhole pressure (corner nx-1,ny-1)
RW_FT = 0.354 # wellbore radius
SKIN = 0.0 # mechanical skin
BETA = 1.127e-3 # Darcy field-unit constant
def _make_helpers(nx, ny, length_ft, width_ft, height_ft, k_md, mu_cp, Bo,
rw_ft, skin):
"""Return (idx, tx, ty, wi, dx) for a homogeneous nx-by-ny grid.
idx(i, j) -> flat unknown index (column i along x, row j along y).
tx(a, b) / ty(a, b) -> harmonic-mean transmissibility between two flat
cell indices across an x- or y-face. wi(c) -> Peaceman well index at cell c.
These mirror the verified 1D simulator (T = BETA*k*A/(mu*B*L)).
"""
dx = length_ft / nx
dy = width_ft / ny
k = (np.full(nx * ny, k_md, float) if np.isscalar(k_md)
else np.asarray(k_md, float).ravel())
area_x = dy * height_ft # face area normal to x-flow
area_y = dx * height_ft # face area normal to y-flow
def idx(i, j):
return j * nx + i
def tx(a, b):
kh = 2.0 * k[a] * k[b] / (k[a] + k[b])
return BETA * kh * area_x / (mu_cp * Bo * dx)
def ty(a, b):
kh = 2.0 * k[a] * k[b] / (k[a] + k[b])
return BETA * kh * area_y / (mu_cp * Bo * dy)
def wi(c):
re = 0.2 * dx
return (2.0 * np.pi * BETA * k[c] * height_ft
/ (mu_cp * Bo * (np.log(re / rw_ft) + skin)))
return idx, tx, ty, wi, dx
def solve_2d(nx, ny, length_ft=LENGTH_FT, width_ft=WIDTH_FT,
height_ft=HEIGHT_FT, k_md=K_MD, mu_cp=MU_CP, Bo=BO,
p_inj_psi=P_INJ_PSI, p_prod_psi=P_PROD_PSI,
rw_ft=RW_FT, skin=SKIN):
"""Steady-state 2D pressure for a quarter-five-spot.
Injector (high BHP) at cell (0,0); producer (low BHP) at (nx-1,ny-1).
Assembles the N x N 5-point stencil (N = nx*ny) and solves A p = b once.
Returns the pressure field as a 2D array p[j, i] (row=y, col=x).
"""
idx, tx, ty, wi, _ = _make_helpers(nx, ny, length_ft, width_ft, height_ft,
k_md, mu_cp, Bo, rw_ft, skin)
N = nx * ny
A = lil_matrix((N, N))
b = np.zeros(N)
# 5-point stencil: each in-grid neighbour adds -t off-diagonal and +t on
# the diagonal. No source / no accumulation -> interior rows sum to zero.
for j in range(ny):
for i in range(nx):
p = idx(i, j)
diag = 0.0
if i > 0:
t = tx(p, idx(i - 1, j)); A[p, idx(i - 1, j)] = -t; diag += t
if i < nx - 1:
t = tx(p, idx(i + 1, j)); A[p, idx(i + 1, j)] = -t; diag += t
if j > 0:
t = ty(p, idx(i, j - 1)); A[p, idx(i, j - 1)] = -t; diag += t
if j < ny - 1:
t = ty(p, idx(i, j + 1)); A[p, idx(i, j + 1)] = -t; diag += t
A[p, p] = diag
# BHP wells (Peaceman): wi on the diagonal, wi*p_bhp on the RHS.
inj = idx(0, 0)
prod = idx(nx - 1, ny - 1)
wi_i = wi(inj); A[inj, inj] += wi_i; b[inj] += wi_i * p_inj_psi
wi_p = wi(prod); A[prod, prod] += wi_p; b[prod] += wi_p * p_prod_psi
p = spsolve(csr_matrix(A), b)
return p.reshape(ny, nx)
P = solve_2d(NX, NY)
print("field shape:", P.shape, "(should be (ny, nx))")
print("injector corner P[0,0] =", round(float(P[0, 0]), 1), "psi")
print("producer corner P[-1,-1]=", round(float(P[-1, -1]), 1), "psi")
print("max-min:", round(float(P.max()), 1), "->", round(float(P.min()), 1), "psi")
print("diagonal symmetry |P - P.T| max:", float(np.abs(P - P.T).max()))
lockCopying code is a Full Access feature.