SPE1 walkthrough

The first SPE comparative solution project (Odeh, 1981) is a 10 x 10 x 3 black-oil box of undersaturated oil: gas is injected at one corner and oil produced at the opposite one. opm-tests/spe1/SPE1CASE1.DATA is its deck. This page states the same case as one ps.Model, checks it, compiles it and runs it. Every Python block runs, in order, in one fresh process (tests/test_docs_run.py runs them).

Semantic decomposition

Problem concept	Where the model states it
Run identity and units	`ps.Model(name=, units="field", start=)`
A 10 x 10 x 3 box, three layers of 20, 30 and 50 ft	`grid=ps.grid.Cartesian(shape=, dx=, dy=, dz=, top=)`
Porosity, layered permeability, rock compressibility	`rock=ps.Rock(...)`
Live oil, dry gas and water, with surface densities	`fluid=ps.fluid.BlackOil(oil=ps.pvt.LiveOil, gas=ps.pvt.DryGas, water=ps.pvt.Water)`
Oil-water and gas-oil relative permeability	`relperm=[ps.relperm.OilWater, ps.relperm.GasOil]`
Pressure at a datum, two contacts, dissolved gas by depth	`initial=ps.initial.Equilibrium(...)`
A producer and a gas injector in one group	`wells=[ps.Well(...)]`, `groups=[ps.Group(...)]`
Their perforations and controls, from the start	`ps.Event("2015-01-01", ps.Perforate, ps.Produce, ps.Inject)`
No re-solution of gas into the oil	`ps.pvt.DissolutionLimits(maximum_rate=0.0)` in the start’s event
Ten years of report steps	`ps.Schedule(report=[ps.Every(...)])`
Field and well results	`summary=ps.Summary(...)`

The model

Bare numbers are in field units: psia, ft, mD, stb/day and Mscf/day; the PVT tables’ solution gas-oil ratio is in Mscf/stb and the gas formation volume factor in rb/Mscf.

import numpy as np

import petroscript as ps

layers = np.ones((10, 10, 3))                     # one value per cell, [i-1, j-1, k-1]
thickness = layers * [20.0, 30.0, 50.0]           # ft, the top layer first
permeability = layers * [500.0, 50.0, 200.0]      # mD, the same along x, y and z
sw = [0.12, 0.18, 0.24, 0.30, 0.36, 0.42, 0.48, 0.54, 0.60, 0.66, 0.72, 0.78, 0.84, 0.91, 1.0]
krw = [0.0, 4.64876033057851e-08, 1.86e-07, 4.18388429752066e-07, 7.43801652892562e-07,
       1.16219008264463e-06, 1.67355371900826e-06, 2.27789256198347e-06, 2.97520661157025e-06,
       3.7654958677686e-06, 4.64876033057851e-06, 5.625e-06, 6.69421487603306e-06,
       8.05914256198347e-06, 1e-05]
sg = [0.0, 0.001, 0.02, 0.05, 0.12, 0.2, 0.25, 0.3, 0.4, 0.45, 0.5, 0.6, 0.7, 0.85, 0.88]
krg = [0.0, 0.0, 0.0, 0.005, 0.025, 0.075, 0.125, 0.19, 0.41, 0.6, 0.72, 0.87, 0.94, 0.98, 0.984]
kro = [1.0, 1.0, 0.997, 0.98, 0.7, 0.35, 0.2, 0.09, 0.021, 0.01, 0.001, 0.0001, 0.0, 0.0, 0.0]
no_pc = [0.0] * 15
g1 = ps.Group(name="G1")

