An aquifer under an oil column

examples/aquifer.py

"""An aquifer under an oil column: the same producer with and without water support.

A Fetkovich aquifer connects to the bottom face of the grid. The model is run twice, once as
stated and once with the aquifer removed by its address (``m.set("aquifers", [])``), and the
field pressure shows what the support is worth.

Flow keywords it exercises: AQUFETP and AQUANCON, plus the oil-water cards (PVDO, PVTW, SWOF,
EQUIL, WCONPROD).

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

import petroscript as ps

support = ps.aquifers.Fetkovich(
    datum_depth=2100.0,
    initial_water_volume=5.0e8,            # sm3
    compressibility=1.0e-4,                # 1/bar: water plus rock
    productivity_index=500.0,              # sm3/day/bar
    connections=[ps.aquifers.Connection(cells=ps.Box(i=(1, 10), j=(1, 10), k=(3, 3)),
                                        face="plus_k")],        # the grid's bottom face
)

m = ps.Model(
    name="AQUIFER",
    units="metric",
    start="2025-01-01",
    grid=ps.grid.Cartesian(shape=(10, 10, 3), dx=100.0, dy=100.0, dz=10.0, top=2000.0),
    rock=ps.Rock(porosity=0.2, permeability=200.0, vertical_ratio=0.1,
                 compressibility=5.0e-5, reference_pressure=200.0),
    fluid=ps.fluid.OilWater(
        oil=ps.pvt.DeadOil(pressure=[50.0, 200.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=200.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=200.0,
                                   oil_water_contact=2100.0),
    aquifers=[support],
    wells=[ps.Well(name="PROD", i=5, j=5)],
    schedule=ps.Schedule(
        end="2027-01-01",
        report="monthly",
        events=[ps.Event("2025-01-01",
                         ps.Perforate(well="PROD", k=(1, 1), diameter=0.2),
                         ps.Produce(well="PROD", oil_rate=300.0, bhp_limit=50.0))],
    ),
    summary=ps.Summary(field=["oil_production_rate", "water_production_rate", "pressure"]),
)

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

without = m.set("aquifers", [])
for name, model in (("with", m), ("without", without)):
    deck = ps.compile(model)
    deck.write(f"out/aquifer/{name}")
    if ps.simulators():
        run = deck.run(folder=f"out/aquifer/{name}/run")
        pressure = run.field("pressure")
        first, last, unit = pressure["value"][0], pressure["value"][-1], pressure["unit"][-1]
        print(f"{name} the aquifer: field pressure {first:.1f} -> {last:.1f} {unit}")

if not ps.simulators():
    print("no simulator found (set PETROSCRIPT_FLOW, or put flow on PATH); the decks are in "
          "out/aquifer")