.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "gallery/fluvialreservoir.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_gallery_fluvialreservoir.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_gallery_fluvialreservoir.py:


2D Fluvial Reservoir ESMDA example
==================================

This example uses fluvial models containing different facies.

This in contrast to the basic reservoir example
:ref:`sphx_glr_gallery_basicreservoir.py`, where a single facies was used. The
example also compares the use of ESMDA with and without localization, as
explained in the example :ref:`sphx_glr_gallery_localization.py`.

The fluvial models were generated with ``FLUVSIM`` through ``geomodpy``, for
more information see towards the end of the example where the code is shown to
reproduce the facies.

.. note::

    To retrieve the data, you need to have ``pooch`` installed:

    .. code-block:: bash

        pip install pooch

    or

    .. code-block:: bash

        conda install -c conda-forge pooch

.. GENERATED FROM PYTHON SOURCE LINES 31-44

.. code-block:: Python

    import json

    import pooch
    import numpy as np
    import matplotlib.pyplot as plt

    import dageo

    # For reproducibility, we instantiate a random number generator with a fixed
    # seed. For production, remove the seed!
    rng = np.random.default_rng(1513)









.. GENERATED FROM PYTHON SOURCE LINES 46-48

Load and plot the facies
------------------------

.. GENERATED FROM PYTHON SOURCE LINES 48-77

.. code-block:: Python


    folder = "data"
    ffacies = "facies.npy"
    finput = "facies.json"
    fpfacies = pooch.retrieve(
        "https://raw.github.com/tuda-geo/data/2024-06-18/resmda/"+ffacies,
        "4bfe56c836bf17ca63453c37e5da91cb97bbef8cc6c08d605f70bd64fe7488b2",
        fname=ffacies,
        path=folder,
    )
    facies = np.load(fpfacies)
    ne, nx, ny = facies.shape

    # Define mean permeability per facies
    perm_means = [0.1, 5.0, 3.0]

    # Plot the facies
    fig, axs = plt.subplots(
        2, 5, figsize=(8, 3), sharex=True, sharey=True, constrained_layout=True)
    axs = axs.ravel()
    fig.suptitle(f"Facies {[f'{i} = {p}' for i, p in enumerate(perm_means)]}")
    for i in range(ne):
        im = axs[i].imshow(
            facies[i, ...], cmap=plt.get_cmap("Accent", 3),
            clim=[-0.5, 2.5], origin="lower"
        )
    fig.colorbar(im, ax=axs, ticks=[0, 1, 2], label="Facies code")





.. image-sg:: /gallery/images/sphx_glr_fluvialreservoir_001.png
   :alt: Facies ['0 = 0.1', '1 = 5.0', '2 = 3.0']
   :srcset: /gallery/images/sphx_glr_fluvialreservoir_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Downloading data from 'https://raw.github.com/tuda-geo/data/2024-06-18/resmda/facies.npy' to file '/home/runner/work/dageo/dageo/examples/data/facies.npy'.




.. GENERATED FROM PYTHON SOURCE LINES 78-80

Assign random permeabilities to the facies
------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 80-104

.. code-block:: Python


    perm_min = 0.05
    perm_max = 10.0

    # Instantiate a random permeability instance
    RP = dageo.RandomPermeability(
        nx, ny, perm_mean=None, perm_min=perm_min, perm_max=perm_max
    )

    # Fill the different facies with their permeabilities
    permeabilities = np.empty(facies.shape)
    for code, mean in enumerate(perm_means):
        mask = facies == code
        permeabilities[mask] = RP(ne, perm_mean=mean)[mask]

    fig, axs = plt.subplots(
        2, 5, figsize=(8, 3), sharex=True, sharey=True, constrained_layout=True)
    axs = axs.ravel()
    fig.suptitle("Permeabilities")
    for i in range(ne):
        im = axs[i].imshow(permeabilities[i, ...], origin="lower")
    fig.colorbar(im, ax=axs, label="Log Permeability (mD)")





.. image-sg:: /gallery/images/sphx_glr_fluvialreservoir_002.png
   :alt: Permeabilities
   :srcset: /gallery/images/sphx_glr_fluvialreservoir_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 105-107

Model parameters
----------------

.. GENERATED FROM PYTHON SOURCE LINES 107-135

.. code-block:: Python


    # We take the first model as "true/reference", and the other for ESMDA.
    perm_true = permeabilities[0, ...][None, ...]
    perm_prior = permeabilities[1:, ...]

    # Time steps
    dt = np.zeros(10)+0.0001
    time = np.r_[0, np.cumsum(dt)]
    nt = time.size

    # Assumed standard deviation of our data
    dstd = 0.5

    # Measurement points
    ox = (5, 15, 24)
    oy = (5, 10, 24)

    # Number of data points
    nd = time.size * len(ox)

    # Wells
    wells = np.array([
        [ox[0], oy[0], 180], [5, 12, 120],
        [ox[1], oy[1], 180], [20, 5, 120],
        [ox[2], oy[2], 180], [24, 17, 120]
    ])









