diff --git a/README.rst b/README.rst index f52b810..a5639cf 100644 --- a/README.rst +++ b/README.rst @@ -10,5 +10,6 @@ presented by Emerick and Reynolds, 2013. - **Bug reports:** https://github.com/tuda-geo/resmda/issues -Available through pip: ``pip install resmda``. +Available through pip and conda: +``pip install resmda`` / ``conda install -c conda-forge resmda``. diff --git a/docs/index.rst b/docs/index.rst index a1f0c59..3c5277f 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -68,7 +68,7 @@ resmda Documentation ^^^^^^^^^^^^^ The API reference of resmda includes almost every function and class. - Some of the underlying theory is also described in the docstring. + Some of the underlying theory is also described in the docstrings. +++ diff --git a/docs/manual/about.rst b/docs/manual/about.rst index acd611d..02d99ac 100644 --- a/docs/manual/about.rst +++ b/docs/manual/about.rst @@ -10,15 +10,15 @@ ES-MDA ------ -In the following an introduction to the ES-MDA (Ensemble Smoother with Multiple -Data Assimilation) algorithm following [EmRe13]_: +In the following an introduction to the Ensemble Smoother with Multiple Data +Assimilation (ES-MDA) algorithm following [EmRe13]_: In history-matching problems, it is common to consider solely the parameter-estimation problem and thereby neglecting model uncertainties. Thus, -unlike EnKF, the parameters and states are always consistent (Thulin et al., -2007). This fact helps to explain the better data matches obtained by ES-MDA -compared to EnKF. The analyzed vector of model parameters :math:`m^a` is given -in that case by +unlike with the ensemble Kalman filter (EnKF), the parameters and states are +always consistent (Thulin et al., 2007). This fact helps to explain the better +data matches obtained by ES-MDA compared to EnKF. The analyzed vector of model +parameters :math:`m^a` is given in that case by .. math:: m_j^a = m_j^f + C_\text{MD}^f \left(C_\text{DD}^f + \alpha C_\text{D} diff --git a/docs/manual/installation.rst b/docs/manual/installation.rst index 3d5e40a..78c5663 100644 --- a/docs/manual/installation.rst +++ b/docs/manual/installation.rst @@ -7,6 +7,12 @@ You can install the latest release of resmda simply via ``pip`` pip install resmda +and ``conda`` + +.. code-block:: console + + conda install -c conda-forge resmda + or clone the repository and install it manually with .. code-block:: console diff --git a/docs/manual/license.rst b/docs/manual/license.rst index ac99a66..ad29127 100644 --- a/docs/manual/license.rst +++ b/docs/manual/license.rst @@ -3,8 +3,8 @@ 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. +Everything of resmda is open-source: The code is released under the Apache 2.0 +license, the text (e.g., the manual) is released under the CC-BY 4.0 license. License of Text and Website diff --git a/examples/fluvialreservoir.py b/examples/fluvialreservoir.py index 14fb3ef..6f330c7 100644 --- a/examples/fluvialreservoir.py +++ b/examples/fluvialreservoir.py @@ -161,7 +161,7 @@ def sim(x): 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)) +loc_mat = resmda.localization_matrix(RP.cov, nd_positions, (nx, ny)) ############################################################################### diff --git a/examples/localization.py b/examples/localization.py index fa6ed6f..28d4673 100644 --- a/examples/localization.py +++ b/examples/localization.py @@ -96,7 +96,7 @@ def sim(x): 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)) +loc_mat = resmda.localization_matrix(RP.cov, nd_positions, (nx, ny)) # QC localization matrix fig, ax = plt.subplots(1, 1, constrained_layout=True) diff --git a/requirements-dev.txt b/requirements-dev.txt index 0b5e3aa..d005eae 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -16,6 +16,7 @@ sphinx_gallery>=0.16 pydata_sphinx_theme sphinx_automodapi ipykernel +ipympl pooch # FOR TESTING diff --git a/resmda/__init__.py b/resmda/__init__.py index e3b2ddc..f50e133 100644 --- a/resmda/__init__.py +++ b/resmda/__init__.py @@ -15,8 +15,8 @@ # the License. from resmda import utils -from resmda.utils import Report -from resmda.data_assimilation import esmda, build_localization_matrix +from resmda.utils import localization_matrix, Report +from resmda.data_assimilation import esmda from resmda.reservoir_simulator import Simulator, RandomPermeability @@ -25,8 +25,7 @@ __all__ = ['reservoir_simulator', 'data_assimilation', 'utils', - 'Simulator', 'RandomPermeability', - 'esmda', 'build_localization_matrix', - 'rng', 'Report'] + 'esmda', 'Simulator', 'RandomPermeability', + 'localization_matrix', 'rng', 'Report'] __version__ = utils.__version__ diff --git a/resmda/data_assimilation.py b/resmda/data_assimilation.py index 86712e5..0f192a2 100644 --- a/resmda/data_assimilation.py +++ b/resmda/data_assimilation.py @@ -16,9 +16,9 @@ import numpy as np -from resmda.utils import rng +from resmda import utils -__all__ = ['esmda', 'build_localization_matrix'] +__all__ = ['esmda'] def __dir__(): @@ -47,23 +47,24 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None, Observed data of shape ``(nd)``. sigma : {float, ndarray} Standard deviation(s) of the observation noise. - alphas : {int, ndarray}, default: 4 + alphas : {int, array-like}, default: 4 Inflation factors for ES-MDA. data_prior : ndarray, default: None - Prior data ensemble, of shape (ne, nd). + Prior data ensemble, of shape ``(ne, nd)``. callback_post : function, default: None 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 also ``forward(model_post)``. return_steps : bool, default: False - If true, returns model (and data) of all ES-MDA steps. - Setting ``return_steps`` to True enforces ``return_post_data=True``. + 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. + Seed or random generator for reproducibility; see + :func:`resmda.utils.rng`. localization_matrix : {ndarray, None}, default: None If provided, apply localization to the Kalman gain matrix, of shape - (model-shape, nd). + ``(model-shape, nd)``. Returns @@ -78,6 +79,9 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None, ne = model_prior.shape[0] nd = data_obs.size + # Get random number generator + rng = utils.rng(random) + # Expand sigma if float if np.asarray(sigma).size == 1: sigma = np.zeros(nd) + sigma @@ -106,7 +110,7 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None, # For each ensemble member, perturb the observation vector using # d_uc = d_obs + sqrt(α_i) * C_D^0.5 z_d; z_d ~ N(0, I_N_d) - zd = rng(random).normal(size=(ne, nd)) + zd = rng.normal(size=(ne, nd)) data_pert = data_obs + np.sqrt(alpha) * sigma * zd # == Step (c) of Emerick & Reynolds, 2013 == @@ -166,37 +170,3 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None, return model_post, data_post else: return model_post - - -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 - data positions. - - Parameters - ---------- - cov_matrix : ndarray - 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, - of size (nd, 2). - shape : tuple - Dimensions of the grid (nx, ny). - - Returns - ------- - loc_matrix : ndarray - Localization matrix of shape (nx, ny, nd). - - """ - # Convert 2D positions of data points to 1D indices suitable for accessing - # the covariance matrix - 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, and reshape. - return cov_matrix[:, indices.astype(int)].reshape((*shape, -1), order='F') diff --git a/resmda/reservoir_simulator.py b/resmda/reservoir_simulator.py index 99bf14f..f8de7a2 100644 --- a/resmda/reservoir_simulator.py +++ b/resmda/reservoir_simulator.py @@ -17,9 +17,9 @@ import numpy as np import scipy as sp -from resmda.utils import rng +from resmda import utils -__all__ = ['Simulator', 'RandomPermeability', 'covariance'] +__all__ = ['Simulator', 'RandomPermeability'] def __dir__(): @@ -29,46 +29,62 @@ def __dir__(): class Simulator: """A small 2D Reservoir Simulator. + 2D connection-based, single-phase, single-component system using the Darcy + assumption, employing backward Euler for time discretization and finite + volume for space discretization. It simulates a single-phase fluid (likely + water) with compressibility in a reservoir with constant porosity and a + heterogeneous permeability field. The simulator utilizes SciPy's + ``sparse.linalg.solve`` for pressure solutions and calculates + transmissibilities for inter-block connections. Well modeling is handled + through the Peaceman well model for well indices, with constant pressure + boundary conditions for injection and production wells. The simulation + operates on a 2D grid with user-defined dimensions (nx, ny), uses flexible + time steps, and starts from a specified initial pressure condition. + Physical processes accounted for include fluid density changes with + pressure (modeling a slightly compressible fluid) while assuming constant + fluid viscosity." + Created by following the course **AESM304A - Flow and Simulation of Subsurface processes** at Delft University of Technology (TUD); this particular part was taught by Dr. D.V. Voskov, https://orcid.org/0000-0002-5399-1755. + + Parameters + ---------- + nx, ny : int + Dimension of field (number of cells). + phi : float, default: 0.2 + Porosity (-). + c_f : float, default: 1e-5 + Formation compressibility (1/kPa). + p0 : float, default: 1.0 + Zero pressure (bar). + rho0 : float, default: 1.0 + Fixed density (kg/m3). + mu_w : float, default: 1.0 + Viscosity (cP). + rw : float, default: 0.15 + Well radius (m). + pres_ini : float, default: 150.0 + Initial well pressure (bar). + wells : {ndarray, None}, default: None + Nd array of shape ``(nwells, 3)``, indicating well locations (with cell + indices) and pressure. If None, the default is used, which is + + np.array([[0, 0, 180], [self.nx-1, self.ny-1, 120]]) + + corresponding to a well in the first and in the last cell, with a + pressure of 180 and 120, respectively. + dx, dz : floats, default: 50.0, 10.0 + Cell dimensions in horizontal (dx) and vertical (dz) + directions (m). + """ def __init__(self, nx, ny, phi=0.2, c_f=1e-5, p0=1.0, rho0=1.0, mu_w=1.0, rw=0.15, pres_ini=150.0, wells=None, dx=50.0, dz=10.0): - """Initialize a Simulation instance. - - Parameters - ---------- - nx, ny : int - Dimension of field (number of cells). - phi : float, default: 0.2 - Porosity (-). - c_f : float, default: 1e-5 - Formation compressibility (1/kPa). - p0 : float, default: 1.0 - Initial pressure (bar or kPa?). - rho0 : float, default: 1.0 - Fixed density (kg/m3). - mu_w : float, default: 1.0 - Viscosity (cP - Pa s). - rw : float, default: 0.15 - Well radius (m). - pres_ini : float, default: 150.0 - Initial pressure [?]. - wells : {ndarray, None}, default: None - Nd array of shape (nwells, 3), indicating well locations (with cell - indices) and pressure. If None, the default is used, which is - np.array([[0, 0, 180], [self.nx-1, self.ny-1, 120]]) - corresponding to a well in the first and in the last cell, with a - pressure of 180 and 120, respectively. - dx, dz : floats, default: 50.0, 10.0 - Cell dimensions in horizontal (dx) and vertical (dz) - directions (m). - - """ + """Initialize a Simulation instance.""" self.size = nx*ny self.shape = (nx, ny) @@ -92,9 +108,9 @@ def __init__(self, nx, ny, phi=0.2, c_f=1e-5, p0=1.0, rho0=1.0, mu_w=1.0, if wells is None: # Default wells setup if none provided. Each well is specified by # its grid coordinates followed by its pressure. The first well - # ([0, 0, 180]) is placed at the top-left corner with a pressure of - # 180 units, representing an injection pressure. The second well - # ([self.nx-1, self.ny-1, 120]), is located at the bottom-right + # ([0, 0, 180]) is placed at the bottom-left corner with a pressure + # of 180 units, representing an injection pressure. The second well + # ([self.nx-1, self.ny-1, 120]), is located at the top-right # corner, and has a pressure of 120 units, possibly a lower # pressure or production scenario. self.wells = np.array([[0, 0, 180], [self.nx-1, self.ny-1, 120]]) @@ -127,7 +143,7 @@ def solve(self, pressure, dt): Current pressure state of the reservoir of size ``self.size``. dt : float - Time step for the simulation. + Time step for the simulation (days). Returns ------- @@ -172,8 +188,12 @@ def solve(self, pressure, dt): d[self.locs] += self._add_wells_d # Bring to sparse matrix - offsets = np.array([-self.nx, -1, 0, 1, self.nx]) - data = np.array([mn, m1, d, p1, pn]) + if self.nx > 1: + offsets = np.array([-self.nx, -1, 0, 1, self.nx]) + data = np.array([mn, m1, d, p1, pn]) + else: + offsets = np.array([-1, 0, 1]) + data = np.array([mn+m1, d, p1+pn]) K = sp.sparse.dia_array((data, offsets), shape=(self.size, self.size)) # Solve the system @@ -188,15 +208,15 @@ def __call__(self, perm_fields, dt=np.ones(10)*0.0001, data=False): Parameters ---------- perm_fields : ndarray - Permeability fields to simulate, either of dimension + Permeability fields (mD) to simulate, either of dimension (ne, nx, ny), or of dimension (nx, ny). dt : ndarray, default: np.ones(10)*0.0001 - Time steps to use for simulation. + Time steps to use for simulation (days). data : {False, [ndarray, ndarray]}, default: False - Specific indices [nx, ny] to output data for; if False, return all - data + Specific indices [nx, ny] to output data for; + if False, return all data Returns ------- @@ -232,38 +252,39 @@ def __call__(self, perm_fields, dt=np.ones(10)*0.0001, data=False): class RandomPermeability: - """Return random permeability fields with certain statistical props.""" + """Return random permeability fields with certain statistical props. - def __init__(self, nx, ny, perm_mean, perm_min, perm_max, - length=(10.0, 10.0), theta=45.0, sigma_pr2=1.0, - dtype='float32'): - """Initialize parameters for generating random permeability fields. - Parameters - ---------- - nx, ny : int - Dimensions of the grid. - perm_mean : float - Mean permeability. - perm_min, perm_max : float - Minimum and maximum values for permeability. - length : tuple of two floats, default: (10.0, 10.0) - Length scales for the correlation of permeability. - theta : float, default: 45.0 - Rotation angle for the anisotropy in the permeability field. - sigma_pr2 : float, default: 1.0 - Variance scale for the permeability. - dtype : str, default: 'float32' - Data type for computations, for precision and performance tuning. + Parameters + ---------- + nx, ny : int + Dimensions of the grid. + perm_mean : float + Mean permeability (mD). + perm_min, perm_max : float + Minimum and maximum values for permeability (mD). + length : tuple of two floats, default: (10.0, 10.0) + Length scales for the correlation of permeability (-). + theta : float, default: 45.0 + Rotation angle for the anisotropy in the permeability field (°). + variance : float, default: 1.0 + Variance scale for the permeability. + dtype : str, default: 'float32' + Data type for computations, for precision and performance tuning. - """ - self.nx, self.ny = nx, ny # Grid dimensions - self.nc = nx * ny # Total number of cells - self.perm_mean = perm_mean # Permeability statistics + """ + + def __init__(self, nx, ny, perm_mean, perm_min, perm_max, + length=(10.0, 10.0), theta=45.0, variance=1.0, + dtype='float32'): + """Initialize parameters for generating random permeability fields.""" + self.nx, self.ny = nx, ny # Grid dimensions + self.nc = nx * ny # Total number of cells + self.perm_mean = perm_mean # Permeability statistics self.perm_min, self.perm_max = perm_min, perm_max self.length, self.theta = length, theta # Anisotropy parameters - self.sigma_pr2 = sigma_pr2 # Variance - self.dtype = dtype # Data type + self.variance = variance # Variance + self.dtype = dtype # Data type @property def cov(self): @@ -273,9 +294,9 @@ def cov(self): statistical parameters. """ if not hasattr(self, '_cov'): - self._cov = covariance( + self._cov = utils.gaussian_covariance( nx=self.nx, ny=self.ny, length=self.length, - theta=self.theta, sigma_pr2=self.sigma_pr2, dtype=self.dtype + theta=self.theta, variance=self.variance, dtype=self.dtype ) return self._cov @@ -292,7 +313,7 @@ def lcho(self): def __call__(self, n, perm_mean=None, perm_min=None, perm_max=None, random=None): - """Gerenate n random permeadility fields + """Generate n random permeability fields Generate n random permeability fields using the specified statistical parameters and random seed. @@ -302,9 +323,9 @@ def __call__(self, n, perm_mean=None, perm_min=None, perm_max=None, n : int Number of fields to generate perm_mean : {float, None}, default: None - Mean permeability to override the initialized value. + Mean permeability to override the initialized value (mD). perm_min, perm_max : {float, None}, default: None - Min and max permeability values to clip the fields. + Min and max permeability values to clip the fields (mD). random : {None, int, np.random.Generator}, default: None Seed or random generator for reproducibility. @@ -325,82 +346,10 @@ def __call__(self, n, perm_mean=None, perm_min=None, perm_max=None, # Initialize fields with mean permeability out = np.full((n, self.nx, self.ny), perm_mean, order='F') for i in range(n): - z = rng(random).normal(size=self.nc) # Generate random numbers + z = utils.rng(random).normal(size=self.nc) # Random numbers # Apply the Cholesky transform out[i, ...] += (self.lcho @ z).reshape( (self.nx, self.ny), order='F') # Clip the results to stay within specified bounds return out.clip(perm_min, perm_max) - - -def covariance(nx, ny, length, theta, sigma_pr2, dtype='float32'): - """Return covariance matrix - - Generate covariance matrix based on grid size, anisotropy, and statistical - parameters. - - - Parameters - ---------- - nx, ny : int - Dimensions of the grid. - length : float - Length scales for the correlation of permeability. - theta : float - Rotation angle for the anisotropy in the permeability field. - sigma_pr2 : float - Variance scale for the permeability. - dtype : str, default: 'float32' - Data type for computations. - - - Returns - ------- - cov : ndarray - Covariance matrix for the permeability field. - - """ - nc = nx * ny # Total number of cells - # Precompute cosine and sine of the rotation angle - cost, sint = np.cos(theta), np.sin(theta) - - # 1. Fill the first row of the covariance matrix - tmp1 = np.zeros([nx, nc], dtype=dtype) - for i in range(nx): - tmp1[i, 0] = 1.0 # Set diagonal - for j in range(i+1, nc): - # Distance in the x and y directions - d0 = (j % nx) - i - d1 = (j // nx) - # Rotate coordinates - rot0 = cost*d0 - sint*d1 - rot1 = sint*d0 + cost*d1 - # Calculate the scaled distance - hl = np.sqrt((rot0/length[0])**2 + (rot1/length[1])**2) - - # Sphere formula for covariance, modified for anisotropy - if sigma_pr2: # Non-zero variance scale - if hl <= 1: - tmp1[i, j-i] = sigma_pr2 * (1 - 1.5*hl + hl**3/2) - - else: # Gaspari-Cohn function for smoothness - if hl < 1: - tmp1[i, j-i] = (-(hl**5)/4 + (hl**4)/2 + (hl**3)*5/8 - - (hl**2)*5/3 + 1) - elif hl >= 1 and hl < 2: - tmp1[i, j-i] = ((hl**5)/12 - (hl**4)/2 + (hl**3)*5/8 + - (hl**2)*5/3 - hl*5 + 4 - (1/hl)*2/3) - - # 2. Get the indices of the non-zero columns - ind = np.where(tmp1.sum(axis=0))[0] - - # 3. Expand the non-zero colums ny-times - tmp2 = np.zeros([nc, ind.size], dtype=dtype) - for i, j in enumerate(ind): - n = j//nx - tmp2[:nc-n*nx, i] = np.tile(tmp1[:, j], ny-n) - - # 4. Construct array through sparse diagonal array - cov = sp.sparse.dia_array((tmp2.T, -ind), shape=(nc, nc)) - return cov.toarray() diff --git a/resmda/utils.py b/resmda/utils.py index 281ae74..ad28085 100644 --- a/resmda/utils.py +++ b/resmda/utils.py @@ -17,6 +17,7 @@ from datetime import datetime import numpy as np +import scipy as sp from scooby import Report as ScoobyReport try: @@ -25,13 +26,126 @@ __version__ = 'unknown-'+datetime.today().strftime('%Y%m%d') -__all__ = ['Report', 'rng'] +__all__ = ['gaussian_covariance', 'localization_matrix', 'Report', 'rng'] def __dir__(): return __all__ +def gaussian_covariance(nx, ny, length, theta, variance, dtype='float32'): + """Return covariance matrix with Gaussian properties + + Generate covariance matrix based on grid size, anisotropy, and statistical + parameters. + + + Parameters + ---------- + nx, ny : int + Dimensions of the grid. + length : float + Length scales for the correlation of the property. + theta : float + Rotation angle for the anisotropy in the property field. + variance : float + Variance of the property. + dtype : str, default: 'float32' + Data type for computations. + + + Returns + ------- + cov : ndarray + Covariance matrix for the property field. + + """ + nc = nx * ny # Total number of cells + # Precompute cosine and sine of the rotation angle + cost, sint = np.cos(theta), np.sin(theta) + + # 1. Fill the first row of the covariance matrix + tmp1 = np.zeros([nx, nc], dtype=dtype) + for i in range(nx): + tmp1[i, 0] = 1.0 # Set diagonal + for j in range(i+1, nc): + # Distance in the x and y directions + d0 = (j % nx) - i + d1 = (j // nx) + # Rotate coordinates + rot0 = cost*d0 - sint*d1 + rot1 = sint*d0 + cost*d1 + # Calculate the scaled distance + hl = np.sqrt((rot0/length[0])**2 + (rot1/length[1])**2) + + # Sphere formula for covariance, modified for anisotropy + if variance: # Non-zero variance scale + if hl <= 1: + tmp1[i, j-i] = variance * (1 - 1.5*hl + hl**3/2) + + else: # Gaspari-Cohn function for smoothness + if hl < 1: + tmp1[i, j-i] = (-(hl**5)/4 + (hl**4)/2 + (hl**3)*5/8 - + (hl**2)*5/3 + 1) + elif hl >= 1 and hl < 2: + tmp1[i, j-i] = ((hl**5)/12 - (hl**4)/2 + (hl**3)*5/8 + + (hl**2)*5/3 - hl*5 + 4 - (1/hl)*2/3) + + # 2. Get the indices of the non-zero columns + ind = np.where(tmp1.sum(axis=0))[0] + + # 3. Expand the non-zero colums ny-times + tmp2 = np.zeros([nc, ind.size], dtype=dtype) + for i, j in enumerate(ind): + n = j//nx + tmp2[:nc-n*nx, i] = np.tile(tmp1[:, j], ny-n) + + # 4. Construct array through sparse diagonal array + cov = sp.sparse.dia_array((tmp2.T, -ind), shape=(nc, nc)) + return cov.toarray() + + +def localization_matrix(covariance, data_positions, shape, cov_type='lower'): + """Return a localization matrix + + Build a localization matrix from a full covariance matrix based on specific + data positions. + + Parameters + ---------- + covariance : ndarray + 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, + of size ``(nd, 2)``. + shape : tuple + Dimensions of the grid ``(nx, ny)``. + cov_type : {'lower', 'upper', 'full'}; default: 'lower' + Matrix type of the provided covariance matrix. + + Returns + ------- + loc_matrix : ndarray + Localization matrix of shape ``(nx, ny, nd)``. + + """ + # Convert 2D positions of data points to 1D indices suitable for accessing + # the covariance matrix + indices = np.array( + data_positions[:, 1] * shape[0] + data_positions[:, 0], + ).astype(int) + + # Extract the corresponding columns from the covariance matrix + loc_mat = covariance[:, indices] + if cov_type == 'lower': + loc_mat += np.tril(covariance, -1).T[:, indices] + elif cov_type == 'upper': + loc_mat += np.triu(covariance, 1).T[:, indices] + + # Reshape and return + return loc_mat.reshape((*shape, -1), order='F') + + def rng(random=None): """Module-wide Random Number Generator. @@ -44,7 +158,7 @@ def rng(random=None): - If ``None`` it returns a :func:`numpy.random.default_rng()` instance instantiated on a module level. - If ``int``, it returns a newly created - :func:`numpy.random.default_rng()`` instance, instantiated with + :func:`numpy.random.default_rng()` instance, instantiated with ``int`` as seed. - If it is already a :class:`numpy.random.Generator` instance, it simply returns it. @@ -67,7 +181,10 @@ def rng(random=None): class Report(ScoobyReport): - """Print a Scooby report; see ``scooby.Report()`` for info.""" + """Print a Scooby report. + + For more info, consult https://github.com/banesullivan/scooby/. + """ def __init__(self, **kwargs): """Initiate a scooby.Report instance.""" diff --git a/tests/test_data_assimilation.py b/tests/test_data_assimilation.py new file mode 100644 index 0000000..f136047 --- /dev/null +++ b/tests/test_data_assimilation.py @@ -0,0 +1,94 @@ +import numpy as np +from numpy.testing import assert_allclose + +import resmda +from resmda import data_assimilation + + +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 forward(x, beta): + """Simple model: y = x (1 + β x²) (linear if beta=0).""" + return np.atleast_1d(x * (1 + beta * x**2)) + + +def test_esmda_linear(): + # Use the simple linear ES-MDA example from the gallery as simple test. + xlocation = -1.0 + ne = int(1e7) + obs_std = 1.0 + rng = resmda.utils.rng(1234) # fixed seed for testing + mprior = rng.normal(loc=1.0, scale=obs_std, size=(ne, 1)) + + def lin_fwd(x): + """Linear forward model.""" + return forward(x, beta=0.0) + + l_dobs = lin_fwd(xlocation) + + # Only return final model and data + lm_post, ld_post = resmda.esmda( + model_prior=mprior, + forward=lin_fwd, + data_obs=l_dobs, + sigma=obs_std, + alphas=4, + random=3333, + ) + assert lm_post.shape == (1e7, 1) + assert ld_post.shape == (1e7, 1) + + x, p = pseudopdf(ld_post[:, 0]) + assert_allclose(0.563, np.max(p), atol=0.001) + assert_allclose(0.012, x[np.argmax(p)], atol=0.001) + + # Also return steps + lm_post2, ld_post2 = resmda.esmda( + model_prior=mprior, + forward=lin_fwd, + data_obs=l_dobs, + sigma=obs_std, + alphas=4, + return_steps=True, + random=3333, + ) + assert lm_post2.shape == (5, 1e7, 1) + assert ld_post2.shape == (5, 1e7, 1) + + assert_allclose(lm_post2[-1, :, 0], lm_post[:, 0]) + assert_allclose(ld_post2[-1, :, 0], ld_post[:, 0]) + + def cbp(x): + x[:] /= 100 # Gets the model much narrowed around 0. + + # alpha-array, localization_matrix, callback_post, return only model + lm_post3 = resmda.esmda( + model_prior=mprior, + forward=lin_fwd, + data_obs=l_dobs, + sigma=obs_std, + alphas=[4, 4, 4, 4], + localization_matrix=np.array([[0.5]]), + callback_post=cbp, + return_post_data=False, + random=3333, + ) + assert lm_post3.shape == (1e7, 1) + + ld_post3 = lin_fwd(lm_post3) + x, p = pseudopdf(ld_post3[:, 0]) + assert_allclose(43420257, np.max(p), atol=1) + assert_allclose(0.0, x[np.argmax(p)], atol=1e-8) + + +def test_all_dir(): + assert set(data_assimilation.__all__) == set(dir(data_assimilation)) diff --git a/tests/test_reservoir_simulator.py b/tests/test_reservoir_simulator.py new file mode 100644 index 0000000..5e6bfcb --- /dev/null +++ b/tests/test_reservoir_simulator.py @@ -0,0 +1,151 @@ +import numpy as np +import scipy as sp +from numpy.testing import assert_allclose + +import resmda +from resmda import reservoir_simulator + + +class TestSimulator: + def test_input_and_defaults(self): + RS = reservoir_simulator.Simulator(3, 2) + assert RS.size == 6 + assert RS.shape == (3, 2) + assert RS.nx == 3 + assert RS.ny == 2 + + assert RS.phi == 0.2 + assert RS.c_f == 1e-5 + assert RS.p0 == 1.0 + assert RS.rho0 == 1.0 + assert RS.mu_w == 1.0 + assert RS.rw == 0.15 + assert RS.dx == 50.0 + assert RS.dz == 10.0 + assert RS.pres_ini == 150.0 + assert_allclose([0, 0, 180, 2, 1, 120], RS.wells.ravel()) + volumes = np.ones(6) * 25000 + assert_allclose(RS.volumes, volumes) + assert_allclose(RS._vol_phi_cf, volumes * 0.2e-5) + assert_allclose(RS.locs, [0, 5]) + + wells = np.array([[1, 1, 130], [0, 3, 200]]) + RS = reservoir_simulator.Simulator( + nx=2, ny=4, phi=0.3, c_f=1e-4, p0=0.8, rho0=0.9, mu_w=0.7, + rw=0.13, pres_ini=140.0, dx=30.0, dz=20., wells=wells, + ) + assert RS.size == 8 + assert RS.shape == (2, 4) + assert RS.nx == 2 + assert RS.ny == 4 + + assert RS.phi == 0.3 + assert RS.c_f == 1e-4 + assert RS.p0 == 0.8 + assert RS.rho0 == 0.9 + assert RS.mu_w == 0.7 + assert RS.rw == 0.13 + assert RS.dx == 30.0 + assert RS.dz == 20.0 + assert RS.pres_ini == 140.0 + assert_allclose(wells, RS.wells) + assert_allclose(RS.locs, [3, 6]) + + def test_check_status_quo(self): + # Here we are only checking the Status Quo. + # It would be great to also have some checks with analytical solutions! + nx, ny = 3, 2 + RS = reservoir_simulator.Simulator(nx, ny) + + result = np.array([[ + [150.00000000, 150.00000000], + [150.00000000, 150.00000000], + [150.00000000, 150.00000000]], [ + [154.18781512, 150.91103108], + [150.54640165, 149.45359835], + [149.08896892, 145.81218488]], [ + [158.73203532, 153.71247787], + [151.25155563, 148.74874726], + [146.28778058, 141.26794637] + ]]) + + perm_fields = np.ones((nx, ny)) + dt = np.array([0.001, 0.1]) + out1 = RS(perm_fields, dt=dt) + + assert_allclose(out1, result) + + out2 = RS(np.array([perm_fields, perm_fields]), dt=dt) + assert_allclose(out1, out2[0]) + assert_allclose(out1, out2[1]) + + out3 = RS(perm_fields, dt=dt, data=[1, 1]) + assert_allclose(out1[:, 1, 1], out3) + + def test_check_1d(self): + nx, ny = 1, 10 + wells1 = np.array([[0, 0, 180], [0, 9, 120]]) + wells2 = np.array([[0, 0, 180], [9, 0, 120]]) + perm_fields = np.ones((nx, ny)) + dt = np.array([0.001, 0.1]) + + RS1 = reservoir_simulator.Simulator(nx, ny, wells=wells1) + RS2 = reservoir_simulator.Simulator(ny, nx, wells=wells2) + + out1 = RS1(perm_fields, dt=dt) + out2 = RS2(perm_fields, dt=dt) + + assert_allclose( + out1.reshape(dt.size+1, -1), + out2.reshape(dt.size+1, -1, order='F') + ) + + +class TestRandomPermeability: + def test_input_and_defaults(self): + RP = reservoir_simulator.RandomPermeability(3, 2, 0.5, 0.2, 0.8) + assert RP.nx == 3 + assert RP.ny == 2 + assert RP.nc == 6 + assert RP.perm_mean == 0.5 + assert RP.perm_min == 0.2 + assert RP.perm_max == 0.8 + assert RP.length == (10.0, 10.0) + assert RP.theta == 45.0 + assert RP.variance == 1.0 + assert RP.dtype == 'float32' + + RP = reservoir_simulator.RandomPermeability( + 3, 2, 0.5, 0.2, 0.8, (3.0, 4.0), 0.0, 2.1, 'float64') + assert RP.length == (3.0, 4.0) + assert RP.theta == 0.0 + assert RP.variance == 2.1 + assert RP.dtype == 'float64' + + def test_cov_lcho(self): + # cov is just a call to utils.gaussian_covariance - check. + RP = reservoir_simulator.RandomPermeability(3, 2, 0.5, 0.2, 0.8) + assert_allclose( + RP.cov, + resmda.utils.gaussian_covariance( + 3, 2, RP.length, RP.theta, RP.variance + ), + ) + # lcho is just a call to sp.linalg.cholesky - check. + assert_allclose(RP.lcho, sp.linalg.cholesky(RP.cov, lower=True)) + + def test_call(self): + RP = reservoir_simulator.RandomPermeability(3, 2, 0.5, 0.0, 1.0) + assert_allclose(RP(1, 0.5, 0.5, 0.5), 0.5) + + rng = resmda.utils.rng(4) + result = np.array([[ + [0.00000000, 0.29639514], + [0.00000000, 0.00000000], + [0.83044794, 0.27682457] + ]]) + assert_allclose(RP(1, random=rng), result) + + +def test_all_dir(): + assert set(reservoir_simulator.__all__) == set(dir(reservoir_simulator)) diff --git a/tests/test_utils.py b/tests/test_utils.py index dc5c3d0..84206fc 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,10 +1,90 @@ +import scooby import numpy as np +from numpy.testing import assert_allclose from resmda import utils +def test_gaussian_covariance(): + # Small length, no variance => eye + assert_allclose(utils.gaussian_covariance(4, 4, [0.1, 0.1], 0, 0), + np.eye(4*4)) + + # Small length, some variance => still eye + assert_allclose( + utils.gaussian_covariance( + nx=4, ny=4, length=[0.1, 0.1], theta=0, variance=1 + ), + np.eye(4*4) + ) + + # Simply check some values with variance + x = utils.gaussian_covariance(4, 4, [2, 1], 30, 2) + assert_allclose(np.diag(x), 1.) + assert_allclose(np.diag(x, -1)[3::4], 0) + for i in range(3): + assert_allclose(np.diag(x, -1)[i::4], 0.00024027168) + assert_allclose(np.diag(x, -4), 0.58600724) + for i in [2, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]: + assert_allclose(np.diag(x, -i), 0) + + # Simply check some values without variance + x = utils.gaussian_covariance(4, 4, [1, 2], 0, 0) + assert_allclose(np.diag(x), 1.) + assert_allclose(np.diag(x, -1)[0], 0.20833333) + assert_allclose(np.diag(x, -3)[1], 0.13466999) + assert_allclose(np.diag(x, -4)[0], 0.6848958) + assert_allclose(np.diag(x, -5)[0], 0.13466999) + assert_allclose(np.diag(x, -7)[1], 0.030032475) + assert_allclose(np.diag(x, -8)[0], 0.20833333) + assert_allclose(np.diag(x, -9)[0], 0.030032475) + assert_allclose(np.diag(x, -11)[1], 0.0004445008) + assert_allclose(np.diag(x, -12)[0], 0.016493056) + assert_allclose(np.diag(x, -13)[0], 0.0004445008) + for i in [2, 6, 10, 14, 15]: + assert_allclose(np.diag(x, -i), 0) + + +def test_localization_matrix(): + # [[ 0, 0, 0, 0], + # [ 4, 5, 0, 0], + # [ 8, 9, 10, 0], + # [12, 13, 14, 15]] + tril = np.tril(np.arange(16).reshape(4, 4)) + full = tril + np.tril(tril, -1).T + triu = np.triu(full) + + data_positions = np.array([[1, 0]], dtype=int) + solution = np.array([[[4]], [[5]], [[9]], [[13]]]) + outtril = utils.localization_matrix(tril, data_positions, (4, 1)) + assert_allclose(solution, outtril) + outtriu = utils.localization_matrix(triu, data_positions, (4, 1), 'upper') + assert_allclose(solution, outtriu) + outfull = utils.localization_matrix(full, data_positions, (4, 1), 'full') + assert_allclose(solution, outfull) + + def test_random(): assert isinstance(utils.rng(), np.random.Generator) assert isinstance(utils.rng(11), np.random.Generator) rng = np.random.default_rng() assert rng == utils.rng(rng) + + +def test_Report(capsys): + out, _ = capsys.readouterr() # Empty capsys + + # Reporting is done by the external package scooby. + # We just ensure the shown packages do not change (core and optional). + out1 = utils.Report() + out2 = scooby.Report( + core=['numpy', 'scipy', 'numba', 'resmda'], + optional=['matplotlib', 'IPython'], + ncol=3) + + # Ensure they're the same; exclude time to avoid errors. + assert out1.__repr__()[115:] == out2.__repr__()[115:] + + +def test_all_dir(): + assert set(utils.__all__) == set(dir(utils))