diff --git a/setup.py b/setup.py index e1bfbf0..4c9e9f9 100644 --- a/setup.py +++ b/setup.py @@ -45,5 +45,6 @@ 'libecalc', 'scikit-learn', 'pylops' + ] + EXTRAS['doc'], ) diff --git a/simulator/flow_rock.py b/simulator/flow_rock.py index 31aa975..ef7f03c 100644 --- a/simulator/flow_rock.py +++ b/simulator/flow_rock.py @@ -18,12 +18,256 @@ # 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 +import sys +#sys.path.append("/home/AD.NORCERESEARCH.NO/mlie/") +from PyGRDECL.GRDECL_Parser import GRDECL_Parser # https://github.com/BinWang0213/PyGRDECL/tree/master from scipy.interpolate import interp1d +from scipy.interpolate import griddata from pipt.misc_tools.analysis_tools import store_ensemble_sim_information from geostat.decomp import Cholesky from simulator.eclipse import ecl_100 +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 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, 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]) + + + #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) + + # 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 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 + 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) + + self.calc_pem(time) + + # mask the bulk imp. to get proper dimensions + tmp_value = np.zeros(self.ecl_case.init.shape) + tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp + self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float, + mask=deepcopy(self.ecl_case.init.mask)) + # run filter + self.pem._filter() + vintage.append(deepcopy(self.pem.bulkimp)) + + if hasattr(self.pem, 'baseline'): # 4D measurement + base_time = dt.datetime(self.startDate['year'], self.startDate['month'], + self.startDate['day']) + dt.timedelta(days=self.pem.baseline) + # + self.calc_pem(base_time) + + # mask the bulk imp. to get proper dimensions + tmp_value = np.zeros(self.ecl_case.init.shape) + + tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp + # kill if values are inf or nan + assert not np.isnan(tmp_value).any() + assert not np.isinf(tmp_value).any() + self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float, + mask=deepcopy(self.ecl_case.init.mask)) + self.pem._filter() + + # 4D response + self.pem_result = [] + for i, elem in enumerate(vintage): + self.pem_result.append(elem - deepcopy(self.pem.bulkimp)) + else: + for i, elem in enumerate(vintage): + self.pem_result.append(elem) + + return success + + def extract_data(self, member): + # start by getting the data from the flow simulator + super().extract_data(member) + + # get the sim2seis from file + for prim_ind in self.l_prim: + # Loop over all keys in pred_data (all data types) + for key in self.all_data_types: + if key in ['bulkimp']: + if self.true_prim[1][prim_ind] in self.pem_input['vintage']: + v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = self.pem_result[v].data.flatten() + + + + class flow_sim2seis(flow): """ Couple the OPM-flow simulator with a sim2seis simulator such that both reservoir quantities and petro-elastic @@ -106,7 +350,7 @@ def call_sim(self, folder=None, wait_for_proc=False): success = super().call_sim(folder, wait_for_proc) if success: - # need a if to check that we have correct sim2seis + # need an if to check that we have correct sim2seis # copy relevant sim2seis files into folder. for file in glob.glob('sim2seis_config/*'): shutil.copy(file, 'En_' + str(self.ensemble_member) + os.sep) @@ -123,7 +367,7 @@ def call_sim(self, folder=None, wait_for_proc=False): time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ dt.timedelta(days=assim_time) - self.calc_pem(time) #mali: update class inhere flor_rock. Include calc_pem as method in flow_rock + self.calc_pem(time) #mali: update class inherent in flow_rock. Include calc_pem as method in flow_rock grdecl.write(f'En_{str(self.ensemble_member)}/Vs{v+1}.grdecl', { 'Vs': self.pem.getShearVel()*.1, 'DIMENS': grid['DIMENS']}, multi_file=False) @@ -160,10 +404,12 @@ def extract_data(self, member): self.pred_data[prim_ind][key] = np.sum( np.abs(result[:, :, :, v]), axis=0).flatten() -class flow_rock(flow): +class flow_barycenter(flow): """ Couple the OPM-flow simulator with a rock-physics simulator such that both reservoir quantities and petro-elastic - quantities can be calculated. Inherit the flow class, and use super to call similar functions. + quantities can be calculated. Inherit the flow class, and use super to call similar functions. In the end, the + barycenter and moment of interia for the bulkimpedance objects, are returned as observations. The objects are + identified using k-means clustering, and the number of objects are determined using and elbow strategy. """ def __init__(self, input_dict=None, filename=None, options=None): @@ -187,6 +433,8 @@ def __init__(self, input_dict=None, filename=None, options=None): self.saveinfo = [input_dict['savesiminfo']] self.scale = [] + self.pem_result = [] + self.bar_result = [] def _getpeminfo(self, input_dict): """ @@ -217,6 +465,8 @@ def _getpeminfo(self, input_dict): self.pem_input['overburden'] = elem[1] if elem[0] == 'percentile': # use for scaling self.pem_input['percentile'] = elem[1] + if elem[0] == 'clusters': # number of clusters for each barycenter + self.pem_input['clusters'] = elem[1] pem = getattr(import_module('simulator.rockphysics.' + self.pem_input['model'].split()[0]), self.pem_input['model'].split()[1]) @@ -358,263 +608,37 @@ def call_sim(self, folder=None, wait_for_proc=False): self.pem._filter() # 4D response - self.pem_result = [] for i, elem in enumerate(vintage): self.pem_result.append(elem - deepcopy(self.pem.bulkimp)) else: for i, elem in enumerate(vintage): self.pem_result.append(elem) + # Extract k-means centers and interias for each element in pem_result + if 'clusters' in self.pem_input: + npzfile = np.load(self.pem_input['clusters'], allow_pickle=True) + n_clusters_list = npzfile['n_clusters_list'] + npzfile.close() + else: + n_clusters_list = len(self.pem_result)*[2] + kmeans_kwargs = {"init": "random", "n_init": 10, "max_iter": 300, "random_state": 42} + for i, bulkimp in enumerate(self.pem_result): + std = np.std(bulkimp) + features = np.argwhere(np.squeeze(np.reshape(np.abs(bulkimp), self.ecl_case.init.shape,)) > 3 * std) + scaler = StandardScaler() + scaled_features = scaler.fit_transform(features) + kmeans = KMeans(n_clusters=n_clusters_list[i], **kmeans_kwargs) + kmeans.fit(scaled_features) + kmeans_center = np.squeeze(scaler.inverse_transform(kmeans.cluster_centers_)) # data / measurements + self.bar_result.append(np.append(kmeans_center, kmeans.inertia_)) + return success def extract_data(self, member): # start by getting the data from the flow simulator super().extract_data(member) - # get the sim2seis from file - 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 + # 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: @@ -622,9 +646,10 @@ 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): + +class flow_avo(flow_rock): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): super().__init__(input_dict, filename, options) @@ -634,7 +659,7 @@ def __init__(self, input_dict=None, filename=None, options=None, **kwargs): def setup_fwd_run(self, **kwargs): self.__dict__.update(kwargs) - super().setup_fwd_run() + super().setup_fwd_run(redund_sim=None) def run_fwd_sim(self, state, member_i, del_folder=True): """ @@ -702,7 +727,6 @@ def run_fwd_sim(self, state, member_i, del_folder=True): 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 @@ -719,7 +743,7 @@ def call_sim(self, folder=None, wait_for_proc=False, run_reservoir_model=None, s # 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) + success = super(flow_rock, 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) @@ -727,91 +751,114 @@ def call_sim(self, folder=None, wait_for_proc=False, run_reservoir_model=None, s 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() + self.get_avo_result(folder, save_folder) - #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" + return success - #with open(file_name, "wb") as f: - # dump(**save_dic, f) - np.savez(file_name, **save_dic) + def get_avo_result(self, folder, save_folder): + 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 + self.calc_velocities(folder, save_folder, grid, v) + + # avo data + self._calc_avo_props() + + avo = self.avo_data.flatten(order="F") + + # MLIE: implement 4D avo + if 'baseline' in self.pem_input: # 4D measurement + base_time = dt.datetime(self.startDate['year'], self.startDate['month'], + self.startDate['day']) + dt.timedelta(days=self.pem_input['baseline']) + self.calc_pem(base_time) + # vp, vs, density in reservoir + self.calc_velocities(folder, save_folder, grid, -1) # avo data self._calc_avo_props() + + avo_baseline = self.avo_data.flatten(order="F") + avo = avo - avo_baseline + + + # 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 - 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" - 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) + + def calc_velocities(self, folder, save_folder, grid, v): + # The properties in pem are 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.full(grid['DIMENS'], self.avo_config['vp_shale']) + self.vp[true_indices] = (self.pem.getBulkVel()) + self.vs = np.full(grid['DIMENS'], self.avo_config['vs_shale']) + self.vs[true_indices] = (self.pem.getShearVel()) + self.rho = np.full(grid['DIMENS'], self.avo_config['den_shale']) + self.rho[true_indices] = (self.pem.getDens()) + else: + self.vp = (self.pem.getBulkVel()).reshape((self.NX, self.NY, self.NZ), order='F') + self.vs = (self.pem.getShearVel()).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) + # 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) + super(flow_rock, self).extract_data(member) # get the sim2seis from file for prim_ind in self.l_prim: @@ -825,105 +872,6 @@ def extract_data(self, member): 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 @@ -982,18 +930,27 @@ def _calc_avo_props(self, dt=0.0005): # 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 + # 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, top_res[:, :, np.newaxis]), axis=2) + # 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)) * 0.001, # kg/m^3 -> k/cm^3 + # self.rho, rho_shale * np.ones((self.NX, self.NY, 1)) * 0.001), axis=2) + rho = np.concatenate((rho_shale * np.ones((self.NX, self.NY, 1)), + self.rho, rho_shale * np.ones((self.NX, self.NY, 1))), axis=2) # search for the lowest grid cell thickness and sample the time according to # that grid thickness to preserve the thin layer effect @@ -1017,40 +974,20 @@ 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] - - # 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 + + # Travel time corresponds to reflectivity series t = time_sample[:, :, 0:-1] # interpolation time t_interp = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], self.avo_config['t_sampling']) trace_interp = np.zeros((self.NX, self.NY, len(t_interp))) - # number of pp reflection coefficients in the vertial direction + # number of pp reflection coefficients in the vertical direction + nz_rpp = vp_sample.shape[2] - 1 for i in range(len(self.avo_config['angle'])): @@ -1096,3 +1033,469 @@ def _reformat3D_then_flatten(cls, array, flatten=True, order="F"): return new_array else: return array + +class flow_grav(flow_rock): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): + super().__init__(input_dict, filename, options) + + self.grav_input = {} + assert 'grav' in input_dict, 'To do GRAV simulation, please specify an "GRAV" section in the "FWDSIM" part' + self._get_grav_info() + + def setup_fwd_run(self, **kwargs): + self.__dict__.update(kwargs) + + super().setup_fwd_run(redund_sim=None) + + def run_fwd_sim(self, state, member_i, del_folder=True): + # The inherited simulator also has a run_fwd_sim. Call this. + #self.ensemble_member = member_i + #self.pred_data = super().run_fwd_sim(state, member_i, del_folder=True) + + #return self.pred_data + + 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, save_folder=None): + # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. + # Then, get the pem. + if folder is None: + folder = self.folder + + # run flow simulator + success = super(flow_rock, self).call_sim(folder, True) + # + # use output from flow simulator to forward model gravity response + if success: + self.get_grav_result(folder, save_folder) + + return success + + def get_grav_result(self, folder, save_folder): + 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() + + # cell centers + self.find_cell_centers(grid) + + # receiver locations + self.measurement_locations(grid) + + # loop over vintages with gravity acquisitions + grav_struct = {} + + for v, assim_time in enumerate(self.grav_config['vintage']): + time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ + dt.timedelta(days=assim_time) + + # porosity, saturation, densities, and fluid mass at individual time-steps + grav_struct[v] = self.calc_mass(time) # calculate the mass of each fluid in each grid cell + + # TODO save densities, saturation and mass for each vintage for plotting? + # 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) + if 'baseline' in self.grav_config: # 4D measurement + base_time = dt.datetime(self.startDate['year'], self.startDate['month'], + self.startDate['day']) + dt.timedelta(days=self.grav_config['baseline']) + # porosity, saturation, densities, and fluid mass at time of baseline survey + grav_base = self.calc_mass(base_time) + + + else: + # seafloor gravity only work in 4D mode + print('Need to specify Baseline survey in pipt file') + + vintage = [] + + for v, assim_time in enumerate(self.grav_config['vintage']): + dg = self.calc_grav(grid, grav_base, grav_struct[v]) + vintage.append(deepcopy(dg)) + + save_dic = {'grav': dg, **self.grav_config} + if save_folder is not None: + file_name = save_folder + os.sep + f"grav_vint{v}.npz" if save_folder[-1] != os.sep \ + else save_folder + f"grav_vint{v}.npz" + else: + file_name = folder + os.sep + f"grav_vint{v}.npz" if folder[-1] != os.sep \ + else folder + f"grav_vint{v}.npz" + + # with open(file_name, "wb") as f: + # dump(**save_dic, f) + np.savez(file_name, **save_dic) + + # 4D response + self.grav_result = [] + for i, elem in enumerate(vintage): + self.grav_result.append(elem) + + def calc_mass(self, time): + # fluid phases written to restart file from simulator run + phases = self.ecl_case.init.phases + + grav_input = {} + # get active porosity + # pore volumes at each assimilation step + tmp = self.ecl_case.cell_data('RPORV', time) + grav_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float) + + # extract saturation + if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: # This should be extended + for var in phases: + if var in ['WAT', 'GAS']: + tmp = self.ecl_case.cell_data('S{}'.format(var), time) + grav_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) + + grav_input['SOIL'] = 1 - (grav_input['SWAT'] + grav_input['SGAS']) + + 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) + grav_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) + + grav_input['SOIL'] = 1 - (grav_input['SGAS']) + + else: + print('Type and number of fluids are unspecified in calc_mass') + + + + # fluid densities + for var in phases: + dens = var + '_DEN' + tmp = self.ecl_case.cell_data(dens, time) + grav_input[dens] = np.array(tmp[~tmp.mask], dtype=float) + + + #fluid masses + for var in phases: + mass = var + '_mass' + grav_input[mass] = grav_input[var + '_DEN'] * grav_input['S' + var] * grav_input['PORO'] + + return grav_input + + def calc_grav(self, grid, grav_base, grav_repeat): + + + + #cell_centre = [x, y, z] + cell_centre = self.grav_config['cell_centre'] + x = cell_centre[0] + y = cell_centre[1] + z = cell_centre[2] + + pos = self.grav_config['meas_location'] + + # Initialize dg as a zero array, with shape depending on the condition + # assumes the length of each vector gives the total number of measurement points + N_meas = (len(pos['x'])) + dg = np.zeros(N_meas) # 1D array for dg + + # total fluid mass at this time + phases = self.ecl_case.init.phases + if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: + dm = grav_repeat['OIL_mass'] + grav_repeat['WAT_mass'] + grav_repeat['GAS_mass'] - (grav_base['OIL_mass'] + grav_base['WAT_mass'] + grav_base['GAS_mass']) + + elif 'OIL' in phases and 'GAS' in phases: # Smeaheia model + dm = grav_repeat['OIL_mass'] + grav_repeat['GAS_mass'] - (grav_base['OIL_mass'] + grav_base['GAS_mass']) + + else: + print('Type and number of fluids are unspecified in calc_grav') + + + for j in range(N_meas): + + # Calculate dg for the current measurement location (j, i) + dg_tmp = (z - pos['z'][j]) / ((x[j] - pos['x'][j]) ** 2 + (y[j] - pos['y'][j]) ** 2 + ( + z - pos['z'][j]) ** 2) ** (3 / 2) + + dg[j] = np.dot(dg_tmp, dm) + print(f'Progress: {j + 1}/{N_meas}') # Mimicking waitbar + + # Scale dg by the constant + dg *= 6.67e-3 + + return dg + + def measurement_locations(self, grid): + # Determine the size of the target area as defined by the reservoir area + + #cell_centre = [x, y, z] + cell_centre = self.grav_config['cell_centre'] + xmin = np.min(cell_centre[0]) + xmax = np.max(cell_centre[0]) + ymin = np.min(cell_centre[1]) + ymax = np.max(cell_centre[1]) + + # Make a mesh of the area + pad = self.grav_config.get('padding_reservoir', 3000) # 3 km padding around the reservoir + if 'padding_reservoir' not in self.grav_config: + print('Please specify extent of measurement locations (Padding in pipt file), using 3 km as default') + + xmin -= pad + xmax += pad + ymin -= pad + ymax += pad + + xspan = xmax - xmin + yspan = ymax - ymin + + dxy = self.grav_config.get('grid_spacing', 1500) # + if 'grid_spacing' not in self.grav_config: + print('Please specify grid spacing in pipt file, using 1.5 km as default') + + Nx = int(np.ceil(xspan / dxy)) + Ny = int(np.ceil(yspan / dxy)) + + xvec = np.linspace(xmin, xmax, Nx) + yvec = np.linspace(ymin, ymax, Ny) + + x, y = np.meshgrid(xvec, yvec) + + pos = {'x': x.flatten(), 'y': y.flatten()} + + # Handle seabed map or water depth scalar if defined in pipt + if 'seabed' in self.grav_config and self.grav_config['seabed'] is not None: + pos['z'] = griddata((self.grav_config['seabed']['x'], self.grav_config['seabed']['y']), + self.grav_config['seabed']['z'], (pos['x'], pos['y']), method='nearest') + else: + pos['z'] = np.ones_like(pos['x']) * self.grav_config.get('water_depth', 300) + + if 'water_depth' not in self.grav_config: + print('Please specify water depths in pipt file, using 300 m as default') + + #return pos + self.grav_config['meas_location'] = pos + + def find_cell_centers(self, grid): + + # Find indices where the boolean array is True + indices = np.where(grid['ACTNUM']) + + # `indices` will be a tuple of arrays: (x_indices, y_indices, z_indices) + #nactive = len(actind) # Number of active cells + + coord = grid['COORD'] + zcorn = grid['ZCORN'] + + # Unpack dimensions + #N1, N2, N3 = grid['DIMENS'] + + + c, a, b = indices + # Calculate xt, yt, zt + xb = 0.25 * (coord[a, b, 0, 0] + coord[a, b + 1, 0, 0] + coord[a + 1, b, 0, 0] + coord[a + 1, b + 1, 0, 0]) + yb = 0.25 * (coord[a, b, 0, 1] + coord[a, b + 1, 0, 1] + coord[a + 1, b, 0, 1] + coord[a + 1, b + 1, 0, 1]) + zb = 0.25 * (coord[a, b, 0, 2] + coord[a, b + 1, 0, 2] + coord[a + 1, b, 0, 2] + coord[a + 1, b + 1, 0, 2]) + + xt = 0.25 * (coord[a, b, 1, 0] + coord[a, b + 1, 1, 0] + coord[a + 1, b, 1, 0] + coord[a + 1, b + 1, 1, 0]) + yt = 0.25 * (coord[a, b, 1, 1] + coord[a, b + 1, 1, 1] + coord[a + 1, b, 1, 1] + coord[a + 1, b + 1, 1, 1]) + zt = 0.25 * (coord[a, b, 1, 2] + coord[a, b + 1, 1, 2] + coord[a + 1, b, 1, 2] + coord[a + 1, b + 1, 1, 2]) + + # Calculate z, x, and y positions + z = (zcorn[c, 0, a, 0, b, 0] + zcorn[c, 0, a, 1, b, 0] + zcorn[c, 0, a, 0, b, 1] + zcorn[c, 0, a, 1, b, 1] + + zcorn[c, 1, a, 0, b, 0] + zcorn[c, 1, a, 1, b, 0] + zcorn[c, 1, a, 0, b, 1] + zcorn[c, 1, a, 1, b, 1]) / 8 + + x = xb + (xt - xb) * (z - zb) / (zt - zb) + y = yb + (yt - yb) * (z - zb) / (zt - zb) + + + cell_centre = [x, y, z] + self.grav_config['cell_centre'] = cell_centre + + def _get_grav_info(self, grav_config=None): + """ + GRAV configuration + """ + # list of configuration parameters in the "Grav" section of teh pipt file + config_para_list = ['baseline', 'vintage', 'water_depth', 'padding', 'grid_spacing'] + + if 'grav' in self.input_dict: + self.grav_config = {} + for elem in self.input_dict['grav']: + assert elem[0] in config_para_list, f'Property {elem[0]} not supported' + self.grav_config[elem[0]] = elem[1] + + + else: + self.grav_config = None + + def extract_data(self, member): + # start by getting the data from the flow simulator + super(flow_rock, self).extract_data(member) + + # get the gravity data from results + for prim_ind in self.l_prim: + # Loop over all keys in pred_data (all data types) + for key in self.all_data_types: + if 'grav' in key: + if self.true_prim[1][prim_ind] in self.grav_config['vintage']: + v = self.grav_config['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = self.grav_result[v].flatten() + +class flow_grav_and_avo(flow_avo, flow_grav): + def __init__(self, input_dict=None, filename=None, options=None, **kwargs): + super().__init__(input_dict, filename, options) + + self.grav_input = {} + assert 'grav' in input_dict, 'To do GRAV simulation, please specify an "GRAV" section in the "FWDSIM" part' + self._get_grav_info() + + assert 'avo' in input_dict, 'To do AVO simulation, please specify an "AVO" section in the "FWDSIM" part' + self._get_avo_info() + + def setup_fwd_run(self, **kwargs): + self.__dict__.update(kwargs) + + super().setup_fwd_run(redund_sim=None) + + def run_fwd_sim(self, state, member_i, del_folder=True): + # The inherited simulator also has a run_fwd_sim. Call this. + #self.ensemble_member = member_i + #self.pred_data = super().run_fwd_sim(state, member_i, del_folder=True) + + #return self.pred_data + + 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, save_folder=None): + # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. + # Then, get the pem. + if folder is None: + folder = self.folder + + # run flow simulator + success = super(flow_rock, self).call_sim(folder, True) + + # use output from flow simulator to forward model gravity response + if success: + # calculate gravity data based on flow simulation output + self.get_grav_result(folder, save_folder) + # calculate avo data based on flow simulation output + self.get_avo_result(folder, save_folder) + + return success + + + def extract_data(self, member): + # start by getting the data from the flow simulator i.e. prod. and inj. data + super(flow_rock, self).extract_data(member) + + # get the gravity data from results + for prim_ind in self.l_prim: + # Loop over all keys in pred_data (all data types) + for key in self.all_data_types: + if 'grav' in key: + if self.true_prim[1][prim_ind] in self.grav_config['vintage']: + v = self.grav_config['vintage'].index(self.true_prim[1][prim_ind]) + self.pred_data[prim_ind][key] = self.grav_result[v].flatten() + + if '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] + diff --git a/simulator/rockphysics/softsandrp.py b/simulator/rockphysics/softsandrp.py new file mode 100644 index 0000000..974e2e5 --- /dev/null +++ b/simulator/rockphysics/softsandrp.py @@ -0,0 +1,616 @@ +"""Descriptive description.""" + +__author__ = {'TM', 'TB', 'ML'} + +# standardrp.py +import numpy as np +import sys +import multiprocessing as mp + +from numpy.random import poisson + +# internal load +from misc.system_tools.environ_var import OpenBlasSingleThread # Single threaded OpenBLAS runs + + +class elasticproperties: + """ + Calculate elastic properties from standard + rock-physics models, specifically following Batzle + and Wang, Geophysics, 1992, for fluid properties, and + Report 1 in Abul Fahimuddin's thesis at Universty of + Bergen (2010) for other properties. + + Example + ------- + >>> porosity = 0.2 + ... pressure = 5 + ... phases = ["Oil","Water"] + ... saturations = [0.3, 0.5] + ... + ... satrock = Elasticproperties() + ... satrock.calc_props(phases, saturations, pressure, porosity) + """ + + def __init__(self, input_dict): + self.dens = None + self.bulkmod = None + self.shearmod = None + self.bulkvel = None + self.shearvel = None + self.bulkimp = None + self.shearimp = None + # The overburden for each grid cell must be + # specified as values on an .npz-file whose + # name is given in input_dict. + self.input_dict = input_dict + self._extInfoInputDict() + + def _extInfoInputDict(self): + # The key word for the file name in the + # dictionary must read "overburden" + if 'overburden' in self.input_dict: + obfile = self.input_dict['overburden'] + npzfile = np.load(obfile) + # The values of overburden must have been + # stored on file using: + # np.savez(, + # obvalues=) + self.overburden = npzfile['obvalues'] + npzfile.close() + #else: + # # Norne litho pressure equation in Bar + # P_litho = -49.6 + 0.2027 * Z + 6.127e-6 * Z ** 2 # Using e-6 for scientific notation + # # Convert reservoir pore pressure from Bar to MPa + # P_litho *= 0.1 + # self.overburden = P_litho + + if 'baseline' in self.input_dict: + self.baseline = self.input_dict['baseline'] # 4D baseline + if 'parallel' in self.input_dict: + self.parallel = self.input_dict['parallel'] + + def _filter(self): + bulkmod = self.bulkimp + self.bulkimp = bulkmod.flatten() + + def setup_fwd_run(self, state): + """ + Setup the input parameters to be used in the PEM simulator. Parameters can be an ensemble or a single array. + State is set as an attribute of the simulator, and the correct value is determined in self.pem.calc_props() + + Parameters + ---------- + state : dict + Dictionary of input parameters or states. + + Changelog + --------- + - KF 11/12-2018 + """ + # self.inv_state = {} + # list_pem_param =[el for el in [foo for foo in self.pem['garn'].keys()] + [foo for foo in self.filter.keys()] + + # [foo for foo in self.__dict__.keys()]] + + # list_tot_param = state.keys() + # for param in list_tot_param: + # if param in list_pem_param or (param.split('_')[-1] in ['garn', 'rest']): + # self.inv_state[param] = state[param] + + pass + + def calc_props(self, phases, saturations, pressure, + porosity, dens = None, wait_for_proc=None, ntg=None, Rs=None, press_init=None, ensembleMember=None): + ### + # TODO add fluid densities here -needs to be added as optional input + # + if not isinstance(phases, list): + phases = [phases] + if not isinstance(saturations, list): + saturations = [saturations] + if not isinstance(pressure, list) and \ + type(pressure).__module__ != 'numpy': + pressure = [pressure] + if not isinstance(porosity, list) and \ + type(porosity).__module__ != 'numpy': + porosity = [porosity] + # + # Load "overburden" pressures into local variable to + # comply with remaining code parts + poverburden = self.overburden + if press_init is None: + p_init = self.p_init + else: + p_init = press_init + + # Average number of contacts that each grain has with surrounding grains + coordnumber = self._coordination_number() + + # porosity value separating the porous media's mechanical and acoustic behaviour + phicritical = self._critical_porosity() + + + # Check that no. of phases is equal to no. of + # entries in saturations list + # + assert (len(saturations) == len(phases)) + # + # Make saturation a Numpy array (so that we + # can easily access the values for each + # phase at one grid cell) + # + # Transpose makes it a no. grid cells x phases + # array + saturations = np.array(saturations).T + # + # Check if we actually inputted saturation values + # for a single grid cell. If yes, we redefine + # saturations to get it on the correct form (no. + # grid cells x phases array). + # + if saturations.ndim == 1: + saturations = \ + np.array([[x] for x in saturations]).T + # + # Loop over all grid cells and calculate the + # various saturated properties + # + self.phases = phases + + self.dens = np.zeros(len(saturations[:, 0])) + self.bulkmod = np.zeros(len(saturations[:, 0])) + self.shearmod = np.zeros(len(saturations[:, 0])) + self.bulkvel = np.zeros(len(saturations[:, 0])) + self.shearvel = np.zeros(len(saturations[:, 0])) + self.bulkimp = np.zeros(len(saturations[:, 0])) + self.shearimp = np.zeros(len(saturations[:, 0])) + + if ntg is None: + ntg = [None for _ in range(len(saturations[:, 0]))] + if Rs is None: + Rs = [None for _ in range(len(saturations[:, 0]))] + if p_init is None: + p_init = [None for _ in range(len(saturations[:, 0]))] + + for i in range(len(saturations[:, 0])): + # + # Calculate fluid properties + # + # set Rs if needed + densf, bulkf = \ + self._fluidprops(self.phases, + saturations[i, :], pressure[i], Rs[i]) + # + #denss, bulks, shears = \ + # self._solidprops(porosity[i], ntg[i], i) + + # + denss, bulks, shears = \ + self._solidprops_Johansen() + # + # Calculate dry rock moduli + # + + #bulkd, sheard = \ + # self._dryrockmoduli(porosity[i], + # overburden[i], + # pressure[i], bulks, + # shears, i, ntg[i], p_init[i], denss, Rs[i], self.phases) + # + peff = self._effective_pressure(poverburden[i], pressure[i]) + + + bulkd, sheard = \ + self._dryrockmoduli_Smeaheia(coordnumber, phicritical, porosity[i], peff, bulks, shears) + + # ------------------------------- + # Calculate saturated properties + # ------------------------------- + # + # Density (kg/m3) + # + self.dens[i] = (porosity[i]*densf + + (1-porosity[i])*denss) + # + # Moduli (MPa) + # + self.bulkmod[i] = \ + bulkd + (1 - bulkd/bulks)**2 / \ + (porosity[i]/bulkf + + (1-porosity[i])/bulks - + bulkd/(bulks**2)) + self.shearmod[i] = sheard + # + # Velocities (km/s) + # + self.bulkvel[i] = \ + np.sqrt((abs(self.bulkmod[i]) + + 4*self.shearmod[i]/3)/(self.dens[i])) + self.shearvel[i] = \ + np.sqrt(self.shearmod[i] / + (self.dens[i])) + # + # convert from (km/s) to (m/s) + # + self.bulkvel[i] *= 1000 + self.shearvel[i] *= 1000 + # + # Impedance (m/s)*(kg/m3) + # + self.bulkimp[i] = self.dens[i] * \ + self.bulkvel[i] + self.shearimp[i] = self.dens[i] * \ + self.shearvel[i] + + def getMatchProp(self, petElProp): + if petElProp.lower() == 'density': + self.match_prop = self.getDens() + elif petElProp.lower() == 'bulk_modulus': + self.match_prop = self.getBulkMod() + elif petElProp.lower() == 'shear_modulus': + self.match_prop = self.getShearMod() + elif petElProp.lower() == 'bulk_velocity': + self.match_prop = self.getBulkVel() + elif petElProp.lower() == 'shear_velocity': + self.match_prop = self.getShearVel() + elif petElProp.lower() == "bulk_impedance": + self.match_prop = self.getBulkImp() + elif petElProp.lower() == 'shear_impedance': + self.match_prop = self.getShearImp() + else: + print("\nError in getMatchProp method") + print("No model output type selected for " + "data match.") + print("Legal model output types are " + "(case insensitive):") + print("Density, bulk modulus, shear " + "modulus, bulk velocity,") + print("shear velocity, bulk impedance, " + "shear impedance") + sys.exit(1) + return self.match_prop + + def getDens(self): + return self.dens + + def getBulkMod(self): + return self.bulkmod + + def getShearMod(self): + return self.shearmod + + def getBulkVel(self): + return self.bulkvel + + def getShearVel(self): + return self.shearvel + + def getBulkImp(self): + return self.bulkimp + + def getShearImp(self): + return self.shearimp + + # + # =================================================== + # Fluid properties start + # =================================================== + # + def _fluidprops(self, fphases, fsats, fpress, Rs=None): + # + # Calculate fluid density and bulk modulus + # + # + # Input + # fphases - fluid phases present; Oil + # and/or Water and/or Gas + # fsats - fluid saturation values for + # fluid phases in "fphases" + # fpress - fluid pressure value (MPa) + # Rs - Gas oil ratio. Default value None + + # + # Output + # fdens - density of fluid mixture for + # pressure value "fpress" inherited + # from phaseprops) + # fbulk - bulk modulus of fluid mixture for + # pressure value "fpress" (unit + # inherited from phaseprops) + # + # ----------------------------------------------- + # + fdens = 0.0 + fbinv = 0.0 + + for i in range(len(fphases)): + # + # Calculate mixture properties by summing + # over individual phase properties + # + pdens, pbulk = self._phaseprops(fphases[i], + fpress, Rs) + fdens = fdens + fsats[i]*abs(pdens) + fbinv = fbinv + fsats[i]/abs(pbulk) + fbulk = 1.0/fbinv + # + return fdens, fbulk + # + # --------------------------------------------------- + # + + def _phaseprops(self, fphase, press, Rs=None): + # + # Calculate properties for a single fluid phase + # + # + # Input + # fphase - fluid phase; Oil, Water or Gas + # press - fluid pressure value (MPa) + # + # Output + # pdens - phase density of fluid phase + # "fphase" for pressure value + # "press" (kg/m³) + # pbulk - bulk modulus of fluid phase + # "fphase" for pressure value + # "press" (MPa) + # + # ----------------------------------------------- + # + if fphase.lower() == "oil": + coeffsrho = np.array([0.8, 829.9]) + coeffsbulk = np.array([10.42, 995.79]) + elif fphase.lower() == "wat": + coeffsrho = np.array([0.3, 1067.3]) + coeffsbulk = np.array([9.0, 2807.6]) + elif fphase.lower() == "gas": + coeffsrho = np.array([4.7, 13.4]) + coeffsbulk = np.array([2.75, 0.0]) + else: + print("\nError in phaseprops method") + print("Illegal fluid phase name.") + print("Legal fluid phase names are (case " + "insensitive): Oil, Wat, and Gas.") + sys.exit(1) + # + # Assume simple linear pressure dependencies. + # Coefficients are inferred from + # plots in Batzle and Wang, Geophysics, 1992, + # (where I set the temperature to be 100 degrees + # Celsius, Note also that they give densities in + # g/cc). The resulting straight lines do not fit + # the data extremely well, but they should + # be sufficiently accurate for the purpose of + # this project. + # + pdens = coeffsrho[0]*press + coeffsrho[1] + pbulk = coeffsbulk[0]*press + coeffsbulk[1] + # + return pdens, pbulk + + # + # ======================= + # Fluid properties end + # ======================= + # + + # + # ========================= + # Solid properties start + # ========================= + # + def _solidprops_Johansen(self): + # + # Calculate bulk and shear solid rock (mineral) + # moduli by averaging Hashin-Shtrikman bounds + # + # + # Input + # poro -porosity + # + # Output + # denss - solid rock density (kg/m³) + # bulks - solid rock bulk modulus (unit MPa) + # shears - solid rock shear modulus (unit MPa) + # + # ----------------------------------------------- + # + # + # Solid rock (mineral) density. (Note + # that this is often termed \rho_dry, and not + # \rho_s) + + denss = 2650 # Density of mineral/solid rock kg/m3 + + # + bulks = 37 # (GPa) + shears = 44 # (GPa) + bulks *= 1000 # Convert from GPa to MPa + shears *= 1000 + # + return denss, bulks, shears + # + # + # ======================= + # Solid properties end + # ======================= + # + def _coordination_number(self): + # Applies for granular media + # Average number of contacts that each grain has with surrounding grains + # Coordnumber = 6; simple cubic packing + # Coordnumber = 12; hexagonal close packing + # Needed for Hertz-Mindlin model + # Smeaheia number (Tuhin) + coordnumber = 9 + + return coordnumber + # + def _critical_porosity(self): + # For most porous media there exists a critical porosity + # phi_critical, that seperates their mechanical and acoustic behaviour into two domains. + # For porosities below phi_critical the mineral grains are oad bearing, for values above the grains are + # suspended in the fluids which are load-bearing + # Needed for Hertz-Mindlin model + # Smeaheia number (Tuhin) + phicritical = 0.36 + + return phicritical + # + def _effective_pressure(self, poverb, pfluid): + + # Input + # poverb - overburden pressure (MPa) + # pfluid - fluid pressure (MPa) + + peff = poverb - pfluid + + if peff < 0: + print("\nError in _hertzmindlin method") + print("Negative effective pressure (" + str(peff) + + "). Setting effective pressure to 0.01") + peff = 0.01 + + + + return peff + + # ============================ + # Dry rock properties start + # ============================ + # + def _dryrockmoduli_Smeaheia(self, coordnumber, phicritical, poro, peff, bulks, shears): + # + # + # Calculate bulk and shear dry rock moduli, + + # + # Input + # poro - porosity + # peff - effective pressure overburden - fluid pressure (MPa) + # bulks - bulk solid (mineral) rock bulk + # modulus (MPa) + # shears - solid rock (mineral) shear + # modulus (MPa) + # + # Output + # bulkd - dry rock bulk modulus (unit + # inherited from hertzmindlin and + # variable; bulks) + # sheard - dry rock shear modulus (unit + # inherited from hertzmindlin and + # variable; shears) + # + # ----------------------------------------------- + # + # Calculate Hertz-Mindlin moduli + # + bulkhm, shearhm = self._hertzmindlin_Mavko(peff, bulks, shears, coordnumber, phicritical) + # + bulkd = 1 / ((poro / phicritical) / (bulkhm + 4 / 3 * shearhm) + + (1 - poro / phicritical) / (bulks + 4 / 3 * shearhm)) - 4 / 3 * shearhm + + psi = (9 * bulkhm + 8 * shearhm) / (bulkhm + 2 * shearhm) + + sheard = 1 / ((poro / phicritical) / (shearhm + 1 / 6 * psi * shearhm) + + (1 - poro / phicritical) / (shears + 1 / 6 * psi * shearhm)) - 1 / 6 * psi * shearhm + + #return K_dry, G_dry + return bulkd, sheard + + + # + # --------------------------------------------------- + # + + def _hertzmindlin_Mavko(self, peff, bulks, shears, coordnumber, phicritical): + # + # Calculate bulk and shear Hertz-Mindlin moduli + # adapted from Tuhins kode and "The rock physics handbook", pp247 + # + # + # Input + # p_eff - effective pressure + # bulks - bulk solid (mineral) rock bulk + # modulus (MPa) + # shears - solid rock (mineral) shear + # modulus (MPa) + # coordnumber - average number of contacts that each grain has with surrounding grains + # phicritical - critical porosity + # + # Output + # bulkhm - Hertz-Mindlin bulk modulus + # (MPa) + # shearhm - Hertz-Mindlin shear modulus + # (MPa) + # + # ----------------------------------------------- + # + + + poisson = (3 * bulks - 2 * shears) / (6 * bulks + 2 * shears) + + bulkhm = ((coordnumber ** 2 * (1 - phicritical) ** 2 * shears ** 2 * peff) / + (18 * np.pi ** 2 * (1 - poisson) ** 2)) ** (1 / 3) + shearhm = (5 - 4 * poisson) / (10 - 5 * poisson) * \ + ((3 * coordnumber ** 2 * (1 - phicritical) ** 2 * shears ** 2 * peff) / + (2 * np.pi ** 2 * (1 - poisson) ** 2)) ** (1 / 3) + + + + # + return bulkhm, shearhm + + + # =========================== + # Dry rock properties end + # =========================== + + +if __name__ == '__main__': + # + # Example input with two phases and three grid cells + # + porosity = [0.34999999, 0.34999999, 0.34999999] +# pressure = [ 29.29150963, 29.14003944, 28.88845444] + pressure = [29.3558, 29.2625, 29.3558] +# pressure = [ 25.0, 25.0, 25.0] + phases = ["Oil", "Wat"] +# saturations = [[0.72783828, 0.66568458, 0.58033288], +# [0.27216172, 0.33431542, 0.41966712]] + saturations = [[0.6358, 0.5755, 0.6358], + [0.3641, 0.4245, 0.3641]] +# saturations = [[0.4, 0.5, 0.6], +# [0.6, 0.5, 0.4]] + petElProp = "bulk velocity" + input_dict = {} + input_dict['overburden'] = 'overb.npz' + + print("\nInput:") + print("porosity, pressure:", porosity, pressure) + print("phases, saturations:", phases, saturations) + print("petElProp:", petElProp) + print("input_dict:", input_dict) + + satrock = elasticproperties(input_dict) + + print("overburden:", satrock.overburden) + + satrock.calc_props(phases, saturations, pressure, + porosity) + + print("\nOutput from calc_props:") + print("Density:", satrock.getDens()) + print("Bulk modulus:", satrock.getBulkMod()) + print("Shear modulus:", satrock.getShearMod()) + print("Bulk velocity:", satrock.getBulkVel()) + print("Shear velocity:", satrock.getShearVel()) + print("Bulk impedance:", satrock.getBulkImp()) + print("Shear impedance:", satrock.getShearImp()) + + satrock.getMatchProp(petElProp) + + print("\nOutput from getMatchProp:") + print("Model output selected for data match:", + satrock.match_prop) diff --git a/simulator/rockphysics/standardrp.py b/simulator/rockphysics/standardrp.py index f847ed0..316fd04 100644 --- a/simulator/rockphysics/standardrp.py +++ b/simulator/rockphysics/standardrp.py @@ -135,6 +135,7 @@ def calc_props(self, phases, saturations, pressure, # Load "overburden" into local variable to # comply with remaining code parts overburden = self.overburden + if press_init is None: p_init = self.p_init else: @@ -209,7 +210,7 @@ def calc_props(self, phases, saturations, pressure, # Density # self.dens[i] = (porosity[i]*densf + - (1-porosity[i])*denss)*0.001 + (1-porosity[i])*denss) # # Moduli # @@ -225,10 +226,10 @@ def calc_props(self, phases, saturations, pressure, # instead of km/s) # self.bulkvel[i] = \ - 100*np.sqrt((abs(self.bulkmod[i]) + + 1000*np.sqrt((abs(self.bulkmod[i]) + 4*self.shearmod[i]/3)/(self.dens[i])) self.shearvel[i] = \ - 100*np.sqrt(self.shearmod[i] / + 1000*np.sqrt(self.shearmod[i] / (self.dens[i])) # # Impedances (m/s)*(Kg/m3)