#!/usr/bin/env python3 """ pi_toolkit.py — Productivity Index toolkit for the WPP Karama KRM-4 dataset. Course 01 · Module 02 · Topic 2.2 — The Productivity Index. Reference implementation behind Lecture 2.2 and Field Demonstration L2-2F. All worked numbers are anchored on the CANONICAL KRM-4 dataset from the module technical brief (TECH-BRIEF-m02.md): Two-rate PI test (canonical, Intro-hub data pack): Q1 = 350 STB/d @ Pwf = 4,267 psia (drawdown 583 psi) Q2 = 620 STB/d @ Pwf = 3,817 psia (drawdown 1,033 psi) P-bar = 4,850 psia (build-up anchor) => J = 0.600 STB/d/psi (consistent across rates => linear IPR) Theoretical (PSS, S = 0): Jideal = 0.926 STB/d/psi Flow efficiency: FE = 0.600 / 0.926 = 0.648 (~0.65) Linear AOFP (Jideal): AOFP = 0.926 * 4,850 = 4,491 STB/d NOTE ON THE TWO-RATE CONFLICT (carried as published; see README.md): Topic 2.2's own worked example and Simulator 2 use a DIFFERENT, non-canonical test (Q1 = 750 @ 4,167, Q2 = 1,400 @ 3,586 -> J ~ 1.108). The decks and this toolkit use the canonical 350/620 dataset (J = 0.600) throughout. The non-canonical set is provided only as a documented alternative in this module (it is NOT the KRM-4 PI). Units: oilfield. Q [STB/d], pressure [psia], J [STB/d/psi], k [mD], h [ft], mu [cp], B [RB/STB], skin S dimensionless. """ from __future__ import annotations import math from dataclasses import dataclass # -------------------------------------------------------------------------- # Canonical KRM-4 parameters (TECH-BRIEF-m02 §1.2–1.6) # -------------------------------------------------------------------------- @dataclass(frozen=True) class KRM4: k: float = 18.0 # mD (build-up confirmed) h: float = 95.0 # ft net pay, fully perforated mu: float = 1.4 # cp at P-bar B: float = 1.25 # RB/STB at P-bar (use Bo, NOT Bob = 1.28) re: float = 1320.0 # ft drainage radius (160-acre spacing) rw: float = 0.354 # ft wellbore radius (8.5-in hole) S: float = 5.0 # canonical measured skin (+5) Pbar: float = 4850.0 # psia current average reservoir pressure Pb: float = 3650.0 # psia bubble point C: float = 0.00708 # Darcy radial-flow field constant KRM4_DEFAULT = KRM4() # -------------------------------------------------------------------------- # Core functions # -------------------------------------------------------------------------- def j_from_two_rate(q1: float, pwf1: float, q2: float, pwf2: float) -> float: """Two-rate (slope) PI: J = (Q2 - Q1) / (Pwf1 - Pwf2). Independent of average reservoir pressure. Order-agnostic in the pair. """ return (q2 - q1) / (pwf1 - pwf2) def j_single_point(q: float, pbar: float, pwf: float) -> float: """Single-point PI: J = Q / (P-bar - Pwf). Depends on the P-bar anchor.""" return q / (pbar - pwf) def j_theoretical(k=KRM4_DEFAULT.k, h=KRM4_DEFAULT.h, mu=KRM4_DEFAULT.mu, B=KRM4_DEFAULT.B, re=KRM4_DEFAULT.re, rw=KRM4_DEFAULT.rw, S=0.0, C=KRM4_DEFAULT.C) -> float: """Theoretical PSS productivity index: J = C * k * h / [mu * B * (ln(re/rw) - 0.75 + S)] With S = 0 this returns Jideal (0.926 for canonical KRM-4 inputs). """ geom = math.log(re / rw) - 0.75 + S return C * k * h / (mu * B * geom) def skin_from_j(jmeas: float, k=KRM4_DEFAULT.k, h=KRM4_DEFAULT.h, mu=KRM4_DEFAULT.mu, B=KRM4_DEFAULT.B, re=KRM4_DEFAULT.re, rw=KRM4_DEFAULT.rw, C=KRM4_DEFAULT.C) -> float: """Back-calculate skin implied by a measured J: ln(re/rw) - 0.75 + S = C * k * h / (Jmeas * mu * B) S = that bracket - (ln(re/rw) - 0.75) """ bracket = C * k * h / (jmeas * mu * B) return bracket - (math.log(re / rw) - 0.75) def flow_efficiency(jmeas: float, jideal: float) -> float: """FE = Jactual / Jideal. Equivalent to (ln-0.75) / (ln-0.75 + S).""" return jmeas / jideal def aofp(j: float, pbar: float) -> float: """Absolute open-flow potential of the LINEAR IPR: AOFP = J * P-bar. Valid only if the whole drawdown range stays above the bubble point; it is an index for comparing wells, not a deliverable rate. """ return j * pbar def required_pwf(j: float, pbar: float, q_target: float) -> float: """Flowing pressure needed to deliver a target rate: Pwf = P-bar - Q/J.""" return pbar - q_target / j def rate_at_pwf(j: float, pbar: float, pwf: float) -> float: """Linear IPR rate at a given Pwf: Q = J * (P-bar - Pwf).""" return j * (pbar - pwf) def ipr_linear(j: float, pbar: float, n: int = 11): """Linear IPR table from Pwf = P-bar (Q = 0) down to Pwf = 0 (Q = AOFP). Returns a list of (Pwf, Q) pairs, n evenly spaced Pwf points. """ rows = [] for i in range(n): pwf = pbar * (1 - i / (n - 1)) rows.append((pwf, j * (pbar - pwf))) return rows def stimulation_prize(jpre: float, jpost: float, pbar: float, pwf: float): """Rate uplift from removing damage at fixed operating Pwf. Returns (q_pre, q_post, delta_q, pct_uplift). """ q_pre = rate_at_pwf(jpre, pbar, pwf) q_post = rate_at_pwf(jpost, pbar, pwf) dq = q_post - q_pre return q_pre, q_post, dq, 100.0 * dq / q_pre def stimulation_economics(delta_q: float, price: float = 75.0, opex: float = 12.0, job_cost: float = 800_000.0): """Headline stimulation economics. Returns (net_per_bbl, revenue_per_day, revenue_per_year, payback_days). """ net = price - opex rev_day = delta_q * net rev_yr = rev_day * 365.0 payback = job_cost / rev_day return net, rev_day, rev_yr, payback # -------------------------------------------------------------------------- # __main__ demo — reproduces the brief's worked numbers # -------------------------------------------------------------------------- def _approx(a: float, b: float, tol: float) -> str: return "OK" if abs(a - b) <= tol else f"** MISMATCH (got {a:.4f}, want {b}) **" if __name__ == "__main__": k = KRM4_DEFAULT print("=" * 68) print(" KRM-4 PRODUCTIVITY INDEX — pi_toolkit.py verification") print("=" * 68) # --- Canonical two-rate test --- q1, pwf1 = 350.0, 4267.0 q2, pwf2 = 620.0, 3817.0 dp1, dp2 = k.Pbar - pwf1, k.Pbar - pwf2 j_slope = j_from_two_rate(q1, pwf1, q2, pwf2) j_p1 = j_single_point(q1, k.Pbar, pwf1) j_p2 = j_single_point(q2, k.Pbar, pwf2) print("\n[1] CANONICAL TWO-RATE TEST (350/620 -> J = 0.600)") print(f" rate 1 : {q1:.0f} STB/d @ {pwf1:.0f} psia (drawdown {dp1:.0f} psi)") print(f" rate 2 : {q2:.0f} STB/d @ {pwf2:.0f} psia (drawdown {dp2:.0f} psi)") print(f" slope J = (620-350)/(4267-3817) = 270/450 = {j_slope:.4f} {_approx(j_slope, 0.600, 1e-6)}") print(f" point 1 : 350/{dp1:.0f} = {j_p1:.4f} {_approx(j_p1, 0.600, 5e-4)}") print(f" point 2 : 620/{dp2:.0f} = {j_p2:.4f} {_approx(j_p2, 0.600, 5e-4)}") jmeas = 0.600 # --- Theoretical J --- geom = math.log(k.re / k.rw) numer = k.C * k.k * k.h denom = k.mu * k.B * (geom - 0.75) jideal = j_theoretical(S=0.0) print("\n[2] THEORETICAL Jideal (PSS, S = 0)") print(f" ln(re/rw) = {geom:.4f} (brief 8.224)") print(f" ln(re/rw) - 0.75 = {geom - 0.75:.4f} (brief 7.474)") print(f" numerator 0.00708kh = {numer:.4f} (brief 12.11)") print(f" denom mu*B*7.474 = {denom:.4f} (brief 13.08)") print(f" Jideal = {jideal:.4f} {_approx(jideal, 0.926, 2e-3)}") # --- Skin + FE --- s_implied = skin_from_j(jmeas) fe = flow_efficiency(jmeas, jideal) print("\n[3] DIAGNOSIS FROM Jmeas = 0.600") print(f" implied skin S = {s_implied:.2f} (brief +4.06 ~ +4; canonical +5)") print(f" flow efficiency FE = {fe:.3f} {_approx(fe, 0.648, 5e-3)} (~0.65)") # --- AOFP --- aof_ideal = aofp(jideal, k.Pbar) aof_meas = aofp(jmeas, k.Pbar) print("\n[4] LINEAR AOFP (= J * P-bar)") print(f" AOFP (Jideal 0.926) = {aof_ideal:.0f} STB/d {_approx(aof_ideal, 4491, 8)}") print(f" AOFP (Jmeas 0.600) = {aof_meas:.0f} STB/d {_approx(aof_meas, 2910, 5)}") # --- Target & stimulation prize --- pwf_op = 3100.0 q_target = 1200.0 q_dam = rate_at_pwf(jmeas, k.Pbar, pwf_op) q_clean = rate_at_pwf(jideal, k.Pbar, pwf_op) qpre, qpost, dq, pct = stimulation_prize(jmeas, jideal, k.Pbar, pwf_op) print(f"\n[5] OPERATING POINT (Pwf = {pwf_op:.0f} psia, target {q_target:.0f} STB/d)") print(f" Q damaged (J 0.600) = {q_dam:.0f} STB/d {_approx(q_dam, 1050, 1)} (misses target by 150)") print(f" Q clean (J 0.926) = {q_clean:.0f} STB/d {_approx(q_clean, 1621, 2)} (meets target)") print(f" stimulation prize dQ = +{dq:.0f} STB/d ({pct:+.1f}%) {_approx(dq, 571, 2)}") print(f" required Pwf for target (J 0.600) = {required_pwf(jmeas, k.Pbar, q_target):.0f} psia") print(f" required Pwf for target (J 0.926) = {required_pwf(jideal, k.Pbar, q_target):.0f} psia") # --- Economics --- # Brief carries the canonical stimulation prize as +571 STB/d (rounded # 0.326*1,750); the headline economics ($35,973/day, $13.13 MM/yr) are # quoted on that figure, so reproduce them on +571 rather than the # slightly-rounded +570 that falls out of Jideal = 0.9257. dq_canonical = 571.0 net, rev_day, rev_yr, payback = stimulation_economics(dq_canonical) print(f"\n[6] STIMULATION ECONOMICS (prize +571 STB/d, $75/STB, $12 opex, $800k job)") print(f" net per bbl = ${net:.0f}") print(f" revenue / day = ${rev_day:,.0f} {_approx(rev_day, 35973, 50)}") print(f" revenue / year = ${rev_yr/1e6:.2f} MM (brief ~$13.13 MM)") print(f" payback = {payback:.1f} days {_approx(payback, 22.2, 0.5)}") # --- Linear IPR table head --- print("\n[7] LINEAR IPR (Jideal = 0.926, P-bar = 4,850) — selected points") for pwf in (4850, 4500, 4000, 3650, 3100, 2000, 0): print(f" Pwf {pwf:>5.0f} psia -> Q = {rate_at_pwf(jideal, k.Pbar, pwf):>5.0f} STB/d") print("\n" + "=" * 68) print(" All canonical KRM-4 PI numbers reproduced.") print("=" * 68)