Shut a well on water cut

examples/water_cut_rule.py

"""A rule that fires during the run: shut a producer once its water cut passes 50%.

An injector floods towards two producers. P1 sits close to the injector and waters out first;
the rule ``ps.actions.When(ps.actions.well(p1).water_cut > 0.5, then=[ps.Shut(well="P1")])``
shuts it, and P2 keeps producing. Flow tests the condition at the end of every time step, not
only at report dates, and the summary holds every time step, so the step whose water cut first
passes 0.5 is when the rule fires. The rule is a dated change like any other: it applies from its
``ps.Event``'s date on.

Flow keywords it exercises: ACTIONX with a WWCT condition, WELOPEN SHUT inside it, plus the
waterflood's cards (PVDO, PVTW, SWOF, EQUIL, WCONPROD, WCONINJE).

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

import petroscript as ps

p1 = ps.Well(name="P1", i=3, j=3)
p2 = ps.Well(name="P2", i=10, j=10)
shut_p1 = ps.actions.When(ps.actions.well(p1).water_cut > 0.5,
                          then=[ps.Shut(well="P1")], name="SHUT_P1")

m = ps.Model(
    name="WATER_CUT_RULE",
    units="metric",
    start="2025-01-01",
    grid=ps.grid.Cartesian(shape=(10, 10, 1), dx=50.0, dy=50.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),
    wells=[ps.Well(name="INJ", i=1, j=1), p1, p2],
    schedule=ps.Schedule(
        end="2028-01-01",
        report="monthly",
        events=[
            ps.Event("2025-01-01",
                     ps.Perforate(well="INJ", k=(1, 1), diameter=0.2),
                     ps.Perforate(well="P1", k=(1, 1), diameter=0.2),
                     ps.Perforate(well="P2", k=(1, 1), diameter=0.2),
                     ps.Inject(well="INJ", phase="water", rate=200.0, bhp_limit=400.0),
                     ps.Produce(well="P1", liquid_rate=100.0, bhp_limit=100.0),
                     ps.Produce(well="P2", liquid_rate=100.0, bhp_limit=100.0),
                     shut_p1),
        ],
    ),
    summary=ps.Summary(field=["oil_production_rate"],
                       wells={"*": ["oil_production_rate", "water_cut"]}),
)

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

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

if ps.simulators():
    run = deck.run(folder="out/water_cut_rule/run")
    cut = run.wells("water_cut").filter(well="P1")
    fired = cut.filter(cut["value"] > 0.5)      # the condition, at the end of each time step
    if len(fired):
        print(f"the rule shuts P1 on {fired['date'][0].date()}, "
              f"when its water cut reaches {fired['value'][0]:.2f}")
    else:
        print("the rule did not fire: P1's water cut stays at or below 0.5")
    oil = run.wells("oil_production_rate")
    p1, p2 = (oil.filter(well=name) for name in ("P1", "P2"))
    shut_rate = abs(p1["value"][-1])                    # a shut well's rate is written as -0.0
    print(f"oil rate at the end: P1 {shut_rate:.0f} {p1['unit'][-1]}, "
          f"P2 {p2['value'][-1]:.0f} {p2['unit'][-1]}")
else:
    print("no simulator found (set PETROSCRIPT_FLOW, or put flow on PATH); the deck is in "
          "out/water_cut_rule")