A hand-built lift table

examples/lift_table.py

"""A producer on a lift table, held to a tubing-head pressure limit as the reservoir depletes.

The table gives the bottom-hole pressure needed to lift a liquid rate to the surface at a given
tubing-head pressure and water cut: the hydrostatic head of the column plus friction. It is built
here with NumPy; in practice it comes from a well-performance program. The producer targets
400 sm3/day of liquid with a 20 barsa tubing-head limit; as the pressure falls, the limit takes
over and the rate declines.

The table's axes are named fields. ``bhp`` is shaped (thp, water fraction, gas fraction, artificial
lift, rate), and the control names the table itself (``lift_table=table``), never its number.

Flow keywords it exercises: VFPPROD, WCONPROD with a THP limit and its table number, and the
well's THP summary.

Run it with ``python examples/lift_table.py``; it writes ``out/lift_table/``.
"""

import numpy as np
import polars as pl

import petroscript as ps

thp = np.array([10.0, 20.0, 40.0, 80.0])                 # barsa: up past the THP it starts at
water_cut = np.array([0.0, 0.5, 0.9])
rate = np.array([50.0, 200.0, 400.0, 800.0])             # sm3/day of liquid
head = 2000.0 * 9.81e-5 * (850.0 + (1020.0 - 850.0) * water_cut)   # barsa over 2000 m
friction = 2.0e-5 * rate**2                              # barsa
bhp = (thp[:, None, None] + head[None, :, None] + friction[None, None, :])[:, :, None, None, :]
table = ps.lift.ProductionTable(
    number=1, bhp_reference_depth=2000.0, rate_kind="liquid", rate=rate, thp=thp,
    water_fraction=water_cut, gas_fraction=np.array([0.0]), artificial_lift=np.array([0.0]),
    bhp=bhp,
)

m = ps.Model(
    name="LIFT_TABLE",
    units="metric",
    start="2025-01-01",
    grid=ps.grid.Cartesian(shape=(10, 10, 3), dx=300.0, dy=300.0, dz=20.0, top=2000.0),
    rock=ps.Rock(porosity=0.2, permeability=300.0, vertical_ratio=0.1,
                 compressibility=5.0e-5, reference_pressure=230.0),
    fluid=ps.fluid.OilWater(
        oil=ps.pvt.DeadOil(pressure=[50.0, 230.0, 400.0], oil_fvf=[1.12, 1.10, 1.08],
                           oil_viscosity=[1.8, 2.0, 2.3]),
        water=ps.pvt.Water(reference_pressure=230.0, fvf=1.01, compressibility=4.5e-5,
                           viscosity=0.4, viscosibility=0.0),
        oil_surface_density=850.0,
        water_surface_density=1020.0,
    ),
    relperm=[ps.relperm.OilWater(sw=[0.2, 0.4, 0.6, 0.8], krw=[0.0, 0.05, 0.2, 0.5],
                                 krow=[0.9, 0.4, 0.1, 0.0], pcow=[0.0, 0.0, 0.0, 0.0])],
    initial=ps.initial.Equilibrium(datum_depth=2000.0, datum_pressure=230.0,
                                   oil_water_contact=2100.0),
    wells=[ps.Well(name="PROD", i=5, j=5, bhp_reference_depth=2000.0)],
    lift_tables=[table],
    schedule=ps.Schedule(
        end="2028-01-01",
        report="monthly",
        events=[ps.Event("2025-01-01",
                         ps.Perforate(well="PROD", k=(1, 3), diameter=0.2),
                         ps.Produce(well="PROD", liquid_rate=400.0, bhp_limit=50.0,
                                    thp_limit=20.0, lift_table=table))],
    ),
    summary=ps.Summary(field=["oil_production_rate", "pressure"],
                       wells={"PROD": ["liquid_production_rate", "tubing_head_pressure",
                                       "bottom_hole_pressure"]}),
)

for issue in m.check():
    print(issue.severity, issue.address, issue)

deck = ps.compile(m)
deck.write("out/lift_table")

if ps.simulators():
    run = deck.run(folder="out/lift_table/run")
    liquid = run.wells("liquid_production_rate")
    held = run.wells("tubing_head_pressure").filter(pl.col("value") <= 20.0 + 1e-3)
    print("the THP limit takes over on", held["date"][0].date() if len(held) else "no date yet")
    print(f"liquid rate {liquid['value'][0]:.0f} -> {liquid['value'][-1]:.0f} {liquid['unit'][-1]}")
else:
    print("no simulator found (set PETROSCRIPT_FLOW, or put flow on PATH); the deck is in "
          "out/lift_table")