Grid and rock

Reservoir meaning

The grid divides the reservoir volume into simulation cells. Rock properties then describe how much fluid each cell can store and how easily phases can move through it. The current PetroScript implementation supports simple Cartesian grids and the corner-point grid subset used by the OPM SPE9/SPE10 black-oil benchmarks. Rock properties are attached directly to full-grid cells as NumPy-shaped tensors.

PetroScript API

Cartesian

ps.grid.Cartesian(...) describes a corner-aligned Cartesian grid. shape is (nx, ny, nz). dx, dy, and dz are each a finite scalar, an explicit Quantity scalar, one value per column, row or layer along their own axis (a list), or a NumPy (nx, ny, nz) tensor. top is the depth/elevation of the top surface for each (i, j) column and accepts a finite scalar, explicit Quantity, or NumPy (nx, ny) tensor. Optional active masks use a scalar, one value per layer, or an explicit (nx, ny, nz) tensor containing only 0/1 or boolean active-cell flags.

import numpy as np
import petroscript as ps

shape = (10, 10, 3)
dz = np.broadcast_to(np.array([20.0, 30.0, 50.0]), shape).copy()
grid = ps.grid.Cartesian(shape=shape, dx=1000.0, dy=1000.0, dz=dz, top=8325.0)

When cell widths are constant along each axis but differ by column/layer, state one width per column, row or layer as a list instead of expanding a full 3D tensor:

grid = ps.grid.Cartesian(
    shape=(5, 4, 3),
    dx=[100.0, 100.0, 120.0, 120.0, 150.0],
    dy=[200.0, 200.0, 200.0, 200.0],
    dz=[15.0, 25.0, 45.0],
    top=1000.0,
)

Optional grid-policy fields stay semantic:

coordinate_unit="metres" (or "feet", "centimetres") records the coordinate-length unit for generated grid geometry.
map_axes=ps.grid.MapAxes(y_axis=..., origin=..., x_axis=...) records export map-origin metadata for generated grid files.
minimum_pore_volume=ps.grid.MinPoreVolume(threshold=...) records a minimum pore-volume activation policy.
pinch_policy=ps.grid.Pinch(minimum_thickness=..., gap_mode="gap", transmissibility_mode="all_pinched_cells", multiplier_scope="all_pinched_cells") records a bounded pinch-out connectivity policy.

These objects are policies on the grid, not user-facing keyword rows.

Static depth and pore-volume corrections, and pore-volume scaling, are statements in the model’s edits=, each naming the array it sets or scales and the ps.Box of cells it reaches:

edits = (
    ps.grid.Assign(target="depth", value=2511.0, cells=ps.Box(i=(1, 2), j=(1, 2), k=(1, 1))),
    ps.grid.Assign(target="pore_volume", value=0.0, cells=ps.Box(k=(1, 1))),
    ps.grid.Scale(target="pore_volume", by=1.05, cells=ps.Box(k=(1, 1))),
)
# model = model.apply(*edits)

model.apply(*edits) puts the statements after the ones the model already holds, so a loaded deck’s own EQUALS, MULTIPLY and the rest stay in place and the new statements act on the arrays they leave. model.replace(edits=...) replaces the whole program instead: on a loaded deck it drops every statement the deck stated, so use it only to state a program from scratch.

The edits are a program the model holds as stated, in the order stated, and the render writes it back card for card (EQUALS, MULTIPLY, COPY, OPERATE, ADDREG and the rest, one card per statement, after the array cards of the statement’s section). A reader that wants the field after the program asks the model for it: model.array("permeability_z") is the cell array after the program, what the simulator will see.

The physics checks read the arrays after the program, so an authored edit reaches the same plausibility gate an imported deck’s edits do. The render never evaluates. Arrays the simulator derives (DEPTH, PORV, TRAN*), the face multiplier policies and the grid’s cell-size vectors have no cell field on the model, so an edit on them changes nothing a reader of the model can see.

Reservoir region maps are domain concepts under ps.Regions, the model’s regions=. PetroScript supports transmissibility-region maps and operation-region maps without exposing the backend row names:

multnum_array = np.ones(shape, dtype=np.int64)
opernum_array = np.ones(shape, dtype=np.int64)
regions = ps.Regions(transmissibility=multnum_array, operation=opernum_array)
# model.replace(regions=regions)

Corner-point grids

