Quickstart

PetroScript holds a reservoir model as one typed, frozen value, ps.Model. You build it by keyword or load it from a deck, check it, compile it for a simulator, run it and read the results by name. This page follows Rosa and the Merle gas field through all of it. Every block below runs, in order, in one fresh Python process (tests/test_quickstart.py runs them).

1. What you will not guess

One frozen value. Nothing is edited in place: every verb (m.at, m.set, x.replace) returns a new value, checked at once. Facts that only hold together (a new well and its control) go in one m.set({...}).
Keywords. The model’s parts take keywords: ps.Model(name="MERLE", ...), ps.Well(name="M-1", i=3, j=3); ps.Model("MERLE") is a TypeError. ps.Event, ps.Every and ps.Quantity take their values in order.
Where a name lives. The root holds the model and its top-level parts: ps.Model, ps.Rock, ps.Regions, ps.Well, ps.Group, ps.Schedule, ps.Event, ps.Summary, ps.Output, ps.Numerics and the dated changes. A Model keyword that takes one of several kinds finds them in the namespace of its own word: grid=ps.grid.Cartesian(...), fluid=ps.fluid.GasWater(...), initial=ps.initial.Equilibrium(...). Tables live in ps.pvt, ps.relperm and ps.lift; ps.wells, ps.groups, ps.aquifers and ps.actions hold the rest by topic.
A well is an identity; what it does is dated. ps.Well(name=, i=, j=) only declares it. Its perforations, controls and status are changes the schedule dates: ps.Event(date, *changes) in ps.Schedule(events=[...]), or the verb m.at(date, *changes).
The control rule. In ps.Produce and ps.Inject the one quantity given bare is the target and every other quantity is a *_limit: ps.Produce(well="M-1", gas_rate=500_000.0, bhp_limit=50.0) holds the gas rate while the bottom-hole pressure stays above 50.
Addresses. Everything a model states has an address: wells/M-1, rock/porosity, fluid/1/gas, initial/1/datum_pressure, schedule/2011-03-01/M-2/control. m.get, m.set, m.find and m.diff read and write by address.
Names, never deck words. Fields and results are spelled out (gas_production_rate, bottom_hole_pressure), the lighter phase first in a pair (gas_water_contact, ps.relperm.GasOil). Textbook symbols stay (dx, krg, pc).
Units. A bare number is in the model’s unit system: units="metric", "field" or "lab". A FIELD gas-oil ratio is Mscf/stb, a METRIC one sm3/sm3; a water cut is a fraction. The units table lists each quantity.
Look it up; never guess.

import petroscript as ps

ps.help(ps.Produce)      # a class: its fields, defaults, the card it writes, its address
ps.help("WCONINJE")      # a deck keyword: the class that writes it
ps.help("FGPR")          # a summary mnemonic: the name its result is read by
ps.help("gas_production_total")  # a result's name: its OPM Flow vector (FGPT) in each scope
ps.help("addresses")     # a topic ("units" and "facts" are the others); ps.help() is the index

A wrong name is refused with the nearest one: [field.unknown] Produce: no field 'bhp_limt'; did you mean 'bhp_limit'?. Most refusals are ps.PetroScriptErrors: ps.ModelError (a model that contradicts itself), ps.CompileError, ps.DeckImportError and ps.RunError among them. Some are plain Python errors (errors.md): AttributeError for a name that does not exist, TypeError for arguments by position, a bare ValueError from ps.lift.generate’s checks. A model’s refusal prints [code] address: message and then fix: with what to type.

2. Rosa builds the model

Dry gas over water in a 10 x 10 x 5 box from 2400 to 2500 m. M-1 produces from June 2009 and M-2 from March 2011.

merle = ps.Model(
    name="MERLE",                      # the written file's stem
    units="metric",
    start="2009-06-01",                # a date's one text, as it prints
    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,   # kx = ky, kz = 0.1 kx
                 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.0244, 0.0082, 0.0052, 0.0039],         # rm3/sm3
                          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 tables' temperature (no card; the checks read it)
    ),
    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=3626.0,
                                   gas_water_contact=2450.0),
    wells=[ps.Well(name="M-1", i=3, j=3),
           ps.Well(name="M-2", i=8, j=8, introduced="2011-03-01")],
    schedule=ps.Schedule(
        end="2011-06-01",
        report="monthly",              # or a list of dates
        events=[
            ps.Event("2009-06-01",
                     ps.Perforate(well="M-1", k=(1, 2), diameter=0.2),
                     ps.Produce(well="M-1", gas_rate=500_000.0, bhp_limit=50.0)),
            ps.Event("2011-03-01",
                     ps.Perforate(well="M-2", k=(1, 2), diameter=0.2),
                     ps.Produce(well="M-2", gas_rate=400_000.0, bhp_limit=50.0)),
        ],
    ),
    summary=ps.Summary(field=["gas_production_rate", "gas_production_total", "pressure"],
                       wells={"*": ["gas_production_rate", "bottom_hole_pressure"]}),
)

