Build the grid and the rock

The grid is the model’s grid=: its shape (nx, ny, nz), the cells’ sizes and the depth of the top. The rock is rock=ps.Rock(...): porosity, permeability and compressibility. A model can start uniform and gain detail without changing shape, but each field takes its own shapes:

Field	One number	A list	An array
dx	every cell	one per column: nx values	(nx, ny, nz)
dy	every cell	one per row: ny values	(nx, ny, nz)
dz	every cell	one per layer: nz values	(nx, ny, nz)
top	one depth	no	(nx, ny): one depth per column
active, net_to_gross	every cell	one per layer	(nx, ny, nz)
porosity, the permeabilities	every cell	one per layer	(nx, ny, nz)
compressibility, reference_pressure	the whole rock	one per PVT region	no

A dx list is per column even when nx equals nz, so a list meant per layer silently varies the widths across the grid instead. Cell indices are 1-based and inclusive, as in a deck: ps.Box(i=(1, 3), k=(2, 2)) is columns 1 to 3 of layer 2.

import numpy as np

import petroscript as ps

shape = (10, 10, 3)
top = 2000.0 + 5.0 * np.arange(10)[:, None] + np.zeros((10, 10))   # m: dips 5 m per column
active = np.ones(shape, dtype=int)
active[9, 9, :] = 0                                               # a dead corner (i=10, j=10)

grid = ps.grid.Cartesian(
    shape=shape,
    dx=100.0, dy=100.0,                    # one value: every cell
    dz=[10.0, 15.0, 20.0],                 # one value per layer
    top=top,                               # one depth per column, (nx, ny)
    active=active,                         # 1 active, 0 not
)
permeability_x = np.full(shape, 150.0)
permeability_x[:5, :, :] = 300.0           # the western half is better rock
rock = ps.Rock(
    porosity=[0.22, 0.18, 0.15],           # per layer
    permeability_x=permeability_x,         # mD, a full array
    permeability_y=permeability_x,
    permeability_z=0.1 * permeability_x,
    compressibility=5.0e-5, reference_pressure=200.0,     # 1/bar at bar
)

permeability= with vertical_ratio= is the short form: one horizontal permeability for x and y, and z as a share of it. The two forms mix across directions, but one direction stated in both is refused (rock.permeability.two_forms). Leave the compressibility out and OPM Flow runs incompressible rock, which m.check() warns about.

Change part of a grid

A grid edit changes an array over a box (or a region), after what the model already states. It is applied with m.apply(...), and m.array(name) reads the result:

m = ps.Model(
    name="GRID", units="metric", start="2025-01-01", grid=grid, rock=rock,
    fluid=ps.fluid.OilWater(
        oil=ps.pvt.DeadOil(pressure=[50.0, 300.0], oil_fvf=[1.1, 1.08], oil_viscosity=[2.0, 2.2]),
        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.8], krw=[0.0, 1.0], krow=[1.0, 0.0],
                                 pcow=[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="P", i=3, j=3)],
    schedule=ps.Schedule(end="2025-03-01", report="monthly", events=[ps.Event(
        "2025-01-01", ps.Perforate(well="P", k=(1, 3), diameter=0.2),
        ps.Produce(well="P", oil_rate=200.0, bhp_limit=50.0))]),
)
m = m.apply(ps.grid.Scale(target="permeability_x", by=0.1,         # a tight streak
                          cells=ps.Box(i=(1, 10), j=(4, 4), k=(2, 2))))
m = m.apply(ps.grid.Assign(target="porosity", value=0.05, cells=ps.Box(k=(3, 3))))
streak = m.array("permeability_x")[ps.Box(i=(1, 10), j=(4, 4), k=(2, 2)).index]
assert float(streak.max()) == 30.0                                  # 300 x 0.1

The edits: ps.grid.Assign (set), Scale (multiply), Add, Clip (hold between bounds), Copy (one array into another) and Operate. m.apply keeps the edits a loaded deck already states and adds yours after them; m.replace(edits=[...]) would replace them all.

Other grids

• A dipping or faulted reservoir from a surface: ps.grid.CornerPoint.from_surface(top, dx=, dy=, nz=, thickness=) builds conformable layers under a top surface.
• A deck’s corner-point grid loads as it is: ps.load("CASE.DATA").grid.
• A well-centred radial grid: ps.grid.Radial(...).

ps.help(ps.grid.Cartesian), ps.help(ps.Rock) and the grid and rock reference have every field and edit.