Skip to content

Commit

Permalink
deploy: bfee468
Browse files Browse the repository at this point in the history
  • Loading branch information
@prisae
prisae committed May 29, 2024
0 parents commit 24e013a
Show file tree
Hide file tree
Showing 111 changed files with 18,539 additions and 0 deletions.
4 changes: 4 additions & 0 deletions .buildinfo
Original file line number Diff line number Diff line change
@@ -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: e5475ec1680e119d49c72b93b3032173
tags: 645f666f9bcd5a90fca523b33c5a78b7
Empty file added .nojekyll
Empty file.
173 changes: 173 additions & 0 deletions _downloads/240167853b3ef667eee41137db914fec/basicreservoir.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,173 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib notebook"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# 2D Reservoir ESMDA example\n\n## Data Assimilation using ESMDA in Reservoir Simulation\n\n*Advanced Data Assimilation using Ensemble Smoother Multiple Data Assimilation\n(ESMDA) in Reservoir Simulation.*\n\n\\begin{align}m_j^a = m_j^f + C_{MD}^f (C_{DD}^f + \\alpha_i C_D)^{-1} (d_{uc,j} -\n d_j^f);\\quad \\text{for} \\quad j=1, 2, \\dots, N_e\\end{align}\n\n- Prior model: M (Ne, Nx, Ny)\n- Prior data: D (Ne, Nt)\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\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 = 25\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\n# dt = np.logspace(-5, -3, 10)\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# ox, oy = (1, 10, 20), (1, 20, 10)\n\n# Wells (if None, first and last cells are used with pressure 180 and 120)\n# wells = np.array([[15, 10, 180], [55, 25, 120], [30, 7, 140]])\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# TODO: change scale in imshow to represent meters\n\n# QC covariance, reference model, and first two random models\nfig, axs = plt.subplots(2, 2, constrained_layout=True)\naxs[0, 0].set_title('Model')\nim = axs[0, 0].imshow(perm_true.T, vmin=perm_min, vmax=perm_max)\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, vmin=perm_min, vmax=perm_max)\naxs[1, 1].set_title('Random Model 2')\naxs[1, 1].imshow(perm_prior[1, ...].T, vmin=perm_min, vmax=perm_max)\nfig.colorbar(im, ax=axs[1, :], orientation='horizontal',\n label='Log of Permeability (mD)')\nfor ax in axs[1, :]:\n ax.set_xlabel('x-direction')\nfor ax in axs[:, 0]:\n ax.set_ylabel('y-direction')\nfig.show()"
]
},
{
"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)\nax.set_title('Observed and prior data')\nax.plot(time, data_prior.T, color='.6', alpha=0.5)\nax.plot(time, data_true, 'ko', label='True data')\nax.plot(time, data_obs, 'C3o', label='Obs. data')\nax.legend()\nax.set_xlabel('Time (???)')\nax.set_ylabel('Pressure (???)')\nfig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ESMDA\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def restrict_permeability(x):\n x[:] = np.clip(x, perm_min, perm_max)\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, 3, figsize=(12, 5))\nax[0].set_title('Prior Mean')\nax[0].imshow(perm_prior.mean(axis=0).T)\nax[1].set_title('Post Mean')\nax[1].imshow(perm_post.mean(axis=0).T)\nax[2].set_title('\"Truth\"')\nax[2].imshow(perm_true.T)\nfig.show()"
]
},
{
"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)\nax.set_title('Prior and posterior data')\nax.plot(time, data_prior.T, color='.6', alpha=0.5)\nax.plot(time, data_post.T, color='C0', alpha=0.5)\nax.plot(time, data_true, 'ko')\nax.plot(time, data_obs, 'C3o')\nax.set_xlabel('Time (???)')\nax.set_ylabel('Pressure (???)')\nfig.show()"
]
},
{
"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
}
207 changes: 207 additions & 0 deletions _downloads/3a680ada310b2576fcfb5110c76a1304/basicESMDA.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,207 @@
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]}: $\\beta$ = {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.utils.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()
fig.show()


###############################################################################
# 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.utils.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()
Binary file added _downloads/46b4cb42d5bb56cc39e2b5b2b520b38d/gallery_python.zip
Binary file not shown.
Loading

0 comments on commit 24e013a

Please sign in to comment.