The model counts its own sizes. Merle states no phase cards and no sizes: the phases are its fluid’s (ps.fluid.GasWater flows gas and water), and the wells, connections and table rows the run allocates are counted. capacity=ps.Capacity(wells=50) (or m.set("capacity/wells", 50)) reserves more than the model holds; a loaded deck keeps the sizes it states (m.capacity).
The fluid’s class is its formulation. ps.fluid.GasWater and ps.fluid.OilWater are the two two-phase ones. ps.fluid.BlackOil is every other phase set: three phases, oil and gas, water only (BlackOil(water=...)), and CO2 or H2 storage, where the ps.storage package in physics=[...] supplies the PVT and the phases stay the deck’s own.
Tables are named columns, each a list or an array, one value per row.
ps.Summary asks for results by the names they are read by: field=[names], then wells=, groups= and regions= as {member: [names]}, "*" for every well or group.
State the rock compressibility. PetroScript never supplies one, and says so when it is missing.

3. Check it: read the rows by severity

for issue in merle.check():      # PetroScript's own checks: instant, no simulator, never raises
    print(issue.severity, issue.address, issue)

warning fluid/1/gas [physics.pvt.gas_range_below_initial] The gas PVT table ends at 350, below
  the initial datum pressure (3626); ...
warning initial/1 [physics.pvt.gas_heavier_than_water] At the datum pressure 3626 the gas weighs
  1838 against 1136.84 for the water; ...
info rock/reference_pressure [physics.rock.reference_pressure_remote] The rock's reference
  pressure lies 3376 below the datum pressure: ...
warning initial/1/datum_pressure [physics.pressure_gradient.other_unit] The datum pressure 3626
  is 13.4 times the model's own water column at the datum depth 2420 ...; read in psi it is
  0.923 times that column, nearer than the pressure as stated

Rosa typed the datum pressure in psi. She sets the fact at its address and reads the rows again:

merle = merle.set("initial/1/datum_pressure", 250.0)
rows = merle.check()
errors = [i for i in rows if i.severity == "error"]        # ps.compile refuses these
warned = [i for i in rows if i.severity == "warning"]      # facts to read, each with its address
assert not errors and not warned
print(rows)                 # ends: PetroScript's checks only; the simulator has not seen this model

Severity, never emptiness. An info row says what a check could not judge (the geometry checks stand down on a corner-point grid, and say so), so a sound model can hold rows.
A warning never blocks. ps.help(issue.code) explains one; a model that means it states acknowledged_warnings=[code].
A value may carry its own unit where a bare number would be misread:

in_psi = merle.set("initial/1/datum_pressure", ps.Quantity(3626.0, "psia"))

4. Compile, run and read the results

deck = ps.compile(merle)     # a Deck for "opm-flow", the default: in memory, nothing written or run
for row in deck.check():     # PetroScript's rows, then Flow's: a dry run and one short step
    print(row)               # last: [preflight.passed], or Flow's refusal in its own words
deck.write("merle")          # merle/MERLE.DATA

print(ps.simulators())       # what deck.run() finds: $PETROSCRIPT_FLOW, then OPM Flow on PATH
try:
    run = deck.run(folder="merle/run")     # deck.check() first, then Flow runs into the folder
except ps.CompileError as refused:         # Flow refused the deck before the run
    print(refused)
    raise
except ps.RunError as error:               # no simulator found, or the run stopped short
    print(error.exit_code, error.folder)
    raise

