Reservoir simulation, in Python
Describe the reservoir. PetroScript writes the deck.
PetroScript holds a reservoir model as one typed, frozen value, ps.Model. Build it by keyword or load an Eclipse or Flow deck, check it, compile it for OPM Flow or JutulDarcy, run it and read the results by name.
Drag a well
The whole path
Five calls from a word problem to results.
1. Build or load
ps.Model(...)·ps.load("CASE.DATA")One typed, frozen value.
ps.loadreads an Eclipse or Flow deck with its includes, or refuses it whole, naming every gap.2. Check
m.check()PetroScript's own checks: instant, no simulator, never raising. Every row has a code and an address; most say what to type.
3. Compile
ps.compile(m)A deck in memory for
"opm-flow"(the default) or"jutuldarcy".deck.write(folder)writes it.4. Pre-flight, run
deck.check()·deck.run(folder=...)The simulator's dry run and one short step, then the run.
5. Read
run.field("oil_production_rate")Results by PetroScript's names, as polars frames.
ps.results.open(folder)reads a run back later.
Python in, deck out
The model reads like the word problem.
A 10 × 10 × 3 waterflood: one injector, one producer, five years of monthly reports. Every part of the case is stated by keyword, in the model's units.
- One frozen value. Nothing is edited in place: every verb (
m.at,m.set,x.replace) returns a new value. - The control rule. In
ps.Produceandps.Injectthe one quantity given bare is the target and every other quantity is a*_limit. - Compile, then write.
ps.compile(m)makes the deck in memory;.write(folder)writes it.
import petroscript as ps
m = ps.Model(
name="WATERFLOOD", units="metric", start="2025-01-01",
grid=ps.grid.Cartesian(shape=(10, 10, 3), dx=50.0, dy=50.0,
dz=[4.0, 6.0, 10.0], top=2000.0),
rock=ps.Rock(porosity=0.2, permeability=100.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), ps.Well(name="PROD", i=10, j=10)],
schedule=ps.Schedule(end="2030-01-01", report="monthly", events=[
ps.Event("2025-01-01",
ps.Perforate(well="INJ", k=(1, 3), diameter=0.2),
ps.Perforate(well="PROD", k=(1, 3), diameter=0.2),
ps.Inject(well="INJ", phase="water", rate=350.0, bhp_limit=400.0),
ps.Produce(well="PROD", liquid_rate=350.0, bhp_limit=100.0))]),
summary=ps.Summary(
field=["oil_production_rate", "water_production_rate", "pressure"],
wells={"PROD": ["water_cut", "bottom_hole_pressure"]}),
)
for issue in m.check(): # instant, no simulator
print(issue.severity, issue.address, issue)
deck = ps.compile(m) # OPM Flow, the default target
deck.write("out/waterflood") # out/waterflood/WATERFLOOD.DATA
run = deck.run(folder="out/waterflood/run") # when OPM Flow is found
cut = run.wells("water_cut").filter(well="PROD") # date, well, value, unit-- Compiled with PetroScript 2.0.2 by Blauweiss AI
RUNSPEC
DIMENS
10 10 3 /
...
OIL
WATER
METRIC
START
1 'JAN' 2025 /
...
SCHEDULE
...
WELSPECS
'INJ' 'FIELD' 1 1 1* 'WAT' /
'PROD' 'FIELD' 10 10 1* 'OIL' /
/
COMPDAT
'INJ' 1 1 1 3 'OPEN' 1* 1* 0.2 /
'PROD' 10 10 1 3 'OPEN' 1* 1* 0.2 /
/
WCONINJE
'INJ' 'WATER' 'OPEN' 'RATE' 350 1* 400 /
/
WCONPROD
'PROD' 'OPEN' 'LRAT' 1* 1* 1* 350 1* 100 /
/
DATES
1 'FEB' 2025 /
/The phase cards come from the fluid; the sizes are counted from the model.
The same wells, in a grid
Every fact has an address.
The wells you placed above, now in a 26 × 18 × 6 grid with a fault: water runs ahead in the layers with the most pore space. It is an illustration, not a run. The model behind the page works the same way: each fact it states has a path, and the verbs read and write by that path, each returning a new model.
m.find("wells/*") # ['wells/INJ', 'wells/PROD']
m.get("rock/porosity") # 0.2
m.get("schedule/2025-01-01/PROD/control")
# ps.Produce(well='PROD', liquid_rate=350.0, bhp_limit=100.0)
later = m.at("2027-01-01",
ps.Adjust(well="PROD", liquid_rate=250.0))
later = later.set("rock/porosity", 0.18)
[c.address for c in m.diff(later)]
# ['rock/porosity', 'schedule/2027-01-01']
later.control("PROD", at="2028-01-01")
# ps.Produce(well='PROD', liquid_rate=250.0, bhp_limit=100.0)- Dated changes.
m.at(date, *changes)dates a change,m.set(address, value)sets any fact (Noneremoves one), andm.until(date)moves the run's end. - Adjust or replace.
ps.Adjustchanges one target and keeps the rest of the control; a newps.Producereplaces the whole control. - Loaded decks. A loaded deck is edited the same way:
ps.load("CASE.DATA")gives a model like any other, andm.difflists what changed: address, before, after.
Look it up, never guess
The library answers back.
A refusal names what you typed and what to type instead, at the step where the mistake is made, before anything runs. A model's refusal prints [code] address: message, then fix:.
Tests assert the code, never the prose: ps.testing.assert_issue(lambda: m.replace(rock=ps.Rock(porosity=0.2, permeability=100.0)), code="rock.permeability_z.unstated").
As you type a part: a misspelt field names the right one
ps.Produce(well="PROD", liquid_rate=350.0, bhp_limt=100.0)[field.unknown] Produce: no field 'bhp_limt'; did you mean 'bhp_limit'? Fields: well, oil_rate, ...A word it does not take says what it does take
ps.Inject(well="INJ", phase="gass", rate=350.0, bhp_limit=400.0)Inject.phase takes one of 'oil', 'water', 'gas', not 'gass'; fix: phase='gas'As you build: a value it will not guess
m.replace(rock=ps.Rock(porosity=0.2, permeability=100.0))[rock.permeability_z.unstated] rock.permeability_z is neither stated nor defined by the grid program; state the array or define it with a grid edit fix: state vertical_ratio= with permeability= (permeability_z over it; 1.0 for the same permeability along z), or permeability_z=
As you check: a number in the wrong unit, at its address
m.set("initial/1/datum_pressure", 2900.0).check() # psi, in a metric modelwarning initial/1/datum_pressure [physics.pressure_gradient.other_unit] The datum pressure 2900 is 12.9 times the model's own water column at the datum depth 2000 (224.616: atmospheric pressure plus the water gradient 0.1118 times the depth); read in psi it is 0.89 times that column, nearer than the pressure as statedOne of the three rows it prints. A warning never blocks;
ps.help(issue.code)explains one.As you compile: a card the target ignores is written, and named
ps.compile(m.replace(numerics=ps.Numerics(stack_depth=25)))[target.keyword_ignored] keyword/NSTACK: The opm-flow target reads and then ignores NSTACK; ...As you read: a Flow mnemonic points at the name
run.field("FGPR")[results.unknown] 'FGPR' is a Flow mnemonic: run.field("gas_production_rate")Any time: a name, either way
ps.help("WBHP")WBHP: the well quantity bottom_hole_pressure ps.Summary(wells={"*": ["bottom_hole_pressure"]}); run.wells("bottom_hole_pressure")
One model, more than one engine
Compile it for the simulator that fits.
Each reservoir simulator is a target behind ps.compile(m, target). A target refuses, at the address that states it, what its simulator cannot honour exactly, and every run is read back by the same names.
OPM Flow
defaulttested with OPM Flow 2026.04
- You call
ps.compile(m), the target"opm-flow".- Found by
PETROSCRIPT_FLOW, thenflowon PATH.- Before the run
deck.check()gives PetroScript's rows, then Flow's dry run of the whole deck and one short step. Flow's first report is then set against the model's datum pressure and contacts.- Reads back
run.field,run.wells,run.groupsandrun.regionsas polars frames (date, value, unit); cell arrays at a report time withrun.restart(at=...); andrun.static.- Refuses
- What Flow refuses, in Flow's own words, before the run; and an equation-of-state fluid, by name.
JutulDarcy
tested with JutulDarcy 0.4.0
- You call
ps.compile(m, "jutuldarcy"): the same model, black-oil, as JutulDarcy reads it.- Found by
PETROSCRIPT_JUTULDARCY: a program that runs Julia with JutulDarcy and GeoEnergyIO, or a relay into a container that holds them.- Before the run
deck.check()asks JutulDarcy to parse and set up the deck.- Reads back
- The field's and the wells' rates, totals, ratios, pressure and volumes, by the same names and units as a Flow run. Cell arrays, the grid and the log are not read yet, and each says so.
- Refuses
[target.unsupported]for a keyword it cannot honour exactly,[target.item_dropped]for an item it would read past (a THP limit, a group control), and an equation-of-state fluid. A group or region summary quantity is a warning.
NeqSim
tested with NeqSim 3.25.0
- Its role
- Not a reservoir target: it makes the lift tables a Flow model's producers flow through.
- You call
ps.lift.generate(tubing, model=m, well=..., fluid=..., composition=..., correlation=..., number=..., rate=..., rate_kind=...).- Found by
PETROSCRIPT_NEQSIM, orsimulator=for one call.- Gives back
- An ordinary
ps.lift.ProductionTablewhosesourcenames NeqSim, its version, the correlation, the composition, the datum and where the THP axis came from. A control uses it withlift_table=. - Refuses
- A point NeqSim's driver refuses is named by its axis values, and the table is refused whole. An unknown or mixed-case word is refused before anything runs.
Today: models compile and run on OPM Flow and JutulDarcy, and NeqSim makes a producer's lift table. Compositional EOS fluids are modelled (ps.fluid.Compositional with ps.pvt.Components): a compositional deck loads and writes back; neither target runs one yet.
In development: injector lift tables, PVT tables and flowlines with NeqSim; compositional and thermal models on JutulDarcy; and GEOS.
m = ps.load("opm-tests/spe1/SPE1CASE2.DATA") # SPE1, field units
decks = {t: ps.compile(m, t) for t in ("opm-flow", "jutuldarcy")}
runs = {t: d.run(folder=f"out/two_targets/{t}/run") for t, d in decks.items()}
oil = {t: run.field("oil_production_rate") for t, run in runs.items()}
both = oil["jutuldarcy"].join(oil["opm-flow"], on="date", suffix="_flow")The example runs both and reports the difference by date, never judging it: two solvers on one model differ a little, most where a well switches control.
tubing = ps.lift.Tubing(
measured_depth=[0.0, 2000.0], true_vertical_depth=[0.0, 2000.0],
inner_diameter=[0.124], roughness=[1.5e-5],
wellhead_temperature=40.0, temperature_gradient=0.03)
table = ps.lift.generate(
tubing, model=m, well="PROD", fluid=oil, composition="RES",
correlation="beggs_brill", number=1,
rate=[50.0, 200.0, 400.0, 800.0], rate_kind="liquid",
thp=[10.0, 20.0, 40.0, 80.0], water_fraction=[0.0, 0.5], gas_fraction=[0.0])
print(table.source) # NeqSim's version, the correlation, the composition, ...Start here
From install to your own deck.
Install
Python 3.10 or newer; numpy, polars and pydantic come with it.
A simulator is a separate program. PetroScript finds OPM Flow through PETROSCRIPT_FLOW, then flow on PATH. Without one, everything up to the deck works; deck.run raises a ps.RunError that says none was found.
Learn
The quickstart: one gas field built by keyword, checked, compiled, run and read, then loaded back and edited by address.
The examples: small scripts, each run end to end.
Bring a deck
m = ps.load("CASE.DATA")
m.to_python("case")
$ petroscript load CASE.DATA case
$ petroscript check case/CASE.pyA deck the model cannot hold is refused whole with a ps.DeckImportError: the first gap, then one line per gap.