The technical reference for engineers and reviewers.

What this engine is

haql-engine is a from-scratch, correctness-first three-phase black-oil simulator written in Python and NumPy. It is built to make a reservoir result inspectable: every run writes a machine-readable record, a provenance manifest and a material-balance ledger.

Quick start

Run the smallest useful checks first. These commands assume the engine directory is the current directory.

cd engine

# 1. Fast internal checks: Jacobian, conservation, determinism.
python3 selftest.py

# 2. One short SPE1 run.
python3 run_spe1.py --T 30 --quiet --json /tmp/spe1_30.json

# 3. A small ensemble. Use --workers only when you have CPU/RAM to spare.
python3 run_uncertainty.py \
  --input uncertainty_input.example.json \
  --out /tmp/haql-ensemble \
  --T 30 --workers 2

# 4. Aggregate the ensemble into P10/P50/P90 and CSV exports.
python3 run_forecast.py \
  --input /tmp/haql-ensemble/forecast_input.json \
  --json /tmp/haql-forecast.json \
  --csv-dir /tmp/haql-forecast-csv

If a command fails, read the last log under /tmp/haql-ensemble/logs/. A failed realization is recorded in experiment.json; it is not silently removed.

Architecture

The code is deliberately split into a numerical core, a case driver and two surrounding layers: geological realizations and forecast aggregation.

deckparse.py ──► fieldgrid.py ──► blackoil_vs.py ──► run_spe1.py
     │                │                  │                 │
     │                │                  │                 ├── provenance.py
     │                │                  │                 └── JSON / material balance
     │                │                  │
     │                │                  └── PVT + relperm + wells + Newton
     │                └── Cartesian grid, depth, transmissibility
     └── SPE1 deck reader

geology.py ──► run_uncertainty.py ──► run_spe1.py (one process per realization)
                     │
                     └── experiment.json + realization_summary.csv + forecast_input.json
                                                   │
                                                   ▼
                                          forecast.py / run_forecast.py

Design rules

• No hidden state. Inputs, solver settings and source identity are hashed into a run manifest.
• No silent failures. A failed or incomplete realization remains in the experiment record.
• No fake precision. Every published number is tied to a fixture, a gate or a named limitation.
• No claim without a test. The regression suites are part of the product, not an afterthought.

How a run works

1. Parse the deck. deckparse.py tokenises the SPE1-style input. The current driver is intentionally explicit about its supported subset.
2. Build the grid. fieldgrid.py creates the Cartesian cells, layer depths, volumes, transmissibilities and well indices.
3. Initialise the state. Pressure, solution GOR and saturations are set from the datum pressure, contacts and optional scalar overrides.
4. Assemble the residual. blackoil_vs.py writes the mass-balance equations for oil, water and gas, including well source terms.
5. Solve Newton. The Jacobian is assembled and the linear system is solved until the residual gate is met or the timestep is cut.
6. Switch variables when needed. Each cell carries either free-gas saturation or dissolved gas, depending on its phase state.
7. Record the result. The driver writes report-step records, material-balance drift, a terminal status and a provenance block.

The full numerical detail is in engine/README.md and engine/VALIDATION.md; this page is the map, not a substitute for those documents.

Uncertainty workflow

geology.py turns one base grid into reproducible realizations. A realization is data, not a hidden random seed: it is written to JSON and hashed before the solver sees it.

The generated ensemble is driven by run_uncertainty.py. It writes one run per realization, records failures, and produces a forecast_input.json for the aggregation layer.

Fault barriers

Each fault in structure.faults may carry a transmissibility_multiplier in [0, 1]. The default is 1.0, which preserves the pre-barrier behaviour. Aliases are transmissibility, transmissibility_factor and barrier_multiplier; specifying more than one alias is rejected.

"structure": {
  "enabled": true,
  "sigma_ft": 0.0,
  "faults": [
    {"axis": "x", "position_ft": 3000.0, "throw_ft": 12.0,
     "side": "positive", "transmissibility_multiplier": 0.0}
  ]
}

"scalars": {
  "fault_transmissibility_multiplier": {
    "distribution": "uniform", "min": 0.0, "max": 1.0
  }
}

A fault normal to x or y acts on the cross-fault faces between the last cell centre below position_ft and the first cell centre at or above it. This is the same boundary the structural throw uses. side selects the displaced block, not the barrier face. Multiple faults on one face multiply. The effective factor is the per-fault baseline times the sampled global scalar. Set structure.enabled to true; a disabled structure emits no fault barrier.

The realization summary records summary.fault_barrier. The applied grid exposes grid.face_multiplier with one value per grid.faces entry, and a run artifact records fault_summary. The run configuration also stores the fault block and scalar override for provenance.

Full public explanation: Fault-flow barriers. Source schema: engine/GEOLOGICAL_UNCERTAINTY.md.

Forecast layer

forecast.py is analysis-only. It reads completed runs, aligns them on a common time horizon, and reports per-case statistics. It never calls the solver and never mutates a run.

• Quantile or exceedance P10/P50/P90 conventions, selectable in the input.
• Common-horizon interpolation with explicit no-extrapolation behavior.
• Mean, median, sample standard deviation, min and max for each metric.
• Failed-run reporting and parent experiment links, not silent exclusion.
• JSON and CSV outputs for downstream review.

If the ensemble is small, the quantiles are descriptive, not a certified reserve estimate. If the horizon differs between runs, the forecast truncates and says so.

Parallel ensembles new

run_uncertainty.py --workers N now runs independent realizations through separate subprocesses. The solver itself remains single-threaded; the parallelism is at the ensemble level.

python3 run_uncertainty.py \
  --input uncertainty_input.example.json \
  --out /tmp/haql-ensemble \
  --T 365 --workers 4

# deterministic experiment order is always real-0000, real-0001, ...
# failed runs are still written and linked, not removed.

Read the dedicated page for determinism, resource guidance and failure behavior: Parallel ensemble execution →

Evidence, not adjectives

Full detail: the public evidence page and engine/VALIDATION.md.

Limits you must state before showing it

• Not a full-field simulator. Cartesian SPE1-class grids; no corner-point/EGRID, no SPE9, no real field grids.
• Not a fault-seal calculator. Fault-face transmissibility is a user-supplied multiplier. There is no juxtaposition, shale-gouge, capillary-entry-pressure or geomechanical seal model, and faults are vertical planes on a regular x/y grid.
• Not fast. It is slower than OPM Flow; the value is inspection and cross-checking, not throughput.
• Not history-matched. The uncertainty ensemble is not conditioned to observed production.
• Not a reserves certification. P10/P50/P90 are descriptive outputs of the supplied model and distributions.
• Not a commercial support product. There is no SLA, no enterprise security certification and no client field audit yet.

File map

Glossary

Feature log

Each shipped feature gets a page here so the documentation and the code move together.