Corner-point grids describe pillars and corner depths rather than a rectilinear cell size. This is the common shape for field-scale legacy datasets because cell tops, bottoms, and pillars can vary independently. A ps.grid.CornerPoint holds pillars (six numbers per pillar: the x, y and z of its top point, then of its bottom point) and corner_depths (eight per cell), with active marking the active cells. ps.io.read_grdecl(path) reads a GRDECL file’s geometry into one and ps.io.write_grdecl(grid, path=...) writes it back; ps.grid.CornerPoint.from_box restates a Cartesian grid as the corner-point grid it describes.

PetroScript currently supports the regular OPM benchmark subset needed by SPE9 and SPE10 MOD01: grid dimensions, flat float64 pillar coordinates, flat float64 corner depths, and optional integer (nx, ny, nz) active-cell masks containing only 0/1 flags. A corner-point grid is the model’s grid= as a Cartesian one is:

box = ps.grid.CornerPoint.from_box(
    ps.grid.Cartesian(shape=(2, 2, 1), dx=100.0, dy=100.0, dz=20.0, top=1000.0)
)
geometry = ps.grid.CornerPoint(
    shape=(2, 2, 1),
    pillars=box.pillars,
    corner_depths=box.corner_depths,
    active=np.ones((2, 2, 1), dtype=np.int64),
)
# model = model.replace(grid=geometry)

A grid between two surfaces

Hand-computing pillars and corner_depths for a structure is where corner-point authoring goes wrong: the pillar bookkeeping is easy to transpose, and a transposed pillar produces inverted cells that no amount of fluid description will rescue. ps.grid.CornerPoint.from_surface computes them for a grid between two surfaces. It returns an ordinary CornerPoint, and the compiler emits the same geometry keywords; it adds no deck keyword and no grid type of its own.

ps.grid.CornerPoint.from_surface takes a top surface of node depths indexed [j][i] — that is, (ny + 1) rows of (nx + 1) values, one per pillar — and infers nx and ny from its shape. Pass either a uniform thickness or an explicit bottom surface of the same shape; a varying thickness is expressed as a base surface.

from_surface stands vertical pillars on a regular dx by dy lattice and checks that node depths increase strictly downwards along each pillar before the grid is constructed, so it cannot return an inverted or zero-volume cell. A bottom surface that crosses above top is refused.

Rock

ps.Rock(...) stores pore-volume and directional flow-capacity inputs. Porosity is a fraction. The permeability is stated as permeability with a vertical_ratio, or as the directional permeability_x, permeability_y, and permeability_z. Each property is a finite scalar, explicit Quantity, one value per layer, or explicit 3D NumPy tensor; other 1D/2D arrays fail at object construction.

rock = ps.Rock(porosity=[0.20, 0.25, 0.30], permeability=150.0, vertical_ratio=0.1,
               compressibility=4.0e-6, reference_pressure=3000.0)

Rock compressibility

ps.Rock(compressibility=..., reference_pressure=...) describes how pore volume changes around a pressure anchor, the two fields stated together. Both values are finite scalar semantic values (float or explicit Quantity), or one value each per PVT region as two lists of one length (per rock region with table_regions="rock_regions"). It is optional and belongs with rock properties rather than fluid PVT.

RockCompaction

ps.grid.RockCompaction describes water-induced pore-volume change for one or more rock regions. Use it when the rock response depends on pressure and water saturation instead of a single compressibility coefficient.

The API is intentionally region-oriented, one table per rock region in each field:

loading (ROCK2D, the pore-volume multipliers) and unloading (ROCK2DTR, the transmissibility multipliers; neither is a loading or an unloading branch): rows of a pressure, then one multiplier per water-saturation node, (pressure, *multipliers);
water_saturation_nodes: the region’s water saturations;
overburden: rows of (depth, overburden_pressure);
model: the compaction hysteresis model name as a semantic selection.

Attach the resulting ps.grid.RockCompaction to ps.Rock(compaction=...) and use ps.Regions(rock=...) when cells should select different compaction families. The current compiler lowers this to the implemented OPM Flow rock-compaction keyword family.

Regions

Region maps store integer region IDs over the grid. Implemented fields are saturation, the directional relperm_x, relperm_y and relperm_z, pvt, equilibration, rock, and Regions.flux. A region map can be a scalar positive integer or an explicit positive-integer NumPy tensor that later matches the model grid. Region IDs are one-based and select the corresponding family/table defined by the model; strings, booleans, zero, fractional IDs, negative IDs, and quantity wrappers are rejected locally. Flux regions are used by bounded transmissibility-region policies rather than by stringly keyword cards.