ps.compile prints each warning once, through ps.PetroScriptWarning, and raises ps.CompileError on an error row. Most decks Flow would refuse are refused while the model is built: SPE1 without its RSVD table raises ps.ModelError [physics.initial.depth_table_missing]. What only Flow finds, deck.run() finds first: it runs deck.check(), so a deck Flow refuses raises ps.CompileError before the run, [preflight.refused] and Flow’s last lines (deck.run(preflight=False) skips the pre-flight). A run that stops short later raises ps.RunError, with Flow’s lines and the run’s folder.

Results are read by PetroScript’s names, never Flow’s mnemonics:

import numpy as np

gas = run.field("gas_production_rate")       # a polars frame: date, value, unit
last = gas["value"][-1]                      # the value at the run's last step
bhp = run.wells("bottom_hole_pressure")      # date, well, value, unit; run.groups, run.regions too
try:
    run.field("FGPR")
except ps.PetroScriptError as refused:
    print(refused)           # [results.unknown] 'FGPR' is a Flow mnemonic: run.field("gas_...")

average = run.field("pressure")              # Flow's average, weighted by hydrocarbon pore volume
pressure = run.restart(at="2011-06-01")["pressure"]    # a cell array (nx, ny, nz) at a report time
pore_volume = run.static["pore_volume"]      # run.initial is the state at the start
mean = np.nansum(pressure * pore_volume) / pore_volume.sum()
i, j, k = np.unravel_index(np.nanargmin(pressure), pressure.shape)    # the lowest-pressure cell
again = ps.results.open("merle/run")         # the same Run, read back from its folder at any time

A cell the run did not simulate is NaN in a result array and 0 in pore_volume, so pressure.mean() would be NaN: use the summary’s average, or np.nansum over the pore volume, and np.nanargmin to find a cell.

5. Load a deck and edit it with verbs

loaded = ps.load("merle/MERLE.DATA")     # any Eclipse or Flow deck, with its includes

A deck the model cannot hold is refused whole with a ps.DeckImportError: the first gap, then N more gaps: and one line per gap, each led by its code. A loaded model is edited like any other:

m = loaded.until("2012-06-01", report="monthly")        # move the run's end
m = m.at("2011-09-01", ps.Shut(well="M-2"))             # a dated change
m = m.at("2011-12-15", ps.Open(well="M-2"), ps.Adjust(well="M-1", gas_rate=450_000.0))
m = m.set("rock/porosity", 0.2)                         # any fact by its address; None removes one
print(m.find("schedule/*/M-1/*", since="2011-06-01"))   # addresses, by pattern and date
for change in loaded.diff(m):                           # what changed: address, before, after
    print(change.address)
assert loaded.set({c.address: c.after for c in loaded.diff(m)}) == m

ps.Adjust changes one target or limit and keeps the rest of the control. A new ps.Produce replaces the whole control: a limit it does not restate is gone.
ps.Shut closes a well at the surface, ps.Stop leaves it open to the reservoir and ps.Open reopens it.
A date that is no report time becomes one. A date past the end is refused: m.until first.

Two edits come up again and again. One quantity from a date on, the later plan kept:

now = m.control("M-1", at="2012-03-01")                 # the control in force on a date
assert isinstance(now, ps.Produce)                      # it may be an Inject, a History or None
m = m.at("2012-03-01", now.replace(bhp_limit=60.0))     # typed: a misspelt field is refused
m = m.set({a: m.get(a).replace(bhp_limit=60.0)
           for a in m.find("schedule/*/M-1/control", since="2012-03-01")})

x.replace(...) is each statement’s own copy with changes. Python’s dataclasses.replace refuses a PetroScript statement, which is not a dataclass. An injector’s target is rate=: ps.Adjust(well="I-1", rate=...), or now.replace(rate=...) on the ps.Inject in force.

A well’s later plan superseded, then restated from a date:

m = m.set({a: None for a in m.find("schedule/*/M-2/*", since="2011-09-01")})
m = m.at("2011-09-01", ps.Produce(well="M-2", gas_rate=250_000.0, bhp_limit=50.0))

A new well mid-run is one m.set of the facts that change together: the well, its perforation and its control, each at its address:

m = m.set({
    "wells/M-3": ps.Well(name="M-3", i=5, j=5, introduced="2012-01-01"),
    "schedule/2012-01-01/M-3/perforation/1": ps.Perforate(well="M-3", k=(1, 2), diameter=0.2),
    "schedule/2012-01-01/M-3/control": ps.Produce(well="M-3", gas_rate=300_000.0, bhp_limit=50.0),
})

