deploy: 55ff177

.buildinfo

@@ -0,0 +1,4 @@
+# Sphinx build info version 1
+# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
+config: 27f628a3846decce8a9f1557573596e8
+tags: 645f666f9bcd5a90fca523b33c5a78b7

_downloads/240167853b3ef667eee41137db914fec/basicreservoir.ipynb

@@ -0,0 +1,173 @@
+{
+  "cells": [
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "%matplotlib widget"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n# 2D Reservoir ESMDA example\n\nEnsemble Smoother Multiple Data Assimilation (ES-MDA) in Reservoir Simulation.\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import numpy as np\nimport matplotlib.pyplot as plt\n\nimport resmda\n\n# For reproducibility, we instantiate a random number generator with a fixed\n# seed. For production, remove the seed!\nrng = np.random.default_rng(1848)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Model parameters\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Grid extension\nnx = 30\nny = 25\nnc = nx*ny\n\n# Permeabilities\nperm_mean = 3.0\nperm_min = 0.5\nperm_max = 5.0\n\n# ESMDA parameters\nne = 100                  # Number of ensembles\ndt = np.zeros(10)+0.0001  # Time steps (could be irregular, e.g., increasing!)\ntime = np.r_[0, np.cumsum(dt)]\nnt = time.size\n\n# Assumed sandard deviation of our data\ndstd = 0.5\n\n# Observation location indices (should be well locations)\nox, oy = 1, 1\n\n# Wells (if None, first and last cells are used with pressure 180 and 120)\nwells = None"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Create permeability maps for ESMDA\n\nWe will create a set of permeability maps that will serve as our initial\nguess (prior). These maps are generated using a Gaussian random field and are\nconstrained by certain statistical properties.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Get the model and ne prior models\nRP = resmda.RandomPermeability(nx, ny, perm_mean, perm_min, perm_max)\nperm_true = RP(1, random=rng)\nperm_prior = RP(ne, random=rng)\n\n\n# QC covariance, reference model, and first two random models\npinp1 = {\"origin\": \"lower\", \"vmin\": perm_min, \"vmax\": perm_max}\nfig, axs = plt.subplots(2, 2, figsize=(6, 6), constrained_layout=True)\naxs[0, 0].set_title(\"Model\")\nim = axs[0, 0].imshow(perm_true.T, **pinp1)\naxs[0, 1].set_title(\"Lower Covariance Matrix\")\nim2 = axs[0, 1].imshow(RP.cov, cmap=\"plasma\")\naxs[1, 0].set_title(\"Random Model 1\")\naxs[1, 0].imshow(perm_prior[0, ...].T, **pinp1)\naxs[1, 1].set_title(\"Random Model 2\")\naxs[1, 1].imshow(perm_prior[1, ...].T, **pinp1)\nfig.colorbar(im, ax=axs[1, :], orientation=\"horizontal\",\n             label=\"Log of Permeability (mD)\")\nfor ax in axs[1, :].ravel():\n    ax.set_xlabel(\"x-direction\")\nfor ax in axs[:, 0].ravel():\n    ax.set_ylabel(\"y-direction\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Run the prior models and the reference case\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Instantiate reservoir simulator\nRS = resmda.Simulator(nx, ny, wells=wells)\n\n\ndef sim(x):\n    \"\"\"Custom fct to use exp(x), and specific dt & location.\"\"\"\n    return RS(np.exp(x), dt=dt, data=(ox, oy))\n\n\n# Simulate data for the prior and true fields\ndata_prior = sim(perm_prior)\ndata_true = sim(perm_true)\ndata_obs = rng.normal(data_true, dstd)\n\n# QC data and priors\nfig, ax = plt.subplots(1, 1, constrained_layout=True)\nax.set_title(\"Observed and prior data\")\nax.plot(time*24*60*60, data_prior.T, color=\".6\", alpha=0.5)\nax.plot(time*24*60*60, data_true, \"ko\", label=\"True data\")\nax.plot(time*24*60*60, data_obs, \"C3o\", label=\"Obs. data\")\nax.legend()\nax.set_xlabel(\"Time (s)\")\nax.set_ylabel(\"Pressure (bar)\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## ESMDA\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "def restrict_permeability(x):\n    \"\"\"Restrict possible permeabilities.\"\"\"\n    np.clip(x, perm_min, perm_max, out=x)\n\n\nperm_post, data_post = resmda.esmda(\n    model_prior=perm_prior,\n    forward=sim,\n    data_obs=data_obs,\n    sigma=dstd,\n    alphas=4,\n    data_prior=data_prior,\n    callback_post=restrict_permeability,\n    random=rng,\n)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Posterior Analysis\n\nAfter running ESMDA, it's crucial to analyze the posterior ensemble of\nmodels. Here, we visualize the first three realizations from both the prior\nand posterior ensembles to see how the models have been updated.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Plot posterior\nfig, ax = plt.subplots(1, 2, figsize=(8, 5), constrained_layout=True)\npinp2 = {\"origin\": \"lower\", \"vmin\": 2.5, \"vmax\": 3.5}\nax[0].set_title(\"Prior Mean\")\nim = ax[0].imshow(perm_prior.mean(axis=0).T, **pinp2)\nax[1].set_title(\"Post Mean\")\nax[1].imshow(perm_post.mean(axis=0).T, **pinp2)\nfig.colorbar(im, ax=ax, label=\"Log of Permeability (mD)\",\n             orientation=\"horizontal\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Observing the monitored pressure at cell (1,1) for all realizations and the\nreference case, we can see that the ensemble of models after the assimilation\nsteps (in blue) is closer to the reference case (in red) than the prior\nensemble (in gray). This indicates that the ESMDA method is effectively\nupdating the models to better represent the observed data.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Compare posterior to prior and observed data\nfig, ax = plt.subplots(1, 1, constrained_layout=True)\nax.set_title(\"Prior and posterior data\")\nax.plot(time*24*60*60, data_prior.T, color=\".6\", alpha=0.5)\nax.plot(time*24*60*60, data_post.T, color=\"C0\", alpha=0.5)\nax.plot(time*24*60*60, data_true, \"ko\")\nax.plot(time*24*60*60, data_obs, \"C3o\")\nax.set_xlabel(\"Time (s)\")\nax.set_ylabel(\"Pressure (bar)\")"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "resmda.Report()"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.11.9"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}

_downloads/3a680ada310b2576fcfb5110c76a1304/basicESMDA.py

@@ -0,0 +1,206 @@
+r"""
+Linear and non-linear ES-MDA examples
+=====================================
+
+A basic example of ES-MDA using a simple 1D equation.
+
+Geir Evensen gave a talk on *Properties of Iterative Ensemble Smoothers and
+Strategies for Conditioning on Production Data* at the IPAM in May 2017.
+
+Here we reproduce the examples he showed on pages 34 and 38. The material can
+be found at:
+
+- PDF: http://helper.ipam.ucla.edu/publications/oilws3/oilws3_14079.pdf
+- Video can be found here:
+  https://www.ipam.ucla.edu/programs/workshops/workshop-iii-data-assimilation-uncertainty-reduction-and-optimization-for-subsurface-flow/?tab=schedule
+
+Geir gives the ES-MDA equations as
+
+.. math::
+    x_{j,i+1} &= x_{j,i} + (C^e_{xy})_i \left((C^e_{yy})_i +
+                 \alpha_iC^e_{dd}\right)^{-1} \left(d + \sqrt{\alpha_i}
+                 \varepsilon_j - g(x_{j,i})\right) \ , \\
+    y_{j,i+1} &= g(x_{j,i+1}) \ .
+
+The model used for this example is
+
+.. math::
+    y = x(1+\beta x^2) \ ,
+
+which is a linear model if :math:`\beta=0`.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import resmda
+
+# sphinx_gallery_thumbnail_number = 3
+
+###############################################################################
+# Forward model
+# -------------
+
+
+def forward(x, beta):
+    """Simple model: y = x (1 + β x²) (linear if beta=0)."""
+    return np.atleast_1d(x * (1 + beta * x**2))
+
+
+fig, axs = plt.subplots(
+        1, 2, figsize=(8, 3), sharex=True, constrained_layout=True)
+fig.suptitle("Forward Model:  y = x (1 + β x²)")
+px = np.linspace(-5, 5, 301)
+for i, b in enumerate([0.0, 0.2]):
+    axs[i].set_title(
+            f"{['Linear model', 'Nonlinear model'][i]}: β = {b}")
+    axs[i].plot(px, forward(px, b))
+    axs[i].set_xlabel('x')
+    axs[i].set_ylabel('y')
+
+
+###############################################################################
+# Plotting functions
+# ------------------
+
+def pseudopdf(data, bins=200, density=True, **kwargs):
+    """Return the pdf from a simple bin count.
+
+    If the data contains a lot of samples, this should be "smooth" enough - and
+    much faster than estimating the pdf using, e.g.,
+    `scipy.stats.gaussian_kde`.
+    """
+    x, y = np.histogram(data, bins=bins, density=density, **kwargs)
+    return (y[:-1]+y[1:])/2, x
+
+
+def plot_result(mpost, dpost, dobs, title, ylim):
+    """Wrapper to use the same plotting for the linear and non-linear case."""
+
+    fig, (ax1, ax2) = plt.subplots(
+            1, 2, figsize=(10, 4), sharey=True, constrained_layout=True)
+    fig.suptitle(title)
+
+    # Plot Likelihood
+    ax2.plot(
+        *pseudopdf(resmda.rng.normal(dobs, size=(ne, dobs.size))),
+        'C2', lw=2, label='Datum'
+    )
+
+    # Plot steps
+    na = mpost.shape[0]-1
+    for i in range(na+1):
+        params = {
+            'color': 'C0' if i == na else 'C3',    # Last blue, rest red
+            'lw': 2 if i in [0, na] else 1,        # First/last thick
+            'alpha': 1 if i in [0, na] else i/na,  # start faint
+            'label': ['Initial', *((na-2)*('',)), 'MDA steps', 'MDA'][i],
+        }
+        ax1.plot(*pseudopdf(mpost[i, :, 0], range=(-3, 5)), **params)
+        ax2.plot(*pseudopdf(dpost[i, :, 0], range=(-5, 8)), **params)
+
+    # Axis and labels
+    ax1.set_title('Model Parameter Domain')
+    ax1.set_xlabel('x')
+    ax1.set_ylim(ylim)
+    ax1.set_xlim([-3, 5])
+    ax1.legend()
+    ax2.set_title('Data Domain')
+    ax2.set_xlabel('y')
+    ax2.set_xlim([-5, 8])
+    ax2.legend()
+
+
+###############################################################################
+# Linear case
+# -----------
+#
+# Prior model parameters and ES-MDA parameters
+# ''''''''''''''''''''''''''''''''''''''''''''
+#
+# In reality, the prior would be :math:`j` models provided by, e.g., the
+# geologists. Here we create $j$ realizations using a normal distribution of a
+# defined mean and standard deviation.
+
+# Point of our "observation"
+xlocation = -1.0
+
+# Ensemble size
+ne = int(1e7)
+
+# Data standard deviation: ones (for this scenario)
+obs_std = 1.0
+
+# Prior: Let's start with ones as our prior guess
+mprior = resmda.rng.normal(loc=1.0, scale=obs_std, size=(ne, 1))
+
+
+###############################################################################
+# Run ES-MDA and plot
+# '''''''''''''''''''
+
+def lin_fwd(x):
+    """Linear forward model."""
+    return forward(x, beta=0.0)
+
+
+# Sample an "observation".
+l_dobs = lin_fwd(xlocation)
+
+lm_post, ld_post = resmda.esmda(
+    model_prior=mprior,
+    forward=lin_fwd,
+    data_obs=l_dobs,
+    sigma=obs_std,
+    alphas=10,
+    return_steps=True,  # To get intermediate results
+)
+
+
+###############################################################################
+
+plot_result(lm_post, ld_post, l_dobs, title='Linear Case', ylim=[0, 0.6])
+
+
+###############################################################################
+# Original figure from Geir's presentation
+# ''''''''''''''''''''''''''''''''''''''''
+#
+# .. image:: ../_static/figures/Geir-IrisTalk-2017-34.png
+
+
+###############################################################################
+# Nonlinear case
+# --------------
+
+def nonlin_fwd(x):
+    """Nonlinear forward model."""
+    return forward(x, beta=0.2)
+
+
+# Sample a nonlinear observation; the rest of the parameters remains the same.
+n_dobs = nonlin_fwd(xlocation)
+nm_post, nd_post = resmda.esmda(
+    model_prior=mprior,
+    forward=nonlin_fwd,
+    data_obs=n_dobs,
+    sigma=obs_std,
+    alphas=10,
+    return_steps=True,
+)
+
+###############################################################################
+
+plot_result(nm_post, nd_post, n_dobs, title='Nonlinear Case', ylim=[0, 0.7])
+
+
+###############################################################################
+# Original figure from Geir's presentation
+# ''''''''''''''''''''''''''''''''''''''''
+#
+# .. image:: ../_static/figures/Geir-IrisTalk-2017-38.png
+
+
+###############################################################################
+
+resmda.Report()