Faults and transmissibility policies

ps.Fault(name=..., segments=[...]) names a fault, each segment a row (i, j, k, face): one-based inclusive (first, last) ranges of cells and the cell face the fault runs along ("plus_i", "minus_i", "plus_j" and so on). The model’s faults= holds them.

fault = ps.Fault(name="F1", segments=[((1, 1), (1, 4), (1, 3), "plus_i")])
fault_multiplier = ps.grid.FaultMultiplier(fault="F1", factor=0.1)
# model.replace(faults=[fault, fault_multiplier])

ps.grid.FaultMultiplier(fault=..., factor=...) applies a transmissibility multiplier to a named fault. The structure validator checks that the fault exists and that segment ranges fit the grid.

ps.grid.RegionPairMultiplier(from_region=..., to_region=..., factor=..., directions=...) applies a factor between two positive flux-region IDs. Use ps.Regions(flux=...) to assign flux regions over cells. Directions are the words "x", "y" and "z".

ps.grid.RegionPairMultiplier takes from_region=, to_region=, factor=, the region map= it reads, directions=, and connections= ("all", "neighbor_only", or "non_neighbor_only") selecting which connection subset the factor applies to. Diffusivity and thermal-conductivity multipliers between regions (MULTREGD, MULTREGH) have no PetroScript statement: OPM Flow stops on both cards.

A cell is a one-based (i, j, k) tuple. ps.grid.NonNeighbor(first=, second=, transmissibility=) states a connection between two cells, and ps.grid.ConnectionTransmissibilityEdit(first=, second=, mode=, value=) multiplies or overrides a connection’s transmissibility. ps.grid.NonNeighborConnectionPolicy(ignore_simulator_connections=True) in the model’s non_neighbor_connections= disables simulator-generated non-neighbor links. These names are domain concepts; the compiler owns any target keyword spelling.

Pinch and minimum pore-volume policy

ps.grid.MinPoreVolume(threshold=...) stores the minimum pore-volume cutoff as a positive scalar semantic policy on the grid. ps.grid.Pinch stores the minimum thickness, its gap_mode ("gap" or "no_gap"), an optional transmissibility_mode and an optional multiplier_scope for the bounded pinch-out behavior OPM Flow runs.

Radial/spider geometry

Radial and Spider describe radial/spider geometry from semantic dimensions: inner_radius, radial_widths, angular_widths, layer thickness, top depth, and optional active-cell intent: ps.grid.Radial(...) and ps.grid.Spider(...). Users do not spell target keyword rows; the compiler owns RADIAL, SPIDER, INRAD, DRV, DTHETAV, PERMR, and PERMTHT lowering.

Units and shapes

Under FIELD, grid lengths and depths are deck-native feet, permeability is md, and porosity is unitless. Shapes follow NumPy conventions:

scalar: one finite value applied everywhere;
(nx, ny, nz): one finite value per full-grid cell for dx, dy, dz, porosity, permeability_x, permeability_y, and permeability_z;
(nx, ny): one finite value per surface column for top.
Cartesian active: scalar or (nx, ny, nz) 0/1 or boolean active-cell flags.
corner-point coord: flat float64 pillar coordinate values for (nx + 1) * (ny + 1) pillars;
corner-point zcorn: flat float64 eight-corner depths per full-grid cell.
corner-point actnum: (nx, ny, nz) integer active-cell tensor with 0/1 flags.
region maps: scalar or (nx, ny, nz) positive integer arrays.
coordinate_unit: one of the words "metres", "feet" and "centimetres".
minimum pore-volume and pinch thresholds: finite positive scalar values.
fault segments: one-based inclusive (start, end) ranges over i, j, and k, with the face a word such as "plus_i".
connection cells: one-based positive (i, j, k) handles that must fit within the model grid dimensions.

A layered property may be stated as one value per layer (porosity=[0.20, 0.25, 0.30]), which the model spreads over each layer’s cells, or as an explicit broadcasted tensor where that improves readability, for example np.broadcast_to(np.array([500.0, 50.0, 200.0]), shape). For a model rebuilt from a deck’s arrays (SPE9’s permeability fields, for example), compute directional permeability_x, permeability_y, and permeability_z tensors from the reservoir intent and pass them to Rock; do not replay materialized fixture arrays as an authoring shortcut.

