From 76a470aed5e52af0d4936cb78db78845b2990571 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Fri, 31 May 2024 14:38:44 +0200 Subject: [PATCH 1/3] Add localization example; publish PR's --- .github/workflows/documentation.yml | 10 ++ examples/localization.py | 192 ++++++++++++++++++++++++++++ 2 files changed, 202 insertions(+) create mode 100644 examples/localization.py diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 63a5177..386ac29 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -45,3 +45,13 @@ jobs: github_token: ${{ secrets.GITHUB_TOKEN }} publish_dir: ./docs/_build/html/ force_orphan: true + + - name: Deploy to GitHub Pages + uses: peaceiris/actions-gh-pages@v4 + if: success() && github.event_name == 'pull_request' + with: + publish_branch: gh-pull-request + github_token: ${{ secrets.GITHUB_TOKEN }} + publish_dir: ./docs/_build/html/ + force_orphan: true + destination_dir: pr diff --git a/examples/localization.py b/examples/localization.py new file mode 100644 index 0000000..d1e2322 --- /dev/null +++ b/examples/localization.py @@ -0,0 +1,192 @@ +r""" +Localization +========================== + +This example follows contextually :ref:`sphx_glr_gallery_basicreservoir.py`, +but uses several well doublets and compares ES-MDA with and without +localization. +""" +import numpy as np +import matplotlib.pyplot as plt + +import resmda + +# For reproducibility, we instantiate a random number generator with a fixed +# seed. For production, remove the seed! +rng = np.random.default_rng(2020) + +# sphinx_gallery_thumbnail_number = 1 + +############################################################################### +# Model parameters +# ---------------- + +# Grid extension +nx = 30 +ny = 30 +nc = nx*ny + +# Permeabilities +perm_mean = 3.0 +perm_min = 0.5 +perm_max = 5.0 + +# ESMDA parameters +ne = 100 # Number of ensembles +dt = np.zeros(10)+0.0001 # Time steps (could be irregular, e.g., increasing!) +time = np.r_[0, np.cumsum(dt)] + +# Assumed sandard deviation of our data +dstd = 0.5 + +# Observation location indices (should be well locations) +ox = (5, 15, 24) +oy = (5, 10, 24) + +# Number of data points +nd = time.size * len(ox) + +# Wells +wells = np.array([ + [5, 5, 180], [5, 12, 120], + [15, 10, 180], [20, 5, 120], + [24, 24, 180], [24, 17, 120] +]) + + +############################################################################### +# Create permeability maps for ESMDA +# ---------------------------------- +# +# We will create a set of permeability maps that will serve as our initial +# guess (prior). These maps are generated using a Gaussian random field and are +# constrained by certain statistical properties. + +# Get the model and ne prior models +RP = resmda.RandomPermeability(nx, ny, perm_mean, perm_min, perm_max) +perm_true = RP(1, random=rng) +perm_prior = RP(ne, random=rng) + + +############################################################################### +# Run the prior models and the reference case +# ------------------------------------------- + +# Instantiate reservoir simulator +RS = resmda.Simulator(nx, ny, wells=wells) + + +def sim(x): + """Custom fct to use exp(x), and specific dt & location.""" + return RS(np.exp(x), dt=dt, data=(ox, oy)).reshape((x.shape[0], -1)) + + +# Simulate data for the prior and true fields +data_prior = sim(perm_prior) +data_true = sim(perm_true) +data_obs = rng.normal(data_true, dstd) +data_obs[0, :3] = data_true[0, :3] +print(data_obs.shape) + + +############################################################################### +# Localization Matrix +# ------------------- + +# Vector of nd length with the well x and y position for each nd data point +nd_positions = np.tile(np.array([ox, oy]), time.size).T + +# Create matrix +loc_mat = resmda.build_localization_matrix(RP.cov, nd_positions, (nx, ny)) + +# QC localization matrix +fig, ax = plt.subplots(1, 1, sharex=True, sharey=True, constrained_layout=True) +ax.imshow(loc_mat.sum(axis=2).T) +ax.contour(loc_mat.sum(axis=2).T, levels=[2.0, ], colors='w') +for well in wells: + ax.plot(well[0], well[1], ['C3v', 'C1^'][int(well[2] == 120)]) + + +############################################################################### +# ES-MDA +# ------ + + +def restrict_permeability(x): + x[:] = np.clip(x, perm_min, perm_max) + + +inp = { + 'model_prior': perm_prior, + 'forward': sim, + 'data_obs': data_obs, + 'sigma': dstd, + 'alphas': 4, + 'data_prior': data_prior, + 'callback_post': restrict_permeability, + 'random': rng, +} + + +############################################################################### +# Without localization +# '''''''''''''''''''' + +nl_perm_post, nl_data_post = resmda.esmda(**inp) + + +############################################################################### +# With localization +# ''''''''''''''''' + +wl_perm_post, wl_data_post = resmda.esmda(**inp, localization_matrix=loc_mat) + + +############################################################################### +# Compare permeabilities +# ---------------------- + +# Plot posterior +fig, axs = plt.subplots( + 2, 2, figsize=(6, 6), sharex=True, sharey=True, constrained_layout=True) +axs[0, 0].set_title('Prior Mean') +axs[0, 0].imshow(perm_prior.mean(axis=0).T) +axs[0, 1].set_title('"Truth"') +axs[0, 1].imshow(perm_true.T) + + +axs[1, 0].set_title('Post Mean without localization') +axs[1, 0].imshow(nl_perm_post.mean(axis=0).T) +axs[1, 1].set_title('Post Mean with localization') +axs[1, 1].imshow(wl_perm_post.mean(axis=0).T) +axs[1, 1].contour(loc_mat.sum(axis=2).T, levels=[2.0, ], colors='w') + +for ax in axs.ravel(): + for well in wells: + ax.plot(well[0], well[1], ['C3v', 'C1^'][int(well[2] == 120)]) + + +############################################################################### +# Compare data +# ------------ + +# QC data and priors +fig, axs = plt.subplots( + 2, 3, figsize=(8, 5), sharex=True, sharey=True, constrained_layout=True) +fig.suptitle('Prior and posterior data') +for i, ax in enumerate(axs[0, :]): + ax.plot(time, data_prior[..., i::3].T, color='.6', alpha=0.5) + ax.plot(time, nl_data_post[..., i::3].T, color='C0', alpha=0.5) + ax.plot(time, data_obs[0, i::3], 'C3o') + ax.set_ylabel('Pressure') +for i, ax in enumerate(axs[1, :]): + ax.plot(time, data_prior[..., i::3].T, color='.6', alpha=0.5) + ax.plot(time, wl_data_post[..., i::3].T, color='C0', alpha=0.5) + ax.plot(time, data_obs[0, i::3], 'C3o') + ax.set_xlabel('Time') + ax.set_ylabel('Pressure') + + +############################################################################### + +resmda.Report() From 1be8283698b431ffdbd286417f3c04069adfa163 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Fri, 31 May 2024 14:51:31 +0200 Subject: [PATCH 2/3] Second try --- .github/workflows/documentation.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 386ac29..8f7a170 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -50,7 +50,6 @@ jobs: uses: peaceiris/actions-gh-pages@v4 if: success() && github.event_name == 'pull_request' with: - publish_branch: gh-pull-request github_token: ${{ secrets.GITHUB_TOKEN }} publish_dir: ./docs/_build/html/ force_orphan: true From 8545a8e51b0cce9e4b98cf819224708f57a8b45c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Fri, 31 May 2024 14:54:23 +0200 Subject: [PATCH 3/3] Give up for now --- .github/workflows/documentation.yml | 9 --------- 1 file changed, 9 deletions(-) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 8f7a170..63a5177 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -45,12 +45,3 @@ jobs: github_token: ${{ secrets.GITHUB_TOKEN }} publish_dir: ./docs/_build/html/ force_orphan: true - - - name: Deploy to GitHub Pages - uses: peaceiris/actions-gh-pages@v4 - if: success() && github.event_name == 'pull_request' - with: - github_token: ${{ secrets.GITHUB_TOKEN }} - publish_dir: ./docs/_build/html/ - force_orphan: true - destination_dir: pr