From 5a3e7e6fdf527a53346d4fffbfd4e915d054d5c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Thu, 30 May 2024 14:16:02 +0200 Subject: [PATCH] Localization II (#7) --- docs/_static/style.css | 33 +++++++++++++++++ docs/conf.py | 6 +++ docs/index.rst | 64 +++++++++++++++++++++++++++++--- docs/manual/about.rst | 3 ++ docs/manual/index.rst | 1 + docs/manual/installation.rst | 10 ++++- docs/manual/license.rst | 34 +++++++++++++++++ examples/README.rst | 6 ++- examples/basicESMDA.py | 4 +- examples/basicreservoir.py | 14 ++----- resmda/__init__.py | 6 ++- resmda/data_assimilation.py | 72 ++++++++++++++---------------------- resmda/utils.py | 2 +- 13 files changed, 188 insertions(+), 67 deletions(-) create mode 100644 docs/_static/style.css create mode 100644 docs/manual/license.rst diff --git a/docs/_static/style.css b/docs/_static/style.css new file mode 100644 index 0000000..5d65031 --- /dev/null +++ b/docs/_static/style.css @@ -0,0 +1,33 @@ +/* "emsig"-colours + blue : #447fc1 + blue-lighter : #61a2d8 + blue-darker : #285985 + red : #ef413d +*/ + +/* Button colours */ + +.sd-btn-info { + font-weight: bold; + background: #447fc1 !important; + border-color: #447fc1 !important; +} + +.sd-btn-info:hover { + background: #ef413d !important; + border-color: #ef413d !important; +} + +/* Card headers */ +.sd-card-header { + font-size: var(--pst-font-size-h5); + font-weight: bold; + padding: 1rem 2rem; + text-align: center; +} + +/* fa stuff */ +.fa { + margin: 0 1rem 0 0; + color: #ef413d; +} diff --git a/docs/conf.py b/docs/conf.py index a1c509f..eb39077 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -102,6 +102,7 @@ "logo": { "text": "resmda", }, + 'navigation_with_keys': True, "github_url": "https://github.com/tuda-geo/resmda", # "use_edit_page_button": True, } @@ -116,3 +117,8 @@ html_use_modindex = True html_file_suffix = '.html' htmlhelp_basename = 'resmda' +html_css_files = [ + "style.css", + "https://cdnjs.cloudflare.com/ajax/libs/font-awesome/4.7.0/" + + "css/font-awesome.min.css" +] diff --git a/docs/index.rst b/docs/index.rst index 2e2831d..a1f0c59 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -18,11 +18,63 @@ resmda Documentation api/index -**TODOs** +.. grid:: 1 1 3 3 + :gutter: 2 -- Add tests -- Finish docstrings & code documentation -- Compare with existing esmda's -- docs: + .. grid-item-card:: - - only deploy for versions + :fa:`book;fa-4x` + + User Guide + ^^^^^^^^^^ + + The manual contains installation instructions, some background theory, + important references, license information, and more. + + +++ + + .. button-ref:: manual + :expand: + :color: info + :click-parent: + + To the user guide + + .. grid-item-card:: + + :fa:`image;fa-4x` + + Gallery + ^^^^^^^ + + The gallery contains examples on the usage of resmda, and is generally + the best way to get started. Download them and modify them to your + needs! + + +++ + + .. button-ref:: resmda_gallery + :expand: + :color: info + :click-parent: + + To the gallery + + .. grid-item-card:: + + :fa:`code;fa-4x` + + API reference + ^^^^^^^^^^^^^ + + The API reference of resmda includes almost every function and class. + Some of the underlying theory is also described in the docstring. + + +++ + + .. button-ref:: api + :expand: + :color: info + :click-parent: + + To the API reference diff --git a/docs/manual/about.rst b/docs/manual/about.rst index eccf751..e501e7a 100644 --- a/docs/manual/about.rst +++ b/docs/manual/about.rst @@ -4,9 +4,12 @@ About A simple 2D reservoir simulator and a straight-forward implementation of the basic *Ensemble Smoother with Multiple Data Assimilation* (ES-MDA) algorithm. +.. _esmda: + ES-MDA ------ + In the following an introduction to the ES-MDA (Ensemble Smoother with Multiple Data Assimilation) algorithm following [EmRe13]_: diff --git a/docs/manual/index.rst b/docs/manual/index.rst index 5590f88..105c225 100644 --- a/docs/manual/index.rst +++ b/docs/manual/index.rst @@ -16,4 +16,5 @@ User Guide about installation references + license diff --git a/docs/manual/installation.rst b/docs/manual/installation.rst index e31fe49..2c002e4 100644 --- a/docs/manual/installation.rst +++ b/docs/manual/installation.rst @@ -1,8 +1,16 @@ Installation ============ -You can install resmda simply via ``pip``: +You can install the latest release of resmda simply via ``pip``: .. code-block:: console pip install resmda + +or clone the repository and run within the command + +.. code-block:: console + + make install + +to get the latest version. diff --git a/docs/manual/license.rst b/docs/manual/license.rst new file mode 100644 index 0000000..ac99a66 --- /dev/null +++ b/docs/manual/license.rst @@ -0,0 +1,34 @@ +.. _license: + +License +======= + +Everything of resmda is open-source: Code is released under the Apache 2.0 +license, text (e.g., the manual) is released under the CC-BY 4.0 license. + + +License of Text and Website +--------------------------- + +This work is licensed under Attribution 4.0 International (CC-BY-4.0). + + https://creativecommons.org/licenses/by/4.0/ + + +License of Code +--------------- + +Copyright 2024 Dieter Werthmüller, Gabriel Serrao Seabra + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + https://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. + diff --git a/examples/README.rst b/examples/README.rst index a1496fe..434dae0 100644 --- a/examples/README.rst +++ b/examples/README.rst @@ -1,4 +1,8 @@ -.. _sphx_glr_gallery: +.. _resmda_gallery: +======= Gallery ======= + +:Release: |version| +:Date: |today| diff --git a/examples/basicESMDA.py b/examples/basicESMDA.py index 7550539..f7e525a 100644 --- a/examples/basicESMDA.py +++ b/examples/basicESMDA.py @@ -83,7 +83,7 @@ def plot_result(mpost, dpost, dobs, title, ylim): # Plot Likelihood ax2.plot( - *pseudopdf(resmda.utils.rng().normal(dobs, size=(ne, dobs.size))), + *pseudopdf(resmda.rng.normal(dobs, size=(ne, dobs.size))), 'C2', lw=2, label='Datum' ) @@ -133,7 +133,7 @@ def plot_result(mpost, dpost, dobs, title, ylim): 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)) +mprior = resmda.rng.normal(loc=1.0, scale=obs_std, size=(ne, 1)) ############################################################################### diff --git a/examples/basicreservoir.py b/examples/basicreservoir.py index debfebe..998f73c 100644 --- a/examples/basicreservoir.py +++ b/examples/basicreservoir.py @@ -2,24 +2,16 @@ 2D Reservoir ESMDA example ========================== -Data Assimilation using ESMDA in Reservoir Simulation ------------------------------------------------------ +Ensemble Smoother Multiple Data Assimilation (ES-MDA) in Reservoir Simulation. -*Advanced Data Assimilation using Ensemble Smoother Multiple Data Assimilation -(ESMDA) in Reservoir Simulation.* - -.. math:: - m_j^a = m_j^f + C_{MD}^f (C_{DD}^f + \alpha_i C_D)^{-1} (d_{uc,j} - - d_j^f);\quad \text{for} \quad j=1, 2, \dots, N_e - -- Prior model: M (Ne, Nx, Ny) -- Prior data: D (Ne, Nt) """ 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(1848) # sphinx_gallery_thumbnail_number = 4 diff --git a/resmda/__init__.py b/resmda/__init__.py index 730e088..e3b2ddc 100644 --- a/resmda/__init__.py +++ b/resmda/__init__.py @@ -20,9 +20,13 @@ from resmda.reservoir_simulator import Simulator, RandomPermeability +# Initialize a random number generator. +rng = utils.rng() + + __all__ = ['reservoir_simulator', 'data_assimilation', 'utils', 'Simulator', 'RandomPermeability', 'esmda', 'build_localization_matrix', - 'Report'] + 'rng', 'Report'] __version__ = utils.__version__ diff --git a/resmda/data_assimilation.py b/resmda/data_assimilation.py index 16e3050..86712e5 100644 --- a/resmda/data_assimilation.py +++ b/resmda/data_assimilation.py @@ -18,7 +18,7 @@ from resmda.utils import rng -__all__ = ['esmda'] +__all__ = ['esmda', 'build_localization_matrix'] def __dir__(): @@ -26,12 +26,12 @@ def __dir__(): def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None, - callback_post=None, return_post_data=True, return_steps=False, - random=None, localization_matrix=None): - """ESMDA algorithm (Emerick and Reynolds, 2013) with optional localization. - - ES-MDA as presented by [EmRe13]_. + localization_matrix=None, callback_post=None, return_post_data=True, + return_steps=False, random=None): + """ES-MDA algorithm ([EmRe13]_) with optional localization. + Consult the section :ref:`esmda` in the manual for the theory and more + information about ES-MDA. Parameters ---------- @@ -52,16 +52,18 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None, data_prior : ndarray, default: None Prior data ensemble, of shape (ne, nd). callback_post : function, default: None - Function to be executed after each ESMDA iteration. + Function to be executed after each ES-MDA iteration to the posterior + model, ``callback_post(model_post)``. return_post_data : bool, default: True - If true, returns data + If true, returns also ``forward(model_post)``. return_steps : bool, default: False - If true, returns model (and data) of all ESMDA steps. - Setting ``return_steps`` to True wil enforce ``return_post_data``. + If true, returns model (and data) of all ES-MDA steps. + Setting ``return_steps`` to True enforces ``return_post_data=True``. random : {None, int, np.random.Generator}, default: None Seed or random generator for reproducibility. localization_matrix : {ndarray, None}, default: None - If provided, apply localization to the Kalman gain matrix. + If provided, apply localization to the Kalman gain matrix, of shape + (model-shape, nd). Returns @@ -91,7 +93,7 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None, # Loop over alphas for i, alpha in enumerate(alphas): - print(f"ESMDA step {i+1: 3d}; α={alpha}") + print(f"ES-MDA step {i+1: 3d}; α={alpha}") # == Step (a) of Emerick & Reynolds, 2013 == # Run the ensemble from time zero. @@ -132,7 +134,7 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None, # Apply localization if provided if localization_matrix is not None: - K *= localization_matrix[..., None] + K *= localization_matrix # Update the ensemble parameters model_post += np.moveaxis(K @ (data_pert - data_prior).T, -1, 0) @@ -166,29 +168,7 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None, return model_post -def convert_positions_to_indices(positions, grid_dimensions): - """Convert 2D grid positions to 1D indices assuming zero-indexed positions. - - Parameters - ---------- - positions : ndarray - Array of (x, y) positions in the grid. - grid_dimensions : tuple - Dimensions of the grid (nx, ny). - - Returns - ------- - indices : ndarray - Array of indices corresponding to the positions in a flattened array. - - """ - nx, ny = grid_dimensions - # Ensure the positions are zero-indexed and correctly calculated for - # row-major order - return positions[:, 1] * nx + positions[:, 0] - - -def build_localization_matrix(cov_matrix, data_positions, grid_dimensions): +def build_localization_matrix(cov_matrix, data_positions, shape): """Build a localization matrix Build a localization matrix from a full covariance matrix based on specific @@ -197,22 +177,26 @@ def build_localization_matrix(cov_matrix, data_positions, grid_dimensions): Parameters ---------- cov_matrix : ndarray - The full (nx*ny) x (nx*ny) covariance matrix. + The lower triangular covariance matrix (nx*ny, nx*ny). data_positions : ndarray - Positions in the grid for each data point (e.g., wells), zero-indexed. - grid_dimensions : tuple + Positions in the grid for each data point (e.g., wells), zero-indexed, + of size (nd, 2). + shape : tuple Dimensions of the grid (nx, ny). Returns ------- loc_matrix : ndarray - Localization matrix of shape (nx*ny, number of data positions). + Localization matrix of shape (nx, ny, nd). """ # Convert 2D positions of data points to 1D indices suitable for accessing # the covariance matrix - indices = convert_positions_to_indices( - data_positions, grid_dimensions).astype(int) + indices = data_positions[:, 1] * shape[0] + data_positions[:, 0] + + # Create full matrix from lower triangular matrix + cov_matrix = cov_matrix + np.tril(cov_matrix, -1).T + # Extract the columns from the covariance matrix corresponding to each data - # point's position - return cov_matrix[:, indices] + # point's position, and reshape. + return cov_matrix[:, indices.astype(int)].reshape((*shape, -1), order='F') diff --git a/resmda/utils.py b/resmda/utils.py index b617ef2..281ae74 100644 --- a/resmda/utils.py +++ b/resmda/utils.py @@ -35,7 +35,7 @@ def __dir__(): def rng(random=None): """Module-wide Random Number Generator. - Mainly meant for internal use. + Instantiate a random number generator. Parameters