Validation and common mistakes

Model constructors validate local numeric type, rank, dtype, and finiteness. Cross-object checks still validate full grid compatibility before rendering. Common mistakes are swapping (nx, ny, nz) order, providing a vector whose length is no axis of the grid, using percentages instead of porosity fractions, passing string region IDs, or accidentally assigning the vertical permeability equal to horizontal permeability when anisotropy matters.

Lowers to

Current status: emitted now for implemented Cartesian, corner-point, and radial/spider geometry paths.

DIMENS records (nx, ny, nz) in RUNSPEC; radial/spider grids also emit RADIAL or SPIDER.
DX, DY, DZ, and TOPS describe Cartesian geometry in GRID.
DXV, DYV, and DZV describe Cartesian axis-vector geometry in GRID when a width is stated as one value per column, row or layer (a list).
INRAD, DRV, DTHETAV, DZ, and TOPS describe radial/spider geometry in GRID.
SPECGRID, COORD, ZCORN, and ACTNUM describe implemented corner-point geometry and activity masks.
MAPAXES, GRIDUNIT, MINPV, MINPVV, and PINCH describe map-export metadata, coordinate-unit, minimum pore-volume, and bounded pinch policy when supplied. MINPV and MINPVV are two statements and both survive: the manual’s own MINPVV entry (6.3.124) says the array “is different to the MINPV and MINPORV keywords … that set a constant minimum threshold pore volume for all cells”, and the simulator applies each. A deck stating the scalar and then assigning the array over boxes — opm-tests/minpv/MINPV300 — keeps both cards; dropping the scalar leaves every cell outside those boxes on MINPVV’s own 1.0e-6 default and activates cells the deck’s run deactivates. The array goes through the box-edit workspace like any other stated array, so EQUALS MINPVV over a box composes with it.
The active-cell mask carries both directions. EQUALS ACTNUM states one flag per cell — 6.3.1 spells it “integers equal to either 0 or 1” — and the mask takes the flag the deck wrote — 1 for active as readily as 0 for inactive — applying the boxes in deck order so a later card wins where two overlap. A whole-field ACTNUM 1 opening the mask and boxed ACTNUM 0 cards carving out of it is the corpus’s ordinary shape; reading the opening card as the closing one is not a narrower import, it is the opposite grid.
MULTPV, PORV, and DEPTH describe pore-volume scaling, pore-volume override, and cell-depth override: a loaded deck’s grid-stated multiplier, and the program’s ps.grid.Scale and ps.grid.Assign statements on pore_volume and depth.
MULTNUM and OPERNUM describe implemented transmissibility and operation region maps from ps.Regions.
PORO, PERMX, PERMY, and PERMZ describe Cartesian/corner-point rock properties in GRID; radial/spider geometry lowers horizontal rock properties as PERMR and PERMTHT plus PERMZ. The model keeps one kind per axis: on a radial grid permeability_x is the radial permeability and permeability_y the angular one, and every card that names them, the base arrays and the program’s COPY, EQUALS, MULTIPLY and OPERATE rows alike, spells them PERMR and PERMTHT.
FLUXNUM, FAULTS, MULTFLT, and MULTREGT describe implemented flux regions, semantic faults, fault multipliers, and region transmissibility multipliers. A MULTFLT is restated where the deck’s own clock put it. Stated in GRID or EDIT it is a property of the model and is restated in that section; stated in SCHEDULE it is an EVENT, and it is restated in SCHEDULE at the position it was stated, after the report steps that preceded it. The distinction is not cosmetic: a grid edit standing before the first report step is byte-identical to the same card in GRID, which is why the GRID and EDIT populations cannot tell the difference — but once a report step precedes the card, applying it from time zero opens the fault early, which moves transmissibility, the non-neighbor-connection set the simulator derives from it, and the time axis that follows. The faults themselves stay declared in GRID: a scheduled card names a fault, it does not redeclare one, so FAULTS must still stand before any MULTFLT that refers to it. VAPPARS is restated the same way, for the same reason: stated in SOLUTION it names the model’s starting condition and is restated there, and stated in SCHEDULE it is an event restated in SCHEDULE. Read into the initial state instead, a schedule-stated card would be hoisted to time zero: on opm-tests/norne/NORNE_ATW2013 that moves the first diverging restart record from the seventh to the first after time zero, on 32,555 PRESSURE cells.
MULTX, MULTX-, MULTY, MULTY-, MULTZ, and MULTZ- describe directional face transmissibility multipliers, with BOX/ENDBOX around any that name less than the whole grid. An imported multiplier is restated in the section the deck stated it in — GRID multipliers in GRID, EDIT multipliers in EDIT. The two sections are separate statements that the simulator composes by multiplication; they are not two spellings of one number, and folding them into a single composite discarded an attribution that matters wherever an EDIT-section TRANX/TRANY/TRANZ mounts an explicit transmissibility the two stand on opposite sides of. A model that states no section — every authored model — still leaves the placement to the compiler, which puts a whole-grid card in GRID and a boxed one in EDIT.