.. GENERATED FROM PYTHON SOURCE LINES 136-138

Run the prior models and the reference case
-------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 138-155

.. code-block:: Python


    # Instantiate reservoir simulator
    RS = dageo.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]









.. GENERATED FROM PYTHON SOURCE LINES 156-158

Localization Matrix
-------------------

.. GENERATED FROM PYTHON SOURCE LINES 158-166

.. code-block:: Python


    # 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 = dageo.localization_matrix(RP.cov, nd_positions, (nx, ny))









.. GENERATED FROM PYTHON SOURCE LINES 167-169

ESMDA
-----

.. GENERATED FROM PYTHON SOURCE LINES 169-188

.. code-block:: Python



    def restrict_permeability(x):
        """Restrict possible permeabilities."""
        np.clip(x, perm_min, perm_max, out=x)


    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,
    }









.. GENERATED FROM PYTHON SOURCE LINES 189-191

Without localization
''''''''''''''''''''

.. GENERATED FROM PYTHON SOURCE LINES 191-195

.. code-block:: Python


    nl_perm_post, nl_data_post = dageo.esmda(**inp)






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    ESMDA step   1; α=4.0
    ESMDA step   2; α=4.0
    ESMDA step   3; α=4.0
    ESMDA step   4; α=4.0




.. GENERATED FROM PYTHON SOURCE LINES 196-198

