Simulator: MiniRes, a pure-Python two-phase forward model

docs/tutorials/README.md

@@ -5,6 +5,7 @@ Here are some tutorials.
 ## Running PIPT and POPT
 
 - [`tutorial_pipt.ipynb`](pipt/TinyBox/tutorial_pipt): Tutorial for running PIPT
+- [`tutorial_minires.ipynb`](pipt/MiniRes/tutorial_minires): The same, with a forward model that needs nothing installed -- a MiniRes waterflood
 - [`tutorial_popt.ipynb`](popt/5Spot/tutorial_popt): Tutorial for running POPT
 
 ## Data structures

docs/tutorials/pipt/MiniRes/CONFIG_ESMDA.toml

@@ -0,0 +1,70 @@
+[ensemble]
+    ne    = 100
+    state = "permx"   # log-permeability, one value per grid cell
+
+    [ensemble.prior_permx]
+        vario = "sph"
+        mean  = "priormean.npz"
+        var   = 1.0
+        range = 5.0
+        aniso = 1.0
+        angle = 0.0
+        grid  = [16, 16, 1]
+
+[dataassim]
+    savefolder = "Results"
+    scheme     = "esmda"
+    analysis   = "approx"
+    energy     = 99.0
+    obsname    = "steps"
+    data       = "data.csv"
+    datavar    = "var.csv"
+    savedata   = ["ensemble_misfit"]
+
+    # ESMDA settings
+    [dataassim.mda]
+        tot_assim_steps = 3
+        inflation_param = [3, 3, 3]
+
+[simulator]
+    # The report points are *step indices*: MiniRes takes a uniform dt.
+    dt          = 0.025
+    reporttype  = "steps"
+    reportpoint = [2, 4, 6, 8, 10, 12]
+    parallel    = 1
+    datatype    = [
+        "WWCT:PRD1", "WWCT:PRD2", "WWCT:PRD3", "WWCT:PRD4",
+        "WOPR:PRD1", "WOPR:PRD2", "WOPR:PRD3", "WOPR:PRD4",
+    ]
+
+    # Passed to minires.ResSim as it stands -- less the permeability,
+    # which is what the ensemble state supplies, one field per member.
+    [simulator.model]
+        Nx  = 16
+        Ny  = 16
+        por = 0.2
+
+        [[simulator.model.wells]]
+            name = "INJ1"
+            xy   = [0.5, 0.5]
+            rate = 1.0
+
+        [[simulator.model.wells]]
+            name = "PRD1"
+            xy   = [0.05, 0.05]
+            rate = -0.25
+
+        [[simulator.model.wells]]
+            name = "PRD2"
+            xy   = [0.95, 0.05]
+            rate = -0.25
+
+        [[simulator.model.wells]]
+            name = "PRD3"
+            xy   = [0.05, 0.95]
+            rate = -0.25
+
+        [[simulator.model.wells]]
+            name = "PRD4"
+            xy   = [0.95, 0.95]
+            rate = -0.25

docs/tutorials/pipt/MiniRes/data.csv

@@ -0,0 +1,7 @@
+steps,WWCT:PRD1,WWCT:PRD2,WWCT:PRD3,WWCT:PRD4,WOPR:PRD1,WOPR:PRD2,WOPR:PRD3,WOPR:PRD4
+2,0.015235853987721568,-0.05199920531202478,0.03752255979032287,0.0470282358195607,0.15244824056730816,0.1848910246568841,0.25639202015836426,0.23418787038282088
+4,-0.0006616316175194097,-0.04265219637867901,0.04397906690243332,0.37986283178125513,0.25325692831363705,0.30636206034840163,0.2733731750727793,0.12179206810338616
+6,0.7121617362686763,-0.04794413004144995,0.7997165437846246,0.8058189238661145,0.06732583255659913,0.21595352277980293,0.12217855975388625,0.0401947210427031
+8,0.7975451683203961,0.10294571650886156,0.9101860547384159,0.916591579588462,0.06589624072266154,0.2414029516420759,0.13618973117833594,0.005099405088204681
+10,0.8651076650607917,0.6923317350907747,0.956140474899063,0.9865003950070657,0.021622676738664103,0.024737083276081456,-0.022559436727134008,0.05001669429847978
+12,0.9665251856623154,0.8433864077844226,0.9132882369878853,0.9576057423100545,0.023493671181500596,0.05687725624415992,0.0569305083093295,0.024680358399554493

docs/tutorials/pipt/MiniRes/make_case.py

@@ -0,0 +1,51 @@
+"""Make this tutorial's synthetic case: a truth, and observations of it.
+
+A twin experiment: the "true" permeability is itself a draw from the prior that
+``CONFIG_ESMDA.toml`` specifies, MiniRes is run on it, and the well reports are
+perturbed by the observation noise. Nothing external is needed -- which is the
+point of this tutorial's simulator.
+
+    python make_case.py
+"""
+
+from copy import deepcopy
+
+import numpy as np
+import pandas as pd
+
+from input_output import read_config
+from misc.sampling import random_stream
+from misc.structures.layout import StateLayout
+from pipt.misc_tools.extract_tools import extract_prior_info
+from simulator.minires import MiniRes
+
+SIGMA = 0.05  # observation noise (absolute, on both the water cut and the oil rate)
+
+kwda, kwsim, kwens = read_config.read("CONFIG_ESMDA.toml")
+nx, ny, _ = kwens["prior_permx"]["grid"]
+
+# The prior's mean: a flat log-permeability field. Also fixes the state's length.
+np.savez("priormean.npz", permx=np.zeros(nx * ny))
+
+# The truth: one draw from that prior
+truth, _ = StateLayout.from_prior_info(
+    extract_prior_info(deepcopy(kwens)), 1, rng=random_stream(7), save=False
+)
+truth = truth[:, 0]
+np.savez("truth.npz", permx=truth)
+
+# The observations: MiniRes on the truth, plus noise
+records = MiniRes(kwsim).run_fwd_sim({"permx": truth}, 0)
+report = kwsim["reportpoint"]
+
+data = pd.DataFrame([{k: float(v[0]) for k, v in row.items()} for row in records], index=report)
+data += SIGMA * np.random.default_rng(42).standard_normal(data.shape)
+data.index.name = kwsim["reporttype"]
+data.to_csv("data.csv")
+
+var = pd.DataFrame({col: [f"['abs', {SIGMA**2}]"] * len(report) for col in data.columns}, index=report)
+var.index.name = data.index.name
+var.to_csv("var.csv")
+
+print(f"Wrote priormean.npz, truth.npz, data.csv, var.csv ({data.shape[0]} report points"
+      f" x {data.shape[1]} data types)")