A unit slip caught by m.check()

examples/gas_water_check.py

"""A gas field over water, where ``m.check()`` catches a unit slip before anything runs.

The gas formation volume factors below were computed with the field-units formula, so they are
in rb/scf, and typed in as rm3/sm3: 5.6146 times too small. The model is well-formed (the table is
positive and falls with pressure), so construction accepts it, but the physics checks read what it
implies: a z-factor near 0.18, outside what any real gas has, and the row names the suspected
unit pair. The row carries its address, so the fix is one ``m.set`` there.

Flow keywords it exercises: GAS WATER, PVDG, PVTW, SGFN, SWFN, EQUIL with a gas-water contact,
WCONPROD on a gas rate, and the field and well gas summaries.

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

import petroscript as ps

m = ps.Model(
    name="GAS_WATER",
    units="metric",
    start="2025-01-01",
    grid=ps.grid.Cartesian(shape=(10, 10, 5), dx=200.0, dy=200.0, dz=20.0, top=2400.0),
    rock=ps.Rock(porosity=0.22, permeability=150.0, vertical_ratio=0.1,
                 compressibility=4.5e-5, reference_pressure=250.0),
    fluid=ps.fluid.GasWater(
        gas=ps.pvt.DryGas(pressure=[50.0, 150.0, 250.0, 350.0],                # barsa
                          gas_fvf=[0.004346, 0.001460, 0.000926, 0.000695],    # rb/scf, by slip
                          gas_viscosity=[0.0135, 0.0155, 0.0185, 0.0215]),     # cP
        water=ps.pvt.Water(reference_pressure=250.0, fvf=1.02, compressibility=4.4e-5,
                           viscosity=0.3, viscosibility=0.0),
        gas_surface_density=0.78,
        water_surface_density=1000.0,
        temperature=90.0,                  # degC: the checks read it to judge the z-factor
    ),
    relperm=[
        ps.relperm.Gas(sg=[0.0, 0.4, 0.8], krg=[0.0, 0.35, 1.0], pc=[0.0, 0.0, 0.0]),
        ps.relperm.Water(sw=[0.2, 0.6, 1.0], krw=[0.0, 0.2, 1.0], pc=[0.0, 0.0, 0.0]),
    ],
    initial=ps.initial.Equilibrium(datum_depth=2420.0, datum_pressure=250.0,
                                   gas_water_contact=2450.0),
    wells=[ps.Well(name="G1", i=5, j=5)],
    schedule=ps.Schedule(
        end="2027-01-01",
        report="monthly",
        events=[ps.Event("2025-01-01",
                         ps.Perforate(well="G1", k=(1, 2), diameter=0.2),
                         ps.Produce(well="G1", gas_rate=500_000.0, bhp_limit=50.0))],
    ),
    summary=ps.Summary(field=["gas_production_rate", "pressure"],
                       wells={"G1": ["bottom_hole_pressure"]}),
)

for issue in m.check():                    # the slip shows up here, with its address
    print(issue.severity, issue.address, issue)

gas = m.get("fluid/1/gas")
m = m.set("fluid/1/gas", gas.replace(gas_fvf=[0.0244, 0.0082, 0.0052, 0.0039]))   # rm3/sm3
assert [issue for issue in m.check() if issue.severity != "info"] == []

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

if ps.simulators():
    run = deck.run(folder="out/gas_water_check/run")
    rate, pressure = run.field("gas_production_rate"), run.field("pressure")
    print(f"gas rate at the end: {rate['value'][-1]:.0f} {rate['unit'][-1]}; pressure "
          f"{pressure['value'][0]:.1f} -> {pressure['value'][-1]:.1f} {pressure['unit'][-1]}")
else:
    print("no simulator found (set PETROSCRIPT_FLOW, or put flow on PATH); the deck is in "
          "out/gas_water_check")