With localization
'''''''''''''''''

.. GENERATED FROM PYTHON SOURCE LINES 198-202

.. code-block:: Python


    wl_perm_post, wl_data_post = dageo.esmda(**inp, localization_matrix=loc_mat)






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    ESMDA step   1; α=4.0
    ESMDA step   2; α=4.0
    ESMDA step   3; α=4.0
    ESMDA step   4; α=4.0




.. GENERATED FROM PYTHON SOURCE LINES 203-205

Compare permeabilities
----------------------

.. GENERATED FROM PYTHON SOURCE LINES 205-234

.. code-block:: Python


    # Plot posterior
    fig, axs = plt.subplots(
        1, 3, figsize=(8, 4), sharex=True, sharey=True, constrained_layout=True)

    par = {"vmin": perm_min, "vmax": perm_max, "origin": "lower"}

    axs[0].set_title("Prior Mean")
    im = axs[0].imshow(perm_prior.mean(axis=0).T, **par)


    axs[1].set_title("Post Mean; No localization")
    axs[1].imshow(nl_perm_post.mean(axis=0).T, **par)

    axs[2].set_title("Post Mean: Localization")
    axs[2].imshow(wl_perm_post.mean(axis=0).T, **par)
    axs[2].contour(loc_mat.sum(axis=2).T, levels=[2.0, ], colors="w")

    fig.colorbar(im, ax=axs, label="Log Permeabilities (mD)",
                 orientation="horizontal")

    for ax in axs:
        for well in wells:
            ax.plot(well[0], well[1], ["C3v", "C1^"][int(well[2] == 120)])
    for ax in axs:
        ax.set_xlabel('x-direction')
    axs[0].set_ylabel('y-direction')





.. image-sg:: /gallery/images/sphx_glr_fluvialreservoir_003.png
   :alt: Prior Mean, Post Mean; No localization, Post Mean: Localization
   :srcset: /gallery/images/sphx_glr_fluvialreservoir_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 235-237

Compare data
------------

.. GENERATED FROM PYTHON SOURCE LINES 237-259

.. code-block:: Python


    # 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.set_title(f"Well ({ox[i]}; {oy[i]})")
        ax.plot(time*24*60*60, data_prior[..., i::3].T, color=".6", alpha=0.5)
        ax.plot(time*24*60*60, nl_data_post[..., i::3].T, color="C0", alpha=0.5)
        ax.plot(time*24*60*60, data_obs[0, i::3], "C3o")
    for i, ax in enumerate(axs[1, :]):
        ax.plot(time*24*60*60, data_prior[..., i::3].T, color=".6", alpha=0.5)
        ax.plot(time*24*60*60, wl_data_post[..., i::3].T, color="C0", alpha=0.5)
        ax.plot(time*24*60*60, data_obs[0, i::3], "C3o")
        ax.set_xlabel("Time (s)")
    for i, ax in enumerate(axs[:, 0]):
        ax.set_ylabel("Pressure (bar)")
    for i, txt in enumerate(["No l", "L"]):
        axs[i, 2].yaxis.set_label_position("right")
        axs[i, 2].set_ylabel(f"{txt}ocalization")





.. image-sg:: /gallery/images/sphx_glr_fluvialreservoir_004.png
   :alt: Prior and posterior data, Well (5; 5), Well (15; 10), Well (24; 24)
   :srcset: /gallery/images/sphx_glr_fluvialreservoir_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 260-372

Reproduce the facies
--------------------

.. note::

    The following cell is about how to reproduce the facies data loaded in
    ``facies.npy``. This was created using ``geomodpy``.

    ``geomodpy`` (Guillaume Rongier, 2023) is not open-source yet. The
    functionality of geomodpy that we use here is a python wrapper for the
    Fortran library ``FLUVSIM`` published in:

        **Deutsch, C. V., and T. T. Tran**, 2002, FLUVSIM: a program for
        object-based stochastic modeling of fluvial depositional systems:
        Computers & Geosciences, 28(4), 525--535.

        DOI: `10.1016/S0098-3004(01)00075-9
        <https://doi.org/10.1016/S0098-3004(01)00075-9>`_.


.. code-block:: python

    # ==== Required imports ====
    import json
    import numpy as np

    # FLUVSIM Version used: 2.900
    from geomodpy.wrapper.fluvsim import FLUVSIM

    # For reproducibility, we instantiate a random number generator with a
    # fixed seed. For production, remove the seed!
    rng = np.random.default_rng(1848)


    # ==== Define the geological parameters ====

    # Here we define the geological parameters by means of their normal
    # distribution parameters

    # Each tuple stands for (mean, std); lists contain several of them.
    geol_distributions = {
        "channel_orientation": (60, 20),
        "channel_amplitude": [(100, 1), (250, 1), (400, 1)],
        "channel_wavelength": [(1000, 5), (2000, 5), (3000, 5)],
        "channel_thickness": [(4, 0.1), (8, 0.1), (11, 0.1)],
        "channel_thickness_undulation": (1, 0.02),
        "channel_thickness_undulation_wavelength": [
            (250, 1), (400, 1), (450, 1)
        ],
        "channel_width_thickness_ratio": [(40, 0.5), (50, 0.5), (60, 0.5)],
        "channel_width_undulation": (1, 0.02),
        "channel_width_undulation_wavelength": (250, 1),
        "channel_prop": (0.4, 0.005),
    }


    def generate_geol_params(geol_dists):
        """Generate geological parameters from normal distributions.

        Expects for each parameter a tuple of two values, or a list
        containing tuples of two values each.
        """
        geol_params = {}
        for param, dist in geol_dists.items():
            if isinstance(dist, list):
                geol_params[param] = tuple(
                    rng.normal(mean, std) for mean, std in dist
                )
            else:
                geol_params[param] = rng.normal(*dist)
        return geol_params


    # ==== Create the facies ====

    # Number of sets and realizations
    nsets = 2
    nreal = 5

    # Model dimension
    nx, ny, nz = 64, 64, 1

    # Pre-allocate containers to store all realizations and their
    # corresponding parameters
    all_params = {}
    facies = np.zeros((nsets * nreal, nz, nx, ny), dtype="i4")

    for i in range(nsets):  # We create two sets
        print(f"Generating realization {i+1} of {nsets}")

        params = generate_geol_params(geol_distributions)
        all_params[f"set-{i}"] = params

        fluvsim = FLUVSIM(
            shape=(nx, ny, nz),
            spacing=(50, 50, 1),
            origin=(25, 25, 0.5),
            n_realizations=nreal,
            **params
        )

        realizations = fluvsim.run().data_vars["facies code"].values
        facies[i*nreal:(i+1)*nreal, ...] = realizations.astype("i4")


    # ==== Save the outputs ====

    # Save the input parameters to FLUVSIM as a json.
    with open("facies.json", "w") as f:
        json.dump(all_params, f, indent=2)
    # Save the facies values as a compressed npy-file.
    np.save("facies.npy", facies.squeeze(), allow_pickle=False)

.. GENERATED FROM PYTHON SOURCE LINES 375-380

Input parameters to ``FLUVSIM``
-------------------------------
For reproducibility purposes we print here the used input values to FLUVSIM.
These are, just as the data themselves, online at
https://github.com/tuda-geo/data/resmda.

.. GENERATED FROM PYTHON SOURCE LINES 380-390

.. code-block:: Python


    fpinput = pooch.retrieve(
        "https://raw.github.com/tuda-geo/data/2024-06-18/resmda/"+finput,
        "db2cb8a620775c68374c24a4fa811f6350381c7fc98a823b9571136d307540b4",
        fname=finput,
        path=folder,
    )
    with open(fpinput, "r") as f:
        print(json.dumps(json.load(f), indent=2))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Downloading data from 'https://raw.github.com/tuda-geo/data/2024-06-18/resmda/facies.json' to file '/home/runner/work/dageo/dageo/examples/data/facies.json'.
    {
      "set-0": {
        "channel_orientation": 63.93933122172756,
        "channel_amplitude": [
          99.72003201597424,
          250.09738942599537,
          399.99645679311095
        ],
        "channel_wavelength": [
          998.3771346070469,
          2006.0222525108188,
          3002.013904724877
        ],
        "channel_thickness": [
          4.015660555164119,
          8.02125812391539,
          10.807448731138143
        ],
        "channel_thickness_undulation": 1.006385274247584,
        "channel_thickness_undulation_wavelength": [
          250.36201979679302,
          398.4759608983245,
          448.44486265914315
        ],
        "channel_width_thickness_ratio": [
          40.545201273906024,
          50.53224723965168,
          60.21913791333122
        ],
        "channel_width_undulation": 1.021202467139173,
        "channel_width_undulation_wavelength": 249.2684837598085,
        "channel_prop": 0.40189379634538847
      },
      "set-1": {
        "channel_orientation": 45.94027637621462,
        "channel_amplitude": [
          98.16203534405199,
          251.1448878929924,
          400.3245308818956
        ],
        "channel_wavelength": [
          1009.8934540579546,
          2008.4295824697977,
          3003.9136764267087
        ],
        "channel_thickness": [
          4.111706626685296,
          8.002547496480005,
          10.775694168554176
        ],
        "channel_thickness_undulation": 1.0147354365869505,
        "channel_thickness_undulation_wavelength": [
          249.9803344958115,
          398.8015770326054,
          448.3657322896644
        ],
        "channel_width_thickness_ratio": [
          39.87761151661253,
          49.4521315931951,
          58.02613825217046
        ],
        "channel_width_undulation": 0.9945245219176657,
        "channel_width_undulation_wavelength": 250.56212687685527,
        "channel_prop": 0.4022891542425987
      }
    }




.. GENERATED FROM PYTHON SOURCE LINES 391-392

.. code-block:: Python

    dageo.Report()





.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <table style='border: 1.5px solid;'>
      <tr>
         <td style='text-align: center; font-weight: bold; font-size: 1.2em; border: 1px solid;' colspan='6'>Fri Dec 06 15:58:54 2024 UTC</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>OS</td>
        <td style='text-align: left; border: 1px solid;'>Linux (Ubuntu 22.04)</td>
        <td style='text-align: right; border: 1px solid;'>CPU(s)</td>
        <td style='text-align: left; border: 1px solid;'>4</td>
        <td style='text-align: right; border: 1px solid;'>Machine</td>
        <td style='text-align: left; border: 1px solid;'>x86_64</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>Architecture</td>
        <td style='text-align: left; border: 1px solid;'>64bit</td>
        <td style='text-align: right; border: 1px solid;'>RAM</td>
        <td style='text-align: left; border: 1px solid;'>15.6 GiB</td>
        <td style='text-align: right; border: 1px solid;'>Environment</td>
        <td style='text-align: left; border: 1px solid;'>Python</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>File system</td>
        <td style='text-align: left; border: 1px solid;'>ext4</td>
      </tr>
      <tr>
         <td style='text-align: center; border: 1px solid;' colspan='6'>Python 3.12.7 (main, Oct  1 2024, 15:17:55) [GCC 11.4.0]</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>dageo</td>
        <td style='text-align: left; border: 1px solid;'>1.1.0</td>
        <td style='text-align: right; border: 1px solid;'>numpy</td>
        <td style='text-align: left; border: 1px solid;'>2.1.3</td>
        <td style='text-align: right; border: 1px solid;'>scipy</td>
        <td style='text-align: left; border: 1px solid;'>1.14.1</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>matplotlib</td>
        <td style='text-align: left; border: 1px solid;'>3.9.3</td>
        <td style='text-align: right; border: 1px solid;'>IPython</td>
        <td style='text-align: left; border: 1px solid;'>8.30.0</td>
        <td style= border: 1px solid;'></td>
        <td style= border: 1px solid;'></td>
      </tr>
    </table>
    </div>
    <br />
    <br />


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 16.255 seconds)

**Estimated memory usage:**  349 MB


.. _sphx_glr_download_gallery_fluvialreservoir.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: fluvialreservoir.ipynb <fluvialreservoir.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: fluvialreservoir.py <fluvialreservoir.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: fluvialreservoir.zip <fluvialreservoir.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_