Build a model from a word problem

A model starts as a description in words. This page turns one into a PetroScript model, one question at a time, then checks it, writes the deck and runs it. Do not start from deck keywords: start from what the reservoir is, and say each fact with the keyword argument that holds it.

The problem

A square oil reservoir, 1,500 m on a side and 30 m thick in three equal layers, its top at 2,000 m. Porosity 0.22, horizontal permeability 200 mD, vertical a tenth of it; rock compressibility 5e-5 1/bar at 200 bar. Dead oil: formation volume factor 1.10 rm3/sm3 and viscosity 2 cP at 200 bar, surface density 850 kg/m3; water: 0.4 cP, 1020 kg/m3. Connate water 0.2, residual oil 0.25, Corey exponents 2, end-point water relative permeability 0.4. The pressure is 210 bar at 2,000 m, and the oil-water contact is at 2,040 m. A producer in the centre makes 1,000 sm3/day of oil, its bottom-hole pressure kept above 100 bar; a water injector in a corner takes 1,200 sm3/day below 300 bar. Run three years from 2025, report monthly, and record the field’s oil rate, water cut and pressure, and the producer’s bottom-hole pressure.

The questions to answer first

Question	What it fixes	Where it goes
What is the run?	name, unit system, start	ps.Model(name=, units=, start=)
What grid?	geometry	grid=ps.grid.Cartesian(...)
What rock?	porosity, permeability, compressibility	rock=ps.Rock(...)
Which fluids?	phases and their PVT	fluid=ps.fluid.OilWater(...) with ps.pvt tables
How do they flow together?	relative permeability	relperm=[ps.relperm.OilWater(...)]
What is the initial state?	datum pressure, contacts	initial=ps.initial.Equilibrium(...)
Which wells, doing what, when?	wells, perforations, controls	wells=[ps.Well(...)], schedule=ps.Schedule(events=[ps.Event(...)])
What do you want to read after?	results	summary=ps.Summary(...)

Fill the table before writing code. If the problem leaves out a value that changes what the reservoir does (a contact, a compressibility, a relative permeability), ask for it rather than invent it. This problem is in METRIC units: m, barsa, sm3/day, rm3/sm3, cP, kg/m3, mD.

Each part, by its keyword

The grid: 15 x 15 cells of 100 m, three layers of 10 m. The model’s parts take keywords.

import numpy as np

import petroscript as ps

grid = ps.grid.Cartesian(shape=(15, 15, 3), dx=100.0, dy=100.0, dz=10.0, top=2000.0)

The rock is three facts: porosity, permeability and compressibility (with the pressure it is measured at). Leave the compressibility out and OPM Flow runs incompressible rock, which m.check() reports as a warning; PetroScript never supplies the number.

rock = ps.Rock(porosity=0.22, permeability=200.0, vertical_ratio=0.1,   # kz = 0.1 kx
               compressibility=5.0e-5, reference_pressure=200.0)        # 1/bar at 200 bar

Two phases, oil and water, make a ps.fluid.OilWater. A table is named columns, one value per row; the dead oil’s formation volume factor and viscosity are stated against pressure:

fluid = ps.fluid.OilWater(
    oil=ps.pvt.DeadOil(pressure=[100.0, 200.0, 300.0],      # barsa
                       oil_fvf=[1.11, 1.10, 1.09],          # rm3/sm3
                       oil_viscosity=[1.9, 2.0, 2.1]),      # cP
    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,                              # kg/m3
    water_surface_density=1020.0,
)

The relative permeability comes from the Corey parameters. Build the columns with NumPy; a table takes lists, tuples or arrays:

sw = np.linspace(0.2, 0.75, 8)                               # connate water to 1 - Sor
s = (sw - 0.2) / (0.75 - 0.2)                                # normalised saturation
relperm = ps.relperm.OilWater(sw=sw, krw=0.4 * s**2, krow=(1.0 - s) ** 2,
                              pcow=np.zeros_like(sw))

The initial state: the pressure at one depth and the contact. PetroScript integrates the rest.

initial = ps.initial.Equilibrium(datum_depth=2000.0, datum_pressure=210.0,
                                 oil_water_contact=2040.0)

A well is an identity: a name and a column. What it does is dated: its perforations and its control are changes a ps.Event states. In ps.Produce and ps.Inject the one quantity given bare is the target, every other a *_limit:

wells = [ps.Well(name="PROD", i=8, j=8), ps.Well(name="INJ", i=1, j=1)]
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.Perforate(well="INJ", k=(1, 3), diameter=0.2),
                     ps.Produce(well="PROD", oil_rate=1000.0, bhp_limit=100.0),
                     ps.Inject(well="INJ", phase="water", rate=1200.0, bhp_limit=300.0))],
)

Results are asked for, never inferred, by the names you will read them by:

summary = ps.Summary(field=["oil_production_rate", "water_cut", "pressure"],
                     wells={"PROD": ["bottom_hole_pressure"]})

The model, checked

m = ps.Model(name="SQUARE", units="metric", start="2025-01-01", grid=grid, rock=rock,
             fluid=fluid, relperm=[relperm], initial=initial, wells=wells,
             schedule=schedule, summary=summary)

for issue in m.check():               # instant, no simulator: errors, warnings, info
    print(issue.severity, issue.address, issue)
assert not [issue for issue in m.check() if issue.severity == "error"]

A model that contradicts itself is refused as it is built, naming the address and the fix: a well with no control ([well.control.missing]), a missing initial state, a table out of order. A warning never blocks; read it, and fix the number at its address (m = m.set("rock/compressibility", 4.0e-5)) or state acknowledged_warnings=[code] when the model means it.

Write it and run it

deck = ps.compile(m)                  # the deck for OPM Flow, in memory
deck.write("square")                  # square/SQUARE.DATA
if ps.simulators():                   # OPM Flow found: PETROSCRIPT_FLOW or flow on PATH
    run = deck.run(folder="square/run")
    print(run.field("oil_production_rate").tail(3))       # date, value, unit

Read a finished run says what else the run holds. To change the model later, keep using verbs: m.at(date, ...), m.set(address, value), x.replace(...) (Edit a loaded deck).