diff --git a/ensemble/ensemble.py b/ensemble/ensemble.py index 982f78c..48cd716 100644 --- a/ensemble/ensemble.py +++ b/ensemble/ensemble.py @@ -15,6 +15,7 @@ from tqdm.auto import tqdm from p_tqdm import p_map import logging +from geostat.decomp import Cholesky # Making realizations # Internal imports import pipt.misc_tools.analysis_tools as at @@ -26,7 +27,6 @@ from misc.system_tools.environ_var import OpenBlasSingleThread # Single threaded OpenBLAS runs - class Ensemble: """ Class for organizing misc. variables and simulator for an ensemble-based inversion run. Here, the forecast step @@ -139,15 +139,24 @@ def __init__(self, keys_en, sim, redund_sim=None): # individually). self.state = {key: val for key, val in tmp_load.items()} - # Find the number of ensemble members from state variable + # Find the number of ensemble members from loaded state variables tmp_ne = [] for tmp_state in self.state.keys(): tmp_ne.extend([self.state[tmp_state].shape[1]]) - if max(tmp_ne) != min(tmp_ne): - print('\033[1;33mInput states have different ensemble size\033[1;m') - sys.exit(1) - self.ne = min(tmp_ne) - + + if 'ne' not in self.keys_en: # NE not specified in input file + if max(tmp_ne) != min(tmp_ne): #Check loaded ensembles are the same size (if more than one state variable) + print('\033[1;33mInput states have different ensemble size\033[1;m') + sys.exit(1) + self.ne = min(tmp_ne) # Use the number of ensemble members in loaded ensemble + else: + # Use the number of ensemble members specified in input file (may be fewer than loaded) + self.ne = int(self.keys_en['ne']) + if self.ne <= min(tmp_ne): + # pick correct number of ensemble members + self.state = {key: val[:,:self.ne] for key, val in self.state.items()} + else: + print('\033[1;33mInput states are smaller than NE\033[1;m') if 'multilevel' in self.keys_en: ml_info = extract.extract_multilevel_info(self.keys_en) self.multilevel, self.tot_level, self.ml_ne, self.ML_error_corr, self.error_comp_scheme, self.ML_corr_done = ml_info @@ -338,6 +347,20 @@ def calc_prediction(self, input_state=None, save_prediction=None): # Index list of ensemble members list_member_index = list(range(self.ne)) + # modified by xluo, for including the simulation of the mean reservoir model + # as used in the RLM-MAC algorithm + if 'daalg' in self.keys_en and self.keys_en['daalg'][1] == 'gies': + list_state.append({}) + list_member_index.append(self.ne) + + for key in self.state.keys(): + tmp_state = np.zeros(list_state[0][key].shape[0]) + + for i in range(self.ne): + tmp_state += list_state[i][key] + + list_state[self.ne][key] = tmp_state / self.ne + if no_tot_run==1: # if not in parallel we use regular loop en_pred = [self.sim.run_fwd_sim(state, member_index) for state, member_index in tqdm(zip(list_state, list_member_index), total=len(list_state))] @@ -392,6 +415,7 @@ def calc_prediction(self, input_state=None, save_prediction=None): else: # Run prediction in parallel using p_map en_pred = p_map(self.sim.run_fwd_sim, list_state, list_member_index, num_cpus=no_tot_run, disable=self.disable_tqdm) + # List successful runs and crashes list_crash = [indx for indx, el in enumerate(en_pred) if el is False] list_success = [indx for indx, el in enumerate(en_pred) if el is not False] diff --git a/pipt/loop/assimilation.py b/pipt/loop/assimilation.py index 6803d27..be9272b 100644 --- a/pipt/loop/assimilation.py +++ b/pipt/loop/assimilation.py @@ -529,8 +529,12 @@ def post_process_forecast(self): for k in pred_data_tmp[i]: # DATATYPE if vintage < len(self.ensemble.sparse_info['mask']) and \ len(pred_data_tmp[i][k]) == int(np.sum(self.ensemble.sparse_info['mask'][vintage])): - self.ensemble.pred_data[i][k] = np.zeros( - (len(self.ensemble.obs_data[i][k]), self.ensemble.ne)) + if self.ensemble.keys_da['daalg'][1] == 'gies': + self.ensemble.pred_data[i][k] = np.zeros( + (len(self.ensemble.obs_data[i][k]), self.ensemble.ne+1)) + else: + self.ensemble.pred_data[i][k] = np.zeros( + (len(self.ensemble.obs_data[i][k]), self.ensemble.ne)) for m in range(pred_data_tmp[i][k].shape[1]): data_array = self.ensemble.compress(pred_data_tmp[i][k][:, m], vintage, self.ensemble.sparse_info['use_ensemble']) diff --git a/pipt/loop/ensemble.py b/pipt/loop/ensemble.py index 5a44ba3..b986a76 100644 --- a/pipt/loop/ensemble.py +++ b/pipt/loop/ensemble.py @@ -540,7 +540,24 @@ def _ext_scaling(self): self.state_scaling = at.calc_scaling( self.prior_state, self.list_states, self.prior_info) - self.Am = None + delta_scaled_prior = self.state_scaling[:, None] * \ + np.dot(at.aug_state(self.prior_state, self.list_states), self.proj) + + u_d, s_d, v_d = np.linalg.svd(delta_scaled_prior, full_matrices=False) + + # remove the last singular value/vector. This is because numpy returns all ne values, while the last is actually + # zero. This part is a good place to include eventual additional truncation. + energy = 0 + trunc_index = len(s_d) - 1 # inititallize + for c, elem in enumerate(s_d): + energy += elem + if energy / sum(s_d) >= self.trunc_energy: + trunc_index = c # take the index where all energy is preserved + break + u_d, s_d, v_d = u_d[:, :trunc_index + + 1], s_d[:trunc_index + 1], v_d[:trunc_index + 1, :] + self.Am = np.dot(u_d, np.eye(trunc_index+1) * + ((s_d**(-1))[:, None])) # notation from paper def save_temp_state_assim(self, ind_save): """ @@ -694,7 +711,7 @@ def compress(self, data=None, vintage=0, aug_coeff=None): data_array = None - elif aug_coeff is None: + elif aug_coeff is None: # compress predicted data data_array, wdec_rec = self.sparse_data[vintage].compress(data) rec = self.sparse_data[vintage].reconstruct( @@ -703,7 +720,7 @@ def compress(self, data=None, vintage=0, aug_coeff=None): self.data_rec.append([]) self.data_rec[vintage].append(rec) - elif not aug_coeff: + elif not aug_coeff: # compress true data, aug_coeff = false options = copy(self.sparse_info) # find the correct mask for the vintage diff --git a/pipt/misc_tools/analysis_tools.py b/pipt/misc_tools/analysis_tools.py index 59f748f..de0140a 100644 --- a/pipt/misc_tools/analysis_tools.py +++ b/pipt/misc_tools/analysis_tools.py @@ -542,8 +542,9 @@ def calc_objectivefun(pert_obs, pred_data, Cd): data_misfit : array-like Nex1 array containing objective function values. """ - ne = pred_data.shape[1] - r = (pred_data - pert_obs) + #ne = pred_data.shape[1] + ne = pert_obs.shape[1] + r = (pred_data[:, :ne] - pert_obs) # This is necessary to use to gies code that xilu has implemented if len(Cd.shape) == 1: precission = Cd**(-1) data_misfit = np.diag(r.T.dot(r*precission[:, None])) diff --git a/pipt/misc_tools/extract_tools.py b/pipt/misc_tools/extract_tools.py index 8f25616..c59d2d5 100644 --- a/pipt/misc_tools/extract_tools.py +++ b/pipt/misc_tools/extract_tools.py @@ -271,8 +271,12 @@ def organize_sparse_representation(info: Union[dict,list]) -> dict: sparse['dim'] = [dim[2], dim[1], dim[0]] # Read mask_files - sparse['mask'] = [] - for idx, filename in enumerate(info['mask'], start=1): + sparse['mask'] = [] + m_info = info['mask'] + # allow for one mask with filename given as string + if isinstance(m_info, str): + m_info = [m_info] + for idx, filename in enumerate(m_info, start=1): if not os.path.exists(filename): mask = np.ones(sparse['dim'], dtype=bool) np.savez(f'mask_{idx}.npz', mask=mask) diff --git a/pipt/misc_tools/qaqc_tools.py b/pipt/misc_tools/qaqc_tools.py index e97f86c..85e2001 100644 --- a/pipt/misc_tools/qaqc_tools.py +++ b/pipt/misc_tools/qaqc_tools.py @@ -73,8 +73,10 @@ def __init__(self, keys, obs_data, datavar, logger=None, prior_info=None, sim=No assim_step = 0 # Assume simultaneous assimiation assim_ind = [keys['obsname'], keys['assimindex'][assim_step]] + #assim_ind = [keys['obsname'], keys['assimindex']] if isinstance(assim_ind[1], list): # Check if prim. ind. is a list self.l_prim = [int(x) for x in assim_ind[1]] + #self.l_prim = [int(x[0]) for x in assim_ind[1]] else: # Float self.l_prim = [int(assim_ind[1])] @@ -133,16 +135,19 @@ def set(self, pred_data, state=None, lam=None): else: self.en_fcst[typ] = np.array( [self.pred_data[ind][typ].flatten() for ind in self.l_prim if + self.obs_data[ind][typ] is not None and sum(np.isnan(self.obs_data[ind][typ])) == 0 and self.obs_data[ind][typ].shape == (1,)]) - l = [self.pred_data[ind][typ] for ind in self.l_prim if sum(np.isnan(self.obs_data[ind][typ])) == 0 + l = [self.pred_data[ind][typ] for ind in self.l_prim if + self.obs_data[ind][typ] is not None + and sum(np.isnan(self.obs_data[ind][typ])) == 0 and self.obs_data[ind][typ].shape[0] > 1] if l: self.en_fcst_vec[typ] = np.concatenate(l) self.state = state self.lam = lam - def calc_coverage(self, line=None, field_dim=None): + def calc_coverage(self, line=None, field_dim=None, uxl = None, uil = None, contours = None, uxl_c = None, uil_c = None): """ Calculate the Data coverage for production and seismic data. For seismic data the plotting is based on the importance-scaled coverage developed by Espen O. Lie from GeoCore. @@ -280,7 +285,7 @@ def _plot_coverage_1D(line, field_dim): plt.savefig(filename) os.system('convert ' + filename + '.png' + ' -trim ' + filename + '.png') - for typ in [dat for dat in self.data_types if not dat in ['bulkimp', 'sim2seis']]: # Only well data + for typ in [dat for dat in self.data_types if not dat in ['bulkimp', 'sim2seis', 'avo', 'grav']]: # Only well data if hasattr(self, 'multilevel'): # calc for each level plt.figure() cover_low = [True for _ in self.en_obs[typ]] @@ -330,12 +335,13 @@ def _plot_coverage_1D(line, field_dim): # Plot the seismic data data_sim = [] data = [] - supported_data = ['sim2seis', 'bulkimp'] + supported_data = ['sim2seis', 'bulkimp', 'avo', 'grav'] my_data = [dat for dat in supported_data if dat in self.data_types] if len(my_data) == 0: return else: my_data = my_data[0] + #my_data = my_data[1] # get the data seis_scaling = 1.0 @@ -391,15 +397,23 @@ def _plot_coverage_1D(line, field_dim): rgb.append(f(attr)) rgb = np.dstack(rgb) - try: - uxl = loadmat('seglines.mat')['uxl'].flatten() - uil = loadmat('seglines.mat')['uil'].flatten() - except: - uxl = [0, field_dim[0]] - uil = [0, field_dim[1]] + if uxl is None and uil is None: + try: + uxl = loadmat('seglines.mat')['uxl'].flatten() + uil = loadmat('seglines.mat')['uil'].flatten() + except: + uxl = [0, field_dim[0]] + uil = [0, field_dim[1]] + extent = (uxl[0], uxl[-1], uil[-1], uil[0]) plt.figure() plt.imshow(rgb, extent=extent) + if contours is not None and uil_c is not None and uxl_c is not None: + plt.contour(uxl_c, uil_c, contours[::-1, :], levels=1, colors='black') + plt.xlim(uxl[0], uxl[-1]) + plt.ylim(uil[-1], uil[0]) + plt.xlabel('Easting (km)') + plt.ylabel('Northing (km)') plt.title('Coverage - not scaled by Importance - epsilon=' + str(nl)) filename = self.folder + 'coverage_vint_' + str(vint) plt.savefig(filename) @@ -415,6 +429,12 @@ def _plot_coverage_1D(line, field_dim): rgb_scaled = cv2.cvtColor(hls, cv2.COLOR_HLS2RGB) rgb = rgb_scaled / 255 plt.imshow(rgb, extent=extent) + if contours is not None and uil_c is not None and uxl_c is not None: + plt.contour(uxl_c, uil_c, contours[::-1, :], levels=1, colors='black', extent=extent) + plt.xlim(uxl[0], uxl[-1]) + plt.ylim(uil[-1], uil[0]) + plt.xlabel('Easting (km)') + plt.ylabel('Northing (km)') plt.title('Coverage - scaled by Importance - epsilon=' + str(nl)) filename = self.folder + 'coverage_importance_vint_' + str(vint) plt.savefig(filename) @@ -422,7 +442,13 @@ def _plot_coverage_1D(line, field_dim): plt.close() plt.figure() - plt.imshow(sat, extent=extent) + plt.imshow(sat[::-1,:], extent=extent) + if contours is not None and uil_c is not None and uxl_c is not None: + plt.contour(uxl_c, uil_c, contours[::-1, :], levels=1, colors='black', extent=extent) + plt.xlim(uxl[0], uxl[-1]) + plt.ylim(uil[-1], uil[0]) + plt.xlabel('Easting (km)') + plt.ylabel('Northing (km)') plt.title('Importance - epsilon=' + str(nl)) filename = self.folder + 'importance_vint_' + str(vint) plt.savefig(filename) @@ -594,7 +620,8 @@ def _plot_kg(_field=None): else: idx = actnum kg = np.ma.array(data=tmp, mask=~idx) - dim = (self.prior_info[param]['nx'], self.prior_info[param]['ny'], self.prior_info[param]['nz']) + #dim = (self.prior_info[param]['nx'], self.prior_info[param]['ny'], self.prior_info[param]['nz']) + dim = next((item[1] for item in self.prior_info[param] if item[0] == 'grid'), None) input_time = None if write_to_resinsight: if time is None: @@ -643,7 +670,7 @@ def _plot_kg(_field=None): :] mean_residual = (en_obs[typ][ind] - en_ml_fcst[typ][l][ind, :]).mean() mean_residual = mean_residual[np.newaxis, np.newaxis].flatten() - delta_d = (en_obs[typ][ind] - en_ml_fcst[typ][l][ind, :])[np.newaxis, :] + delta_d = (en_obs[typ][ind] - en_ml_fcst[typ][l][ind, :self.ne])[np.newaxis, :] X2 = _calc_proj() self.state = self.ML_state[l] num_cell = self.state[param].shape[0] @@ -656,10 +683,10 @@ def _plot_kg(_field=None): self.state = copy.deepcopy(self.ML_state) delattr(self, 'ML_state') else: - pert_pred = (en_fcst[typ][ind, :] - en_fcst[typ][ind, :].mean())[np.newaxis, :] - mean_residual = (en_obs[typ][ind] - en_fcst[typ][ind, :]).mean() + pert_pred = (en_fcst[typ][ind, :self.ne] - en_fcst[typ][ind, :self.ne].mean())[np.newaxis, :] + mean_residual = (en_obs[typ][ind] - en_fcst[typ][ind, :self.ne]).mean() mean_residual = mean_residual[np.newaxis, np.newaxis].flatten() - delta_d = (en_obs[typ][ind] - en_fcst[typ][ind, :])[np.newaxis, :] + delta_d = (en_obs[typ][ind] - en_fcst[typ][ind, :self.ne])[np.newaxis, :] X2 = _calc_proj() num_cell = self.state[param].shape[0] tmp = _calc_kalman_gain() @@ -688,7 +715,7 @@ def _plot_kg(_field=None): if len(en_ml_fcst[typ][l].shape) == 2: pert_pred = en_ml_fcst[typ][l] - np.dot(en_ml_fcst[typ][l].mean(axis=1)[:, np.newaxis], np.ones((1, self.ml_ne[l]))) - delta_d = en_obs[typ] - en_ml_fcst[typ][l] + delta_d = en_obs[typ] - en_ml_fcst[typ][l][:,:self.ne] mean_residual = (en_obs[typ] - en_ml_fcst[typ][l]).mean(axis=1) X2 = _calc_proj() self.state = self.ML_state[l] @@ -709,10 +736,10 @@ def _plot_kg(_field=None): else: # combine time instances if len(en_fcst[typ].shape) == 2: - pert_pred = en_fcst[typ] - np.dot(en_fcst[typ].mean(axis=1)[:, np.newaxis], + pert_pred = en_fcst[typ][:, :self.ne] - np.dot(en_fcst[typ][:, :self.ne].mean(axis=1)[:, np.newaxis], np.ones((1, self.ne))) - delta_d = en_obs[typ] - en_fcst[typ] - mean_residual = (en_obs[typ] - en_fcst[typ]).mean(axis=1) + delta_d = en_obs[typ] - en_fcst[typ][:, :self.ne] + mean_residual = (en_obs[typ] - en_fcst[typ][:, :self.ne]).mean(axis=1) X2 = _calc_proj() for param in self.list_state: num_cell = self.state[param].shape[0] @@ -750,12 +777,12 @@ def _plot_kg(_field=None): param = el[1] time = len(self.l_prim) time_str = '-' - pert_pred = en_fcst[typ] - np.dot(en_fcst[typ].mean(axis=1)[:, np.newaxis], + pert_pred = en_fcst[typ][:, :self.ne] - np.dot(en_fcst[typ][:, :self.ne].mean(axis=1)[:, np.newaxis], np.ones((1, self.ne))) mean_residual = (en_obs[typ] - en_fcst[typ]).mean(axis=1) t_var = [self.datavar[ind][typ] for ind in en_time[typ] if self.datavar[ind][typ] is not None] X2 = _calc_proj() - delta_d = en_obs[typ] - en_fcst[typ] + delta_d = en_obs[typ] - en_fcst[typ][:, :self.ne] num_cell = self.state[param].shape[0] tmp = _calc_kalman_gain() if el_ind < len(kg_max_mean): @@ -805,7 +832,7 @@ def calc_mahalanobis(self, combi_list=(1, None)): en_fcst_pert = en_fcst + np.sqrt(en_var[:, 0])[:, np.newaxis] * \ np.random.randn(en_fcst.shape[0], en_fcst.shape[1]) - else: # some data should be defines as blocks. To get the correct measure we project the data onto the subspace + else: # some data should be defined as blocks. To get the correct measure we project the data onto the subspace # spanned by the first principal component. The level 1, 2 and 3. Difference is then calculated in # similar fashion as for the full data-space. have simple rules for generating combinations. All data are # aquired at some time, at some position, and there might be multiple data types at the same time and @@ -814,7 +841,7 @@ def calc_mahalanobis(self, combi_list=(1, None)): if 'time' in combi_type or 'vector' in combi_type: tmp_fcst = [] for typ in self.data_types: - tmp_fcst.append([self.en_fcst[typ][ind, :][np.newaxis, :] for ind in self.l_prim + tmp_fcst.append([self.en_fcst[typ][ind, :self.ne][np.newaxis, :self.ne] for ind in self.l_prim if self.obs_data[ind][typ] is not None and sum( np.isnan(self.obs_data[ind][typ])) == 0]) filt_fcst = [x for x in tmp_fcst if len(x)] # remove all empty lists @@ -1001,7 +1028,7 @@ def calc_da_stat(self, options=None): idx = np.ones(self.state[key].shape[0], dtype=bool) else: idx = actnum - if M.size == np.sum(idx): # we have a grid parameter + if M.size == np.sum(idx) and M.size > 1: # we have a grid parameter tmp = np.zeros(M.shape) tmp[M > S] = 1 tmp[M > 2 * S] = 2 @@ -1013,6 +1040,7 @@ def calc_da_stat(self, options=None): data[idx] = tmp field = np.ma.array(data=data, mask=~idx) dim = (self.prior_info[key]['nx'], self.prior_info[key]['ny'], self.prior_info[key]['nz']) + #dim = next((item[1] for item in self.prior_info[key] if item[0] == 'grid'), None) input_time = None if hasattr(self.sim, 'write_to_grid'): self.sim.write_to_grid(field, f'da_stat_{key}', self.folder, dim, input_time) diff --git a/pipt/misc_tools/wavelet_tools.py b/pipt/misc_tools/wavelet_tools.py index a907cb6..ac206de 100644 --- a/pipt/misc_tools/wavelet_tools.py +++ b/pipt/misc_tools/wavelet_tools.py @@ -16,7 +16,6 @@ class SparseRepresentation: def __init__(self, options): # options: dim, actnum, level, wname, colored_noise, threshold_rule, th_mult, # use_hard_th, keep_ca, inactive_value - # dim must be given as (nz,ny,nx) self.options = options self.num_grid = np.prod(self.options['dim']) @@ -30,7 +29,7 @@ def __init__(self, options): self.ca_leading_index = None self.ca_leading_coeff = None - # Function to doing image compression. If the function is called without threshold, then the leading indices must + # Function for image compression. If the function is called without threshold, then the leading indices must # be defined in the class. Typically, this is done by running the compression on true data with a given threshold. def compress(self, data, th_mult=None): if ('inactive_value' not in self.options) or (self.options['inactive_value'] is None): @@ -43,13 +42,10 @@ def compress(self, data, th_mult=None): if 'min_noise' not in self.options: self.options['min_noise'] = 1.0e-9 signal = signal.reshape(self.options['dim'], order=self.options['order']) - # get the signal back into its original shape (nx,ny,nz) - signal = signal.transpose((2, 1, 0)) - # pywt throws a warning in case of single-dimentional entries in the shape of the signal. - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - wdec = pywt.wavedecn( - signal, self.options['wname'], 'symmetric', int(self.options['level'])) + + # Wavelet decomposition + wdec = pywt.wavedecn(signal, self.options['wname'], 'symmetric', int(self.options['level'])) + wdec_rec = deepcopy(wdec) # Perform thresholding if the threshold is given as input. @@ -65,7 +61,12 @@ def compress(self, data, th_mult=None): # Initialize # Note: the keys below are organized the same way as in Matlab. - keys = ['daa', 'ada', 'dda', 'aad', 'dad', 'add', 'ddd'] + if signal.ndim == 3: + keys = ['daa', 'ada', 'dda', 'aad', 'dad', 'add', 'ddd'] + details = 'ddd' + elif signal.ndim == 2: + keys = ['da', 'ad', 'dd'] + details = 'dd' if true_data: for level in range(0, int(self.options['level'])+1): num_subband = 1 if level == 0 else len(keys) @@ -85,7 +86,7 @@ def compress(self, data, th_mult=None): # In the white noise case estimated std is based on the high (hhh) subband only if true_data and not self.options['colored_noise']: - subband_hhh = wdec[-1]['ddd'].flatten() + subband_hhh = wdec[-1][details].flatten() est_noise_level = np.median( np.abs(subband_hhh - np.median(subband_hhh))) / 0.6745 # estimated noise std est_noise_level = np.maximum(est_noise_level, self.options['min_noise']) @@ -224,10 +225,10 @@ def reconstruct(self, wdec_rec): print('No signal to reconstruct') sys.exit(1) + # reconstruct from wavelet coefficients data_rec = pywt.waverecn(wdec_rec, self.options['wname'], 'symmetric') - data_rec = data_rec.transpose((2, 1, 0)) # flip the axes - dim = self.options['dim'] - data_rec = data_rec[0:dim[0], 0:dim[1], 0:dim[2]] # severe issure here + data_rec = data_rec[tuple(slice(0, s) for s in self.options['dim'])] + data_rec = data_rec.flatten(order=self.options['order']) data_rec = data_rec[self.options['mask']] diff --git a/pipt/update_schemes/gies/gies_base.py b/pipt/update_schemes/gies/gies_base.py new file mode 100644 index 0000000..fc96ee5 --- /dev/null +++ b/pipt/update_schemes/gies/gies_base.py @@ -0,0 +1,284 @@ +""" +EnRML type schemes +""" +# External imports +import pipt.misc_tools.analysis_tools as at +from geostat.decomp import Cholesky +from pipt.loop.ensemble import Ensemble +from pipt.update_schemes.update_methods_ns.subspace_update import subspace_update +from pipt.update_schemes.update_methods_ns.full_update import full_update +from pipt.update_schemes.update_methods_ns.approx_update import approx_update +import sys +import pkgutil +import inspect +import numpy as np +import copy as cp +from scipy.linalg import cholesky, solve + +# Internal imports + + +class GIESMixIn(Ensemble): + """ + This is a base template for implementating the generalized iterative ensemble smoother (GIES) in the following papers: + Luo, Xiaodong. "Novel iterative ensemble smoothers derived from a class of generalized cost functions." + Computational Geosciences 25.3 (2021): 1159-1189. + Luo, Xiaodong, and William C. Cruz. "Data assimilation with soft constraints (DASC) through a generalized iterative + ensemble smoother." Computational Geosciences 26.3 (2022): 571-594. + """ + + def __init__(self, keys_da, keys_fwd, sim): + """ + The class is initialized by passing the PIPT init. file upwards in the hierarchy to be read and parsed in + `pipt.input_output.pipt_init.ReadInitFile`. + + Parameters + ---------- + init_file: str + PIPT init. file containing info. to run the inversion algorithm + """ + # Pass the init_file upwards in the hierarchy + super().__init__(keys_da, keys_fwd, sim) + + if self.restart is False: + # Save prior state in separate variable + self.prior_state = cp.deepcopy(self.state) + + # Extract parameters like conv. tol. and damping param. from ITERATION keyword in DATAASSIM + self._ext_iter_param() + + # Within variables + self.prev_data_misfit = None # Data misfit at previous iteration + if 'actnum' in self.keys_da.keys(): + try: + self.actnum = np.load(self.keys_da['actnum'])['actnum'] + except: + print('ACTNUM file cannot be loaded!') + else: + self.actnum = None + # At the moment, the iterative loop is threated as an iterative smoother and thus we check if assim. indices + # are given as in the Simultaneous loop. + self.check_assimindex_simultaneous() + # define the assimilation index + self.assim_index = [self.keys_da['obsname'], self.keys_da['assimindex'][0]] + # define the list of states + self.list_states = list(self.state.keys()) + # define the list of datatypes + self.list_datatypes, self.list_act_datatypes = at.get_list_data_types( + self.obs_data, self.assim_index) + # Get the perturbed observations and observation scaling + self.data_random_state = cp.deepcopy(np.random.get_state()) + self._ext_obs() + # Get state scaling and svd of scaled prior + self._ext_scaling() + self.current_state = cp.deepcopy(self.state) + + def calc_analysis(self): + """ + Calculate the update step in LM-EnRML, which is just the Levenberg-Marquardt update algorithm with + the sensitivity matrix approximated by the ensemble. + """ + + # reformat predicted data + _, self.aug_pred_data = at.aug_obs_pred_data(self.obs_data, self.pred_data, self.assim_index, + self.list_datatypes) + + if self.iteration == 1: # first iteration + data_misfit = at.calc_objectivefun( + self.real_obs_data, self.aug_pred_data, self.cov_data) + + # Store the (mean) data misfit (also for conv. check) + self.data_misfit = np.mean(data_misfit) + self.prior_data_misfit = np.mean(data_misfit) + self.data_misfit_std = np.std(data_misfit) + + #self.logger.info( + # f'Prior run complete with data misfit: {self.prior_data_misfit:0.1f}. Lambda for initial analysis: {self.lam}') + + if 'localanalysis' in self.keys_da: + self.local_analysis_update() + else: + # Mean pred_data and perturbation matrix with scaling + if len(self.scale_data.shape) == 1: + #self.pert_preddata = np.dot(np.expand_dims(self.scale_data ** (-1), axis=1), + # np.ones((1, self.ne))) * np.dot(self.aug_pred_data[:, 0:self.ne], self.proj) + self.pert_preddata = np.dot(np.expand_dims(self.scale_data ** (-1), axis=1), + np.ones((1, self.ne))) * (self.aug_pred_data[:, 0:self.ne] - self.aug_pred_data[:, self.ne, None]) + else: + #self.pert_preddata = solve( + # self.scale_data, np.dot(self.aug_pred_data[:, 0:self.ne], self.proj)) + self.pert_preddata = solve( + self.scale_data, self.aug_pred_data[:, 0:self.ne] - self.aug_pred_data[:, self.ne, None]) + + aug_state = at.aug_state(self.current_state, self.list_states) + self.update() # run ordinary analysis + if hasattr(self, 'step'): + aug_state_upd = aug_state + self.step + if hasattr(self, 'w_step'): + self.W = self.current_W + self.w_step + aug_prior_state = at.aug_state(self.prior_state, self.list_states) + aug_state_upd = np.dot(aug_prior_state, (np.eye( + self.ne) + self.W / np.sqrt(self.ne - 1))) + + # Extract updated state variables from aug_update + self.state = at.update_state(aug_state_upd, self.state, self.list_states) + self.state = at.limits(self.state, self.prior_info) + + def check_convergence(self): + """ + Check if LM-EnRML have converged based on evaluation of change sizes of objective function, state and damping + parameter. + + Returns + ------- + conv: bool + Logic variable telling if algorithm has converged + why_stop: dict + Dict. with keys corresponding to conv. criteria, with logical variable telling which of them that has been + met + """ + + _, pred_data = at.aug_obs_pred_data(self.obs_data, self.pred_data, self.assim_index, + self.list_datatypes) + # Initialize the initial success value + success = False + + # if inital conv. check, there are no prev_data_misfit + self.prev_data_misfit = self.data_misfit + self.prev_data_misfit_std = self.data_misfit_std + + # Calc. std dev of data misfit (used to update lamda) + # mat_obs = np.dot(obs_data_vector.reshape((len(obs_data_vector),1)), np.ones((1, self.ne))) # use the perturbed + # data instead. + + data_misfit = at.calc_objectivefun(self.real_obs_data, pred_data, self.cov_data) + + self.data_misfit = np.mean(data_misfit) + self.data_misfit_std = np.std(data_misfit) + + # # Calc. mean data misfit for convergence check, using the updated state variable + # self.data_misfit = np.dot((mean_preddata - obs_data_vector).T, + # solve(cov_data, (mean_preddata - obs_data_vector))) + + # Convergence check: Relative step size of data misfit or state change less than tolerance + if abs(1 - (self.data_misfit / self.prev_data_misfit)) < self.data_misfit_tol: + #or self.lam >= self.lam_max: + # Logical variables for conv. criteria + why_stop = {'data_misfit_stop': 1 - (self.data_misfit / self.prev_data_misfit) < self.data_misfit_tol, + 'data_misfit': self.data_misfit, + 'prev_data_misfit': self.prev_data_misfit, + 'lambda': self.lam} + if hasattr(self, 'lam_max'): + why_stop['lambda_stop'] = (self.lam >= self.lam_max) + + if self.data_misfit >= self.prev_data_misfit: + success = False + self.logger.info( + f'Iterations have converged after {self.iteration} iterations. Objective function reduced ' + f'from {self.prior_data_misfit:0.2f} to {self.prev_data_misfit:0.2f}') + else: + self.logger.info( + f'Iterations have converged after {self.iteration} iterations. Objective function reduced ' + f'from {self.prior_data_misfit:0.2f} to {self.data_misfit:0.2f}') + # Return conv = True, why_stop var. + return True, success, why_stop + + else: # conv. not met + # Logical variables for conv. criteria + why_stop = {'data_misfit_stop': 1 - (self.data_misfit / self.prev_data_misfit) < self.data_misfit_tol, + 'data_misfit': self.data_misfit, + 'prev_data_misfit': self.prev_data_misfit, + 'lambda': self.lam} + if hasattr(self, 'lam_max'): + why_stop['lambda_stop'] = (self.lam >= self.lam_max) + + ############################################### + ##### update Lambda step-size values ########## + ############################################### + # If reduction in mean data misfit, reduce damping param + if self.data_misfit < self.prev_data_misfit: + # Reduce damping parameter (divide calculations for ANALYSISDEBUG purpose) + if not hasattr(self, 'lam_min'): + self.lam = self.lam / self.gamma + else: + if self.lam > self.lam_min: + self.lam = self.lam / self.gamma + + success = True + self.current_state = cp.deepcopy(self.state) + if hasattr(self, 'W'): + self.current_W = cp.deepcopy(self.W) + + else: # Reject iteration, and increase lam + # Increase damping parameter (divide calculations for ANALYSISDEBUG purpose) + self.lam = self.lam * self.gamma + success = False + + self.logger.info(f'Iter {self.iteration}: ') + if success: + #self.logger.info(f'Successfull iteration number {self.iteration}! Objective function reduced from ' + # f'{self.prev_data_misfit:0.1f} to {self.data_misfit:0.1f}. New Lamba for next analysis: ' + # f'{self.lam}') + pass + + # self.prev_data_misfit = self.data_misfit + # self.prev_data_misfit_std = self.data_misfit_std + else: + #self.logger.info(f'Failed iteration number {self.iteration}! Objective function increased from ' + # f'{self.prev_data_misfit:0.1f} to {self.data_misfit:0.1f}. New Lamba for repeated analysis: ' + # f'{self.lam}') + # Reset the objective function after report + self.data_misfit = self.prev_data_misfit + self.data_misfit_std = self.prev_data_misfit_std + + # Return conv = False, why_stop var. + return False, success, why_stop + + def _ext_iter_param(self): + """ + Extract parameters needed in LM-EnRML from the ITERATION keyword given in the DATAASSIM part of PIPT init. + file. These parameters include convergence tolerances and parameters for the damping parameter. Default + values for these parameters have been given here, if they are not provided in ITERATION. + """ + + # Predefine all the default values + self.data_misfit_tol = 0.01 + self.step_tol = 0.01 + self.lam = 1.0 + #self.lam_max = 1e10 + #self.lam_min = 0.01 + self.gamma = 2 + self.trunc_energy = 0.95 + self.iteration = 0 + + # Loop over options in ITERATION and extract the parameters we want + for i, opt in enumerate(list(zip(*self.keys_da['iteration']))[0]): + if opt == 'data_misfit_tol': + self.data_misfit_tol = self.keys_da['iteration'][i][1] + if opt == 'step_tol': + self.step_tol = self.keys_da['iteration'][i][1] + if opt == 'lambda': + self.lam = self.keys_da['iteration'][i][1] + if opt == 'lambda_max': + self.lam_max = self.keys_da['iteration'][i][1] + if opt == 'lambda_min': + self.lam_min = self.keys_da['iteration'][i][1] + if opt == 'lambda_factor': + self.gamma = self.keys_da['iteration'][i][1] + + if 'energy' in self.keys_da: + # initial energy (Remember to extract this) + self.trunc_energy = self.keys_da['energy'] + if self.trunc_energy > 1: # ensure that it is given as percentage + self.trunc_energy /= 100. + + + + diff --git a/pipt/update_schemes/gies/gies_rlmmac.py b/pipt/update_schemes/gies/gies_rlmmac.py new file mode 100644 index 0000000..80d03fa --- /dev/null +++ b/pipt/update_schemes/gies/gies_rlmmac.py @@ -0,0 +1,7 @@ + +from pipt.update_schemes.gies.gies_base import GIESMixIn +from pipt.update_schemes.gies.rlmmac_update import rlmmac_update + + +class gies_rlmmac(GIESMixIn, rlmmac_update): + pass \ No newline at end of file diff --git a/pipt/update_schemes/gies/rlmmac_update.py b/pipt/update_schemes/gies/rlmmac_update.py new file mode 100644 index 0000000..ba70301 --- /dev/null +++ b/pipt/update_schemes/gies/rlmmac_update.py @@ -0,0 +1,214 @@ +"""EnRML (IES) without the prior increment term.""" + +import numpy as np +from copy import deepcopy +import copy as cp +from scipy.linalg import solve, solve_banded, cholesky, lu_solve, lu_factor, inv +import pickle +import pipt.misc_tools.analysis_tools as at +from pipt.misc_tools.cov_regularization import _calc_loc + +class rlmmac_update(): + """ + Regularized Levenburg-Marquardt algorithm for Minimum Average Cost (RLM-MAC) problem, following the paper: + Luo, Xiaodong, et al. "Iterative ensemble smoother as an approximate solution to a regularized + minimum-average-cost problem: theory and applications." SPE Journal 2015: 962-982. + """ + + def update(self): + # calc the svd of the scaled data pertubation matrix + u_d, s_d, v_d = np.linalg.svd(self.pert_preddata, full_matrices=False) + aug_state = at.aug_state(self.current_state, self.list_states, self.cell_index) + + # remove the last singular value/vector. This is because numpy returns all ne values, while the last is actually + # zero. This part is a good place to include eventual additional truncation. + if self.trunc_energy < 1: + ti = (np.cumsum(s_d) / sum(s_d)) <= self.trunc_energy + u_d, s_d, v_d = u_d[:, ti].copy(), s_d[ti].copy(), v_d[ti, :].copy() + if 'localization' in self.keys_da: + if 'emp_cov' in self.keys_da and self.keys_da['emp_cov'] == 'yes': + if len(self.scale_data.shape) == 1: + E_hat = np.dot(np.expand_dims(self.scale_data ** (-1), + axis=1), np.ones((1, self.ne))) * self.E + x_0 = np.dot(np.diag(s_d[:] ** (-1)), np.dot(u_d[:, :].T, E_hat)) + Lam, z = np.linalg.eig(np.dot(x_0, x_0.T)) + else: + E_hat = solve(self.scale_data, self.E) + # x_0 = np.diag(s_d[:] ** (-1)) @ u_d[:, :].T @ E_hat + x_0 = np.dot(np.diag(s_d[:] ** (-1)), np.dot(u_d[:, :].T, E_hat)) + Lam, z = np.linalg.eig(np.dot(x_0, x_0.T)) + + # assume S_d = U S V^T; then Kalman gain K = S_m S_d^T (S_d S_d^T + gamma E E^T)^{-1} (d^o - g(m)) + # For this particular part A + # = S_d^T (S_d S_d^T + gamma E E^T)^{-1} + # = ... + # = V (I + gamma (S^-1 U^T E) (S^-1 U^T E)^T)^{-1} (S^-1 U^T) + # let S^-1 U^T E = Z L Z^T + # = V Z (I + gamma L)^{-1} Z^T S^-1 U^T + alpha = self.lam * np.sum(Lam**2) / len(Lam) + X = v_d.T @ z @ np.diag(1 / (alpha * Lam + np.eye(len(Lam)))) @ z.T @ np.diag(s_d[:] ** (-1)) @ u_d.T + #X = np.dot(np.dot(v_d.T, z), solve((self.lam + 1) * np.diag(Lam) + np.eye(len(Lam)), + # np.dot(u_d[:, :], np.dot(np.diag(s_d[:] ** (-1)).T, z)).T)) + + else: + alpha = self.lam * np.sum(s_d**2) / len(s_d) + X = v_d.T @ np.diag(s_d / (s_d ** 2 + alpha)) @ u_d.T + #X = np.dot(np.dot(v_d.T, np.diag(s_d)), solve(((self.lam + 1) * np.eye(len(s_d)) + np.diag(s_d ** 2)), u_d.T)) + + # we must perform localization + # store the size of all data + data_size = [[self.obs_data[int(time)][data].size if self.obs_data[int(time)][data] is not None else 0 + for data in self.list_datatypes] for time in self.assim_index[1]] + + f = self.keys_da['localization'] + + if f[1][0] == 'autoadaloc': + + # Mean state and perturbation matrix + mean_state = np.mean(aug_state, 1) + if 'emp_cov' in self.keys_da and self.keys_da['emp_cov'] == 'yes': + pert_state = (self.state_scaling**(-1))[:, None] * (aug_state - np.dot(np.resize(mean_state, (len(mean_state), 1)), + np.ones((1, self.ne)))) + else: + pert_state = (self.state_scaling**(-1) + )[:, None] * np.dot(aug_state, self.proj) + if len(self.scale_data.shape) == 1: + scaled_delta_data = np.dot(np.expand_dims(self.scale_data ** (-1), axis=1), + np.ones((1, pert_state.shape[1]))) * ( + self.real_obs_data - self.aug_pred_data[:, 0:self.ne]) + else: + scaled_delta_data = solve( + self.scale_data, (self.real_obs_data - self.aug_pred_data[:, 0:self.ne])) + + self.step = self.localization.auto_ada_loc(self.state_scaling[:, None] * pert_state, np.dot(X, scaled_delta_data), + self.list_states, + **{'prior_info': self.prior_info}) + elif 'localanalysis' in self.localization.loc_info and self.localization.loc_info['localanalysis']: + if 'distance' in self.localization.loc_info: + weight = _calc_loc(self.localization.loc_info['range'], self.localization.loc_info['distance'], + self.prior_info[self.list_states[0]], self.localization.loc_info['type'], self.ne) + else: + # if no distance, do full update + weight = np.ones((aug_state.shape[0], X.shape[1])) + mean_state = np.mean(aug_state, 1) + if 'emp_cov' in self.keys_da and self.keys_da['emp_cov'] == 'yes': + pert_state = (aug_state - np.dot(np.resize(mean_state, (len(mean_state), 1)), + np.ones((1, self.ne)))) + else: + pert_state = (aug_state - np.dot(np.resize(mean_state, (len(mean_state), 1)), + np.ones((1, self.ne)))) / (np.sqrt(self.ne - 1)) + + if len(self.scale_data.shape) == 1: + scaled_delta_data = np.dot(np.expand_dims(self.scale_data ** (-1), axis=1), + np.ones((1, pert_state.shape[1]))) * ( + self.real_obs_data - self.aug_pred_data) + else: + scaled_delta_data = solve( + self.scale_data, (self.real_obs_data - self.aug_pred_data)) + try: + self.step = weight.multiply( + np.dot(pert_state, X)).dot(scaled_delta_data) + except: + self.step = (weight*(np.dot(pert_state, X))).dot(scaled_delta_data) + + elif sum(['dist_loc' in el for el in f]) >= 1: + local_mask = self.localization.localize(self.list_datatypes, [self.keys_da['truedataindex'][int(elem)] + for elem in self.assim_index[1]], + self.list_states, self.ne, self.prior_info, data_size) + mean_state = np.mean(aug_state, 1) + if 'emp_cov' in self.keys_da and self.keys_da['emp_cov'] == 'yes': + pert_state = (aug_state - np.dot(np.resize(mean_state, (len(mean_state), 1)), + np.ones((1, self.ne)))) + else: + pert_state = (aug_state - np.dot(np.resize(mean_state, (len(mean_state), 1)), + np.ones((1, self.ne)))) / (np.sqrt(self.ne - 1)) + + if len(self.scale_data.shape) == 1: + scaled_delta_data = np.dot(np.expand_dims(self.scale_data ** (-1), axis=1), + np.ones((1, pert_state.shape[1]))) * ( + self.real_obs_data - self.aug_pred_data) + else: + scaled_delta_data = solve( + self.scale_data, (self.real_obs_data - self.aug_pred_data)) + + self.step = local_mask.multiply( + np.dot(pert_state, X)).dot(scaled_delta_data) + + else: + act_data_list = {} + count = 0 + for i in self.assim_index[1]: + for el in self.list_datatypes: + if self.real_obs_data[int(i)][el] is not None: + act_data_list[( + el, float(self.keys_da['truedataindex'][int(i)]))] = count + count += 1 + + well = [w for w in + set([el[0] for el in self.localization.loc_info.keys() if type(el) == tuple])] + times = [t for t in set( + [el[1] for el in self.localization.loc_info.keys() if type(el) == tuple])] + tot_dat_index = {} + for uniq_well in well: + tmp_index = [] + for t in times: + if (uniq_well, t) in act_data_list: + tmp_index.append(act_data_list[(uniq_well, t)]) + tot_dat_index[uniq_well] = tmp_index + if 'emp_cov' in self.keys_da and self.keys_da['emp_cov'] == 'yes': + emp_cov = True + else: + emp_cov = False + + self.step = at.parallel_upd(self.list_states, self.prior_info, self.current_state, X, + self.localization.loc_info, self.real_obs_data, self.aug_pred_data, + int(self.keys_fwd['parallel']), + actnum=self.localization.loc_info['actnum'], + field_dim=self.localization.loc_info['field'], + act_data_list=tot_dat_index, + scale_data=self.scale_data, + num_states=len( + [el for el in self.list_states]), + emp_d_cov=emp_cov) + self.step = at.aug_state(self.step, self.list_states) + + else: + # Mean state and perturbation matrix + mean_state = np.mean(aug_state, 1) + if 'emp_cov' in self.keys_da and self.keys_da['emp_cov'] == 'yes': + pert_state = (self.state_scaling**(-1))[:, None] * (aug_state - np.dot(np.resize(mean_state, (len(mean_state), 1)), + np.ones((1, self.ne)))) + else: + pert_state = (self.state_scaling**(-1) + )[:, None] * np.dot(aug_state, self.proj) + if 'emp_cov' in self.keys_da and self.keys_da['emp_cov'] == 'yes': + if len(self.scale_data.shape) == 1: + E_hat = np.dot(np.expand_dims(self.scale_data ** (-1), + axis=1), np.ones((1, self.ne))) * self.E + x_0 = np.dot(np.diag(s_d[:] ** (-1)), np.dot(u_d[:, :].T, E_hat)) + Lam, z = np.linalg.eig(np.dot(x_0, x_0.T)) + x_1 = np.dot(np.dot(u_d[:, :], np.dot(np.diag(s_d[:] ** (-1)).T, z)).T, + np.dot(np.expand_dims(self.scale_data ** (-1), axis=1), np.ones((1, self.ne))) * + (self.real_obs_data - self.aug_pred_data)) + else: + E_hat = solve(self.scale_data, self.E) + x_0 = np.dot(np.diag(s_d[:] ** (-1)), np.dot(u_d[:, :].T, E_hat)) + Lam, z = np.linalg.eig(np.dot(x_0, x_0.T)) + x_1 = np.dot(np.dot(u_d[:, :], np.dot(np.diag(s_d[:] ** (-1)).T, z)).T, + solve(self.scale_data, (self.real_obs_data - self.aug_pred_data))) + + x_2 = solve((self.lam + 1) * np.diag(Lam) + np.eye(len(Lam)), x_1) + x_3 = np.dot(np.dot(v_d.T, z), x_2) + delta_1 = np.dot(self.state_scaling[:, None] * pert_state, x_3) + self.step = delta_1 + else: + # Compute the approximate update (follow notation in paper) + if len(self.scale_data.shape) == 1: + x_1 = np.dot(u_d.T, np.dot(np.expand_dims(self.scale_data ** (-1), axis=1), np.ones((1, self.ne))) * + (self.real_obs_data - self.aug_pred_data[:, 0:self.ne])) + else: + x_1 = np.dot(u_d.T, solve(self.scale_data, + (self.real_obs_data - self.aug_pred_data[:, 0:self.ne]))) + x_2 = solve(((self.lam + 1) * np.eye(len(s_d)) + np.diag(s_d ** 2)), x_1) + x_3 = np.dot(np.dot(v_d.T, np.diag(s_d)), x_2) + self.step = np.dot(self.state_scaling[:, None] * pert_state, x_3) diff --git a/setup.py b/setup.py index 8f59cb8..2078a82 100644 --- a/setup.py +++ b/setup.py @@ -43,6 +43,8 @@ 'tomli-w', 'pyyaml', 'libecalc==8.23.1', # pin version to avoid frequent modifications - 'scikit-learn' + 'scikit-learn', + 'pylops' + ] + EXTRAS['doc'], ) diff --git a/simulator/calc_pem.py b/simulator/calc_pem.py new file mode 100644 index 0000000..6186f42 --- /dev/null +++ b/simulator/calc_pem.py @@ -0,0 +1,105 @@ +from simulator.opm import flow +from importlib import import_module +import datetime as dt +import numpy as np +import os +from misc import ecl, grdecl +import shutil +import glob +from subprocess import Popen, PIPE +import mat73 +from copy import deepcopy +from sklearn.cluster import KMeans +from sklearn.preprocessing import StandardScaler +from mako.lookup import TemplateLookup +from mako.runtime import Context + +# from pylops import avo +from pylops.utils.wavelets import ricker +from pylops.signalprocessing import Convolve1D +from misc.PyGRDECL.GRDECL_Parser import GRDECL_Parser # https://github.com/BinWang0213/PyGRDECL/tree/master +from scipy.interpolate import interp1d +from pipt.misc_tools.analysis_tools import store_ensemble_sim_information +from geostat.decomp import Cholesky +from simulator.eclipse import ecl_100 + + + +def calc_pem(self, time): + # fluid phases written to restart file from simulator run + phases = self.ecl_case.init.phases + + pem_input = {} + # get active porosity + tmp = self.ecl_case.cell_data('PORO') + if 'compaction' in self.pem_input: + multfactor = self.ecl_case.cell_data('PORV_RC', time) + + pem_input['PORO'] = np.array(multfactor[~tmp.mask] * tmp[~tmp.mask], dtype=float) + else: + pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float) + + # get active NTG if needed + if 'ntg' in self.pem_input: + if self.pem_input['ntg'] == 'no': + pem_input['NTG'] = None + else: + tmp = self.ecl_case.cell_data('NTG') + pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) + else: + tmp = self.ecl_case.cell_data('NTG') + pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) + + for var in phases: + tmp = self.ecl_case.cell_data(var, time) + pem_input[var] = np.array(tmp[~tmp.mask], dtype=float) # only active, and conv. to float + + if 'RS' in self.ecl_case.cell_data: + tmp = self.ecl_case.cell_data('RS', time) + pem_input['RS'] = np.array(tmp[~tmp.mask], dtype=float) + else: + pem_input['RS'] = None + print('RS is not a variable in the ecl_case') + + # extract pressure + tmp = self.ecl_case.cell_data('PRESSURE', time) + pem_input['PRESSURE'] = np.array(tmp[~tmp.mask], dtype=float) + + if 'press_conv' in self.pem_input: + pem_input['PRESSURE'] = pem_input['PRESSURE'] * self.pem_input['press_conv'] + + tmp = self.ecl_case.cell_data('PRESSURE', 1) + + if hasattr(self.pem, 'p_init'): + P_init = self.pem.p_init * np.ones(tmp.shape)[~tmp.mask] + else: + P_init = np.array(tmp[~tmp.mask], dtype=float) # initial pressure is first + + if 'press_conv' in self.pem_input: + P_init = P_init * self.pem_input['press_conv'] + + # extract saturations + if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: # This should be extended + saturations = [1 - (pem_input['SWAT'] + pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)] + for ph in phases] + elif 'OIL' in phases and 'GAS' in phases: # Smeaheia model + saturations = [pem_input['S{}'.format(ph)] for ph in phases] + else: + print('Type and number of fluids are unspecified in calc_pem') + + # fluid saturations in dictionary + tmp_s = {f'S{ph}': saturations[i] for i, ph in enumerate(phases)} + self.sats.extend([tmp_s]) + + # Get elastic parameters + if hasattr(self, 'ensemble_member') and (self.ensemble_member is not None) and \ + (self.ensemble_member >= 0): + self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'], + ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init, + ensembleMember=self.ensemble_member) + else: + self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'], + ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init) + + + diff --git a/simulator/flow_rock.py b/simulator/flow_rock.py index 3f34466..d2b6857 100644 --- a/simulator/flow_rock.py +++ b/simulator/flow_rock.py @@ -1,9 +1,12 @@ """Descriptive description.""" +from selectors import SelectSelector + from simulator.opm import flow from importlib import import_module import datetime as dt import numpy as np import os +import pandas as pd from misc import ecl, grdecl import shutil import glob @@ -12,7 +15,463 @@ from copy import deepcopy from sklearn.cluster import KMeans from sklearn.preprocessing import StandardScaler +from scipy.optimize import fsolve +from scipy.special import jv # Bessel function of the first kind +from scipy.integrate import quad +from scipy.special import j0 +from mako.lookup import TemplateLookup +from mako.runtime import Context +#import cProfile +#import pstats + +# from pylops import avo +from pylops.utils.wavelets import ricker +from pylops.signalprocessing import Convolve1D +import sys +#from PyGRDECL.GRDECL_Parser import GRDECL_Parser # https://github.com/BinWang0213/PyGRDECL/tree/master +from scipy.interpolate import interp1d +from scipy.interpolate import griddata +from pipt.misc_tools.analysis_tools import store_ensemble_sim_information +from geostat.decomp import Cholesky +#from simulator.eclipse import ecl_100 +from CoolProp.CoolProp import PropsSI # http://coolprop.org/#high-level-interface-example + +class mixIn_multi_data(): + + def find_cell_centre(self, grid): + # Find indices where the boolean array is True + indices = np.where(grid['ACTNUM']) + + coord = grid['COORD'] + zcorn = grid['ZCORN'] + + c, a, b = indices + # Calculate xt, yt, zt + xb = 0.25 * (coord[a, b, 0, 0] + coord[a, b + 1, 0, 0] + coord[a + 1, b, 0, 0] + coord[a + 1, b + 1, 0, 0]) + yb = 0.25 * (coord[a, b, 0, 1] + coord[a, b + 1, 0, 1] + coord[a + 1, b, 0, 1] + coord[a + 1, b + 1, 0, 1]) + zb = 0.25 * (coord[a, b, 0, 2] + coord[a, b + 1, 0, 2] + coord[a + 1, b, 0, 2] + coord[a + 1, b + 1, 0, 2]) + + xt = 0.25 * (coord[a, b, 1, 0] + coord[a, b + 1, 1, 0] + coord[a + 1, b, 1, 0] + coord[a + 1, b + 1, 1, 0]) + yt = 0.25 * (coord[a, b, 1, 1] + coord[a, b + 1, 1, 1] + coord[a + 1, b, 1, 1] + coord[a + 1, b + 1, 1, 1]) + zt = 0.25 * (coord[a, b, 1, 2] + coord[a, b + 1, 1, 2] + coord[a + 1, b, 1, 2] + coord[a + 1, b + 1, 1, 2]) + + # Calculate z, x, and y positions + z = (zcorn[c, 0, a, 0, b, 0] + zcorn[c, 0, a, 1, b, 0] + zcorn[c, 0, a, 0, b, 1] + zcorn[c, 0, a, 1, b, 1] + + zcorn[c, 1, a, 0, b, 0] + zcorn[c, 1, a, 1, b, 0] + zcorn[c, 1, a, 0, b, 1] + zcorn[c, 1, a, 1, b, 1]) / 8 + + x = xb + (xt - xb) * (z - zb) / (zt - zb) + y = yb + (yt - yb) * (z - zb) / (zt - zb) + + cell_centre = [x, y, z] + return cell_centre + + + def get_seabed_depths(self, file_path): + # Read the data while skipping the header comments + # We'll assume the header data ends before the numerical data + # The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead + water_depths = pd.read_csv(file_path, comment='#', sep=r'\s+', + header=None) # delim_whitespace=True, header=None) + + # Give meaningful column names: + water_depths.columns = ['x', 'y', 'z', 'column', 'row'] + + return water_depths + + def measurement_locations(self, grid, water_depth, pad=1500, dxy=3000, well_coord = None, dxy_fine = 1500, r0 = 5000): + + # Determine the size of the measurement area as defined by the field extent + cell_centre = self.find_cell_centre(grid) + x_min = np.min(cell_centre[0]) + x_max = np.max(cell_centre[0]) + y_min = np.min(cell_centre[1]) + y_max = np.max(cell_centre[1]) + + x_min -= pad + x_max += pad + y_min -= pad + y_max += pad + + x_span = x_max - x_min + y_span = y_max - y_min + + nx = int(np.ceil(x_span / dxy)) + ny = int(np.ceil(y_span / dxy)) + + x_vec = np.linspace(x_min, x_max, nx) + y_vec = np.linspace(y_min, y_max, ny) + x, y = np.meshgrid(x_vec, y_vec) + + # allow for finer measurement grid around injection well + if well_coord is not None: + # choose center point and radius for area with finer measurement grid + # y = 64, x = 34 # alpha well position Smeaheia + xc = cell_centre[well_coord[0]] + yc = cell_centre[well_coord[1]] + + # Fine grid covering bounding box of the circle (clamped to domain) + fx_min = max(x_min, xc - r0) + fx_max = min(x_max, xc + r0) + fy_min = max(y_min, yc - r0) + fy_max = min(y_max, yc + r0) + + pts_coarse = np.column_stack((x.ravel(), y.ravel())) + + if fx_max > fx_min and fy_max > fy_min: + nfx = int(np.ceil((fx_max - fx_min) / dxy_fine)) + 1 + nfy = int(np.ceil((fy_max - fy_min) / dxy_fine)) + 1 + x_fine = np.linspace(fx_min, fx_max, nfx) + y_fine = np.linspace(fy_min, fy_max, nfy) + xf, yf = np.meshgrid(x_fine, y_fine) + pts_fine = np.column_stack((xf.ravel(), yf.ravel())) + + # Keep only fine points inside the circle + d2 = (pts_fine[:, 0] - xc) ** 2 + (pts_fine[:, 1] - yc) ** 2 + mask_inside = d2 <= r0 ** 2 + pts_fine_inside = pts_fine[mask_inside] + + # remove the coarse points inside the circle + d2 = (pts_coarse[:, 0] - xc) ** 2 + (pts_coarse[:, 1] - yc) ** 2 + mask_inside = d2 <= r0 ** 2 + pts_coarse = pts_coarse[~mask_inside] + + # Combine and remove duplicates by rounding to a tolerance or using a structured array + # Use tolerance based on the smaller spacing + tol = min(dxy, dxy_fine) * 1e-3 + all_pts = np.vstack((pts_coarse, pts_fine_inside)) + + # Round coordinates to avoid floating point duplicates then use np.unique + # Determine digits to round so differences smaller than tol collapse + digits = max(0, int(-np.floor(np.log10(tol)))) + all_pts_rounded = np.round(all_pts, digits) + uniq_pts = np.unique(all_pts_rounded, axis=0) + x = uniq_pts[:, 0] + y = uniq_pts[:, 1] + + + pos = {'x': x.flatten(), 'y': y.flatten()} + + # Seabed map or water depth scalar depending on input + if isinstance(water_depth, float): + pos['z'] = np.ones_like(pos['x']) * water_depth + else: + pos['z'] = griddata((water_depth['x'], water_depth['y']), + np.abs(water_depth['z']), (pos['x'], pos['y']), + method='nearest') # z is positive downwards + return pos + + +class flow_rock(flow): + """ + Couple the OPM-flow simulator with a rock-physics simulator such that both reservoir quantities and petro-elastic + quantities can be calculated. Inherit the flow class, and use super to call similar functions. + """ + + def __init__(self, input_dict=None, filename=None, options=None): + super().__init__(input_dict) + self._getpeminfo(input_dict) + + self.date_slack = None + if 'date_slack' in input_dict: + self.date_slack = int(input_dict['date_slack']) + + # If we want to extract, or evaluate, something uniquely from the ensemble specific run we can + # run a user defined code to do this. + self.saveinfo = None + if 'savesiminfo' in input_dict: + # Make sure "ANALYSISDEBUG" gives a list + if isinstance(input_dict['savesiminfo'], list): + self.saveinfo = input_dict['savesiminfo'] + else: + self.saveinfo = [input_dict['savesiminfo']] + + self.scale = [] + + # Store dynamic variables in case they are provided in the state + self.state = None + self.no_flow = False + + def _getpeminfo(self, input_dict): + """ + Get, and return, flow and PEM modules + """ + if 'pem' in input_dict: + self.pem_input = {} + for elem in input_dict['pem']: + if elem[0] == 'model': # Set the petro-elastic model + self.pem_input['model'] = elem[1] + if elem[0] == 'depth': # provide the npz of depth values + self.pem_input['depth'] = elem[1] + if elem[0] == 'actnum': # the npz of actnum values + self.pem_input['actnum'] = elem[1] + if elem[0] == 'baseline': # the time for the baseline 4D measurement + self.pem_input['baseline'] = elem[1] + if elem[0] == 'vintage': + self.pem_input['vintage'] = elem[1] + if not type(self.pem_input['vintage']) == list: + self.pem_input['vintage'] = [elem[1]] + if elem[0] == 'ntg': + self.pem_input['ntg'] = elem[1] + if elem[0] == 'press_conv': + self.pem_input['press_conv'] = elem[1] + if elem[0] == 'compaction': + self.pem_input['compaction'] = True + if elem[0] == 'overburden': # the npz of overburden values + self.pem_input['overburden'] = elem[1] + if elem[0] == 'percentile': # use for scaling + self.pem_input['percentile'] = elem[1] + if elem[0] == 'phases': # get the fluid phases + self.pem_input['phases'] = elem[1] + if elem[0] == 'grid': # get the model grid + self.pem_input['grid'] = elem[1] + if elem[0] == 'param_file': # get model parameters required for pem + self.pem_input['param_file'] = elem[1] + + + pem = getattr(import_module('simulator.rockphysics.' + + self.pem_input['model'].split()[0]), self.pem_input['model'].split()[1]) + + self.pem = pem(self.pem_input) + + else: + self.pem = None + + def _get_pem_input(self, type, time=None): + if self.no_flow: # get variable from state + if any(type.lower() in key for key in self.state.keys()) and time > 0: + data = self.state[type.lower()+'_'+str(time)] + mask = np.zeros(data.shape, dtype=bool) + return np.ma.array(data=data, dtype=data.dtype, + mask=mask) + else: # read parameter from file + param_file = self.pem_input['param_file'] + npzfile = np.load(param_file) + parameter = npzfile[type] + npzfile.close() + data = parameter[:,self.ensemble_member] + mask = np.zeros(data.shape, dtype=bool) + return np.ma.array(data=data, dtype=data.dtype, + mask=mask) + else: # get variable of parameter from flow simulation + return self.ecl_case.cell_data(type,time) + + def calc_pem(self, time, time_index=None): + + if self.no_flow: + time_input = time_index + else: + time_input = time + + # fluid phases written given as input + phases = str.upper(self.pem_input['phases']) + phases = phases.split() + + pem_input = {} + tmp_dyn_var = {} + # get active porosity + tmp = self._get_pem_input('PORO') # self.ecl_case.cell_data('PORO') + if 'compaction' in self.pem_input: + multfactor = self._get_pem_input('PORV_RC', time_input) + pem_input['PORO'] = np.array(multfactor[~tmp.mask] * tmp[~tmp.mask], dtype=float) + else: + pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float) + + # get active NTG if needed + if 'ntg' in self.pem_input: + if self.pem_input['ntg'] == 'no': + pem_input['NTG'] = None + else: + tmp = self._get_pem_input('NTG') + pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) + else: + tmp = self._get_pem_input('NTG') + pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) + + if 'RS' in self.pem_input: #ecl_case.cell_data: # to be more robust! + tmp = self._get_pem_input('RS', time_input) + pem_input['RS'] = np.array(tmp[~tmp.mask], dtype=float) + else: + pem_input['RS'] = None + print('RS is not a variable in the ecl_case') + + # extract pressure + tmp = self._get_pem_input('PRESSURE', time_input) + pem_input['PRESSURE'] = np.array(tmp[~tmp.mask], dtype=float) + + # convert pressure from Bar to MPa + if 'press_conv' in self.pem_input and time_input == time: + pem_input['PRESSURE'] = pem_input['PRESSURE'] * self.pem_input['press_conv'] + + if hasattr(self.pem, 'p_init'): + P_init = self.pem.p_init * np.ones(tmp.shape)[~tmp.mask] + else: + P_init = np.array(tmp[~tmp.mask], dtype=float) # initial pressure is first + + if 'press_conv' in self.pem_input and time_input == time: + P_init = P_init * self.pem_input['press_conv'] + + # extract saturations + if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: # This should be extended + for var in phases: + if var in ['WAT', 'GAS']: + tmp = self._get_pem_input('S{}'.format(var), time_input) + pem_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) + pem_input['S{}'.format(var)] = np.clip(pem_input['S{}'.format(var)], 0, 1) + + pem_input['SOIL'] = np.clip(1 - (pem_input['SWAT'] + pem_input['SGAS']), 0, 1) + saturations = [ np.clip(1 - (pem_input['SWAT'] + pem_input['SGAS']), 0, 1) if ph == 'OIL' else pem_input['S{}'.format(ph)] + for ph in phases] + elif 'WAT' in phases and 'GAS' in phases: # Smeaheia model using OPM CO2Store + for var in phases: + if var in ['GAS']: + tmp = self._get_pem_input('S{}'.format(var), time_input) + pem_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) + pem_input['S{}'.format(var)] = np.clip(pem_input['S{}'.format(var)] , 0, 1) + pem_input['SWAT'] = 1 - pem_input['SGAS'] + saturations = [1 - (pem_input['SGAS']) if ph == 'WAT' else pem_input['S{}'.format(ph)] for ph in phases] + + elif 'OIL' in phases and 'GAS' in phases: # Original Smeaheia model + for var in phases: + if var in ['GAS']: + tmp = self._get_pem_input('S{}'.format(var), time_input) + pem_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) + pem_input['S{}'.format(var)] = np.clip(pem_input['S{}'.format(var)], 0, 1) + pem_input['SOIL'] = 1 - pem_input['SGAS'] + saturations = [1 - (pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)] for ph in phases] + + else: + + print('Type and number of fluids are unspecified in calc_pem') + + # fluid saturations in dictionary + # tmp_dyn_var = {f'S{ph}': saturations[i] for i, ph in enumerate(phases)} + for var in phases: + tmp_dyn_var[f'S{var}'] = pem_input[f'S{var}'] + + tmp_dyn_var['PRESSURE'] = pem_input['PRESSURE'] + self.dyn_var.extend([tmp_dyn_var]) + + if not self.no_flow: + keywords = self.ecl_case.arrays(time) + keywords = [s.strip() for s in keywords] # Remove leading/trailing spaces + #for key in self.all_data_types: + #if 'grav' in key: + densities = [] + for var in phases: + # fluid densities + dens = var + '_DEN' + if dens in keywords: + tmp = self._get_pem_input(dens, time_input) + pem_input[dens] = np.array(tmp[~tmp.mask], dtype=float) + # extract densities + densities.append(pem_input[dens]) + else: + densities = None + # pore volumes at each assimilation step + if 'RPORV' in keywords: + tmp = self._get_pem_input('RPORV', time_input) + pem_input['RPORV'] = np.array(tmp[~tmp.mask], dtype=float) + else: + densities = None + + # Get elastic parameters + if hasattr(self, 'ensemble_member') and (self.ensemble_member is not None) and \ + (self.ensemble_member >= 0): + self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'], + dens = densities, ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init, + ensembleMember=self.ensemble_member) + else: + self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'], + dens = densities, ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init) + + def setup_fwd_run(self, redund_sim): + super().setup_fwd_run(redund_sim=redund_sim) + + def run_fwd_sim(self, state, member_i, del_folder=True): + # The inherited simulator also has a run_fwd_sim. Call this. + self.ensemble_member = member_i + + # Check if dynamic variables are provided in the state. If that is the case, do not run flow simulator + if any('sgas' in key for key in state.keys()) or any('swat' in key for key in state.keys()) or any('pressure' in key for key in state.keys()): + self.state = {} + for key in state.keys(): + self.state[key] = state[key] + self.no_flow = True + #self.pred_data = self.extract_data(member_i) + #else: + self.pred_data = super().run_fwd_sim(state, member_i, del_folder=del_folder) + + return self.pred_data + + def call_sim(self, folder=None, wait_for_proc=False): + # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. + # Then, get the pem. + if not self.no_flow: + success = super().call_sim(folder, wait_for_proc) + else: + success = True + + if success: + self.ecl_case = ecl.EclipseCase( + 'En_' + str(self.ensemble_member) + os.sep + self.file + '.DATA') + phases = self.ecl_case.init.phases + self.dyn_var = [] + vintage = [] + # loop over seismic vintages + for v, assim_time in enumerate(self.pem_input['vintage']): + time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ + dt.timedelta(days=assim_time) + + self.calc_pem(time, v+1) + + # mask the bulk imp. to get proper dimensions + tmp_value = np.zeros(self.ecl_case.init.shape) + tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp + self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float, + mask=deepcopy(self.ecl_case.init.mask)) + # run filter + self.pem._filter() + vintage.append(deepcopy(self.pem.bulkimp)) + + if hasattr(self.pem, 'baseline'): # 4D measurement + base_time = dt.datetime(self.startDate['year'], self.startDate['month'], + self.startDate['day']) + dt.timedelta(days=self.pem.baseline) + + self.calc_pem(base_time, 0) + + # mask the bulk imp. to get proper dimensions + tmp_value = np.zeros(self.ecl_case.init.shape) + + tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp + # kill if values are inf or nan + assert not np.isnan(tmp_value).any() + assert not np.isinf(tmp_value).any() + self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float, + mask=deepcopy(self.ecl_case.init.mask)) + self.pem._filter() + + # 4D response + self.pem_result = [] + for i, elem in enumerate(vintage): + self.pem_result.append(elem - deepcopy(self.pem.bulkimp)) + else: + for i, elem in enumerate(vintage): + self.pem_result.append(elem) + + return success + + def extract_data(self, member): + # start by getting the data from the flow simulator + super().extract_data(member) + # get the sim2seis from file + for prim_ind in self.l_prim: + # Loop over all keys in pred_data (all data types) + for key in self.all_data_types: + if key in ['bulkimp']: + if self.true_prim[1][prim_ind] in self.pem_input['vintage']: + v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = self.pem_result[v].data.flatten() class flow_sim2seis(flow): """ @@ -55,7 +514,7 @@ def _getpeminfo(self, input_dict): self.pem_input['depth'] = elem[1] if elem[0] == 'actnum': # the npz of actnum values self.pem_input['actnum'] = elem[1] - if elem[0] == 'baseline': # the time for the baseline 4D measurment + if elem[0] == 'baseline': # the time for the baseline 4D measurement self.pem_input['baseline'] = elem[1] if elem[0] == 'vintage': self.pem_input['vintage'] = elem[1] @@ -96,7 +555,7 @@ def call_sim(self, folder=None, wait_for_proc=False): success = super().call_sim(folder, wait_for_proc) if success: - # need a if to check that we have correct sim2seis + # need an if to check that we have correct sim2seis # copy relevant sim2seis files into folder. for file in glob.glob('sim2seis_config/*'): shutil.copy(file, 'En_' + str(self.ensemble_member) + os.sep) @@ -106,64 +565,20 @@ def call_sim(self, folder=None, wait_for_proc=False): grid = self.ecl_case.grid() phases = self.ecl_case.init.phases - if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: # This should be extended - vintage = [] - # loop over seismic vintages - for v, assim_time in enumerate(self.pem_input['vintage']): - time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ + self.dyn_var = [] + vintage = [] + # loop over seismic vintages + for v, assim_time in enumerate(self.pem_input['vintage']): + time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ dt.timedelta(days=assim_time) - pem_input = {} - # get active porosity - tmp = self.ecl_case.cell_data('PORO') - if 'compaction' in self.pem_input: - multfactor = self.ecl_case.cell_data('PORV_RC', time) - - pem_input['PORO'] = np.array( - multfactor[~tmp.mask]*tmp[~tmp.mask], dtype=float) - else: - pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float) - # get active NTG if needed - if 'ntg' in self.pem_input: - if self.pem_input['ntg'] == 'no': - pem_input['NTG'] = None - else: - tmp = self.ecl_case.cell_data('NTG') - pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) - else: - tmp = self.ecl_case.cell_data('NTG') - pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) - - for var in ['SWAT', 'SGAS', 'PRESSURE', 'RS']: - tmp = self.ecl_case.cell_data(var, time) - # only active, and conv. to float - pem_input[var] = np.array(tmp[~tmp.mask], dtype=float) - - if 'press_conv' in self.pem_input: - pem_input['PRESSURE'] = pem_input['PRESSURE'] * \ - self.pem_input['press_conv'] - - tmp = self.ecl_case.cell_data('PRESSURE', 1) - if hasattr(self.pem, 'p_init'): - P_init = self.pem.p_init*np.ones(tmp.shape)[~tmp.mask] - else: - # initial pressure is first - P_init = np.array(tmp[~tmp.mask], dtype=float) - - if 'press_conv' in self.pem_input: - P_init = P_init*self.pem_input['press_conv'] - saturations = [1 - (pem_input['SWAT'] + pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)] - for ph in phases] - # Get the pressure - self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'], - ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init, - ensembleMember=self.ensemble_member) + self.calc_pem(time) #mali: update class inherent in flow_rock. Include calc_pem as method in flow_rock - grdecl.write(f'En_{str(self.ensemble_member)}/Vs{v+1}.grdecl', { + grdecl.write(f'En_{str(self.ensemble_member)}/Vs{v+1}.grdecl', { 'Vs': self.pem.getShearVel()*.1, 'DIMENS': grid['DIMENS']}, multi_file=False) - grdecl.write(f'En_{str(self.ensemble_member)}/Vp{v+1}.grdecl', { + grdecl.write(f'En_{str(self.ensemble_member)}/Vp{v+1}.grdecl', { 'Vp': self.pem.getBulkVel()*.1, 'DIMENS': grid['DIMENS']}, multi_file=False) - grdecl.write(f'En_{str(self.ensemble_member)}/rho{v+1}.grdecl', + grdecl.write(f'En_{str(self.ensemble_member)}/rho{v+1}.grdecl', {'rho': self.pem.getDens(), 'DIMENS': grid['DIMENS']}, multi_file=False) current_folder = os.getcwd() @@ -194,10 +609,12 @@ def extract_data(self, member): self.pred_data[prim_ind][key] = np.sum( np.abs(result[:, :, :, v]), axis=0).flatten() -class flow_rock(flow): +class flow_barycenter(flow): """ Couple the OPM-flow simulator with a rock-physics simulator such that both reservoir quantities and petro-elastic - quantities can be calculated. Inherit the flow class, and use super to call similar functions. + quantities can be calculated. Inherit the flow class, and use super to call similar functions. In the end, the + barycenter and moment of interia for the bulkimpedance objects, are returned as observations. The objects are + identified using k-means clustering, and the number of objects are determined using and elbow strategy. """ def __init__(self, input_dict=None, filename=None, options=None): @@ -221,6 +638,8 @@ def __init__(self, input_dict=None, filename=None, options=None): self.saveinfo = [input_dict['savesiminfo']] self.scale = [] + self.pem_result = [] + self.bar_result = [] def _getpeminfo(self, input_dict): """ @@ -251,6 +670,8 @@ def _getpeminfo(self, input_dict): self.pem_input['overburden'] = elem[1] if elem[0] == 'percentile': # use for scaling self.pem_input['percentile'] = elem[1] + if elem[0] == 'clusters': # number of clusters for each barycenter + self.pem_input['clusters'] = elem[1] pem = getattr(import_module('simulator.rockphysics.' + self.pem_input['model'].split()[0]), self.pem_input['model'].split()[1]) @@ -392,267 +813,1992 @@ def call_sim(self, folder=None, wait_for_proc=False): self.pem._filter() # 4D response - self.pem_result = [] for i, elem in enumerate(vintage): self.pem_result.append(elem - deepcopy(self.pem.bulkimp)) else: for i, elem in enumerate(vintage): self.pem_result.append(elem) + # Extract k-means centers and interias for each element in pem_result + if 'clusters' in self.pem_input: + npzfile = np.load(self.pem_input['clusters'], allow_pickle=True) + n_clusters_list = npzfile['n_clusters_list'] + npzfile.close() + else: + n_clusters_list = len(self.pem_result)*[2] + kmeans_kwargs = {"init": "random", "n_init": 10, "max_iter": 300, "random_state": 42} + for i, bulkimp in enumerate(self.pem_result): + std = np.std(bulkimp) + features = np.argwhere(np.squeeze(np.reshape(np.abs(bulkimp), self.ecl_case.init.shape,)) > 3 * std) + scaler = StandardScaler() + scaled_features = scaler.fit_transform(features) + kmeans = KMeans(n_clusters=n_clusters_list[i], **kmeans_kwargs) + kmeans.fit(scaled_features) + kmeans_center = np.squeeze(scaler.inverse_transform(kmeans.cluster_centers_)) # data / measurements + self.bar_result.append(np.append(kmeans_center, kmeans.inertia_)) + return success def extract_data(self, member): # start by getting the data from the flow simulator super().extract_data(member) - # get the sim2seis from file + # get the barycenters and inertias for prim_ind in self.l_prim: # Loop over all keys in pred_data (all data types) for key in self.all_data_types: - if key in ['bulkimp']: + if key in ['barycenter']: if self.true_prim[1][prim_ind] in self.pem_input['vintage']: v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) - self.pred_data[prim_ind][key] = self.pem_result[v].data.flatten() - -class flow_barycenter(flow): - """ - Couple the OPM-flow simulator with a rock-physics simulator such that both reservoir quantities and petro-elastic - quantities can be calculated. Inherit the flow class, and use super to call similar functions. In the end, the - barycenter and moment of interia for the bulkimpedance objects, are returned as observations. The objects are - identified using k-means clustering, and the number of objects are determined using and elbow strategy. - """ + self.pred_data[prim_ind][key] = self.bar_result[v].flatten() - def __init__(self, input_dict=None, filename=None, options=None): +class flow_avo(flow_rock, mixIn_multi_data): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): super().__init__(input_dict, filename, options) - self._getpeminfo(input_dict) - self.dum_file_root = 'dummy.txt' - self.dum_entry = str(0) - self.date_slack = None - if 'date_slack' in input_dict: - self.date_slack = int(input_dict['date_slack']) + assert 'avo' in input_dict, 'To do AVO simulation, please specify an "AVO" section in the "FWDSIM" part' + self._get_avo_info() - # If we want to extract, or evaluate, something uniquely from the ensemble specific run we can - # run a user defined code to do this. - self.saveinfo = None - if 'savesiminfo' in input_dict: - # Make sure "ANALYSISDEBUG" gives a list - if isinstance(input_dict['savesiminfo'], list): - self.saveinfo = input_dict['savesiminfo'] + def setup_fwd_run(self, **kwargs): + self.__dict__.update(kwargs) + + super().setup_fwd_run(redund_sim=None) + + def run_fwd_sim(self, state, member_i, del_folder=True): + """ + Setup and run the AVO forward simulator. + + Parameters + ---------- + state : dict + Dictionary containing the ensemble state. + + member_i : int + Index of the ensemble member. any index < 0 (e.g., -1) means the ground truth in synthetic case studies + + del_folder : bool, optional + Boolean to determine if the ensemble folder should be deleted. Default is False. + """ + + # The inherited simulator also has a run_fwd_sim. Call this. + self.ensemble_member = member_i + #return super().run_fwd_sim(state, member_i, del_folder=del_folder) + + + self.pred_data = super().run_fwd_sim(state, member_i, del_folder=del_folder) + return self.pred_data + + def call_sim(self, folder=None, wait_for_proc=False, run_reservoir_model=None, save_folder=None): + # replace the sim2seis part (which is unusable) by avo based on Pylops + + if folder is None: + folder = self.folder + else: + self.folder = folder + + if not self.no_flow: + # call call_sim in flow class (skip flow_rock, go directly to flow which is a parent of flow_rock) + success = super(flow_rock, self).call_sim(folder, wait_for_proc) + else: + success = True + + if success: + self.get_avo_result(folder, save_folder) + + return success + + def get_avo_result(self, folder, save_folder): + + if self.no_flow: + grid_file = self.pem_input['grid'] + grid = np.load(grid_file) + zcorn = grid['ZCORN'] + dz = np.diff(zcorn[:, 0, :, 0, :, 0], axis=0) + # Extract the last layer + last_layer = dz[-1, :, :] + # Reshape to ensure it has the same number of dimensions + last_layer = last_layer.reshape(1, dz.shape[1], dz.shape[2]) + # Concatenate to the original array along the first axis + dz = np.concatenate([dz, last_layer], axis=0) + f_dim = [grid['DIMENS'][2], grid['DIMENS'][1], grid['DIMENS'][0]] + else: + self.ecl_case = ecl.EclipseCase(folder + os.sep + self.file + '.DATA') if folder[-1] != os.sep \ + else ecl.EclipseCase(folder + self.file + '.DATA') + grid = self.ecl_case.grid() + zcorn = grid['ZCORN'] + ecl_init = ecl.EclipseInit(folder + os.sep + self.file + '.DATA') if folder[-1] != os.sep \ + else ecl.EclipseCase(folder + self.file + '.DATA') + dz = ecl_init.cell_data('DZ') + f_dim = [ecl_init.init.nk, ecl_init.init.nj, ecl_init.init.ni] + + + # ecl_init = ecl.EclipseInit(ecl_case) + # f_dim = [self.ecl_case.init.nk, self.ecl_case.init.nj, self.ecl_case.init.ni] + #f_dim = [self.NZ, self.NY, self.NX] + # phases = self.ecl_case.init.phases + self.dyn_var = [] + # coarsening of avo data - should be given as input in pipt + step_x = 1 + step_y = 1 + if 'baseline' in self.pem_input: # 4D measurement + base_time = dt.datetime(self.startDate['year'], self.startDate['month'], + self.startDate['day']) + dt.timedelta(days=self.pem_input['baseline']) + + + self.calc_pem(base_time,0) + # vp, vs, density in reservoir + self.calc_velocities(folder, save_folder, grid, -1, f_dim) + + if not self.no_flow: + # vp, vs, density in reservoir + vp, vs, rho = self.calc_velocities(folder, save_folder, grid, 0, f_dim) + + # avo data + # self._calc_avo_props() + avo_data_baseline, Rpp_baseline, vp_baseline, vs_baseline = self._calc_avo_props_active_cells(grid, vp, vs, rho, dz, zcorn) + kept_data = avo_data_baseline[::step_x, ::step_y, :].copy() + avo_data_baseline[:] = np.nan + avo_data_baseline[::step_x, ::step_y, :] = kept_data + # TODO: check which order to use, + # need to correlate with pipt/toml input-file and compression code + avo_baseline = avo_data_baseline.flatten(order="C") + avo_baseline = avo_baseline[~np.isnan(avo_baseline)] + #rho_baseline = rho_sample + tmp = self._get_pem_input('PRESSURE', base_time) + PRESSURE_baseline = np.array(tmp[~tmp.mask], dtype=float) + tmp = self._get_pem_input('SGAS', base_time) + SGAS_baseline = np.array(tmp[~tmp.mask], dtype=float) + print('OPM flow is used') else: - self.saveinfo = [input_dict['savesiminfo']] + file_name = f"avo_vint0_{folder}.npz" if folder[-1] != os.sep \ + else f"avo_vint0_{folder[:-1]}.npz" + + avo_baseline = np.load(file_name, allow_pickle=True)['avo_bl'] + #Rpp_baseline = np.load(file_name, allow_pickle=True)['Rpp_bl'] + #vs_baseline = np.load(file_name, allow_pickle=True)['Vs_bl'] + #vp_baseline = np.load(file_name, allow_pickle=True)['Vp_bl'] + #rho_baseline = np.load(file_name, allow_pickle=True)['Rho_bl'] + + vintage = [] + # loop over seismic vintages + for v, assim_time in enumerate(self.pem_input['vintage']): + time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ + dt.timedelta(days=assim_time) + + # extract dynamic variables from simulation run + self.calc_pem(time, v+1) + + # vp, vs, density in reservoir + vp, vs, rho = self.calc_velocities(folder, save_folder, grid, v+1, f_dim) + + # avo data + #self._calc_avo_props() + avo_data, Rpp, vp_sample, vs_sample = self._calc_avo_props_active_cells(grid, vp, vs, rho, dz, zcorn) + #avo_data = avo_data[::step_x, ::step_y, :] + # make mask for wavelet decomposition + kept_data = avo_data[::step_x, ::step_y, :].copy() + avo_data[:] = np.nan + avo_data[::step_x, ::step_y, :] = kept_data + mask = np.ones(np.shape(avo_data), dtype=bool) + mask[np.isnan(avo_data)]=False + np.savez(f'mask_{v}.npz', mask=mask) + #TODO: check order of flattening as above + avo = avo_data.flatten(order="C") + avo = avo[~np.isnan(avo)] + + tmp = self._get_pem_input('PRESSURE', time) + PRESSURE = np.array(tmp[~tmp.mask], dtype=float) + tmp = self._get_pem_input('SGAS', time) + SGAS = np.array(tmp[~tmp.mask], dtype=float) + + + # MLIE: implement 4D avo + if 'baseline' in self.pem_input: # 4D measurement + avo = avo - avo_baseline + Rpp = Rpp - Rpp_baseline + vs_sample = vs_sample - vs_baseline + vp_sample = vp_sample - vp_baseline + PRESSURE = PRESSURE - PRESSURE_baseline + SGAS = SGAS - SGAS_baseline + #rho = self.rho_sample - rho_baseline + print('Time-lapse avo') + #else: + # Rpp = self.Rpp + # Vs = self.vs_sample + # Vp = self.vp_sample + # rho = self.rho_sample + + + + # XLUO: self.ensemble_member < 0 => reference reservoir model in synthetic case studies + # the corresonding (noisy) data are observations in data assimilation + if 'add_synthetic_noise' in self.input_dict and self.ensemble_member < 0: + non_nan_idx = np.argwhere(~np.isnan(avo)) + data_std = np.std(avo[non_nan_idx]) + if self.input_dict['add_synthetic_noise'][0] == 'snr': + noise_std = np.sqrt(self.input_dict['add_synthetic_noise'][1]) * data_std + avo[non_nan_idx] += noise_std * np.random.randn(avo[non_nan_idx].size, 1) + else: + noise_std = 0.0 # simulated data don't contain noise + + vintage.append(deepcopy(avo)) + + if v == 0: + save_dic = {'avo': avo, 'noise_std': noise_std, 'Rpp': Rpp, 'Vs': vs_sample, 'Vp': vp_sample, 'PRESSURE': PRESSURE, 'SGAS': SGAS, **self.avo_config} + #save_dic = {'avo': avo, 'noise_std': noise_std, 'Rpp': Rpp, 'Vs': vs_sample, 'Vp': vp_sample, 'rho': rho_sample, #**self.avo_config, + # 'Vs_bl': vs_baseline, 'Vp_bl': vp_baseline, 'avo_bl': avo_baseline, 'Rpp_bl': Rpp_baseline, 'rho_bl': rho_baseline, **self.avo_config} + #save_dic = {'avo': avo, 'noise_std': noise_std, 'Rpp': Rpp, 'Vs': vs_sample, 'Vp': vp_sample, + # **self.avo_config} + else: + save_dic = {'avo': avo, 'noise_std': noise_std, 'Rpp': Rpp, 'Vs': vs_sample, 'Vp': vp_sample, 'PRESSURE': PRESSURE, 'SGAS': SGAS}#, 'Rpp': Rpp, 'Vs': vs_sample, 'Vp': vp_sample, + #'rho': rho_sample, **self.avo_config} + + if save_folder is not None: + file_name = save_folder + os.sep + f"avo_vint{v}.npz" if save_folder[-1] != os.sep \ + else save_folder + f"avo_vint{v}.npz" + #np.savez(file_name, **save_dic) + else: + file_name = folder + os.sep + f"avo_vint{v}.npz" if folder[-1] != os.sep \ + else folder + f"avo_vint{v}.npz" + file_name_rec = 'Ensemble_results/' + f"avo_vint{v}_{folder}.npz" if folder[-1] != os.sep \ + else 'Ensemble_results/' + f"avo_vint{v}_{folder[:-1]}.npz" + np.savez(file_name_rec, **save_dic) + # with open(file_name, "wb") as f: + # dump(**save_dic, f) + np.savez(file_name, **save_dic) + # 4D response + self.avo_result = [] + for i, elem in enumerate(vintage): + self.avo_result.append(elem) + + def calc_velocities(self, folder, save_folder, grid, v, f_dim): + # The properties in pem are only given in the active cells + # indices of active cells: + + true_indices = np.where(grid['ACTNUM']) + + vp = np.full(f_dim, np.nan) + vp[true_indices] = (self.pem.getBulkVel()) + vs = np.full(f_dim, np.nan) + vs[true_indices] = (self.pem.getShearVel()) + rho = np.full(f_dim, np.nan) + rho[true_indices] = (self.pem.getDens()) + + + + ## Debug + #bulkmod = np.full(f_dim, np.nan) + #bulkmod[true_indices] = self.pem.getBulkMod() + #self.shearmod = np.full(f_dim, np.nan) + #self.shearmod[true_indices] = self.pem.getShearMod() + #self.poverburden = np.full(f_dim, np.nan) + #self.poverburden[true_indices] = self.pem.getOverburdenP() + #self.pressure = np.full(f_dim, np.nan) + #self.pressure[true_indices] = self.pem.getPressure() + #self.peff = np.full(f_dim, np.nan) + #self.peff[true_indices] = self.pem.getPeff() + #porosity = np.full(f_dim, np.nan) + #porosity[true_indices] = self.pem.getPorosity() + #if self.dyn_var: + # sgas = np.full(f_dim, np.nan) + # sgas[true_indices] = self.dyn_var[v]['SGAS'] + # soil = np.full(f_dim, np.nan) + # soil[true_indices] = self.dyn_var[v]['SOIL'] + # pdyn = np.full(f_dim, np.nan) + # pdyn[true_indices] = self.dyn_var[v]['PRESSURE'] + # + #if self.dyn_var is None: + # save_dic = {'vp': vp, 'vs': vs, 'rho': rho}#, 'bulkmod': self.bulkmod, 'shearmod': self.shearmod, + # #'Pov': self.poverburden, 'P': self.pressure, 'Peff': self.peff, 'por': porosity} # for debugging + #else: + # save_dic = {'vp': vp, 'vs': vs, 'rho': rho}#, 'por': porosity, 'sgas': sgas, 'Pd': pdyn} + + #if save_folder is not None: + # file_name = save_folder + os.sep + f"vp_vs_rho_vint{v}.npz" if save_folder[-1] != os.sep \ + # else save_folder + f"vp_vs_rho_vint{v}.npz" + # np.savez(file_name, **save_dic) + #else: + # file_name_rec = 'Ensemble_results/' + f"vp_vs_rho_vint{v}_{folder}.npz" if folder[-1] != os.sep \ + # else 'Ensemble_results/' + f"vp_vs_rho_vint{v}_{folder[:-1]}.npz" + # np.savez(file_name_rec, **save_dic) + # for debugging + return vp, vs, rho + + def extract_data(self, member): + # start by getting the data from the flow simulator + super(flow_rock, self).extract_data(member) + + # get the avo from file + for prim_ind in self.l_prim: + # Loop over all keys in pred_data (all data types) + for key in self.all_data_types: + if 'avo' in key: + if self.true_prim[1][prim_ind] in self.pem_input['vintage']: + idx = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) + filename = self.folder + os.sep + key + '_vint' + str(idx) + '.npz' if self.folder[-1] != os.sep \ + else self.folder + key + '_vint' + str(idx) + '.npz' + with np.load(filename) as f: + self.pred_data[prim_ind][key] = f[key] + # + #v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) + #self.pred_data[prim_ind][key] = self.avo_result[v].flatten() + + def _get_avo_info(self, avo_config=None): + """ + AVO configuration + """ + # list of configuration parameters in the "AVO" section + config_para_list = ['dz', 'tops', 'angle', 'frequency', 'wave_len', 'vp_shale', 'vs_shale', + 'den_shale', 't_min', 't_max', 't_sampling', 'pp_func'] + if 'avo' in self.input_dict: + self.avo_config = {} + for elem in self.input_dict['avo']: + assert elem[0] in config_para_list, f'Property {elem[0]} not supported' + if elem[0] == 'vintage' and not isinstance(elem[1], list): + elem[1] = [elem[1]] + self.avo_config[elem[0]] = elem[1] + + # if only one angle is considered, convert self.avo_config['angle'] into a list, as required later + if isinstance(self.avo_config['angle'], float): + self.avo_config['angle'] = [self.avo_config['angle']] + + # self._get_DZ(file=self.avo_config['dz']) # =>self.DZ + kw_file = {'DZ': self.avo_config['dz'], 'TOPS': self.avo_config['tops']} + self._get_props(kw_file) + self.overburden = self.pem_input['overburden'] + + # make sure that the "pylops" package is installed + # See https://github.com/PyLops/pylops + self.pp_func = getattr(import_module('pylops.avo.avo'), self.avo_config['pp_func']) + + else: + self.avo_config = None + + def _get_props(self, kw_file): + # extract properties (specified by keywords) in (possibly) different files + # kw_file: a dictionary contains "keyword: file" pairs + # Note that all properties are reshaped into the reservoir model dimension (NX, NY, NZ) + # using the "F" order + for kw in kw_file: + file = kw_file[kw] + if file.endswith('.npz'): + with np.load(file) as f: + exec(f'self.{kw} = f[ "{kw}" ]') + self.NX, self.NY, self.NZ = f['NX'], f['NY'], f['NZ'] + else: + try: + self.NX = int(self.input_dict['dimension'][0]) + self.NY = int(self.input_dict['dimension'][1]) + self.NZ = int(self.input_dict['dimension'][2]) + except: + for item in self.input_dict['pem']: + if item[0] == 'dimension': + dimension = item[1] + break + self.NX = int(dimension[0]) + self.NY = int(dimension[1]) + self.NZ = int(dimension[2]) + # reader = GRDECL_Parser(filename=file) + # reader.read_GRDECL() + # exec(f"self.{kw} = reader.{kw}.reshape((reader.NX, reader.NY, reader.NZ), order='F')") + # self.NX, self.NY, self.NZ = reader.NX, reader.NY, reader.NZ + # eval(f'np.savez("./{kw}.npz", {kw}=self.{kw}, NX=self.NX, NY=self.NY, NZ=self.NZ)') + + def _calc_avo_props(self, dt=0.0005): + # dt is the fine resolution sampling rate + # convert properties in reservoir model to time domain + vp_shale = self.avo_config['vp_shale'] # scalar value (code may not work for matrix value) + vs_shale = self.avo_config['vs_shale'] # scalar value + rho_shale = self.avo_config['den_shale'] # scalar value + + # Two-way travel time of the top of the reservoir + # TOPS[:, :, 0] corresponds to the depth profile of the reservoir top on the first layer + top_res = 2 * self.TOPS[:, :, 0] / vp_shale + + # Cumulative traveling time through the reservoir in vertical direction + cum_time_res = np.cumsum(2 * self.DZ / self.vp, axis=2) + top_res[:, :, np.newaxis] + + # assumes underburden to be constant. No reflections from underburden. Hence set traveltime to underburden very large + underburden = top_res + np.max(cum_time_res) + + # total travel time + # cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res), axis=2) + cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res, underburden[:, :, np.newaxis]), axis=2) + + + # add overburden and underburden of Vp, Vs and Density + vp = np.concatenate((vp_shale * np.ones((self.NX, self.NY, 1)), + self.vp, vp_shale * np.ones((self.NX, self.NY, 1))), axis=2) + vs = np.concatenate((vs_shale * np.ones((self.NX, self.NY, 1)), + self.vs, vs_shale * np.ones((self.NX, self.NY, 1))), axis=2) + + #rho = np.concatenate((rho_shale * np.ones((self.NX, self.NY, 1)) * 0.001, # kg/m^3 -> k/cm^3 + # self.rho, rho_shale * np.ones((self.NX, self.NY, 1)) * 0.001), axis=2) + rho = np.concatenate((rho_shale * np.ones((self.NX, self.NY, 1)), + self.rho, rho_shale * np.ones((self.NX, self.NY, 1))), axis=2) + + # search for the lowest grid cell thickness and sample the time according to + # that grid thickness to preserve the thin layer effect + time_sample = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], dt) + if time_sample.shape[0] == 1: + time_sample = time_sample.reshape(-1) + time_sample = np.tile(time_sample, (self.NX, self.NY, 1)) + + vp_sample = np.tile(vp[:, :, 1][..., np.newaxis], (1, 1, time_sample.shape[2])) + vs_sample = np.tile(vs[:, :, 1][..., np.newaxis], (1, 1, time_sample.shape[2])) + rho_sample = np.tile(rho[:, :, 1][..., np.newaxis], (1, 1, time_sample.shape[2])) + + for m in range(self.NX): + for l in range(self.NY): + for k in range(time_sample.shape[2]): + # find the right interval of time_sample[m, l, k] belonging to, and use + # this information to allocate vp, vs, rho + idx = np.searchsorted(cum_time[m, l, :], time_sample[m, l, k], side='left') + idx = idx if idx < len(cum_time[m, l, :]) else len(cum_time[m, l, :]) - 1 + vp_sample[m, l, k] = vp[m, l, idx] + vs_sample[m, l, k] = vs[m, l, idx] + rho_sample[m, l, k] = rho[m, l, idx] + + + + + # Ricker wavelet + wavelet, t_axis, wav_center = ricker(np.arange(0, self.avo_config['wave_len'], dt), + f0=self.avo_config['frequency']) + + + # Travel time corresponds to reflectivity series + t = time_sample[:, :, 0:-1] + + # interpolation time + t_interp = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], self.avo_config['t_sampling']) + trace_interp = np.zeros((self.NX, self.NY, len(t_interp))) + + # number of pp reflection coefficients in the vertical direction + + nz_rpp = vp_sample.shape[2] - 1 + + for i in range(len(self.avo_config['angle'])): + angle = self.avo_config['angle'][i] + Rpp = self.pp_func(vp_sample[:, :, :-1], vs_sample[:, :, :-1], rho_sample[:, :, :-1], + vp_sample[:, :, 1:], vs_sample[:, :, 1:], rho_sample[:, :, 1:], angle) + + for m in range(self.NX): + for l in range(self.NY): + # convolution with the Ricker wavelet + conv_op = Convolve1D(nz_rpp, h=wavelet, offset=wav_center, dtype="float32") + w_trace = conv_op * Rpp[m, l, :] + + # Sample the trace into regular time interval + f = interp1d(np.squeeze(t[m, l, :]), np.squeeze(w_trace), + kind='nearest', fill_value='extrapolate') + trace_interp[m, l, :] = f(t_interp) + + if i == 0: + avo_data = trace_interp # 3D + elif i == 1: + avo_data = np.stack((avo_data, trace_interp), axis=-1) # 4D + else: + avo_data = np.concatenate((avo_data, trace_interp[:, :, :, np.newaxis]), axis=3) # 4D + + self.avo_data = avo_data + + def _calc_avo_props_active_cells(self, grid, vp, vs, rho, dz, zcorn, dt=0.0005): + # dt is the fine resolution sampling rate + # convert properties in reservoir model to time domain + vp_shale = self.avo_config['vp_shale'] # scalar value (code may not work for matrix value) + vs_shale = self.avo_config['vs_shale'] # scalar value + rho_shale = self.avo_config['den_shale'] # scalar value + + + actnum = grid['ACTNUM'] + # Find indices where the boolean array is True + active_indices = np.where(actnum) + c, a, b = active_indices + + # Two-way travel time tp the top of the reservoir + top_res = 2 * zcorn[0, 0, :, 0, :, 0] / vp_shale + + # depth difference between cells in z-direction: + depth_differences = dz#(zcorn[:, 0, :, 0, :, 0] , axis=0) + + + # Cumulative traveling time through the reservoir in vertical direction + #cum_time_res = 2 * zcorn[:, 0, :, 0, :, 0] / self.vp + top_res[np.newaxis, :, :] + cum_time_res = np.cumsum(2 * depth_differences / vp, axis=0) + top_res[np.newaxis, :, :] + # assumes under burden to be constant. No reflections from under burden. Hence set travel time to under burden very large + underburden = top_res + np.nanmax(cum_time_res) + + # total travel time + # cum_time = np.concat enate((top_res[:, :, np.newaxis], cum_time_res), axis=2) + cum_time = np.concatenate((top_res[np.newaxis, :, :], cum_time_res, underburden[np.newaxis, :, :]), axis=0) + + # add overburden and underburden values for Vp, Vs and Density + vp = np.concatenate((vp_shale * np.ones((1, self.NY, self.NX)), + vp, vp_shale * np.ones((1, self.NY, self.NX))), axis=0) + vs = np.concatenate((vs_shale * np.ones((1, self.NY, self.NX)), + vs, vs_shale * np.ones((1, self.NY, self.NX))), axis=0) + rho = np.concatenate((rho_shale * np.ones((1, self.NY, self.NX)), + rho, rho_shale * np.ones((1, self.NY, self.NX))), axis=0) + + + # Combine a and b into a 2D array (each column represents a vector) + ab = np.column_stack((a, b)) + + # Extract unique rows and get the indices of those unique rows + unique_rows, indices = np.unique(ab, axis=0, return_index=True) + + # search for the lowest grid cell thickness and sample the time according to + # that grid thickness to preserve the thin layer effect + time_sample = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], dt) + if time_sample.shape[0] == 1: + time_sample = time_sample.reshape(-1) + #time_sample = np.tile(time_sample, (len(indices), 1)) + time_sample = np.tile(time_sample, (self.NX, self.NY, 1)) + + #vp_sample = vp_shale * np.ones((self.NX, self.NY, time_sample.shape[2])) + #vs_sample = vs_shale * np.ones((self.NX, self.NY, time_sample.shape[2])) + #rho_sample = rho_shale * np.ones((self.NX, self.NY, time_sample.shape[2])) + vp_sample = np.full([self.NX, self.NY, time_sample.shape[2]], np.nan) + vs_sample = np.full([self.NX, self.NY, time_sample.shape[2]], np.nan) + rho_sample = np.full([self.NX, self.NY, time_sample.shape[2]], np.nan) + + + for ind in range(len(indices)): + for k in range(time_sample.shape[2]): + # find the right interval of time_sample[m, l, k] belonging to, and use + # this information to allocate vp, vs, rho + idx = np.searchsorted(cum_time[:, a[indices[ind]], b[indices[ind]]], time_sample[b[indices[ind]], a[indices[ind]], k], side='left') + idx = idx if idx < len(cum_time[:, a[indices[ind]], b[indices[ind]]]) else len( + cum_time[:,a[indices[ind]], b[indices[ind]]]) - 1 + vp_sample[b[indices[ind]], a[indices[ind]], k] = vp[idx, a[indices[ind]], b[indices[ind]]] + vs_sample[b[indices[ind]], a[indices[ind]], k] = vs[idx, a[indices[ind]], b[indices[ind]]] + rho_sample[b[indices[ind]], a[indices[ind]], k] = rho[idx, a[indices[ind]], b[indices[ind]]] + + # Ricker wavelet + wavelet, t_axis, wav_center = ricker(np.arange(0, self.avo_config['wave_len']-dt, dt), + f0=self.avo_config['frequency']) + + # Travel time corresponds to reflectivity series + #t = time_sample[:, 0:-1] + t = time_sample[:, :, 0:-1] + + # interpolation time + t_interp = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], self.avo_config['t_sampling']) + #trace_interp = np.zeros((len(indices), len(t_interp))) + #trace_interp = np.zeros((self.NX, self.NY, len(t_interp))) + trace_interp = np.full([self.NX, self.NY, len(t_interp)], np.nan) + + # number of pp reflection coefficients in the vertical direction + nz_rpp = vp_sample.shape[2] - 1 + conv_op = Convolve1D(nz_rpp, h=wavelet, offset=wav_center, dtype="float32") + + avo_data = [] + Rpp = [] + for i in range(len(self.avo_config['angle'])): + angle = self.avo_config['angle'][i] + Rpp = self.pp_func(vp_sample[:, :, :-1], vs_sample[:, :, :-1], rho_sample[:, :, :-1], + vp_sample[:, :, 1:], vs_sample[:, :, 1:], rho_sample[:, :, 1:], angle) + + for ind in range(len(indices)): + # convolution with the Ricker wavelet + + w_trace = conv_op * Rpp[b[indices[ind]], a[indices[ind]], :] + + # Sample the trace into regular time interval + f = interp1d(np.squeeze(t[b[indices[ind]], a[indices[ind]], :]), np.squeeze(w_trace), + kind='nearest', fill_value='extrapolate') + trace_interp[b[indices[ind]], a[indices[ind]], :] = f(t_interp) + + if i == 0: + avo_data = trace_interp # 3D + elif i == 1: + avo_data = np.stack((avo_data, trace_interp), axis=-1) # 4D + else: + avo_data = np.concatenate((avo_data, trace_interp[:, :, np.newaxis]), axis=2) # 4D + + return avo_data, Rpp, vp_sample, vs_sample + #self.avo_data = avo_data + #self.Rpp = Rpp + #self.vp_sample = vp_sample + #self.vs_sample = vs_sample + #self.rho_sample = rho_sample + + + + def _calc_avo_props_active_cells_org(self, grid, vp, vs, rho, dz, zcorn, dt=0.0005): + # dt is the fine resolution sampling rate + # convert properties in reservoir model to time domain + vp_shale = self.avo_config['vp_shale'] # scalar value (code may not work for matrix value) + vs_shale = self.avo_config['vs_shale'] # scalar value + rho_shale = self.avo_config['den_shale'] # scalar value + + + actnum = grid['ACTNUM'] + # Find indices where the boolean array is True + active_indices = np.where(actnum) + # # # + + # Two-way travel time tp the top of the reservoir + + + c, a, b = active_indices + + # Two-way travel time tp the top of the reservoir + top_res = 2 * zcorn[0, 0, :, 0, :, 0] / vp_shale + + # depth difference between cells in z-direction: + depth_differences = dz#(zcorn[:, 0, :, 0, :, 0] , axis=0) + + + # Cumulative traveling time through the reservoir in vertical direction + #cum_time_res = 2 * zcorn[:, 0, :, 0, :, 0] / self.vp + top_res[np.newaxis, :, :] + cum_time_res = np.cumsum(2 * depth_differences / vp, axis=0) + top_res[np.newaxis, :, :] + # assumes under burden to be constant. No reflections from under burden. Hence set travel time to under burden very large + underburden = top_res + np.nanmax(cum_time_res) + + # total travel time + # cum_time = np.concat enate((top_res[:, :, np.newaxis], cum_time_res), axis=2) + cum_time = np.concatenate((top_res[np.newaxis, :, :], cum_time_res, underburden[np.newaxis, :, :]), axis=0) + + # add overburden and underburden values for Vp, Vs and Density + vp = np.concatenate((vp_shale * np.ones((1, self.NY, self.NX)), + vp, vp_shale * np.ones((1, self.NY, self.NX))), axis=0) + vs = np.concatenate((vs_shale * np.ones((1, self.NY, self.NX)), + vs, vs_shale * np.ones((1, self.NY, self.NX))), axis=0) + rho = np.concatenate((rho_shale * np.ones((1, self.NY, self.NX)), + rho, rho_shale * np.ones((1, self.NY, self.NX))), axis=0) + + + # Combine a and b into a 2D array (each column represents a vector) + ab = np.column_stack((a, b)) + + # Extract unique rows and get the indices of those unique rows + unique_rows, indices = np.unique(ab, axis=0, return_index=True) + + # search for the lowest grid cell thickness and sample the time according to + # that grid thickness to preserve the thin layer effect + time_sample = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], dt) + if time_sample.shape[0] == 1: + time_sample = time_sample.reshape(-1) + time_sample = np.tile(time_sample, (len(indices), 1)) + + vp_sample = np.zeros((len(indices), time_sample.shape[1])) + vs_sample = np.zeros((len(indices), time_sample.shape[1])) + rho_sample = np.zeros((len(indices), time_sample.shape[1])) + + for ind in range(len(indices)): + for k in range(time_sample.shape[1]): + # find the right interval of time_sample[m, l, k] belonging to, and use + # this information to allocate vp, vs, rho + idx = np.searchsorted(cum_time[:, a[indices[ind]], b[indices[ind]]], time_sample[ind, k], side='left') + idx = idx if idx < len(cum_time[:, a[indices[ind]], b[indices[ind]]]) else len( + cum_time[:,a[indices[ind]], b[indices[ind]]]) - 1 + vp_sample[ind, k] = vp[idx, a[indices[ind]], b[indices[ind]]] + vs_sample[ind, k] = vs[idx, a[indices[ind]], b[indices[ind]]] + rho_sample[ind, k] = rho[idx, a[indices[ind]], b[indices[ind]]] + + # Ricker wavelet + wavelet, t_axis, wav_center = ricker(np.arange(0, self.avo_config['wave_len']-dt, dt), + f0=self.avo_config['frequency']) + + # Travel time corresponds to reflectivity series + t = time_sample[:, 0:-1] + + # interpolation time + t_interp = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], self.avo_config['t_sampling']) + trace_interp = np.zeros((len(indices), len(t_interp))) + + # number of pp reflection coefficients in the vertical direction + nz_rpp = vp_sample.shape[1] - 1 + conv_op = Convolve1D(nz_rpp, h=wavelet, offset=wav_center, dtype="float32") + + avo_data = [] + Rpp = [] + for i in range(len(self.avo_config['angle'])): + angle = self.avo_config['angle'][i] + Rpp = self.pp_func(vp_sample[:, :-1], vs_sample[:, :-1], rho_sample[:, :-1], + vp_sample[:, 1:], vs_sample[:, 1:], rho_sample[:, 1:], angle) + + for ind in range(len(indices)): + # convolution with the Ricker wavelet + + w_trace = conv_op * Rpp[ind, :] + + # Sample the trace into regular time interval + f = interp1d(np.squeeze(t[ind, :]), np.squeeze(w_trace), + kind='nearest', fill_value='extrapolate') + trace_interp[ind, :] = f(t_interp) + + if i == 0: + avo_data = trace_interp # 3D + elif i == 1: + avo_data = np.stack((avo_data, trace_interp), axis=-1) # 4D + else: + avo_data = np.concatenate((avo_data, trace_interp[:, :, np.newaxis]), axis=2) # 4D + + return avo_data, Rpp, vp_sample, vs_sample, rho_sample + #self.avo_data = avo_data + #self.Rpp = Rpp + #self.vp_sample = vp_sample + #self.vs_sample = vs_sample + #self.rho_sample = rho_sample + + def _calc_avo_props_active_cells_org(self, grid, dt=0.0005): + # dt is the fine resolution sampling rate + # convert properties in reservoir model to time domain + vp_shale = self.avo_config['vp_shale'] # scalar value (code may not work for matrix value) + vs_shale = self.avo_config['vs_shale'] # scalar value + rho_shale = self.avo_config['den_shale'] # scalar value + + # check if Nz, is at axis = 0, then transpose to dimensions, Nx, ny, Nz + if grid['ACTNUM'].shape[0] == self.NZ: + vp = np.transpose(self.vp, (2, 1, 0)) + vs = np.transpose(self.vs, (2, 1, 0)) + rho = np.transpose(self.rho, (2, 1, 0)) + actnum = np.transpose(grid['ACTNUM'], (2, 1, 0)) + else: + actnum = grid['ACTNUM'] + vp = self.vp + vs = self.vs + rho = self.rho + # # # + + # Two-way travel time of the top of the reservoir + # TOPS[:, :, 0] corresponds to the depth profile of the reservoir top on the first layer + top_res = 2 * self.TOPS[:, :, 0] / vp_shale + + # Cumulative traveling time through the reservoir in vertical direction + cum_time_res = np.nancumsum(2 * self.DZ / vp, axis=2) + top_res[:, :, np.newaxis] + + # assumes under burden to be constant. No reflections from under burden. Hence set travel time to under burden very large + underburden = top_res + np.max(cum_time_res) + + # total travel time + # cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res), axis=2) + cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res, underburden[:, :, np.newaxis]), axis=2) + + # add overburden and underburden of Vp, Vs and Density + vp = np.concatenate((vp_shale * np.ones((self.NX, self.NY, 1)), + vp, vp_shale * np.ones((self.NX, self.NY, 1))), axis=2) + vs = np.concatenate((vs_shale * np.ones((self.NX, self.NY, 1)), + vs, vs_shale * np.ones((self.NX, self.NY, 1))), axis=2) + #rho = np.concatenate((rho_shale * np.ones((self.NX, self.NY, 1)) * 0.001, # kg/m^3 -> k/cm^3 + # self.rho, rho_shale * np.ones((self.NX, self.NY, 1)) * 0.001), axis=2) + rho = np.concatenate((rho_shale * np.ones((self.NX, self.NY, 1)), + rho, rho_shale * np.ones((self.NX, self.NY, 1))), axis=2) + + # get indices of active cells + + indices = np.where(actnum) + a, b, c = indices + # Combine a and b into a 2D array (each column represents a vector) + ab = np.column_stack((a, b)) + + # Extract unique rows and get the indices of those unique rows + unique_rows, indices = np.unique(ab, axis=0, return_index=True) + + + # search for the lowest grid cell thickness and sample the time according to + # that grid thickness to preserve the thin layer effect + time_sample = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], dt) + if time_sample.shape[0] == 1: + time_sample = time_sample.reshape(-1) + time_sample = np.tile(time_sample, (len(indices), 1)) + + vp_sample = np.zeros((len(indices), time_sample.shape[1])) + vs_sample = np.zeros((len(indices), time_sample.shape[1])) + rho_sample = np.zeros((len(indices), time_sample.shape[1])) + + + for ind in range(len(indices)): + for k in range(time_sample.shape[1]): + # find the right interval of time_sample[m, l, k] belonging to, and use + # this information to allocate vp, vs, rho + idx = np.searchsorted(cum_time[a[indices[ind]], b[indices[ind]], :], time_sample[ind, k], side='left') + idx = idx if idx < len(cum_time[a[indices[ind]], b[indices[ind]], :]) else len(cum_time[a[indices[ind]], b[indices[ind]], :]) - 1 + vp_sample[ind, k] = vp[a[indices[ind]], b[indices[ind]], idx] + vs_sample[ind, k] = vs[a[indices[ind]], b[indices[ind]], idx] + rho_sample[ind, k] = rho[a[indices[ind]], b[indices[ind]], idx] + + + # Ricker wavelet + wavelet, t_axis, wav_center = ricker(np.arange(0, self.avo_config['wave_len'], dt), + f0=self.avo_config['frequency']) + + # Travel time corresponds to reflectivity series + t = time_sample[:, 0:-1] + + # interpolation time + t_interp = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], self.avo_config['t_sampling']) + trace_interp = np.zeros((len(indices), len(t_interp))) + + # number of pp reflection coefficients in the vertical direction + nz_rpp = vp_sample.shape[1] - 1 + conv_op = Convolve1D(nz_rpp, h=wavelet, offset=wav_center, dtype="float32") + + avo_data = [] + Rpp = [] + for i in range(len(self.avo_config['angle'])): + angle = self.avo_config['angle'][i] + Rpp = self.pp_func(vp_sample[:, :-1], vs_sample[:, :-1], rho_sample[:, :-1], + vp_sample[:, 1:], vs_sample[:, 1:], rho_sample[:, 1:], angle) + + + + for ind in range(len(indices)): + # convolution with the Ricker wavelet + + w_trace = conv_op * Rpp[ind, :] + + # Sample the trace into regular time interval + f = interp1d(np.squeeze(t[ind, :]), np.squeeze(w_trace), + kind='nearest', fill_value='extrapolate') + trace_interp[ind, :] = f(t_interp) + + if i == 0: + avo_data = trace_interp # 3D + elif i == 1: + avo_data = np.stack((avo_data, trace_interp), axis=-1) # 4D + else: + avo_data = np.concatenate((avo_data, trace_interp[:, :, :, np.newaxis]), axis=3) # 4D + + self.avo_data = avo_data + self.Rpp = Rpp + self.vp_sample = vp_sample + self.vs_sample = vs_sample + self.rho_sample = rho_sample + + + @classmethod + def _reformat3D_then_flatten(cls, array, flatten=True, order="F"): + """ + XILU: Quantities read by "EclipseData.cell_data" are put in the axis order of [nz, ny, nx]. To be consisent with + ECLIPSE/OPM custom, we need to change the axis order. We further flatten the array according to the specified order + """ + array = np.array(array) + if len(array.shape) != 1: # if array is a 1D array, then do nothing + assert isinstance(array, np.ndarray) and len(array.shape) == 3, "Only 3D numpy array are supported" + + # axis [0 (nz), 1 (ny), 2 (nx)] -> [2 (nx), 1 (ny), 0 (nz)] + new_array = np.transpose(array, axes=[2, 1, 0]) + if flatten: + new_array = new_array.flatten(order=order) + + return new_array + else: + return array + +class flow_grav(flow_rock, mixIn_multi_data): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): + super().__init__(input_dict, filename, options) + + self.grav_input = {} + assert 'grav' in input_dict, 'To do GRAV simulation, please specify an "GRAV" section in the "FWDSIM" part' + self._get_grav_info() + + def setup_fwd_run(self, **kwargs): + self.__dict__.update(kwargs) + + super().setup_fwd_run(redund_sim=None) + + def run_fwd_sim(self, state, member_i, del_folder=True): + # The inherited simulator also has a run_fwd_sim. Call this. + self.ensemble_member = member_i + #return super().run_fwd_sim(state, member_i, del_folder=del_folder) + self.pred_data = super().run_fwd_sim(state, member_i, del_folder) + return self.pred_data + + def call_sim(self, folder=None, wait_for_proc=False, save_folder=None): + # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. + # Then, get the pem. + if folder is None: + folder = self.folder + + if not self.no_flow: + # call call_sim in flow class (skip flow_rock, go directly to flow which is a parent of flow_rock) + success = super(flow_rock, self).call_sim(folder, True) + else: + success = True + # + # use output from flow simulator to forward model gravity response + if success: + self.get_grav_result(folder, save_folder) + + return success + + def get_grav_result(self, folder, save_folder): + if self.no_flow: + grid_file = self.pem_input['grid'] + grid = np.load(grid_file) + else: + self.ecl_case = ecl.EclipseCase(folder + os.sep + self.file + '.DATA') if folder[-1] != os.sep \ + else ecl.EclipseCase(folder + self.file + '.DATA') + grid = self.ecl_case.grid() + + + self.dyn_var = [] + + # cell centers + #cell_centre = self.find_cell_centre(grid) + + # receiver locations + # Make a mesh of the area + pad = self.grav_config.get('padding', 1500) # 3 km padding around the reservoir + if 'padding' not in self.grav_config: + print('Please specify extent of measurement locations (Padding in input file), using 1.5 km as default') + dxy = self.grav_config.get('grid_spacing', 1500) # + if 'grid_spacing' not in self.grav_config: + print('Please specify grid spacing in input file, using 1.5 km as default') + if 'seabed' in self.grav_config and self.grav_config['seabed'] is not None: + file_path = self.grav_config['seabed'] + water_depth = self.get_seabed_depths(file_path) + else: + water_depth = self.grav_config.get('water_depth', 300) + if 'water_depth' not in self.grav_config: + print('Please specify water depths in input file, using 300 m as default') + pos = self.measurement_locations(grid, water_depth, pad, dxy) + + # loop over vintages with gravity acquisitions + grav_struct = {} + + if 'baseline' in self.grav_config: # 4D measurement + base_time = dt.datetime(self.startDate['year'], self.startDate['month'], + self.startDate['day']) + dt.timedelta(days=self.grav_config['baseline']) + # porosity, saturation, densities, and fluid mass at time of baseline survey + grav_base = self.calc_mass(base_time, 0) + + + else: + # seafloor gravity only works in 4D mode + grav_base = None + print('Need to specify Baseline survey for gravity in input file') + + for v, assim_time in enumerate(self.grav_config['vintage']): + time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ + dt.timedelta(days=assim_time) + + # porosity, saturation, densities, and fluid mass at individual time-steps + grav_struct[v] = self.calc_mass(time, v+1) # calculate the mass of each fluid in each grid cell + + + + vintage = [] + + + for v, assim_time in enumerate(self.grav_config['vintage']): + dg = self.calc_grav(grid, grav_base, grav_struct[v], pos) + vintage.append(deepcopy(dg)) + + #save_dic = {'grav': dg, **self.grav_config} + save_dic = { + 'grav': dg, 'P_vint': grav_struct[v]['PRESSURE'], 'rho_gas_vint':grav_struct[v]['GAS_DEN'], + 'meas_location': pos, **self.grav_config, + **{key: grav_struct[v][key] - grav_base[key] for key in grav_struct[v].keys()} + } + if save_folder is not None: + file_name = save_folder + os.sep + f"grav_vint{v}.npz" if save_folder[-1] != os.sep \ + else save_folder + f"grav_vint{v}.npz" + else: + file_name = folder + os.sep + f"grav_vint{v}.npz" if folder[-1] != os.sep \ + else folder + f"grav_vint{v}.npz" + prior_folder = 'Prior_ensemble_results' + try: + files = os.listdir(prior_folder) + filename_to_check = f"grav_vint{v}_{folder}.npz" + + if filename_to_check in files: + file_name_rec = 'Ensemble_results/' + f"grav_vint{v}_{folder}.npz" if folder[-1] != os.sep \ + else 'Ensemble_results/' + f"grav_vint{v}_{folder[:-1]}.npz" + else: + file_name_rec = 'Prior_ensemble_results/' + f"grav_vint{v}_{folder}.npz" if folder[-1] != os.sep \ + else 'Prior_ensemble_results/' + f"grav_vint{v}_{folder[:-1]}.npz" + + except: + file_name_rec = 'Ensemble_results/' + f"grav_vint{v}_{folder}.npz" if folder[-1] != os.sep \ + else 'Ensemble_results/' + f"grav_vint{v}_{folder[:-1]}.npz" + np.savez(file_name_rec, **save_dic) + + np.savez(file_name, **save_dic) + + + # 4D response + self.grav_result = [] + for i, elem in enumerate(vintage): + self.grav_result.append(elem) + + def calc_mass(self, time, time_index = None): + + if self.no_flow: + time_input = time_index + else: + time_input = time + + # fluid phases given as input + phases = str.upper(self.pem_input['phases']) + phases = phases.split() + # + grav_input = {} + tmp_dyn_var = {} + + + tmp = self._get_pem_input('RPORV', time_input) + grav_input['RPORV'] = np.array(tmp[~tmp.mask], dtype=float) + + tmp = self._get_pem_input('PRESSURE', time_input) + #if time_input == time_index and time_index > 0: # to be activiated in case on inverts for Delta Pressure + # # Inverts for changes in dynamic variables using time-lapse data + # tmp_baseline = self._get_pem_input('PRESSURE', 0) + # tmp = tmp + tmp_baseline + grav_input['PRESSURE'] = np.array(tmp[~tmp.mask], dtype=float) + # convert pressure from Bar to MPa + if 'press_conv' in self.pem_input and time_input == time: + grav_input['PRESSURE'] = grav_input['PRESSURE'] * self.pem_input['press_conv'] + #else: + # print('Keyword RPORV missing from simulation output, need updated pore volumes at each assimilation step') + # extract saturation + if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: # This should be extended + for var in phases: + if var in ['WAT', 'GAS']: + tmp = self._get_pem_input('S{}'.format(var), time_input) + #if time_input == time_index and time_index > 0: # to be activated in case on inverts for Delta S + # # Inverts for changes in dynamic variables using time-lapse data + # tmp_baseline = self._get_pem_input('S{}'.format(var), 0) + # tmp = tmp + tmp_baseline + #tmp = self.ecl_case.cell_data('S{}'.format(var), time) + grav_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) + grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] > 1] = 1 + grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] < 0] = 0 + + grav_input['SOIL'] = 1 - (grav_input['SWAT'] + grav_input['SGAS']) + grav_input['SOIL'][grav_input['SOIL'] > 1] = 1 + grav_input['SOIL'][grav_input['SOIL'] < 0] = 0 + + + tmp_dyn_var['SWAT'] = grav_input['SWAT'] # = {f'S{ph}': saturations[i] for i, ph in enumerate(phases)} + tmp_dyn_var['SGAS'] = grav_input['SGAS'] + tmp_dyn_var['SOIL'] = grav_input['SOIL'] + + + elif 'WAT' in phases and 'GAS' in phases: # Smeaheia model + for var in phases: + if var in ['GAS']: + tmp = self._get_pem_input('S{}'.format(var), time_input) + #if time_input == time_index and time_index > 0: # to be activated in case on inverts for Delta S + # Inverts for changes in dynamic variables using time-lapse data + # tmp_baseline = self._get_pem_input('S{}'.format(var), 0) + # tmp = tmp + tmp_baseline + #tmp = self.ecl_case.cell_data('S{}'.format(var), time) + grav_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) + grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] > 1] = 1 + grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] < 0] = 0 + + grav_input['SWAT'] = 1 - (grav_input['SGAS']) + + # fluid saturation + tmp_dyn_var['SWAT'] = grav_input['SWAT'] #= {f'S{ph}': saturations[i] for i, ph in enumerate(phases)} + tmp_dyn_var['SGAS'] = grav_input['SGAS'] + + elif 'OIL' in phases and 'GAS' in phases: # Original Smeaheia model + for var in phases: + if var in ['GAS']: + tmp = self._get_pem_input('S{}'.format(var), time_input) + #if time_input == time_index and time_index > 0: # to be activated in case on inverts for Delta S + # Inverts for changes in dynamic variables using time-lapse data + # tmp_baseline = self._get_pem_input('S{}'.format(var), 0) + # tmp = tmp + tmp_baseline + #tmp = self.ecl_case.cell_data('S{}'.format(var), time) + grav_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) + grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] > 1] = 1 + grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] < 0] = 0 + + grav_input['SOIL'] = 1 - (grav_input['SGAS']) + + # fluid saturation + tmp_dyn_var['SOIL'] = grav_input['SOIL'] #= {f'S{ph}': saturations[i] for i, ph in enumerate(phases)} + tmp_dyn_var['SGAS'] = grav_input['SGAS'] + + else: + print('Type and number of fluids are unspecified in calc_mass') + + + # fluid densities + for var in phases: + dens = var + '_DEN' + #tmp = self.ecl_case.cell_data(dens, time) + if self.no_flow: + if any('pressure' in key for key in self.state.keys()): + if 'press_conv' in self.pem_input: + conv2pa = 1e6 #MPa to Pa + else: + conv2pa = 1e5 # Bar to Pa + + if var == 'GAS': + if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: + tmp = PropsSI('D', 'T', 298.15, 'P', grav_input['PRESSURE']*conv2pa, 'Methane') + elif 'WAT' in phases and 'GAS' in grav_input['PRESSURE']: # Smeaheia model T = 37 C + tmp = PropsSI('D', 'T', 310.15, 'P', grav_input['PRESSURE']*conv2pa, 'CO2') + mask = np.zeros(tmp.shape, dtype=bool) + tmp = np.ma.array(data=tmp, dtype=tmp.dtype, mask=mask) + elif var == 'WAT': + tmp = PropsSI('D', 'T|liquid', 298.15, 'P', grav_input['PRESSURE']*conv2pa, 'Water') + mask = np.zeros(tmp.shape, dtype=bool) + tmp = np.ma.array(data=tmp, dtype=tmp.dtype, mask=mask) + else: + tmp = self._get_pem_input(dens, time_input) + else: + tmp = self._get_pem_input(dens, time_input) + grav_input[dens] = np.array(tmp[~tmp.mask], dtype=float) + tmp_dyn_var[dens] = grav_input[dens] + else: + tmp = self._get_pem_input(dens, time_input) + grav_input[dens] = np.array(tmp[~tmp.mask], dtype=float) + tmp_dyn_var[dens] = grav_input[dens] + + + tmp_dyn_var['PRESSURE'] = grav_input['PRESSURE'] + tmp_dyn_var['RPORV'] = grav_input['RPORV'] + self.dyn_var.extend([tmp_dyn_var]) + + #fluid masses + for var in phases: + mass = var + '_mass' + grav_input[mass] = grav_input[var + '_DEN'] * grav_input['S' + var] * grav_input['RPORV'] + + return grav_input + + def calc_grav(self, grid, grav_base, grav_repeat, pos): + + cell_centre = self.find_cell_centre(grid) + x = cell_centre[0] + y = cell_centre[1] + z = cell_centre[2] + + # Initialize dg as a zero array, with shape depending on the condition + # assumes the length of each vector gives the total number of measurement points + n_meas = (len(pos['x'])) + dg = np.zeros(n_meas) # 1D array for dg + dg[:] = np.nan + + # fluid phases given as input + phases = str.upper(self.pem_input['phases']) + phases = phases.split() + if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: + dm = grav_repeat['OIL_mass'] + grav_repeat['WAT_mass'] + grav_repeat['GAS_mass'] - (grav_base['OIL_mass'] + grav_base['WAT_mass'] + grav_base['GAS_mass']) + + elif 'OIL' in phases and 'GAS' in phases: # Original Smeaheia model + dm = grav_repeat['OIL_mass'] + grav_repeat['GAS_mass'] - (grav_base['OIL_mass'] + grav_base['GAS_mass']) + # dm = grav_repeat['WAT_mass'] + grav_repeat['GAS_mass'] - (grav_base['WAT_mass'] + grav_base['GAS_mass']) + + elif 'WAT' in phases and 'GAS' in phases: # Smeaheia model + dm = grav_repeat['WAT_mass'] + grav_repeat['GAS_mass'] - (grav_base['WAT_mass'] + grav_base['GAS_mass']) + #dm = grav_repeat['WAT_mass'] + grav_repeat['GAS_mass'] - (grav_base['WAT_mass'] + grav_base['GAS_mass']) + + else: + dm = None + print('Type and number of fluids are unspecified in calc_grav') + + + for j in range(n_meas): + + # Calculate dg for the current measurement location (j, i) + dg_tmp = (z - pos['z'][j]) / ((x - pos['x'][j]) ** 2 + (y - pos['y'][j]) ** 2 + ( + z - pos['z'][j]) ** 2) ** (3 / 2) + + dg[j] = np.dot(dg_tmp, dm) + #print(f'Progress: {j + 1}/{n_meas}') # Mimicking wait bar + + # Scale dg by the constant + dg *= 6.67e-3 + + return dg + + def _get_grav_info(self, grav_config=None): + """ + GRAV configuration + """ + # list of configuration parameters in the "Grav" section of teh pipt file + config_para_list = ['baseline', 'vintage', 'water_depth', 'padding', 'grid_spacing', 'seabed'] + + if 'grav' in self.input_dict: + self.grav_config = {} + for elem in self.input_dict['grav']: + assert elem[0] in config_para_list, f'Property {elem[0]} not supported' + if elem[0] == 'vintage' and not isinstance(elem[1], list): + elem[1] = [elem[1]] + self.grav_config[elem[0]] = elem[1] + else: + self.grav_config = None + + def extract_data(self, member): + # start by getting the data from the flow simulator + super(flow_rock, self).extract_data(member) + + # get the gravity data from results + for prim_ind in self.l_prim: + # Loop over all keys in pred_data (all data types) + for key in self.all_data_types: + if 'grav' in key: + if self.true_prim[1][prim_ind] in self.grav_config['vintage']: + v = self.grav_config['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = self.grav_result[v].flatten() + +class flow_seafloor_disp(flow_rock, mixIn_multi_data): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): + super().__init__(input_dict, filename, options) + + assert 'sea_disp' in input_dict, 'To do subsidence/uplift simulation, please specify an "SEA_DISP" section in the pipt file' + self._get_disp_info() + + + def setup_fwd_run(self, **kwargs): + self.__dict__.update(kwargs) + + super().setup_fwd_run(redund_sim=None) + + def run_fwd_sim(self, state, member_i, del_folder=True): + # The inherited simulator also has a run_fwd_sim. Call this. + self.ensemble_member = member_i + self.pred_data = super().run_fwd_sim(state, member_i, del_folder) + + return self.pred_data + + def call_sim(self, folder=None, wait_for_proc=False, save_folder=None): + # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. + # Then, get the pem. + if folder is None: + folder = self.folder + + # run flow simulator + # success = True + success = super(flow_rock, self).call_sim(folder, True) + + # use output from flow simulator to forward model gravity response + if success: + # calculate gravity data based on flow simulation output + self.get_displacement_result(folder, save_folder) + + + return success + + def get_displacement_result(self, folder, save_folder): + if self.no_flow: + grid_file = self.pem_input['grid'] + grid = np.load(grid_file) + else: + self.ecl_case = ecl.EclipseCase(folder + os.sep + self.file + '.DATA') if folder[-1] != os.sep \ + else ecl.EclipseCase(folder + self.file + '.DATA') + grid = self.ecl_case.grid() + + self.dyn_var = [] + + # receiver locations + pad = self.disp_config.get('padding', 1500) # 3 km padding around the reservoir + if 'padding' not in self.disp_config: + print('Please specify extent of measurement locations, padding in input file, using 1.5 km as default') + dxy = self.disp_config.get('grid_spacing', 1500) # + if 'grid_spacing' not in self.disp_config: + print('Please specify grid spacing in input file, using 1.5 km as default') + if 'seabed' in self.disp_config and self.disp_config['seabed'] is not None: + file_path = self.disp_config['seabed'] + water_depth = self.get_seabed_depths(file_path) + else: + water_depth = self.disp_config.get('water_depth', 300) + if 'water_depth' not in self.disp_config: + print('Please specify water depths in input file, using 300 m as default') + pos = self.measurement_locations(grid, water_depth, pad, dxy) + + # loop over vintages with gravity acquisitions + disp_struct = {} + + if 'baseline' in self.disp_config: # 4D measurement + base_time = dt.datetime(self.startDate['year'], self.startDate['month'], + self.startDate['day']) + dt.timedelta(days=self.disp_config['baseline']) + # pore volume at time of baseline survey + disp_base = self.get_pore_volume(base_time, 0) + + else: + # seafloor displacement only work in 4D mode + disp_base = None + print('Need to specify Baseline survey for displacement modelling in input file') + + for v, assim_time in enumerate(self.disp_config['vintage']): + time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ + dt.timedelta(days=assim_time) + + # pore volume and pressure at individual time-steps + disp_struct[v] = self.get_pore_volume(time, v+1) # calculate the mass of each fluid in each grid cell + + vintage = [] + + for v, assim_time in enumerate(self.disp_config['vintage']): + # calculate subsidence and uplift + dz_seafloor = self.map_z_response(disp_base, disp_struct[v], grid, pos) + vintage.append(deepcopy(dz_seafloor)) + + save_dic = {'sea_disp': dz_seafloor, 'meas_location': pos, **self.disp_config} + if save_folder is not None: + file_name = save_folder + os.sep + f"sea_disp_vint{v}.npz" if save_folder[-1] != os.sep \ + else save_folder + f"sea_disp_vint{v}.npz" + else: + file_name = folder + os.sep + f"sea_disp_vint{v}.npz" if folder[-1] != os.sep \ + else folder + f"sea_disp_vint{v}.npz" + file_name_rec = 'Ensemble_results/' + f"sea_disp_vint{v}_{folder}.npz" if folder[-1] != os.sep \ + else 'Ensemble_results/' + f"sea_disp_vint{v}_{folder[:-1]}.npz" + np.savez(file_name_rec, **save_dic) + + np.savez(file_name, **save_dic) + + + # 4D response + self.disp_result = [] + for i, elem in enumerate(vintage): + self.disp_result.append(elem) + + def get_pore_volume(self, time, time_index = None): + + if self.no_flow: + time_input = time_index + else: + time_input = time + + # + disp_input = {} + tmp_dyn_var = {} + + + tmp = self._get_pem_input('RPORV', time_input) + disp_input['RPORV'] = np.array(tmp[~tmp.mask], dtype=float) + + tmp = self._get_pem_input('PRESSURE', time_input) + #if time_input == time_index and time_index > 0: # to be activiated in case on inverts for Delta Pressure + # # Inverts for changes in dynamic variables using time-lapse data + # tmp_baseline = self._get_pem_input('PRESSURE', 0) + # tmp = tmp + tmp_baseline + disp_input['PRESSURE'] = np.array(tmp[~tmp.mask], dtype=float) + # convert pressure from Bar to MPa + if 'press_conv' in self.pem_input and time_input == time: + disp_input['PRESSURE'] = disp_input['PRESSURE'] * self.pem_input['press_conv'] + #else: + # print('Keyword RPORV missing from simulation output, need pdated porevolumes at each assimilation step') + + + tmp_dyn_var['PRESSURE'] = disp_input['PRESSURE'] + tmp_dyn_var['RPORV'] = disp_input['RPORV'] + self.dyn_var.extend([tmp_dyn_var]) + + return disp_input + + def compute_horizontal_distance(self, pos, x, y): + dx = pos['x'][:, np.newaxis] - x + dy = pos['y'][:, np.newaxis] - y + rho = np.sqrt(dx ** 2 + dy ** 2).flatten() + return rho + + def map_z_response(self, base, repeat, grid, pos): + """ + Maps out subsidence and uplift based either on the simulation + model pressure drop (method = 'pressure') or simulated change in pore volume + using either the van Opstal or Geertsma forward model + + Arguments: + base -- A dictionary containing baseline pressures and pore volumes. + repeat -- A dictionary containing pressures and pore volumes at repeat measurements. + + compute subsidence at position 'pos', b + + Output is modeled subsidence in cm. + + """ + + # Method to compute pore volume change + method = self.disp_config['method'].lower() + + # Forward model to compute subsidence/uplift response + model = self.disp_config['model'].lower() + + if self.disp_config['poisson'] > 0.5: + poisson = 0.5 + print('Poisson\'s ratio exceeds physical limits, setting it to 0.5') + else: + poisson = self.disp_config['poisson'] + + # Depth of rigid basement + z_base = self.disp_config['z_base'] + + compressibility = self.disp_config['compressibility'] # 1/MPa + + E = ((1 + poisson) * (1 - 2 * poisson)) / ((1 - poisson) * compressibility) + + # coordinates of cell centres + cell_centre = self.find_cell_centre(grid) + + # compute pore volume change between baseline and repeat survey + # based on the reservoir pore volumes in the individual vintages + if method == 'pressure': + dV = base['RPORV'] * (base['PRESSURE'] - repeat['PRESSURE']) * compressibility + else: + dV = base['RPORV'] - repeat['RPORV'] + + # coordinates of active cell centres + x = cell_centre[0] + y = cell_centre[1] + z = cell_centre[2] + + + if model == 'van_Opstal': + # Represents a signal change for subsidence/uplift. + #trans_func = t_van_opstal + component = ["Geertsma_vertical", "System_3_vertical"] + else: # Use Geertsma + #trans_func = t_geertsma + component = ["Geertsma_vertical"] + # Initialization + dz_1_2 = 0 + dz_3 = 0 + + # indices of active cells: + true_indices = np.where(grid['ACTNUM']) + # number of active gridcells + n_nucleous = len(true_indices[0]) + assert n_nucleous == len(x) + #pr = cProfile.Profile() + #pr.enable() + + for j in range(n_nucleous): + rho = self.compute_horizontal_distance(pos, x[j], y[j]) + THH, TRB = self.compute_deformation_transfer(pos['z'], z[j], z_base, rho, poisson, E, dV[j], component) + dz_1_2 = dz_1_2 + THH + dz_3 = dz_3 + TRB + + #pr.disable() + + # Print profiling results + #stats = pstats.Stats(pr) + #stats.strip_dirs() + #stats.sort_stats('cumulative') + #stats.print_stats() + + if model == 'van_Opstal': + # Represents a signal change for subsidence/uplift. + dz = dz_1_2 + dz_3 + else: # Use Geertsma + dz = dz_1_2 + + # Convert from meters to centimeters + dz *= 100 + + return dz + + def compute_van_opstal_transfer_function(self, z_res, z_base, rho, poisson): + """ + Compute the Van Opstal transfer function. + + Args: + z_res -- Numpy array of depths to reservoir cells [m]. + z_base -- Distance to the basement [m]. + rho -- Numpy array of horizontal distances in the field [m]. + poisson -- Poisson's ratio. + + Returns: + T -- Numpy array of the transfer function values. + T_geertsma -- Numpy array of the Geertsma transfer function values. + """ + + # Change to km scale + rho = rho / 1e3 + z_res = z_res / 1e3 + z_base = z_base / 1e3 + + + # Find lambda max (to optimize Hilbert transform) + cutoff = 1e-10 # Function value at max lambda + try: + lambda_max = fsolve(lambda x: 4 * (2 * x * z_base + 1) / (3 - 4 * poisson) * np.exp( + x * (np.max(z_res) - 2 * z_base)) - cutoff, 10)[0] + except: + lambda_max = 10 # Default value if unable to solve for max lambda + + lambda_vals = np.linspace(0, lambda_max, 100) + # range of lateral distances between measurement location and reservoir cells + nj = len(rho) + # range of vertical distances between measurement location and reservoir cells + ni = len(z_res) + # initialize + t_van_opstal = np.zeros((ni, nj)) + + # input function to make a hankel transform of order 0 of + c_t = self.van_opstal(lambda_vals, z_res[0], z_base, poisson) + + h_t, i_t = self.h_t(c_t, lambda_vals, rho) # Extract integrand + t_van_opstal[0, :] = (2 * z_res[0] / (rho ** 2 + z_res[0] ** 2) ** (3 / 2)) + h_t / (2 * np.pi) + + for i in range(1, ni): + C = self.van_opstal(lambda_vals, z_res[i], z_base, poisson) + h_t = self.h_t(C, lambda_vals, rho, i_t) + t_van_opstal[i, :] = (2 * z_res[i] / (rho ** 2 + z_res[i] ** 2) ** (3 / 2)) + h_t / (2 * np.pi) + + t_van_opstal *= 1e-6 # Convert back to meters + + t_geertsma = (2 * z_res[:, np.newaxis] / ((np.ones((ni, 1)) * rho) ** 2 + (z_res[:, np.newaxis]) ** 2) ** ( + 3 / 2)) * 1e-6 + + return t_van_opstal, t_geertsma - self.scale = [] - self.pem_result = [] - self.bar_result = [] + def van_opstal(self, lambda_vals, z_res, z_base, poisson): + """ + Compute the Van Opstal transfer function. - def _getpeminfo(self, input_dict): + Args: + lambda_vals -- Numpy array of lambda values. + z_res -- Depth to reservoir [m]. + z_base -- Distance to the basement [m]. + poisson -- Poisson's ratio. + + Returns: + value -- Numpy array of computed values. """ - Get, and return, flow and PEM modules + + term1 = np.exp(lambda_vals * z_res) * (2 * lambda_vals * z_base + 1) + term2 = np.exp(-lambda_vals * z_res) * ( + 4 * lambda_vals ** 2 * z_base ** 2 + 2 * lambda_vals * z_base + (3 - 4 * poisson) ** 2) + + term3_numer = (3 - 4 * poisson) * ( + np.exp(-lambda_vals * (2 * z_base + z_res)) - np.exp(-lambda_vals * (2 * z_base - z_res))) + term3_denom = 2 * ((1 - 2 * poisson) ** 2 + lambda_vals ** 2 * z_base ** 2 + (3 - 4 * poisson) * np.cosh( + lambda_vals * z_base) ** 2) + + value = term1 - term2 - (term3_numer / term3_denom) + + return value + + def hankel_transform_order_0(self, f, r_max, num_points=1000): """ - if 'pem' in input_dict: - self.pem_input = {} - for elem in input_dict['pem']: - if elem[0] == 'model': # Set the petro-elastic model - self.pem_input['model'] = elem[1] - if elem[0] == 'depth': # provide the npz of depth values - self.pem_input['depth'] = elem[1] - if elem[0] == 'actnum': # the npz of actnum values - self.pem_input['actnum'] = elem[1] - if elem[0] == 'baseline': # the time for the baseline 4D measurment - self.pem_input['baseline'] = elem[1] - if elem[0] == 'vintage': - self.pem_input['vintage'] = elem[1] - if not type(self.pem_input['vintage']) == list: - self.pem_input['vintage'] = [elem[1]] - if elem[0] == 'ntg': - self.pem_input['ntg'] = elem[1] - if elem[0] == 'press_conv': - self.pem_input['press_conv'] = elem[1] - if elem[0] == 'compaction': - self.pem_input['compaction'] = True - if elem[0] == 'overburden': # the npz of overburden values - self.pem_input['overburden'] = elem[1] - if elem[0] == 'percentile': # use for scaling - self.pem_input['percentile'] = elem[1] - if elem[0] == 'clusters': # number of clusters for each barycenter - self.pem_input['clusters'] = elem[1] + Computes the Hankel transform of order 0 of a function f(r). - pem = getattr(import_module('simulator.rockphysics.' + - self.pem_input['model'].split()[0]), self.pem_input['model'].split()[1]) + Parameters: + - f: callable, the function to transform, f(r) + - r_max: float, upper limit of the integral (approximate infinity) + - num_points: int, number of points for numerical integration - self.pem = pem(self.pem_input) + Returns: + - k_values: array of k values + - H_k: array of Hankel transform evaluated at k_values + """ + r = np.linspace(0, r_max, num_points) + dr = r[1] - r[0] + f_r = f(r) + + def integrand(r, k): + return f(r) * j0(k * r) * r + + # Define a range of k values to evaluate + k_min, k_max = 0, 10 # adjust as needed + k_values = np.linspace(k_min, k_max, 100) + + H_k = [] + + for k in k_values: + # Perform numerical integration over r + result, _ = quad(integrand, 0, r_max, args=(k,)) + H_k.append(result) + + return k_values, np.array(H_k) + + def makeL(self, poisson, k, c, A_g, eps, lambda_): + L = A_g * ( + (4 * poisson - 3 + 2 * k * lambda_) * np.exp(-lambda_ * (k + c)) + - np.exp(lambda_ * eps * (k - c)) + ) + return L + + def makeM(self, poisson, k, c, A_g, eps, lambda_): + M = A_g * ( + (4 * poisson - 3 - 2 * k * lambda_) * np.exp(-lambda_ * (k + c)) + - eps * np.exp(lambda_ * eps * (k - c)) + ) + return M + + def makeDelta(self, poisson, k, lambda_): + Delta = ( + (4 * poisson - 3) * np.cosh(k * lambda_) ** 2 + - (k * lambda_) ** 2 + - (1 - 2 * poisson) ** 2 + ) + return Delta + + def makeB(self, poisson, k, c, A_g, eps, lambda_): + L = self.makeL(poisson, k, c, A_g, eps, lambda_) + M = self.makeM(poisson, k, c, A_g, eps, lambda_) + Delta = self.makeDelta(poisson, k, lambda_) + + numerator = ( + lambda_ * L * (2 * (1 - poisson) * np.cosh(k * lambda_) - lambda_ * k * np.sinh(k * lambda_)) + + lambda_ * M * ((1 - 2 * poisson) * np.sinh(k * lambda_) + k * lambda_ * np.cosh(k * lambda_)) + ) + + B = numerator / Delta + return B + + def makeC(self,poisson, k, c, A_g, eps, lambda_): + L = self.makeL(poisson, k, c, A_g, eps, lambda_) + M = self.makeM(poisson, k, c, A_g, eps, lambda_) + Delta = self.makeDelta(poisson, k, lambda_) + + numerator = ( + lambda_ * L * ((1 - 2 * poisson) * np.sinh(k * lambda_) - lambda_ * k * np.cosh(k * lambda_)) + + lambda_ * M * (2 * (1 - poisson) * np.cosh(k * lambda_) + k * lambda_ * np.sinh(k * lambda_)) + ) + + C = numerator / Delta + return C + + def uHH_integrand(self, lambda_, z, rho, eps, c, poisson): + val = lambda_ * (eps * np.exp(lambda_ * eps * (z - c)) + + (3 - 4 * poisson + 2 * z * lambda_) * + np.exp(-lambda_ * (z + c))) + return val * j0(lambda_ * rho) + + def compute_deformation_transfer(self, z, c, k, rho, poisson, E, dV, component): + scale = 1000 # convert to km + # depth of receiver positions + z = np.array(z)/scale + rho = np.array(rho)/scale + # depth of reservoir cell + c = c/scale + k = k/scale + component = list(component) + # number of measurement locations + n_rec = len(z) + assert len(rho) == n_rec + THH = np.zeros(n_rec) + TRB = np.zeros(n_rec) + + # Constants + A_g = -dV * E / (4 * np.pi * (1 + poisson)) + uHH_outside_intregral = -(A_g * (1 + poisson)) / E + uRB_outside_intregral = (1 + poisson) / E + + for c_n in component: + if c_n == 'Geertsma_vertical': + lambda_max = 15 / np.max(rho).item() + for i in range(n_rec): + if rho[i] > np.abs(c-z[i])*3: + THH[i] = 0 + else: + eps = np.sign(c - z[i]) + THH[i] = quad(lambda lambda_var: self.uHH_integrand(lambda_var, z[i], rho[i], eps, c, poisson), 0, lambda_max)[0] * uHH_outside_intregral + THH[i] = THH[i]*scale**-2 + + elif c_n == 'System_3_vertical': + lambda_max = 30 / np.max(rho).item() + num_points = 500 + lambda_grid = np.linspace(0, lambda_max, num_points) + + sinh_z = np.sinh(z[:, np.newaxis] * lambda_grid) + cosh_z = np.cosh(z[:, np.newaxis] * lambda_grid) + J0_rho = j0(lambda_grid * rho) + + # + for i in range(n_rec): + if rho[i] > np.abs(c-z[i])*3: + TRB[i] = 0 + else: + z_i = z[i] + sinh_z_i = sinh_z[i] + cosh_z_i = cosh_z[i] + + b_values = self.makeB(poisson, k, c, A_g, -1, lambda_grid) + c_values = self.makeC(poisson, k, c, A_g, -1, lambda_grid) + + part1 = b_values * (lambda_grid * z_i * cosh_z_i - (1 - 2 * poisson) * sinh_z_i) + part2 = c_values * ((2 * (1 - poisson) * cosh_z_i) - lambda_grid * z_i * sinh_z_i) + + values = (part1 + part2) * J0_rho[:, i]#J0_rho_j + integral_result = np.trapz(values, lambda_grid) + TRB[i] = integral_result * uRB_outside_intregral + TRB[i] = TRB[i]*scale**-2 + + return THH, TRB + + def h_t(self, h, r=None, k=None, i_k=None): + """ + Hankel transform of order 0. + + Args: + h -- Signal h(r). + r -- Radial positions [m] (optional). + k -- Spatial frequencies [rad/m] (optional). + I -- Integration kernel (optional). + + Returns: + h_t -- Spectrum H(k). + I -- Integration kernel. + """ + + # Check if h is a vector + if h.ndim > 1: + raise ValueError('Signal must be a vector.') + + if r is None or len(r) == 0: + r = np.arange(len(h)) # Default to 0:numel(h)-1 else: - self.pem = None + r = np.sort(r) + h = h[np.argsort(r)] # Sort h according to sorted r + + if k is None or len(k) == 0: + k = np.pi / len(h) * np.arange(len(h)) # Default spatial frequencies + + if i_k is None: + # Create integration kernel I + r = np.concatenate([(r[:-1] + r[1:]) / 2, [r[-1]]]) # Midpoints plus last point + i_k = (2 * np.pi / k[:, np.newaxis]) * r * jv(1, k[:, np.newaxis] * r) # Bessel function + i_k[k == 0, :] = np.pi * r * r + i_k = i_k - np.hstack([np.zeros((len(k), 1)), i_k[:, :-1]]) # Shift integration kernel + else: + # Ensure I is sorted based on r + i_k = i_k[:, np.argsort(r)] - def setup_fwd_run(self, redund_sim): - super().setup_fwd_run(redund_sim=redund_sim) + # Compute Hankel Transform + h_t = np.reshape(i_k @ h.flatten(), k.shape) + + + + return h_t, i_k + + def _get_disp_info(self, disp_config=None): + """ + seafloor displacement (uplift/subsidence) configuration + """ + # list of configuration parameters in the "Grav" section of teh pipt file + config_para_list = ['baseline', 'vintage', 'method', 'model', 'poisson', 'compressibility', + 'z_base', 'grid_spacing', 'padding', 'seabed', 'water_depth'] + + if 'sea_disp' in self.input_dict: + self.disp_config = {} + for elem in self.input_dict['sea_disp']: + assert elem[0] in config_para_list, f'Property {elem[0]} not supported' + if elem[0] == 'vintage' and not isinstance(elem[1], list): + elem[1] = [elem[1]] + self.disp_config[elem[0]] = elem[1] + else: + self.disp_config = None + + def extract_data(self, member): + # start by getting the data from the flow simulator i.e. prod. and inj. data + super(flow_rock, self).extract_data(member) + + # get the gravity data from results + for prim_ind in self.l_prim: + # Loop over all keys in pred_data (all data types) + for key in self.all_data_types: + if 'sea_disp' in key: + if self.true_prim[1][prim_ind] in self.disp_config['vintage']: + v = self.disp_config['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = self.disp_result[v].flatten() + +class flow_grav_and_avo(flow_avo, flow_grav): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): + super().__init__(input_dict, filename, options) + + self.grav_input = {} + assert 'grav' in input_dict, 'To do GRAV simulation, please specify an "GRAV" section in the "FWDSIM" part' + self._get_grav_info() + + assert 'avo' in input_dict, 'To do AVO simulation, please specify an "AVO" section in the "FWDSIM" part' + self._get_avo_info() + + def setup_fwd_run(self, **kwargs): + self.__dict__.update(kwargs) + + super().setup_fwd_run(redund_sim=None) def run_fwd_sim(self, state, member_i, del_folder=True): # The inherited simulator also has a run_fwd_sim. Call this. self.ensemble_member = member_i - self.pred_data = super().run_fwd_sim(state, member_i, del_folder=True) + self.pred_data = super().run_fwd_sim(state, member_i, del_folder) return self.pred_data - def call_sim(self, folder=None, wait_for_proc=False): + def call_sim(self, folder=None, wait_for_proc=False, save_folder=None): # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. # Then, get the pem. - success = super().call_sim(folder, wait_for_proc) + if folder is None: + folder = self.folder + else: + self.folder = folder - if success: - self.ecl_case = ecl.EclipseCase( - 'En_' + str(self.ensemble_member) + os.sep + self.file + '.DATA') - phases = self.ecl_case.init.phases - #if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: # This should be extended - if 'WAT' in phases and 'GAS' in phases: - vintage = [] - # loop over seismic vintages - for v, assim_time in enumerate(self.pem_input['vintage']): - time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ - dt.timedelta(days=assim_time) - pem_input = {} - # get active porosity - tmp = self.ecl_case.cell_data('PORO') - if 'compaction' in self.pem_input: - multfactor = self.ecl_case.cell_data('PORV_RC', time) + # run flow simulator + # success = True + success = super(flow_rock, self).call_sim(folder, True) - pem_input['PORO'] = np.array( - multfactor[~tmp.mask]*tmp[~tmp.mask], dtype=float) - else: - pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float) - # get active NTG if needed - if 'ntg' in self.pem_input: - if self.pem_input['ntg'] == 'no': - pem_input['NTG'] = None - else: - tmp = self.ecl_case.cell_data('NTG') - pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) - else: - tmp = self.ecl_case.cell_data('NTG') - pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) + # use output from flow simulator to forward model gravity response + if success: + # calculate gravity data based on flow simulation output + self.get_grav_result(folder, save_folder) + # calculate avo data based on flow simulation output + self.get_avo_result(folder, save_folder) - pem_input['RS'] = None - for var in ['SWAT', 'SGAS', 'PRESSURE', 'RS']: - try: - tmp = self.ecl_case.cell_data(var, time) - except: - pass - # only active, and conv. to float - pem_input[var] = np.array(tmp[~tmp.mask], dtype=float) + return success - if 'press_conv' in self.pem_input: - pem_input['PRESSURE'] = pem_input['PRESSURE'] * \ - self.pem_input['press_conv'] - tmp = self.ecl_case.cell_data('PRESSURE', 1) - if hasattr(self.pem, 'p_init'): - P_init = self.pem.p_init*np.ones(tmp.shape)[~tmp.mask] - else: - # initial pressure is first - P_init = np.array(tmp[~tmp.mask], dtype=float) + def extract_data(self, member): + # start by getting the data from the flow simulator i.e. prod. and inj. data + super(flow_rock, self).extract_data(member) - if 'press_conv' in self.pem_input: - P_init = P_init*self.pem_input['press_conv'] + # get the gravity data from results + for prim_ind in self.l_prim: + # Loop over all keys in pred_data (all data types) + for key in self.all_data_types: + if 'grav' in key: + if self.true_prim[1][prim_ind] in self.grav_config['vintage']: + v = self.grav_config['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = self.grav_result[v].flatten() - saturations = [1 - (pem_input['SWAT'] + pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)] - for ph in phases] - # Get the pressure - self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'], - ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init, - ensembleMember=self.ensemble_member) - # mask the bulkimp to get proper dimensions - tmp_value = np.zeros(self.ecl_case.init.shape) - tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp - self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float, - mask=deepcopy(self.ecl_case.init.mask)) - # run filter - self.pem._filter() - vintage.append(deepcopy(self.pem.bulkimp)) + if 'avo' in key: + if self.true_prim[1][prim_ind] in self.pem_input['vintage']: + idx = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) + filename = self.folder + os.sep + key + '_vint' + str(idx) + '.npz' if self.folder[-1] != os.sep \ + else self.folder + key + '_vint' + str(idx) + '.npz' + with np.load(filename) as f: + self.pred_data[prim_ind][key] = f[key] + #v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) + #self.pred_data[prim_ind][key] = self.avo_result[v].flatten() + +class flow_grav_seafloor_disp_and_avo(flow_avo, flow_grav, flow_seafloor_disp): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): + super().__init__(input_dict, filename, options) - if hasattr(self.pem, 'baseline'): # 4D measurement - base_time = dt.datetime(self.startDate['year'], self.startDate['month'], - self.startDate['day']) + dt.timedelta(days=self.pem.baseline) - # pem_input = {} - # get active porosity - tmp = self.ecl_case.cell_data('PORO') + assert 'grav' in input_dict, 'To do GRAV simulation, please specify an "GRAV" section in the "FWDSIM" part' + self._get_grav_info() - if 'compaction' in self.pem_input: - multfactor = self.ecl_case.cell_data('PORV_RC', base_time) + assert 'avo' in input_dict, 'To do AVO simulation, please specify an "AVO" section in the "FWDSIM" part' + self._get_avo_info() - pem_input['PORO'] = np.array( - multfactor[~tmp.mask] * tmp[~tmp.mask], dtype=float) - else: - pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float) + assert 'sea_disp' in input_dict, 'To do subsidence/uplift simulation, please specify an "SEA_DISP" section in the "FWDSIM" part' + self._get_disp_info() - pem_input['RS'] = None - for var in ['SWAT', 'SGAS', 'PRESSURE', 'RS']: - try: - tmp = self.ecl_case.cell_data(var, base_time) - except: - pass - # only active, and conv. to float - pem_input[var] = np.array(tmp[~tmp.mask], dtype=float) + def setup_fwd_run(self, **kwargs): + self.__dict__.update(kwargs) - if 'press_conv' in self.pem_input: - pem_input['PRESSURE'] = pem_input['PRESSURE'] * \ - self.pem_input['press_conv'] + super().setup_fwd_run(redund_sim=None) - saturations = [1 - (pem_input['SWAT'] + pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)] - for ph in phases] - # Get the pressure - self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'], - ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init, - ensembleMember=None) + def run_fwd_sim(self, state, member_i, del_folder=True): + # The inherited simulator also has a run_fwd_sim. Call this. + self.ensemble_member = member_i + self.pred_data = super().run_fwd_sim(state, member_i, del_folder) - # mask the bulkimp to get proper dimensions - tmp_value = np.zeros(self.ecl_case.init.shape) + return self.pred_data - tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp - # kill if values are inf or nan - assert not np.isnan(tmp_value).any() - assert not np.isinf(tmp_value).any() - self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float, - mask=deepcopy(self.ecl_case.init.mask)) - self.pem._filter() + def call_sim(self, folder=None, wait_for_proc=False, save_folder=None): + # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. + # Then, get the pem. + if folder is None: + folder = self.folder + else: + self.folder = folder - # 4D response - for i, elem in enumerate(vintage): - self.pem_result.append(elem - deepcopy(self.pem.bulkimp)) - else: - for i, elem in enumerate(vintage): - self.pem_result.append(elem) + # run flow simulator + # success = True + success = super(flow_rock, self).call_sim(folder, True) - # Extract k-means centers and interias for each element in pem_result - if 'clusters' in self.pem_input: - npzfile = np.load(self.pem_input['clusters'], allow_pickle=True) - n_clusters_list = npzfile['n_clusters_list'] - npzfile.close() - else: - n_clusters_list = len(self.pem_result)*[2] - kmeans_kwargs = {"init": "random", "n_init": 10, "max_iter": 300, "random_state": 42} - for i, bulkimp in enumerate(self.pem_result): - std = np.std(bulkimp) - features = np.argwhere(np.squeeze(np.reshape(np.abs(bulkimp), self.ecl_case.init.shape,)) > 3 * std) - scaler = StandardScaler() - scaled_features = scaler.fit_transform(features) - kmeans = KMeans(n_clusters=n_clusters_list[i], **kmeans_kwargs) - kmeans.fit(scaled_features) - kmeans_center = np.squeeze(scaler.inverse_transform(kmeans.cluster_centers_)) # data / measurements - self.bar_result.append(np.append(kmeans_center, kmeans.inertia_)) + # use output from flow simulator to forward model gravity response + if success: + # calculate gravity data based on flow simulation output + self.get_grav_result(folder, save_folder) + # calculate avo data based on flow simulation output + self.get_avo_result(folder, save_folder) + # calculate gravity data based on flow simulation output + self.get_displacement_result(folder, save_folder) return success def extract_data(self, member): - # start by getting the data from the flow simulator - super().extract_data(member) + # start by getting the data from the flow simulator i.e. prod. and inj. data + super(flow_rock, self).extract_data(member) - # get the barycenters and inertias + # get the gravity data from results for prim_ind in self.l_prim: # Loop over all keys in pred_data (all data types) for key in self.all_data_types: - if key in ['barycenter']: + if 'grav' in key: + if self.true_prim[1][prim_ind] in self.grav_config['vintage']: + v = self.grav_config['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = self.grav_result[v].flatten() + + if 'avo' in key: if self.true_prim[1][prim_ind] in self.pem_input['vintage']: - v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) - self.pred_data[prim_ind][key] = self.bar_result[v].flatten() + idx = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) + filename = self.folder + os.sep + key + '_vint' + str(idx) + '.npz' if self.folder[-1] != os.sep \ + else self.folder + key + '_vint' + str(idx) + '.npz' + with np.load(filename) as f: + self.pred_data[prim_ind][key] = f[key] + + if 'sea_disp' in key: + if self.true_prim[1][prim_ind] in self.disp_config['vintage']: + v = self.disp_config['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = self.disp_result[v].flatten() + +class flow_grav_seafloor_disp(flow_grav, flow_seafloor_disp): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): + super().__init__(input_dict, filename, options) + + assert 'grav' in input_dict, 'To do GRAV simulation, please specify an "GRAV" section in the "FWDSIM" part' + self._get_grav_info() + + assert 'sea_disp' in input_dict, 'To do subsidence/uplift simulation, please specify an "SEA_DISP" section in the "FWDSIM" part' + self._get_disp_info() + + def setup_fwd_run(self, **kwargs): + self.__dict__.update(kwargs) + + super().setup_fwd_run(redund_sim=None) + + def run_fwd_sim(self, state, member_i, del_folder=True): + # The inherited simulator also has a run_fwd_sim. Call this. + self.ensemble_member = member_i + self.pred_data = super().run_fwd_sim(state, member_i, del_folder) + + return self.pred_data + + def call_sim(self, folder=None, wait_for_proc=False, save_folder=None): + # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. + # Then, get the pem. + if folder is None: + folder = self.folder + else: + self.folder = folder + + # run flow simulator + # success = True + success = super(flow_rock, self).call_sim(folder, True) + + # use output from flow simulator to forward model gravity response + if success: + # calculate gravity data based on flow simulation output + self.get_grav_result(folder, save_folder) + # calculate avo data based on flow simulation output + #self.get_avo_result(folder, save_folder) + # calculate gravity data based on flow simulation output + self.get_displacement_result(folder, save_folder) + + return success + + def extract_data(self, member): + # start by getting the data from the flow simulator i.e. prod. and inj. data + super(flow_rock, self).extract_data(member) + + # get the gravity data from results + for prim_ind in self.l_prim: + # Loop over all keys in pred_data (all data types) + for key in self.all_data_types: + if 'grav' in key: + if self.true_prim[1][prim_ind] in self.grav_config['vintage']: + v = self.grav_config['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = self.grav_result[v].flatten() + + if 'sea_disp' in key: + if self.true_prim[1][prim_ind] in self.disp_config['vintage']: + v = self.disp_config['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = self.disp_result[v].flatten() diff --git a/simulator/rockphysics/softsandrp.py b/simulator/rockphysics/softsandrp.py new file mode 100644 index 0000000..dac1a5a --- /dev/null +++ b/simulator/rockphysics/softsandrp.py @@ -0,0 +1,881 @@ +"""Descriptive description.""" + +__author__ = {'TM', 'TB', 'ML'} + +# standardrp.py +import numpy as np +import sys +import multiprocessing as mp +from CoolProp.CoolProp import PropsSI # http://coolprop.org/#high-level-interface-example +import CoolProp.CoolProp as CP +# Density of carbon dioxide at 100 bar and 25C # Smeaheia 37 degrees C +#rho_co2 = PropsSI('D', 'T', 298.15, 'P', 100e5, 'CO2') +from numpy.random import poisson + +# internal load +from misc.system_tools.environ_var import OpenBlasSingleThread # Single threaded OpenBLAS runs + + +class elasticproperties: + """ + Calculate elastic properties from standard + rock-physics models, specifically following Batzle + and Wang, Geophysics, 1992, for fluid properties, and + Report 1 in Abul Fahimuddin's thesis at Universty of + Bergen (2010) for other properties. + + Example + ------- + >>> porosity = 0.2 + ... pressure = 5 + ... phases = ["Oil","Water"] + ... saturations = [0.3, 0.5] + ... + ... satrock = Elasticproperties() + ... satrock.calc_props(phases, saturations, pressure, porosity) + """ + + def __init__(self, input_dict): + self.dens = None + self.bulkmod = None + self.shearmod = None + self.bulkvel = None + self.shearvel = None + self.bulkimp = None + self.shearimp = None + # The overburden for each grid cell must be + # specified as values on an .npz-file whose + # name is given in input_dict. + self.input_dict = input_dict + self._extInfoInputDict() + + def _extInfoInputDict(self): + # The key word for the file name in the + # dictionary must read "overburden" + if 'overburden' in self.input_dict: + obfile = self.input_dict['overburden'] + npzfile = np.load(obfile) + # The values of overburden must have been + # stored on file using: + # np.savez(, + # obvalues=) + self.overburden = npzfile['obvalues'] + npzfile.close() + #else: + # # Norne litho pressure equation in Bar + # P_litho = -49.6 + 0.2027 * Z + 6.127e-6 * Z ** 2 # Using e-6 for scientific notation + # # Convert reservoir pore pressure from Bar to MPa + # P_litho *= 0.1 + # self.overburden = P_litho + + if 'baseline' in self.input_dict: + self.baseline = self.input_dict['baseline'] # 4D baseline + if 'parallel' in self.input_dict: + self.parallel = self.input_dict['parallel'] + + def _filter(self): + bulkmod = self.bulkimp + self.bulkimp = bulkmod.flatten() + + def setup_fwd_run(self, state): + """ + Setup the input parameters to be used in the PEM simulator. Parameters can be an ensemble or a single array. + State is set as an attribute of the simulator, and the correct value is determined in self.pem.calc_props() + + Parameters + ---------- + state : dict + Dictionary of input parameters or states. + + Changelog + --------- + - KF 11/12-2018 + """ + # self.inv_state = {} + # list_pem_param =[el for el in [foo for foo in self.pem['garn'].keys()] + [foo for foo in self.filter.keys()] + + # [foo for foo in self.__dict__.keys()]] + + # list_tot_param = state.keys() + # for param in list_tot_param: + # if param in list_pem_param or (param.split('_')[-1] in ['garn', 'rest']): + # self.inv_state[param] = state[param] + + pass + + def calc_props(self, phases, saturations, pressure, + porosity, dens = None, wait_for_proc=None, ntg=None, Rs=None, press_init=None, ensembleMember=None): + ### + + # + if not isinstance(phases, list): + phases = [phases] + if not isinstance(saturations, list): + saturations = [saturations] + if not isinstance(pressure, list) and \ + type(pressure).__module__ != 'numpy': + pressure = [pressure] + if not isinstance(porosity, list) and \ + type(porosity).__module__ != 'numpy': + porosity = [porosity] + # + # Load "overburden" pressures into local variable to + # comply with remaining code parts + poverburden = self.overburden + + # debug + self.pressure = pressure + self.peff = poverburden - pressure + self.porosity = porosity + + if press_init is None: + p_init = self.p_init + else: + p_init = press_init + + # Average number of contacts that each grain has with surrounding grains + coordnumber = self._coordination_number() + + # porosity value separating the porous media's mechanical and acoustic behaviour + phicritical = self._critical_porosity() + + + # Check that no. of phases is equal to no. of + # entries in saturations list + # + assert (len(saturations) == len(phases)) + # + # Make saturation a Numpy array (so that we + # can easily access the values for each + # phase at one grid cell) + # + # Transpose makes it a no. grid cells x phases + # array + saturations = np.array(saturations).T + # + # Check if we actually inputted saturation values + # for a single grid cell. If yes, we redefine + # saturations to get it on the correct form (no. + # grid cells x phases array). + # + if saturations.ndim == 1: + saturations = \ + np.array([[x] for x in saturations]).T + # + # Loop over all grid cells and calculate the + # various saturated properties + # + self.phases = phases + + self.dens = np.zeros(len(saturations[:, 0])) + self.bulkmod = np.zeros(len(saturations[:, 0])) + self.shearmod = np.zeros(len(saturations[:, 0])) + self.bulkvel = np.zeros(len(saturations[:, 0])) + self.shearvel = np.zeros(len(saturations[:, 0])) + self.bulkimp = np.zeros(len(saturations[:, 0])) + self.shearimp = np.zeros(len(saturations[:, 0])) + + if ntg is None: + ntg = [None for _ in range(len(saturations[:, 0]))] + if Rs is None: + Rs = [None for _ in range(len(saturations[:, 0]))] + if p_init is None: + p_init = [None for _ in range(len(saturations[:, 0]))] + + + if dens is not None: + assert (len(dens) == len(phases)) + # Transpose makes it a no. grid cells x phases array + dens = np.array(dens).T + + # + denss, bulks, shears = self._solidprops_Johansen() + + for i in range(len(saturations[:, 0])): + # + # Calculate fluid properties + # + if dens is None: + densf_SI = self._fluid_densSIprop(self.phases, + saturations[i, :], pressure[i]) + bulkf_Brie = self._fluidprops_Brie(self.phases, saturations[i, :], pressure[i], densf_SI) + densf, bulkf = \ + self._fluidprops_Wood(self.phases, + saturations[i, :], pressure[i], Rs[i]) + else: + densf = self._fluid_dens(saturations[i, :], dens[i, :]) + + bulkf = self._fluidprops_Brie(self.phases, saturations[i, :], pressure[i], densf) + # + #denss, bulks, shears = \ + # self._solidprops(porosity[i], ntg[i], i) + + # + # Calculate dry rock moduli + # + + #bulkd, sheard = \ + # self._dryrockmoduli(porosity[i], + # overburden[i], + # pressure[i], bulks, + # shears, i, ntg[i], p_init[i], denss, Rs[i], self.phases) + # + peff = self._effective_pressure(poverburden[i], pressure[i]) + + + bulkd, sheard = \ + self._dryrockmoduli_Smeaheia(coordnumber, phicritical, porosity[i], peff, bulks, shears) + + # ------------------------------- + # Calculate saturated properties + # ------------------------------- + # + # Density (kg/m3) + # + self.dens[i] = (porosity[i]*densf + + (1-porosity[i])*denss) + # + # Moduli (MPa) + # + self.bulkmod[i] = \ + bulkd + (1 - bulkd/bulks)**2 / \ + (porosity[i]/bulkf + + (1-porosity[i])/bulks - + bulkd/(bulks**2)) + self.shearmod[i] = sheard + # + # Velocities (km/s) + # + self.bulkvel[i] = \ + np.sqrt((abs(self.bulkmod[i]) + + 4*self.shearmod[i]/3)/(self.dens[i])) + self.shearvel[i] = \ + np.sqrt(self.shearmod[i] / + (self.dens[i])) + # + # convert from (km/s) to (m/s) + # + self.bulkvel[i] *= 1000 + self.shearvel[i] *= 1000 + # + # Impedance (m/s)*(kg/m3) + # + self.bulkimp[i] = self.dens[i] * \ + self.bulkvel[i] + self.shearimp[i] = self.dens[i] * \ + self.shearvel[i] + + + + def getMatchProp(self, petElProp): + if petElProp.lower() == 'density': + self.match_prop = self.getDens() + elif petElProp.lower() == 'bulk_modulus': + self.match_prop = self.getBulkMod() + elif petElProp.lower() == 'shear_modulus': + self.match_prop = self.getShearMod() + elif petElProp.lower() == 'bulk_velocity': + self.match_prop = self.getBulkVel() + elif petElProp.lower() == 'shear_velocity': + self.match_prop = self.getShearVel() + elif petElProp.lower() == "bulk_impedance": + self.match_prop = self.getBulkImp() + elif petElProp.lower() == 'shear_impedance': + self.match_prop = self.getShearImp() + else: + print("\nError in getMatchProp method") + print("No model output type selected for " + "data match.") + print("Legal model output types are " + "(case insensitive):") + print("Density, bulk modulus, shear " + "modulus, bulk velocity,") + print("shear velocity, bulk impedance, " + "shear impedance") + sys.exit(1) + return self.match_prop + + def getDens(self): + return self.dens + + def getBulkMod(self): + return self.bulkmod + + def getShearMod(self): + return self.shearmod + + def getBulkVel(self): + return self.bulkvel + + def getShearVel(self): + return self.shearvel + + def getBulkImp(self): + return self.bulkimp + + def getShearImp(self): + return self.shearimp + + def getOverburdenP(self): + return self.overburden + + def getPressure(self): + return self.pressure + + def getPeff(self): + return self.peff + + def getPorosity(self): + return self.porosity + # + # =================================================== + # Fluid properties start + # =================================================== + # + def _fluid_densSIprop(self, phases, fsats, press, t= 37, CO2 = None): + + conv2Pa = 1e6 # MPa to Pa + ta = t + 273.15 # absolute temp in K + # fluid densities + fdens = 0.0 + + for i in range(len(phases)): + # + # Calculate mixture properties by summing + # over individual phase properties + # + + var = phases[i] + if var == 'GAS' and CO2 is None: + pdens = PropsSI('D', 'T', ta, 'P', press * conv2Pa, 'Methane') + elif var == 'GAS' and CO2 is True: + pdens = PropsSI('D', 'T', ta, 'P', press * conv2Pa, 'CO2') + elif var == 'OIL': + CP.get_global_param_string('predefined_mixtures').split(',')[0:6] + #pdens = CP.PropsSI('D', 'T', ta, 'P', press * conv2Pa, 'Ekofisk.mix') + pdens = CP.PropsSI('D', 'T', ta, 'P', press * conv2Pa, 'butane') + elif var == 'WAT': + pdens = PropsSI('D', 'T|liquid', ta, 'P', press * conv2Pa, 'Water') + + fdens = fdens + fsats[i] * abs(pdens) + + return fdens + + def _fluidprops_Wood(self, fphases, fsats, fpress, Rs=None): + # + # Calculate fluid density and bulk modulus + # + # + # Input + # fphases - fluid phases present; Oil + # and/or Water and/or Gas + # fsats - fluid saturation values for + # fluid phases in "fphases" + # fpress - fluid pressure value (MPa) + # Rs - Gas oil ratio. Default value None + + # + # Output + # fdens - density of fluid mixture for + # pressure value "fpress" inherited + # from phaseprops) + # fbulk - bulk modulus of fluid mixture for + # pressure value "fpress" (unit + # inherited from phaseprops) + # + # ----------------------------------------------- + # + fdens = 0.0 + fbinv = 0.0 + + for i in range(len(fphases)): + # + # Calculate mixture properties by summing + # over individual phase properties + # + pdens, pbulk = self._phaseprops(fphases[i], + fpress, Rs) + fdens = fdens + fsats[i]*abs(pdens) + fbinv = fbinv + fsats[i]/abs(pbulk) + fbulk = 1.0/fbinv + # + return fdens, fbulk + # + # --------------------------------------------------- + # + + def _fluid_dens(self, fsatsp, fdensp): + fdens = sum(fsatsp * fdensp) + return fdens + + def _fluidprops_Brie(self, fphases, fsats, fpress, fdens, Rs=None, e = 5): + # + # Calculate fluid density and bulk modulus BRIE et al. 1995 + # Assumes two phases liquid and gas + # + # Input + # fphases - fluid phases present; Oil + # and/or Water and/or Gas + # fsats - fluid saturation values for + # fluid phases in "fphases" + # fdens - fluid density for given pressure and temperature + # fpress - fluid pressure value (MPa) + # Rs - Gas oil ratio. Default value None + # e - Brie's exponent (e= 5 Utsira sand filled with brine and CO2 + # Figure 7 in Carcione et al. 2006 "Physics and Seismic Modeling + # for Monitoring CO 2 Storage" + + # + # Output + # fbulk - bulk modulus of fluid mixture for + # pressure value "fpress" (unit + # inherited from phaseprops) + # + # ----------------------------------------------- + # + + + for i in range(len(fphases)): + # + if fphases[i].lower() in ["oil", "wat"]: + fsatsl = fsats[i] + pbulkl = self._phaseprops_Smeaheia(fphases[i], fpress, fdens, Rs) + elif fphases[i].lower() in ["gas"]: + pbulkg = self._phaseprops_Smeaheia(fphases[i], fpress, fdens, Rs) + + + fbulk = (pbulkl - pbulkg) * (fsatsl)**e + pbulkg + + # + return fbulk + # + # --------------------------------------------------- + # + @staticmethod + def pseudo_p_t(pres, t, gs): + """Calculate the pseudoreduced temperature and pressure according to Thomas et al. 1970. + + Parameters + ---------- + pres : float or array-like + Pressure in MPa + t : float or array-like + Temperature in °C + gs : float + Gas gravity + + Returns + ------- + float or array-like + Ta: absolute temperature + Ppr:pseudoreduced pressure + Tpr:pseudoreduced temperature + """ + + # convert the temperature to absolute temperature + ta = t + 273.15 + p_pr = pres / (4.892 - 0.4048 * gs) + t_pr = ta / (94.72 + 170.75 * gs) + return ta, p_pr, t_pr + # + # --------------------------------------------------- + # + @staticmethod + def dz_dp(p_pr, t_pr): + """Values for dZ/dPpr obtained from equation 10b in Batzle and Wang (1992). + """ + # analytic + dz_dp = (0.03 + 0.00527 * (3.5 - t_pr) ** 3) + 0.109 * (3.85 - t_pr) ** 2 * 1.2 * p_pr ** 0.2 * -( + 0.45 + 8 * (0.56 - 1 / t_pr) ** 2) / t_pr * np.exp( + -(0.45 + 8 * (0.56 - 1 / t_pr) ** 2) * p_pr ** 1.2 / t_pr) + + # numerical approximation + # dzdp= 1.938783*P_pr**0.2*(1 - 0.25974025974026*T_pr)**2*(-8*(0.56 - 1/T_pr)**2 - 0.45)* + # np.exp(P_pr**1.2*(-8*(0.56 - 1/T_pr)**2 - 0.45)/T_pr)/T_pr + 0.22595125*(1 - 0.285714285714286*T_pr)**3 + # + 0.03 + return dz_dp + # + #----------------------------------------------------------- + # + def _phaseprops_Smeaheia(self, fphase, press, fdens, Rs=None, t = 37, CO2 = True): + # + # Calculate properties for a single fluid phase + # + # + # Input + # fphase - fluid phase; Oil, Water or Gas + # press - fluid pressure value (MPa) + # fdens - fluid density (kg/m3) + # t - temperature in degrees C + # + # Output + # pbulk - bulk modulus of fluid phase + # "fphase" for pressure value + # "press" (MPa) + # + # ----------------------------------------------- + # References + # ---------- + # Xu, H. (2006). Calculation of CO2 acoustic properties using Batzle-Wang equations. Geophysics, 71(2), F21-F23. + # """ + + if fphase.lower() == "wat": # refers to pure water or brine + #Compute the bulk modulus of pure water as a function of temperature and pressure + #using Batzle and Wang (1992). + if np.any(press > 100): + print('pressures above about 100 MPa-> inaccurate estimations of water velocity') + w = np.array([[1.40285e+03, 1.52400e+00, 3.43700e-03, -1.19700e-05], + [4.87100e+00, -1.11000e-02, 1.73900e-04, -1.62800e-06], + [-4.78300e-02, 2.74700e-04, -2.13500e-06, 1.23700e-08], + [1.48700e-04, -6.50300e-07, -1.45500e-08, 1.32700e-10], + [-2.19700e-07, 7.98700e-10, 5.23000e-11, -4.61400e-13]]) + v_w = sum(w[i, j] * t ** i * press ** j for i in range(5) for j in range(4)) # m/s + K_w = fdens * v_w ** 2 * 1e-6 + if CO2 is True: # refers to brine + salinity = 35000 / 1000000 + s1 = 1170 - 9.6 * t + 0.055 * t ** 2 - 8.5e-5 * t ** 3 + 2.6 * press - 0.0029 * t * press - 0.0476 * press ** 2 + s15 = 780 - 10 * press + 0.16 * press ** 2 + s2 = -820 + v_b = v_w + s1 * salinity + s15 * salinity ** 1.5 + s2 * salinity ** 2 + x = 300 * press - 2400 * press * salinity + t * (80 + 3 * t - 3300 * salinity - 13 * press + 47 * press * salinity) + rho_b = fdens + salinity * (0.668 + 0.44 * salinity + 1e-6 * x) + pbulk = rho_b * v_b ** 2 * 1e-6 + else: + pbulk = K_w + + elif fphase.lower() == "gas" and CO2 is True: # refers to CO2 + R = 8.3145 # J.mol-1K-1 gas constant for CO2 + gs = 1.5189 # Specific gravity #https://www.engineeringtoolbox.com/specific-gravities-gases-d_334.html + ta, p_pr, t_pr = self.pseudo_p_t(press, t, gs) + + E = 0.109 * (3.85 - t_pr) ** 2 * np.exp(-(0.45 + 8 * (0.56 - 1 / t_pr) ** 2) * p_pr ** 1.2 / t_pr) + Z = (0.03 + 0.00527 * (3.5 - t_pr) ** 3) * p_pr + (0.642 * t_pr - 0.007 * t_pr ** 4 - 0.52) + E + rho = 28.8 * gs * press / (Z * R * ta) # g/cm3 + + r_0 = 0.85 + 5.6 / (p_pr + 2) + 27.1 / (p_pr + 3.5) ** 2 - 8.7 * np.exp(-0.65 * (p_pr + 1)) + dz_dp = self.dz_dp(p_pr, t_pr) + pbulk = press / (1 - p_pr * dz_dp / Z) * r_0 + + #pbulk_test = self.test_new_implementation(press) + #print(np.max(pbulk-pbulk_test)) + + elif fphase.lower() == "gas": # refers to Methane + gs = 0.5537 #https://www.engineeringtoolbox.com/specific-gravities-gases-d_334.html + R = 8.3145 # J.mol-1K-1 gas constant + ta, p_pr, t_pr = self.pseudo_p_t(press, t, gs) + E = 0.109 * (3.85 - t_pr) ** 2 * np.exp(-(0.45 + 8 * (0.56 - 1 / t_pr) ** 2) * p_pr ** 1.2 / t_pr) + Z = (0.03 + 0.00527 * (3.5 - t_pr) ** 3) * p_pr + (0.642 * t_pr - 0.007 * t_pr ** 4 - 0.52) + E + rho = 28.8 * gs * press / (Z * R * ta) # g/cm3 + + r_0 = 0.85 + 5.6 / (p_pr + 2) + 27.1 / (p_pr + 3.5) ** 2 - 8.7 * np.exp(-0.65 * (p_pr + 1)) + dz_dp = self.dz_dp(p_pr, t_pr) + pbulk = press / (1 - p_pr * dz_dp / Z) * r_0 + + elif fphase.lower() == "oil": #pure oil + # Estimate the oil bulk modulus at specific temperature and pressure. + v = 2096 * (fdens / (2600 - fdens)) ** 0.5 - 3.7 * t + 4.64 * press + 0.0115 * ( + 4.12 * (1080 / fdens - 1) ** 0.5 - 1) * t * press + pbulk = fdens * v ** 2 + + + # + return pbulk + + # + def test_new_implementation(self, press): + # Values from .DATA file for Smeaheia (converted to MPa) + press_range = np.array( + [0.101, 0.885, 1.669, 2.453, 3.238, 4.022, 4.806, 5.590, 6.2098, 7.0899, 7.6765, 8.2630, 8.8495, 9.4359, + 10.0222, 10.6084, 11.1945, 14.7087, 17.6334, 20.856, 23.4695, 27.5419]) # Example pressures in MPa + Bo_values = np.array( + [1.07365, 0.11758, 0.05962, 0.03863, 0.02773, 0.02100, 0.01639, 0.01298, 0.010286, 0.007578, 0.005521, + 0.003314, 0.003034, 0.002919, 0.002851, 0.002802, 0.002766, 0.002648, 0.002599, 0.002566, 0.002546, + 0.002525]) # Example formation volume factors in m^3/kg + + # Calculate numerical derivative of Bo with respect to Pressure + dBo_dP = - np.gradient(Bo_values, press_range) + # Calculate isothermal compressibility (van der Waals) + compressibility = (1 / Bo_values) * dBo_dP # Resulting array of compressibility values + bulk_mod = 1 / compressibility + + # Find the index of the closest pressure value in b + closest_index = (np.abs(press_range - press)).argmin() + + # Extract the corresponding value from a + pbulk_test = bulk_mod[closest_index] + return pbulk_test + + def _phaseprops(self, fphase, press, Rs=None): + # + # Calculate properties for a single fluid phase + # + # + # Input + # fphase - fluid phase; Oil, Water or Gas + # press - fluid pressure value (MPa) + # + # Output + # pdens - phase density of fluid phase + # "fphase" for pressure value + # "press" (kg/m³) + # pbulk - bulk modulus of fluid phase + # "fphase" for pressure value + # "press" (MPa) + # + # ----------------------------------------------- + # + if fphase.lower() == "oil": + coeffsrho = np.array([0.8, 829.9]) + coeffsbulk = np.array([10.42, 995.79]) + elif fphase.lower() == "wat": + coeffsrho = np.array([0.3, 1067.3]) + coeffsbulk = np.array([9.0, 2807.6]) + elif fphase.lower() == "gas": + coeffsrho = np.array([4.7, 13.4]) + coeffsbulk = np.array([2.75, 0.0]) + else: + print("\nError in phaseprops method") + print("Illegal fluid phase name.") + print("Legal fluid phase names are (case " + "insensitive): Oil, Wat, and Gas.") + sys.exit(1) + # + # Assume simple linear pressure dependencies. + # Coefficients are inferred from + # plots in Batzle and Wang, Geophysics, 1992, + # (where I set the temperature to be 100 degrees + # Celsius, Note also that they give densities in + # g/cc). The resulting straight lines do not fit + # the data extremely well, but they should + # be sufficiently accurate for the purpose of + # this project. + # + pdens = coeffsrho[0]*press + coeffsrho[1] + pbulk = coeffsbulk[0]*press + coeffsbulk[1] + # + return pdens, pbulk + + # + # ======================= + # Fluid properties end + # ======================= + # + + # + # ========================= + # Solid properties start + # ========================= + # + def _solidprops_Johansen(self): + # + # Calculate bulk and shear solid rock (mineral) + # moduli by averaging Hashin-Shtrikman bounds + # + # + # Input + # poro -porosity + # + # Output + # denss - solid rock density (kg/m³) + # bulks - solid rock bulk modulus (unit MPa) + # shears - solid rock shear modulus (unit MPa) + # + # ----------------------------------------------- + # + # + # Solid rock (mineral) density. (Note + # that this is often termed \rho_dry, and not + # \rho_s) + + denss = 2650 # Density of mineral/solid rock kg/m3 + + # + bulks = 37 # (GPa) + shears = 44 # (GPa) + bulks *= 1000 # Convert from GPa to MPa + shears *= 1000 + # + return denss, bulks, shears + # + # + # ======================= + # Solid properties end + # ======================= + # + def _coordination_number(self): + # Applies for granular media + # Average number of contacts that each grain has with surrounding grains + # Coordnumber = 6; simple cubic packing + # Coordnumber = 12; hexagonal close packing + # Needed for Hertz-Mindlin model + # Smeaheia number (Tuhin) + coordnumber = 9 + + return coordnumber + # + def _critical_porosity(self): + # For most porous media there exists a critical porosity + # phi_critical, that seperates their mechanical and acoustic behaviour into two domains. + # For porosities below phi_critical the mineral grains are oad bearing, for values above the grains are + # suspended in the fluids which are load-bearing + # Needed for Hertz-Mindlin model + # Smeaheia number (Tuhin) + phicritical = 0.36 + + return phicritical + # + def _effective_pressure(self, poverb, pfluid): + + # Input + # poverb - overburden pressure (MPa) + # pfluid - fluid pressure (MPa) + + peff = poverb - pfluid + + if peff < 0: + print("\nError in _hertzmindlin method") + print("Negative effective pressure (" + str(peff) + + "). Setting effective pressure to 0.01") + peff = 0.01 + + + + return peff + + # ============================ + # Dry rock properties start + # ============================ + # + def _dryrockmoduli_Smeaheia(self, coordnumber, phicritical, poro, peff, bulks, shears): + # + # + # Calculate bulk and shear dry rock moduli, + + # + # Input + # poro - porosity + # peff - effective pressure overburden - fluid pressure (MPa) + # bulks - bulk solid (mineral) rock bulk + # modulus (MPa) + # shears - solid rock (mineral) shear + # modulus (MPa) + # + # Output + # bulkd - dry rock bulk modulus (unit + # inherited from hertzmindlin and + # variable; bulks) + # sheard - dry rock shear modulus (unit + # inherited from hertzmindlin and + # variable; shears) + # + # ----------------------------------------------- + # + # Calculate Hertz-Mindlin moduli + # + bulkhm, shearhm = self._hertzmindlin_Mavko(peff, bulks, shears, coordnumber, phicritical) + # + bulkd = 1 / ((poro / phicritical) / (bulkhm + 4 / 3 * shearhm) + + (1 - poro / phicritical) / (bulks + 4 / 3 * shearhm)) - 4 / 3 * shearhm + + psi = (9 * bulkhm + 8 * shearhm) / (bulkhm + 2 * shearhm) + + sheard = 1 / ((poro / phicritical) / (shearhm + 1 / 6 * psi * shearhm) + + (1 - poro / phicritical) / (shears + 1 / 6 * psi * shearhm)) - 1 / 6 * psi * shearhm + + #return K_dry, G_dry + return bulkd, sheard + + + # + # --------------------------------------------------- + # + + def _hertzmindlin_Mavko(self, peff, bulks, shears, coordnumber, phicritical): + # + # Calculate bulk and shear Hertz-Mindlin moduli + # adapted from Tuhins kode and "The rock physics handbook", pp247 + # + # + # Input + # p_eff - effective pressure + # bulks - bulk solid (mineral) rock bulk + # modulus (MPa) + # shears - solid rock (mineral) shear + # modulus (MPa) + # coordnumber - average number of contacts that each grain has with surrounding grains + # phicritical - critical porosity + # + # Output + # bulkhm - Hertz-Mindlin bulk modulus + # (MPa) + # shearhm - Hertz-Mindlin shear modulus + # (MPa) + # + # ----------------------------------------------- + # + + + poisson = (3 * bulks - 2 * shears) / (6 * bulks + 2 * shears) + + bulkhm = ((coordnumber ** 2 * (1 - phicritical) ** 2 * shears ** 2 * peff) / + (18 * np.pi ** 2 * (1 - poisson) ** 2)) ** (1 / 3) + shearhm = (5 - 4 * poisson) / (10 - 5 * poisson) * \ + ((3 * coordnumber ** 2 * (1 - phicritical) ** 2 * shears ** 2 * peff) / + (2 * np.pi ** 2 * (1 - poisson) ** 2)) ** (1 / 3) + + + + # + return bulkhm, shearhm + + + # =========================== + # Dry rock properties end + # =========================== + + +if __name__ == '__main__': + # + # Example input with two phases and three grid cells + # + porosity = [0.34999999, 0.34999999, 0.34999999] +# pressure = [ 29.29150963, 29.14003944, 28.88845444] + pressure = [29.3558, 29.2625, 29.3558] +# pressure = [ 25.0, 25.0, 25.0] + phases = ["Oil", "Wat"] +# saturations = [[0.72783828, 0.66568458, 0.58033288], +# [0.27216172, 0.33431542, 0.41966712]] + saturations = [[0.6358, 0.5755, 0.6358], + [0.3641, 0.4245, 0.3641]] +# saturations = [[0.4, 0.5, 0.6], +# [0.6, 0.5, 0.4]] + petElProp = "bulk velocity" + input_dict = {} + input_dict['overburden'] = 'overb.npz' + + print("\nInput:") + print("porosity, pressure:", porosity, pressure) + print("phases, saturations:", phases, saturations) + print("petElProp:", petElProp) + print("input_dict:", input_dict) + + satrock = elasticproperties(input_dict) + + print("overburden:", satrock.overburden) + + satrock.calc_props(phases, saturations, pressure, + porosity) + + print("\nOutput from calc_props:") + print("Density:", satrock.getDens()) + print("Bulk modulus:", satrock.getBulkMod()) + print("Shear modulus:", satrock.getShearMod()) + print("Bulk velocity:", satrock.getBulkVel()) + print("Shear velocity:", satrock.getShearVel()) + print("Bulk impedance:", satrock.getBulkImp()) + print("Shear impedance:", satrock.getShearImp()) + + satrock.getMatchProp(petElProp) + + print("\nOutput from getMatchProp:") + print("Model output selected for data match:", + satrock.match_prop) diff --git a/simulator/rockphysics/standardrp.py b/simulator/rockphysics/standardrp.py index f847ed0..dde7b09 100644 --- a/simulator/rockphysics/standardrp.py +++ b/simulator/rockphysics/standardrp.py @@ -90,7 +90,7 @@ def setup_fwd_run(self, state): pass def calc_props(self, phases, saturations, pressure, - porosity, wait_for_proc=None, ntg=None, Rs=None, press_init=None, ensembleMember=None): + porosity, dens = None, wait_for_proc=None, ntg=None, Rs=None, press_init=None, ensembleMember=None): ### # This doesn't initialize for models with no uncertainty ### @@ -135,11 +135,19 @@ def calc_props(self, phases, saturations, pressure, # Load "overburden" into local variable to # comply with remaining code parts overburden = self.overburden + if press_init is None: p_init = self.p_init else: p_init = press_init # + poverburden = self.overburden + + # debug + self.pressure = pressure + self.peff = poverburden - pressure + self.porosity = porosity + # Check that no. of phases is equal to no. of # entries in saturations list # @@ -208,8 +216,7 @@ def calc_props(self, phases, saturations, pressure, # # Density # - self.dens[i] = (porosity[i]*densf + - (1-porosity[i])*denss)*0.001 + self.dens[i] = (porosity[i]*densf + (1-porosity[i])*denss) # # Moduli # @@ -225,10 +232,10 @@ def calc_props(self, phases, saturations, pressure, # instead of km/s) # self.bulkvel[i] = \ - 100*np.sqrt((abs(self.bulkmod[i]) + + 1000*np.sqrt((abs(self.bulkmod[i]) + 4*self.shearmod[i]/3)/(self.dens[i])) self.shearvel[i] = \ - 100*np.sqrt(self.shearmod[i] / + 1000*np.sqrt(self.shearmod[i] / (self.dens[i])) # # Impedances (m/s)*(Kg/m3) @@ -287,6 +294,17 @@ def getBulkImp(self): def getShearImp(self): return self.shearimp + def getOverburdenP(self): + return self.overburden + + def getPressure(self): + return self.pressure + + def getPeff(self): + return self.peff + + def getPorosity(self): + return self.porosity # # =================================================== # Fluid properties start