The composing reading is now measured rather than argued. opm-tests/mult/MULTXYZ11 states all six MULT* arrays in both sections, and the target’s own INIT writes each resolved multiplier as the PRODUCT of the two — 0.5 × 5.0 = 2.5, 2.1435 × 0.43451435 = 0.931381509, six for six to float32 storage precision, with neither stage’s value alone appearing anywhere.

A statement made in both stages: what the target does, per family

GRID and EDIT are two statements rather than two spellings of one, and what the simulator does with the pair is not the same for every family. Each rule below is pinned to a controlled run of the pinned target rather than to a reading, because the three families disagree and an analogy from one to another states numbers the run does not use.

Face multipliers COMPOSE. Product of the two stages, measured above.
Region transmissibility multipliers REPLACE — the later stage wins. MULTXYZ11 also states MULTREGT on the pair 2→1 twice, GRID 0.45678 and EDIT 1.09876. Every face crossing the pair carries 1.09876 relative to a comparable face that does not, in all three directions to eight significant figures. The product, 0.50189, appears nowhere. The controlled arm separates replacement from a card going unread: with the EDIT card stripped and the deck re-run, the same faces carry 0.45678. So each card is restated in its own section and the simulator does the replacing — merging the factors would state a transmissibility the run does not use, by a factor of 2.19.
A transform on an array the simulator READS is fidelity-only. dispersion/DISPERC_INJECTION scales PORO with an EDIT-section OPERATE; in the target’s INIT the porosity is the scaled value and PORV / (PORO × DX × DY × DZ × NTG) is 1.00000000 in all 9,000 cells, on both sides of the edited box. The pore volume is derived from the post-edit porosity, so the card computes the same numbers wherever it stands and its stage decides only which section it is restated in.

The other direction — a GRID-section card on an EDIT-stage array — keeps its refusal. Nothing has measured it, and moving a statement about a derived quantity earlier is not the same statement.

TRANX, TRANY, and TRANZ describe explicit EDIT-stage transmissibility across an axis’s positive faces, which REPLACES what the simulator computes from permeability and geometry rather than scaling it (7.3.69-7.3.71). A uniform override and a varying field are two statements. One value over the whole grid is a scalar the edit algebra composes with, and it renders as an EDIT-section EQUALS. A field that varies cell by cell says something no scalar can — on opm-tests/mult/MULT2D_TRANS a column of sealed faces inside an otherwise uniform sheet — and is carried as a cell field in the model’s transmissibility=, restated as its own TRAN* array card with the deck’s own cell count. Collapsing the field to any one of its values states a different grid.
MULTREGT region ids may be left unstated, and the absence is carried as an absence. Manual 6.3.138 gives items 1 and 2 no default value, so a 1* there is the card declining to name a region — coupling every one of them — rather than shorthand for a number. Resolving it into the pairs a region map holds would state couplings the deck never wrote, and would state them against whatever the map holds at PetroScript’s reading rather than the simulator’s.
ROCK describes rock compressibility in PROPS when supplied.
ROCKCOMP, ROCKOPTS, ROCKTAB, ROCK2D, ROCK2DTR, ROCKWNOD, OVERBURD, and ROCKNUM describe implemented rock-compaction semantics when supplied (see Rock compaction).
INIT is emitted for grid output when the grid asks for it. The compiler writes no NOECHO or ECHO: Flow 2026.04 does not support either and warns on both.

radial/spider geometry is implemented through ps.grid.Radial and ps.grid.Spider. LGR geometry, local grid refinement, and unstructured grids remain future/unsupported. Non-neighbor connection editing and bounded region-set materialization are implemented through the semantic APIs documented above.