Time-Step Sensitivity of the Implicit Simulator

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 diags
from scipy.sparse.linalg import spsolve


# ── Embedded 1D implicit single-phase simulator (do not edit) ────────────
class Reservoir1D:
    """1D single-phase implicit simulator (harmonic-mean T, real area, Peaceman wells)."""
    def __init__(self, nx, length_ft, width_ft, height_ft, k_md, phi, ct_psi, mu_cp, Pi_psi, Bo=1.0):
        self.nx = nx; self.dx = length_ft / nx
        self.area = width_ft * height_ft; self.vb = self.area * self.dx; self.height = height_ft
        self.k = np.full(nx, k_md, float) if np.isscalar(k_md) else np.asarray(k_md, float)
        self.phi = phi; self.ct = ct_psi; self.mu = mu_cp; self.Bo = Bo
        self.P = np.full(nx, Pi_psi, float); self.Pi = Pi_psi; self.beta = 1.127e-3; self.wells = []
    def _T(self, i, j):
        kh = 2*self.k[i]*self.k[j]/(self.k[i]+self.k[j])
        return self.beta*kh*self.area/(self.mu*self.dx)
    def _wi(self, w):
        re = 0.2*self.dx
        return 2*np.pi*self.beta*self.k[w["cell"]]*self.height/(self.mu*(np.log(re/w["rw"])+w["skin"]))
    def add_well(self, cell, rate_stbd=None, bhp_psi=None, rw_ft=0.354, skin=0.0):
        self.wells.append(dict(cell=cell, rate=rate_stbd, bhp=bhp_psi, rw=rw_ft, skin=skin))
    def step(self, dt):
        n = self.nx; C = self.vb*self.phi*self.ct/(5.615*dt)
        lo = np.zeros(n); di = np.zeros(n); up = np.zeros(n); b = np.zeros(n)
        for i in range(n):
            tl = self._T(i, i-1) if i > 0 else 0.0; tr = self._T(i, i+1) if i < n-1 else 0.0
            di[i] = C + tl + tr
            if i > 0: lo[i] = -tl
            if i < n-1: up[i] = -tr
            b[i] = C*self.P[i]
        for w in self.wells:
            i = w["cell"]
            if w["rate"] is not None: b[i] += w["rate"]*self.Bo
            elif w["bhp"] is not None: wi = self._wi(w); di[i] += wi; b[i] += wi*w["bhp"]
        self.P = spsolve(diags([lo[1:], di, up[:-1]], [-1, 0, 1], format="csc"), b)
    def run(self, days, dt):
        for _ in range(int(days/dt)): self.step(dt)


# ── OML 58 depletion problem (do not edit) ───────────────────────────────
NX = 50               # grid cells
LENGTH_FT = 5000.0    # reservoir length
WIDTH_FT = 1000.0     # cross-section width
HEIGHT_FT = 50.0      # net pay
K_MD = 20.0           # permeability (md)
PHI = 0.20            # porosity
CT_PSI = 1.0e-5       # total compressibility (1/psi)
MU_CP = 2.0           # oil viscosity (cp)
PI_PSI = 3000.0       # initial reservoir pressure
BO = 1.2              # formation volume factor

WELL_CELL = 0         # producer on the edge cell
BHP_PSI = 1000.0      # constant bottom-hole pressure
SIM_DAYS = 180.0      # depletion period
DT_LIST = [1, 10, 30, 90]  # time steps to compare (days)


def build_model():
    """Fresh Reservoir1D at initial pressure PI_PSI (call once per run)."""
    return Reservoir1D(NX, LENGTH_FT, WIDTH_FT, HEIGHT_FT, K_MD, PHI,
                       CT_PSI, MU_CP, PI_PSI, Bo=BO)


def final_pavg(dt_days):
    """Mean reservoir pressure (psi) after SIM_DAYS, run at time step dt_days."""
    m = build_model()
    m.add_well(cell=WELL_CELL, bhp_psi=BHP_PSI)
    m.run(SIM_DAYS, dt_days)
    return float(np.mean(m.P))


def sweep(dt_list):
    """Dict mapping each dt_days -> its final mean pressure (psi)."""
    return {dt: final_pavg(dt) for dt in dt_list}


pavg_by_dt = sweep(DT_LIST)

print("mean pressure by dt (days):", pavg_by_dt)

lockCopying code is a Full Access feature.