From 463a5f31bad3c43c09197da04c64769a93001fa1 Mon Sep 17 00:00:00 2001 From: gabrielserrao Date: Wed, 12 Jun 2024 16:04:27 +0200 Subject: [PATCH 01/13] adding post perm plot --- examples/fluvialExample.py | 257 +++++++++++++++++++++++++++++++++++++ 1 file changed, 257 insertions(+) create mode 100644 examples/fluvialExample.py diff --git a/examples/fluvialExample.py b/examples/fluvialExample.py new file mode 100644 index 0000000..19918a8 --- /dev/null +++ b/examples/fluvialExample.py @@ -0,0 +1,257 @@ +#%% +r""" +2D Fluvial Reservoir ESMDA example +========================== + +Ensemble Smoother Multiple Data Assimilation (ES-MDA) in Reservoir Simulation. + +""" +import numpy as np +import matplotlib.pyplot as plt +import xarray as xr + +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 + +############################################################################### +# Model parameters +# ---------------- +# Permeabilities +#read xarray from file +''' +xarray with 10000 realizations of facies for channelized reservoirs +all the models are 64x64x1 with 100x100x1 spacing (dx,dy,dz) which means they are pseudo 3D models + +We will read the first 1001 realizations and generate the permeability maps for the first 1000 realizations and +the first will be the reference case + +The property of interest is called facies code and we will assign a fixed permeability to each facies code + +Later, we can distribute the permeability across the reservoir using the facies code and the variance approach as in the basis case + +''' +facies_ensemble = xr.open_dataset('realizations.nc') + +reference_case = facies_ensemble.isel(Realization=0) +prior_ensemble = facies_ensemble.isel(Realization=slice(1,1000)) + +fac_array = reference_case['facies code'].values.astype(float) +fac_prior_array = prior_ensemble['facies code'].values.astype(float) +fac_prior_array = np.squeeze(fac_prior_array) #remove an extra dimension + +#assign bounds for the permeability +perm_min = 0.1 +perm_max =5 +# Grid extension +nx = fac_array.shape[1] +ny = fac_array.shape[2] +nc = nx * ny +#assigning permeability for the facies code + +perm_means = np.array([0.1, 5, 3.0]) + +def assign_permeability(facies_array: np.array, perm_means: np.array, nx: int, ny: int, rng: np.random.Generator) -> np.array: + ''' + This function receives an array of any shape and converts the facies code to permeability using different + means for each facies provided in the perm_means vector. The final permeability map is a combination of the permeabilities assigned to each facies. + ''' + # Initialize the final permeability map with zeros, ensure it is float + permeability_map = np.zeros_like(facies_array, dtype=float) + + # Iterate over each facies code and assign permeability + for facies_code, mean in enumerate(perm_means): + # Create a RandomPermeability instance for the current facies + RP = resmda.RandomPermeability(nx, ny, mean, perm_min, perm_max, dtype='float64') + + + # Generate permeability for the current facies + facies_perm = RP(n=int(facies_array.shape[0]), random=rng) + + # Overlay the permeability of the current facies onto the final permeability map + mask = facies_array == facies_code + permeability_map[mask] = facies_perm[mask] + + return permeability_map + +# Now generate the permeability maps for the reference case and the prior ensemble +perm_true = assign_permeability(fac_array, perm_means, nx, ny, rng) +perm_prior = assign_permeability(fac_prior_array, perm_means, nx, ny, rng) + + + +#%% + + +# ESMDA parameters +ne = perm_prior.shape[0] # 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 +#%% +# Assumed sandard deviation of our data +dstd = 0.5 + +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] +]) + + +# TODO: change scale in imshow to represent meters + +# QC reference model and first two random models +fig, axs = plt.subplots(1, 3, constrained_layout=True) +axs[0].set_title('Model') +im = axs[0].imshow(perm_true.T, vmin=perm_min, vmax=perm_max) +axs[1].set_title('Random Model 1') +axs[1].imshow(perm_prior[0, ...].T, vmin=perm_min, vmax=perm_max) +axs[2].set_title('Random Model 2') +axs[2].imshow(perm_prior[1, ...].T, vmin=perm_min, vmax=perm_max) +fig.colorbar(im, ax=axs, orientation='horizontal', + label='Log of Permeability (mD)') +for ax in axs: + ax.set_xlabel('x-direction') +axs[0].set_ylabel('y-direction') +fig.show() +#%% + + +############################################################################### +# 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] + + +############################################################################### +# 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 +RP = resmda.RandomPermeability(nx, ny, 1, perm_min, perm_max) #had to include it here just for the localization +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): + """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, +} + + +############################################################################### +# 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 f296d27884f724a82fde82b8554fb2ec7a670879 Mon Sep 17 00:00:00 2001 From: gabrielserrao Date: Wed, 12 Jun 2024 17:50:37 +0200 Subject: [PATCH 02/13] adding post perm map with and without loc --- examples/fluvialExample.py | 66 ++++++++++++++++++++++++++++++++------ 1 file changed, 57 insertions(+), 9 deletions(-) diff --git a/examples/fluvialExample.py b/examples/fluvialExample.py index 19918a8..baded7e 100644 --- a/examples/fluvialExample.py +++ b/examples/fluvialExample.py @@ -35,10 +35,10 @@ Later, we can distribute the permeability across the reservoir using the facies code and the variance approach as in the basis case ''' -facies_ensemble = xr.open_dataset('realizations.nc') - +facies_ensemble = xr.open_dataset('/samoa/data/smrserraoseabr/GenerateModels/realizations.nc') +ne =100 #number of prior models reference_case = facies_ensemble.isel(Realization=0) -prior_ensemble = facies_ensemble.isel(Realization=slice(1,1000)) +prior_ensemble = facies_ensemble.isel(Realization=slice(1,ne+1)) fac_array = reference_case['facies code'].values.astype(float) fac_prior_array = prior_ensemble['facies code'].values.astype(float) @@ -82,7 +82,15 @@ def assign_permeability(facies_array: np.array, perm_means: np.array, nx: int, n perm_true = assign_permeability(fac_array, perm_means, nx, ny, rng) perm_prior = assign_permeability(fac_prior_array, perm_means, nx, ny, rng) - +#%% +#plot perm_true +plt.imshow(perm_true.T) +plt.colorbar(label='Log of Permeability (mD)') +plt.xlabel('x-direction') +plt.ylabel('y-direction') +plt.title('True Permeability Map') +plt.savefig(f'True_Permeability_Map_{ne}_Ensembles.png') +plt.show() #%% @@ -126,6 +134,7 @@ def assign_permeability(facies_array: np.array, perm_means: np.array, nx: int, n for ax in axs: ax.set_xlabel('x-direction') axs[0].set_ylabel('y-direction') +fig.savefig(f'QC_Models_{ne}_Ensembles.png') fig.show() #%% @@ -167,7 +176,7 @@ def sim(x): 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)]) - +fig.savefig(f'Localization_Matrix_{ne}_Ensembles.png') ############################################################################### # ES-MDA @@ -227,7 +236,7 @@ def restrict_permeability(x): for ax in axs.ravel(): for well in wells: ax.plot(well[0], well[1], ['C3v', 'C1^'][int(well[2] == 120)]) - +fig.savefig(f'Compare_Permeabilities_{ne}_Ensembles.png') ############################################################################### # Compare data @@ -248,10 +257,49 @@ def restrict_permeability(x): ax.plot(time, data_obs[0, i::3], 'C3o') ax.set_xlabel('Time') ax.set_ylabel('Pressure') - +fig.savefig(f'Compare_Data_{ne}_Ensembles.png') +plt.show() ############################################################################### -resmda.Report() - # %% +#Plot posterior models individually + +# Create separate figures for with and without localization +fig_nl, axs_nl = plt.subplots(3, 4, figsize=(12, 9), constrained_layout=True) +fig_nl.suptitle('Without Localization') +fig_wl, axs_wl = plt.subplots(3, 4, figsize=(12, 9), constrained_layout=True) +fig_wl.suptitle('With Localization') +model_indices = [0, 1, 2, 3] # Indices of the models to display + +for i, model_idx in enumerate(model_indices): + # Without localization + im_nl_prior = axs_nl[0, i].imshow(perm_prior[model_idx, ...].T) + axs_nl[0, i].set_title(f'Prior Model {model_idx+1}') + + im_nl_post = axs_nl[1, i].imshow(nl_perm_post[model_idx, ...].T) + axs_nl[1, i].set_title(f'Post Model {model_idx+1}') + + diff_nl = nl_perm_post[model_idx, ...] - perm_prior[model_idx, ...] + im_nl_diff = axs_nl[2, i].imshow(diff_nl.T) + axs_nl[2, i].set_title(f'Difference Model {model_idx+1}') + + # With localization + im_wl_prior = axs_wl[0, i].imshow(perm_prior[model_idx, ...].T) + axs_wl[0, i].set_title(f'Prior Model {model_idx+1}') + + im_wl_post = axs_wl[1, i].imshow(wl_perm_post[model_idx, ...].T) + axs_wl[1, i].set_title(f'Post Model {model_idx+1}') + + diff_wl = wl_perm_post[model_idx, ...] - perm_prior[model_idx, ...] + im_wl_diff = axs_wl[2, i].imshow(diff_wl.T) + axs_wl[2, i].set_title(f'Difference Model {model_idx+1}') + +# Add colorbars to the last column of each row for both figures +for row in range(3): + fig_nl.colorbar(axs_nl[row, -1].get_images()[0], ax=axs_nl[row, :], location='right', aspect=20) + fig_wl.colorbar(axs_wl[row, -1].get_images()[0], ax=axs_wl[row, :], location='right', aspect=20) +fig_nl.savefig(f'Posterior_Models_Without_Localization_{ne}_Ensembles.png') +fig_wl.savefig(f'Posterior_Models_With_Localization_{ne}_Ensembles.png') +# %% +resmda.Report() \ No newline at end of file From cd91143a758fef10492336e9b789aab1fbdc92f7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Tue, 18 Jun 2024 11:48:21 +0200 Subject: [PATCH 03/13] Spelling etc (#10) --- .gitignore | 2 ++ docs/manual/about.rst | 40 +++++++++++++++++++----------------- docs/manual/installation.rst | 6 ++++-- 3 files changed, 27 insertions(+), 21 deletions(-) diff --git a/.gitignore b/.gitignore index a53eba5..e4bc802 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,7 @@ # Directories and file types __pycache__/ +*.nc +*.pdf # Sphinx docs/_build/ diff --git a/docs/manual/about.rst b/docs/manual/about.rst index e501e7a..acd611d 100644 --- a/docs/manual/about.rst +++ b/docs/manual/about.rst @@ -13,9 +13,12 @@ ES-MDA In the following an introduction to the ES-MDA (Ensemble Smoother with Multiple Data Assimilation) algorithm following [EmRe13]_: -In history-matching problems, it is common to only consider the -parameter-estimation problem (neglecting model uncertainties). In that case, -the analyzed vector of model parameters :math:`m^a` is given by +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 .. math:: m_j^a = m_j^f + C_\text{MD}^f \left(C_\text{DD}^f + \alpha C_\text{D} @@ -25,24 +28,25 @@ for ensembles :math:`j=1, 2, \dots, N_e`. Here, - :math:`^a`: analysis; - :math:`^f`: forecast; -- :math:`m^f`: prior vector of model parameters; -- :math:`d^f`: vector of predicted data; +- :math:`m^f`: prior vector of model parameters (:math:`N_m`); - :math:`C_\text{MD}^f`: cross-covariance matrix between :math:`m^f` and - :math:`d^f`; -- :math:`C_\text{DD}^f`: :math:`N_d \times N_d` auto-covariance matrix of - predicted data; -- :math:`d_\text{obs}`: :math:`N_d`-dimensional vector of observed data; -- :math:`d_\text{uc} = d_\text{obs} + \sqrt{\alpha}C_\text{D}^{1/2} z_d, \ z_d - \sim \mathcal{N}(0, I_{N_d})`; -- :math:`C_\text{D}`: :math:`N_d \times N_d` covariance matrix of observed data - measurement errors; -- :math:`\alpha`: ES-MDA coefficient. + :math:`d^f` (:math:`N_m \times N_d`); +- :math:`C_\text{DD}^f`: auto-covariance matrix of predicted data + (:math:`N_d \times N_d`); +- :math:`C_\text{D}`: covariance matrix of observed data measurement errors + (:math:`N_d \times N_d`); +- :math:`\alpha`: ES-MDA coefficient; +- :math:`d_\text{uc}` : vector of perturbed data, obtained from the + vector of observed data, :math:`d_\text{obs}` (:math:`N_d`); +- :math:`d^f`: vector of predicted data (:math:`N_d`). The prior vector of model parameters, :math:`m^f_j`, can in reality be :math:`j` possible models :math:`m^f` given from an analyst (e.g., the geologist). In theoretical tests, these are usually created by perturbing the prior :math:`m^f` by, e.g., adding random Gaussian noise. +The ES-MDA algorithm follows [EmRe13]_: + 1. Choose the number of data assimilations, :math:`N_a`, and the coefficients :math:`\alpha_i` for :math:`i = 1, \dots, N_a`. 2. For :math:`i = 1` to :math:`N_a`: @@ -73,13 +77,11 @@ method. In this case, we start assimilating data with a large value of :math:`\alpha`, which reduces the magnitude of the initial updates; then, we gradually decrease :math:`\alpha`. -For ES-MDA, we only consider the parameter-estimation problem. 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. - Reservoir Model --------------- The implemented small 2D Reservoir Simulator was 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. +course material of **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. diff --git a/docs/manual/installation.rst b/docs/manual/installation.rst index 2c002e4..3d5e40a 100644 --- a/docs/manual/installation.rst +++ b/docs/manual/installation.rst @@ -1,16 +1,18 @@ Installation ============ -You can install the latest release of 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 +or clone the repository and install it manually with .. code-block:: console + git clone git@github.com:tuda-geo/resmda + cd resmda make install to get the latest version. From 2989fe417e38ff1971181132973207055da0378a Mon Sep 17 00:00:00 2001 From: gabrielserrao Date: Wed, 12 Jun 2024 16:04:27 +0200 Subject: [PATCH 04/13] adding post perm plot --- examples/fluvialExample.py | 257 +++++++++++++++++++++++++++++++++++++ 1 file changed, 257 insertions(+) create mode 100644 examples/fluvialExample.py diff --git a/examples/fluvialExample.py b/examples/fluvialExample.py new file mode 100644 index 0000000..19918a8 --- /dev/null +++ b/examples/fluvialExample.py @@ -0,0 +1,257 @@ +#%% +r""" +2D Fluvial Reservoir ESMDA example +========================== + +Ensemble Smoother Multiple Data Assimilation (ES-MDA) in Reservoir Simulation. + +""" +import numpy as np +import matplotlib.pyplot as plt +import xarray as xr + +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 + +############################################################################### +# Model parameters +# ---------------- +# Permeabilities +#read xarray from file +''' +xarray with 10000 realizations of facies for channelized reservoirs +all the models are 64x64x1 with 100x100x1 spacing (dx,dy,dz) which means they are pseudo 3D models + +We will read the first 1001 realizations and generate the permeability maps for the first 1000 realizations and +the first will be the reference case + +The property of interest is called facies code and we will assign a fixed permeability to each facies code + +Later, we can distribute the permeability across the reservoir using the facies code and the variance approach as in the basis case + +''' +facies_ensemble = xr.open_dataset('realizations.nc') + +reference_case = facies_ensemble.isel(Realization=0) +prior_ensemble = facies_ensemble.isel(Realization=slice(1,1000)) + +fac_array = reference_case['facies code'].values.astype(float) +fac_prior_array = prior_ensemble['facies code'].values.astype(float) +fac_prior_array = np.squeeze(fac_prior_array) #remove an extra dimension + +#assign bounds for the permeability +perm_min = 0.1 +perm_max =5 +# Grid extension +nx = fac_array.shape[1] +ny = fac_array.shape[2] +nc = nx * ny +#assigning permeability for the facies code + +perm_means = np.array([0.1, 5, 3.0]) + +def assign_permeability(facies_array: np.array, perm_means: np.array, nx: int, ny: int, rng: np.random.Generator) -> np.array: + ''' + This function receives an array of any shape and converts the facies code to permeability using different + means for each facies provided in the perm_means vector. The final permeability map is a combination of the permeabilities assigned to each facies. + ''' + # Initialize the final permeability map with zeros, ensure it is float + permeability_map = np.zeros_like(facies_array, dtype=float) + + # Iterate over each facies code and assign permeability + for facies_code, mean in enumerate(perm_means): + # Create a RandomPermeability instance for the current facies + RP = resmda.RandomPermeability(nx, ny, mean, perm_min, perm_max, dtype='float64') + + + # Generate permeability for the current facies + facies_perm = RP(n=int(facies_array.shape[0]), random=rng) + + # Overlay the permeability of the current facies onto the final permeability map + mask = facies_array == facies_code + permeability_map[mask] = facies_perm[mask] + + return permeability_map + +# Now generate the permeability maps for the reference case and the prior ensemble +perm_true = assign_permeability(fac_array, perm_means, nx, ny, rng) +perm_prior = assign_permeability(fac_prior_array, perm_means, nx, ny, rng) + + + +#%% + + +# ESMDA parameters +ne = perm_prior.shape[0] # 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 +#%% +# Assumed sandard deviation of our data +dstd = 0.5 + +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] +]) + + +# TODO: change scale in imshow to represent meters + +# QC reference model and first two random models +fig, axs = plt.subplots(1, 3, constrained_layout=True) +axs[0].set_title('Model') +im = axs[0].imshow(perm_true.T, vmin=perm_min, vmax=perm_max) +axs[1].set_title('Random Model 1') +axs[1].imshow(perm_prior[0, ...].T, vmin=perm_min, vmax=perm_max) +axs[2].set_title('Random Model 2') +axs[2].imshow(perm_prior[1, ...].T, vmin=perm_min, vmax=perm_max) +fig.colorbar(im, ax=axs, orientation='horizontal', + label='Log of Permeability (mD)') +for ax in axs: + ax.set_xlabel('x-direction') +axs[0].set_ylabel('y-direction') +fig.show() +#%% + + +############################################################################### +# 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] + + +############################################################################### +# 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 +RP = resmda.RandomPermeability(nx, ny, 1, perm_min, perm_max) #had to include it here just for the localization +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): + """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, +} + + +############################################################################### +# 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 64b9a14be214903a2307c271518551d837bcdca2 Mon Sep 17 00:00:00 2001 From: gabrielserrao Date: Wed, 12 Jun 2024 17:50:37 +0200 Subject: [PATCH 05/13] adding post perm map with and without loc --- examples/fluvialExample.py | 66 ++++++++++++++++++++++++++++++++------ 1 file changed, 57 insertions(+), 9 deletions(-) diff --git a/examples/fluvialExample.py b/examples/fluvialExample.py index 19918a8..baded7e 100644 --- a/examples/fluvialExample.py +++ b/examples/fluvialExample.py @@ -35,10 +35,10 @@ Later, we can distribute the permeability across the reservoir using the facies code and the variance approach as in the basis case ''' -facies_ensemble = xr.open_dataset('realizations.nc') - +facies_ensemble = xr.open_dataset('/samoa/data/smrserraoseabr/GenerateModels/realizations.nc') +ne =100 #number of prior models reference_case = facies_ensemble.isel(Realization=0) -prior_ensemble = facies_ensemble.isel(Realization=slice(1,1000)) +prior_ensemble = facies_ensemble.isel(Realization=slice(1,ne+1)) fac_array = reference_case['facies code'].values.astype(float) fac_prior_array = prior_ensemble['facies code'].values.astype(float) @@ -82,7 +82,15 @@ def assign_permeability(facies_array: np.array, perm_means: np.array, nx: int, n perm_true = assign_permeability(fac_array, perm_means, nx, ny, rng) perm_prior = assign_permeability(fac_prior_array, perm_means, nx, ny, rng) - +#%% +#plot perm_true +plt.imshow(perm_true.T) +plt.colorbar(label='Log of Permeability (mD)') +plt.xlabel('x-direction') +plt.ylabel('y-direction') +plt.title('True Permeability Map') +plt.savefig(f'True_Permeability_Map_{ne}_Ensembles.png') +plt.show() #%% @@ -126,6 +134,7 @@ def assign_permeability(facies_array: np.array, perm_means: np.array, nx: int, n for ax in axs: ax.set_xlabel('x-direction') axs[0].set_ylabel('y-direction') +fig.savefig(f'QC_Models_{ne}_Ensembles.png') fig.show() #%% @@ -167,7 +176,7 @@ def sim(x): 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)]) - +fig.savefig(f'Localization_Matrix_{ne}_Ensembles.png') ############################################################################### # ES-MDA @@ -227,7 +236,7 @@ def restrict_permeability(x): for ax in axs.ravel(): for well in wells: ax.plot(well[0], well[1], ['C3v', 'C1^'][int(well[2] == 120)]) - +fig.savefig(f'Compare_Permeabilities_{ne}_Ensembles.png') ############################################################################### # Compare data @@ -248,10 +257,49 @@ def restrict_permeability(x): ax.plot(time, data_obs[0, i::3], 'C3o') ax.set_xlabel('Time') ax.set_ylabel('Pressure') - +fig.savefig(f'Compare_Data_{ne}_Ensembles.png') +plt.show() ############################################################################### -resmda.Report() - # %% +#Plot posterior models individually + +# Create separate figures for with and without localization +fig_nl, axs_nl = plt.subplots(3, 4, figsize=(12, 9), constrained_layout=True) +fig_nl.suptitle('Without Localization') +fig_wl, axs_wl = plt.subplots(3, 4, figsize=(12, 9), constrained_layout=True) +fig_wl.suptitle('With Localization') +model_indices = [0, 1, 2, 3] # Indices of the models to display + +for i, model_idx in enumerate(model_indices): + # Without localization + im_nl_prior = axs_nl[0, i].imshow(perm_prior[model_idx, ...].T) + axs_nl[0, i].set_title(f'Prior Model {model_idx+1}') + + im_nl_post = axs_nl[1, i].imshow(nl_perm_post[model_idx, ...].T) + axs_nl[1, i].set_title(f'Post Model {model_idx+1}') + + diff_nl = nl_perm_post[model_idx, ...] - perm_prior[model_idx, ...] + im_nl_diff = axs_nl[2, i].imshow(diff_nl.T) + axs_nl[2, i].set_title(f'Difference Model {model_idx+1}') + + # With localization + im_wl_prior = axs_wl[0, i].imshow(perm_prior[model_idx, ...].T) + axs_wl[0, i].set_title(f'Prior Model {model_idx+1}') + + im_wl_post = axs_wl[1, i].imshow(wl_perm_post[model_idx, ...].T) + axs_wl[1, i].set_title(f'Post Model {model_idx+1}') + + diff_wl = wl_perm_post[model_idx, ...] - perm_prior[model_idx, ...] + im_wl_diff = axs_wl[2, i].imshow(diff_wl.T) + axs_wl[2, i].set_title(f'Difference Model {model_idx+1}') + +# Add colorbars to the last column of each row for both figures +for row in range(3): + fig_nl.colorbar(axs_nl[row, -1].get_images()[0], ax=axs_nl[row, :], location='right', aspect=20) + fig_wl.colorbar(axs_wl[row, -1].get_images()[0], ax=axs_wl[row, :], location='right', aspect=20) +fig_nl.savefig(f'Posterior_Models_Without_Localization_{ne}_Ensembles.png') +fig_wl.savefig(f'Posterior_Models_With_Localization_{ne}_Ensembles.png') +# %% +resmda.Report() \ No newline at end of file From 7d6f066e026031ec238e04ccd65ced8969d0f521 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Tue, 18 Jun 2024 17:58:44 +0200 Subject: [PATCH 06/13] Working on fluvial example --- .gitignore | 1 + examples/fluvialExample.py | 305 ----------------------------- examples/fluvialreservoir.py | 360 +++++++++++++++++++++++++++++++++++ requirements-dev.txt | 1 + 4 files changed, 362 insertions(+), 305 deletions(-) delete mode 100644 examples/fluvialExample.py create mode 100644 examples/fluvialreservoir.py diff --git a/.gitignore b/.gitignore index e4bc802..71761dc 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,7 @@ __pycache__/ *.nc *.pdf +*.zip # Sphinx docs/_build/ diff --git a/examples/fluvialExample.py b/examples/fluvialExample.py deleted file mode 100644 index baded7e..0000000 --- a/examples/fluvialExample.py +++ /dev/null @@ -1,305 +0,0 @@ -#%% -r""" -2D Fluvial Reservoir ESMDA example -========================== - -Ensemble Smoother Multiple Data Assimilation (ES-MDA) in Reservoir Simulation. - -""" -import numpy as np -import matplotlib.pyplot as plt -import xarray as xr - -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 - -############################################################################### -# Model parameters -# ---------------- -# Permeabilities -#read xarray from file -''' -xarray with 10000 realizations of facies for channelized reservoirs -all the models are 64x64x1 with 100x100x1 spacing (dx,dy,dz) which means they are pseudo 3D models - -We will read the first 1001 realizations and generate the permeability maps for the first 1000 realizations and -the first will be the reference case - -The property of interest is called facies code and we will assign a fixed permeability to each facies code - -Later, we can distribute the permeability across the reservoir using the facies code and the variance approach as in the basis case - -''' -facies_ensemble = xr.open_dataset('/samoa/data/smrserraoseabr/GenerateModels/realizations.nc') -ne =100 #number of prior models -reference_case = facies_ensemble.isel(Realization=0) -prior_ensemble = facies_ensemble.isel(Realization=slice(1,ne+1)) - -fac_array = reference_case['facies code'].values.astype(float) -fac_prior_array = prior_ensemble['facies code'].values.astype(float) -fac_prior_array = np.squeeze(fac_prior_array) #remove an extra dimension - -#assign bounds for the permeability -perm_min = 0.1 -perm_max =5 -# Grid extension -nx = fac_array.shape[1] -ny = fac_array.shape[2] -nc = nx * ny -#assigning permeability for the facies code - -perm_means = np.array([0.1, 5, 3.0]) - -def assign_permeability(facies_array: np.array, perm_means: np.array, nx: int, ny: int, rng: np.random.Generator) -> np.array: - ''' - This function receives an array of any shape and converts the facies code to permeability using different - means for each facies provided in the perm_means vector. The final permeability map is a combination of the permeabilities assigned to each facies. - ''' - # Initialize the final permeability map with zeros, ensure it is float - permeability_map = np.zeros_like(facies_array, dtype=float) - - # Iterate over each facies code and assign permeability - for facies_code, mean in enumerate(perm_means): - # Create a RandomPermeability instance for the current facies - RP = resmda.RandomPermeability(nx, ny, mean, perm_min, perm_max, dtype='float64') - - - # Generate permeability for the current facies - facies_perm = RP(n=int(facies_array.shape[0]), random=rng) - - # Overlay the permeability of the current facies onto the final permeability map - mask = facies_array == facies_code - permeability_map[mask] = facies_perm[mask] - - return permeability_map - -# Now generate the permeability maps for the reference case and the prior ensemble -perm_true = assign_permeability(fac_array, perm_means, nx, ny, rng) -perm_prior = assign_permeability(fac_prior_array, perm_means, nx, ny, rng) - -#%% -#plot perm_true -plt.imshow(perm_true.T) -plt.colorbar(label='Log of Permeability (mD)') -plt.xlabel('x-direction') -plt.ylabel('y-direction') -plt.title('True Permeability Map') -plt.savefig(f'True_Permeability_Map_{ne}_Ensembles.png') -plt.show() - -#%% - - -# ESMDA parameters -ne = perm_prior.shape[0] # 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 -#%% -# Assumed sandard deviation of our data -dstd = 0.5 - -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] -]) - - -# TODO: change scale in imshow to represent meters - -# QC reference model and first two random models -fig, axs = plt.subplots(1, 3, constrained_layout=True) -axs[0].set_title('Model') -im = axs[0].imshow(perm_true.T, vmin=perm_min, vmax=perm_max) -axs[1].set_title('Random Model 1') -axs[1].imshow(perm_prior[0, ...].T, vmin=perm_min, vmax=perm_max) -axs[2].set_title('Random Model 2') -axs[2].imshow(perm_prior[1, ...].T, vmin=perm_min, vmax=perm_max) -fig.colorbar(im, ax=axs, orientation='horizontal', - label='Log of Permeability (mD)') -for ax in axs: - ax.set_xlabel('x-direction') -axs[0].set_ylabel('y-direction') -fig.savefig(f'QC_Models_{ne}_Ensembles.png') -fig.show() -#%% - - -############################################################################### -# 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] - - -############################################################################### -# 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 -RP = resmda.RandomPermeability(nx, ny, 1, perm_min, perm_max) #had to include it here just for the localization -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)]) -fig.savefig(f'Localization_Matrix_{ne}_Ensembles.png') - -############################################################################### -# ES-MDA -# ------ - - -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, -} - - -############################################################################### -# 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)]) -fig.savefig(f'Compare_Permeabilities_{ne}_Ensembles.png') - -############################################################################### -# 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') -fig.savefig(f'Compare_Data_{ne}_Ensembles.png') -plt.show() - -############################################################################### - -# %% -#Plot posterior models individually - -# Create separate figures for with and without localization -fig_nl, axs_nl = plt.subplots(3, 4, figsize=(12, 9), constrained_layout=True) -fig_nl.suptitle('Without Localization') -fig_wl, axs_wl = plt.subplots(3, 4, figsize=(12, 9), constrained_layout=True) -fig_wl.suptitle('With Localization') -model_indices = [0, 1, 2, 3] # Indices of the models to display - -for i, model_idx in enumerate(model_indices): - # Without localization - im_nl_prior = axs_nl[0, i].imshow(perm_prior[model_idx, ...].T) - axs_nl[0, i].set_title(f'Prior Model {model_idx+1}') - - im_nl_post = axs_nl[1, i].imshow(nl_perm_post[model_idx, ...].T) - axs_nl[1, i].set_title(f'Post Model {model_idx+1}') - - diff_nl = nl_perm_post[model_idx, ...] - perm_prior[model_idx, ...] - im_nl_diff = axs_nl[2, i].imshow(diff_nl.T) - axs_nl[2, i].set_title(f'Difference Model {model_idx+1}') - - # With localization - im_wl_prior = axs_wl[0, i].imshow(perm_prior[model_idx, ...].T) - axs_wl[0, i].set_title(f'Prior Model {model_idx+1}') - - im_wl_post = axs_wl[1, i].imshow(wl_perm_post[model_idx, ...].T) - axs_wl[1, i].set_title(f'Post Model {model_idx+1}') - - diff_wl = wl_perm_post[model_idx, ...] - perm_prior[model_idx, ...] - im_wl_diff = axs_wl[2, i].imshow(diff_wl.T) - axs_wl[2, i].set_title(f'Difference Model {model_idx+1}') - -# Add colorbars to the last column of each row for both figures -for row in range(3): - fig_nl.colorbar(axs_nl[row, -1].get_images()[0], ax=axs_nl[row, :], location='right', aspect=20) - fig_wl.colorbar(axs_wl[row, -1].get_images()[0], ax=axs_wl[row, :], location='right', aspect=20) -fig_nl.savefig(f'Posterior_Models_Without_Localization_{ne}_Ensembles.png') -fig_wl.savefig(f'Posterior_Models_With_Localization_{ne}_Ensembles.png') -# %% -resmda.Report() \ No newline at end of file diff --git a/examples/fluvialreservoir.py b/examples/fluvialreservoir.py new file mode 100644 index 0000000..d01ee2b --- /dev/null +++ b/examples/fluvialreservoir.py @@ -0,0 +1,360 @@ +r""" +2D Fluvial Reservoir ESMDA example +================================== + +In contrast to the basic reservoir example +:ref:`sphx_glr_gallery_basicreservoir.py`, where a single facies was used, this +example uses fluvial models. It also makes use of localization as explained in +: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. + + +""" +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) + +# Adjust this path to a folder of your choice. +data_path = os.path.join('..', 'download', '') + +# sphinx_gallery_thumbnail_number = 3 + +############################################################################### +# Load and plot the facies +# ------------------------ + +fname = 'facies.npy' +pooch.retrieve( + 'https://raw.github.com/emsig/data/2021-05-21/emg3d/surveys/'+fname, + 'hashTODO', + fname=fname, + path=data_path, +) +facies = np.load(data_path + fname) +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() +code_to_perm = [str(i)+" = "+str(p) for i, p in enumerate(perm_means)] +fig.suptitle(f'Facies: {code_to_perm}') +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') + + +############################################################################### +# Assign random permeabilities to the facies +# ------------------------------------------ + +perm_min = 0.05 +perm_max = 10.0 + +# Instantiate a random permeability instance +RP = resmda.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) + + +############################################################################### +# Model parameters +# ---------------- + +# We take the first model as "true/reference", and the other for ES-MDA. +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 + +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] +]) + + +############################################################################### +# 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] + + +############################################################################### +# 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)) + + +############################################################################### +# ES-MDA +# ------ + + +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, +} + + +############################################################################### +# 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)]) +fig.savefig(f'Compare_Permeabilities_{ne}_Ensembles.png') + + +############################################################################### +# 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') +fig.savefig(f'Compare_Data_{ne}_Ensembles.png') +fig.show() + +############################################################################### +# Reproduce the facies +# -------------------- +# +# .. note:: +# +# The following sections are about how to reproduce the facies data +# ``facies.npz``. 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 +# open-source Fortran library FLUVSIM: +# +# **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 +# `_. +# +# +# Required imports +# ---------------- +# +# .. code-block:: python +# +# 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 +# +# +# .. code-block:: python +# +# # 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. +# """ +# 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 +# ----------------- +# +# .. code-block:: python +# +# # 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 output +# --------------- +# +# .. code-block:: python +# +# with open('facies.json', "w") as f: +# json.dump(all_params, f, indent=2) +# np.save('facies', facies.squeeze(), allow_pickle=False) + + +############################################################################### +resmda.Report() diff --git a/requirements-dev.txt b/requirements-dev.txt index 230467f..84611cb 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -16,6 +16,7 @@ sphinx_gallery pydata_sphinx_theme sphinx_automodapi ipykernel +pooch # FOR TESTING flake8 From 2a2c49f25cf11cc6582933df9686bf5f4de5e2ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Tue, 18 Jun 2024 20:58:11 +0200 Subject: [PATCH 07/13] Make data download work --- .gitignore | 1 + examples/fluvialreservoir.py | 85 ++++++++++++++++++++---------------- 2 files changed, 48 insertions(+), 38 deletions(-) diff --git a/.gitignore b/.gitignore index 71761dc..2ff0282 100644 --- a/.gitignore +++ b/.gitignore @@ -9,6 +9,7 @@ docs/_build/ docs/api/resmda* docs/savefig/ docs/gallery/* +download/ # Pytest and coverage related htmlcov diff --git a/examples/fluvialreservoir.py b/examples/fluvialreservoir.py index d01ee2b..3e1d951 100644 --- a/examples/fluvialreservoir.py +++ b/examples/fluvialreservoir.py @@ -13,6 +13,12 @@ """ +import os +import json + +# To retrieve the data, you need to have pooch installed: +# ``pip install pooch`` or ``conda install -c conda-forge pooch`` +import pooch import numpy as np import matplotlib.pyplot as plt @@ -33,8 +39,8 @@ fname = 'facies.npy' pooch.retrieve( - 'https://raw.github.com/emsig/data/2021-05-21/emg3d/surveys/'+fname, - 'hashTODO', + 'https://raw.github.com/tuda-geo/data/2024-06-18/resmda/'+fname, + '4bfe56c836bf17ca63453c37e5da91cb97bbef8cc6c08d605f70bd64fe7488b2', fname=fname, path=data_path, ) @@ -204,7 +210,6 @@ def restrict_permeability(x): for ax in axs.ravel(): for well in wells: ax.plot(well[0], well[1], ['C3v', 'C1^'][int(well[2] == 120)]) -fig.savefig(f'Compare_Permeabilities_{ne}_Ensembles.png') ############################################################################### @@ -226,7 +231,6 @@ def restrict_permeability(x): ax.plot(time, data_obs[0, i::3], 'C3o') ax.set_xlabel('Time') ax.set_ylabel('Pressure') -fig.savefig(f'Compare_Data_{ne}_Ensembles.png') fig.show() ############################################################################### @@ -235,26 +239,24 @@ def restrict_permeability(x): # # .. note:: # -# The following sections are about how to reproduce the facies data -# ``facies.npz``. 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 -# open-source Fortran library FLUVSIM: +# The following cell is about how to reproduce the facies data loaded in +# ``facies.npy``. This was created using ``geomodpy``. # -# **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. +# ``geomodpy`` (Guillaume Rongier, 2023) is not open-source yet. The +# functionality of geomodpy that we use here is a python wrapper for the +# open-source Fortran library ``FLUVSIM``: # -# DOI: `10.1016/S0098-3004(01)00075-9 -# `_. +# **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 +# `_. # -# Required imports -# ---------------- # # .. code-block:: python # +# # ==== Required imports ==== # import json # import numpy as np # @@ -264,16 +266,12 @@ def restrict_permeability(x): # # 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 # # -# .. code-block:: python +# # ==== 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 = { @@ -296,7 +294,7 @@ def restrict_permeability(x): # """Generate geological parameters from normal distributions. # # Expects for each parameter a tuple of two values, or a list -# containing tuples of two values. +# containing tuples of two values each. # """ # geol_params = {} # for param, dist in geol_dists.items(): @@ -307,13 +305,9 @@ def restrict_permeability(x): # else: # geol_params[param] = rng.normal(*dist) # return geol_params - - -############################################################################### -# Create the facies -# ----------------- # -# .. code-block:: python +# +# # ==== Create the facies ==== # # # Number of sets and realizations # nsets = 2 @@ -343,18 +337,33 @@ def restrict_permeability(x): # # realizations = fluvsim.run().data_vars['facies code'].values # facies[i*nreal:(i+1)*nreal, ...] = realizations.astype('i4') - - -############################################################################### -# Save the output -# --------------- # -# .. code-block:: python # +# # ==== Save the output ==== +# +# # 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', facies.squeeze(), allow_pickle=False) +############################################################################### +# 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. + +fname = 'facies.json' +pooch.retrieve( + 'https://raw.github.com/tuda-geo/data/2024-06-18/resmda/'+fname, + 'db2cb8a620775c68374c24a4fa811f6350381c7fc98a823b9571136d307540b4', + fname=fname, + path=data_path, +) +with open(data_path + fname, 'r') as f: + print(json.dumps(json.load(f), indent=2)) + ############################################################################### resmda.Report() From 2a34ae1d60297a3107e8fdba4c923311ca3a9a08 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Wed, 19 Jun 2024 12:26:18 +0200 Subject: [PATCH 08/13] small edits --- examples/basicESMDA.py | 1 + examples/fluvialreservoir.py | 21 ++++++++++++++++++--- examples/localization.py | 3 +++ 3 files changed, 22 insertions(+), 3 deletions(-) diff --git a/examples/basicESMDA.py b/examples/basicESMDA.py index f7e525a..43c8839 100644 --- a/examples/basicESMDA.py +++ b/examples/basicESMDA.py @@ -57,6 +57,7 @@ def forward(x, beta): axs[i].plot(px, forward(px, b)) axs[i].set_xlabel('x') axs[i].set_ylabel('y') +fig.show() ############################################################################### diff --git a/examples/fluvialreservoir.py b/examples/fluvialreservoir.py index 3e1d951..acbbcb8 100644 --- a/examples/fluvialreservoir.py +++ b/examples/fluvialreservoir.py @@ -4,20 +4,32 @@ In contrast to the basic reservoir example :ref:`sphx_glr_gallery_basicreservoir.py`, where a single facies was used, this -example uses fluvial models. It also makes use of localization as explained in +example uses fluvial models containing different facies. It also compares the +use of ES-MDA with and without localization, as explained in :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 """ import os import json -# To retrieve the data, you need to have pooch installed: -# ``pip install pooch`` or ``conda install -c conda-forge pooch`` import pooch import numpy as np import matplotlib.pyplot as plt @@ -62,6 +74,7 @@ clim=[-0.5, 2.5], origin='lower' ) fig.colorbar(im, ax=axs, ticks=[0, 1, 2], label='Facies code') +fig.show() ############################################################################### @@ -89,6 +102,7 @@ for i in range(ne): im = axs[i].imshow(permeabilities[i, ...], origin='lower') fig.colorbar(im, ax=axs) +fig.show() ############################################################################### @@ -210,6 +224,7 @@ def restrict_permeability(x): for ax in axs.ravel(): for well in wells: ax.plot(well[0], well[1], ['C3v', 'C1^'][int(well[2] == 120)]) +fig.show() ############################################################################### diff --git a/examples/localization.py b/examples/localization.py index 8d42bf7..1cb321e 100644 --- a/examples/localization.py +++ b/examples/localization.py @@ -104,6 +104,7 @@ def sim(x): 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)]) +fig.show() ############################################################################### @@ -164,6 +165,7 @@ def restrict_permeability(x): for ax in axs.ravel(): for well in wells: ax.plot(well[0], well[1], ['C3v', 'C1^'][int(well[2] == 120)]) +fig.show() ############################################################################### @@ -185,6 +187,7 @@ def restrict_permeability(x): ax.plot(time, data_obs[0, i::3], 'C3o') ax.set_xlabel('Time') ax.set_ylabel('Pressure') +fig.show() ############################################################################### From a50c89d45368c84a3709496ad4ca1e23f6cf2132 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Wed, 19 Jun 2024 17:35:21 +0200 Subject: [PATCH 09/13] WIP; need to check: lables, c-scales, axis --- examples/basicreservoir.py | 26 ++++++++++++++------------ examples/fluvialreservoir.py | 10 +++++----- examples/localization.py | 20 ++++++++++++++------ 3 files changed, 33 insertions(+), 23 deletions(-) diff --git a/examples/basicreservoir.py b/examples/basicreservoir.py index 767052a..ba654f5 100644 --- a/examples/basicreservoir.py +++ b/examples/basicreservoir.py @@ -21,7 +21,7 @@ # ---------------- # Grid extension -nx = 25 +nx = 30 ny = 25 nc = nx*ny @@ -61,23 +61,22 @@ 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) +pinp = {'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, vmin=perm_min, vmax=perm_max) +im = axs[0, 0].imshow(perm_true.T, **pinp) 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, 0].imshow(perm_prior[0, ...].T, **pinp) axs[1, 1].set_title('Random Model 2') -axs[1, 1].imshow(perm_prior[1, ...].T, vmin=perm_min, vmax=perm_max) +axs[1, 1].imshow(perm_prior[1, ...].T, **pinp) fig.colorbar(im, ax=axs[1, :], orientation='horizontal', label='Log of Permeability (mD)') -for ax in axs[1, :]: +for ax in axs[1, :].ravel(): ax.set_xlabel('x-direction') -for ax in axs[:, 0]: +for ax in axs[:, 0].ravel(): ax.set_ylabel('y-direction') fig.show() @@ -142,14 +141,17 @@ def restrict_permeability(x): # models. Here, we visualize the first three realizations from both the prior # and posterior ensembles to see how the models have been updated. +# TODO: change from mean to a particular example + # 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) +im = ax[0].imshow(perm_prior.mean(axis=0).T, **pinp) ax[1].set_title('Post Mean') -ax[1].imshow(perm_post.mean(axis=0).T) +ax[1].imshow(perm_post.mean(axis=0).T, **pinp) ax[2].set_title('"Truth"') -ax[2].imshow(perm_true.T) +ax[2].imshow(perm_true.T, **pinp) +fig.colorbar(im, ax=ax[2], label='Log of Permeability (mD)') fig.show() diff --git a/examples/fluvialreservoir.py b/examples/fluvialreservoir.py index acbbcb8..e983ad5 100644 --- a/examples/fluvialreservoir.py +++ b/examples/fluvialreservoir.py @@ -210,15 +210,15 @@ def restrict_permeability(x): 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, 0].imshow(perm_prior.mean(axis=0).T, vmin=perm_min, vmax=perm_max) axs[0, 1].set_title('"Truth"') -axs[0, 1].imshow(perm_true.T) +axs[0, 1].imshow(perm_true.T, vmin=perm_min, vmax=perm_max) axs[1, 0].set_title('Post Mean without localization') -axs[1, 0].imshow(nl_perm_post.mean(axis=0).T) +axs[1, 0].imshow(nl_perm_post.mean(axis=0).T, vmin=perm_min, vmax=perm_max) axs[1, 1].set_title('Post Mean with localization') -axs[1, 1].imshow(wl_perm_post.mean(axis=0).T) +axs[1, 1].imshow(wl_perm_post.mean(axis=0).T, vmin=perm_min, vmax=perm_max) axs[1, 1].contour(loc_mat.sum(axis=2).T, levels=[2.0, ], colors='w') for ax in axs.ravel(): @@ -259,7 +259,7 @@ def restrict_permeability(x): # # ``geomodpy`` (Guillaume Rongier, 2023) is not open-source yet. The # functionality of geomodpy that we use here is a python wrapper for the -# open-source Fortran library ``FLUVSIM``: +# 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: diff --git a/examples/localization.py b/examples/localization.py index 1cb321e..7750052 100644 --- a/examples/localization.py +++ b/examples/localization.py @@ -22,7 +22,7 @@ # ---------------- # Grid extension -nx = 30 +nx = 35 ny = 30 nc = nx*ny @@ -100,8 +100,10 @@ def sim(x): # 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.imshow(loc_mat.sum(axis=2).T, origin='lower') ax.contour(loc_mat.sum(axis=2).T, levels=[2.0, ], colors='w') +ax.set_xlabel('x-direction') +ax.set_ylabel('y-direction') for well in wells: ax.plot(well[0], well[1], ['C3v', 'C1^'][int(well[2] == 120)]) fig.show() @@ -151,20 +153,26 @@ def restrict_permeability(x): 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) +im = axs[0, 0].imshow(perm_prior.mean(axis=0).T, origin='lower') axs[0, 1].set_title('"Truth"') -axs[0, 1].imshow(perm_true.T) +axs[0, 1].imshow(perm_true.T, origin='lower') axs[1, 0].set_title('Post Mean without localization') -axs[1, 0].imshow(nl_perm_post.mean(axis=0).T) +axs[1, 0].imshow(nl_perm_post.mean(axis=0).T, origin='lower') axs[1, 1].set_title('Post Mean with localization') -axs[1, 1].imshow(wl_perm_post.mean(axis=0).T) +axs[1, 1].imshow(wl_perm_post.mean(axis=0).T, origin='lower') axs[1, 1].contour(loc_mat.sum(axis=2).T, levels=[2.0, ], colors='w') +fig.colorbar(im, ax=axs, orientation='horizontal', + label='Log of Permeability (mD)') for ax in axs.ravel(): for well in wells: ax.plot(well[0], well[1], ['C3v', 'C1^'][int(well[2] == 120)]) +for ax in axs[1, :].ravel(): + ax.set_xlabel('x-direction') +for ax in axs[:, 0].ravel(): + ax.set_ylabel('y-direction') fig.show() From d77829d77436a2c6368bce603e82cc74235c750a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Mon, 1 Jul 2024 17:32:57 +0200 Subject: [PATCH 10/13] working on docs --- .gitignore | 1 + docs/conf.py | 3 +-- requirements-dev.txt | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.gitignore b/.gitignore index 2ff0282..9de7d08 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ docs/api/resmda* docs/savefig/ docs/gallery/* download/ +docs/sg_execution_times.rst # Pytest and coverage related htmlcov diff --git a/docs/conf.py b/docs/conf.py index eb39077..20d0e12 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,7 +1,6 @@ import time import warnings from resmda import __version__ -from sphinx_gallery.sorting import FileNameSortKey # ==== 1. Extensions ==== @@ -43,7 +42,7 @@ # 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 diff --git a/requirements-dev.txt b/requirements-dev.txt index 84611cb..0b5e3aa 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -8,11 +8,11 @@ matplotlib setuptools_scm # FOR DOCUMENTATION -sphinx +sphinx>=7.3 numpydoc sphinx_design sphinx_numfig -sphinx_gallery +sphinx_gallery>=0.16 pydata_sphinx_theme sphinx_automodapi ipykernel From 886a44a27d2411558d9fe55e684bf995d5a4d1ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Thu, 11 Jul 2024 12:01:34 +0200 Subject: [PATCH 11/13] Clean up fluvial example --- docs/conf.py | 2 +- examples/fluvialreservoir.py | 132 ++++++++++++++++++----------------- 2 files changed, 68 insertions(+), 66 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 20d0e12..9267b86 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -38,7 +38,7 @@ 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) diff --git a/examples/fluvialreservoir.py b/examples/fluvialreservoir.py index e983ad5..b525a21 100644 --- a/examples/fluvialreservoir.py +++ b/examples/fluvialreservoir.py @@ -5,7 +5,7 @@ In contrast to the basic reservoir example :ref:`sphx_glr_gallery_basicreservoir.py`, where a single facies was used, this example uses fluvial models containing different facies. It also compares the -use of ES-MDA with and without localization, as explained in +use of ES-MDA 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 @@ -38,10 +38,10 @@ # For reproducibility, we instantiate a random number generator with a fixed # seed. For production, remove the seed! -rng = np.random.default_rng(1848) +rng = np.random.default_rng(1324) # Adjust this path to a folder of your choice. -data_path = os.path.join('..', 'download', '') +data_path = os.path.join("..", "download", "") # sphinx_gallery_thumbnail_number = 3 @@ -49,10 +49,10 @@ # Load and plot the facies # ------------------------ -fname = 'facies.npy' +fname = "facies.npy" pooch.retrieve( - 'https://raw.github.com/tuda-geo/data/2024-06-18/resmda/'+fname, - '4bfe56c836bf17ca63453c37e5da91cb97bbef8cc6c08d605f70bd64fe7488b2', + "https://raw.github.com/tuda-geo/data/2024-06-18/resmda/"+fname, + "4bfe56c836bf17ca63453c37e5da91cb97bbef8cc6c08d605f70bd64fe7488b2", fname=fname, path=data_path, ) @@ -63,17 +63,16 @@ 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) +fig, axs = plt.subplots( + 2, 5, figsize=(8, 3), sharex=True, sharey=True, constrained_layout=True) axs = axs.ravel() -code_to_perm = [str(i)+" = "+str(p) for i, p in enumerate(perm_means)] -fig.suptitle(f'Facies: {code_to_perm}') +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' + clim=[-0.5, 2.5], origin="lower" ) -fig.colorbar(im, ax=axs, ticks=[0, 1, 2], label='Facies code') +fig.colorbar(im, ax=axs, ticks=[0, 1, 2], label="Facies code") fig.show() @@ -95,12 +94,12 @@ 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) +fig, axs = plt.subplots( + 2, 5, figsize=(8, 3), sharex=True, sharey=True, constrained_layout=True) axs = axs.ravel() -fig.suptitle('Permeabilities') +fig.suptitle("Permeabilities") for i in range(ne): - im = axs[i].imshow(permeabilities[i, ...], origin='lower') + im = axs[i].imshow(permeabilities[i, ...], origin="lower") fig.colorbar(im, ax=axs) fig.show() @@ -121,6 +120,7 @@ # Assumed standard deviation of our data dstd = 0.5 +# Measurement points ox = (5, 15, 24) oy = (5, 10, 24) @@ -129,9 +129,9 @@ # Wells wells = np.array([ - [5, 5, 180], [5, 12, 120], - [15, 10, 180], [20, 5, 120], - [24, 24, 180], [24, 17, 120] + [ox[0], oy[0], 180], [5, 12, 120], + [ox[1], oy[1], 180], [20, 5, 120], + [ox[2], oy[2], 180], [24, 17, 120] ]) @@ -177,14 +177,14 @@ def restrict_permeability(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, + "model_prior": perm_prior, + "forward": sim, + "data_obs": data_obs, + "sigma": dstd, + "alphas": 4, + "data_prior": data_prior, + "callback_post": restrict_permeability, + "random": rng, } @@ -209,22 +209,21 @@ def restrict_permeability(x): # 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].set_title("Prior Mean") axs[0, 0].imshow(perm_prior.mean(axis=0).T, vmin=perm_min, vmax=perm_max) -axs[0, 1].set_title('"Truth"') +axs[0, 1].set_title("'Truth'") axs[0, 1].imshow(perm_true.T, vmin=perm_min, vmax=perm_max) -axs[1, 0].set_title('Post Mean without localization') +axs[1, 0].set_title("Post Mean without localization") axs[1, 0].imshow(nl_perm_post.mean(axis=0).T, vmin=perm_min, vmax=perm_max) -axs[1, 1].set_title('Post Mean with localization') +axs[1, 1].set_title("Post Mean with localization") axs[1, 1].imshow(wl_perm_post.mean(axis=0).T, vmin=perm_min, vmax=perm_max) -axs[1, 1].contour(loc_mat.sum(axis=2).T, levels=[2.0, ], colors='w') +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)]) -fig.show() + ax.plot(well[0], well[1], ["C3v", "C1^"][int(well[2] == 120)]) ############################################################################### @@ -234,19 +233,22 @@ def restrict_permeability(x): # 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') +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') + 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, 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') -fig.show() + 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") + ############################################################################### # Reproduce the facies @@ -290,18 +292,18 @@ def restrict_permeability(x): # # # 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': [ +# "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), +# "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), # } # # @@ -334,13 +336,13 @@ def restrict_permeability(x): # # Pre-allocate containers to store all realizations and their # # corresponding parameters # all_params = {} -# facies = np.zeros((nsets * nreal, nz, nx, ny), dtype='i4') +# 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 +# all_params[f"set-{i}"] = params # # fluvsim = FLUVSIM( # shape=(nx, ny, nz), @@ -350,17 +352,17 @@ def restrict_permeability(x): # **params # ) # -# realizations = fluvsim.run().data_vars['facies code'].values -# facies[i*nreal:(i+1)*nreal, ...] = realizations.astype('i4') +# realizations = fluvsim.run().data_vars["facies code"].values +# facies[i*nreal:(i+1)*nreal, ...] = realizations.astype("i4") # # # # ==== Save the output ==== # # # Save the input parameters to FLUVSIM as a json. -# with open('facies.json', "w") as f: +# 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', facies.squeeze(), allow_pickle=False) +# np.save("facies", facies.squeeze(), allow_pickle=False) ############################################################################### @@ -370,14 +372,14 @@ def restrict_permeability(x): # These are, just as the data themselves, online at # https://github.com/tuda-geo/data/resmda. -fname = 'facies.json' +fname = "facies.json" pooch.retrieve( - 'https://raw.github.com/tuda-geo/data/2024-06-18/resmda/'+fname, - 'db2cb8a620775c68374c24a4fa811f6350381c7fc98a823b9571136d307540b4', + "https://raw.github.com/tuda-geo/data/2024-06-18/resmda/"+fname, + "db2cb8a620775c68374c24a4fa811f6350381c7fc98a823b9571136d307540b4", fname=fname, path=data_path, ) -with open(data_path + fname, 'r') as f: +with open(data_path + fname, "r") as f: print(json.dumps(json.load(f), indent=2)) ############################################################################### From 7664cb11ef558be190b6a260435e699723946561 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Thu, 11 Jul 2024 12:59:34 +0200 Subject: [PATCH 12/13] Clean up fluvial example II --- examples/fluvialreservoir.py | 34 +++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/examples/fluvialreservoir.py b/examples/fluvialreservoir.py index b525a21..fbbbdf1 100644 --- a/examples/fluvialreservoir.py +++ b/examples/fluvialreservoir.py @@ -38,7 +38,7 @@ # For reproducibility, we instantiate a random number generator with a fixed # seed. For production, remove the seed! -rng = np.random.default_rng(1324) +rng = np.random.default_rng(1513) # Adjust this path to a folder of your choice. data_path = os.path.join("..", "download", "") @@ -73,7 +73,6 @@ clim=[-0.5, 2.5], origin="lower" ) fig.colorbar(im, ax=axs, ticks=[0, 1, 2], label="Facies code") -fig.show() ############################################################################### @@ -100,8 +99,7 @@ fig.suptitle("Permeabilities") for i in range(ne): im = axs[i].imshow(permeabilities[i, ...], origin="lower") -fig.colorbar(im, ax=axs) -fig.show() +fig.colorbar(im, ax=axs, label="Log10 Permeability (mD)") ############################################################################### @@ -208,20 +206,25 @@ def restrict_permeability(x): # 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, vmin=perm_min, vmax=perm_max) -axs[0, 1].set_title("'Truth'") -axs[0, 1].imshow(perm_true.T, vmin=perm_min, vmax=perm_max) + 1, 3, figsize=(8, 4), sharex=True, sharey=True, constrained_layout=True) +par = {"vmin": perm_min, "vmax": perm_max, "origin": "lower"} -axs[1, 0].set_title("Post Mean without localization") -axs[1, 0].imshow(nl_perm_post.mean(axis=0).T, vmin=perm_min, vmax=perm_max) -axs[1, 1].set_title("Post Mean with localization") -axs[1, 1].imshow(wl_perm_post.mean(axis=0).T, vmin=perm_min, vmax=perm_max) -axs[1, 1].contour(loc_mat.sum(axis=2).T, levels=[2.0, ], colors="w") +axs[0].set_title("Prior Mean") +im = axs[0].imshow(perm_prior.mean(axis=0).T, **par) -for ax in axs.ravel(): + +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="Log10 Permeabilities (mD)", + orientation="horizontal") + +for ax in axs: for well in wells: ax.plot(well[0], well[1], ["C3v", "C1^"][int(well[2] == 120)]) @@ -235,6 +238,7 @@ def restrict_permeability(x): 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") From 85ff9f2bfa33b0f70274264de410e8f1ee6b2ab1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Thu, 11 Jul 2024 13:44:36 +0200 Subject: [PATCH 13/13] Clean up fluvial example III --- docs/conf.py | 2 +- examples/basicESMDA.py | 8 ++-- examples/basicreservoir.py | 79 ++++++++++++++++------------------ examples/fluvialreservoir.py | 7 ++- examples/localization.py | 82 +++++++++++++++++++----------------- 5 files changed, 88 insertions(+), 90 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 9267b86..e74345b 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -48,7 +48,7 @@ # 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 diff --git a/examples/basicESMDA.py b/examples/basicESMDA.py index 43c8839..5562435 100644 --- a/examples/basicESMDA.py +++ b/examples/basicESMDA.py @@ -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 @@ -53,11 +53,10 @@ 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') -fig.show() ############################################################################### @@ -110,7 +109,6 @@ def plot_result(mpost, dpost, dobs, title, ylim): ax2.set_xlabel('y') ax2.set_xlim([-5, 8]) ax2.legend() - fig.show() ############################################################################### diff --git a/examples/basicreservoir.py b/examples/basicreservoir.py index ba654f5..d786718 100644 --- a/examples/basicreservoir.py +++ b/examples/basicreservoir.py @@ -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 @@ -62,23 +61,22 @@ # QC covariance, reference model, and first two random models -pinp = {'origin': 'lower', 'vmin': perm_min, 'vmax': perm_max} +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, **pinp) -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, **pinp) -axs[1, 1].set_title('Random Model 2') -axs[1, 1].imshow(perm_prior[1, ...].T, **pinp) -fig.colorbar(im, ax=axs[1, :], orientation='horizontal', - label='Log of Permeability (mD)') +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') + ax.set_xlabel("x-direction") for ax in axs[:, 0].ravel(): - ax.set_ylabel('y-direction') -fig.show() + ax.set_ylabel("y-direction") ############################################################################### @@ -100,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)") ############################################################################### @@ -141,18 +138,15 @@ def restrict_permeability(x): # models. Here, we visualize the first three realizations from both the prior # and posterior ensembles to see how the models have been updated. -# TODO: change from mean to a particular example - # Plot posterior -fig, ax = plt.subplots(1, 3, figsize=(12, 5)) -ax[0].set_title('Prior Mean') -im = ax[0].imshow(perm_prior.mean(axis=0).T, **pinp) -ax[1].set_title('Post Mean') -ax[1].imshow(perm_post.mean(axis=0).T, **pinp) -ax[2].set_title('"Truth"') -ax[2].imshow(perm_true.T, **pinp) -fig.colorbar(im, ax=ax[2], label='Log of Permeability (mD)') -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") ############################################################################### @@ -164,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)") ############################################################################### diff --git a/examples/fluvialreservoir.py b/examples/fluvialreservoir.py index fbbbdf1..14fb3ef 100644 --- a/examples/fluvialreservoir.py +++ b/examples/fluvialreservoir.py @@ -99,7 +99,7 @@ fig.suptitle("Permeabilities") for i in range(ne): im = axs[i].imshow(permeabilities[i, ...], origin="lower") -fig.colorbar(im, ax=axs, label="Log10 Permeability (mD)") +fig.colorbar(im, ax=axs, label="Log Permeability (mD)") ############################################################################### @@ -221,12 +221,15 @@ def restrict_permeability(x): 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="Log10 Permeabilities (mD)", +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') ############################################################################### diff --git a/examples/localization.py b/examples/localization.py index 7750052..fa6ed6f 100644 --- a/examples/localization.py +++ b/examples/localization.py @@ -39,7 +39,7 @@ # Assumed sandard deviation of our data dstd = 0.5 -# Observation location indices (should be well locations) +# Measurement points ox = (5, 15, 24) oy = (5, 10, 24) @@ -48,9 +48,9 @@ # Wells wells = np.array([ - [5, 5, 180], [5, 12, 120], - [15, 10, 180], [20, 5, 120], - [24, 24, 180], [24, 17, 120] + [ox[0], oy[0], 180], [5, 12, 120], + [ox[1], oy[1], 180], [20, 5, 120], + [ox[2], oy[2], 180], [24, 17, 120] ]) @@ -99,14 +99,14 @@ def sim(x): 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, origin='lower') -ax.contour(loc_mat.sum(axis=2).T, levels=[2.0, ], colors='w') +fig, ax = plt.subplots(1, 1, constrained_layout=True) +im = ax.imshow(loc_mat.sum(axis=2).T/time.size, origin='lower', cmap="plasma") +ax.contour(loc_mat.sum(axis=2).T/time.size, levels=[0.2, ], colors='w') ax.set_xlabel('x-direction') ax.set_ylabel('y-direction') for well in wells: ax.plot(well[0], well[1], ['C3v', 'C1^'][int(well[2] == 120)]) -fig.show() +fig.colorbar(im, ax=ax, label="Localization Weight (-)") ############################################################################### @@ -151,29 +151,30 @@ def restrict_permeability(x): # 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') -im = axs[0, 0].imshow(perm_prior.mean(axis=0).T, origin='lower') -axs[0, 1].set_title('"Truth"') -axs[0, 1].imshow(perm_true.T, origin='lower') - - -axs[1, 0].set_title('Post Mean without localization') -axs[1, 0].imshow(nl_perm_post.mean(axis=0).T, origin='lower') -axs[1, 1].set_title('Post Mean with localization') -axs[1, 1].imshow(wl_perm_post.mean(axis=0).T, origin='lower') -axs[1, 1].contour(loc_mat.sum(axis=2).T, levels=[2.0, ], colors='w') -fig.colorbar(im, ax=axs, orientation='horizontal', - label='Log of Permeability (mD)') - -for ax in axs.ravel(): + 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[1, :].ravel(): + ax.plot(well[0], well[1], ["C3v", "C1^"][int(well[2] == 120)]) +for ax in axs: ax.set_xlabel('x-direction') -for ax in axs[:, 0].ravel(): - ax.set_ylabel('y-direction') -fig.show() +axs[0].set_ylabel('y-direction') ############################################################################### @@ -185,17 +186,20 @@ def restrict_permeability(x): 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') + 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, 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') -fig.show() + 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") ###############################################################################