spe1 = ps.Model(
    name="SPE1CASE1",
    units="field",
    start="2015-01-01",
    grid=ps.grid.Cartesian(shape=(10, 10, 3), dx=1000.0, dy=1000.0, dz=thickness, top=8325.0),
    rock=ps.Rock(porosity=0.3, permeability_x=permeability, permeability_y=permeability,
                 permeability_z=permeability, compressibility=3e-06, reference_pressure=14.7),
    fluid=ps.fluid.BlackOil(
        oil=ps.pvt.LiveOil(        # long form: each solution GOR's saturated row first
            solution_gor=[0.001, 0.0905, 0.18, 0.371, 0.636, 0.775, 0.93, 1.27, 1.27, 1.618,
                          1.618],
            pressure=[14.7, 264.7, 514.7, 1014.7, 2014.7, 2514.7, 3014.7, 4014.7, 9014.7,
                      5014.7, 9014.7],
            oil_fvf=[1.062, 1.15, 1.207, 1.295, 1.435, 1.5, 1.565, 1.695, 1.579, 1.827, 1.737],
            oil_viscosity=[1.04, 0.975, 0.91, 0.83, 0.695, 0.641, 0.594, 0.51, 0.74, 0.449,
                           0.631],
        ),
        gas=ps.pvt.DryGas(
            pressure=[14.7, 264.7, 514.7, 1014.7, 2014.7, 2514.7, 3014.7, 4014.7, 5014.7, 9014.7],
            gas_fvf=[166.666, 12.093, 6.274, 3.197, 1.614, 1.294, 1.08, 0.811, 0.649, 0.386],
            gas_viscosity=[0.008, 0.0096, 0.0112, 0.014, 0.0189, 0.0208, 0.0228, 0.0268, 0.0309,
                           0.047],
        ),
        water=ps.pvt.Water(reference_pressure=4017.55, fvf=1.038, compressibility=3.22e-06,
                           viscosity=0.318, viscosibility=0.0),
        oil_surface_density=53.66,      # lb/ft3
        gas_surface_density=0.0533,
        water_surface_density=64.49,
    ),
    relperm=[
        ps.relperm.OilWater(sw=sw, krw=krw, krow=kro, pcow=no_pc),
        ps.relperm.GasOil(sg=sg, krg=krg, krog=kro, pcog=no_pc),
    ],
    initial=ps.initial.Equilibrium(
        datum_depth=8400.0, datum_pressure=4800.0,
        oil_water_contact=8450.0, gas_oil_contact=8300.0,
        solution_gor_by_depth=ps.initial.RsVsDepth(depth=[8300.0, 8450.0],
                                                   solution_gor=[1.27, 1.27]),
        accuracy=0.0,     # EQUIL item 9 as SPE1 states it: each cell from its centre
    ),
    wells=[ps.Well(name="PROD", i=10, j=10, group=g1, bhp_reference_depth=8400.0),
           ps.Well(name="INJ", i=1, j=1, group=g1, bhp_reference_depth=8335.0)],
    groups=[g1],
    schedule=ps.Schedule(
        report=[ps.Every(*[31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] * 10)],   # days
        events=[ps.Event("2015-01-01",
                         ps.Perforate(well="PROD", k=(3, 3), diameter=0.5),
                         ps.Perforate(well="INJ", k=(1, 1), diameter=0.5),
                         ps.Produce(well="PROD", oil_rate=20000.0, bhp_limit=1000.0),
                         ps.Inject(well="INJ", phase="gas", rate=100000.0, bhp_limit=9014.0),
                         ps.pvt.DissolutionLimits(maximum_rate=0.0))],
    ),
    summary=ps.Summary(field=["oil_production_rate", "gas_oil_ratio"],
                       wells={"*": ["bottom_hole_pressure", "oil_production_rate",
                                    "gas_injection_rate"]}),
)

A layered field is an array. dz= and the permeabilities take one value per cell, indexed [i - 1, j - 1, k - 1]; a uniform one is a bare number (porosity=0.3).
The report steps are the deck’s. SPE1 reports at steps of 31, 28, 31, … days, a 28-day February every year, so ps.Every states the lengths rather than "monthly".
The control rule. In ps.Produce the bare quantity is the target and every other one a limit: the producer holds 20,000 stb/day of oil while its bottom-hole pressure stays above 1,000 psia, stated at its reference depth of 8,400 ft.
The preferred phase comes from the control. ps.Well states no phase here; the compile writes each well’s from its first control (oil for PROD, gas for INJ).
EQUIL item 9. accuracy=0.0 initializes each cell from its centre, as the deck states. Left unstated, the card defaults it and Flow applies its own default, -5: each cell’s initial pressures and saturations are averaged over ten horizontal slices of it (2 × 5).

Check and compile

rows = spe1.check()
assert [row for row in rows if row.severity == "error"] == []
assert spe1.get("wells/PROD").bhp_reference_depth == 8400.0

deck = ps.compile(spe1)
text = deck.files["SPE1CASE1.DATA"]
assert "\nDISGAS\n" in text           # a live oil: the gas dissolves
assert "DRSDT" in text                # the re-solution limit
deck.write("spe1")                    # spe1/SPE1CASE1.DATA

The deck in the corpus is the same reservoir. ps.load reads it into a model whose grid, rock, fluid, saturation functions and groups equal the ones typed above:

benchmark = ps.load("opm-tests/spe1/SPE1CASE1.DATA")
for part in ("grid", "rock", "fluid", "relperm", "groups"):
    assert getattr(spe1, part) == getattr(benchmark, part), part
rest = {change.address.split("/")[0] for change in spe1.diff(benchmark)}
assert rest == {"title", "capacity", "dimensions", "summary", "wells", "schedule"}

The rest of spe1.diff(benchmark) is what the deck states and the model above leaves to its defaults: the title, the sizes, its wider summary request, and the values that equal the defaults (the preferred phases, and each perforation’s column and open status).

Run it and read the results

run = deck.run(folder="spe1/run")            # OPM Flow, found as ps.simulators() lists it
oil = run.field("oil_production_rate")       # polars: date, value, unit
assert oil["unit"][0] == "stb/day"
assert oil["value"][0] == 20000.0            # the producer holds its target at first
bhp = run.wells("bottom_hole_pressure")       # date, well, value, unit
assert set(bhp["well"]) == {"INJ", "PROD"}
assert oil["value"][-1] < 20000.0            # ten years on, the target is no longer met

The oil rate holds at 20,000 stb/day while the producer’s pressure allows it. The injected gas breaks through, the producing gas-oil ratio climbs, and once the bottom-hole pressure reaches its 1,000 psia limit the well holds that pressure and the oil rate falls.