Read a finished run

A run leaves its files in a folder. deck.run() returns a Run that reads them, and ps.results.open(folder) reads the same folder at any time later. Everything is read by PetroScript’s names, the names a model asks for results by (gas_production_rate, never FGPR), and each value comes with its unit.

The run read here is the quickstart’s gas field, Merle (METRIC), run in OPM Flow:

import datetime

import numpy as np
import polars as pl

import petroscript as ps

merle = ps.Model(
    name="MERLE", units="metric", start="2009-06-01",
    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,
                 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],
                          gas_fvf=[0.0244, 0.0082, 0.0052, 0.0039],
                          gas_viscosity=[0.0135, 0.0155, 0.0185, 0.0215]),
        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,
    ),
    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=250.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", 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"]}),
)
run = ps.compile(merle).run(folder="merle/run")      # needs OPM Flow: docs/install.md
again = ps.results.open("merle/run")                 # the same run, read back from its folder

Time series: the summary

Each scope has its reader, and each returns a polars frame:

gas = run.field("gas_production_rate")        # date, value, unit
bhp = run.wells("bottom_hole_pressure")       # date, well, value, unit
everything = run.frame()                      # date, scope, name, quantity, value, unit
m1 = bhp.filter(pl.col("well") == "M-1")      # one well
print(gas["unit"][0], gas["value"][-1])       # sm3/day, and the last value

run.groups(name) and run.regions(name) read the other scopes the same way. Only what the model’s ps.Summary asked for is there: a name the run did not write is refused with the nearest it did ([results.unknown] ... nearest: gas_production_rate), and a Flow mnemonic with its name (run.field("FGPR") answers run.field("gas_production_rate")). ps.help("FGPR") does the same lookup.

A row is a time step, not a report date. Flow writes the summary at every step it takes, so the first rows of a monthly run sit at days 1, 4 and 13. run.times holds every step, in days since the start. To read a value on a date, filter to it:

on_new_year = gas.filter(pl.col("date").dt.date() == datetime.date(2010, 1, 1))
up_to = gas.filter(pl.col("date").dt.date() <= datetime.date(2010, 1, 1)).tail(1)   # or before
total = run.field("gas_production_total")["value"][-1]                              # at the end

The date a well reached a level. A producer’s bottom-hole pressure falls toward its limit and an injector’s rises, so test <= for a producer (>= for an injector), with a small tolerance: the simulator holds a pressure at its limit only to its solver’s accuracy. The first row that passes is the first step at or past the level; print it.

m1 = run.wells("bottom_hole_pressure").filter(pl.col("well") == "M-1")
reached = m1.filter(pl.col("value") <= 240.0 + 1e-3)      # barsa: M-1 is a producer
print(reached["date"][0].date().isoformat() if len(reached) else "not reached")

Test the pressure against its limit, not the rate against its target: a rate can read a hair under its target long before a limit takes over (SPE1’s producer reads 19,999.998 stb/day against its 20,000 target in September 2016, its BHP still at 3,509 psia, and reaches its 1,000 psia limit only in September 2017).

Cell arrays: the state in every cell

start = run.initial["pressure"]                            # (nx, ny, nz) at the start
later = run.restart(at="2011-06-01")["pressure"]           # at a restart report
dates = list(run.restarts)                                 # every restart report's date
final = run.restarts[dates[-1]]["pressure"]                # keyed by date: list(...)[-1] is the last
pore_volume = run.static["pore_volume"]                    # arrays that never change
average = np.nansum(later * pore_volume) / pore_volume.sum()   # pore-volume weighted
i, j, k = np.unravel_index(np.nanargmin(later), later.shape)   # the lowest-pressure cell

A cell the run did not simulate is NaN in a result array and 0 in pore_volume, so a plain later.mean() comes back NaN rather than quietly wrong: weight by the pore volume with np.nansum, or read the field average the summary already holds (run.field("pressure"), weighted by hydrocarbon pore volume). The arrays are read-only; copy one (later.copy()) to change it. The names are pressure, the saturations (water_saturation, gas_saturation, oil_saturation), solution_gor, and the static arrays list(run.static) lists (pore_volume, porosity, permeability_x, depth, the region maps, …). A static name asked of run.initial is answered with run.static["..."].

What else a run says

• run.units is the run’s unit system and run.phases the phases it simulated.
• run.grid is the grid the run wrote.
• run.log is what the simulator said while it ran: run.log.warnings, run.log.errors.
• A file the run did not keep is refused when a read needs it, naming the file: a folder with only the summary answers the summary and refuses the cell arrays.