Use regions and faults

Real reservoirs differ from place to place: rock types with their own relative permeability, compartments with their own PVT or initial pressure, faults that slow flow across them. In PetroScript each of these is a region map (an integer per cell, regions numbered from 1) plus the tables or statements keyed by region number.

What differs	The map (ps.Regions(...))	What is keyed by region
Relative permeability	saturation=	relperm={1: (...), 2: (...)}
PVT	pvt=	fluid={1: ..., 2: ...}
Initial state	equilibration=	initial={1: ..., 2: ...}
Reporting (in-place, region rates)	fluid_in_place=	ps.Summary(regions={1: [...]})
Grid edits	operation=, transmissibility=, flux=	cells=ps.grid.InRegion("operation", 2)

The example: a 4 x 4 x 2 oil field (METRIC) whose bottom layer is a second rock type.

import numpy as np

import petroscript as ps

shape = (4, 4, 2)
rock_type = np.ones(shape, dtype=int)
rock_type[:, :, 1] = 2                                   # layer 2 is rock type 2

top_rock = ps.relperm.OilWater(sw=[0.2, 0.8], krw=[0.0, 1.0], krow=[1.0, 0.0],
                               pcow=[0.0, 0.0])
tight_rock = ps.relperm.OilWater(sw=[0.3, 0.7], krw=[0.0, 0.5], krow=[1.0, 0.0],
                                 pcow=[0.0, 0.0])

m = ps.Model(
    name="REGIONS", units="metric", start="2025-01-01",
    grid=ps.grid.Cartesian(shape=shape, dx=100.0, dy=100.0, dz=10.0, top=2000.0),
    rock=ps.Rock(porosity=[0.22, 0.15], permeability=[200.0, 20.0],   # one value per layer
                 vertical_ratio=0.1, compressibility=5.0e-5, reference_pressure=200.0),
    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,
    ),
    regions=ps.Regions(saturation=rock_type, fluid_in_place=rock_type),
    relperm={1: (top_rock,), 2: (tight_rock,)},          # each region's tables
    initial=ps.initial.Equilibrium(datum_depth=2000.0, datum_pressure=200.0,
                                   oil_water_contact=2030.0),
    wells=[ps.Well(name="P", i=2, j=2)],
    schedule=ps.Schedule(end="2025-06-01", report="monthly", events=[
        ps.Event("2025-01-01", ps.Perforate(well="P", k=(1, 2), diameter=0.2),
                 ps.Produce(well="P", oil_rate=200.0, bhp_limit=50.0))]),
    summary=ps.Summary(field=["oil_production_rate"], regions={1: ["oil_in_place"],
                                                               2: ["oil_in_place"]}),
)
assert "SATNUM" in ps.compile(m).files["REGIONS.DATA"]

relperm=[top_rock, tight_rock], a flat list of two oil-water tables, is refused: relperm states water_oil more than once: state each once. A flat list is one region’s tables (an oil-water and a gas-oil table, say); several regions are a mapping by region number, each region’s tables a tuple (a list runs too).

The same holds for the fluid (fluid={1: ..., 2: ...} with ps.Regions(pvt=...)) and the initial state (initial={1: ..., 2: ...} with ps.Regions(equilibration=...)).

Faults

A fault is a named set of cell faces, stated on the model’s faults=. Each segment is (i, j, k, face): three 1-based, inclusive (first, last) ranges and the face the fault runs along ("plus_i", "minus_i", "plus_j", …). Its transmissibility multiplier sits beside it from the start, or in a ps.Event from a date:

fault = ps.Fault(name="F1", segments=[((2, 2), (1, 4), (1, 2), "plus_i")])  # east of i=2
m = m.replace(faults=[fault, ps.grid.FaultMultiplier(fault="F1", factor=0.01)])
m = m.at("2025-04-01", ps.grid.FaultMultiplier(fault="F1", factor=0.5))    # x0.5 from then
deck = ps.compile(m).files["REGIONS.DATA"]
assert "FAULTS" in deck and "MULTFLT" in deck

Multipliers stated on successive dates multiply, as Flow reads them: 0.5 then 0.5 is 0.25.

Edits by region

A grid edit takes a region of an operation, transmissibility or flux map as its cells:

m = m.set("regions/operation", rock_type)
m = m.apply(ps.grid.Scale(target="permeability_x", by=0.5,
                          cells=ps.grid.InRegion("operation", 2)))
layer_2 = m.array("permeability_x")[ps.Box(k=(2, 2)).index]
assert float(layer_2.max()) == 10.0                                   # 20 mD halved

m.array(name) reads a cell array after the grid edits. A transmissibility multiplier (transmissibility_multiplier_z, Flow’s MULTZ) loads, edits and writes back, but m.array cannot read it yet ([array.unknown]): read the statements that set it with m.get("edits").

ps.help(ps.Regions) lists every map, and the grid and rock reference every grid edit.