diff --git a/simulator/flow_rock.py b/simulator/flow_rock.py index ef7f03c..92b63f5 100644 --- a/simulator/flow_rock.py +++ b/simulator/flow_rock.py @@ -27,6 +27,7 @@ from geostat.decomp import Cholesky from simulator.eclipse import ecl_100 + class flow_rock(flow): """ Couple the OPM-flow simulator with a rock-physics simulator such that both reservoir quantities and petro-elastic @@ -167,16 +168,21 @@ def calc_pem(self, time): #for key in self.all_data_types: - # if 'grav' in key: - # for var in phases: - # # fluid densities - # dens = [var + '_DEN'] - # tmp = self.ecl_case.cell_data(dens, time) - # pem_input[dens] = np.array(tmp[~tmp.mask], dtype=float) - - # # pore volumes at each assimilation step - # tmp = self.ecl_case.cell_data('RPORV', time) - # pem_input['RPORV'] = np.array(tmp[~tmp.mask], dtype=float) + #if 'grav' in key: + for var in phases: + # fluid densities + dens = [var + '_DEN'] + try: + tmp = self.ecl_case.cell_data(dens, time) + pem_input[dens] = np.array(tmp[~tmp.mask], dtype=float) + except: + pem_input[dens] = None + # pore volumes at each assimilation step + try: + tmp = self.ecl_case.cell_data('RPORV', time) + pem_input['RPORV'] = np.array(tmp[~tmp.mask], dtype=float) + except: + pem_input['RPORV'] = None # Get elastic parameters if hasattr(self, 'ensemble_member') and (self.ensemble_member is not None) and \ @@ -649,7 +655,6 @@ def extract_data(self, member): class flow_avo(flow_rock): - def __init__(self, input_dict=None, filename=None, options=None, **kwargs): super().__init__(input_dict, filename, options) @@ -774,7 +779,8 @@ def get_avo_result(self, folder, save_folder): self.calc_velocities(folder, save_folder, grid, v) # avo data - self._calc_avo_props() + #self._calc_avo_props() + self._calc_avo_props_active_cells(grid) avo = self.avo_data.flatten(order="F") @@ -787,8 +793,9 @@ def get_avo_result(self, folder, save_folder): self.calc_velocities(folder, save_folder, grid, -1) # avo data - self._calc_avo_props() - + #self._calc_avo_props() + self._calc_avo_props_active_cells(grid) + avo_baseline = self.avo_data.flatten(order="F") avo = avo - avo_baseline @@ -974,6 +981,7 @@ def _calc_avo_props(self, dt=0.0005): vs_sample[m, l, k] = vs[m, l, idx] rho_sample[m, l, k] = rho[m, l, idx] + # Ricker wavelet wavelet, t_axis, wav_center = ricker(np.arange(0, self.avo_config['wave_len'], dt), f0=self.avo_config['frequency']) @@ -1015,6 +1023,114 @@ def _calc_avo_props(self, dt=0.0005): self.avo_data = avo_data + + def _calc_avo_props_active_cells(self, grid, dt=0.005): + # dt is the fine resolution sampling rate + # convert properties in reservoir model to time domain + vp_shale = self.avo_config['vp_shale'] # scalar value (code may not work for matrix value) + vs_shale = self.avo_config['vs_shale'] # scalar value + rho_shale = self.avo_config['den_shale'] # scalar value + + # Two-way travel time of the top of the reservoir + # TOPS[:, :, 0] corresponds to the depth profile of the reservoir top on the first layer + top_res = 2 * self.TOPS[:, :, 0] / vp_shale + + # Cumulative traveling time through the reservoir in vertical direction + cum_time_res = np.cumsum(2 * self.DZ / self.vp, axis=2) + top_res[:, :, np.newaxis] + + # assumes underburden to be constant. No reflections from underburden. Hence set traveltime to underburden very large + underburden = top_res + np.max(cum_time_res) + + # total travel time + # cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res), axis=2) + cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res, underburden[:, :, np.newaxis]), axis=2) + + # add overburden and underburden of Vp, Vs and Density + vp = np.concatenate((vp_shale * np.ones((self.NX, self.NY, 1)), + self.vp, vp_shale * np.ones((self.NX, self.NY, 1))), axis=2) + vs = np.concatenate((vs_shale * np.ones((self.NX, self.NY, 1)), + self.vs, vs_shale * np.ones((self.NX, self.NY, 1))), axis=2) + #rho = np.concatenate((rho_shale * np.ones((self.NX, self.NY, 1)) * 0.001, # kg/m^3 -> k/cm^3 + # self.rho, rho_shale * np.ones((self.NX, self.NY, 1)) * 0.001), axis=2) + rho = np.concatenate((rho_shale * np.ones((self.NX, self.NY, 1)), + self.rho, rho_shale * np.ones((self.NX, self.NY, 1))), axis=2) + + # get indices of active cells + actnum = np.transpose(grid['ACTNUM'], (2, 1, 0)) + indices = np.where(actnum) + a, b, c = indices + # Combine a and b into a 2D array (each column represents a vector) + ab = np.column_stack((a, b)) + + # Extract unique rows and get the indices of those unique rows + unique_rows, indices = np.unique(ab, axis=0, return_index=True) + + + # search for the lowest grid cell thickness and sample the time according to + # that grid thickness to preserve the thin layer effect + time_sample = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], dt) + if time_sample.shape[0] == 1: + time_sample = time_sample.reshape(-1) + time_sample = np.tile(time_sample, (len(indices), 1)) + + vp_sample = np.zeros((len(indices), time_sample.shape[1])) + vs_sample = np.zeros((len(indices), time_sample.shape[1])) + rho_sample = np.zeros((len(indices), time_sample.shape[1])) + + + for ind in range(len(indices)): + for k in range(time_sample.shape[1]): + # find the right interval of time_sample[m, l, k] belonging to, and use + # this information to allocate vp, vs, rho + idx = np.searchsorted(cum_time[a[indices[ind]], b[indices[ind]], :], time_sample[ind, k], side='left') + idx = idx if idx < len(cum_time[a[indices[ind]], b[indices[ind]], :]) else len(cum_time[a[indices[ind]], b[indices[ind]], :]) - 1 + vp_sample[ind, k] = vp[a[indices[ind]], b[indices[ind]], idx] + vs_sample[ind, k] = vs[a[indices[ind]], b[indices[ind]], idx] + rho_sample[ind, k] = rho[a[indices[ind]], b[indices[ind]], idx] + + + # Ricker wavelet + wavelet, t_axis, wav_center = ricker(np.arange(0, self.avo_config['wave_len'], dt), + f0=self.avo_config['frequency']) + + # Travel time corresponds to reflectivity series + t = time_sample[:, 0:-1] + + # interpolation time + t_interp = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], self.avo_config['t_sampling']) + trace_interp = np.zeros((len(indices), len(t_interp))) + + # number of pp reflection coefficients in the vertical direction + nz_rpp = vp_sample.shape[1] - 1 + conv_op = Convolve1D(nz_rpp, h=wavelet, offset=wav_center, dtype="float32") + + for i in range(len(self.avo_config['angle'])): + angle = self.avo_config['angle'][i] + Rpp = self.pp_func(vp_sample[:, :-1], vs_sample[:, :-1], rho_sample[:, :-1], + vp_sample[:, 1:], vs_sample[:, 1:], rho_sample[:, 1:], angle) + + + + for ind in range(len(indices)): + # convolution with the Ricker wavelet + + w_trace = conv_op * Rpp[ind, :] + + # Sample the trace into regular time interval + f = interp1d(np.squeeze(t[ind, :]), np.squeeze(w_trace), + kind='nearest', fill_value='extrapolate') + trace_interp[ind, :] = f(t_interp) + + if i == 0: + avo_data = trace_interp # 3D + elif i == 1: + avo_data = np.stack((avo_data, trace_interp), axis=-1) # 4D + else: + avo_data = np.concatenate((avo_data, trace_interp[:, :, :, np.newaxis]), axis=3) # 4D + + self.avo_data = avo_data + + @classmethod def _reformat3D_then_flatten(cls, array, flatten=True, order="F"): """ @@ -1323,8 +1439,9 @@ def measurement_locations(self, grid): def find_cell_centers(self, grid): # Find indices where the boolean array is True - indices = np.where(grid['ACTNUM']) - + #indices = np.where(grid['ACTNUM']) + actnum = np.transpose(grid['ACTNUM'], (2, 1, 0)) + indices = np.where(actnum) # `indices` will be a tuple of arrays: (x_indices, y_indices, z_indices) #nactive = len(actind) # Number of active cells @@ -1335,7 +1452,7 @@ def find_cell_centers(self, grid): #N1, N2, N3 = grid['DIMENS'] - c, a, b = indices + b, a, c = indices # Calculate xt, yt, zt xb = 0.25 * (coord[a, b, 0, 0] + coord[a, b + 1, 0, 0] + coord[a + 1, b, 0, 0] + coord[a + 1, b + 1, 0, 0]) yb = 0.25 * (coord[a, b, 0, 1] + coord[a, b + 1, 0, 1] + coord[a + 1, b, 0, 1] + coord[a + 1, b + 1, 0, 1]) @@ -1466,6 +1583,7 @@ def call_sim(self, folder=None, wait_for_proc=False, save_folder=None): folder = self.folder # run flow simulator + #success = True #super(flow_rock, self).call_sim(folder, True) success = super(flow_rock, self).call_sim(folder, True) # use output from flow simulator to forward model gravity response