Chapter 3 — The numbers that will not fit in the module

An authoring module is source a person reads. A million-cell porosity array, an observed rate history, a lift table of seventy-five thousand pressures: none of those belong in source. They sit in files beside the module, and the module reads them with the tools that read numbers, numpy and polars, into the fields that take them.

from pathlib import Path

import numpy as np
import polars as pl

here = Path(__file__).parent                    # the module's own folder, not the cwd
model = model.set("rock/porosity", np.load(here / "poro.npy"))   # one value per cell
porosity = model.rock.porosity
assert isinstance(porosity, np.ndarray) and not porosity.flags.writeable  # its own read-only copy

A cell array is an array. Every field that holds one value per cell takes a numpy array of the grid’s shape, indexed [i - 1, j - 1, k - 1], or a single value for every cell; a few, such as the grid’s dz and the initial water saturation, also take one value per layer. The model holds its own read-only copy, so rewriting the file or the array later changes no built model.

A table is named columns. A PVT or saturation table, a lift table’s axes, a depth table: each column is a field taking a list or an array, and from_frame reads them by name from what polars reads:

oil = ps.pvt.DeadOil.from_frame(pl.read_csv(here / "pvdo.csv"))   # pressure, oil_fvf, oil_viscosity
model = model.set("fluid/1/oil", oil)

A history is statements. Observed rates are ps.History statements, one per well and date, each in the step of its date; a date inside the run that is no report time becomes one:

history = pl.read_csv(here / "history.csv", try_parse_dates=True)   # date, well, oil_rate
for row in history.iter_rows(named=True):
    model = model.at(row["date"], ps.History(well=row["well"], oil_rate=row["oil_rate"]))
assert len(model.find("schedule/*/PROD1/control")) == 1 + history.height

A body of numbers is an array too. A lift table’s bottom-hole pressures are one array whose axes are the tubing-head pressure, the water fraction, the gas fraction, the artificial lift and the rate:

table = ps.lift.ProductionTable(
    number=1, bhp_reference_depth=2500.0, rate_kind="oil",
    rate=[100.0, 500.0], thp=[20.0, 40.0],
    bhp=np.load(here / "lift.npy"),             # shape (2, 1, 1, 1, 2): thp, ..., rate
)

Geometry has its own file. A corner-point grid comes from a GRDECL file: ps.io.read_grdecl(path) reads its SPECGRID, COORD, ZCORN and, when stated, ACTNUM and MAPAXES into a ps.grid.CornerPoint, and refuses any other keyword in the file by name: such a file is numbers, never a deck. ps.io.write_grdecl writes the file it reads back:

assert isinstance(model.grid, ps.grid.Cartesian)   # write_grdecl takes a box or corner-point grid
ps.io.write_grdecl(model.grid, path=here / "grid.grdecl")
corner_point = ps.io.read_grdecl(here / "grid.grdecl")
assert corner_point.shape == model.grid.shape

Every failure happens at the line that reads. A missing file, a missing column, an array of the wrong shape: numpy and polars raise at the load, and the model refuses a wrong shape or a missing column at the field that takes it, naming the field.

The model writes its own arrays back beside itself. m.to_python(folder) writes the model’s module, HELLO.py here, and beside it, in HELLO_arrays/, each array of more than 64 values as a .npy the module loads from its own folder. A module that resolves its data through a folder parameter (def build(data=Path(__file__).parent): ...) runs against another realization when called with another folder: that is a convention, deliberately not machinery.