One model, two simulators

examples/two_targets.py

"""One model, two simulators: compile SPE1 for OPM Flow and for JutulDarcy, run both, compare.

SPE1 (Odeh, 1981) is the industry's first comparative-solution benchmark: gas injected into a
three-phase black oil. This is its FIELD-unit case, SPE1CASE2. Both targets take it whole:
JutulDarcy's capability table refuses nothing in it, and JutulDarcy runs PetroScript's deck exactly
as it runs the deck's own source. One model compiles to two decks. Each run is read by the same
names, in the same units, so comparing them is a join on the date. The difference between the
engines is reported, never judged here: two solvers on one model differ a little, most where a
well switches control.

Flow keywords it exercises: the whole of SPE1CASE2 (PVTO, PVDG, SGOF, EQUIL, WCONPROD, WCONINJE),
plus FOPT in the summary.

Run it with ``python examples/two_targets.py [DECK]``. DECK defaults to SPE1CASE2 in the opm-tests
corpus at ``opm-tests/``; it writes ``out/two_targets/``. Flow is found through
``PETROSCRIPT_FLOW`` or ``flow`` on PATH, JutulDarcy through ``PETROSCRIPT_JUTULDARCY``.
"""

import sys

import petroscript as ps

path = sys.argv[1] if len(sys.argv) > 1 else "opm-tests/spe1/SPE1CASE2.DATA"
m = ps.load(path)
m = m.replace(summary=m.summary.replace(field=[*m.summary.field, "oil_production_total"]))

decks = {target: ps.compile(m, target) for target in ("opm-flow", "jutuldarcy")}
for target, deck in decks.items():
    deck.write(f"out/two_targets/{target}")
    print(target, "deck:", ", ".join(deck.files))

found = {simulator.name for simulator in ps.simulators()}
runs = {t: d.run(folder=f"out/two_targets/{t}/run") for t, d in decks.items() if t in found}
if not runs:
    print("no simulator found (set PETROSCRIPT_FLOW and PETROSCRIPT_JUTULDARCY); the decks are in "
          "out/two_targets")
    sys.exit(0)
for target in decks.keys() - runs.keys():
    print(f"{target}: no runner found, so there is nothing to compare its run with")


def series(run, quantity, well=None):
    """One quantity by date: the field's, or one well's."""

    frame = run.field(quantity) if well is None else run.wells(quantity).filter(well=well)
    return frame.select("date", "value", "unit")


asked = [
    ("oil rate", "oil_production_rate", None),
    ("cumulative oil", "oil_production_total", None),
    ("PROD BHP", "bottom_hole_pressure", "PROD"),
]
if len(runs) == 1:
    (target, run), = runs.items()
    for label, quantity, well in asked:
        last = series(run, quantity, well).row(-1)
        print(f"{target} {label} at the end: {last[1]:.1f} {last[2]}")
    sys.exit(0)

flow, jutul = runs["opm-flow"], runs["jutuldarcy"]
print(f"{'date':10}  {'what':14}  {'OPM Flow':>14}  {'JutulDarcy':>14}  {'diff':>7}")
for label, quantity, well in asked:
    # Flow writes every time step, JutulDarcy each report: compare on JutulDarcy's report dates
    both = series(jutul, quantity, well).join(series(flow, quantity, well), on="date",
                                               suffix="_flow")
    both = both.with_columns(diff=(both["value"] - both["value_flow"]) / both["value_flow"].abs())
    shown = both.select("date", "value", "unit", "value_flow", "diff").gather_every(24)
    for date, value, unit, flowed, diff in shown.iter_rows():
        print(f"{date:%Y-%m-%d}  {label:14}  {flowed:14.1f}  {value:14.1f}  {diff:7.2%}  {unit}")
    worst = both["diff"].abs().max()
    print(f"largest difference, {label}: {worst:.2%} over {both.height} report dates")