Set alone, a well there from the start is refused at once: [well.control.missing] Well 'M-4' has no reservoir control .... A well introduced later, such as M-3, is refused when the model is checked or compiled: [well.introduced.unanchored] Well 'M-3' is introduced 2012-01-01 and no perforation or control at that report time names it.

A day number is a date: the model’s start is a datetime.date, so day 365 is m.start + datetime.timedelta(days=365). A grid edit goes after the program the deck already states with m.apply; m.replace(edits=...) would replace that program, dropping every statement the deck made:

import datetime

m = m.at(m.start + datetime.timedelta(days=365), ps.Adjust(well="M-1", gas_rate=400_000.0))
m = m.apply(ps.grid.Scale(target="permeability_x", by=0.1, cells=ps.Box(i=(1, 3), j=(1, 3))))

6. The rest of the root, by example

history = ps.History(well="M-1", gas_rate=480_000.0, water_rate=2.0, control="gas_rate")
# observed rates, as a history match holds them; control= names the one the well is held to
injector = ps.Inject(well="I-1", phase="water", rate=500.0, bhp_limit=300.0)
observed = ps.InjectionHistory(well="I-1", phase="water", rate=480.0)
platform = ps.Group(name="PLATFORM")                     # ps.Well(..., group=platform)
steps = ps.Schedule(report=[ps.Every(31, 28, 31)])       # report steps in days: a loaded TSTEP card
fault = ps.Fault(name="F1", segments=[((5, 5), (1, 10), (1, 5), "plus_i")])  # Model(faults=[fault])
top_two = m.array("porosity")[ps.Box(k=(1, 2)).index]    # a cell array, and a 1-based inclusive box
layers = np.ones(m.grid.shape, dtype=int)
layers[:, :, 3:] = 2
m = m.set("regions", ps.Regions(fluid_in_place=layers))  # region maps: integers per cell
m = m.at("2012-04-01", ps.Numerics(next_step=1.0))       # numerics: Model(numerics=), or dated
m = m.set("output/static_file", True)                    # ps.Output: what a run writes, by address
rule = ps.actions.When(ps.actions.field().water_cut > 0.9, then=[ps.Shut(well="M-1")], name="WET")
m = m.at(m.start, rule)                                  # a rule the run evaluates, from day one
assert "ACTIONX" in ps.compile(m).files["MERLE.DATA"]   # ps.actions: rules written as ACTIONX

A lift table is a ps.lift table, and a control names it by the table itself. NeqSim can make one from the tubing and the fluid (how-to):

table = ps.lift.ProductionTable(
    number=1, bhp_reference_depth=2400.0, rate_kind="oil",
    rate=[100.0, 500.0, 1000.0], thp=[20.0, 40.0],
    bhp=np.array([[[[[120.0, 130.0, 150.0]]]], [[[[140.0, 150.0, 170.0]]]]]),
)                            # bhp's axes: thp, water fraction, gas fraction, artificial lift, rate
on_lift = ps.Produce(well="P1", oil_rate=800.0, bhp_limit=100.0, thp_limit=30.0, lift_table=table)

A check of your own is any object with a code that returns ps.Issue rows, warnings or info. Once registered it runs inside every m.check(), compile and deck.check():

class BhpFloor:
    code = "merle.bhp_floor"

    def __call__(self, model: ps.Model) -> list[ps.Issue]:
        controls = {a: model.get(a) for a in model.find("schedule/*/*/control")}
        return [
            ps.Issue(code=self.code, severity="warning", address=address,
                     message="the BHP limit is below the field's floor of 55 barsa",
                     fix="state bhp_limit=55.0 or more")
            for address, control in controls.items()
            if isinstance(control, ps.Produce) and (control.bhp_limit or 0.0) < 55.0
        ]


ps.checks.register(BhpFloor())
print([issue.address for issue in m.check() if issue.code == BhpFloor.code])

A model also describes and prints itself:

facts = m.describe()["facts"]            # JSON-safe: what the model states and what follows from it
print(facts["in_place"]["field"]["gas_in_place_upper_bound"]["value"])
m.to_python("merle_script")              # merle_script/MERLE.py rebuilds the model

ps.help() lists every class and namespace. The topics this page leaves out are a ps.help away: ps.aquifers, ps.groups (group targets), ps.actions (rules that fire during the run), ps.wells (the dated well statements beyond the root’s) and ps.bundle (a model from sheets).