From 87f7f6a9c471ac1d183ed27f35a6394d7d776b66 Mon Sep 17 00:00:00 2001 From: mlie Date: Thu, 6 Nov 2025 14:30:40 +0100 Subject: [PATCH] Move simulators to SimulatorWrap --- docs/gen_ref_pages.py | 6 +- simulator/calc_pem.py | 105 - simulator/eclipse.py | 1131 ----------- simulator/flow_rock.py | 2804 --------------------------- simulator/opm.py | 378 ---- simulator/rockphysics/__init__.py | 1 - simulator/rockphysics/softsandrp.py | 881 --------- simulator/rockphysics/standardrp.py | 795 -------- 8 files changed, 5 insertions(+), 6096 deletions(-) delete mode 100644 simulator/calc_pem.py delete mode 100644 simulator/eclipse.py delete mode 100644 simulator/flow_rock.py delete mode 100644 simulator/opm.py delete mode 100644 simulator/rockphysics/__init__.py delete mode 100644 simulator/rockphysics/softsandrp.py delete mode 100644 simulator/rockphysics/standardrp.py diff --git a/docs/gen_ref_pages.py b/docs/gen_ref_pages.py index fd166860..b57b056f 100644 --- a/docs/gen_ref_pages.py +++ b/docs/gen_ref_pages.py @@ -19,7 +19,11 @@ # Skip files that don't have docstrings (avoid build abortion) # More elaborate solution: https://github.com/mkdocstrings/mkdocstrings/discussions/412 - if (txt := path.read_text()) and txt.splitlines()[0][0] not in ["'", '"']: + if (txt := path.read_text()) and (not txt.splitlines()[0] or txt.splitlines()[0][0] not in ["'", '"']): + print(path,' will not be included in docs because of unknown issue') + continue + except UnicodeDecodeError: + print(path, ' will not be included in docs because of unknown issue') continue parts = tuple(path.relative_to(src).with_suffix("").parts) diff --git a/simulator/calc_pem.py b/simulator/calc_pem.py deleted file mode 100644 index 6186f420..00000000 --- a/simulator/calc_pem.py +++ /dev/null @@ -1,105 +0,0 @@ -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/eclipse.py b/simulator/eclipse.py deleted file mode 100644 index 32d83d06..00000000 --- a/simulator/eclipse.py +++ /dev/null @@ -1,1131 +0,0 @@ -"""Wrap Eclipse""" -# External imports -import numpy as np -import sys -import os -from copy import deepcopy -from mako.lookup import TemplateLookup -from mako.runtime import Context -from multiprocessing import Process -import datetime as dt -from scipy import interpolate -from subprocess import call, DEVNULL -from misc import ecl, grdecl -from shutil import rmtree, copytree # rmtree for removing folders -import time -# import rips -from glob import glob - -# Internal imports -from misc.system_tools.environ_var import EclipseRunEnvironment -from pipt.misc_tools.analysis_tools import store_ensemble_sim_information -from pipt.misc_tools.extract_tools import list_to_dict - - -class eclipse: - """ - Class for running the Schlumberger eclipse 100 black oil reservoir simulator. For more information see GeoQuest: - ECLIPSE reference manual 2009.1. Schlumberger, GeoQuest (2009). - - To run this class, eclipse must be installed and elcrun must be in the system path! - """ - - def __init__(self, input_dict=None, filename=None, options=None): - """ - The inputs are all optional, but in the same fashion as the other simulators a system must be followed. - The input_dict can be utilized as a single input. Here all nescessary info is stored. Alternatively, - if input_dict is not defined, all the other input variables must be defined. - - Parameters - ---------- - input_dict : dict, optional - Dictionary containing all information required to run the simulator. - - - parallel: number of forward simulations run in parallel - - simoptions: options for the simulations - - mpi: option to use mpi (always use > 2 cores) - - sim_path: Path to the simulator - - sim_flag: Flags sent to the simulator (see simulator documentation for all possibilities) - - sim_limit: maximum number of seconds a simulation can run before being killed - - runfile: name of the simulation input file - - reportpoint: these are the dates the simulator reports results - - reporttype: this key states that the report poins are given as dates - - datatype: the data types the simulator reports - - replace: replace failed simulations with randomly selected successful ones - - rerun: in case of failure, try to rerun the simulator the given number of times - - startdate: simulaton start date - - saveforecast: save the predicted measurements for each iteration - - filename : str, optional - Name of the .mako file utilized to generate the ECL input .DATA file. Must be in uppercase for the - ECL simulator. - - options : dict, optional - Dictionary with options for the simulator. - - Returns - ------- - initial_object : object - Initial object from the class ecl_100. - """ - - # IN - self.input_dict = input_dict - self.file = filename - self.options = options - - self.upscale = None - # If input option 1 is selected - if self.input_dict is not None: - self._extInfoInputDict() - - # Allocate for internal use - self.inv_stat = None - self.static_state = None - - # Multilevel default value - self.level = -1 - - def _extInfoInputDict(self): - """ - Extract the manditory and optional information from the input_dict dictionary. - - Parameters - ---------- - input_dict : dict - Dictionary containing all information required to run the simulator (defined in self). - - Returns - ------- - filename : str - Name of the .mako file utilized to generate the ECL input .DATA file. - """ - # Chech for mandatory keys - assert 'reporttype' in self.input_dict, 'Reporttype is missing, please specify this' - assert 'reportpoint' in self.input_dict, 'Reportpoint is missing, please specify this' - - self.true_prim = [self.input_dict['reporttype'], self.input_dict['reportpoint']] - - self.true_order = [self.input_dict['reporttype'], self.input_dict['reportpoint']] - self.all_data_types = self.input_dict['datatype'] - self.l_prim = [int(i) for i in range(len(self.true_prim[1]))] - - # In the ecl framework, all reference to the filename should be uppercase - self.file = self.input_dict['runfile'].upper() - - # Extract sim options - simoptions = self.input_dict.get('simoptions', {}) - if isinstance(simoptions, list): - simoptions = list_to_dict(simoptions) - - self.options = {} - self.options['sim_path'] = simoptions.get('sim_path', '') - self.options['sim_flag'] = simoptions.get('sim_flag', '') - self.options['mpi'] = simoptions.get('mpi', '') - self.options['mpiarray'] = simoptions.get('mpiarray', '') - self.options['parsing-strictness'] = simoptions.get('parsing-strictness', '') - - if 'sim_limit' in self.input_dict: - self.options['sim_limit'] = self.input_dict['sim_limit'] - - if 'reportdates' in self.input_dict: - self.reportdates = [ - x * 30 for x in range(1, int(self.input_dict['reportdates'][1]))] - - if 'read_sch' in self.input_dict: - # self.read_sch = self.input_dict['read_sch'][0] - load_file = np.load(self.input_dict['read_sch'], allow_pickle=True) - self.reportdates = load_file[load_file.files[0]] - - if 'startdate' in self.input_dict: - self.startDate = {} - # assume date is on form day/month/year - tmpDate = [int(elem) for elem in self.input_dict['startdate'].split('/')] - - self.startDate['day'] = tmpDate[0] - self.startDate['month'] = tmpDate[1] - self.startDate['year'] = tmpDate[2] - - if 'realizations' in self.input_dict: - self.realizations = self.input_dict['realizations'] - - if 'trunc_level' in self.input_dict: - self.trunc_level = self.input_dict['trunc_level'] - - if 'rerun' in self.input_dict: - self.rerun = int(self.input_dict['rerun']) - else: - self.rerun = 0 - - # 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 self.input_dict: - # Make sure "ANALYSISDEBUG" gives a list - if isinstance(self.input_dict['savesiminfo'], list): - self.saveinfo = self.input_dict['savesiminfo'] - else: - self.saveinfo = [self.input_dict['savesiminfo']] - - if 'upscale' in self.input_dict: - self.upscale = {} - for i in range(0, len(self.input_dict['upscale'])): - if self.input_dict['upscale'][i][0] == 'state': - # Set the parameter we upscale with regards to, must be the same as one free parameter - self.upscale['state'] = self.input_dict['upscale'][i][1] - # Set the dimension of the parameterfield to be upscaled (x and y) - self.upscale['dim'] = self.input_dict['upscale'][i][2] - if self.input_dict['upscale'][i][0] == 'maxtrunc': - # Set the ratio for the maximum truncation value - self.upscale['maxtrunc'] = self.input_dict['upscale'][i][1] - if self.input_dict['upscale'][i][0] == 'maxdiff': - # Set the ratio for the maximum truncation of differences - self.upscale['maxdiff'] = self.input_dict['upscale'][i][1] - if self.input_dict['upscale'][i][0] == 'wells': - # Set the list of well indexes, this is a list of lists where each element in the outer list gives - # a well coordinate (x and y) as elements in the inner list. - self.upscale['wells'] = [] - for j in range(1, len(self.input_dict['upscale'][i])): - self.upscale['wells'].append( - [int(elem) for elem in self.input_dict['upscale'][i][j]]) - if self.input_dict['upscale'][i][0] == 'radius': - # List of radius lengths - self.upscale['radius'] = [int(elem) - for elem in self.input_dict['upscale'][i][1]] - if self.input_dict['upscale'][i][0] == 'us_type': - self.upscale['us_type'] = self.input_dict['upscale'][i][1] - - # Check that we have as many radius elements as wells - if 'radius' in self.upscale: - if len(self.upscale['radius']) != len(self.upscale['wells']): - sys.exit('ERROR: Missmatch between number of well inputs and number of radius elements. Please check ' - 'the input file') - else: - self.upscale['radius'] = [] - self.upscale['wells'] = [] - - # The simulator should run on different levels - if 'multilevel' in self.input_dict: - # extract list of levels - self.multilevel = self.input_dict['multilevel'] - else: - # if not initiallize as list with one element - self.multilevel = [False] - - def setup_fwd_run(self, **kwargs): - """ - Setup the simulator. - - Attributes - ---------- - assimIndex : int - Gives the index-type (e.g. step,time,etc.) and the index for the - data to be assimilated - trueOrder : - Gives the index-type (e.g. step,time,etc.) and the index of the true data - """ - - self.level = -1 # default value - self.__dict__.update(kwargs) # parse kwargs input into class attributes - - if hasattr(self, 'reportdates'): - self.report = {'dates': self.reportdates} - elif 'reportmonths' in self.input_dict: # for optimization - self.report = { - 'days': [30 * i for i in range(1, int(self.input_dict['reportmonths'][1]))]} - else: - assimIndex = [i for i in range(len(self.l_prim))] - trueOrder = self.true_order - - self.pred_data = [deepcopy({}) for _ in range(len(assimIndex))] - for ind in self.l_prim: - for key in self.all_data_types: - self.pred_data[ind][key] = np.zeros((1, 1)) - - if isinstance(trueOrder[1], list): # Check if true data prim. ind. is a list - self.true_prim = [trueOrder[0], [x for x in trueOrder[1]]] - else: # Float - self.true_prim = [trueOrder[0], [trueOrder[1]]] - # self.all_data_types = list(pred_data[0].keys()) - - # Initiallise space to store the number of active cells. This is only for the upscaling option. - if 'upscale' in self.input_dict: - self.num_act = [] - - # Initialise error summary - self.error_smr = [] - - # Initiallize run time summary - self.run_time = [] - - # Check that the .mako file is in the current working directory - if not os.path.isfile('%s.mako' % self.file): - sys.exit( - 'ERROR: .mako file is not in the current working directory. This file must be defined') - - def run_fwd_sim(self, state, member_i, del_folder=True,nosim=False): - """ - Setup and run the ecl_100 forward simulator. All the parameters are defined as attributes, and the name of the - parameters are initialized in setupFwdRun. This method sets up and runs all the individual ensemble members. - This method is based on writing .DATA file for each ensemble member using the mako template, for more info - regarding mako see http://www.makotemplates.org/ - - Parameters - ---------- - state : dict - Dictionary containing the ensemble state. - - member_i : int - Index of the ensemble member. - - del_folder : bool, optional - Boolean to determine if the ensemble folder should be deleted. Default is False. - - nosim : bool, optional - Boolean to determine if the simulation should be run. Default is False. - """ - if hasattr(self, 'level'): - state['level'] = self.level - else: - state['level'] = -1 # default value - os.mkdir('En_' + str(member_i)) - folder = 'En_' + str(member_i) + os.sep - - state['member'] = member_i - # If the run is upscaled, run the upscaling procedure - if self.upscale is not None: - if hasattr(self, 'level'): # if this is a multilevel run, we must set the level - self.upscale['maxtrunc'] = self.trunc_level[self.level] - if self.upscale['maxtrunc'] > 0: # if the truncation level is 0, we do not perform upscaling - self.coarsen(folder, state) - # if the level is 0, and upscaling has been performed earlier this must be - elif hasattr(self, 'coarse'): - # removed - del self.coarse - - # start by generating the .DATA file, using the .mako template situated in ../folder - self._runMako(folder, state) - if nosim: # for hpc we only want to generate the .DATA file - return # exit the function - else: - 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 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 remove_folder(self, member): - folder = 'En_' + str(member) + os.sep - try: - rmtree(folder) # Try to delete folder - except: # If deleting fails, just rename to 'En__date' - os.rename(folder, folder + '_' + - dt.datetime.now().strftime("%m-%d-%Y_%H-%M-%S")) - - def extract_data(self, member): - # get the formated data - 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 self.pred_data[prim_ind][key] is not None: # Obs. data at assim. step - true_data_info = [self.true_prim[0], self.true_prim[1][prim_ind]] - try: - data_array = self.get_sim_results(key, true_data_info, member) - self.pred_data[prim_ind][key] = data_array - except: - print(f'Failed to extract {key} at {prim_ind} for member {member}') - pass - - def coarsen(self, folder, ensembleMember=None): - """ - This method utilized one field parameter to upscale the computational grid. A coarsening file is written to the - ensemble folder, and the eclipse calculates the upscaled permeabilities, porosities and transmissibilities - based on the new grid and the original parameter values. - - Parameters - ---------- - folder : str, optional - Path to the ecl_100 run folder. - - ensembleMember : int, optional - Index of the ensemble member to run. - - Changelog - --------- - - KF 17/9-2015 Added uniform upscaling as an option - - KF 6/01-17 - """ - - # Get the parameter values for the current ensemble member. Note that we select only one parameter - if ensembleMember is not None: - coarsenParam = self.inv_state[self.upscale['state']][:, ensembleMember] - else: - coarsenParam = self.inv_state[self.upscale['state']] - - if self.upscale['us_type'].lower() == 'haar': # the upscaling is based on a Haar wavelet - # Do not add any dead-cells - # TODO: Make a better method for determining dead or inactive cells. - - orig_wght = np.ones((self.upscale['dim'][0], self.upscale['dim'][1])) - # Get the new cell structure by a 2-D unbalanced Haar transform. - wght = self._Haar(coarsenParam.reshape( - (self.upscale['dim'][0], self.upscale['dim'][1])), orig_wght) - - elif self.upscale['us_type'].lower() == 'unif': - # This option generates a grid where the cells are upscaled uniformly in a dyadic manner. We utilize the - # same framework to define the parameters to be coarsened as in the Haar case. Hence, we only need to - # calculate the wght variable. The level of upscaling can be determined by a fraction, however, for this - # case the fraction gives the number of upscaling levels, 1 is all possible levels, 0 is no upscaling. - wght = self._unif((self.upscale['dim'][0], self.upscale['dim'][1])) - - self.write_coarse(folder, wght, coarsenParam.reshape( - (int(self.upscale['dim'][0]), int(self.upscale['dim'][1])))) - - def write_coarse(self, folder, whgt, image): - """ - This function writes the include file coarsen to the ecl run. This file tels ECL to coarsen the grid. - - Parameters - ---------- - folder : str - Path to the ecl_100 run. - - whgt : float - Weight of the transformed cells. - - image : array-like - Original image. - - Changelog - --------- - - KF 17/9-2015 - """ - - well = self.upscale['wells'] - radius = self.upscale['radius'] - - well_cells = self._nodeIndex(image.shape[1], image.shape[0], well, radius) - coarse = np.array([[True] * image.shape[1]] * image.shape[0]) - tmp = coarse.flatten() - tmp[well_cells] = False - coarse = tmp.reshape(image.shape) - - # f = open(folder + 'coarsen.dat', 'w+') - # f.write('COARSEN \n') - ecl_coarse = [] - - for level in range(len(whgt) - 1, -1, -1): - weight_at_level = whgt[level] - x_dim = weight_at_level.shape[0] - y_dim = weight_at_level.shape[1] - - merged_at_level = weight_at_level[int(x_dim / 2):, :int(y_dim / 2)] - - level_dim = 2 ** (level + 1) - - for i in range(0, merged_at_level.shape[0]): - for j in range(0, merged_at_level.shape[1]): - if merged_at_level[i, j] == 1: - # Do this in two stages to avoid errors with the dimensions - # - # only merging square cells - if (i + 1) * level_dim <= image.shape[0] and (j + 1) * level_dim <= image.shape[1]: - if coarse[i * level_dim:(i + 1) * level_dim, j * level_dim:(j + 1) * level_dim].all(): - coarse[i * level_dim:(i + 1) * level_dim, - j * level_dim:(j + 1) * level_dim] = False - ecl_coarse.append([j * level_dim + 1, (j + 1) * level_dim, - i * level_dim + 1, (i + 1) * level_dim, 1, 1, 1, 1, 1]) - # f.write('%i %i %i %i 1 1 1 1 1 / \n' % (j * level_dim + 1, (j + 1) * level_dim, - # i * level_dim + 1, (i + 1) * level_dim)) - # cells at first edge, non-square - if (i + 1)*level_dim > image.shape[0] and (j + 1) * level_dim <= image.shape[1]: - if coarse[i * level_dim::, j * level_dim:(j + 1) * level_dim].all(): - coarse[i * level_dim::, j * - level_dim:(j + 1) * level_dim] = False - ecl_coarse.append([j * level_dim + 1, (j + 1) * level_dim, - i * level_dim + 1, image.shape[0], 1, 1, 1, 1, 1]) - # f.write('%i %i %i %i 1 1 1 1 1 / \n' % (j * level_dim + 1, (j + 1) * level_dim, - # i * level_dim + 1, image.shape[0])) - # cells at second edge, non-square - if (j + 1)*level_dim > image.shape[1] and (i + 1) * level_dim <= image.shape[0]: - if coarse[i * level_dim:(i + 1) * level_dim, j * level_dim::].all(): - coarse[i * level_dim:(i + 1) * level_dim, - j * level_dim::] = False - ecl_coarse.append([j * level_dim + 1, image.shape[1], - i * level_dim + 1, (i + 1) * level_dim, 1, 1, 1, 1, 1]) - # f.write('%i %i %i %i 1 1 1 1 1 / \n' % (j * level_dim + 1, image.shape[1], - # i * level_dim + 1, (i + 1) * level_dim)) - # cells at intersection between first and second edge, non-square - if (i + 1) * level_dim > image.shape[0] and (j + 1) * level_dim > image.shape[1]: - if coarse[i * level_dim::, j * level_dim::].all(): - coarse[i * level_dim::, j * level_dim::] = False - ecl_coarse.append([j * level_dim + 1, image.shape[1], - i * level_dim + 1, image.shape[0], 1, 1, 1, 1, 1]) - # f.write('%i %i %i %i 1 1 1 1 1 / \n' % (j * level_dim + 1, image.shape[1], - # i * level_dim + 1, image.shape[0])) - # f.write('/') - # f.close() - self.coarse = ecl_coarse - - def _nodeIndex(self, x_dir, y_dir, cells, listRadius): - # Find the node index for the cells, and the cells in a radius around. - index = [] - for k in range(0, len(cells)): - cell_x = cells[k][0] - 1 # Remove one to make python equivalent to ecl format - cell_y = cells[k][1] - 1 # Remove one to make python equivalent to ecl format - radius = listRadius[k] - # tmp_index = [cell_x + cell_y*x_dir] - tmp_index = [] - for i in range(-radius, radius + 1): - for j in range(-radius, radius + 1): - x = cell_x + i - y = cell_y + j - if (0 <= x < x_dir and 0 <= y < y_dir) and (np.sqrt(i ** 2 + j ** 2) <= radius) \ - and (x + y * x_dir < y_dir * x_dir): - tmp_index.append(x + y * x_dir) - index.extend(tmp_index) - - return index - - def _Haar(self, image, weights): - """ - The 2D unconditional Haar transform. - - Parameters - ---------- - image : array-like - The original image where the Haar-transform is performed. This could be the permeability or the porosity field. - - orig_weights : array-like - The original weights of the image. If some cells are inactive or dead, their weight can be set to zero. - - Returns - ------- - tot_weights : array-like - A vector/matrix of the same size as the image, with the weights of the new cells. This is utilized for writing - the coarsening file but is not sufficient to recreate the image. - """ - # Todo: make an option which stores the transformed field. - - # Define the two truncation levels as a fraction of the larges value in the original image. The fraction is defined - # through the inputs. - max_value = self.upscale['maxtrunc'] * np.max(image) - max_diff = self.upscale['maxdiff'] * (np.max(image) - np.min(image)) - - tot_transf = [] - tot_weight = [] - tot_alpha = [] - count = 0 - allow_merge = np.ones(image.shape) - while image.shape[0] > 1: - level_alpha = [] - if isinstance(image.shape, tuple): - # the image has two axes. - # Perform the transformation for all the coulmns, and then for the transpose of the results - - for axis in range(0, 2): - # Initialise the matrix for storing the details and smoothing - level_diff_columns = np.empty( - (int(np.ceil(len(image[:, 0]) / 2)), len(image[0, :]))) - # level_diff_rows = np.empty((len(image[:, 0]), int(np.ceil(len(image[0, :])/2)))) - level_smooth_columns = np.empty( - (int(np.ceil(len(image[:, 0]) / 2)), len(image[0, :]))) - level_weight_columns = np.empty( - (int(np.ceil(len(image[:, 0]) / 2)), len(image[0, :]))) - level_alpha_columns = np.empty( - (int(np.ceil(len(image[:, 0]) / 2)), len(image[0, :]))) - - for i in range(0, image.shape[1]): - if len(image[:, i]) % 2 == 0: - diff = image[1::2, i] - image[::2, i] - new_weights = weights[::2, i] + weights[1::2, i] - alpha = (weights[1::2, i] / new_weights) - smooth = image[::2, i] + alpha * diff - - else: - tmp_img = image[:-1, i] - diff = tmp_img[1::2] - tmp_img[::2] - - tmp_weight = weights[:-1, i] - new_weights = tmp_weight[::2] + tmp_weight[1::2] - new_weights = np.append(new_weights, weights[-1, i]) - - tmp_img = image[:-1, i] - alpha = (weights[1::2, i] / new_weights[:-1]) - smooth = tmp_img[::2] + alpha * diff - smooth = np.append(smooth, image[-1, i]) - alpha = np.append(alpha, 0) - diff = np.append(diff, 0) - - level_diff_columns[:, i] = diff - level_weight_columns[:, i] = new_weights - level_smooth_columns[:, i] = smooth - level_alpha_columns[:, i] = alpha - - # image = level_smooth_columns.T - image = np.vstack((level_smooth_columns, level_diff_columns)).T - weights = np.vstack((level_weight_columns, level_weight_columns)).T - level_alpha.append(level_alpha_columns) - - image, weights, allow_merge, end_ctrl = self._haarTrunc(image, weights, max_value, max_diff, - allow_merge) - - tot_transf.append(image) - tot_weight.append(weights) - tot_alpha.append(level_alpha) - - if end_ctrl: - break - image = image[:image.shape[0] / 2, :image.shape[1] / 2] - weights = weights[:weights.shape[0] / 2, :weights.shape[1] / 2] - - else: - if len(image) % 2 == 0: - diff = image[1::2] - image[::2] - new_weights = weights[::2] + weights[1::2] - smooth = image[::2] + (weights[1::2] / new_weights) * diff - - else: - tmp_img = image[:-1] - diff = tmp_img[1::2] - tmp_img[::2] - - tmp_weight = weights[:-1] - new_weights = tmp_weight[::2] + tmp_weight[1::2] - new_weights = np.append(new_weights, weights[-1]) - - tmp_img = image[:-1] - smooth = tmp_img[::2] + (weights[1::2] / new_weights[:-1]) * diff - smooth = np.append(smooth, image[-1]) - - tot_transf.append(smooth) - tot_weight.append(new_weights) - - # weights = new_weights - # image = smooth - count += 1 - - return tot_weight - - def _unif(self, dim): - # calculate the uniform upscaling steps - - min_dim = min(dim) # find the lowest dimension - max_levels = 0 - while min_dim/2 >= 1: - min_dim = min_dim/2 - max_levels += 1 - - # conservative choice - trunc_level = int(np.floor(max_levels*self.upscale['maxtrunc'])) - - # all the initial cells are upscaled - wght = [np.ones(tuple([int(elem) for elem in dim]))] - - for i in range(trunc_level): - new_dim = [int(np.ceil(elem/2)) - for elem in wght[i].shape] # always take a larger dimension - wght.append(np.ones(tuple(new_dim))) - - return wght - - def _haarTrunc(self, image, weights, max_val, max_diff, merge): - """ - Function for truncating the wavelets. Based on the max_val and max_diff values this function set the detail - coaffiecient to zero if the value of this coefficient is below max_diff, and the value of the smooth coefficient - is below max_val. - - Parameters - ---------- - image : array-like - The transformed image. - - weights : array-like - The weights of the transformed image. - - max_val : float - Smooth values above this value are not allowed. - - max_diff : float - Detail coefficients higher than this value are not truncated. - - merge : array-like - Matrix/vector of booleans defining whether a merge is allowed. - - Returns - ------- - image : array-like - Transformed image with truncated coefficients. - - weights : array-like - Updated weights. - - allow_merge : array-like - Booleans keeping track of allowed merges. - - end_ctrl : bool - Boolean to control whether further upscaling is possible. - """ - # the image will always contain even elements in x and y dir - x_dir = int(image.shape[0] / 2) - y_dir = int(image.shape[1] / 2) - - smooth = image[:x_dir, :y_dir] - smooth_diff = image[x_dir:, :y_dir] - diff_smooth = image[:x_dir, y_dir:] - diff_diff = image[x_dir:, y_dir:] - - weights_diff = weights[x_dir:, :y_dir] - - allow_merge = np.array( - [[False] * y_dir] * x_dir) - - # ensure that last row does not be left behind during non-dyadic upscaling - if image.shape[0] > merge.shape[0]: - merge = np.insert(merge, -1, merge[-1, :], axis=0) - if image.shape[1] > merge.shape[1]: - merge = np.insert(merge, -1, merge[:, -1], axis=1) - - cand1 = zip(*np.where(smooth < max_val)) - - end_ctrl = True - - for elem in cand1: - # If the method wants to merge cells outside the box the if sentence gives an error, make an exception for - # this hence do not merge these cells. - try: - if abs(smooth_diff[elem]) <= max_diff and abs(diff_smooth[elem]) <= max_diff and \ - merge[elem[0] * 2, elem[1] * 2] and merge[elem[0] * 2 + 1, elem[1] * 2] and \ - merge[elem[0] * 2 + 1, elem[1] * 2 + 1] and merge[elem[0] * 2, elem[1] * 2 + 1]: - diff_smooth[elem] = 0 - smooth_diff[elem] = 0 - diff_diff[elem] = 0 - weights_diff[elem] = True - allow_merge[elem] = True - end_ctrl = False - except: - pass - # for elem in cand2: - # if smooth[elem] <= max_val: - # diff_smooth[elem] = 0 - # smooth_diff[elem] = 0 - # diff_diff[elem] = 0 - - image[x_dir:, :y_dir] = smooth_diff - image[:x_dir, y_dir:] = diff_smooth - image[x_dir:, y_dir:] = diff_diff - - weights[x_dir:, :y_dir] = weights_diff - - return image, weights, allow_merge, end_ctrl - - def _runMako(self, folder, state): - """ - Read the mako template (.mako) file from ../folder, and render the correct data file (.DATA) in folder. - - Parameters - ---------- - Folder : str, optional - Folder for the ecl_100 run. - - ensembleMember : int, optional - Index of the ensemble member to run. - - Changelog - --------- - - KF 14/9-2015 - """ - - # Check and add report time - if hasattr(self, 'report'): - for key in self.report: - state[key] = self.report[key] - - if 'mako_kwargs' in self.input_dict: - mako_kwargs = dict(self.input_dict['mako_kwargs']) - state.update(mako_kwargs) - - # Convert drilling order (float) to drilling queue (integer) - drilling order optimization - # if "drillingorder" in en_dict: - # dorder = en_dict['drillingorder'] - # en_dict['drillingqueue'] = np.argsort(dorder)[::-1][:len(dorder)] - - # Add startdate - if hasattr(self, 'startDate'): - state['startdate'] = dt.datetime( - self.startDate['year'], self.startDate['month'], self.startDate['day']) - - # Add the coarsing values - if hasattr(self, 'coarse'): - state['coarse'] = self.coarse - - # Look for the mako file - lkup = TemplateLookup(directories=os.getcwd(), - input_encoding='utf-8') - - # Get template - # If we need the E300 run, define a E300 data file with _E300 added to the end. - if hasattr(self, 'E300'): - if self.E300: - tmpl = lkup.get_template('%s.mako' % (self.file + '_E300')) - else: - tmpl = lkup.get_template('%s.mako' % self.file) - else: - tmpl = lkup.get_template('%s.mako' % self.file) - - # use a context and render onto a file - with open('{0}{1}'.format(folder + self.file, '.DATA'), 'w') as f: - ctx = Context(f, **state) - tmpl.render_context(ctx) - - def get_sim_results(self, whichResponse, ext_data_info=None, member=None): - """ - Read the output from simulator and store as an array. Optionally, if the DA method is based on an ensemble - method, the output must be read inside a folder. - - Parameters - ---------- - whichResponse : str - Which of the responses is to be outputted (e.g., WBHP PRO-1, WOPR, PRESS, OILSAT, etc). - - ext_data_info : tuple, optional - Tuple containing the assimilation step information, including the place of assimilation (e.g., which TIME) and the - index of this assimilation place. - - member : int, optional - Ensemble member that is finished. - - Returns - ------- - yFlow : array-like - Array containing the response from ECL 100. The response type is chosen by the user in options['data_type']. - - Notes - ----- - - Modified the ecl package to allow reading the summary data directly, hence, we get cell, summary, and field data - from the ecl package. - - KF 29/10-2015 - - Modified the ecl package to read RFT files. To obtain, e.g,. RFT pressures form well 'PRO-1" whichResponse - must be rft_PRESSURE PRO-1 - """ - # Check that we have no trailing spaces - whichResponse = whichResponse.strip() - - # if ensemble DA method - if member is not None: - # Get results - if hasattr(self, 'ecl_case'): - # En_XX/YYYY.DATA is the folder setup - rt_mem = int(self.ecl_case.root.split('/')[0].split('_')[1]) - if rt_mem != member: # wrong case - self.ecl_case = ecl.EclipseCase('En_' + str(member) + os.sep + self.file + '.DATA') - else: - self.ecl_case = ecl.EclipseCase('En_' + str(member) + os.sep + self.file + '.DATA') - if ext_data_info[0] == 'days': - time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ - dt.timedelta(days=ext_data_info[1]) - dates = self.ecl_case.by_date - if time not in dates and 'date_slack' in self.input_dict: - slack = int(self.input_dict['date_slack']) - if slack > 0: - v = [el for el in dates if np.abs( - (el-time).total_seconds()) < slack] - if len(v) > 0: - time = v[0] - else: - time = ext_data_info[1] - - # Check if the data is a field or well data, by checking if the well is defined - if len(whichResponse.split(' ')) == 2: - # if rft, search for rft_ - if 'rft_' in whichResponse: - # to speed up when performing the prediction step - if hasattr(self, 'rft_case'): - rt_mem = int(self.rft_case.root.split('/')[0].split('_')[1]) - if rt_mem != member: - self.rft_case = ecl.EclipseRFT('En_' + str(member) + os.sep + self.file) - else: - self.rft_case = ecl.EclipseRFT( - 'En_' + str(member) + os.sep + self.file) - # Get the data. Due to formating we can slice the property. - rft_prop = self.rft_case.rft_data(well=whichResponse.split( - ' ')[1], prop=whichResponse.split(' ')[0][4:]) - # rft_data are collected for open connections. This may vary throughout the simulation, hence we - # must also collect the depth for the rft_data to check if all data is present - rft_depth = self.rft_case.rft_data( - well=whichResponse.split(' ')[1], prop='DEPTH') - # to check this we import the referance depth if this is available. If not we assume that the data - # is ok. - try: - ref_depth_f = np.load(whichResponse.split( - ' ')[1].upper() + '_rft_ref_depth.npz') - ref_depth = ref_depth_f[ref_depth_f.files[0]] - yFlow = np.array([]) - interp = interpolate.interp1d(rft_depth, rft_prop, kind='linear', bounds_error=False, - fill_value=(rft_prop[0], rft_prop[-1])) - for d in ref_depth: - yFlow = np.append(yFlow, interp(d)) - except: - yFlow = rft_prop - else: - # If well, read the rsm file - if ext_data_info is not None: # Get the data at a specific well and time - yFlow = self.ecl_case.summary_data(whichResponse, time) - elif len(whichResponse.split(' ')) == 1: # field data - if whichResponse.upper() in ['FOPT', 'FWPT', 'FGPT', 'FWIT', 'FGIT']: - if ext_data_info is not None: - yFlow = self.ecl_case.summary_data(whichResponse, time) - elif whichResponse.upper() in ['PERMX', 'PERMY', 'PERMZ', 'PORO', 'NTG', 'SATNUM', - 'MULTNUM', 'OPERNUM']: - yFlow = np.array([self.ecl_case.cell_data(whichResponse).flatten()[time]]) # assume that time is the index - else: - yFlow = self.ecl_case.cell_data(whichResponse, time).flatten() - if yFlow is None: - yFlow = self.ecl_case.cell_data(whichResponse).flatten() - - # store the run time. NB: elapsed must be defined in .DATA file for this to work - if 'save_elapsed' in self.input_dict and len(self.run_time) <= member: - self.run_time.extend(self.ecl_case.summary_data('ELAPSED', time)) - - # If we have performed coarsening, we store the number of active grid-cells - if self.upscale is not None: - # Get this number from INIT file - with ecl.EclipseFile('En_' + str(member) + os.sep + self.file, 'INIT') as case: - intHead = case.get('INTEHEAD') - # The active cell is element 12 of this vector, index 11 in python indexing... - active_cells = intHead[11] - if len(self.num_act) <= member: - self.num_act.extend([active_cells]) - - else: - case = ecl.EclipseCase(self.file + '.DATA') - if ext_data_info[0] == 'days': - time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ - dt.timedelta(days=ext_data_info[1]) - else: - time = ext_data_info[1] - - # Check if the data is a field or well data, by checking if the well is defined - if len(whichResponse.split(' ')) == 2: - # if rft, search for rft_ - if 'rft_' in whichResponse: - rft_case = ecl.EclipseRFT(self.file) - # Get the data. Due to formating we can slice the property. - rft_prop = rft_case.rft_data(well=whichResponse.split( - ' ')[1], prop=whichResponse.split(' ')[0][4:]) - # rft_data are collected for open connections. This may vary throughout the simulation, hence we - # must also collect the depth for the rft_data to check if all data is present - rft_depth = rft_case.rft_data( - well=whichResponse.split(' ')[1], prop='DEPTH') - try: - ref_depth_f = np.load(whichResponse.split( - ' ')[1].upper() + '_rft_ref_depth.npz') - ref_depth = ref_depth_f[ref_depth_f.files[0]] - yFlow = np.array([]) - interp = interpolate.interp1d(rft_depth, rft_prop, kind='linear', bounds_error=False, - fill_value=(rft_prop[0], rft_prop[-1])) - for d in ref_depth: - yFlow = np.append(yFlow, interp(d)) - except: - yFlow = rft_prop - else: - # If well, read the rsm file - if ext_data_info is not None: # Get the data at a specific well and time - yFlow = case.summary_data(whichResponse, time) - elif len(whichResponse.split(' ')) == 1: - if whichResponse in ['FOPT', 'FWPT', 'FGPT', 'FWIT', 'FGIT']: - if ext_data_info is not None: - yFlow = case.summary_data(whichResponse, time) - else: - yFlow = case.cell_data(whichResponse, time).flatten() - if yFlow is None: - yFlow = case.cell_data(whichResponse).flatten() - - # If we have performed coarsening, we store the number of active grid-cells - if self.upscale is not None: - # Get this number from INIT file - with ecl.EclipseFile('En_' + str(member) + os.sep + self.file,'INIT') as case: - intHead = case.get('INTEHEAD') - # The active cell is element 12 of this vector, index 11 in python indexing... - active_cells = intHead[11] - if len(self.num_act) <= member: - self.num_act.extend([active_cells]) - - return yFlow - - def store_fwd_debug(self, assimstep): - if 'fwddebug' in self.keys_fwd: - # Init dict. of variables to save - save_dict = {} - - # Make sure "ANALYSISDEBUG" gives a list - if isinstance(self.keys_fwd['fwddebug'], list): - debug = self.keys_fwd['fwddebug'] - else: - debug = [self.keys_fwd['fwddebug']] - - # Loop over variables to store in save list - for var in debug: - # Save with key equal variable name and the actual variable - if isinstance(eval('self.' + var), dict): - # save directly - np.savez('fwd_debug_%s_%i' % (var, assimstep), **eval('self.' + var)) - else: - save_dict[var] = eval('self.' + var) - - np.savez('fwd_debug_%i' % (assimstep), **save_dict) - - def write_to_grid(self, value, propname, path, dim, t_ind=None): - if t_ind == None: - trans_dict = {} - - def _lookup(kw): - return trans_dict[kw] if kw in trans_dict else kw - - # Write a quantity to the grid as a grdecl file - with open(path + propname + '.grdecl', 'wb') as fileobj: - grdecl._write_kw(fileobj, propname, value, _lookup, dim) - else: - pass - # some errors with rips - # p = Process(target=_write_to_resinsight, args=(list(value[~value.mask]),propname, t_ind)) - # Find an open resinsight case - # p.start() - # time.sleep(1) - # p.terminate() - -# def _write_to_resinsight(value, name,t_ind): -# resinsight = rips.Instance.find() -# case = resinsight.project.case(case_id=0) -# case.set_active_cell_property(value, 'GENERATED', name,t_ind) - - -class ecl_100(eclipse): - ''' - ecl_100 class - ''' - - def call_sim(self, path=None, wait_for_proc=False): - """ - Method for calling the ecl_100 simulator. - - Parameters - ---------- - path : str - Alternative folder for the ecl_100.data file. - - wait_for_proc : bool, optional - Logical variable to wait for the simulator to finish. Default is False. - - Returns - ------- - .RSM : str - Run summary file in the standard ECL format. Well data are collected from this file. - - .RST : str - Restart file in the standard ECL format. Pressure and saturation data are collected from this file. - - .PRT : str - Info file to be used for checking errors in the run. - - - Changelog - --------- - - KF 14/9-2015 - """ - - # Filename - if path is not None: - filename = path + self.file - else: - filename = self.file - - # Run the simulator: - success = True - try: - with EclipseRunEnvironment(filename): - com = ['eclrun', '--nocleanup', 'eclipse', filename + '.DATA'] - if 'sim_limit' in self.options: - call(com, stdout=DEVNULL, timeout=self.options['sim_limit']) - else: - call(com, stdout=DEVNULL) - raise ValueError - except: - print('\nError in the eclipse run.') # add rerun? - if not os.path.exists('Crashdump'): - copytree(path, 'Crashdump') - success = False - - return success - - -class ecl_300(eclipse): - ''' - eclipse 300 class - ''' - - def call_sim(self, path=None, wait_for_proc=False): - """ - Method for calling the ecl_300 simulator. - - Parameters - ---------- - path : str - Alternative folder for the ecl_100.data file. - - wait_for_proc : bool, optional - Logical variable to wait for the simulator to finish. Default is False. - - !!! note - For now, this option is only utilized in a single localization option. - - Returns - ------- - RSM : str - Run summary file in the standard ECL format. Well data are collected from this file. - - RST : str - Restart file in the standard ECL format. Pressure and saturation data are collected from this file. - - PRT : str - Info file to be used for checking errors in the run. - - Changelog - --------- - - KF 8/10-2015 - - - """ - # Filename - if path is not None: - filename = path + self.file - else: - filename = self.file - - # Run the simulator: - with EclipseRunEnvironment(filename): - call(['eclrun', '--nocleanup', 'e300', filename + '.DATA'], stdout=DEVNULL) diff --git a/simulator/flow_rock.py b/simulator/flow_rock.py deleted file mode 100644 index d2b6857b..00000000 --- a/simulator/flow_rock.py +++ /dev/null @@ -1,2804 +0,0 @@ -"""Descriptive description.""" -from selectors import SelectSelector - -from simulator.opm import flow -from importlib import import_module -import datetime as dt -import numpy as np -import os -import pandas as pd -from misc import ecl, grdecl -import shutil -import glob -from subprocess import Popen, PIPE -import mat73 -from copy import deepcopy -from sklearn.cluster import KMeans -from sklearn.preprocessing import StandardScaler -from scipy.optimize import fsolve -from scipy.special import jv # Bessel function of the first kind -from scipy.integrate import quad -from scipy.special import j0 -from mako.lookup import TemplateLookup -from mako.runtime import Context -#import cProfile -#import pstats - -# from pylops import avo -from pylops.utils.wavelets import ricker -from pylops.signalprocessing import Convolve1D -import sys -#from PyGRDECL.GRDECL_Parser import GRDECL_Parser # https://github.com/BinWang0213/PyGRDECL/tree/master -from scipy.interpolate import interp1d -from scipy.interpolate import griddata -from pipt.misc_tools.analysis_tools import store_ensemble_sim_information -from geostat.decomp import Cholesky -#from simulator.eclipse import ecl_100 -from CoolProp.CoolProp import PropsSI # http://coolprop.org/#high-level-interface-example - -class mixIn_multi_data(): - - def find_cell_centre(self, grid): - # Find indices where the boolean array is True - indices = np.where(grid['ACTNUM']) - - coord = grid['COORD'] - zcorn = grid['ZCORN'] - - c, a, b = indices - # Calculate xt, yt, zt - xb = 0.25 * (coord[a, b, 0, 0] + coord[a, b + 1, 0, 0] + coord[a + 1, b, 0, 0] + coord[a + 1, b + 1, 0, 0]) - yb = 0.25 * (coord[a, b, 0, 1] + coord[a, b + 1, 0, 1] + coord[a + 1, b, 0, 1] + coord[a + 1, b + 1, 0, 1]) - zb = 0.25 * (coord[a, b, 0, 2] + coord[a, b + 1, 0, 2] + coord[a + 1, b, 0, 2] + coord[a + 1, b + 1, 0, 2]) - - xt = 0.25 * (coord[a, b, 1, 0] + coord[a, b + 1, 1, 0] + coord[a + 1, b, 1, 0] + coord[a + 1, b + 1, 1, 0]) - yt = 0.25 * (coord[a, b, 1, 1] + coord[a, b + 1, 1, 1] + coord[a + 1, b, 1, 1] + coord[a + 1, b + 1, 1, 1]) - zt = 0.25 * (coord[a, b, 1, 2] + coord[a, b + 1, 1, 2] + coord[a + 1, b, 1, 2] + coord[a + 1, b + 1, 1, 2]) - - # Calculate z, x, and y positions - z = (zcorn[c, 0, a, 0, b, 0] + zcorn[c, 0, a, 1, b, 0] + zcorn[c, 0, a, 0, b, 1] + zcorn[c, 0, a, 1, b, 1] + - zcorn[c, 1, a, 0, b, 0] + zcorn[c, 1, a, 1, b, 0] + zcorn[c, 1, a, 0, b, 1] + zcorn[c, 1, a, 1, b, 1]) / 8 - - x = xb + (xt - xb) * (z - zb) / (zt - zb) - y = yb + (yt - yb) * (z - zb) / (zt - zb) - - cell_centre = [x, y, z] - return cell_centre - - - def get_seabed_depths(self, file_path): - # Read the data while skipping the header comments - # We'll assume the header data ends before the numerical data - # The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead - water_depths = pd.read_csv(file_path, comment='#', sep=r'\s+', - header=None) # delim_whitespace=True, header=None) - - # Give meaningful column names: - water_depths.columns = ['x', 'y', 'z', 'column', 'row'] - - return water_depths - - def measurement_locations(self, grid, water_depth, pad=1500, dxy=3000, well_coord = None, dxy_fine = 1500, r0 = 5000): - - # Determine the size of the measurement area as defined by the field extent - cell_centre = self.find_cell_centre(grid) - x_min = np.min(cell_centre[0]) - x_max = np.max(cell_centre[0]) - y_min = np.min(cell_centre[1]) - y_max = np.max(cell_centre[1]) - - x_min -= pad - x_max += pad - y_min -= pad - y_max += pad - - x_span = x_max - x_min - y_span = y_max - y_min - - nx = int(np.ceil(x_span / dxy)) - ny = int(np.ceil(y_span / dxy)) - - x_vec = np.linspace(x_min, x_max, nx) - y_vec = np.linspace(y_min, y_max, ny) - x, y = np.meshgrid(x_vec, y_vec) - - # allow for finer measurement grid around injection well - if well_coord is not None: - # choose center point and radius for area with finer measurement grid - # y = 64, x = 34 # alpha well position Smeaheia - xc = cell_centre[well_coord[0]] - yc = cell_centre[well_coord[1]] - - # Fine grid covering bounding box of the circle (clamped to domain) - fx_min = max(x_min, xc - r0) - fx_max = min(x_max, xc + r0) - fy_min = max(y_min, yc - r0) - fy_max = min(y_max, yc + r0) - - pts_coarse = np.column_stack((x.ravel(), y.ravel())) - - if fx_max > fx_min and fy_max > fy_min: - nfx = int(np.ceil((fx_max - fx_min) / dxy_fine)) + 1 - nfy = int(np.ceil((fy_max - fy_min) / dxy_fine)) + 1 - x_fine = np.linspace(fx_min, fx_max, nfx) - y_fine = np.linspace(fy_min, fy_max, nfy) - xf, yf = np.meshgrid(x_fine, y_fine) - pts_fine = np.column_stack((xf.ravel(), yf.ravel())) - - # Keep only fine points inside the circle - d2 = (pts_fine[:, 0] - xc) ** 2 + (pts_fine[:, 1] - yc) ** 2 - mask_inside = d2 <= r0 ** 2 - pts_fine_inside = pts_fine[mask_inside] - - # remove the coarse points inside the circle - d2 = (pts_coarse[:, 0] - xc) ** 2 + (pts_coarse[:, 1] - yc) ** 2 - mask_inside = d2 <= r0 ** 2 - pts_coarse = pts_coarse[~mask_inside] - - # Combine and remove duplicates by rounding to a tolerance or using a structured array - # Use tolerance based on the smaller spacing - tol = min(dxy, dxy_fine) * 1e-3 - all_pts = np.vstack((pts_coarse, pts_fine_inside)) - - # Round coordinates to avoid floating point duplicates then use np.unique - # Determine digits to round so differences smaller than tol collapse - digits = max(0, int(-np.floor(np.log10(tol)))) - all_pts_rounded = np.round(all_pts, digits) - uniq_pts = np.unique(all_pts_rounded, axis=0) - x = uniq_pts[:, 0] - y = uniq_pts[:, 1] - - - pos = {'x': x.flatten(), 'y': y.flatten()} - - # Seabed map or water depth scalar depending on input - if isinstance(water_depth, float): - pos['z'] = np.ones_like(pos['x']) * water_depth - else: - pos['z'] = griddata((water_depth['x'], water_depth['y']), - np.abs(water_depth['z']), (pos['x'], pos['y']), - method='nearest') # z is positive downwards - return pos - - -class flow_rock(flow): - """ - Couple the OPM-flow simulator with a rock-physics simulator such that both reservoir quantities and petro-elastic - quantities can be calculated. Inherit the flow class, and use super to call similar functions. - """ - - def __init__(self, input_dict=None, filename=None, options=None): - super().__init__(input_dict) - self._getpeminfo(input_dict) - - self.date_slack = None - if 'date_slack' in input_dict: - self.date_slack = int(input_dict['date_slack']) - - # If we want to extract, or evaluate, something uniquely from the ensemble specific run we can - # run a user defined code to do this. - self.saveinfo = None - if 'savesiminfo' in input_dict: - # Make sure "ANALYSISDEBUG" gives a list - if isinstance(input_dict['savesiminfo'], list): - self.saveinfo = input_dict['savesiminfo'] - else: - self.saveinfo = [input_dict['savesiminfo']] - - self.scale = [] - - # Store dynamic variables in case they are provided in the state - self.state = None - self.no_flow = False - - def _getpeminfo(self, input_dict): - """ - Get, and return, flow and PEM modules - """ - if 'pem' in input_dict: - self.pem_input = {} - for elem in input_dict['pem']: - if elem[0] == 'model': # Set the petro-elastic model - self.pem_input['model'] = elem[1] - if elem[0] == 'depth': # provide the npz of depth values - self.pem_input['depth'] = elem[1] - if elem[0] == 'actnum': # the npz of actnum values - self.pem_input['actnum'] = elem[1] - if elem[0] == 'baseline': # the time for the baseline 4D measurement - self.pem_input['baseline'] = elem[1] - if elem[0] == 'vintage': - self.pem_input['vintage'] = elem[1] - if not type(self.pem_input['vintage']) == list: - self.pem_input['vintage'] = [elem[1]] - if elem[0] == 'ntg': - self.pem_input['ntg'] = elem[1] - if elem[0] == 'press_conv': - self.pem_input['press_conv'] = elem[1] - if elem[0] == 'compaction': - self.pem_input['compaction'] = True - if elem[0] == 'overburden': # the npz of overburden values - self.pem_input['overburden'] = elem[1] - if elem[0] == 'percentile': # use for scaling - self.pem_input['percentile'] = elem[1] - if elem[0] == 'phases': # get the fluid phases - self.pem_input['phases'] = elem[1] - if elem[0] == 'grid': # get the model grid - self.pem_input['grid'] = elem[1] - if elem[0] == 'param_file': # get model parameters required for pem - self.pem_input['param_file'] = elem[1] - - - pem = getattr(import_module('simulator.rockphysics.' + - self.pem_input['model'].split()[0]), self.pem_input['model'].split()[1]) - - self.pem = pem(self.pem_input) - - else: - self.pem = None - - def _get_pem_input(self, type, time=None): - if self.no_flow: # get variable from state - if any(type.lower() in key for key in self.state.keys()) and time > 0: - data = self.state[type.lower()+'_'+str(time)] - mask = np.zeros(data.shape, dtype=bool) - return np.ma.array(data=data, dtype=data.dtype, - mask=mask) - else: # read parameter from file - param_file = self.pem_input['param_file'] - npzfile = np.load(param_file) - parameter = npzfile[type] - npzfile.close() - data = parameter[:,self.ensemble_member] - mask = np.zeros(data.shape, dtype=bool) - return np.ma.array(data=data, dtype=data.dtype, - mask=mask) - else: # get variable of parameter from flow simulation - return self.ecl_case.cell_data(type,time) - - def calc_pem(self, time, time_index=None): - - if self.no_flow: - time_input = time_index - else: - time_input = time - - # fluid phases written given as input - phases = str.upper(self.pem_input['phases']) - phases = phases.split() - - pem_input = {} - tmp_dyn_var = {} - # get active porosity - tmp = self._get_pem_input('PORO') # self.ecl_case.cell_data('PORO') - if 'compaction' in self.pem_input: - multfactor = self._get_pem_input('PORV_RC', time_input) - pem_input['PORO'] = np.array(multfactor[~tmp.mask] * tmp[~tmp.mask], dtype=float) - else: - pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float) - - # get active NTG if needed - if 'ntg' in self.pem_input: - if self.pem_input['ntg'] == 'no': - pem_input['NTG'] = None - else: - tmp = self._get_pem_input('NTG') - pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) - else: - tmp = self._get_pem_input('NTG') - pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float) - - if 'RS' in self.pem_input: #ecl_case.cell_data: # to be more robust! - tmp = self._get_pem_input('RS', time_input) - pem_input['RS'] = np.array(tmp[~tmp.mask], dtype=float) - else: - pem_input['RS'] = None - print('RS is not a variable in the ecl_case') - - # extract pressure - tmp = self._get_pem_input('PRESSURE', time_input) - pem_input['PRESSURE'] = np.array(tmp[~tmp.mask], dtype=float) - - # convert pressure from Bar to MPa - if 'press_conv' in self.pem_input and time_input == time: - pem_input['PRESSURE'] = pem_input['PRESSURE'] * self.pem_input['press_conv'] - - if hasattr(self.pem, 'p_init'): - P_init = self.pem.p_init * np.ones(tmp.shape)[~tmp.mask] - else: - P_init = np.array(tmp[~tmp.mask], dtype=float) # initial pressure is first - - if 'press_conv' in self.pem_input and time_input == time: - P_init = P_init * self.pem_input['press_conv'] - - # extract saturations - if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: # This should be extended - for var in phases: - if var in ['WAT', 'GAS']: - tmp = self._get_pem_input('S{}'.format(var), time_input) - pem_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) - pem_input['S{}'.format(var)] = np.clip(pem_input['S{}'.format(var)], 0, 1) - - pem_input['SOIL'] = np.clip(1 - (pem_input['SWAT'] + pem_input['SGAS']), 0, 1) - saturations = [ np.clip(1 - (pem_input['SWAT'] + pem_input['SGAS']), 0, 1) if ph == 'OIL' else pem_input['S{}'.format(ph)] - for ph in phases] - elif 'WAT' in phases and 'GAS' in phases: # Smeaheia model using OPM CO2Store - for var in phases: - if var in ['GAS']: - tmp = self._get_pem_input('S{}'.format(var), time_input) - pem_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) - pem_input['S{}'.format(var)] = np.clip(pem_input['S{}'.format(var)] , 0, 1) - pem_input['SWAT'] = 1 - pem_input['SGAS'] - saturations = [1 - (pem_input['SGAS']) if ph == 'WAT' else pem_input['S{}'.format(ph)] for ph in phases] - - elif 'OIL' in phases and 'GAS' in phases: # Original Smeaheia model - for var in phases: - if var in ['GAS']: - tmp = self._get_pem_input('S{}'.format(var), time_input) - pem_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) - pem_input['S{}'.format(var)] = np.clip(pem_input['S{}'.format(var)], 0, 1) - pem_input['SOIL'] = 1 - pem_input['SGAS'] - saturations = [1 - (pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)] for ph in phases] - - else: - - print('Type and number of fluids are unspecified in calc_pem') - - # fluid saturations in dictionary - # tmp_dyn_var = {f'S{ph}': saturations[i] for i, ph in enumerate(phases)} - for var in phases: - tmp_dyn_var[f'S{var}'] = pem_input[f'S{var}'] - - tmp_dyn_var['PRESSURE'] = pem_input['PRESSURE'] - self.dyn_var.extend([tmp_dyn_var]) - - if not self.no_flow: - keywords = self.ecl_case.arrays(time) - keywords = [s.strip() for s in keywords] # Remove leading/trailing spaces - #for key in self.all_data_types: - #if 'grav' in key: - densities = [] - for var in phases: - # fluid densities - dens = var + '_DEN' - if dens in keywords: - tmp = self._get_pem_input(dens, time_input) - pem_input[dens] = np.array(tmp[~tmp.mask], dtype=float) - # extract densities - densities.append(pem_input[dens]) - else: - densities = None - # pore volumes at each assimilation step - if 'RPORV' in keywords: - tmp = self._get_pem_input('RPORV', time_input) - pem_input['RPORV'] = np.array(tmp[~tmp.mask], dtype=float) - else: - densities = None - - # Get elastic parameters - if hasattr(self, 'ensemble_member') and (self.ensemble_member is not None) and \ - (self.ensemble_member >= 0): - self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'], - dens = densities, ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init, - ensembleMember=self.ensemble_member) - else: - self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'], - dens = densities, ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init) - - def setup_fwd_run(self, redund_sim): - super().setup_fwd_run(redund_sim=redund_sim) - - def run_fwd_sim(self, state, member_i, del_folder=True): - # The inherited simulator also has a run_fwd_sim. Call this. - self.ensemble_member = member_i - - # Check if dynamic variables are provided in the state. If that is the case, do not run flow simulator - if any('sgas' in key for key in state.keys()) or any('swat' in key for key in state.keys()) or any('pressure' in key for key in state.keys()): - self.state = {} - for key in state.keys(): - self.state[key] = state[key] - self.no_flow = True - #self.pred_data = self.extract_data(member_i) - #else: - self.pred_data = super().run_fwd_sim(state, member_i, del_folder=del_folder) - - return self.pred_data - - def call_sim(self, folder=None, wait_for_proc=False): - # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. - # Then, get the pem. - if not self.no_flow: - success = super().call_sim(folder, wait_for_proc) - else: - success = True - - if success: - self.ecl_case = ecl.EclipseCase( - 'En_' + str(self.ensemble_member) + os.sep + self.file + '.DATA') - phases = self.ecl_case.init.phases - self.dyn_var = [] - vintage = [] - # loop over seismic vintages - for v, assim_time in enumerate(self.pem_input['vintage']): - time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ - dt.timedelta(days=assim_time) - - self.calc_pem(time, v+1) - - # mask the bulk imp. to get proper dimensions - tmp_value = np.zeros(self.ecl_case.init.shape) - tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp - self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float, - mask=deepcopy(self.ecl_case.init.mask)) - # run filter - self.pem._filter() - vintage.append(deepcopy(self.pem.bulkimp)) - - if hasattr(self.pem, 'baseline'): # 4D measurement - base_time = dt.datetime(self.startDate['year'], self.startDate['month'], - self.startDate['day']) + dt.timedelta(days=self.pem.baseline) - - self.calc_pem(base_time, 0) - - # mask the bulk imp. to get proper dimensions - tmp_value = np.zeros(self.ecl_case.init.shape) - - tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp - # kill if values are inf or nan - assert not np.isnan(tmp_value).any() - assert not np.isinf(tmp_value).any() - self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float, - mask=deepcopy(self.ecl_case.init.mask)) - self.pem._filter() - - # 4D response - self.pem_result = [] - for i, elem in enumerate(vintage): - self.pem_result.append(elem - deepcopy(self.pem.bulkimp)) - else: - for i, elem in enumerate(vintage): - self.pem_result.append(elem) - - return success - - def extract_data(self, member): - # start by getting the data from the flow simulator - super().extract_data(member) - - # get the sim2seis from file - for prim_ind in self.l_prim: - # Loop over all keys in pred_data (all data types) - for key in self.all_data_types: - if key in ['bulkimp']: - if self.true_prim[1][prim_ind] in self.pem_input['vintage']: - v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) - self.pred_data[prim_ind][key] = self.pem_result[v].data.flatten() - -class flow_sim2seis(flow): - """ - 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 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) - - 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 - self.dyn_var = [] - vintage = [] - # loop over seismic vintages - for v, assim_time in enumerate(self.pem_input['vintage']): - time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ - dt.timedelta(days=assim_time) - - self.calc_pem(time) #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) - 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_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_rock, mixIn_multi_data): - 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(redund_sim=None) - - def run_fwd_sim(self, state, member_i, del_folder=True): - """ - Setup and run the AVO forward simulator. - - Parameters - ---------- - state : dict - Dictionary containing the ensemble state. - - member_i : int - Index of the ensemble member. any index < 0 (e.g., -1) means the ground truth in synthetic case studies - - del_folder : bool, optional - Boolean to determine if the ensemble folder should be deleted. Default is False. - """ - - # The inherited simulator also has a run_fwd_sim. Call this. - self.ensemble_member = member_i - #return super().run_fwd_sim(state, member_i, del_folder=del_folder) - - - self.pred_data = super().run_fwd_sim(state, member_i, del_folder=del_folder) - return self.pred_data - - def call_sim(self, folder=None, wait_for_proc=False, run_reservoir_model=None, save_folder=None): - # replace the sim2seis part (which is unusable) by avo based on Pylops - - if folder is None: - folder = self.folder - else: - self.folder = folder - - if not self.no_flow: - # call call_sim in flow class (skip flow_rock, go directly to flow which is a parent of flow_rock) - success = super(flow_rock, self).call_sim(folder, wait_for_proc) - else: - success = True - - if success: - self.get_avo_result(folder, save_folder) - - return success - - def get_avo_result(self, folder, save_folder): - - if self.no_flow: - grid_file = self.pem_input['grid'] - grid = np.load(grid_file) - zcorn = grid['ZCORN'] - dz = np.diff(zcorn[:, 0, :, 0, :, 0], axis=0) - # Extract the last layer - last_layer = dz[-1, :, :] - # Reshape to ensure it has the same number of dimensions - last_layer = last_layer.reshape(1, dz.shape[1], dz.shape[2]) - # Concatenate to the original array along the first axis - dz = np.concatenate([dz, last_layer], axis=0) - f_dim = [grid['DIMENS'][2], grid['DIMENS'][1], grid['DIMENS'][0]] - else: - self.ecl_case = ecl.EclipseCase(folder + os.sep + self.file + '.DATA') if folder[-1] != os.sep \ - else ecl.EclipseCase(folder + self.file + '.DATA') - grid = self.ecl_case.grid() - zcorn = grid['ZCORN'] - ecl_init = ecl.EclipseInit(folder + os.sep + self.file + '.DATA') if folder[-1] != os.sep \ - else ecl.EclipseCase(folder + self.file + '.DATA') - dz = ecl_init.cell_data('DZ') - f_dim = [ecl_init.init.nk, ecl_init.init.nj, ecl_init.init.ni] - - - # ecl_init = ecl.EclipseInit(ecl_case) - # f_dim = [self.ecl_case.init.nk, self.ecl_case.init.nj, self.ecl_case.init.ni] - #f_dim = [self.NZ, self.NY, self.NX] - # phases = self.ecl_case.init.phases - self.dyn_var = [] - # coarsening of avo data - should be given as input in pipt - step_x = 1 - step_y = 1 - if 'baseline' in self.pem_input: # 4D measurement - base_time = dt.datetime(self.startDate['year'], self.startDate['month'], - self.startDate['day']) + dt.timedelta(days=self.pem_input['baseline']) - - - self.calc_pem(base_time,0) - # vp, vs, density in reservoir - self.calc_velocities(folder, save_folder, grid, -1, f_dim) - - if not self.no_flow: - # vp, vs, density in reservoir - vp, vs, rho = self.calc_velocities(folder, save_folder, grid, 0, f_dim) - - # avo data - # self._calc_avo_props() - avo_data_baseline, Rpp_baseline, vp_baseline, vs_baseline = self._calc_avo_props_active_cells(grid, vp, vs, rho, dz, zcorn) - kept_data = avo_data_baseline[::step_x, ::step_y, :].copy() - avo_data_baseline[:] = np.nan - avo_data_baseline[::step_x, ::step_y, :] = kept_data - # TODO: check which order to use, - # need to correlate with pipt/toml input-file and compression code - avo_baseline = avo_data_baseline.flatten(order="C") - avo_baseline = avo_baseline[~np.isnan(avo_baseline)] - #rho_baseline = rho_sample - tmp = self._get_pem_input('PRESSURE', base_time) - PRESSURE_baseline = np.array(tmp[~tmp.mask], dtype=float) - tmp = self._get_pem_input('SGAS', base_time) - SGAS_baseline = np.array(tmp[~tmp.mask], dtype=float) - print('OPM flow is used') - else: - file_name = f"avo_vint0_{folder}.npz" if folder[-1] != os.sep \ - else f"avo_vint0_{folder[:-1]}.npz" - - avo_baseline = np.load(file_name, allow_pickle=True)['avo_bl'] - #Rpp_baseline = np.load(file_name, allow_pickle=True)['Rpp_bl'] - #vs_baseline = np.load(file_name, allow_pickle=True)['Vs_bl'] - #vp_baseline = np.load(file_name, allow_pickle=True)['Vp_bl'] - #rho_baseline = np.load(file_name, allow_pickle=True)['Rho_bl'] - - vintage = [] - # loop over seismic vintages - for v, assim_time in enumerate(self.pem_input['vintage']): - time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ - dt.timedelta(days=assim_time) - - # extract dynamic variables from simulation run - self.calc_pem(time, v+1) - - # vp, vs, density in reservoir - vp, vs, rho = self.calc_velocities(folder, save_folder, grid, v+1, f_dim) - - # avo data - #self._calc_avo_props() - avo_data, Rpp, vp_sample, vs_sample = self._calc_avo_props_active_cells(grid, vp, vs, rho, dz, zcorn) - #avo_data = avo_data[::step_x, ::step_y, :] - # make mask for wavelet decomposition - kept_data = avo_data[::step_x, ::step_y, :].copy() - avo_data[:] = np.nan - avo_data[::step_x, ::step_y, :] = kept_data - mask = np.ones(np.shape(avo_data), dtype=bool) - mask[np.isnan(avo_data)]=False - np.savez(f'mask_{v}.npz', mask=mask) - #TODO: check order of flattening as above - avo = avo_data.flatten(order="C") - avo = avo[~np.isnan(avo)] - - tmp = self._get_pem_input('PRESSURE', time) - PRESSURE = np.array(tmp[~tmp.mask], dtype=float) - tmp = self._get_pem_input('SGAS', time) - SGAS = np.array(tmp[~tmp.mask], dtype=float) - - - # MLIE: implement 4D avo - if 'baseline' in self.pem_input: # 4D measurement - avo = avo - avo_baseline - Rpp = Rpp - Rpp_baseline - vs_sample = vs_sample - vs_baseline - vp_sample = vp_sample - vp_baseline - PRESSURE = PRESSURE - PRESSURE_baseline - SGAS = SGAS - SGAS_baseline - #rho = self.rho_sample - rho_baseline - print('Time-lapse avo') - #else: - # Rpp = self.Rpp - # Vs = self.vs_sample - # Vp = self.vp_sample - # rho = self.rho_sample - - - - # XLUO: self.ensemble_member < 0 => reference reservoir model in synthetic case studies - # the corresonding (noisy) data are observations in data assimilation - if 'add_synthetic_noise' in self.input_dict and self.ensemble_member < 0: - non_nan_idx = np.argwhere(~np.isnan(avo)) - data_std = np.std(avo[non_nan_idx]) - if self.input_dict['add_synthetic_noise'][0] == 'snr': - noise_std = np.sqrt(self.input_dict['add_synthetic_noise'][1]) * data_std - avo[non_nan_idx] += noise_std * np.random.randn(avo[non_nan_idx].size, 1) - else: - noise_std = 0.0 # simulated data don't contain noise - - vintage.append(deepcopy(avo)) - - if v == 0: - save_dic = {'avo': avo, 'noise_std': noise_std, 'Rpp': Rpp, 'Vs': vs_sample, 'Vp': vp_sample, 'PRESSURE': PRESSURE, 'SGAS': SGAS, **self.avo_config} - #save_dic = {'avo': avo, 'noise_std': noise_std, 'Rpp': Rpp, 'Vs': vs_sample, 'Vp': vp_sample, 'rho': rho_sample, #**self.avo_config, - # 'Vs_bl': vs_baseline, 'Vp_bl': vp_baseline, 'avo_bl': avo_baseline, 'Rpp_bl': Rpp_baseline, 'rho_bl': rho_baseline, **self.avo_config} - #save_dic = {'avo': avo, 'noise_std': noise_std, 'Rpp': Rpp, 'Vs': vs_sample, 'Vp': vp_sample, - # **self.avo_config} - else: - save_dic = {'avo': avo, 'noise_std': noise_std, 'Rpp': Rpp, 'Vs': vs_sample, 'Vp': vp_sample, 'PRESSURE': PRESSURE, 'SGAS': SGAS}#, 'Rpp': Rpp, 'Vs': vs_sample, 'Vp': vp_sample, - #'rho': rho_sample, **self.avo_config} - - if save_folder is not None: - file_name = save_folder + os.sep + f"avo_vint{v}.npz" if save_folder[-1] != os.sep \ - else save_folder + f"avo_vint{v}.npz" - #np.savez(file_name, **save_dic) - else: - file_name = folder + os.sep + f"avo_vint{v}.npz" if folder[-1] != os.sep \ - else folder + f"avo_vint{v}.npz" - file_name_rec = 'Ensemble_results/' + f"avo_vint{v}_{folder}.npz" if folder[-1] != os.sep \ - else 'Ensemble_results/' + f"avo_vint{v}_{folder[:-1]}.npz" - np.savez(file_name_rec, **save_dic) - # with open(file_name, "wb") as f: - # dump(**save_dic, f) - np.savez(file_name, **save_dic) - # 4D response - self.avo_result = [] - for i, elem in enumerate(vintage): - self.avo_result.append(elem) - - def calc_velocities(self, folder, save_folder, grid, v, f_dim): - # The properties in pem are only given in the active cells - # indices of active cells: - - true_indices = np.where(grid['ACTNUM']) - - vp = np.full(f_dim, np.nan) - vp[true_indices] = (self.pem.getBulkVel()) - vs = np.full(f_dim, np.nan) - vs[true_indices] = (self.pem.getShearVel()) - rho = np.full(f_dim, np.nan) - rho[true_indices] = (self.pem.getDens()) - - - - ## Debug - #bulkmod = np.full(f_dim, np.nan) - #bulkmod[true_indices] = self.pem.getBulkMod() - #self.shearmod = np.full(f_dim, np.nan) - #self.shearmod[true_indices] = self.pem.getShearMod() - #self.poverburden = np.full(f_dim, np.nan) - #self.poverburden[true_indices] = self.pem.getOverburdenP() - #self.pressure = np.full(f_dim, np.nan) - #self.pressure[true_indices] = self.pem.getPressure() - #self.peff = np.full(f_dim, np.nan) - #self.peff[true_indices] = self.pem.getPeff() - #porosity = np.full(f_dim, np.nan) - #porosity[true_indices] = self.pem.getPorosity() - #if self.dyn_var: - # sgas = np.full(f_dim, np.nan) - # sgas[true_indices] = self.dyn_var[v]['SGAS'] - # soil = np.full(f_dim, np.nan) - # soil[true_indices] = self.dyn_var[v]['SOIL'] - # pdyn = np.full(f_dim, np.nan) - # pdyn[true_indices] = self.dyn_var[v]['PRESSURE'] - # - #if self.dyn_var is None: - # save_dic = {'vp': vp, 'vs': vs, 'rho': rho}#, 'bulkmod': self.bulkmod, 'shearmod': self.shearmod, - # #'Pov': self.poverburden, 'P': self.pressure, 'Peff': self.peff, 'por': porosity} # for debugging - #else: - # save_dic = {'vp': vp, 'vs': vs, 'rho': rho}#, 'por': porosity, 'sgas': sgas, 'Pd': pdyn} - - #if save_folder is not None: - # file_name = save_folder + os.sep + f"vp_vs_rho_vint{v}.npz" if save_folder[-1] != os.sep \ - # else save_folder + f"vp_vs_rho_vint{v}.npz" - # np.savez(file_name, **save_dic) - #else: - # file_name_rec = 'Ensemble_results/' + f"vp_vs_rho_vint{v}_{folder}.npz" if folder[-1] != os.sep \ - # else 'Ensemble_results/' + f"vp_vs_rho_vint{v}_{folder[:-1]}.npz" - # np.savez(file_name_rec, **save_dic) - # for debugging - return vp, vs, rho - - def extract_data(self, member): - # start by getting the data from the flow simulator - super(flow_rock, self).extract_data(member) - - # get the avo from file - for prim_ind in self.l_prim: - # Loop over all keys in pred_data (all data types) - for key in self.all_data_types: - if 'avo' in key: - if self.true_prim[1][prim_ind] in self.pem_input['vintage']: - idx = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) - filename = self.folder + os.sep + key + '_vint' + str(idx) + '.npz' if self.folder[-1] != os.sep \ - else self.folder + key + '_vint' + str(idx) + '.npz' - with np.load(filename) as f: - self.pred_data[prim_ind][key] = f[key] - # - #v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) - #self.pred_data[prim_ind][key] = self.avo_result[v].flatten() - - def _get_avo_info(self, avo_config=None): - """ - AVO configuration - """ - # list of configuration parameters in the "AVO" section - config_para_list = ['dz', 'tops', 'angle', 'frequency', 'wave_len', 'vp_shale', 'vs_shale', - 'den_shale', 't_min', 't_max', 't_sampling', 'pp_func'] - if 'avo' in self.input_dict: - self.avo_config = {} - for elem in self.input_dict['avo']: - assert elem[0] in config_para_list, f'Property {elem[0]} not supported' - if elem[0] == 'vintage' and not isinstance(elem[1], list): - elem[1] = [elem[1]] - self.avo_config[elem[0]] = elem[1] - - # if only one angle is considered, convert self.avo_config['angle'] into a list, as required later - if isinstance(self.avo_config['angle'], float): - self.avo_config['angle'] = [self.avo_config['angle']] - - # self._get_DZ(file=self.avo_config['dz']) # =>self.DZ - kw_file = {'DZ': self.avo_config['dz'], 'TOPS': self.avo_config['tops']} - self._get_props(kw_file) - self.overburden = self.pem_input['overburden'] - - # make sure that the "pylops" package is installed - # See https://github.com/PyLops/pylops - self.pp_func = getattr(import_module('pylops.avo.avo'), self.avo_config['pp_func']) - - else: - self.avo_config = None - - def _get_props(self, kw_file): - # extract properties (specified by keywords) in (possibly) different files - # kw_file: a dictionary contains "keyword: file" pairs - # Note that all properties are reshaped into the reservoir model dimension (NX, NY, NZ) - # using the "F" order - for kw in kw_file: - file = kw_file[kw] - if file.endswith('.npz'): - with np.load(file) as f: - exec(f'self.{kw} = f[ "{kw}" ]') - self.NX, self.NY, self.NZ = f['NX'], f['NY'], f['NZ'] - else: - try: - self.NX = int(self.input_dict['dimension'][0]) - self.NY = int(self.input_dict['dimension'][1]) - self.NZ = int(self.input_dict['dimension'][2]) - except: - for item in self.input_dict['pem']: - if item[0] == 'dimension': - dimension = item[1] - break - self.NX = int(dimension[0]) - self.NY = int(dimension[1]) - self.NZ = int(dimension[2]) - # reader = GRDECL_Parser(filename=file) - # reader.read_GRDECL() - # exec(f"self.{kw} = reader.{kw}.reshape((reader.NX, reader.NY, reader.NZ), order='F')") - # self.NX, self.NY, self.NZ = reader.NX, reader.NY, reader.NZ - # eval(f'np.savez("./{kw}.npz", {kw}=self.{kw}, NX=self.NX, NY=self.NY, NZ=self.NZ)') - - def _calc_avo_props(self, dt=0.0005): - # dt is the fine resolution sampling rate - # convert properties in reservoir model to time domain - vp_shale = self.avo_config['vp_shale'] # scalar value (code may not work for matrix value) - vs_shale = self.avo_config['vs_shale'] # scalar value - rho_shale = self.avo_config['den_shale'] # scalar value - - # Two-way travel time of the top of the reservoir - # TOPS[:, :, 0] corresponds to the depth profile of the reservoir top on the first layer - top_res = 2 * self.TOPS[:, :, 0] / vp_shale - - # Cumulative traveling time through the reservoir in vertical direction - cum_time_res = np.cumsum(2 * self.DZ / self.vp, axis=2) + top_res[:, :, np.newaxis] - - # assumes underburden to be constant. No reflections from underburden. Hence set traveltime to underburden very large - underburden = top_res + np.max(cum_time_res) - - # total travel time - # cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res), axis=2) - cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res, underburden[:, :, np.newaxis]), axis=2) - - - # add overburden and underburden of Vp, Vs and Density - vp = np.concatenate((vp_shale * np.ones((self.NX, self.NY, 1)), - self.vp, vp_shale * np.ones((self.NX, self.NY, 1))), axis=2) - vs = np.concatenate((vs_shale * np.ones((self.NX, self.NY, 1)), - self.vs, vs_shale * np.ones((self.NX, self.NY, 1))), axis=2) - - #rho = np.concatenate((rho_shale * np.ones((self.NX, self.NY, 1)) * 0.001, # kg/m^3 -> k/cm^3 - # self.rho, rho_shale * np.ones((self.NX, self.NY, 1)) * 0.001), axis=2) - rho = np.concatenate((rho_shale * np.ones((self.NX, self.NY, 1)), - self.rho, rho_shale * np.ones((self.NX, self.NY, 1))), axis=2) - - # search for the lowest grid cell thickness and sample the time according to - # that grid thickness to preserve the thin layer effect - time_sample = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], dt) - if time_sample.shape[0] == 1: - time_sample = time_sample.reshape(-1) - time_sample = np.tile(time_sample, (self.NX, self.NY, 1)) - - vp_sample = np.tile(vp[:, :, 1][..., np.newaxis], (1, 1, time_sample.shape[2])) - vs_sample = np.tile(vs[:, :, 1][..., np.newaxis], (1, 1, time_sample.shape[2])) - rho_sample = np.tile(rho[:, :, 1][..., np.newaxis], (1, 1, time_sample.shape[2])) - - for m in range(self.NX): - for l in range(self.NY): - for k in range(time_sample.shape[2]): - # find the right interval of time_sample[m, l, k] belonging to, and use - # this information to allocate vp, vs, rho - idx = np.searchsorted(cum_time[m, l, :], time_sample[m, l, k], side='left') - idx = idx if idx < len(cum_time[m, l, :]) else len(cum_time[m, l, :]) - 1 - vp_sample[m, l, k] = vp[m, l, idx] - vs_sample[m, l, k] = vs[m, l, idx] - rho_sample[m, l, k] = rho[m, l, idx] - - - - - # Ricker wavelet - wavelet, t_axis, wav_center = ricker(np.arange(0, self.avo_config['wave_len'], dt), - f0=self.avo_config['frequency']) - - - # Travel time corresponds to reflectivity series - t = time_sample[:, :, 0:-1] - - # interpolation time - t_interp = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], self.avo_config['t_sampling']) - trace_interp = np.zeros((self.NX, self.NY, len(t_interp))) - - # number of pp reflection coefficients in the vertical direction - - nz_rpp = vp_sample.shape[2] - 1 - - for i in range(len(self.avo_config['angle'])): - angle = self.avo_config['angle'][i] - Rpp = self.pp_func(vp_sample[:, :, :-1], vs_sample[:, :, :-1], rho_sample[:, :, :-1], - vp_sample[:, :, 1:], vs_sample[:, :, 1:], rho_sample[:, :, 1:], angle) - - for m in range(self.NX): - for l in range(self.NY): - # convolution with the Ricker wavelet - conv_op = Convolve1D(nz_rpp, h=wavelet, offset=wav_center, dtype="float32") - w_trace = conv_op * Rpp[m, l, :] - - # Sample the trace into regular time interval - f = interp1d(np.squeeze(t[m, l, :]), np.squeeze(w_trace), - kind='nearest', fill_value='extrapolate') - trace_interp[m, l, :] = f(t_interp) - - if i == 0: - avo_data = trace_interp # 3D - elif i == 1: - avo_data = np.stack((avo_data, trace_interp), axis=-1) # 4D - else: - avo_data = np.concatenate((avo_data, trace_interp[:, :, :, np.newaxis]), axis=3) # 4D - - self.avo_data = avo_data - - def _calc_avo_props_active_cells(self, grid, vp, vs, rho, dz, zcorn, dt=0.0005): - # dt is the fine resolution sampling rate - # convert properties in reservoir model to time domain - vp_shale = self.avo_config['vp_shale'] # scalar value (code may not work for matrix value) - vs_shale = self.avo_config['vs_shale'] # scalar value - rho_shale = self.avo_config['den_shale'] # scalar value - - - actnum = grid['ACTNUM'] - # Find indices where the boolean array is True - active_indices = np.where(actnum) - c, a, b = active_indices - - # Two-way travel time tp the top of the reservoir - top_res = 2 * zcorn[0, 0, :, 0, :, 0] / vp_shale - - # depth difference between cells in z-direction: - depth_differences = dz#(zcorn[:, 0, :, 0, :, 0] , axis=0) - - - # Cumulative traveling time through the reservoir in vertical direction - #cum_time_res = 2 * zcorn[:, 0, :, 0, :, 0] / self.vp + top_res[np.newaxis, :, :] - cum_time_res = np.cumsum(2 * depth_differences / vp, axis=0) + top_res[np.newaxis, :, :] - # assumes under burden to be constant. No reflections from under burden. Hence set travel time to under burden very large - underburden = top_res + np.nanmax(cum_time_res) - - # total travel time - # cum_time = np.concat enate((top_res[:, :, np.newaxis], cum_time_res), axis=2) - cum_time = np.concatenate((top_res[np.newaxis, :, :], cum_time_res, underburden[np.newaxis, :, :]), axis=0) - - # add overburden and underburden values for Vp, Vs and Density - vp = np.concatenate((vp_shale * np.ones((1, self.NY, self.NX)), - vp, vp_shale * np.ones((1, self.NY, self.NX))), axis=0) - vs = np.concatenate((vs_shale * np.ones((1, self.NY, self.NX)), - vs, vs_shale * np.ones((1, self.NY, self.NX))), axis=0) - rho = np.concatenate((rho_shale * np.ones((1, self.NY, self.NX)), - rho, rho_shale * np.ones((1, self.NY, self.NX))), axis=0) - - - # Combine a and b into a 2D array (each column represents a vector) - ab = np.column_stack((a, b)) - - # Extract unique rows and get the indices of those unique rows - unique_rows, indices = np.unique(ab, axis=0, return_index=True) - - # search for the lowest grid cell thickness and sample the time according to - # that grid thickness to preserve the thin layer effect - time_sample = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], dt) - if time_sample.shape[0] == 1: - time_sample = time_sample.reshape(-1) - #time_sample = np.tile(time_sample, (len(indices), 1)) - time_sample = np.tile(time_sample, (self.NX, self.NY, 1)) - - #vp_sample = vp_shale * np.ones((self.NX, self.NY, time_sample.shape[2])) - #vs_sample = vs_shale * np.ones((self.NX, self.NY, time_sample.shape[2])) - #rho_sample = rho_shale * np.ones((self.NX, self.NY, time_sample.shape[2])) - vp_sample = np.full([self.NX, self.NY, time_sample.shape[2]], np.nan) - vs_sample = np.full([self.NX, self.NY, time_sample.shape[2]], np.nan) - rho_sample = np.full([self.NX, self.NY, time_sample.shape[2]], np.nan) - - - for ind in range(len(indices)): - for k in range(time_sample.shape[2]): - # find the right interval of time_sample[m, l, k] belonging to, and use - # this information to allocate vp, vs, rho - idx = np.searchsorted(cum_time[:, a[indices[ind]], b[indices[ind]]], time_sample[b[indices[ind]], a[indices[ind]], k], side='left') - idx = idx if idx < len(cum_time[:, a[indices[ind]], b[indices[ind]]]) else len( - cum_time[:,a[indices[ind]], b[indices[ind]]]) - 1 - vp_sample[b[indices[ind]], a[indices[ind]], k] = vp[idx, a[indices[ind]], b[indices[ind]]] - vs_sample[b[indices[ind]], a[indices[ind]], k] = vs[idx, a[indices[ind]], b[indices[ind]]] - rho_sample[b[indices[ind]], a[indices[ind]], k] = rho[idx, a[indices[ind]], b[indices[ind]]] - - # Ricker wavelet - wavelet, t_axis, wav_center = ricker(np.arange(0, self.avo_config['wave_len']-dt, dt), - f0=self.avo_config['frequency']) - - # Travel time corresponds to reflectivity series - #t = time_sample[:, 0:-1] - t = time_sample[:, :, 0:-1] - - # interpolation time - t_interp = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], self.avo_config['t_sampling']) - #trace_interp = np.zeros((len(indices), len(t_interp))) - #trace_interp = np.zeros((self.NX, self.NY, len(t_interp))) - trace_interp = np.full([self.NX, self.NY, len(t_interp)], np.nan) - - # number of pp reflection coefficients in the vertical direction - nz_rpp = vp_sample.shape[2] - 1 - conv_op = Convolve1D(nz_rpp, h=wavelet, offset=wav_center, dtype="float32") - - avo_data = [] - Rpp = [] - for i in range(len(self.avo_config['angle'])): - angle = self.avo_config['angle'][i] - Rpp = self.pp_func(vp_sample[:, :, :-1], vs_sample[:, :, :-1], rho_sample[:, :, :-1], - vp_sample[:, :, 1:], vs_sample[:, :, 1:], rho_sample[:, :, 1:], angle) - - for ind in range(len(indices)): - # convolution with the Ricker wavelet - - w_trace = conv_op * Rpp[b[indices[ind]], a[indices[ind]], :] - - # Sample the trace into regular time interval - f = interp1d(np.squeeze(t[b[indices[ind]], a[indices[ind]], :]), np.squeeze(w_trace), - kind='nearest', fill_value='extrapolate') - trace_interp[b[indices[ind]], a[indices[ind]], :] = f(t_interp) - - if i == 0: - avo_data = trace_interp # 3D - elif i == 1: - avo_data = np.stack((avo_data, trace_interp), axis=-1) # 4D - else: - avo_data = np.concatenate((avo_data, trace_interp[:, :, np.newaxis]), axis=2) # 4D - - return avo_data, Rpp, vp_sample, vs_sample - #self.avo_data = avo_data - #self.Rpp = Rpp - #self.vp_sample = vp_sample - #self.vs_sample = vs_sample - #self.rho_sample = rho_sample - - - - def _calc_avo_props_active_cells_org(self, grid, vp, vs, rho, dz, zcorn, dt=0.0005): - # dt is the fine resolution sampling rate - # convert properties in reservoir model to time domain - vp_shale = self.avo_config['vp_shale'] # scalar value (code may not work for matrix value) - vs_shale = self.avo_config['vs_shale'] # scalar value - rho_shale = self.avo_config['den_shale'] # scalar value - - - actnum = grid['ACTNUM'] - # Find indices where the boolean array is True - active_indices = np.where(actnum) - # # # - - # Two-way travel time tp the top of the reservoir - - - c, a, b = active_indices - - # Two-way travel time tp the top of the reservoir - top_res = 2 * zcorn[0, 0, :, 0, :, 0] / vp_shale - - # depth difference between cells in z-direction: - depth_differences = dz#(zcorn[:, 0, :, 0, :, 0] , axis=0) - - - # Cumulative traveling time through the reservoir in vertical direction - #cum_time_res = 2 * zcorn[:, 0, :, 0, :, 0] / self.vp + top_res[np.newaxis, :, :] - cum_time_res = np.cumsum(2 * depth_differences / vp, axis=0) + top_res[np.newaxis, :, :] - # assumes under burden to be constant. No reflections from under burden. Hence set travel time to under burden very large - underburden = top_res + np.nanmax(cum_time_res) - - # total travel time - # cum_time = np.concat enate((top_res[:, :, np.newaxis], cum_time_res), axis=2) - cum_time = np.concatenate((top_res[np.newaxis, :, :], cum_time_res, underburden[np.newaxis, :, :]), axis=0) - - # add overburden and underburden values for Vp, Vs and Density - vp = np.concatenate((vp_shale * np.ones((1, self.NY, self.NX)), - vp, vp_shale * np.ones((1, self.NY, self.NX))), axis=0) - vs = np.concatenate((vs_shale * np.ones((1, self.NY, self.NX)), - vs, vs_shale * np.ones((1, self.NY, self.NX))), axis=0) - rho = np.concatenate((rho_shale * np.ones((1, self.NY, self.NX)), - rho, rho_shale * np.ones((1, self.NY, self.NX))), axis=0) - - - # Combine a and b into a 2D array (each column represents a vector) - ab = np.column_stack((a, b)) - - # Extract unique rows and get the indices of those unique rows - unique_rows, indices = np.unique(ab, axis=0, return_index=True) - - # search for the lowest grid cell thickness and sample the time according to - # that grid thickness to preserve the thin layer effect - time_sample = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], dt) - if time_sample.shape[0] == 1: - time_sample = time_sample.reshape(-1) - time_sample = np.tile(time_sample, (len(indices), 1)) - - vp_sample = np.zeros((len(indices), time_sample.shape[1])) - vs_sample = np.zeros((len(indices), time_sample.shape[1])) - rho_sample = np.zeros((len(indices), time_sample.shape[1])) - - for ind in range(len(indices)): - for k in range(time_sample.shape[1]): - # find the right interval of time_sample[m, l, k] belonging to, and use - # this information to allocate vp, vs, rho - idx = np.searchsorted(cum_time[:, a[indices[ind]], b[indices[ind]]], time_sample[ind, k], side='left') - idx = idx if idx < len(cum_time[:, a[indices[ind]], b[indices[ind]]]) else len( - cum_time[:,a[indices[ind]], b[indices[ind]]]) - 1 - vp_sample[ind, k] = vp[idx, a[indices[ind]], b[indices[ind]]] - vs_sample[ind, k] = vs[idx, a[indices[ind]], b[indices[ind]]] - rho_sample[ind, k] = rho[idx, a[indices[ind]], b[indices[ind]]] - - # Ricker wavelet - wavelet, t_axis, wav_center = ricker(np.arange(0, self.avo_config['wave_len']-dt, dt), - f0=self.avo_config['frequency']) - - # Travel time corresponds to reflectivity series - t = time_sample[:, 0:-1] - - # interpolation time - t_interp = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], self.avo_config['t_sampling']) - trace_interp = np.zeros((len(indices), len(t_interp))) - - # number of pp reflection coefficients in the vertical direction - nz_rpp = vp_sample.shape[1] - 1 - conv_op = Convolve1D(nz_rpp, h=wavelet, offset=wav_center, dtype="float32") - - avo_data = [] - Rpp = [] - for i in range(len(self.avo_config['angle'])): - angle = self.avo_config['angle'][i] - Rpp = self.pp_func(vp_sample[:, :-1], vs_sample[:, :-1], rho_sample[:, :-1], - vp_sample[:, 1:], vs_sample[:, 1:], rho_sample[:, 1:], angle) - - for ind in range(len(indices)): - # convolution with the Ricker wavelet - - w_trace = conv_op * Rpp[ind, :] - - # Sample the trace into regular time interval - f = interp1d(np.squeeze(t[ind, :]), np.squeeze(w_trace), - kind='nearest', fill_value='extrapolate') - trace_interp[ind, :] = f(t_interp) - - if i == 0: - avo_data = trace_interp # 3D - elif i == 1: - avo_data = np.stack((avo_data, trace_interp), axis=-1) # 4D - else: - avo_data = np.concatenate((avo_data, trace_interp[:, :, np.newaxis]), axis=2) # 4D - - return avo_data, Rpp, vp_sample, vs_sample, rho_sample - #self.avo_data = avo_data - #self.Rpp = Rpp - #self.vp_sample = vp_sample - #self.vs_sample = vs_sample - #self.rho_sample = rho_sample - - def _calc_avo_props_active_cells_org(self, grid, dt=0.0005): - # dt is the fine resolution sampling rate - # convert properties in reservoir model to time domain - vp_shale = self.avo_config['vp_shale'] # scalar value (code may not work for matrix value) - vs_shale = self.avo_config['vs_shale'] # scalar value - rho_shale = self.avo_config['den_shale'] # scalar value - - # check if Nz, is at axis = 0, then transpose to dimensions, Nx, ny, Nz - if grid['ACTNUM'].shape[0] == self.NZ: - vp = np.transpose(self.vp, (2, 1, 0)) - vs = np.transpose(self.vs, (2, 1, 0)) - rho = np.transpose(self.rho, (2, 1, 0)) - actnum = np.transpose(grid['ACTNUM'], (2, 1, 0)) - else: - actnum = grid['ACTNUM'] - vp = self.vp - vs = self.vs - rho = self.rho - # # # - - # Two-way travel time of the top of the reservoir - # TOPS[:, :, 0] corresponds to the depth profile of the reservoir top on the first layer - top_res = 2 * self.TOPS[:, :, 0] / vp_shale - - # Cumulative traveling time through the reservoir in vertical direction - cum_time_res = np.nancumsum(2 * self.DZ / vp, axis=2) + top_res[:, :, np.newaxis] - - # assumes under burden to be constant. No reflections from under burden. Hence set travel time to under burden very large - underburden = top_res + np.max(cum_time_res) - - # total travel time - # cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res), axis=2) - cum_time = np.concatenate((top_res[:, :, np.newaxis], cum_time_res, underburden[:, :, np.newaxis]), axis=2) - - # add overburden and underburden of Vp, Vs and Density - vp = np.concatenate((vp_shale * np.ones((self.NX, self.NY, 1)), - vp, vp_shale * np.ones((self.NX, self.NY, 1))), axis=2) - vs = np.concatenate((vs_shale * np.ones((self.NX, self.NY, 1)), - vs, vs_shale * np.ones((self.NX, self.NY, 1))), axis=2) - #rho = np.concatenate((rho_shale * np.ones((self.NX, self.NY, 1)) * 0.001, # kg/m^3 -> k/cm^3 - # self.rho, rho_shale * np.ones((self.NX, self.NY, 1)) * 0.001), axis=2) - rho = np.concatenate((rho_shale * np.ones((self.NX, self.NY, 1)), - rho, rho_shale * np.ones((self.NX, self.NY, 1))), axis=2) - - # get indices of active cells - - indices = np.where(actnum) - a, b, c = indices - # Combine a and b into a 2D array (each column represents a vector) - ab = np.column_stack((a, b)) - - # Extract unique rows and get the indices of those unique rows - unique_rows, indices = np.unique(ab, axis=0, return_index=True) - - - # search for the lowest grid cell thickness and sample the time according to - # that grid thickness to preserve the thin layer effect - time_sample = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], dt) - if time_sample.shape[0] == 1: - time_sample = time_sample.reshape(-1) - time_sample = np.tile(time_sample, (len(indices), 1)) - - vp_sample = np.zeros((len(indices), time_sample.shape[1])) - vs_sample = np.zeros((len(indices), time_sample.shape[1])) - rho_sample = np.zeros((len(indices), time_sample.shape[1])) - - - for ind in range(len(indices)): - for k in range(time_sample.shape[1]): - # find the right interval of time_sample[m, l, k] belonging to, and use - # this information to allocate vp, vs, rho - idx = np.searchsorted(cum_time[a[indices[ind]], b[indices[ind]], :], time_sample[ind, k], side='left') - idx = idx if idx < len(cum_time[a[indices[ind]], b[indices[ind]], :]) else len(cum_time[a[indices[ind]], b[indices[ind]], :]) - 1 - vp_sample[ind, k] = vp[a[indices[ind]], b[indices[ind]], idx] - vs_sample[ind, k] = vs[a[indices[ind]], b[indices[ind]], idx] - rho_sample[ind, k] = rho[a[indices[ind]], b[indices[ind]], idx] - - - # Ricker wavelet - wavelet, t_axis, wav_center = ricker(np.arange(0, self.avo_config['wave_len'], dt), - f0=self.avo_config['frequency']) - - # Travel time corresponds to reflectivity series - t = time_sample[:, 0:-1] - - # interpolation time - t_interp = np.arange(self.avo_config['t_min'], self.avo_config['t_max'], self.avo_config['t_sampling']) - trace_interp = np.zeros((len(indices), len(t_interp))) - - # number of pp reflection coefficients in the vertical direction - nz_rpp = vp_sample.shape[1] - 1 - conv_op = Convolve1D(nz_rpp, h=wavelet, offset=wav_center, dtype="float32") - - avo_data = [] - Rpp = [] - for i in range(len(self.avo_config['angle'])): - angle = self.avo_config['angle'][i] - Rpp = self.pp_func(vp_sample[:, :-1], vs_sample[:, :-1], rho_sample[:, :-1], - vp_sample[:, 1:], vs_sample[:, 1:], rho_sample[:, 1:], angle) - - - - for ind in range(len(indices)): - # convolution with the Ricker wavelet - - w_trace = conv_op * Rpp[ind, :] - - # Sample the trace into regular time interval - f = interp1d(np.squeeze(t[ind, :]), np.squeeze(w_trace), - kind='nearest', fill_value='extrapolate') - trace_interp[ind, :] = f(t_interp) - - if i == 0: - avo_data = trace_interp # 3D - elif i == 1: - avo_data = np.stack((avo_data, trace_interp), axis=-1) # 4D - else: - avo_data = np.concatenate((avo_data, trace_interp[:, :, :, np.newaxis]), axis=3) # 4D - - self.avo_data = avo_data - self.Rpp = Rpp - self.vp_sample = vp_sample - self.vs_sample = vs_sample - self.rho_sample = rho_sample - - - @classmethod - def _reformat3D_then_flatten(cls, array, flatten=True, order="F"): - """ - XILU: Quantities read by "EclipseData.cell_data" are put in the axis order of [nz, ny, nx]. To be consisent with - ECLIPSE/OPM custom, we need to change the axis order. We further flatten the array according to the specified order - """ - array = np.array(array) - if len(array.shape) != 1: # if array is a 1D array, then do nothing - assert isinstance(array, np.ndarray) and len(array.shape) == 3, "Only 3D numpy array are supported" - - # axis [0 (nz), 1 (ny), 2 (nx)] -> [2 (nx), 1 (ny), 0 (nz)] - new_array = np.transpose(array, axes=[2, 1, 0]) - if flatten: - new_array = new_array.flatten(order=order) - - return new_array - else: - return array - -class flow_grav(flow_rock, mixIn_multi_data): - def __init__(self, input_dict=None, filename=None, options=None, **kwargs): - super().__init__(input_dict, filename, options) - - self.grav_input = {} - assert 'grav' in input_dict, 'To do GRAV simulation, please specify an "GRAV" section in the "FWDSIM" part' - self._get_grav_info() - - def setup_fwd_run(self, **kwargs): - self.__dict__.update(kwargs) - - super().setup_fwd_run(redund_sim=None) - - def run_fwd_sim(self, state, member_i, del_folder=True): - # The inherited simulator also has a run_fwd_sim. Call this. - self.ensemble_member = member_i - #return super().run_fwd_sim(state, member_i, del_folder=del_folder) - self.pred_data = super().run_fwd_sim(state, member_i, del_folder) - return self.pred_data - - def call_sim(self, folder=None, wait_for_proc=False, save_folder=None): - # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. - # Then, get the pem. - if folder is None: - folder = self.folder - - if not self.no_flow: - # call call_sim in flow class (skip flow_rock, go directly to flow which is a parent of flow_rock) - success = super(flow_rock, self).call_sim(folder, True) - else: - success = True - # - # use output from flow simulator to forward model gravity response - if success: - self.get_grav_result(folder, save_folder) - - return success - - def get_grav_result(self, folder, save_folder): - if self.no_flow: - grid_file = self.pem_input['grid'] - grid = np.load(grid_file) - else: - self.ecl_case = ecl.EclipseCase(folder + os.sep + self.file + '.DATA') if folder[-1] != os.sep \ - else ecl.EclipseCase(folder + self.file + '.DATA') - grid = self.ecl_case.grid() - - - self.dyn_var = [] - - # cell centers - #cell_centre = self.find_cell_centre(grid) - - # receiver locations - # Make a mesh of the area - pad = self.grav_config.get('padding', 1500) # 3 km padding around the reservoir - if 'padding' not in self.grav_config: - print('Please specify extent of measurement locations (Padding in input file), using 1.5 km as default') - dxy = self.grav_config.get('grid_spacing', 1500) # - if 'grid_spacing' not in self.grav_config: - print('Please specify grid spacing in input file, using 1.5 km as default') - if 'seabed' in self.grav_config and self.grav_config['seabed'] is not None: - file_path = self.grav_config['seabed'] - water_depth = self.get_seabed_depths(file_path) - else: - water_depth = self.grav_config.get('water_depth', 300) - if 'water_depth' not in self.grav_config: - print('Please specify water depths in input file, using 300 m as default') - pos = self.measurement_locations(grid, water_depth, pad, dxy) - - # loop over vintages with gravity acquisitions - grav_struct = {} - - if 'baseline' in self.grav_config: # 4D measurement - base_time = dt.datetime(self.startDate['year'], self.startDate['month'], - self.startDate['day']) + dt.timedelta(days=self.grav_config['baseline']) - # porosity, saturation, densities, and fluid mass at time of baseline survey - grav_base = self.calc_mass(base_time, 0) - - - else: - # seafloor gravity only works in 4D mode - grav_base = None - print('Need to specify Baseline survey for gravity in input file') - - for v, assim_time in enumerate(self.grav_config['vintage']): - time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ - dt.timedelta(days=assim_time) - - # porosity, saturation, densities, and fluid mass at individual time-steps - grav_struct[v] = self.calc_mass(time, v+1) # calculate the mass of each fluid in each grid cell - - - - vintage = [] - - - for v, assim_time in enumerate(self.grav_config['vintage']): - dg = self.calc_grav(grid, grav_base, grav_struct[v], pos) - vintage.append(deepcopy(dg)) - - #save_dic = {'grav': dg, **self.grav_config} - save_dic = { - 'grav': dg, 'P_vint': grav_struct[v]['PRESSURE'], 'rho_gas_vint':grav_struct[v]['GAS_DEN'], - 'meas_location': pos, **self.grav_config, - **{key: grav_struct[v][key] - grav_base[key] for key in grav_struct[v].keys()} - } - if save_folder is not None: - file_name = save_folder + os.sep + f"grav_vint{v}.npz" if save_folder[-1] != os.sep \ - else save_folder + f"grav_vint{v}.npz" - else: - file_name = folder + os.sep + f"grav_vint{v}.npz" if folder[-1] != os.sep \ - else folder + f"grav_vint{v}.npz" - prior_folder = 'Prior_ensemble_results' - try: - files = os.listdir(prior_folder) - filename_to_check = f"grav_vint{v}_{folder}.npz" - - if filename_to_check in files: - file_name_rec = 'Ensemble_results/' + f"grav_vint{v}_{folder}.npz" if folder[-1] != os.sep \ - else 'Ensemble_results/' + f"grav_vint{v}_{folder[:-1]}.npz" - else: - file_name_rec = 'Prior_ensemble_results/' + f"grav_vint{v}_{folder}.npz" if folder[-1] != os.sep \ - else 'Prior_ensemble_results/' + f"grav_vint{v}_{folder[:-1]}.npz" - - except: - file_name_rec = 'Ensemble_results/' + f"grav_vint{v}_{folder}.npz" if folder[-1] != os.sep \ - else 'Ensemble_results/' + f"grav_vint{v}_{folder[:-1]}.npz" - np.savez(file_name_rec, **save_dic) - - np.savez(file_name, **save_dic) - - - # 4D response - self.grav_result = [] - for i, elem in enumerate(vintage): - self.grav_result.append(elem) - - def calc_mass(self, time, time_index = None): - - if self.no_flow: - time_input = time_index - else: - time_input = time - - # fluid phases given as input - phases = str.upper(self.pem_input['phases']) - phases = phases.split() - # - grav_input = {} - tmp_dyn_var = {} - - - tmp = self._get_pem_input('RPORV', time_input) - grav_input['RPORV'] = np.array(tmp[~tmp.mask], dtype=float) - - tmp = self._get_pem_input('PRESSURE', time_input) - #if time_input == time_index and time_index > 0: # to be activiated in case on inverts for Delta Pressure - # # Inverts for changes in dynamic variables using time-lapse data - # tmp_baseline = self._get_pem_input('PRESSURE', 0) - # tmp = tmp + tmp_baseline - grav_input['PRESSURE'] = np.array(tmp[~tmp.mask], dtype=float) - # convert pressure from Bar to MPa - if 'press_conv' in self.pem_input and time_input == time: - grav_input['PRESSURE'] = grav_input['PRESSURE'] * self.pem_input['press_conv'] - #else: - # print('Keyword RPORV missing from simulation output, need updated pore volumes at each assimilation step') - # extract saturation - if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: # This should be extended - for var in phases: - if var in ['WAT', 'GAS']: - tmp = self._get_pem_input('S{}'.format(var), time_input) - #if time_input == time_index and time_index > 0: # to be activated in case on inverts for Delta S - # # Inverts for changes in dynamic variables using time-lapse data - # tmp_baseline = self._get_pem_input('S{}'.format(var), 0) - # tmp = tmp + tmp_baseline - #tmp = self.ecl_case.cell_data('S{}'.format(var), time) - grav_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) - grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] > 1] = 1 - grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] < 0] = 0 - - grav_input['SOIL'] = 1 - (grav_input['SWAT'] + grav_input['SGAS']) - grav_input['SOIL'][grav_input['SOIL'] > 1] = 1 - grav_input['SOIL'][grav_input['SOIL'] < 0] = 0 - - - tmp_dyn_var['SWAT'] = grav_input['SWAT'] # = {f'S{ph}': saturations[i] for i, ph in enumerate(phases)} - tmp_dyn_var['SGAS'] = grav_input['SGAS'] - tmp_dyn_var['SOIL'] = grav_input['SOIL'] - - - elif 'WAT' in phases and 'GAS' in phases: # Smeaheia model - for var in phases: - if var in ['GAS']: - tmp = self._get_pem_input('S{}'.format(var), time_input) - #if time_input == time_index and time_index > 0: # to be activated in case on inverts for Delta S - # Inverts for changes in dynamic variables using time-lapse data - # tmp_baseline = self._get_pem_input('S{}'.format(var), 0) - # tmp = tmp + tmp_baseline - #tmp = self.ecl_case.cell_data('S{}'.format(var), time) - grav_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) - grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] > 1] = 1 - grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] < 0] = 0 - - grav_input['SWAT'] = 1 - (grav_input['SGAS']) - - # fluid saturation - tmp_dyn_var['SWAT'] = grav_input['SWAT'] #= {f'S{ph}': saturations[i] for i, ph in enumerate(phases)} - tmp_dyn_var['SGAS'] = grav_input['SGAS'] - - elif 'OIL' in phases and 'GAS' in phases: # Original Smeaheia model - for var in phases: - if var in ['GAS']: - tmp = self._get_pem_input('S{}'.format(var), time_input) - #if time_input == time_index and time_index > 0: # to be activated in case on inverts for Delta S - # Inverts for changes in dynamic variables using time-lapse data - # tmp_baseline = self._get_pem_input('S{}'.format(var), 0) - # tmp = tmp + tmp_baseline - #tmp = self.ecl_case.cell_data('S{}'.format(var), time) - grav_input['S{}'.format(var)] = np.array(tmp[~tmp.mask], dtype=float) - grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] > 1] = 1 - grav_input['S{}'.format(var)][grav_input['S{}'.format(var)] < 0] = 0 - - grav_input['SOIL'] = 1 - (grav_input['SGAS']) - - # fluid saturation - tmp_dyn_var['SOIL'] = grav_input['SOIL'] #= {f'S{ph}': saturations[i] for i, ph in enumerate(phases)} - tmp_dyn_var['SGAS'] = grav_input['SGAS'] - - else: - print('Type and number of fluids are unspecified in calc_mass') - - - # fluid densities - for var in phases: - dens = var + '_DEN' - #tmp = self.ecl_case.cell_data(dens, time) - if self.no_flow: - if any('pressure' in key for key in self.state.keys()): - if 'press_conv' in self.pem_input: - conv2pa = 1e6 #MPa to Pa - else: - conv2pa = 1e5 # Bar to Pa - - if var == 'GAS': - if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: - tmp = PropsSI('D', 'T', 298.15, 'P', grav_input['PRESSURE']*conv2pa, 'Methane') - elif 'WAT' in phases and 'GAS' in grav_input['PRESSURE']: # Smeaheia model T = 37 C - tmp = PropsSI('D', 'T', 310.15, 'P', grav_input['PRESSURE']*conv2pa, 'CO2') - mask = np.zeros(tmp.shape, dtype=bool) - tmp = np.ma.array(data=tmp, dtype=tmp.dtype, mask=mask) - elif var == 'WAT': - tmp = PropsSI('D', 'T|liquid', 298.15, 'P', grav_input['PRESSURE']*conv2pa, 'Water') - mask = np.zeros(tmp.shape, dtype=bool) - tmp = np.ma.array(data=tmp, dtype=tmp.dtype, mask=mask) - else: - tmp = self._get_pem_input(dens, time_input) - else: - tmp = self._get_pem_input(dens, time_input) - grav_input[dens] = np.array(tmp[~tmp.mask], dtype=float) - tmp_dyn_var[dens] = grav_input[dens] - else: - tmp = self._get_pem_input(dens, time_input) - grav_input[dens] = np.array(tmp[~tmp.mask], dtype=float) - tmp_dyn_var[dens] = grav_input[dens] - - - tmp_dyn_var['PRESSURE'] = grav_input['PRESSURE'] - tmp_dyn_var['RPORV'] = grav_input['RPORV'] - self.dyn_var.extend([tmp_dyn_var]) - - #fluid masses - for var in phases: - mass = var + '_mass' - grav_input[mass] = grav_input[var + '_DEN'] * grav_input['S' + var] * grav_input['RPORV'] - - return grav_input - - def calc_grav(self, grid, grav_base, grav_repeat, pos): - - cell_centre = self.find_cell_centre(grid) - x = cell_centre[0] - y = cell_centre[1] - z = cell_centre[2] - - # Initialize dg as a zero array, with shape depending on the condition - # assumes the length of each vector gives the total number of measurement points - n_meas = (len(pos['x'])) - dg = np.zeros(n_meas) # 1D array for dg - dg[:] = np.nan - - # fluid phases given as input - phases = str.upper(self.pem_input['phases']) - phases = phases.split() - if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases: - dm = grav_repeat['OIL_mass'] + grav_repeat['WAT_mass'] + grav_repeat['GAS_mass'] - (grav_base['OIL_mass'] + grav_base['WAT_mass'] + grav_base['GAS_mass']) - - elif 'OIL' in phases and 'GAS' in phases: # Original Smeaheia model - dm = grav_repeat['OIL_mass'] + grav_repeat['GAS_mass'] - (grav_base['OIL_mass'] + grav_base['GAS_mass']) - # dm = grav_repeat['WAT_mass'] + grav_repeat['GAS_mass'] - (grav_base['WAT_mass'] + grav_base['GAS_mass']) - - elif 'WAT' in phases and 'GAS' in phases: # Smeaheia model - dm = grav_repeat['WAT_mass'] + grav_repeat['GAS_mass'] - (grav_base['WAT_mass'] + grav_base['GAS_mass']) - #dm = grav_repeat['WAT_mass'] + grav_repeat['GAS_mass'] - (grav_base['WAT_mass'] + grav_base['GAS_mass']) - - else: - dm = None - print('Type and number of fluids are unspecified in calc_grav') - - - for j in range(n_meas): - - # Calculate dg for the current measurement location (j, i) - dg_tmp = (z - pos['z'][j]) / ((x - pos['x'][j]) ** 2 + (y - pos['y'][j]) ** 2 + ( - z - pos['z'][j]) ** 2) ** (3 / 2) - - dg[j] = np.dot(dg_tmp, dm) - #print(f'Progress: {j + 1}/{n_meas}') # Mimicking wait bar - - # Scale dg by the constant - dg *= 6.67e-3 - - return dg - - def _get_grav_info(self, grav_config=None): - """ - GRAV configuration - """ - # list of configuration parameters in the "Grav" section of teh pipt file - config_para_list = ['baseline', 'vintage', 'water_depth', 'padding', 'grid_spacing', 'seabed'] - - if 'grav' in self.input_dict: - self.grav_config = {} - for elem in self.input_dict['grav']: - assert elem[0] in config_para_list, f'Property {elem[0]} not supported' - if elem[0] == 'vintage' and not isinstance(elem[1], list): - elem[1] = [elem[1]] - self.grav_config[elem[0]] = elem[1] - else: - self.grav_config = None - - def extract_data(self, member): - # start by getting the data from the flow simulator - super(flow_rock, self).extract_data(member) - - # get the gravity data from results - for prim_ind in self.l_prim: - # Loop over all keys in pred_data (all data types) - for key in self.all_data_types: - if 'grav' in key: - if self.true_prim[1][prim_ind] in self.grav_config['vintage']: - v = self.grav_config['vintage'].index(self.true_prim[1][prim_ind]) - self.pred_data[prim_ind][key] = self.grav_result[v].flatten() - -class flow_seafloor_disp(flow_rock, mixIn_multi_data): - def __init__(self, input_dict=None, filename=None, options=None, **kwargs): - super().__init__(input_dict, filename, options) - - assert 'sea_disp' in input_dict, 'To do subsidence/uplift simulation, please specify an "SEA_DISP" section in the pipt file' - self._get_disp_info() - - - def setup_fwd_run(self, **kwargs): - self.__dict__.update(kwargs) - - super().setup_fwd_run(redund_sim=None) - - def run_fwd_sim(self, state, member_i, del_folder=True): - # The inherited simulator also has a run_fwd_sim. Call this. - self.ensemble_member = member_i - self.pred_data = super().run_fwd_sim(state, member_i, del_folder) - - return self.pred_data - - def call_sim(self, folder=None, wait_for_proc=False, save_folder=None): - # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. - # Then, get the pem. - if folder is None: - folder = self.folder - - # run flow simulator - # success = True - success = super(flow_rock, self).call_sim(folder, True) - - # use output from flow simulator to forward model gravity response - if success: - # calculate gravity data based on flow simulation output - self.get_displacement_result(folder, save_folder) - - - return success - - def get_displacement_result(self, folder, save_folder): - if self.no_flow: - grid_file = self.pem_input['grid'] - grid = np.load(grid_file) - else: - self.ecl_case = ecl.EclipseCase(folder + os.sep + self.file + '.DATA') if folder[-1] != os.sep \ - else ecl.EclipseCase(folder + self.file + '.DATA') - grid = self.ecl_case.grid() - - self.dyn_var = [] - - # receiver locations - pad = self.disp_config.get('padding', 1500) # 3 km padding around the reservoir - if 'padding' not in self.disp_config: - print('Please specify extent of measurement locations, padding in input file, using 1.5 km as default') - dxy = self.disp_config.get('grid_spacing', 1500) # - if 'grid_spacing' not in self.disp_config: - print('Please specify grid spacing in input file, using 1.5 km as default') - if 'seabed' in self.disp_config and self.disp_config['seabed'] is not None: - file_path = self.disp_config['seabed'] - water_depth = self.get_seabed_depths(file_path) - else: - water_depth = self.disp_config.get('water_depth', 300) - if 'water_depth' not in self.disp_config: - print('Please specify water depths in input file, using 300 m as default') - pos = self.measurement_locations(grid, water_depth, pad, dxy) - - # loop over vintages with gravity acquisitions - disp_struct = {} - - if 'baseline' in self.disp_config: # 4D measurement - base_time = dt.datetime(self.startDate['year'], self.startDate['month'], - self.startDate['day']) + dt.timedelta(days=self.disp_config['baseline']) - # pore volume at time of baseline survey - disp_base = self.get_pore_volume(base_time, 0) - - else: - # seafloor displacement only work in 4D mode - disp_base = None - print('Need to specify Baseline survey for displacement modelling in input file') - - for v, assim_time in enumerate(self.disp_config['vintage']): - time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \ - dt.timedelta(days=assim_time) - - # pore volume and pressure at individual time-steps - disp_struct[v] = self.get_pore_volume(time, v+1) # calculate the mass of each fluid in each grid cell - - vintage = [] - - for v, assim_time in enumerate(self.disp_config['vintage']): - # calculate subsidence and uplift - dz_seafloor = self.map_z_response(disp_base, disp_struct[v], grid, pos) - vintage.append(deepcopy(dz_seafloor)) - - save_dic = {'sea_disp': dz_seafloor, 'meas_location': pos, **self.disp_config} - if save_folder is not None: - file_name = save_folder + os.sep + f"sea_disp_vint{v}.npz" if save_folder[-1] != os.sep \ - else save_folder + f"sea_disp_vint{v}.npz" - else: - file_name = folder + os.sep + f"sea_disp_vint{v}.npz" if folder[-1] != os.sep \ - else folder + f"sea_disp_vint{v}.npz" - file_name_rec = 'Ensemble_results/' + f"sea_disp_vint{v}_{folder}.npz" if folder[-1] != os.sep \ - else 'Ensemble_results/' + f"sea_disp_vint{v}_{folder[:-1]}.npz" - np.savez(file_name_rec, **save_dic) - - np.savez(file_name, **save_dic) - - - # 4D response - self.disp_result = [] - for i, elem in enumerate(vintage): - self.disp_result.append(elem) - - def get_pore_volume(self, time, time_index = None): - - if self.no_flow: - time_input = time_index - else: - time_input = time - - # - disp_input = {} - tmp_dyn_var = {} - - - tmp = self._get_pem_input('RPORV', time_input) - disp_input['RPORV'] = np.array(tmp[~tmp.mask], dtype=float) - - tmp = self._get_pem_input('PRESSURE', time_input) - #if time_input == time_index and time_index > 0: # to be activiated in case on inverts for Delta Pressure - # # Inverts for changes in dynamic variables using time-lapse data - # tmp_baseline = self._get_pem_input('PRESSURE', 0) - # tmp = tmp + tmp_baseline - disp_input['PRESSURE'] = np.array(tmp[~tmp.mask], dtype=float) - # convert pressure from Bar to MPa - if 'press_conv' in self.pem_input and time_input == time: - disp_input['PRESSURE'] = disp_input['PRESSURE'] * self.pem_input['press_conv'] - #else: - # print('Keyword RPORV missing from simulation output, need pdated porevolumes at each assimilation step') - - - tmp_dyn_var['PRESSURE'] = disp_input['PRESSURE'] - tmp_dyn_var['RPORV'] = disp_input['RPORV'] - self.dyn_var.extend([tmp_dyn_var]) - - return disp_input - - def compute_horizontal_distance(self, pos, x, y): - dx = pos['x'][:, np.newaxis] - x - dy = pos['y'][:, np.newaxis] - y - rho = np.sqrt(dx ** 2 + dy ** 2).flatten() - return rho - - def map_z_response(self, base, repeat, grid, pos): - """ - Maps out subsidence and uplift based either on the simulation - model pressure drop (method = 'pressure') or simulated change in pore volume - using either the van Opstal or Geertsma forward model - - Arguments: - base -- A dictionary containing baseline pressures and pore volumes. - repeat -- A dictionary containing pressures and pore volumes at repeat measurements. - - compute subsidence at position 'pos', b - - Output is modeled subsidence in cm. - - """ - - # Method to compute pore volume change - method = self.disp_config['method'].lower() - - # Forward model to compute subsidence/uplift response - model = self.disp_config['model'].lower() - - if self.disp_config['poisson'] > 0.5: - poisson = 0.5 - print('Poisson\'s ratio exceeds physical limits, setting it to 0.5') - else: - poisson = self.disp_config['poisson'] - - # Depth of rigid basement - z_base = self.disp_config['z_base'] - - compressibility = self.disp_config['compressibility'] # 1/MPa - - E = ((1 + poisson) * (1 - 2 * poisson)) / ((1 - poisson) * compressibility) - - # coordinates of cell centres - cell_centre = self.find_cell_centre(grid) - - # compute pore volume change between baseline and repeat survey - # based on the reservoir pore volumes in the individual vintages - if method == 'pressure': - dV = base['RPORV'] * (base['PRESSURE'] - repeat['PRESSURE']) * compressibility - else: - dV = base['RPORV'] - repeat['RPORV'] - - # coordinates of active cell centres - x = cell_centre[0] - y = cell_centre[1] - z = cell_centre[2] - - - if model == 'van_Opstal': - # Represents a signal change for subsidence/uplift. - #trans_func = t_van_opstal - component = ["Geertsma_vertical", "System_3_vertical"] - else: # Use Geertsma - #trans_func = t_geertsma - component = ["Geertsma_vertical"] - # Initialization - dz_1_2 = 0 - dz_3 = 0 - - # indices of active cells: - true_indices = np.where(grid['ACTNUM']) - # number of active gridcells - n_nucleous = len(true_indices[0]) - assert n_nucleous == len(x) - #pr = cProfile.Profile() - #pr.enable() - - for j in range(n_nucleous): - rho = self.compute_horizontal_distance(pos, x[j], y[j]) - THH, TRB = self.compute_deformation_transfer(pos['z'], z[j], z_base, rho, poisson, E, dV[j], component) - dz_1_2 = dz_1_2 + THH - dz_3 = dz_3 + TRB - - #pr.disable() - - # Print profiling results - #stats = pstats.Stats(pr) - #stats.strip_dirs() - #stats.sort_stats('cumulative') - #stats.print_stats() - - if model == 'van_Opstal': - # Represents a signal change for subsidence/uplift. - dz = dz_1_2 + dz_3 - else: # Use Geertsma - dz = dz_1_2 - - # Convert from meters to centimeters - dz *= 100 - - return dz - - def compute_van_opstal_transfer_function(self, z_res, z_base, rho, poisson): - """ - Compute the Van Opstal transfer function. - - Args: - z_res -- Numpy array of depths to reservoir cells [m]. - z_base -- Distance to the basement [m]. - rho -- Numpy array of horizontal distances in the field [m]. - poisson -- Poisson's ratio. - - Returns: - T -- Numpy array of the transfer function values. - T_geertsma -- Numpy array of the Geertsma transfer function values. - """ - - # Change to km scale - rho = rho / 1e3 - z_res = z_res / 1e3 - z_base = z_base / 1e3 - - - # Find lambda max (to optimize Hilbert transform) - cutoff = 1e-10 # Function value at max lambda - try: - lambda_max = fsolve(lambda x: 4 * (2 * x * z_base + 1) / (3 - 4 * poisson) * np.exp( - x * (np.max(z_res) - 2 * z_base)) - cutoff, 10)[0] - except: - lambda_max = 10 # Default value if unable to solve for max lambda - - lambda_vals = np.linspace(0, lambda_max, 100) - # range of lateral distances between measurement location and reservoir cells - nj = len(rho) - # range of vertical distances between measurement location and reservoir cells - ni = len(z_res) - # initialize - t_van_opstal = np.zeros((ni, nj)) - - # input function to make a hankel transform of order 0 of - c_t = self.van_opstal(lambda_vals, z_res[0], z_base, poisson) - - h_t, i_t = self.h_t(c_t, lambda_vals, rho) # Extract integrand - t_van_opstal[0, :] = (2 * z_res[0] / (rho ** 2 + z_res[0] ** 2) ** (3 / 2)) + h_t / (2 * np.pi) - - for i in range(1, ni): - C = self.van_opstal(lambda_vals, z_res[i], z_base, poisson) - h_t = self.h_t(C, lambda_vals, rho, i_t) - t_van_opstal[i, :] = (2 * z_res[i] / (rho ** 2 + z_res[i] ** 2) ** (3 / 2)) + h_t / (2 * np.pi) - - t_van_opstal *= 1e-6 # Convert back to meters - - t_geertsma = (2 * z_res[:, np.newaxis] / ((np.ones((ni, 1)) * rho) ** 2 + (z_res[:, np.newaxis]) ** 2) ** ( - 3 / 2)) * 1e-6 - - return t_van_opstal, t_geertsma - - def van_opstal(self, lambda_vals, z_res, z_base, poisson): - """ - Compute the Van Opstal transfer function. - - Args: - lambda_vals -- Numpy array of lambda values. - z_res -- Depth to reservoir [m]. - z_base -- Distance to the basement [m]. - poisson -- Poisson's ratio. - - Returns: - value -- Numpy array of computed values. - """ - - term1 = np.exp(lambda_vals * z_res) * (2 * lambda_vals * z_base + 1) - term2 = np.exp(-lambda_vals * z_res) * ( - 4 * lambda_vals ** 2 * z_base ** 2 + 2 * lambda_vals * z_base + (3 - 4 * poisson) ** 2) - - term3_numer = (3 - 4 * poisson) * ( - np.exp(-lambda_vals * (2 * z_base + z_res)) - np.exp(-lambda_vals * (2 * z_base - z_res))) - term3_denom = 2 * ((1 - 2 * poisson) ** 2 + lambda_vals ** 2 * z_base ** 2 + (3 - 4 * poisson) * np.cosh( - lambda_vals * z_base) ** 2) - - value = term1 - term2 - (term3_numer / term3_denom) - - return value - - def hankel_transform_order_0(self, f, r_max, num_points=1000): - """ - Computes the Hankel transform of order 0 of a function f(r). - - Parameters: - - f: callable, the function to transform, f(r) - - r_max: float, upper limit of the integral (approximate infinity) - - num_points: int, number of points for numerical integration - - Returns: - - k_values: array of k values - - H_k: array of Hankel transform evaluated at k_values - """ - r = np.linspace(0, r_max, num_points) - dr = r[1] - r[0] - f_r = f(r) - - def integrand(r, k): - return f(r) * j0(k * r) * r - - # Define a range of k values to evaluate - k_min, k_max = 0, 10 # adjust as needed - k_values = np.linspace(k_min, k_max, 100) - - H_k = [] - - for k in k_values: - # Perform numerical integration over r - result, _ = quad(integrand, 0, r_max, args=(k,)) - H_k.append(result) - - return k_values, np.array(H_k) - - def makeL(self, poisson, k, c, A_g, eps, lambda_): - L = A_g * ( - (4 * poisson - 3 + 2 * k * lambda_) * np.exp(-lambda_ * (k + c)) - - np.exp(lambda_ * eps * (k - c)) - ) - return L - - def makeM(self, poisson, k, c, A_g, eps, lambda_): - M = A_g * ( - (4 * poisson - 3 - 2 * k * lambda_) * np.exp(-lambda_ * (k + c)) - - eps * np.exp(lambda_ * eps * (k - c)) - ) - return M - - def makeDelta(self, poisson, k, lambda_): - Delta = ( - (4 * poisson - 3) * np.cosh(k * lambda_) ** 2 - - (k * lambda_) ** 2 - - (1 - 2 * poisson) ** 2 - ) - return Delta - - def makeB(self, poisson, k, c, A_g, eps, lambda_): - L = self.makeL(poisson, k, c, A_g, eps, lambda_) - M = self.makeM(poisson, k, c, A_g, eps, lambda_) - Delta = self.makeDelta(poisson, k, lambda_) - - numerator = ( - lambda_ * L * (2 * (1 - poisson) * np.cosh(k * lambda_) - lambda_ * k * np.sinh(k * lambda_)) - + lambda_ * M * ((1 - 2 * poisson) * np.sinh(k * lambda_) + k * lambda_ * np.cosh(k * lambda_)) - ) - - B = numerator / Delta - return B - - def makeC(self,poisson, k, c, A_g, eps, lambda_): - L = self.makeL(poisson, k, c, A_g, eps, lambda_) - M = self.makeM(poisson, k, c, A_g, eps, lambda_) - Delta = self.makeDelta(poisson, k, lambda_) - - numerator = ( - lambda_ * L * ((1 - 2 * poisson) * np.sinh(k * lambda_) - lambda_ * k * np.cosh(k * lambda_)) - + lambda_ * M * (2 * (1 - poisson) * np.cosh(k * lambda_) + k * lambda_ * np.sinh(k * lambda_)) - ) - - C = numerator / Delta - return C - - def uHH_integrand(self, lambda_, z, rho, eps, c, poisson): - val = lambda_ * (eps * np.exp(lambda_ * eps * (z - c)) + - (3 - 4 * poisson + 2 * z * lambda_) * - np.exp(-lambda_ * (z + c))) - return val * j0(lambda_ * rho) - - def compute_deformation_transfer(self, z, c, k, rho, poisson, E, dV, component): - scale = 1000 # convert to km - # depth of receiver positions - z = np.array(z)/scale - rho = np.array(rho)/scale - # depth of reservoir cell - c = c/scale - k = k/scale - component = list(component) - # number of measurement locations - n_rec = len(z) - assert len(rho) == n_rec - THH = np.zeros(n_rec) - TRB = np.zeros(n_rec) - - # Constants - A_g = -dV * E / (4 * np.pi * (1 + poisson)) - uHH_outside_intregral = -(A_g * (1 + poisson)) / E - uRB_outside_intregral = (1 + poisson) / E - - for c_n in component: - if c_n == 'Geertsma_vertical': - lambda_max = 15 / np.max(rho).item() - for i in range(n_rec): - if rho[i] > np.abs(c-z[i])*3: - THH[i] = 0 - else: - eps = np.sign(c - z[i]) - THH[i] = quad(lambda lambda_var: self.uHH_integrand(lambda_var, z[i], rho[i], eps, c, poisson), 0, lambda_max)[0] * uHH_outside_intregral - THH[i] = THH[i]*scale**-2 - - elif c_n == 'System_3_vertical': - lambda_max = 30 / np.max(rho).item() - num_points = 500 - lambda_grid = np.linspace(0, lambda_max, num_points) - - sinh_z = np.sinh(z[:, np.newaxis] * lambda_grid) - cosh_z = np.cosh(z[:, np.newaxis] * lambda_grid) - J0_rho = j0(lambda_grid * rho) - - # - for i in range(n_rec): - if rho[i] > np.abs(c-z[i])*3: - TRB[i] = 0 - else: - z_i = z[i] - sinh_z_i = sinh_z[i] - cosh_z_i = cosh_z[i] - - b_values = self.makeB(poisson, k, c, A_g, -1, lambda_grid) - c_values = self.makeC(poisson, k, c, A_g, -1, lambda_grid) - - part1 = b_values * (lambda_grid * z_i * cosh_z_i - (1 - 2 * poisson) * sinh_z_i) - part2 = c_values * ((2 * (1 - poisson) * cosh_z_i) - lambda_grid * z_i * sinh_z_i) - - values = (part1 + part2) * J0_rho[:, i]#J0_rho_j - - integral_result = np.trapz(values, lambda_grid) - TRB[i] = integral_result * uRB_outside_intregral - TRB[i] = TRB[i]*scale**-2 - - return THH, TRB - - def h_t(self, h, r=None, k=None, i_k=None): - """ - Hankel transform of order 0. - - Args: - h -- Signal h(r). - r -- Radial positions [m] (optional). - k -- Spatial frequencies [rad/m] (optional). - I -- Integration kernel (optional). - - Returns: - h_t -- Spectrum H(k). - I -- Integration kernel. - """ - - # Check if h is a vector - if h.ndim > 1: - raise ValueError('Signal must be a vector.') - - if r is None or len(r) == 0: - r = np.arange(len(h)) # Default to 0:numel(h)-1 - else: - r = np.sort(r) - h = h[np.argsort(r)] # Sort h according to sorted r - - if k is None or len(k) == 0: - k = np.pi / len(h) * np.arange(len(h)) # Default spatial frequencies - - if i_k is None: - # Create integration kernel I - r = np.concatenate([(r[:-1] + r[1:]) / 2, [r[-1]]]) # Midpoints plus last point - i_k = (2 * np.pi / k[:, np.newaxis]) * r * jv(1, k[:, np.newaxis] * r) # Bessel function - i_k[k == 0, :] = np.pi * r * r - i_k = i_k - np.hstack([np.zeros((len(k), 1)), i_k[:, :-1]]) # Shift integration kernel - else: - # Ensure I is sorted based on r - i_k = i_k[:, np.argsort(r)] - - # Compute Hankel Transform - h_t = np.reshape(i_k @ h.flatten(), k.shape) - - - - return h_t, i_k - - def _get_disp_info(self, disp_config=None): - """ - seafloor displacement (uplift/subsidence) configuration - """ - # list of configuration parameters in the "Grav" section of teh pipt file - config_para_list = ['baseline', 'vintage', 'method', 'model', 'poisson', 'compressibility', - 'z_base', 'grid_spacing', 'padding', 'seabed', 'water_depth'] - - if 'sea_disp' in self.input_dict: - self.disp_config = {} - for elem in self.input_dict['sea_disp']: - assert elem[0] in config_para_list, f'Property {elem[0]} not supported' - if elem[0] == 'vintage' and not isinstance(elem[1], list): - elem[1] = [elem[1]] - self.disp_config[elem[0]] = elem[1] - else: - self.disp_config = None - - def extract_data(self, member): - # start by getting the data from the flow simulator i.e. prod. and inj. data - super(flow_rock, self).extract_data(member) - - # get the gravity data from results - for prim_ind in self.l_prim: - # Loop over all keys in pred_data (all data types) - for key in self.all_data_types: - if 'sea_disp' in key: - if self.true_prim[1][prim_ind] in self.disp_config['vintage']: - v = self.disp_config['vintage'].index(self.true_prim[1][prim_ind]) - self.pred_data[prim_ind][key] = self.disp_result[v].flatten() - -class flow_grav_and_avo(flow_avo, flow_grav): - def __init__(self, input_dict=None, filename=None, options=None, **kwargs): - super().__init__(input_dict, filename, options) - - self.grav_input = {} - assert 'grav' in input_dict, 'To do GRAV simulation, please specify an "GRAV" section in the "FWDSIM" part' - self._get_grav_info() - - assert 'avo' in input_dict, 'To do AVO simulation, please specify an "AVO" section in the "FWDSIM" part' - self._get_avo_info() - - def setup_fwd_run(self, **kwargs): - self.__dict__.update(kwargs) - - super().setup_fwd_run(redund_sim=None) - - def run_fwd_sim(self, state, member_i, del_folder=True): - # The inherited simulator also has a run_fwd_sim. Call this. - self.ensemble_member = member_i - self.pred_data = super().run_fwd_sim(state, member_i, del_folder) - - return self.pred_data - - def call_sim(self, folder=None, wait_for_proc=False, save_folder=None): - # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. - # Then, get the pem. - if folder is None: - folder = self.folder - else: - self.folder = folder - - # run flow simulator - # success = True - success = super(flow_rock, self).call_sim(folder, True) - - # use output from flow simulator to forward model gravity response - if success: - # calculate gravity data based on flow simulation output - self.get_grav_result(folder, save_folder) - # calculate avo data based on flow simulation output - self.get_avo_result(folder, save_folder) - - 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] - #v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind]) - #self.pred_data[prim_ind][key] = self.avo_result[v].flatten() - -class flow_grav_seafloor_disp_and_avo(flow_avo, flow_grav, flow_seafloor_disp): - def __init__(self, input_dict=None, filename=None, options=None, **kwargs): - super().__init__(input_dict, filename, options) - - 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() - - assert 'sea_disp' in input_dict, 'To do subsidence/uplift simulation, please specify an "SEA_DISP" section in the "FWDSIM" part' - self._get_disp_info() - - def setup_fwd_run(self, **kwargs): - self.__dict__.update(kwargs) - - super().setup_fwd_run(redund_sim=None) - - def run_fwd_sim(self, state, member_i, del_folder=True): - # The inherited simulator also has a run_fwd_sim. Call this. - self.ensemble_member = member_i - self.pred_data = super().run_fwd_sim(state, member_i, del_folder) - - return self.pred_data - - def call_sim(self, folder=None, wait_for_proc=False, save_folder=None): - # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. - # Then, get the pem. - if folder is None: - folder = self.folder - else: - self.folder = folder - - # run flow simulator - # success = True - success = super(flow_rock, self).call_sim(folder, True) - - # use output from flow simulator to forward model gravity response - if success: - # calculate gravity data based on flow simulation output - self.get_grav_result(folder, save_folder) - # calculate avo data based on flow simulation output - self.get_avo_result(folder, save_folder) - # calculate gravity data based on flow simulation output - self.get_displacement_result(folder, save_folder) - - return success - - def extract_data(self, member): - # start by getting the data from the flow simulator i.e. prod. and inj. data - super(flow_rock, self).extract_data(member) - - # get the gravity data from results - for prim_ind in self.l_prim: - # Loop over all keys in pred_data (all data types) - for key in self.all_data_types: - if 'grav' in key: - if self.true_prim[1][prim_ind] in self.grav_config['vintage']: - v = self.grav_config['vintage'].index(self.true_prim[1][prim_ind]) - self.pred_data[prim_ind][key] = self.grav_result[v].flatten() - - if '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] - - if 'sea_disp' in key: - if self.true_prim[1][prim_ind] in self.disp_config['vintage']: - v = self.disp_config['vintage'].index(self.true_prim[1][prim_ind]) - self.pred_data[prim_ind][key] = self.disp_result[v].flatten() - -class flow_grav_seafloor_disp(flow_grav, flow_seafloor_disp): - def __init__(self, input_dict=None, filename=None, options=None, **kwargs): - super().__init__(input_dict, filename, options) - - assert 'grav' in input_dict, 'To do GRAV simulation, please specify an "GRAV" section in the "FWDSIM" part' - self._get_grav_info() - - assert 'sea_disp' in input_dict, 'To do subsidence/uplift simulation, please specify an "SEA_DISP" section in the "FWDSIM" part' - self._get_disp_info() - - def setup_fwd_run(self, **kwargs): - self.__dict__.update(kwargs) - - super().setup_fwd_run(redund_sim=None) - - def run_fwd_sim(self, state, member_i, del_folder=True): - # The inherited simulator also has a run_fwd_sim. Call this. - self.ensemble_member = member_i - self.pred_data = super().run_fwd_sim(state, member_i, del_folder) - - return self.pred_data - - def call_sim(self, folder=None, wait_for_proc=False, save_folder=None): - # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first. - # Then, get the pem. - if folder is None: - folder = self.folder - else: - self.folder = folder - - # run flow simulator - # success = True - success = super(flow_rock, self).call_sim(folder, True) - - # use output from flow simulator to forward model gravity response - if success: - # calculate gravity data based on flow simulation output - self.get_grav_result(folder, save_folder) - # calculate avo data based on flow simulation output - #self.get_avo_result(folder, save_folder) - # calculate gravity data based on flow simulation output - self.get_displacement_result(folder, save_folder) - - return success - - def extract_data(self, member): - # start by getting the data from the flow simulator i.e. prod. and inj. data - super(flow_rock, self).extract_data(member) - - # get the gravity data from results - for prim_ind in self.l_prim: - # Loop over all keys in pred_data (all data types) - for key in self.all_data_types: - if 'grav' in key: - if self.true_prim[1][prim_ind] in self.grav_config['vintage']: - v = self.grav_config['vintage'].index(self.true_prim[1][prim_ind]) - self.pred_data[prim_ind][key] = self.grav_result[v].flatten() - - if 'sea_disp' in key: - if self.true_prim[1][prim_ind] in self.disp_config['vintage']: - v = self.disp_config['vintage'].index(self.true_prim[1][prim_ind]) - self.pred_data[prim_ind][key] = self.disp_result[v].flatten() diff --git a/simulator/opm.py b/simulator/opm.py deleted file mode 100644 index 125e78f6..00000000 --- a/simulator/opm.py +++ /dev/null @@ -1,378 +0,0 @@ -"""Wrap OPM-flow""" -# External imports -from subprocess import call, DEVNULL, run -import os,sys -import shutil -import re -import time -from datetime import timedelta - -# Internal imports -from simulator.eclipse import eclipse -from misc.system_tools.environ_var import OPMRunEnvironment - - -class flow(eclipse): - """ - Class for running OPM flow with Eclipse input files. Inherits eclipse parent class for setting up and running - simulations, and reading the results. - """ - - def __init__(self,input_file=None,initialize_parent=True): - if initialize_parent: - super().__init__(input_file) - else: - self.file = input_file['filename'] - self.options = input_file - - def call_sim(self, folder=None, wait_for_proc=False): - """ - Call OPM flow simulator via shell. - - Parameters - ---------- - folder : str - Folder with runfiles. - - wait_for_proc : bool - Boolean determining if we wait for the process to be done or not. - - Changelog - --------- - - ST 18/10-18 - """ - # Filename - if folder is not None: - filename = folder + self.file - else: - filename = self.file - - success = True - #print(filename) - try: - with OPMRunEnvironment(filename, 'PRT', ['End of simulation', 'NOSIM']): - com = [] - if self.options['mpi']: - com.extend(self.options['mpi'].split()) - com.append(self.options['sim_path'] + 'flow') - if self.options['parsing-strictness']: - com.extend(['--parsing-strictness=' + self.options['parsing-strictness']]) - com.extend(['--output-dir=' + folder, * - self.options['sim_flag'].split(), filename + '.DATA']) - if 'sim_limit' in self.options: - call(com, stdout=DEVNULL, timeout=self.options['sim_limit']) - else: - call(com, stdout=DEVNULL) - raise ValueError # catch errors in run_sim - except Exception as e: - print('\nError in the OPM run.') # add rerun? - if not os.path.exists('Crashdump'): - shutil.copytree(folder, 'Crashdump') - success = False - - return success - - def check_sim_end(self, finished_member=None): - """ - Check in RPT file for "End of simulation" to see if OPM flow is done. - - Changelog - --------- - - ST 19/10-18 - """ - # Initialize output - # member = None - # - # # Search for output.dat file - # for file in os.listdir('En_' + str(finished_member)): # Search within a specific En_folder - # if file.endswith('PRT'): # look in PRT file - # with open('En_' + str(finished_member) + os.sep + file, 'r') as fid: - # for line in fid: - # if re.search('End of simulation', line): - # # TODO: not do time.sleep() - # # time.sleep(0.1) - # member = finished_member - - return finished_member - - @staticmethod - def SLURM_HPC_run(n_e, venv,filename, **kwargs): - """ - HPC run manager for SLURM. - - This function will start num_runs of sim.call_sim() using job arrays in SLURM. - """ - # Extract the filename from the kwargs - filename_str = f'"{filename.upper()}"' if filename is not None else "" - - # Extract mpi flag from kwargs - mpi = kwargs.get("mpi", None) - mpi_str = f'"{mpi}"' if mpi is not None else "mpirun --bind-to none -np 1" - - # set number of tasks to the number following -np in mpi_str (default is 1) - n_tasks = re.search(r"-np (\d+)", mpi_str).group(1) - - - # extract the sim_limit from kwargs. Default is 1 hour - sim_limit = kwargs.get("sim_limit", None) - sim_limit_str = f'--time={str(timedelta(seconds=sim_limit))}' if sim_limit is not None else "--time=01:00:00" - - diff_ne = n_e[-1] - n_e[0] - - slurm_script = f"""#!/bin/bash -#SBATCH --partition=comp -#SBATCH --job-name=EnDA -#SBATCH --array=0-{diff_ne} -#SBATCH {sim_limit_str} -#SBATCH --mem=4G -#SBATCH --ntasks={n_tasks} -#SBATCH --cpus-per-task=2 -#SBATCH --export=ALL -#SBATCH --output=/dev/null - -# OPTIONAL: load modules here -module load Python -export LMOD_DISABLE_SAME_NAME_AUTOSWAP=no -module load opm-simulators - -source {venv} - -# Set folder based on SLURM_ARRAY_TASK_ID -folder="En_$(( {n_e[0]} + SLURM_ARRAY_TASK_ID ))/" - -python -m simulator.opm "$folder" {filename_str} {mpi_str} -""" - script_name = "submit_test_parallel_mpi.sh" - with open(script_name, "w") as f: - f.write(slurm_script) - - # Make it executable (optional): - os.chmod(script_name, 0o755) - - # print(f"Created SLURM script: {script_name}") - # print(f"Submitting array job with {num_runs} tasks...") - - # Submit the script to SLURM - cmd = ["sbatch", script_name] - result = run(cmd, capture_output=True, text=True) - - # remove script file - os.remove(script_name) - - # Extract the job ID from the output - match = re.search(r"Submitted batch job (\d+)", result.stdout) - if match: - return match.group(1) # Return the main job ID - else: - print("Failed to extract Job ID from sbatch output.") - return None - - @staticmethod - def SLURM_ARRAY_HPC_run(n_e, venv, filename, **kwargs): - """ - HPC run manager for Slurm array jobs. - Each ensemble member runs independently in its own task. - - Parameters - ---------- - n_e : list[int] - Indices of ensemble members to simulate. - venv : str - Path to Python virtual environment activate script. - filename : str - Simulation input file. - kwargs : dict - Extra simulation options. Recognized: - - sim_limit (float seconds or str HH:MM:SS) - - mem (default "4G") - - cpus_per_task (default 2) - """ - - # Start and end indices for the array - start_idx = n_e[0] - end_idx = n_e[-1] - num_tasks = end_idx - start_idx - - # Extract options - sim_limit = kwargs.get("sim_limit", None) - if sim_limit is not None: - if isinstance(sim_limit, (int, float)): - from datetime import timedelta - sim_limit_str = f'--time={str(timedelta(seconds=sim_limit))}' - else: - sim_limit_str = f'--time={sim_limit}' - else: - sim_limit_str = "--time=02:00:00" - - mem = kwargs.get("mem", "4G") - cpus_per_task = kwargs.get("cpus_per_task", 1) - - slurm_script = f"""#!/bin/bash -#SBATCH --job-name=EnDA_array -#SBATCH --partition=comp -#SBATCH --array=0-{num_tasks} -#SBATCH --cpus-per-task={cpus_per_task} -#SBATCH --mem={mem} -#SBATCH {sim_limit_str} -#SBATCH --output=logs/job_%A_%a.out -#SBATCH --error=logs/job_%A_%a.err - -module load Python -export LMOD_DISABLE_SAME_NAME_AUTOSWAP=no -module load opm-simulators -source {venv} - -IDX=$(( {start_idx} + SLURM_ARRAY_TASK_ID )) -FOLDER="En_$IDX/" - -python -m simulator.opm "$FOLDER" {filename} "" -""" - - script_name = "submit_array.sh" - with open(script_name, "w") as f: - f.write(slurm_script) - - os.chmod(script_name, 0o755) - - result = run(["sbatch", script_name], capture_output=True, text=True) - - os.remove(script_name) - - match = re.search(r"Submitted batch job (\d+)", result.stdout) - if match: - return match.group(1) - else: - print("Job submission failed:", result.stderr) - return None - - - def are_jobs_done(self, job_id): - """Check if all job array tasks are completed using sacct.""" - check_cmd = ["sacct", "-j", f"{job_id}", "--format=JobID,State", "--noheader"] - - # print(check_cmd) - - check_result = run(check_cmd, capture_output=True, text=True) - - while not len(check_result.stdout): # if spinning up - time.sleep(1) - check_result = run(check_cmd, capture_output=True, text=True) - - # print(check_result.stdout) - - return_states = [] - - job_states = check_result.stdout.strip().split("\n") - for job in job_states: - parts = job.split() - if len(parts) >= 2: - state = parts[1] - if state not in ["COMPLETED", "FAILED", "CANCELLED"]: - return False # A job is still running or pendin - else: - if state == "FAILED" or state == "CANCELLED": - return_states.append(False) - else: - return_states.append(True) - - return return_states - - def wait_for_jobs(self,job_id,wait_time=10): - """Wait until all job array tasks are completed.""" - #print(f"Waiting for job array {job_id} to complete...") - - val = self.are_jobs_done(job_id) - while not val: - time.sleep(wait_time) # Wait for 10 seconds before checking again - val = self.are_jobs_done(job_id) - - return val - - #print(f"All jobs in array {job_id} are completed.") - - - -class ebos(eclipse): - """ - Class for running OPM ebos with Eclipse input files. Inherits eclipse parent class for setting up and running - simulations, and reading the results. - """ - - def call_sim(self, folder=None, wait_for_proc=False): - """ - Call OPM flow simulator via shell. - - Parameters - ---------- - folder : str - Folder with runfiles. - - wait_for_proc : bool - Determines whether to wait for the process to be done or not. - - Changelog - --------- - - RJL 27/08-19 - """ - # Filename - if folder is not None: - filename = folder + self.file - else: - filename = self.file - - # Run simulator - if 'sim_path' not in self.options.keys(): - self.options['sim_path'] = '' - - with OPMRunEnvironment(filename, 'OUT', 'Timing receipt'): - with open(filename+'.OUT', 'w') as f: - call([self.options['sim_path'] + 'ebos', '--output-dir=' + folder, - *self.options['sim_flag'].split(), filename + '.DATA'], stdout=f) - - def check_sim_end(self, finished_member=None): - """ - Check in RPT file for "End of simulation" to see if OPM ebos is done. - - Changelog - --------- - - RJL 27/08-19 - """ - # Initialize output - # member = None - # - # # Search for output.dat file - # for file in os.listdir('En_' + str(finished_member)): # Search within a specific En_folder - # if file.endswith('OUT'): # look in OUT file - # with open('En_' + str(finished_member) + os.sep + file, 'r') as fid: - # for line in fid: - # if re.search('Timing receipt', line): - # # TODO: not do time.sleep() - # # time.sleep(0.1) - # member = finished_member - - return finished_member - - -if __name__ == "__main__": - import sys - if len(sys.argv) != 4: - print("Usage: python -m simulator.opm ") - sys.exit(1) - - folder = sys.argv[1] - filename = sys.argv[2] - mpi = sys.argv[3] - options = {} - options['sim_path'] = '' - options['sim_flag'] = '' - options['mpi'] = mpi - options['parsing-strictness'] = '' - options['filename'] = filename - - sim = flow(input_file=options,initialize_parent=False) - success = sim.call_sim(folder=folder) - if success: - sys.exit(0) - else: - sys.exit(1) # ensure that slurm catch the error - #print("Success!" if success else "Failed.") diff --git a/simulator/rockphysics/__init__.py b/simulator/rockphysics/__init__.py deleted file mode 100644 index 22dde210..00000000 --- a/simulator/rockphysics/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""Compute elastic properties.""" diff --git a/simulator/rockphysics/softsandrp.py b/simulator/rockphysics/softsandrp.py deleted file mode 100644 index dac1a5a6..00000000 --- a/simulator/rockphysics/softsandrp.py +++ /dev/null @@ -1,881 +0,0 @@ -"""Descriptive description.""" - -__author__ = {'TM', 'TB', 'ML'} - -# standardrp.py -import numpy as np -import sys -import multiprocessing as mp -from CoolProp.CoolProp import PropsSI # http://coolprop.org/#high-level-interface-example -import CoolProp.CoolProp as CP -# Density of carbon dioxide at 100 bar and 25C # Smeaheia 37 degrees C -#rho_co2 = PropsSI('D', 'T', 298.15, 'P', 100e5, 'CO2') -from numpy.random import poisson - -# internal load -from misc.system_tools.environ_var import OpenBlasSingleThread # Single threaded OpenBLAS runs - - -class elasticproperties: - """ - Calculate elastic properties from standard - rock-physics models, specifically following Batzle - and Wang, Geophysics, 1992, for fluid properties, and - Report 1 in Abul Fahimuddin's thesis at Universty of - Bergen (2010) for other properties. - - Example - ------- - >>> porosity = 0.2 - ... pressure = 5 - ... phases = ["Oil","Water"] - ... saturations = [0.3, 0.5] - ... - ... satrock = Elasticproperties() - ... satrock.calc_props(phases, saturations, pressure, porosity) - """ - - def __init__(self, input_dict): - self.dens = None - self.bulkmod = None - self.shearmod = None - self.bulkvel = None - self.shearvel = None - self.bulkimp = None - self.shearimp = None - # The overburden for each grid cell must be - # specified as values on an .npz-file whose - # name is given in input_dict. - self.input_dict = input_dict - self._extInfoInputDict() - - def _extInfoInputDict(self): - # The key word for the file name in the - # dictionary must read "overburden" - if 'overburden' in self.input_dict: - obfile = self.input_dict['overburden'] - npzfile = np.load(obfile) - # The values of overburden must have been - # stored on file using: - # np.savez(, - # obvalues=) - self.overburden = npzfile['obvalues'] - npzfile.close() - #else: - # # Norne litho pressure equation in Bar - # P_litho = -49.6 + 0.2027 * Z + 6.127e-6 * Z ** 2 # Using e-6 for scientific notation - # # Convert reservoir pore pressure from Bar to MPa - # P_litho *= 0.1 - # self.overburden = P_litho - - if 'baseline' in self.input_dict: - self.baseline = self.input_dict['baseline'] # 4D baseline - if 'parallel' in self.input_dict: - self.parallel = self.input_dict['parallel'] - - def _filter(self): - bulkmod = self.bulkimp - self.bulkimp = bulkmod.flatten() - - def setup_fwd_run(self, state): - """ - Setup the input parameters to be used in the PEM simulator. Parameters can be an ensemble or a single array. - State is set as an attribute of the simulator, and the correct value is determined in self.pem.calc_props() - - Parameters - ---------- - state : dict - Dictionary of input parameters or states. - - Changelog - --------- - - KF 11/12-2018 - """ - # self.inv_state = {} - # list_pem_param =[el for el in [foo for foo in self.pem['garn'].keys()] + [foo for foo in self.filter.keys()] + - # [foo for foo in self.__dict__.keys()]] - - # list_tot_param = state.keys() - # for param in list_tot_param: - # if param in list_pem_param or (param.split('_')[-1] in ['garn', 'rest']): - # self.inv_state[param] = state[param] - - pass - - def calc_props(self, phases, saturations, pressure, - porosity, dens = None, wait_for_proc=None, ntg=None, Rs=None, press_init=None, ensembleMember=None): - ### - - # - if not isinstance(phases, list): - phases = [phases] - if not isinstance(saturations, list): - saturations = [saturations] - if not isinstance(pressure, list) and \ - type(pressure).__module__ != 'numpy': - pressure = [pressure] - if not isinstance(porosity, list) and \ - type(porosity).__module__ != 'numpy': - porosity = [porosity] - # - # Load "overburden" pressures into local variable to - # comply with remaining code parts - poverburden = self.overburden - - # debug - self.pressure = pressure - self.peff = poverburden - pressure - self.porosity = porosity - - if press_init is None: - p_init = self.p_init - else: - p_init = press_init - - # Average number of contacts that each grain has with surrounding grains - coordnumber = self._coordination_number() - - # porosity value separating the porous media's mechanical and acoustic behaviour - phicritical = self._critical_porosity() - - - # Check that no. of phases is equal to no. of - # entries in saturations list - # - assert (len(saturations) == len(phases)) - # - # Make saturation a Numpy array (so that we - # can easily access the values for each - # phase at one grid cell) - # - # Transpose makes it a no. grid cells x phases - # array - saturations = np.array(saturations).T - # - # Check if we actually inputted saturation values - # for a single grid cell. If yes, we redefine - # saturations to get it on the correct form (no. - # grid cells x phases array). - # - if saturations.ndim == 1: - saturations = \ - np.array([[x] for x in saturations]).T - # - # Loop over all grid cells and calculate the - # various saturated properties - # - self.phases = phases - - self.dens = np.zeros(len(saturations[:, 0])) - self.bulkmod = np.zeros(len(saturations[:, 0])) - self.shearmod = np.zeros(len(saturations[:, 0])) - self.bulkvel = np.zeros(len(saturations[:, 0])) - self.shearvel = np.zeros(len(saturations[:, 0])) - self.bulkimp = np.zeros(len(saturations[:, 0])) - self.shearimp = np.zeros(len(saturations[:, 0])) - - if ntg is None: - ntg = [None for _ in range(len(saturations[:, 0]))] - if Rs is None: - Rs = [None for _ in range(len(saturations[:, 0]))] - if p_init is None: - p_init = [None for _ in range(len(saturations[:, 0]))] - - - if dens is not None: - assert (len(dens) == len(phases)) - # Transpose makes it a no. grid cells x phases array - dens = np.array(dens).T - - # - denss, bulks, shears = self._solidprops_Johansen() - - for i in range(len(saturations[:, 0])): - # - # Calculate fluid properties - # - if dens is None: - densf_SI = self._fluid_densSIprop(self.phases, - saturations[i, :], pressure[i]) - bulkf_Brie = self._fluidprops_Brie(self.phases, saturations[i, :], pressure[i], densf_SI) - densf, bulkf = \ - self._fluidprops_Wood(self.phases, - saturations[i, :], pressure[i], Rs[i]) - else: - densf = self._fluid_dens(saturations[i, :], dens[i, :]) - - bulkf = self._fluidprops_Brie(self.phases, saturations[i, :], pressure[i], densf) - # - #denss, bulks, shears = \ - # self._solidprops(porosity[i], ntg[i], i) - - # - # Calculate dry rock moduli - # - - #bulkd, sheard = \ - # self._dryrockmoduli(porosity[i], - # overburden[i], - # pressure[i], bulks, - # shears, i, ntg[i], p_init[i], denss, Rs[i], self.phases) - # - peff = self._effective_pressure(poverburden[i], pressure[i]) - - - bulkd, sheard = \ - self._dryrockmoduli_Smeaheia(coordnumber, phicritical, porosity[i], peff, bulks, shears) - - # ------------------------------- - # Calculate saturated properties - # ------------------------------- - # - # Density (kg/m3) - # - self.dens[i] = (porosity[i]*densf + - (1-porosity[i])*denss) - # - # Moduli (MPa) - # - self.bulkmod[i] = \ - bulkd + (1 - bulkd/bulks)**2 / \ - (porosity[i]/bulkf + - (1-porosity[i])/bulks - - bulkd/(bulks**2)) - self.shearmod[i] = sheard - # - # Velocities (km/s) - # - self.bulkvel[i] = \ - np.sqrt((abs(self.bulkmod[i]) + - 4*self.shearmod[i]/3)/(self.dens[i])) - self.shearvel[i] = \ - np.sqrt(self.shearmod[i] / - (self.dens[i])) - # - # convert from (km/s) to (m/s) - # - self.bulkvel[i] *= 1000 - self.shearvel[i] *= 1000 - # - # Impedance (m/s)*(kg/m3) - # - self.bulkimp[i] = self.dens[i] * \ - self.bulkvel[i] - self.shearimp[i] = self.dens[i] * \ - self.shearvel[i] - - - - def getMatchProp(self, petElProp): - if petElProp.lower() == 'density': - self.match_prop = self.getDens() - elif petElProp.lower() == 'bulk_modulus': - self.match_prop = self.getBulkMod() - elif petElProp.lower() == 'shear_modulus': - self.match_prop = self.getShearMod() - elif petElProp.lower() == 'bulk_velocity': - self.match_prop = self.getBulkVel() - elif petElProp.lower() == 'shear_velocity': - self.match_prop = self.getShearVel() - elif petElProp.lower() == "bulk_impedance": - self.match_prop = self.getBulkImp() - elif petElProp.lower() == 'shear_impedance': - self.match_prop = self.getShearImp() - else: - print("\nError in getMatchProp method") - print("No model output type selected for " - "data match.") - print("Legal model output types are " - "(case insensitive):") - print("Density, bulk modulus, shear " - "modulus, bulk velocity,") - print("shear velocity, bulk impedance, " - "shear impedance") - sys.exit(1) - return self.match_prop - - def getDens(self): - return self.dens - - def getBulkMod(self): - return self.bulkmod - - def getShearMod(self): - return self.shearmod - - def getBulkVel(self): - return self.bulkvel - - def getShearVel(self): - return self.shearvel - - def getBulkImp(self): - return self.bulkimp - - def getShearImp(self): - return self.shearimp - - def getOverburdenP(self): - return self.overburden - - def getPressure(self): - return self.pressure - - def getPeff(self): - return self.peff - - def getPorosity(self): - return self.porosity - # - # =================================================== - # Fluid properties start - # =================================================== - # - def _fluid_densSIprop(self, phases, fsats, press, t= 37, CO2 = None): - - conv2Pa = 1e6 # MPa to Pa - ta = t + 273.15 # absolute temp in K - # fluid densities - fdens = 0.0 - - for i in range(len(phases)): - # - # Calculate mixture properties by summing - # over individual phase properties - # - - var = phases[i] - if var == 'GAS' and CO2 is None: - pdens = PropsSI('D', 'T', ta, 'P', press * conv2Pa, 'Methane') - elif var == 'GAS' and CO2 is True: - pdens = PropsSI('D', 'T', ta, 'P', press * conv2Pa, 'CO2') - elif var == 'OIL': - CP.get_global_param_string('predefined_mixtures').split(',')[0:6] - #pdens = CP.PropsSI('D', 'T', ta, 'P', press * conv2Pa, 'Ekofisk.mix') - pdens = CP.PropsSI('D', 'T', ta, 'P', press * conv2Pa, 'butane') - elif var == 'WAT': - pdens = PropsSI('D', 'T|liquid', ta, 'P', press * conv2Pa, 'Water') - - fdens = fdens + fsats[i] * abs(pdens) - - return fdens - - def _fluidprops_Wood(self, fphases, fsats, fpress, Rs=None): - # - # Calculate fluid density and bulk modulus - # - # - # Input - # fphases - fluid phases present; Oil - # and/or Water and/or Gas - # fsats - fluid saturation values for - # fluid phases in "fphases" - # fpress - fluid pressure value (MPa) - # Rs - Gas oil ratio. Default value None - - # - # Output - # fdens - density of fluid mixture for - # pressure value "fpress" inherited - # from phaseprops) - # fbulk - bulk modulus of fluid mixture for - # pressure value "fpress" (unit - # inherited from phaseprops) - # - # ----------------------------------------------- - # - fdens = 0.0 - fbinv = 0.0 - - for i in range(len(fphases)): - # - # Calculate mixture properties by summing - # over individual phase properties - # - pdens, pbulk = self._phaseprops(fphases[i], - fpress, Rs) - fdens = fdens + fsats[i]*abs(pdens) - fbinv = fbinv + fsats[i]/abs(pbulk) - fbulk = 1.0/fbinv - # - return fdens, fbulk - # - # --------------------------------------------------- - # - - def _fluid_dens(self, fsatsp, fdensp): - fdens = sum(fsatsp * fdensp) - return fdens - - def _fluidprops_Brie(self, fphases, fsats, fpress, fdens, Rs=None, e = 5): - # - # Calculate fluid density and bulk modulus BRIE et al. 1995 - # Assumes two phases liquid and gas - # - # Input - # fphases - fluid phases present; Oil - # and/or Water and/or Gas - # fsats - fluid saturation values for - # fluid phases in "fphases" - # fdens - fluid density for given pressure and temperature - # fpress - fluid pressure value (MPa) - # Rs - Gas oil ratio. Default value None - # e - Brie's exponent (e= 5 Utsira sand filled with brine and CO2 - # Figure 7 in Carcione et al. 2006 "Physics and Seismic Modeling - # for Monitoring CO 2 Storage" - - # - # Output - # fbulk - bulk modulus of fluid mixture for - # pressure value "fpress" (unit - # inherited from phaseprops) - # - # ----------------------------------------------- - # - - - for i in range(len(fphases)): - # - if fphases[i].lower() in ["oil", "wat"]: - fsatsl = fsats[i] - pbulkl = self._phaseprops_Smeaheia(fphases[i], fpress, fdens, Rs) - elif fphases[i].lower() in ["gas"]: - pbulkg = self._phaseprops_Smeaheia(fphases[i], fpress, fdens, Rs) - - - fbulk = (pbulkl - pbulkg) * (fsatsl)**e + pbulkg - - # - return fbulk - # - # --------------------------------------------------- - # - @staticmethod - def pseudo_p_t(pres, t, gs): - """Calculate the pseudoreduced temperature and pressure according to Thomas et al. 1970. - - Parameters - ---------- - pres : float or array-like - Pressure in MPa - t : float or array-like - Temperature in °C - gs : float - Gas gravity - - Returns - ------- - float or array-like - Ta: absolute temperature - Ppr:pseudoreduced pressure - Tpr:pseudoreduced temperature - """ - - # convert the temperature to absolute temperature - ta = t + 273.15 - p_pr = pres / (4.892 - 0.4048 * gs) - t_pr = ta / (94.72 + 170.75 * gs) - return ta, p_pr, t_pr - # - # --------------------------------------------------- - # - @staticmethod - def dz_dp(p_pr, t_pr): - """Values for dZ/dPpr obtained from equation 10b in Batzle and Wang (1992). - """ - # analytic - dz_dp = (0.03 + 0.00527 * (3.5 - t_pr) ** 3) + 0.109 * (3.85 - t_pr) ** 2 * 1.2 * p_pr ** 0.2 * -( - 0.45 + 8 * (0.56 - 1 / t_pr) ** 2) / t_pr * np.exp( - -(0.45 + 8 * (0.56 - 1 / t_pr) ** 2) * p_pr ** 1.2 / t_pr) - - # numerical approximation - # dzdp= 1.938783*P_pr**0.2*(1 - 0.25974025974026*T_pr)**2*(-8*(0.56 - 1/T_pr)**2 - 0.45)* - # np.exp(P_pr**1.2*(-8*(0.56 - 1/T_pr)**2 - 0.45)/T_pr)/T_pr + 0.22595125*(1 - 0.285714285714286*T_pr)**3 - # + 0.03 - return dz_dp - # - #----------------------------------------------------------- - # - def _phaseprops_Smeaheia(self, fphase, press, fdens, Rs=None, t = 37, CO2 = True): - # - # Calculate properties for a single fluid phase - # - # - # Input - # fphase - fluid phase; Oil, Water or Gas - # press - fluid pressure value (MPa) - # fdens - fluid density (kg/m3) - # t - temperature in degrees C - # - # Output - # pbulk - bulk modulus of fluid phase - # "fphase" for pressure value - # "press" (MPa) - # - # ----------------------------------------------- - # References - # ---------- - # Xu, H. (2006). Calculation of CO2 acoustic properties using Batzle-Wang equations. Geophysics, 71(2), F21-F23. - # """ - - if fphase.lower() == "wat": # refers to pure water or brine - #Compute the bulk modulus of pure water as a function of temperature and pressure - #using Batzle and Wang (1992). - if np.any(press > 100): - print('pressures above about 100 MPa-> inaccurate estimations of water velocity') - w = np.array([[1.40285e+03, 1.52400e+00, 3.43700e-03, -1.19700e-05], - [4.87100e+00, -1.11000e-02, 1.73900e-04, -1.62800e-06], - [-4.78300e-02, 2.74700e-04, -2.13500e-06, 1.23700e-08], - [1.48700e-04, -6.50300e-07, -1.45500e-08, 1.32700e-10], - [-2.19700e-07, 7.98700e-10, 5.23000e-11, -4.61400e-13]]) - v_w = sum(w[i, j] * t ** i * press ** j for i in range(5) for j in range(4)) # m/s - K_w = fdens * v_w ** 2 * 1e-6 - if CO2 is True: # refers to brine - salinity = 35000 / 1000000 - s1 = 1170 - 9.6 * t + 0.055 * t ** 2 - 8.5e-5 * t ** 3 + 2.6 * press - 0.0029 * t * press - 0.0476 * press ** 2 - s15 = 780 - 10 * press + 0.16 * press ** 2 - s2 = -820 - v_b = v_w + s1 * salinity + s15 * salinity ** 1.5 + s2 * salinity ** 2 - x = 300 * press - 2400 * press * salinity + t * (80 + 3 * t - 3300 * salinity - 13 * press + 47 * press * salinity) - rho_b = fdens + salinity * (0.668 + 0.44 * salinity + 1e-6 * x) - pbulk = rho_b * v_b ** 2 * 1e-6 - else: - pbulk = K_w - - elif fphase.lower() == "gas" and CO2 is True: # refers to CO2 - R = 8.3145 # J.mol-1K-1 gas constant for CO2 - gs = 1.5189 # Specific gravity #https://www.engineeringtoolbox.com/specific-gravities-gases-d_334.html - ta, p_pr, t_pr = self.pseudo_p_t(press, t, gs) - - E = 0.109 * (3.85 - t_pr) ** 2 * np.exp(-(0.45 + 8 * (0.56 - 1 / t_pr) ** 2) * p_pr ** 1.2 / t_pr) - Z = (0.03 + 0.00527 * (3.5 - t_pr) ** 3) * p_pr + (0.642 * t_pr - 0.007 * t_pr ** 4 - 0.52) + E - rho = 28.8 * gs * press / (Z * R * ta) # g/cm3 - - r_0 = 0.85 + 5.6 / (p_pr + 2) + 27.1 / (p_pr + 3.5) ** 2 - 8.7 * np.exp(-0.65 * (p_pr + 1)) - dz_dp = self.dz_dp(p_pr, t_pr) - pbulk = press / (1 - p_pr * dz_dp / Z) * r_0 - - #pbulk_test = self.test_new_implementation(press) - #print(np.max(pbulk-pbulk_test)) - - elif fphase.lower() == "gas": # refers to Methane - gs = 0.5537 #https://www.engineeringtoolbox.com/specific-gravities-gases-d_334.html - R = 8.3145 # J.mol-1K-1 gas constant - ta, p_pr, t_pr = self.pseudo_p_t(press, t, gs) - E = 0.109 * (3.85 - t_pr) ** 2 * np.exp(-(0.45 + 8 * (0.56 - 1 / t_pr) ** 2) * p_pr ** 1.2 / t_pr) - Z = (0.03 + 0.00527 * (3.5 - t_pr) ** 3) * p_pr + (0.642 * t_pr - 0.007 * t_pr ** 4 - 0.52) + E - rho = 28.8 * gs * press / (Z * R * ta) # g/cm3 - - r_0 = 0.85 + 5.6 / (p_pr + 2) + 27.1 / (p_pr + 3.5) ** 2 - 8.7 * np.exp(-0.65 * (p_pr + 1)) - dz_dp = self.dz_dp(p_pr, t_pr) - pbulk = press / (1 - p_pr * dz_dp / Z) * r_0 - - elif fphase.lower() == "oil": #pure oil - # Estimate the oil bulk modulus at specific temperature and pressure. - v = 2096 * (fdens / (2600 - fdens)) ** 0.5 - 3.7 * t + 4.64 * press + 0.0115 * ( - 4.12 * (1080 / fdens - 1) ** 0.5 - 1) * t * press - pbulk = fdens * v ** 2 - - - # - return pbulk - - # - def test_new_implementation(self, press): - # Values from .DATA file for Smeaheia (converted to MPa) - press_range = np.array( - [0.101, 0.885, 1.669, 2.453, 3.238, 4.022, 4.806, 5.590, 6.2098, 7.0899, 7.6765, 8.2630, 8.8495, 9.4359, - 10.0222, 10.6084, 11.1945, 14.7087, 17.6334, 20.856, 23.4695, 27.5419]) # Example pressures in MPa - Bo_values = np.array( - [1.07365, 0.11758, 0.05962, 0.03863, 0.02773, 0.02100, 0.01639, 0.01298, 0.010286, 0.007578, 0.005521, - 0.003314, 0.003034, 0.002919, 0.002851, 0.002802, 0.002766, 0.002648, 0.002599, 0.002566, 0.002546, - 0.002525]) # Example formation volume factors in m^3/kg - - # Calculate numerical derivative of Bo with respect to Pressure - dBo_dP = - np.gradient(Bo_values, press_range) - # Calculate isothermal compressibility (van der Waals) - compressibility = (1 / Bo_values) * dBo_dP # Resulting array of compressibility values - bulk_mod = 1 / compressibility - - # Find the index of the closest pressure value in b - closest_index = (np.abs(press_range - press)).argmin() - - # Extract the corresponding value from a - pbulk_test = bulk_mod[closest_index] - return pbulk_test - - def _phaseprops(self, fphase, press, Rs=None): - # - # Calculate properties for a single fluid phase - # - # - # Input - # fphase - fluid phase; Oil, Water or Gas - # press - fluid pressure value (MPa) - # - # Output - # pdens - phase density of fluid phase - # "fphase" for pressure value - # "press" (kg/m³) - # pbulk - bulk modulus of fluid phase - # "fphase" for pressure value - # "press" (MPa) - # - # ----------------------------------------------- - # - if fphase.lower() == "oil": - coeffsrho = np.array([0.8, 829.9]) - coeffsbulk = np.array([10.42, 995.79]) - elif fphase.lower() == "wat": - coeffsrho = np.array([0.3, 1067.3]) - coeffsbulk = np.array([9.0, 2807.6]) - elif fphase.lower() == "gas": - coeffsrho = np.array([4.7, 13.4]) - coeffsbulk = np.array([2.75, 0.0]) - else: - print("\nError in phaseprops method") - print("Illegal fluid phase name.") - print("Legal fluid phase names are (case " - "insensitive): Oil, Wat, and Gas.") - sys.exit(1) - # - # Assume simple linear pressure dependencies. - # Coefficients are inferred from - # plots in Batzle and Wang, Geophysics, 1992, - # (where I set the temperature to be 100 degrees - # Celsius, Note also that they give densities in - # g/cc). The resulting straight lines do not fit - # the data extremely well, but they should - # be sufficiently accurate for the purpose of - # this project. - # - pdens = coeffsrho[0]*press + coeffsrho[1] - pbulk = coeffsbulk[0]*press + coeffsbulk[1] - # - return pdens, pbulk - - # - # ======================= - # Fluid properties end - # ======================= - # - - # - # ========================= - # Solid properties start - # ========================= - # - def _solidprops_Johansen(self): - # - # Calculate bulk and shear solid rock (mineral) - # moduli by averaging Hashin-Shtrikman bounds - # - # - # Input - # poro -porosity - # - # Output - # denss - solid rock density (kg/m³) - # bulks - solid rock bulk modulus (unit MPa) - # shears - solid rock shear modulus (unit MPa) - # - # ----------------------------------------------- - # - # - # Solid rock (mineral) density. (Note - # that this is often termed \rho_dry, and not - # \rho_s) - - denss = 2650 # Density of mineral/solid rock kg/m3 - - # - bulks = 37 # (GPa) - shears = 44 # (GPa) - bulks *= 1000 # Convert from GPa to MPa - shears *= 1000 - # - return denss, bulks, shears - # - # - # ======================= - # Solid properties end - # ======================= - # - def _coordination_number(self): - # Applies for granular media - # Average number of contacts that each grain has with surrounding grains - # Coordnumber = 6; simple cubic packing - # Coordnumber = 12; hexagonal close packing - # Needed for Hertz-Mindlin model - # Smeaheia number (Tuhin) - coordnumber = 9 - - return coordnumber - # - def _critical_porosity(self): - # For most porous media there exists a critical porosity - # phi_critical, that seperates their mechanical and acoustic behaviour into two domains. - # For porosities below phi_critical the mineral grains are oad bearing, for values above the grains are - # suspended in the fluids which are load-bearing - # Needed for Hertz-Mindlin model - # Smeaheia number (Tuhin) - phicritical = 0.36 - - return phicritical - # - def _effective_pressure(self, poverb, pfluid): - - # Input - # poverb - overburden pressure (MPa) - # pfluid - fluid pressure (MPa) - - peff = poverb - pfluid - - if peff < 0: - print("\nError in _hertzmindlin method") - print("Negative effective pressure (" + str(peff) + - "). Setting effective pressure to 0.01") - peff = 0.01 - - - - return peff - - # ============================ - # Dry rock properties start - # ============================ - # - def _dryrockmoduli_Smeaheia(self, coordnumber, phicritical, poro, peff, bulks, shears): - # - # - # Calculate bulk and shear dry rock moduli, - - # - # Input - # poro - porosity - # peff - effective pressure overburden - fluid pressure (MPa) - # bulks - bulk solid (mineral) rock bulk - # modulus (MPa) - # shears - solid rock (mineral) shear - # modulus (MPa) - # - # Output - # bulkd - dry rock bulk modulus (unit - # inherited from hertzmindlin and - # variable; bulks) - # sheard - dry rock shear modulus (unit - # inherited from hertzmindlin and - # variable; shears) - # - # ----------------------------------------------- - # - # Calculate Hertz-Mindlin moduli - # - bulkhm, shearhm = self._hertzmindlin_Mavko(peff, bulks, shears, coordnumber, phicritical) - # - bulkd = 1 / ((poro / phicritical) / (bulkhm + 4 / 3 * shearhm) + - (1 - poro / phicritical) / (bulks + 4 / 3 * shearhm)) - 4 / 3 * shearhm - - psi = (9 * bulkhm + 8 * shearhm) / (bulkhm + 2 * shearhm) - - sheard = 1 / ((poro / phicritical) / (shearhm + 1 / 6 * psi * shearhm) + - (1 - poro / phicritical) / (shears + 1 / 6 * psi * shearhm)) - 1 / 6 * psi * shearhm - - #return K_dry, G_dry - return bulkd, sheard - - - # - # --------------------------------------------------- - # - - def _hertzmindlin_Mavko(self, peff, bulks, shears, coordnumber, phicritical): - # - # Calculate bulk and shear Hertz-Mindlin moduli - # adapted from Tuhins kode and "The rock physics handbook", pp247 - # - # - # Input - # p_eff - effective pressure - # bulks - bulk solid (mineral) rock bulk - # modulus (MPa) - # shears - solid rock (mineral) shear - # modulus (MPa) - # coordnumber - average number of contacts that each grain has with surrounding grains - # phicritical - critical porosity - # - # Output - # bulkhm - Hertz-Mindlin bulk modulus - # (MPa) - # shearhm - Hertz-Mindlin shear modulus - # (MPa) - # - # ----------------------------------------------- - # - - - poisson = (3 * bulks - 2 * shears) / (6 * bulks + 2 * shears) - - bulkhm = ((coordnumber ** 2 * (1 - phicritical) ** 2 * shears ** 2 * peff) / - (18 * np.pi ** 2 * (1 - poisson) ** 2)) ** (1 / 3) - shearhm = (5 - 4 * poisson) / (10 - 5 * poisson) * \ - ((3 * coordnumber ** 2 * (1 - phicritical) ** 2 * shears ** 2 * peff) / - (2 * np.pi ** 2 * (1 - poisson) ** 2)) ** (1 / 3) - - - - # - return bulkhm, shearhm - - - # =========================== - # Dry rock properties end - # =========================== - - -if __name__ == '__main__': - # - # Example input with two phases and three grid cells - # - porosity = [0.34999999, 0.34999999, 0.34999999] -# pressure = [ 29.29150963, 29.14003944, 28.88845444] - pressure = [29.3558, 29.2625, 29.3558] -# pressure = [ 25.0, 25.0, 25.0] - phases = ["Oil", "Wat"] -# saturations = [[0.72783828, 0.66568458, 0.58033288], -# [0.27216172, 0.33431542, 0.41966712]] - saturations = [[0.6358, 0.5755, 0.6358], - [0.3641, 0.4245, 0.3641]] -# saturations = [[0.4, 0.5, 0.6], -# [0.6, 0.5, 0.4]] - petElProp = "bulk velocity" - input_dict = {} - input_dict['overburden'] = 'overb.npz' - - print("\nInput:") - print("porosity, pressure:", porosity, pressure) - print("phases, saturations:", phases, saturations) - print("petElProp:", petElProp) - print("input_dict:", input_dict) - - satrock = elasticproperties(input_dict) - - print("overburden:", satrock.overburden) - - satrock.calc_props(phases, saturations, pressure, - porosity) - - print("\nOutput from calc_props:") - print("Density:", satrock.getDens()) - print("Bulk modulus:", satrock.getBulkMod()) - print("Shear modulus:", satrock.getShearMod()) - print("Bulk velocity:", satrock.getBulkVel()) - print("Shear velocity:", satrock.getShearVel()) - print("Bulk impedance:", satrock.getBulkImp()) - print("Shear impedance:", satrock.getShearImp()) - - satrock.getMatchProp(petElProp) - - print("\nOutput from getMatchProp:") - print("Model output selected for data match:", - satrock.match_prop) diff --git a/simulator/rockphysics/standardrp.py b/simulator/rockphysics/standardrp.py deleted file mode 100644 index dde7b09d..00000000 --- a/simulator/rockphysics/standardrp.py +++ /dev/null @@ -1,795 +0,0 @@ -"""Descriptive description.""" - -__author__ = 'TM' - -# standardrp.py -import numpy as np -import sys -import multiprocessing as mp -# 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. - - Examples - -------- - >>> 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() - 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 a 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): - ### - # This doesn't initialize for models with no uncertainty - ### - # # if some PEM properties have uncertainty, set the correct value - # if ensembleMember is not None: - # for pem_state in self.__dict__.keys(): # loop over all possible pem vaules - # if pem_state is not 'inv_state': # do not alter the ensemble - # if type(eval('self.{}'.format(pem_state))) is dict: - # for el in eval('self.{}'.format(pem_state)).keys(): - # if type(eval('self.{}'.format(pem_state))[el]) is dict: - # for param in eval('self.{}'.format(pem_state))[el].keys(): - # if param in self.inv_state: - # eval('self.{}'.format(pem_state))[el][param]=\ - # self.inv_state[param][:, ensembleMember] - # elif param + '_' + el in self.inv_state: - # eval('self.{}'.format(pem_state))[el][param] = \ - # self.inv_state[param+'_' + el][:, ensembleMember] - # else: - # if el in self.inv_state: - # eval('self.{}'.format(pem_state))[el] = self.inv_state[el][:,ensembleMember] - # else: - # if pem_state in self.inv_state: - # setattr(self,pem_state, self.inv_state[pem_state][:,ensembleMember]) - - # Check if the inputs are given as a list (more - # than one phase) or a single input (single - # phase). If single phase input, make the input a - # list with a single entry (s.t. it can be used - # directly with the methods below) - # - 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" into local variable to - # comply with remaining code parts - overburden = self.overburden - - if press_init is None: - p_init = self.p_init - else: - p_init = press_init - # - poverburden = self.overburden - - # debug - self.pressure = pressure - self.peff = poverburden - pressure - self.porosity = porosity - - # Check that no. of phases is equal to no. of - # entries in saturations list - # - 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) - # - # 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) - # ------------------------------- - # Calculate saturated properties - # ------------------------------- - # - # Density - # - self.dens[i] = (porosity[i]*densf + (1-porosity[i])*denss) - # - # Moduli - # - self.bulkmod[i] = \ - bulkd + (1 - bulkd/bulks)**2 / \ - (porosity[i]/bulkf + - (1-porosity[i])/bulks - - bulkd/(bulks**2)) - self.shearmod[i] = sheard - - # Velocities (due to bulk/shear modulus being - # in MPa, we multiply by 1000 to get m/s - # instead of km/s) - # - self.bulkvel[i] = \ - 1000*np.sqrt((abs(self.bulkmod[i]) + - 4*self.shearmod[i]/3)/(self.dens[i])) - self.shearvel[i] = \ - 1000*np.sqrt(self.shearmod[i] / - (self.dens[i])) - # - # Impedances (m/s)*(Kg/m3) - # - self.bulkimp[i] = self.dens[i] * \ - self.bulkvel[i] - self.shearimp[i] = self.dens[i] * \ - self.shearvel[i] - - def getMatchProp(self, petElProp): - if petElProp.lower() == 'density': - self.match_prop = self.getDens() - elif petElProp.lower() == 'bulk_modulus': - self.match_prop = self.getBulkMod() - elif petElProp.lower() == 'shear_modulus': - self.match_prop = self.getShearMod() - elif petElProp.lower() == 'bulk_velocity': - self.match_prop = self.getBulkVel() - elif petElProp.lower() == 'shear_velocity': - self.match_prop = self.getShearVel() - elif petElProp.lower() == "bulk_impedance": - self.match_prop = self.getBulkImp() - elif petElProp.lower() == 'shear_impedance': - self.match_prop = self.getShearImp() - else: - print("\nError in getMatchProp method") - print("No model output type selected for " - "data match.") - print("Legal model output types are " - "(case insensitive):") - print("Density, bulk modulus, shear " - "modulus, bulk velocity,") - print("shear velocity, bulk impedance, " - "shear impedance") - sys.exit(1) - return self.match_prop - - def getDens(self): - return self.dens - - def getBulkMod(self): - return self.bulkmod - - def getShearMod(self): - return self.shearmod - - def getBulkVel(self): - return self.bulkvel - - def getShearVel(self): - return self.shearvel - - def getBulkImp(self): - return self.bulkimp - - def getShearImp(self): - return self.shearimp - - def getOverburdenP(self): - return self.overburden - - def getPressure(self): - return self.pressure - - def getPeff(self): - return self.peff - - def getPorosity(self): - return self.porosity - # - # =================================================== - # Fluid properties start - # =================================================== - # - def _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(self, poro, ntg=None, ind=None): - # - # 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 - # inherited from hashinshtr) - # shears - solid rock shear modulus (unit - # inherited from hashinshtr) - # - # ----------------------------------------------- - # - # From PetroWiki (kg/m³) - # - densc = 2540 - densq = 2650 - # - # From "Step 1" of "recipe" in Report 1 in Abul - # Fahimuddin's thesis. - # - vclay = 0.7 - 1.58*poro - # - # Calculate solid rock (mineral) density. (Note - # that this is often termed \rho_dry, and not - # \rho_s) - # - denss = densq + vclay*(densc - densq) - # - # Calculate lower and upper bulk and shear - # Hashin-Shtrikman bounds - # - bulkl, bulku, shearl, shearu = \ - self._hashinshtr(vclay) - # - # Calculate bulk and shear solid rock (mineral) - # moduli as arithmetic means of the respective - # bounds - # - bulkb = np.array([bulkl, bulku]) - shearb = np.array([shearl, shearu]) - bulks = np.mean(bulkb) - shears = np.mean(shearb) - # - return denss, bulks, shears - # - # --------------------------------------------------- - # - - def _hashinshtr(self, vclay): - # - # Calculate lower and upper, bulk and shear, - # Hashin-Shtrikman bounds, utilizing that they - # all have the common mathematical form, - # - # f = a + b/((1/c) + d*(1/e)). - # - # - # Input - # vclay - "volume of clay" - # - # Output - # bulkl - lower bulk Hashin-Shtrikman - # bound (MPa) - # bulku - upper bulk Hashin-Shtrikman - # bound (MPa) - # shearl - lower shear Hashin-Shtrikman - # bound (MPa) - # shearu - upper shear Hashin-Shtrikman - # bound (MPa) - # - # ----------------------------------------------- - # - # From table 1 in Report 1 in Abul Fahimuddin's - # thesis (he used GPa, I use MPa), ("c" for clay - # and "q" for quartz.): - # - bulkc = 14900 - bulkq = 37000 - shearc = 1950 - shearq = 44000 - # - # Calculate quantities common for both bulk and - # shear formulas - # - lb = 1 - vclay - ub = vclay - le = bulkc + 4*shearc/3 - ue = bulkq + 4*shearq/3 - # - # Calculate quantities common only for bulk - # formulas - # - bld = vclay - bud = 1 - vclay - blc = bulkq - bulkc - buc = bulkc - bulkq - # - # Calculate quantities common only for shear - # formulas - # - sld = 2*vclay*(bulkc + 2*shearc)/(5*shearc) - sud = 2*(1 - vclay)*(bulkq + 2*shearq)/(5*shearq) - slc = shearq - shearc - suc = shearc - shearq - # - # Calculate bounds utilizing generic formula; - # f = a + b/((1/c) + d*(1/e)). - # - # Lower bulk - # - bulkl = self._genhashinshtr(bulkc, lb, blc, bld, - le) - # - # Upper bulk - # - bulku = self._genhashinshtr(bulkq, ub, buc, bud, - ue) - # - # Lower shear - # - shearl = self._genhashinshtr(shearc, lb, slc, - sld, le) - # - # Upper shear - # - shearu = self._genhashinshtr(shearq, ub, suc, - sud, ue) - # - return bulkl, bulku, shearl, shearu - - # - # --------------------------------------------------- - # - def _genhashinshtr(self, a, b, c, d, e): - # - # Calculate arbitrary Hashin-Shtrikman bound, - # which has the generic form - # - # f = a + b/((1/c) + d*(1/e)) - # - # both for lower and upper bulk and shear bounds - # - # - # Input - # a - see above formula - # b - see above formula - # c - see above formula - # d - see above formula - # e - see above formula - # - # Output - # f - Bulk or shear Hashin-Shtrikman bound - # value - # - # ----------------------------------------------- - # - cinv = 1/c - einv = 1/e - f = a + b/(cinv + d*einv) - # - return f - # - # ======================= - # Solid properties end - # ======================= - # - - # - # ============================ - # Dry rock properties start - # ============================ - # - def _dryrockmoduli(self, poro, poverb, pfluid, bulks, shears, ind=None, ntg=None, p_init=None, denss=None, Rs=None, phases=None): - # - # - # Calculate bulk and shear dry rock moduli, - # utilizing that they have the common - # mathematical form, - # - # -- -- ^(-1) - # |(poro/poroc) 1 - (poro/poroc)| - # f = |------------ + ----------------| - z. - # | a + z b + z | - # --. -- - # - # - # Input - # poro - porosity - # poverb - overburden pressure (MPa) - # pfluid - 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(poverb, - pfluid) - # - # From table 1 in Report 1 in Abul Fahimuddin's - # thesis (I assume \phi_max in that - # table corresponds to \phi_c in his formulas): - # - poroc = 0.4 - # - # Calculate input common to both bulk and - # shear formulas - # - poratio = poro/poroc - # - # Calculate input to bulk formula - # - ba = bulkhm - bb = bulks - bz = 4*shearhm/3 - # - # Calculate input to shear formula - # - sa = shearhm - sb = shears - sz = (shearhm/6)*((9*bulkhm + 8*shearhm) / - (bulkhm + 2*shearhm)) - # - # Calculate moduli - # - bulkd = self._gendryrock(poratio, ba, bb, bz) - sheard = self._gendryrock(poratio, sa, sb, sz) - # - return bulkd, sheard - # - # --------------------------------------------------- - # - - def _hertzmindlin(self, poverb, pfluid): - # - # Calculate bulk and shear Hertz-Mindlin moduli - # utilizing that they have the common - # mathematical form, - # - # f = max*(peff/pref)^kappa. - # - # - # Input - # poverb - overburden pressure - # pfluid - fluid pressure - # - # Output - # bulkhm - Hertz-Mindlin bulk modulus - # (MPa) - # shearhm - Hertz-Mindlin shear modulus - # (MPa) - # - # ----------------------------------------------- - # - # From table 1 in Report 1 in Abul Fahimuddin's - # thesis (he used GPa for the moduli, I use MPa - # also for them): - # - bulkmax = 3310 - shearmax = 2840 - pref = 8.8 - kappa = 0.233 - # - # Calculate moduli - # - 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 - # sys.exit(1) - common = (peff/pref)**kappa - bulkhm = bulkmax*common - shearhm = shearmax*common - # - return bulkhm, shearhm - # - # --------------------------------------------------- - # - - def _gendryrock(self, q, a, b, z): - # - # Calculate arbitrary dry rock moduli, which has - # the generic form - # - # -- -- ^(-1) - # | q 1 - q | - # f = |------------ + ----------------| - z, - # | a + z b + z | - # --. -- - # - # both for bulk and shear moduli - # - # - # Input - # q - see above formula - # a - see above formula - # b - see above formula - # z - see above formula - # - # Output - # f - Bulk or shear dry rock modulus value - # - # ----------------------------------------------- - # - afrac = q/(a + z) - bfrac = (1 - q)/(b + z) - f = 1/(afrac + bfrac) - z - # - return f - # =========================== - # 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)