Skip to content

Commit

Permalink
Fluvial example (#11)
Browse files Browse the repository at this point in the history
Co-authored-by: gabrielserrao <[email protected]>
  • Loading branch information
prisae and gabrielserrao authored Jul 11, 2024
1 parent 8bcc59c commit 55ff177
Show file tree
Hide file tree
Showing 7 changed files with 489 additions and 84 deletions.
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,15 @@
__pycache__/
*.nc
*.pdf
*.zip

# Sphinx
docs/_build/
docs/api/resmda*
docs/savefig/
docs/gallery/*
download/
docs/sg_execution_times.rst

# Pytest and coverage related
htmlcov
Expand Down
7 changes: 3 additions & 4 deletions docs/conf.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
import time
import warnings
from resmda import __version__
from sphinx_gallery.sorting import FileNameSortKey

# ==== 1. Extensions ====

Expand Down Expand Up @@ -39,17 +38,17 @@
sphinx_gallery_conf = {
'examples_dirs': ['../examples', ],
'gallery_dirs': ['gallery', ],
'capture_repr': ('_repr_html_', '__repr__'),
'capture_repr': ('_repr_html_', ),
# Patter to search for example files
"filename_pattern": r"\.py",
# Sort gallery example by file name instead of number of lines (default)
"within_subsection_order": FileNameSortKey,
"within_subsection_order": "FileNameSortKey",
# Remove the settings (e.g., sphinx_gallery_thumbnail_number)
'remove_config_comments': True,
# Show memory
'show_memory': True,
# Custom first notebook cell
'first_notebook_cell': '%matplotlib notebook',
'first_notebook_cell': '%matplotlib widget',
}

# https://github.com/sphinx-gallery/sphinx-gallery/pull/521/files
Expand Down
7 changes: 3 additions & 4 deletions examples/basicESMDA.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@
.. 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})
\varepsilon_j - g(x_{j,i})\right) \ , \\
y_{j,i+1} &= g(x_{j,i+1}) \ .
The model used for this example is
Expand Down Expand Up @@ -53,7 +53,7 @@ def forward(x, beta):
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}")
f"{['Linear model', 'Nonlinear model'][i]}: β = {b}")
axs[i].plot(px, forward(px, b))
axs[i].set_xlabel('x')
axs[i].set_ylabel('y')
Expand Down Expand Up @@ -109,7 +109,6 @@ def plot_result(mpost, dpost, dobs, title, ylim):
ax2.set_xlabel('y')
ax2.set_xlim([-5, 8])
ax2.legend()
fig.show()


###############################################################################
Expand Down
85 changes: 40 additions & 45 deletions examples/basicreservoir.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
# ----------------

# Grid extension
nx = 25
nx = 30
ny = 25
nc = nx*ny

Expand All @@ -32,7 +32,6 @@

# ESMDA parameters
ne = 100 # Number of ensembles
# dt = np.logspace(-5, -3, 10)
dt = np.zeros(10)+0.0001 # Time steps (could be irregular, e.g., increasing!)
time = np.r_[0, np.cumsum(dt)]
nt = time.size
Expand Down Expand Up @@ -61,25 +60,23 @@
perm_prior = RP(ne, random=rng)


# TODO: change scale in imshow to represent meters

# QC covariance, reference model, and first two random models
fig, axs = plt.subplots(2, 2, constrained_layout=True)
axs[0, 0].set_title('Model')
im = axs[0, 0].imshow(perm_true.T, vmin=perm_min, vmax=perm_max)
axs[0, 1].set_title('Lower Covariance Matrix')
im2 = axs[0, 1].imshow(RP.cov, cmap='plasma')
axs[1, 0].set_title('Random Model 1')
axs[1, 0].imshow(perm_prior[0, ...].T, vmin=perm_min, vmax=perm_max)
axs[1, 1].set_title('Random Model 2')
axs[1, 1].imshow(perm_prior[1, ...].T, vmin=perm_min, vmax=perm_max)
fig.colorbar(im, ax=axs[1, :], orientation='horizontal',
label='Log of Permeability (mD)')
for ax in axs[1, :]:
ax.set_xlabel('x-direction')
for ax in axs[:, 0]:
ax.set_ylabel('y-direction')
fig.show()
pinp1 = {"origin": "lower", "vmin": perm_min, "vmax": perm_max}
fig, axs = plt.subplots(2, 2, figsize=(6, 6), constrained_layout=True)
axs[0, 0].set_title("Model")
im = axs[0, 0].imshow(perm_true.T, **pinp1)
axs[0, 1].set_title("Lower Covariance Matrix")
im2 = axs[0, 1].imshow(RP.cov, cmap="plasma")
axs[1, 0].set_title("Random Model 1")
axs[1, 0].imshow(perm_prior[0, ...].T, **pinp1)
axs[1, 1].set_title("Random Model 2")
axs[1, 1].imshow(perm_prior[1, ...].T, **pinp1)
fig.colorbar(im, ax=axs[1, :], orientation="horizontal",
label="Log of Permeability (mD)")
for ax in axs[1, :].ravel():
ax.set_xlabel("x-direction")
for ax in axs[:, 0].ravel():
ax.set_ylabel("y-direction")


###############################################################################
Expand All @@ -101,15 +98,14 @@ def sim(x):
data_obs = rng.normal(data_true, dstd)

# QC data and priors
fig, ax = plt.subplots(1, 1)
ax.set_title('Observed and prior data')
ax.plot(time, data_prior.T, color='.6', alpha=0.5)
ax.plot(time, data_true, 'ko', label='True data')
ax.plot(time, data_obs, 'C3o', label='Obs. data')
fig, ax = plt.subplots(1, 1, constrained_layout=True)
ax.set_title("Observed and prior data")
ax.plot(time*24*60*60, data_prior.T, color=".6", alpha=0.5)
ax.plot(time*24*60*60, data_true, "ko", label="True data")
ax.plot(time*24*60*60, data_obs, "C3o", label="Obs. data")
ax.legend()
ax.set_xlabel('Time (???)')
ax.set_ylabel('Pressure (???)')
fig.show()
ax.set_xlabel("Time (s)")
ax.set_ylabel("Pressure (bar)")


###############################################################################
Expand Down Expand Up @@ -143,14 +139,14 @@ def restrict_permeability(x):
# and posterior ensembles to see how the models have been updated.

# Plot posterior
fig, ax = plt.subplots(1, 3, figsize=(12, 5))
ax[0].set_title('Prior Mean')
ax[0].imshow(perm_prior.mean(axis=0).T)
ax[1].set_title('Post Mean')
ax[1].imshow(perm_post.mean(axis=0).T)
ax[2].set_title('"Truth"')
ax[2].imshow(perm_true.T)
fig.show()
fig, ax = plt.subplots(1, 2, figsize=(8, 5), constrained_layout=True)
pinp2 = {"origin": "lower", "vmin": 2.5, "vmax": 3.5}
ax[0].set_title("Prior Mean")
im = ax[0].imshow(perm_prior.mean(axis=0).T, **pinp2)
ax[1].set_title("Post Mean")
ax[1].imshow(perm_post.mean(axis=0).T, **pinp2)
fig.colorbar(im, ax=ax, label="Log of Permeability (mD)",
orientation="horizontal")


###############################################################################
Expand All @@ -162,15 +158,14 @@ def restrict_permeability(x):


# Compare posterior to prior and observed data
fig, ax = plt.subplots(1, 1)
ax.set_title('Prior and posterior data')
ax.plot(time, data_prior.T, color='.6', alpha=0.5)
ax.plot(time, data_post.T, color='C0', alpha=0.5)
ax.plot(time, data_true, 'ko')
ax.plot(time, data_obs, 'C3o')
ax.set_xlabel('Time (???)')
ax.set_ylabel('Pressure (???)')
fig.show()
fig, ax = plt.subplots(1, 1, constrained_layout=True)
ax.set_title("Prior and posterior data")
ax.plot(time*24*60*60, data_prior.T, color=".6", alpha=0.5)
ax.plot(time*24*60*60, data_post.T, color="C0", alpha=0.5)
ax.plot(time*24*60*60, data_true, "ko")
ax.plot(time*24*60*60, data_obs, "C3o")
ax.set_xlabel("Time (s)")
ax.set_ylabel("Pressure (bar)")


###############################################################################
Expand Down
Loading

0 comments on commit 55ff177

Please sign in to comment.