diff --git a/ensemble/ensemble.py b/ensemble/ensemble.py index fb984d5..676a059 100644 --- a/ensemble/ensemble.py +++ b/ensemble/ensemble.py @@ -553,6 +553,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 self.keys_da['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))] diff --git a/pipt/loop/assimilation.py b/pipt/loop/assimilation.py index 0ccd6dd..926b1a8 100644 --- a/pipt/loop/assimilation.py +++ b/pipt/loop/assimilation.py @@ -580,8 +580,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 675fcc3..131b3c1 100644 --- a/pipt/loop/ensemble.py +++ b/pipt/loop/ensemble.py @@ -732,7 +732,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( @@ -741,7 +741,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 86e7b4c..04e51c5 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/update_schemes/gies/gies_base.py b/pipt/update_schemes/gies/gies_base.py new file mode 100644 index 0000000..f967f63 --- /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_state() + 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 605c44b..e1bfbf0 100644 --- a/setup.py +++ b/setup.py @@ -43,6 +43,7 @@ 'tomli-w', 'pyyaml', 'libecalc', - '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..31aa975 100644 --- a/simulator/flow_rock.py +++ b/simulator/flow_rock.py @@ -12,7 +12,17 @@ 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 class flow_sim2seis(flow): """ @@ -55,7 +65,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] @@ -106,64 +116,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.sats = [] + 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 inhere flor_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() @@ -656,3 +622,477 @@ def extract_data(self, member): 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() + + +class flow_avo(flow_sim2seis): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): + super().__init__(input_dict, filename, options) + + 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() + + 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. + """ + + if member_i >= 0: + folder = 'En_' + str(member_i) + os.sep + if not os.path.exists(folder): + os.mkdir(folder) + else: # XLUO: any negative member_i is considered as the index for the true model + assert 'truth_folder' in self.input_dict, "ensemble member index is negative, please specify " \ + "the folder containing the true model" + if not os.path.exists(self.input_dict['truth_folder']): + os.mkdir(self.input_dict['truth_folder']) + folder = self.input_dict['truth_folder'] + os.sep if self.input_dict['truth_folder'][-1] != os.sep \ + else self.input_dict['truth_folder'] + del_folder = False # never delete this folder + self.folder = folder + self.ensemble_member = member_i + + state['member'] = member_i + + # start by generating the .DATA file, using the .mako template situated in ../folder + self._runMako(folder, state) + success = False + rerun = self.rerun + while rerun >= 0 and not success: + success = self.call_sim(folder, True) + rerun -= 1 + if success: + self.extract_data(member_i) + if del_folder: + if self.saveinfo is not None: # Try to save information + store_ensemble_sim_information(self.saveinfo, member_i) + self.remove_folder(member_i) + return self.pred_data + else: + if hasattr(self, 'redund_sim') and self.redund_sim is not None: + success = self.redund_sim.call_sim(folder, True) + if success: + self.extract_data(member_i) + if del_folder: + if self.saveinfo is not None: # Try to save information + store_ensemble_sim_information(self.saveinfo, member_i) + self.remove_folder(member_i) + return self.pred_data + else: + if del_folder: + self.remove_folder(member_i) + return False + else: + if del_folder: + self.remove_folder(member_i) + return False + + + 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 + + # The field 'run_reservoir_model' can be passed from the method "setup_fwd_run" + if hasattr(self, 'run_reservoir_model'): + run_reservoir_model = self.run_reservoir_model + + if run_reservoir_model is None: + run_reservoir_model = True + + # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. + # Then, get the pem. + if run_reservoir_model: # in case that simulation has already done (e.g., for the true reservoir model) + success = super(flow_sim2seis, self).call_sim(folder, wait_for_proc) + #ecl = ecl_100(filename=self.file) + #ecl.options = self.options + #success = ecl.call_sim(folder, wait_for_proc) + else: + success = True + + if success: + 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() + + #phases = self.ecl_case.init.phases + self.sats = [] + 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) + # vp, vs, density in reservoir + # The properties in pem is only given in the active cells + # indices of active cells: + if grid['ACTNUM'].shape[0] == self.NX: + true_indices = np.where(grid['ACTNUM']) + elif grid['ACTNUM'].shape[0] == self.NZ: + actnum = np.transpose(grid['ACTNUM'], (2, 1, 0)) + true_indices = np.where(actnum) + else: + print('warning: dimension mismatch in line 750 flow_rock.py') + + if len(self.pem.getBulkVel()) == len(true_indices[0]): + self.vp = np.zeros(grid['DIMENS']) + self.vp[true_indices] = (self.pem.getBulkVel() * .1) + self.vs = np.zeros(grid['DIMENS']) + self.vs[true_indices] = (self.pem.getShearVel() * .1) + self.rho = np.zeros(grid['DIMENS']) + self.rho[true_indices] = (self.pem.getDens()) + else: + self.vp = (self.pem.getBulkVel() * .1).reshape((self.NX, self.NY, self.NZ), order='F') + self.vs = (self.pem.getShearVel() * .1).reshape((self.NX, self.NY, self.NZ), order='F') + self.rho = (self.pem.getDens()).reshape((self.NX, self.NY, self.NZ), order='F') # in the unit of g/cm^3 + + save_dic = {'vp': self.vp, 'vs': self.vs, 'rho': self.vp} + 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" + else: + if hasattr(self, 'ensemble_member') and (self.ensemble_member is not None) and \ + (self.ensemble_member >= 0): + file_name = folder + os.sep + f"vp_vs_rho_vint{v}.npz" if folder[-1] != os.sep \ + else folder + f"vp_vs_rho_vint{v}.npz" + else: + file_name = os.getcwd() + os.sep + f"vp_vs_rho_vint{v}.npz" + + #with open(file_name, "wb") as f: + # dump(**save_dic, f) + np.savez(file_name, **save_dic) + + # avo data + self._calc_avo_props() + + avo = self.avo_data.flatten(order="F") + + # 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 + + save_dic = {'avo': avo, 'noise_std': noise_std, **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" + else: + file_name = folder + os.sep + f"avo_vint{v}.npz" if folder[-1] != os.sep \ + else folder + f"avo_vint{v}.npz" + + #with open(file_name, "wb") as f: + # dump(**save_dic, f) + np.savez(file_name, **save_dic) + + return success + + def extract_data(self, member): + # start by getting the data from the flow simulator + super(flow_sim2seis, self).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 '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] + + def _runMako(self, folder, state, addfiles=['properties']): + """ + Hard coding, maybe a better way possible + addfiles: additional files that need to be included into ECLIPSE/OPM DATA file + """ + super()._runMako(folder, state) + + lkup = TemplateLookup(directories=os.getcwd(), input_encoding='utf-8') + for file in addfiles: + if os.path.exists(file + '.mako'): + tmpl = lkup.get_template('%s.mako' % file) + + # use a context and render onto a file + with open('{0}'.format(folder + file), 'w') as f: + ctx = Context(f, **state) + tmpl.render_context(ctx) + + 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) + + + + if 'RS' in self.pem_input: #ecl_case.cell_data: # to be more robust! + 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'] + + + 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 + for var in phases: + if var in ['WAT', 'GAS']: + tmp = self.ecl_case.cell_data('S{}'.format(var), time) + pem_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) + + 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 + for var in phases: + if var in ['GAS']: + tmp = self.ecl_case.cell_data('S{}'.format(var), time) + pem_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) + 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_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) + + 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' + 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: + 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 trough the reservoir in vertical direction + cum_time_res = np.cumsum(2 * self.DZ / self.vp, axis=2) + top_res[:, :, np.newaxis] + # total travel time + cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res, top_res[:, :, 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) + + # 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] + + + # from matplotlib import pyplot as plt + # plt.plot(vp_sample[0, 0, :]) + # plt.show() + + #vp_avg = 0.5 * (vp_sample[:, :, 1:] + vp_sample[:, :, :-1]) + #vs_avg = 0.5 * (vs_sample[:, :, 1:] + vs_sample[:, :, :-1]) + #rho_avg = 0.5 * (rho_sample[:, :, 1:] + rho_sample[:, :, :-1]) + + #vp_diff = vp_sample[:, :, 1:] - vp_sample[:, :, :-1] + #vs_diff = vs_sample[:, :, 1:] - vs_sample[:, :, :-1] + #rho_diff = rho_sample[:, :, 1:] - rho_sample[:, :, :-1] + + #R0_smith = 0.5 * (vp_diff / vp_avg + rho_diff / rho_avg) + #G_smith = -2.0 * (vs_avg / vp_avg) ** 2 * (2.0 * vs_diff / vs_avg + rho_diff / rho_avg) + 0.5 * vp_diff / vp_avg + + # PP reflection coefficients, see, e.g., + # "https://pylops.readthedocs.io/en/latest/api/generated/pylops.avo.avo.approx_zoeppritz_pp.html" + # So far, it seems that "approx_zoeppritz_pp" is the only available option + # approx_zoeppritz_pp(vp1, vs1, rho1, vp0, vs0, rho0, theta1) + avo_data_list = [] + + # 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 sereis + 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 vertial 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 + + @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 arraies 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 diff --git a/simulator/flow_rock_backup.py b/simulator/flow_rock_backup.py new file mode 100644 index 0000000..7eb2e03 --- /dev/null +++ b/simulator/flow_rock_backup.py @@ -0,0 +1,1120 @@ +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 + +class flow_sim2seis(flow): + """ + Couple the OPM-flow simulator with a sim2seis 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, 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']) + + # 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 = [] + + 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] + + 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 setup_fwd_run(self): + super().setup_fwd_run() + + 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) + + 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. + success = super().call_sim(folder, wait_for_proc) + + if success: + # need a 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) + + self.ecl_case = ecl.EclipseCase( + 'En_' + str(self.ensemble_member) + os.sep + self.file + '.DATA') + 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']) + \ + 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) + + 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', { + 'Vp': self.pem.getBulkVel()*.1, 'DIMENS': grid['DIMENS']}, multi_file=False) + 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() + run_folder = current_folder + os.sep + 'En_' + str(self.ensemble_member) + # The sim2seis is invoked via a shell script. The simulations provides outputs. Run, and get all output. Search + # for Done. If not finished in reasonable time -> kill + p = Popen(['./sim2seis.sh', run_folder], stdout=PIPE) + start = time + while b'done' not in p.stdout.readline(): + pass + + # Todo: handle sim2seis or pem error + + 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 ['sim2seis']: + if self.true_prim[1][prim_ind] in self.pem_input['vintage']: + result = mat73.loadmat(f'En_{member}/Data_conv.mat')['data_conv'] + v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = np.sum( + np.abs(result[:, :, :, v]), axis=0).flatten() + +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, 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']) + + # 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 = [] + + 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 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] + + 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 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 + self.pred_data = super().run_fwd_sim(state, member_i, del_folder=True) + + 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. + success = super().call_sim(folder, wait_for_proc) + + 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) + + 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) + + 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) + + 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) + # 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 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') + + if 'compaction' in self.pem_input: + multfactor = self.ecl_case.cell_data('PORV_RC', base_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) + + 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) + + if 'press_conv' in self.pem_input: + pem_input['PRESSURE'] = pem_input['PRESSURE'] * \ + 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=None) + + # 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 + # 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_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. + """ + + def __init__(self, input_dict=None, filename=None, options=None): + 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']) + + # 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 = [] + self.pem_result = [] + self.bar_result = [] + + 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 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] + + 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 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 + self.pred_data = super().run_fwd_sim(state, member_i, del_folder=True) + + 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. + success = super().call_sim(folder, wait_for_proc) + + 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) + + 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) + + 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) + + 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) + # 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 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') + + if 'compaction' in self.pem_input: + multfactor = self.ecl_case.cell_data('PORV_RC', base_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) + + 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) + + if 'press_conv' in self.pem_input: + pem_input['PRESSURE'] = pem_input['PRESSURE'] * \ + 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=None) + + # 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 + # 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 + 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 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 ['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.bar_result[v].flatten() + +class flow_avo(flow_sim2seis): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): + super().__init__(input_dict, filename, options) + + 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() + + 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. + """ + + if member_i >= 0: + folder = 'En_' + str(member_i) + os.sep + if not os.path.exists(folder): + os.mkdir(folder) + else: # XLUO: any negative member_i is considered as the index for the true model + assert 'truth_folder' in self.input_dict, "ensemble member index is negative, please specify " \ + "the folder containing the true model" + if not os.path.exists(self.input_dict['truth_folder']): + os.mkdir(self.input_dict['truth_folder']) + folder = self.input_dict['truth_folder'] + os.sep if self.input_dict['truth_folder'][-1] != os.sep \ + else self.input_dict['truth_folder'] + del_folder = False # never delete this folder + self.folder = folder + self.ensemble_member = member_i + + state['member'] = member_i + + # start by generating the .DATA file, using the .mako template situated in ../folder + self._runMako(folder, state) + success = False + rerun = self.rerun + while rerun >= 0 and not success: + success = self.call_sim(folder, True) + rerun -= 1 + if success: + self.extract_data(member_i) + if del_folder: + if self.saveinfo is not None: # Try to save information + store_ensemble_sim_information(self.saveinfo, member_i) + self.remove_folder(member_i) + return self.pred_data + else: + if hasattr(self, 'redund_sim') and self.redund_sim is not None: + success = self.redund_sim.call_sim(folder, True) + if success: + self.extract_data(member_i) + if del_folder: + if self.saveinfo is not None: # Try to save information + store_ensemble_sim_information(self.saveinfo, member_i) + self.remove_folder(member_i) + return self.pred_data + else: + if del_folder: + self.remove_folder(member_i) + return False + else: + if del_folder: + self.remove_folder(member_i) + return False + + + 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 + + # The field 'run_reservoir_model' can be passed from the method "setup_fwd_run" + if hasattr(self, 'run_reservoir_model'): + run_reservoir_model = self.run_reservoir_model + + if run_reservoir_model is None: + run_reservoir_model = True + + # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. + # Then, get the pem. + if run_reservoir_model: # in case that simulation has already done (e.g., for the true reservoir model) + success = super(flow_sim2seis, self).call_sim(folder, wait_for_proc) + #ecl = ecl_100(filename=self.file) + #ecl.options = self.options + #success = ecl.call_sim(folder, wait_for_proc) + else: + success = True + + if success: + 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() + + 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']) + \ + 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) + #pem_input['PORO'] = np.array(self._reformat3D_then_flatten(multfactor * tmp), dtype=float) + else: + pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float) + #pem_input['PORO'] = np.array(self._reformat3D_then_flatten(tmp), 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) + #pem_input['NTG'] = np.array(self._reformat3D_then_flatten(tmp), dtype=float) + else: + tmp = self.ecl_case.cell_data('NTG') + pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) + #pem_input['NTG'] = np.array(self._reformat3D_then_flatten(tmp), 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) + #pem_input[var] = np.array(self._reformat3D_then_flatten(tmp), dtype=float) + + if 'press_conv' in self.pem_input: + pem_input['PRESSURE'] = pem_input['PRESSURE'] * \ + self.pem_input['press_conv'] + #pem_input['PRESSURE'] = self._reformat3D_then_flatten(pem_input['PRESSURE'] * + # self.pem_input['press_conv']) + + tmp = self.ecl_case.cell_data('PRESSURE', 0) + if hasattr(self.pem, 'p_init'): + P_init = self.pem.p_init * np.ones(tmp.shape)[~tmp.mask] + #P_init = self._reformat3D_then_flatten(self.pem.p_init.reshape(tmp.shape) * np.ones(tmp.shape)) + else: + # initial pressure is first + P_init = np.array(tmp[~tmp.mask], dtype=float) + #P_init = np.array(self._reformat3D_then_flatten(tmp), dtype=float) + + if 'press_conv' in self.pem_input: + P_init = P_init * self.pem_input['press_conv'] + #P_init = self._reformat3D_then_flatten(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] + #saturations = [self._reformat3D_then_flatten(1 - (pem_input['SWAT'] + pem_input['SGAS'])) + # if ph == 'OIL' else pem_input['S{}'.format(ph)] for ph in phases] + + # Get the pressure + 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) + + #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', { + # 'Vp': self.pem.getBulkVel() * .1, 'DIMENS': grid['DIMENS']}, multi_file=False) + #grdecl.write(f'En_{str(self.ensemble_member)}/rho{v + 1}.grdecl', + # {'rho': self.pem.getDens(), 'DIMENS': grid['DIMENS']}, multi_file=False) + + # vp, vs, density + self.vp = (self.pem.getBulkVel() * .1).reshape((self.NX, self.NY, self.NZ), order='F') + self.vs = (self.pem.getShearVel() * .1).reshape((self.NX, self.NY, self.NZ), order='F') + self.rho = (self.pem.getDens()).reshape((self.NX, self.NY, self.NZ), order='F') # in the unit of g/cm^3 + + save_dic = {'vp': self.vp, 'vs': self.vs, 'rho': self.vp} + 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" + else: + if hasattr(self, 'ensemble_member') and (self.ensemble_member is not None) and \ + (self.ensemble_member >= 0): + file_name = folder + os.sep + f"vp_vs_rho_vint{v}.npz" if folder[-1] != os.sep \ + else folder + f"vp_vs_rho_vint{v}.npz" + else: + file_name = os.getcwd() + os.sep + f"vp_vs_rho_vint{v}.npz" + + #with open(file_name, "wb") as f: + # dump(**save_dic, f) + np.savez(file_name, **save_dic) + + # avo data + self._calc_avo_props() + + avo = self.avo_data.flatten(order="F") + + # 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 + + save_dic = {'avo': avo, 'noise_std': noise_std, **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" + else: + # if hasattr(self, 'ensemble_member') and (self.ensemble_member is not None) and \ + # (self.ensemble_member >= 0): + # file_name = folder + os.sep + f"avo_vint{v}.npz" if folder[-1] != os.sep \ + # else folder + f"avo_vint{v}.npz" + # else: + # file_name = os.getcwd() + os.sep + f"avo_vint{v}.npz" + file_name = folder + os.sep + f"avo_vint{v}.npz" if folder[-1] != os.sep \ + else folder + f"avo_vint{v}.npz" + + #with open(file_name, "wb") as f: + # dump(**save_dic, f) + np.savez(file_name, **save_dic) + + return success + + def extract_data(self, member): + # start by getting the data from the flow simulator + super(flow_sim2seis, self).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 '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] + + def _runMako(self, folder, state, addfiles=['properties']): + """ + Hard coding, maybe a better way possible + addfiles: additional files that need to be included into ECLIPSE/OPM DATA file + """ + super()._runMako(folder, state) + + lkup = TemplateLookup(directories=os.getcwd(), input_encoding='utf-8') + for file in addfiles: + if os.path.exists(file + '.mako'): + tmpl = lkup.get_template('%s.mako' % file) + + # use a context and render onto a file + with open('{0}'.format(folder + file), 'w') as f: + ctx = Context(f, **state) + tmpl.render_context(ctx) + + def calc_pem(self, time): + + + 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' + 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: + 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 trough the reservoir in vertical direction + cum_time_res = np.cumsum(2 * self.DZ / self.vp, axis=2) + top_res[:, :, np.newaxis] + # total travel time + cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res, top_res[:, :, 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) + + # 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] + + + # from matplotlib import pyplot as plt + # plt.plot(vp_sample[0, 0, :]) + # plt.show() + + #vp_avg = 0.5 * (vp_sample[:, :, 1:] + vp_sample[:, :, :-1]) + #vs_avg = 0.5 * (vs_sample[:, :, 1:] + vs_sample[:, :, :-1]) + #rho_avg = 0.5 * (rho_sample[:, :, 1:] + rho_sample[:, :, :-1]) + + #vp_diff = vp_sample[:, :, 1:] - vp_sample[:, :, :-1] + #vs_diff = vs_sample[:, :, 1:] - vs_sample[:, :, :-1] + #rho_diff = rho_sample[:, :, 1:] - rho_sample[:, :, :-1] + + #R0_smith = 0.5 * (vp_diff / vp_avg + rho_diff / rho_avg) + #G_smith = -2.0 * (vs_avg / vp_avg) ** 2 * (2.0 * vs_diff / vs_avg + rho_diff / rho_avg) + 0.5 * vp_diff / vp_avg + + # PP reflection coefficients, see, e.g., + # "https://pylops.readthedocs.io/en/latest/api/generated/pylops.avo.avo.approx_zoeppritz_pp.html" + # So far, it seems that "approx_zoeppritz_pp" is the only available option + # approx_zoeppritz_pp(vp1, vs1, rho1, vp0, vs0, rho0, theta1) + avo_data_list = [] + + # 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 sereis + 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 vertial 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 + + @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 arraies 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 \ No newline at end of file diff --git a/simulator/flow_rock_mali.py b/simulator/flow_rock_mali.py new file mode 100644 index 0000000..6abb95d --- /dev/null +++ b/simulator/flow_rock_mali.py @@ -0,0 +1,1130 @@ +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 + +class flow_sim2seis(flow): + """ + Couple the OPM-flow simulator with a sim2seis 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, 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']) + + # 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 = [] + self.sats = [] + + 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] + + 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 calc_pem(self, v, time, 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 ['PRESSURE', 'RS', 'SWAT', 'SGAS']: + 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 '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'] + + saturations = [1 - (pem_input['SWAT'] + pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)] + for ph in phases] + + # fluid saturations in dictionary + tmp_s = {f'S{ph}': saturations[i] for i, ph in enumerate(phases)} + self.sats.extend([tmp_s]) + + # 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 setup_fwd_run(self): + super().setup_fwd_run() + + 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) + + 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. + success = super().call_sim(folder, wait_for_proc) + + if success: + # need a 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) + + self.ecl_case = ecl.EclipseCase( + 'En_' + str(self.ensemble_member) + os.sep + self.file + '.DATA') + + 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 = [] + self.sats = [] + # loop over seismic vintages + for v, assim_time in enumerate([0.0] + self.pem_input['vintage']): + time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ + dt.timedelta(days=assim_time) + + self.calc_pem(v, time, phases) + + 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', { + 'Vp': self.pem.getBulkVel() * .1, 'DIMENS': grid['DIMENS']}, multi_file=False) + 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() + run_folder = current_folder + os.sep + 'En_' + str(self.ensemble_member) + # The sim2seis is invoked via a shell script. The simulations provides outputs. Run, and get all output. Search + # for Done. If not finished in reasonable time -> kill + p = Popen(['./sim2seis.sh', run_folder], stdout=PIPE) + start = time + while b'done' not in p.stdout.readline(): + pass + + # Todo: handle sim2seis or pem error + + 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 ['sim2seis']: + if self.true_prim[1][prim_ind] in self.pem_input['vintage']: + result = mat73.loadmat(f'En_{member}/Data_conv.mat')['data_conv'] + v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = np.sum( + np.abs(result[:, :, :, v]), axis=0).flatten() + +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, 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']) + + # 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 = [] + + + 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 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] + + 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 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 + self.pred_data = super().run_fwd_sim(state, member_i, del_folder=True) + + 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. + success = super().call_sim(folder, wait_for_proc) + + 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) + + 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) + + 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) + + 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) + # 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 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') + + if 'compaction' in self.pem_input: + multfactor = self.ecl_case.cell_data('PORV_RC', base_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) + + 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) + + if 'press_conv' in self.pem_input: + pem_input['PRESSURE'] = pem_input['PRESSURE'] * \ + 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=None) + + # 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 + # 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_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. + """ + + def __init__(self, input_dict=None, filename=None, options=None): + 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']) + + # 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 = [] + self.pem_result = [] + self.bar_result = [] + + 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 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] + + 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 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 + self.pred_data = super().run_fwd_sim(state, member_i, del_folder=True) + + 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. + success = super().call_sim(folder, wait_for_proc) + + 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) + + 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) + + 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) + + 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) + # 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 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') + + if 'compaction' in self.pem_input: + multfactor = self.ecl_case.cell_data('PORV_RC', base_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) + + 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) + + if 'press_conv' in self.pem_input: + pem_input['PRESSURE'] = pem_input['PRESSURE'] * \ + 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=None) + + # 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 + # 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 + 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 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 ['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.bar_result[v].flatten() + +class flow_avo(flow_sim2seis): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): + super().__init__(input_dict, filename, options) + + 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() + + 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. + """ + + if member_i >= 0: + folder = 'En_' + str(member_i) + os.sep + if not os.path.exists(folder): + os.mkdir(folder) + else: # XLUO: any negative member_i is considered as the index for the true model + assert 'truth_folder' in self.input_dict, "ensemble member index is negative, please specify " \ + "the folder containing the true model" + if not os.path.exists(self.input_dict['truth_folder']): + os.mkdir(self.input_dict['truth_folder']) + folder = self.input_dict['truth_folder'] + os.sep if self.input_dict['truth_folder'][-1] != os.sep \ + else self.input_dict['truth_folder'] + del_folder = False # never delete this folder + self.folder = folder + self.ensemble_member = member_i + + state['member'] = member_i + + # start by generating the .DATA file, using the .mako template situated in ../folder + self._runMako(folder, state) + success = False + rerun = self.rerun + while rerun >= 0 and not success: + success = self.call_sim(folder, True) + rerun -= 1 + if success: + self.extract_data(member_i) + if del_folder: + if self.saveinfo is not None: # Try to save information + store_ensemble_sim_information(self.saveinfo, member_i) + self.remove_folder(member_i) + return self.pred_data + else: + if hasattr(self, 'redund_sim') and self.redund_sim is not None: + success = self.redund_sim.call_sim(folder, True) + if success: + self.extract_data(member_i) + if del_folder: + if self.saveinfo is not None: # Try to save information + store_ensemble_sim_information(self.saveinfo, member_i) + self.remove_folder(member_i) + return self.pred_data + else: + if del_folder: + self.remove_folder(member_i) + return False + else: + if del_folder: + self.remove_folder(member_i) + return False + + + 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 + + # The field 'run_reservoir_model' can be passed from the method "setup_fwd_run" + if hasattr(self, 'run_reservoir_model'): + run_reservoir_model = self.run_reservoir_model + + if run_reservoir_model is None: + run_reservoir_model = True + + # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. + # Then, get the pem. + if run_reservoir_model: # in case that simulation has already done (e.g., for the true reservoir model) + success = super(flow_sim2seis, self).call_sim(folder, wait_for_proc) + #ecl = ecl_100(filename=self.file) + #ecl.options = self.options + #success = ecl.call_sim(folder, wait_for_proc) + else: + success = True + + if success: + self.ecl_case = ecl.EclipseCase(folder + os.sep + self.file + '.DATA') if folder[-1] != os.sep \ + else ecl.EclipseCase(folder + self.file + '.DATA') + + self.calc_pem() + + 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']) + \ + 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) + #pem_input['PORO'] = np.array(self._reformat3D_then_flatten(multfactor * tmp), dtype=float) + else: + pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float) + #pem_input['PORO'] = np.array(self._reformat3D_then_flatten(tmp), 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) + #pem_input['NTG'] = np.array(self._reformat3D_then_flatten(tmp), dtype=float) + else: + tmp = self.ecl_case.cell_data('NTG') + pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) + #pem_input['NTG'] = np.array(self._reformat3D_then_flatten(tmp), 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) + #pem_input[var] = np.array(self._reformat3D_then_flatten(tmp), dtype=float) + + if 'press_conv' in self.pem_input: + pem_input['PRESSURE'] = pem_input['PRESSURE'] * \ + self.pem_input['press_conv'] + #pem_input['PRESSURE'] = self._reformat3D_then_flatten(pem_input['PRESSURE'] * + # self.pem_input['press_conv']) + + tmp = self.ecl_case.cell_data('PRESSURE', 0) + if hasattr(self.pem, 'p_init'): + P_init = self.pem.p_init * np.ones(tmp.shape)[~tmp.mask] + #P_init = self._reformat3D_then_flatten(self.pem.p_init.reshape(tmp.shape) * np.ones(tmp.shape)) + else: + # initial pressure is first + P_init = np.array(tmp[~tmp.mask], dtype=float) + #P_init = np.array(self._reformat3D_then_flatten(tmp), dtype=float) + + if 'press_conv' in self.pem_input: + P_init = P_init * self.pem_input['press_conv'] + #P_init = self._reformat3D_then_flatten(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] + #saturations = [self._reformat3D_then_flatten(1 - (pem_input['SWAT'] + pem_input['SGAS'])) + # if ph == 'OIL' else pem_input['S{}'.format(ph)] for ph in phases] + + # Get the pressure + 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) + + #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', { + # 'Vp': self.pem.getBulkVel() * .1, 'DIMENS': grid['DIMENS']}, multi_file=False) + #grdecl.write(f'En_{str(self.ensemble_member)}/rho{v + 1}.grdecl', + # {'rho': self.pem.getDens(), 'DIMENS': grid['DIMENS']}, multi_file=False) + + # vp, vs, density + self.vp = (self.pem.getBulkVel() * .1).reshape((self.NX, self.NY, self.NZ), order='F') + self.vs = (self.pem.getShearVel() * .1).reshape((self.NX, self.NY, self.NZ), order='F') + self.rho = (self.pem.getDens()).reshape((self.NX, self.NY, self.NZ), order='F') # in the unit of g/cm^3 + + save_dic = {'vp': self.vp, 'vs': self.vs, 'rho': self.vp} + 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" + else: + if hasattr(self, 'ensemble_member') and (self.ensemble_member is not None) and \ + (self.ensemble_member >= 0): + file_name = folder + os.sep + f"vp_vs_rho_vint{v}.npz" if folder[-1] != os.sep \ + else folder + f"vp_vs_rho_vint{v}.npz" + else: + file_name = os.getcwd() + os.sep + f"vp_vs_rho_vint{v}.npz" + + #with open(file_name, "wb") as f: + # dump(**save_dic, f) + np.savez(file_name, **save_dic) + + # avo data + self._calc_avo_props() + + avo = self.avo_data.flatten(order="F") + + # 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 + + save_dic = {'avo': avo, 'noise_std': noise_std, **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" + else: + # if hasattr(self, 'ensemble_member') and (self.ensemble_member is not None) and \ + # (self.ensemble_member >= 0): + # file_name = folder + os.sep + f"avo_vint{v}.npz" if folder[-1] != os.sep \ + # else folder + f"avo_vint{v}.npz" + # else: + # file_name = os.getcwd() + os.sep + f"avo_vint{v}.npz" + file_name = folder + os.sep + f"avo_vint{v}.npz" if folder[-1] != os.sep \ + else folder + f"avo_vint{v}.npz" + + #with open(file_name, "wb") as f: + # dump(**save_dic, f) + np.savez(file_name, **save_dic) + + return success + + def extract_data(self, member): + # start by getting the data from the flow simulator + super(flow_sim2seis, self).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 '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] + + def _runMako(self, folder, state, addfiles=['properties']): + """ + Hard coding, maybe a better way possible + addfiles: additional files that need to be included into ECLIPSE/OPM DATA file + """ + super()._runMako(folder, state) + + lkup = TemplateLookup(directories=os.getcwd(), input_encoding='utf-8') + for file in addfiles: + tmpl = lkup.get_template('%s.mako' % file) + + # use a context and render onto a file + with open('{0}'.format(folder + file), 'w') as f: + ctx = Context(f, **state) + tmpl.render_context(ctx) + + 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' + 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: + 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 trough the reservoir in vertical direction + cum_time_res = np.cumsum(2 * self.DZ / self.vp, axis=2) + top_res[:, :, np.newaxis] + cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res, top_res[:, :, 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) + + # 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] + + + # from matplotlib import pyplot as plt + # plt.plot(vp_sample[0, 0, :]) + # plt.show() + + #vp_avg = 0.5 * (vp_sample[:, :, 1:] + vp_sample[:, :, :-1]) + #vs_avg = 0.5 * (vs_sample[:, :, 1:] + vs_sample[:, :, :-1]) + #rho_avg = 0.5 * (rho_sample[:, :, 1:] + rho_sample[:, :, :-1]) + + #vp_diff = vp_sample[:, :, 1:] - vp_sample[:, :, :-1] + #vs_diff = vs_sample[:, :, 1:] - vs_sample[:, :, :-1] + #rho_diff = rho_sample[:, :, 1:] - rho_sample[:, :, :-1] + + #R0_smith = 0.5 * (vp_diff / vp_avg + rho_diff / rho_avg) + #G_smith = -2.0 * (vs_avg / vp_avg) ** 2 * (2.0 * vs_diff / vs_avg + rho_diff / rho_avg) + 0.5 * vp_diff / vp_avg + + # PP reflection coefficients, see, e.g., + # "https://pylops.readthedocs.io/en/latest/api/generated/pylops.avo.avo.approx_zoeppritz_pp.html" + # So far, it seems that "approx_zoeppritz_pp" is the only available option + # approx_zoeppritz_pp(vp1, vs1, rho1, vp0, vs0, rho0, theta1) + avo_data_list = [] + + # 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 sereis + 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 vertial 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 + + @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 arraies 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 \ No newline at end of file