simulator.eclipse

   1"""Wrap Eclipse"""
   2# External imports
   3import numpy as np
   4import sys
   5import os
   6from copy import deepcopy
   7from mako.lookup import TemplateLookup
   8from mako.runtime import Context
   9from multiprocessing import Process
  10import datetime as dt
  11from scipy import interpolate
  12from subprocess import call, DEVNULL
  13from misc import ecl, grdecl
  14from shutil import rmtree, copytree  # rmtree for removing folders
  15import time
  16# import rips
  17
  18# Internal imports
  19from misc.system_tools.environ_var import EclipseRunEnvironment
  20from pipt.misc_tools.analysis_tools import store_ensemble_sim_information
  21
  22
  23class eclipse:
  24    """
  25    Class for running the Schlumberger eclipse 100 black oil reservoir simulator. For more information see  GeoQuest:
  26    ECLIPSE reference manual 2009.1. Schlumberger, GeoQuest (2009).
  27
  28    To run this class, eclipse must be installed and elcrun must be in the system path!
  29    """
  30
  31    def __init__(self, input_dict=None, filename=None, options=None):
  32        """
  33        The inputs are all optional, but in the same fashion as the other simulators a system must be followed.
  34        The input_dict can be utilized as a single input. Here all nescessary info is stored. Alternatively,
  35        if input_dict is not defined, all the other input variables must be defined.
  36
  37        Parameters
  38        ----------
  39        input_dict : dict, optional
  40            Dictionary containing all information required to run the simulator.
  41
  42                - parallel: number of forward simulations run in parallel
  43                - simoptions: options for the simulations
  44                    - mpi: option to use mpi (always use > 2 cores)
  45                    - sim_path: Path to the simulator
  46                    - sim_flag: Flags sent to the simulator (see simulator documentation for all possibilities)
  47                - sim_limit: maximum number of seconds a simulation can run before being killed
  48                - runfile: name of the simulation input file
  49                - reportpoint: these are the dates the simulator reports results
  50                - reporttype: this key states that the report poins are given as dates
  51                - datatype: the data types the simulator reports
  52                - replace: replace failed simulations with randomly selected successful ones
  53                - rerun: in case of failure, try to rerun the simulator the given number of times
  54                - startdate: simulaton start date
  55                - saveforecast: save the predicted measurements for each iteration
  56
  57        filename : str, optional
  58            Name of the .mako file utilized to generate the ECL input .DATA file. Must be in uppercase for the
  59            ECL simulator.
  60
  61        options : dict, optional
  62            Dictionary with options for the simulator.
  63
  64        Returns
  65        -------
  66        initial_object : object
  67            Initial object from the class ecl_100.
  68        """
  69
  70        # IN
  71        self.input_dict = input_dict
  72        self.file = filename
  73        self.options = options
  74
  75        self.upscale = None
  76        # If input option 1 is selected
  77        if self.input_dict is not None:
  78            self._extInfoInputDict()
  79
  80        # Allocate for internal use
  81        self.inv_stat = None
  82        self.static_state = None
  83
  84    def _extInfoInputDict(self):
  85        """
  86        Extract the manditory and optional information from the input_dict dictionary.
  87
  88        Parameters
  89        ----------
  90        input_dict : dict
  91            Dictionary containing all information required to run the simulator (defined in self).
  92
  93        Returns
  94        -------
  95        filename : str
  96            Name of the .mako file utilized to generate the ECL input .DATA file.
  97        """
  98        # Chech for mandatory keys
  99        assert 'reporttype' in self.input_dict, 'Reporttype is missing, please specify this'
 100        assert 'reportpoint' in self.input_dict, 'Reportpoint is missing, please specify this'
 101
 102        self.true_prim = [self.input_dict['reporttype'], self.input_dict['reportpoint']]
 103
 104        self.true_order = [self.input_dict['reporttype'], self.input_dict['reportpoint']]
 105        self.all_data_types = self.input_dict['datatype']
 106        self.l_prim = [int(i) for i in range(len(self.true_prim[1]))]
 107
 108        # In the ecl framework, all reference to the filename should be uppercase
 109        self.file = self.input_dict['runfile'].upper()
 110        self.options = {}
 111        self.options['sim_path'] = ''
 112        self.options['sim_flag'] = ''
 113        self.options['mpi'] = ''
 114        self.options['parsing-strictness'] = ''
 115        # Loop over options in SIMOPTIONS and extract the parameters we want
 116        if 'simoptions' in self.input_dict:
 117            if type(self.input_dict['simoptions'][0]) == str:
 118                self.input_dict['simoptions'] = [self.input_dict['simoptions']]
 119            for i, opt in enumerate(list(zip(*self.input_dict['simoptions']))[0]):
 120                if opt == 'sim_path':
 121                    self.options['sim_path'] = self.input_dict['simoptions'][i][1]
 122                if opt == 'sim_flag':
 123                    self.options['sim_flag'] = self.input_dict['simoptions'][i][1]
 124                if opt == 'mpi':
 125                    self.options['mpi'] = self.input_dict['simoptions'][i][1]
 126                if opt == 'parsing-strictness':
 127                    self.options['parsing-strictness'] = self.input_dict['simoptions'][i][1]
 128        if 'sim_limit' in self.input_dict:
 129            self.options['sim_limit'] = self.input_dict['sim_limit']
 130
 131        if 'reportdates' in self.input_dict:
 132            self.reportdates = [
 133                x * 30 for x in range(1, int(self.input_dict['reportdates'][1]))]
 134
 135        if 'read_sch' in self.input_dict:
 136            # self.read_sch = self.input_dict['read_sch'][0]
 137            load_file = np.load(self.input_dict['read_sch'], allow_pickle=True)
 138            self.reportdates = load_file[load_file.files[0]]
 139
 140        if 'startdate' in self.input_dict:
 141            self.startDate = {}
 142            # assume date is on form day/month/year
 143            tmpDate = [int(elem) for elem in self.input_dict['startdate'].split('/')]
 144
 145            self.startDate['day'] = tmpDate[0]
 146            self.startDate['month'] = tmpDate[1]
 147            self.startDate['year'] = tmpDate[2]
 148
 149        if 'realizations' in self.input_dict:
 150            self.realizations = self.input_dict['realizations']
 151
 152        if 'trunc_level' in self.input_dict:
 153            self.trunc_level = self.input_dict['trunc_level']
 154
 155        if 'rerun' in self.input_dict:
 156            self.rerun = int(self.input_dict['rerun'])
 157        else:
 158            self.rerun = 0
 159
 160        # If we want to extract, or evaluate, something uniquely from the ensemble specific run we can
 161        # run a user defined code to do this.
 162        self.saveinfo = None
 163        if 'savesiminfo' in self.input_dict:
 164            # Make sure "ANALYSISDEBUG" gives a list
 165            if isinstance(self.input_dict['savesiminfo'], list):
 166                self.saveinfo = self.input_dict['savesiminfo']
 167            else:
 168                self.saveinfo = [self.input_dict['savesiminfo']]
 169
 170        if 'upscale' in self.input_dict:
 171            self.upscale = {}
 172            for i in range(0, len(self.input_dict['upscale'])):
 173                if self.input_dict['upscale'][i][0] == 'state':
 174                    # Set the parameter we upscale with regards to, must be the same as one free parameter
 175                    self.upscale['state'] = self.input_dict['upscale'][i][1]
 176                    # Set the dimension of the parameterfield to be upscaled (x and y)
 177                    self.upscale['dim'] = self.input_dict['upscale'][i][2]
 178                if self.input_dict['upscale'][i][0] == 'maxtrunc':
 179                    # Set the ratio for the maximum truncation value
 180                    self.upscale['maxtrunc'] = self.input_dict['upscale'][i][1]
 181                if self.input_dict['upscale'][i][0] == 'maxdiff':
 182                    # Set the ratio for the maximum truncation of differences
 183                    self.upscale['maxdiff'] = self.input_dict['upscale'][i][1]
 184                if self.input_dict['upscale'][i][0] == 'wells':
 185                    # Set the list of well indexes, this is a list of lists where each element in the outer list gives
 186                    # a well coordinate (x and y) as elements in the inner list.
 187                    self.upscale['wells'] = []
 188                    for j in range(1, len(self.input_dict['upscale'][i])):
 189                        self.upscale['wells'].append(
 190                            [int(elem) for elem in self.input_dict['upscale'][i][j]])
 191                if self.input_dict['upscale'][i][0] == 'radius':
 192                    # List of radius lengths
 193                    self.upscale['radius'] = [int(elem)
 194                                              for elem in self.input_dict['upscale'][i][1]]
 195                if self.input_dict['upscale'][i][0] == 'us_type':
 196                    self.upscale['us_type'] = self.input_dict['upscale'][i][1]
 197
 198            # Check that we have as many radius elements as wells
 199            if 'radius' in self.upscale:
 200                if len(self.upscale['radius']) != len(self.upscale['wells']):
 201                    sys.exit('ERROR: Missmatch between number of well inputs and number of radius elements. Please check '
 202                             'the input file')
 203            else:
 204                self.upscale['radius'] = []
 205                self.upscale['wells'] = []
 206
 207        # The simulator should run on different levels
 208        if 'multilevel' in self.input_dict:
 209            # extract list of levels
 210            self.multilevel = self.input_dict['multilevel']
 211        else:
 212            # if not initiallize as list with one element
 213            self.multilevel = [False]
 214
 215    def setup_fwd_run(self, **kwargs):
 216        """
 217        Setup the simulator.
 218
 219        Parameters
 220        ----------
 221        assimIndex: int
 222            Gives the index-type (e.g. step,time,etc.) and the index for the
 223            data to be assimilated
 224        trueOrder:
 225            Gives the index-type (e.g. step,time,etc.) and the index of the true data
 226        """
 227        self.__dict__.update(kwargs)  # parse kwargs input into class attributes
 228
 229        if hasattr(self, 'reportdates'):
 230            self.report = {'dates': self.reportdates}
 231        elif 'reportmonths' in self.input_dict:  # for optimization
 232            self.report = {
 233                'days': [30 * i for i in range(1, int(self.input_dict['reportmonths'][1]))]}
 234        else:
 235            assimIndex = [i for i in range(len(self.l_prim))]
 236            trueOrder = self.true_order
 237
 238        self.pred_data = [deepcopy({}) for _ in range(len(assimIndex))]
 239        for ind in self.l_prim:
 240            for key in self.all_data_types:
 241                self.pred_data[ind][key] = np.zeros((1, 1))
 242
 243        if isinstance(trueOrder[1], list):  # Check if true data prim. ind. is a list
 244            self.true_prim = [trueOrder[0], [x for x in trueOrder[1]]]
 245        else:  # Float
 246            self.true_prim = [trueOrder[0], [trueOrder[1]]]
 247        # self.all_data_types = list(pred_data[0].keys())
 248
 249        # Initiallise space to store the number of active cells. This is only for the upscaling option.
 250        if 'upscale' in self.input_dict:
 251            self.num_act = []
 252
 253        # Initialise error summary
 254        self.error_smr = []
 255
 256        # Initiallize run time summary
 257        self.run_time = []
 258
 259        # Check that the .mako file is in the current working directory
 260        if not os.path.isfile('%s.mako' % self.file):
 261            sys.exit(
 262                'ERROR: .mako file is not in the current working directory. This file must be defined')
 263
 264    def run_fwd_sim(self, state, member_i, del_folder=True):
 265        """
 266        Setup and run the ecl_100 forward simulator. All the parameters are defined as attributes, and the name of the
 267        parameters are initialized in setupFwdRun. This method sets up and runs all the individual ensemble members.
 268        This method is based on writing .DATA file for each ensemble member using the mako template, for more info
 269        regarding mako see http://www.makotemplates.org/
 270
 271        Parameters
 272        ----------
 273        state : dict
 274            Dictionary containing the ensemble state.
 275
 276        member_i : int
 277            Index of the ensemble member.
 278
 279        del_folder : bool, optional
 280            Boolean to determine if the ensemble folder should be deleted. Default is False.
 281        """
 282        if hasattr(self, 'level'):
 283            state['level'] = self.level
 284        else:
 285            state['level'] = 0  # default value
 286        os.mkdir('En_' + str(member_i))
 287        folder = 'En_' + str(member_i) + os.sep
 288
 289        state['member'] = member_i
 290        # If the run is upscaled, run the upscaling procedure
 291        if self.upscale is not None:
 292            if hasattr(self, 'level'):  # if this is a multilevel run, we must set the level
 293                self.upscale['maxtrunc'] = self.trunc_level[self.level]
 294            if self.upscale['maxtrunc'] > 0:  # if the truncation level is 0, we do not perform upscaling
 295                self.coarsen(folder, state)
 296            # if the level is 0, and upscaling has been performed earlier this must be
 297            elif hasattr(self, 'coarse'):
 298                # removed
 299                del self.coarse
 300
 301        # start by generating the .DATA file, using the .mako template situated in ../folder
 302        self._runMako(folder, state)
 303        success = False
 304        rerun = self.rerun
 305        while rerun >= 0 and not success:
 306            success = self.call_sim(folder, True)
 307            rerun -= 1
 308        if success:
 309            self.extract_data(member_i)
 310            if del_folder:
 311                if self.saveinfo is not None:  # Try to save information
 312                    store_ensemble_sim_information(self.saveinfo, member_i)
 313                self.remove_folder(member_i)
 314            return self.pred_data
 315        else:
 316            if self.redund_sim is not None:
 317                success = self.redund_sim.call_sim(folder, True)
 318                if success:
 319                    self.extract_data(member_i)
 320                    if del_folder:
 321                        if self.saveinfo is not None:  # Try to save information
 322                            store_ensemble_sim_information(self.saveinfo, member_i)
 323                        self.remove_folder(member_i)
 324                    return self.pred_data
 325                else:
 326                    if del_folder:
 327                        self.remove_folder(member_i)
 328                    return False
 329            else:
 330                if del_folder:
 331                    self.remove_folder(member_i)
 332                return False
 333
 334    def remove_folder(self, member):
 335        folder = 'En_' + str(member) + os.sep
 336        try:
 337            rmtree(folder)  # Try to delete folder
 338        except:  # If deleting fails, just rename to 'En_<ensembleMember>_date'
 339            os.rename(folder, folder + '_' +
 340                      dt.datetime.now().strftime("%m-%d-%Y_%H-%M-%S"))
 341
 342    def extract_data(self, member):
 343        # get the formated data
 344        for prim_ind in self.l_prim:
 345            # Loop over all keys in pred_data (all data types)
 346            for key in self.all_data_types:
 347                if self.pred_data[prim_ind][key] is not None:  # Obs. data at assim. step
 348                    true_data_info = [self.true_prim[0], self.true_prim[1][prim_ind]]
 349                    try:
 350                        data_array = self.get_sim_results(key, true_data_info, member)
 351                        self.pred_data[prim_ind][key] = data_array
 352                    except:
 353                        pass
 354
 355    def coarsen(self, folder, ensembleMember=None):
 356        """
 357        This method utilized one field parameter to upscale the computational grid. A coarsening file is written to the
 358        ensemble folder, and the eclipse calculates the upscaled permeabilities, porosities and transmissibilities
 359        based on the new grid and the original parameter values.
 360
 361        Parameters
 362        ----------
 363        Folder : str, optional
 364            Path to the ecl_100 run folder.
 365
 366        ensembleMember : int, optional
 367            Index of the ensemble member to run.
 368
 369        Changelog
 370        ---------
 371        - KF 17/9-2015 Added uniform upscaling as an option
 372        - KF 6/01-17
 373        """
 374
 375        # Get the parameter values for the current ensemble member. Note that we select only one parameter
 376        if ensembleMember is not None:
 377            coarsenParam = self.inv_state[self.upscale['state']][:, ensembleMember]
 378        else:
 379            coarsenParam = self.inv_state[self.upscale['state']]
 380
 381        if self.upscale['us_type'].lower() == 'haar':  # the upscaling is based on a Haar wavelet
 382            # Do not add any dead-cells
 383            # TODO: Make a better method for determining dead or inactive cells.
 384
 385            orig_wght = np.ones((self.upscale['dim'][0], self.upscale['dim'][1]))
 386            # Get the new cell structure by a 2-D unbalanced Haar transform.
 387            wght = self._Haar(coarsenParam.reshape(
 388                (self.upscale['dim'][0], self.upscale['dim'][1])), orig_wght)
 389
 390        elif self.upscale['us_type'].lower() == 'unif':
 391            # This option generates a grid where the cells are upscaled uniformly in a dyadic manner. We utilize the
 392            # same framework to define the parameters to be coarsened as in the Haar case. Hence, we only need to
 393            # calculate the wght variable. The level of upscaling can be determined by a fraction, however, for this
 394            # case the fraction gives the number of upscaling levels, 1 is all possible levels, 0 is no upscaling.
 395            wght = self._unif((self.upscale['dim'][0], self.upscale['dim'][1]))
 396
 397        self.write_coarse(folder, wght, coarsenParam.reshape(
 398            (int(self.upscale['dim'][0]), int(self.upscale['dim'][1]))))
 399
 400    def write_coarse(self, folder, whgt, image):
 401        """
 402        This function writes the include file coarsen to the ecl run. This file tels ECL to coarsen the grid.
 403
 404        Parameters
 405        ----------
 406        folder : str
 407            Path to the ecl_100 run.
 408
 409        whgt : float
 410            Weight of the transformed cells.
 411
 412        image : array-like
 413            Original image.
 414
 415        Changelog
 416        ---------
 417        - KF 17/9-2015
 418        """
 419
 420        well = self.upscale['wells']
 421        radius = self.upscale['radius']
 422
 423        well_cells = self._nodeIndex(image.shape[1], image.shape[0], well, radius)
 424        coarse = np.array([[True] * image.shape[1]] * image.shape[0])
 425        tmp = coarse.flatten()
 426        tmp[well_cells] = False
 427        coarse = tmp.reshape(image.shape)
 428
 429        # f = open(folder + 'coarsen.dat', 'w+')
 430        # f.write('COARSEN \n')
 431        ecl_coarse = []
 432
 433        for level in range(len(whgt) - 1, -1, -1):
 434            weight_at_level = whgt[level]
 435            x_dim = weight_at_level.shape[0]
 436            y_dim = weight_at_level.shape[1]
 437
 438            merged_at_level = weight_at_level[int(x_dim / 2):, :int(y_dim / 2)]
 439
 440            level_dim = 2 ** (level + 1)
 441
 442            for i in range(0, merged_at_level.shape[0]):
 443                for j in range(0, merged_at_level.shape[1]):
 444                    if merged_at_level[i, j] == 1:
 445                        # Do this in two stages to avoid errors with the dimensions
 446                        #
 447                        # only merging square cells
 448                        if (i + 1) * level_dim <= image.shape[0] and (j + 1) * level_dim <= image.shape[1]:
 449                            if coarse[i * level_dim:(i + 1) * level_dim, j * level_dim:(j + 1) * level_dim].all():
 450                                coarse[i * level_dim:(i + 1) * level_dim,
 451                                       j * level_dim:(j + 1) * level_dim] = False
 452                                ecl_coarse.append([j * level_dim + 1, (j + 1) * level_dim,
 453                                                   i * level_dim + 1, (i + 1) * level_dim, 1, 1, 1, 1, 1])
 454                                # f.write('%i %i %i %i 1 1 1 1 1 / \n' % (j * level_dim + 1, (j + 1) * level_dim,
 455                                #                                         i * level_dim + 1, (i + 1) * level_dim))
 456                        # cells at first edge, non-square
 457                        if (i + 1)*level_dim > image.shape[0] and (j + 1) * level_dim <= image.shape[1]:
 458                            if coarse[i * level_dim::, j * level_dim:(j + 1) * level_dim].all():
 459                                coarse[i * level_dim::, j *
 460                                       level_dim:(j + 1) * level_dim] = False
 461                                ecl_coarse.append([j * level_dim + 1, (j + 1) * level_dim,
 462                                                   i * level_dim + 1, image.shape[0], 1, 1, 1, 1, 1])
 463                                # f.write('%i %i %i %i 1 1 1 1 1 / \n' % (j * level_dim + 1, (j + 1) * level_dim,
 464                                #                                         i * level_dim + 1, image.shape[0]))
 465                        # cells at second edge, non-square
 466                        if (j + 1)*level_dim > image.shape[1] and (i + 1) * level_dim <= image.shape[0]:
 467                            if coarse[i * level_dim:(i + 1) * level_dim, j * level_dim::].all():
 468                                coarse[i * level_dim:(i + 1) * level_dim,
 469                                       j * level_dim::] = False
 470                                ecl_coarse.append([j * level_dim + 1, image.shape[1],
 471                                                   i * level_dim + 1, (i + 1) * level_dim, 1, 1, 1, 1, 1])
 472                                # f.write('%i %i %i %i 1 1 1 1 1 / \n' % (j * level_dim + 1, image.shape[1],
 473                                #                                         i * level_dim + 1, (i + 1) * level_dim))
 474                        # cells at intersection between first and second edge, non-square
 475                        if (i + 1) * level_dim > image.shape[0] and (j + 1) * level_dim > image.shape[1]:
 476                            if coarse[i * level_dim::, j * level_dim::].all():
 477                                coarse[i * level_dim::, j * level_dim::] = False
 478                                ecl_coarse.append([j * level_dim + 1, image.shape[1],
 479                                                   i * level_dim + 1, image.shape[0], 1, 1, 1, 1, 1])
 480                                # f.write('%i %i %i %i 1 1 1 1 1 / \n' % (j * level_dim + 1, image.shape[1],
 481                                #                                         i * level_dim + 1, image.shape[0]))
 482        # f.write('/')
 483        # f.close()
 484        self.coarse = ecl_coarse
 485
 486    def _nodeIndex(self, x_dir, y_dir, cells, listRadius):
 487        # Find the node index for the cells, and the cells in a radius around.
 488        index = []
 489        for k in range(0, len(cells)):
 490            cell_x = cells[k][0] - 1  # Remove one to make python equivalent to ecl format
 491            cell_y = cells[k][1] - 1  # Remove one to make python equivalent to ecl format
 492            radius = listRadius[k]
 493            # tmp_index = [cell_x + cell_y*x_dir]
 494            tmp_index = []
 495            for i in range(-radius, radius + 1):
 496                for j in range(-radius, radius + 1):
 497                    x = cell_x + i
 498                    y = cell_y + j
 499                    if (0 <= x < x_dir and 0 <= y < y_dir) and (np.sqrt(i ** 2 + j ** 2) <= radius) \
 500                            and (x + y * x_dir < y_dir * x_dir):
 501                        tmp_index.append(x + y * x_dir)
 502            index.extend(tmp_index)
 503
 504        return index
 505
 506    def _Haar(self, image, weights):
 507        """
 508        The 2D unconditional Haar transform.
 509
 510        Parameters
 511        ----------
 512        image : array-like
 513            The original image where the Haar-transform is performed. This could be the permeability or the porosity field.
 514
 515        orig_weights : array-like
 516            The original weights of the image. If some cells are inactive or dead, their weight can be set to zero.
 517
 518        Returns
 519        -------
 520        tot_weights : array-like
 521            A vector/matrix of the same size as the image, with the weights of the new cells. This is utilized for writing
 522            the coarsening file but is not sufficient to recreate the image.
 523        """
 524        # Todo: make an option which stores the transformed field.
 525
 526        # Define the two truncation levels as a fraction of the larges value in the original image. The fraction is defined
 527        # through the inputs.
 528        max_value = self.upscale['maxtrunc'] * np.max(image)
 529        max_diff = self.upscale['maxdiff'] * (np.max(image) - np.min(image))
 530
 531        tot_transf = []
 532        tot_weight = []
 533        tot_alpha = []
 534        count = 0
 535        allow_merge = np.ones(image.shape)
 536        while image.shape[0] > 1:
 537            level_alpha = []
 538            if isinstance(image.shape, tuple):
 539                # the image has two axes.
 540                # Perform the transformation for all the coulmns, and then for the transpose of the results
 541
 542                for axis in range(0, 2):
 543                    # Initialise the matrix for storing the details and smoothing
 544                    level_diff_columns = np.empty(
 545                        (int(np.ceil(len(image[:, 0]) / 2)), len(image[0, :])))
 546                    # level_diff_rows = np.empty((len(image[:, 0]), int(np.ceil(len(image[0, :])/2))))
 547                    level_smooth_columns = np.empty(
 548                        (int(np.ceil(len(image[:, 0]) / 2)), len(image[0, :])))
 549                    level_weight_columns = np.empty(
 550                        (int(np.ceil(len(image[:, 0]) / 2)), len(image[0, :])))
 551                    level_alpha_columns = np.empty(
 552                        (int(np.ceil(len(image[:, 0]) / 2)), len(image[0, :])))
 553
 554                    for i in range(0, image.shape[1]):
 555                        if len(image[:, i]) % 2 == 0:
 556                            diff = image[1::2, i] - image[::2, i]
 557                            new_weights = weights[::2, i] + weights[1::2, i]
 558                            alpha = (weights[1::2, i] / new_weights)
 559                            smooth = image[::2, i] + alpha * diff
 560
 561                        else:
 562                            tmp_img = image[:-1, i]
 563                            diff = tmp_img[1::2] - tmp_img[::2]
 564
 565                            tmp_weight = weights[:-1, i]
 566                            new_weights = tmp_weight[::2] + tmp_weight[1::2]
 567                            new_weights = np.append(new_weights, weights[-1, i])
 568
 569                            tmp_img = image[:-1, i]
 570                            alpha = (weights[1::2, i] / new_weights[:-1])
 571                            smooth = tmp_img[::2] + alpha * diff
 572                            smooth = np.append(smooth, image[-1, i])
 573                            alpha = np.append(alpha, 0)
 574                            diff = np.append(diff, 0)
 575
 576                        level_diff_columns[:, i] = diff
 577                        level_weight_columns[:, i] = new_weights
 578                        level_smooth_columns[:, i] = smooth
 579                        level_alpha_columns[:, i] = alpha
 580
 581                    # image = level_smooth_columns.T
 582                    image = np.vstack((level_smooth_columns, level_diff_columns)).T
 583                    weights = np.vstack((level_weight_columns, level_weight_columns)).T
 584                    level_alpha.append(level_alpha_columns)
 585
 586                image, weights, allow_merge, end_ctrl = self._haarTrunc(image, weights, max_value, max_diff,
 587                                                                        allow_merge)
 588
 589                tot_transf.append(image)
 590                tot_weight.append(weights)
 591                tot_alpha.append(level_alpha)
 592
 593                if end_ctrl:
 594                    break
 595                image = image[:image.shape[0] / 2, :image.shape[1] / 2]
 596                weights = weights[:weights.shape[0] / 2, :weights.shape[1] / 2]
 597
 598            else:
 599                if len(image) % 2 == 0:
 600                    diff = image[1::2] - image[::2]
 601                    new_weights = weights[::2] + weights[1::2]
 602                    smooth = image[::2] + (weights[1::2] / new_weights) * diff
 603
 604                else:
 605                    tmp_img = image[:-1]
 606                    diff = tmp_img[1::2] - tmp_img[::2]
 607
 608                    tmp_weight = weights[:-1]
 609                    new_weights = tmp_weight[::2] + tmp_weight[1::2]
 610                    new_weights = np.append(new_weights, weights[-1])
 611
 612                    tmp_img = image[:-1]
 613                    smooth = tmp_img[::2] + (weights[1::2] / new_weights[:-1]) * diff
 614                    smooth = np.append(smooth, image[-1])
 615
 616                tot_transf.append(smooth)
 617                tot_weight.append(new_weights)
 618
 619            # weights = new_weights
 620            # image = smooth
 621            count += 1
 622
 623        return tot_weight
 624
 625    def _unif(self, dim):
 626        # calculate the uniform upscaling steps
 627
 628        min_dim = min(dim)  # find the lowest dimension
 629        max_levels = 0
 630        while min_dim/2 >= 1:
 631            min_dim = min_dim/2
 632            max_levels += 1
 633
 634        # conservative choice
 635        trunc_level = int(np.floor(max_levels*self.upscale['maxtrunc']))
 636
 637        # all the initial cells are upscaled
 638        wght = [np.ones(tuple([int(elem) for elem in dim]))]
 639
 640        for i in range(trunc_level):
 641            new_dim = [int(np.ceil(elem/2))
 642                       for elem in wght[i].shape]  # always take a larger dimension
 643            wght.append(np.ones(tuple(new_dim)))
 644
 645        return wght
 646
 647    def _haarTrunc(self, image, weights, max_val, max_diff, merge):
 648        """
 649        Function for truncating the wavelets. Based on the max_val and max_diff values this function set the detail
 650        coaffiecient to zero if the value of this coefficient is below max_diff, and the value of the smooth coefficient
 651        is below max_val.
 652
 653        Parameters
 654        ----------
 655        image : array-like
 656            The transformed image.
 657
 658        weights : array-like
 659            The weights of the transformed image.
 660
 661        max_val : float
 662            Smooth values above this value are not allowed.
 663
 664        max_diff : float
 665            Detail coefficients higher than this value are not truncated.
 666
 667        merge : array-like
 668            Matrix/vector of booleans defining whether a merge is allowed.
 669
 670        Returns
 671        -------
 672        image : array-like
 673            Transformed image with truncated coefficients.
 674
 675        weights : array-like
 676            Updated weights.
 677
 678        allow_merge : array-like
 679            Booleans keeping track of allowed merges.
 680
 681        end_ctrl : bool
 682            Boolean to control whether further upscaling is possible.
 683        """
 684        # the image will always contain even elements in x and y dir
 685        x_dir = int(image.shape[0] / 2)
 686        y_dir = int(image.shape[1] / 2)
 687
 688        smooth = image[:x_dir, :y_dir]
 689        smooth_diff = image[x_dir:, :y_dir]
 690        diff_smooth = image[:x_dir, y_dir:]
 691        diff_diff = image[x_dir:, y_dir:]
 692
 693        weights_diff = weights[x_dir:, :y_dir]
 694
 695        allow_merge = np.array(
 696            [[False] * y_dir] * x_dir)
 697
 698        # ensure that last row does not be left behind during non-dyadic upscaling
 699        if image.shape[0] > merge.shape[0]:
 700            merge = np.insert(merge, -1, merge[-1, :], axis=0)
 701        if image.shape[1] > merge.shape[1]:
 702            merge = np.insert(merge, -1, merge[:, -1], axis=1)
 703
 704        cand1 = zip(*np.where(smooth < max_val))
 705
 706        end_ctrl = True
 707
 708        for elem in cand1:
 709            # If the method wants to merge cells outside the box the if sentence gives an error, make an exception for
 710            # this hence do not merge these cells.
 711            try:
 712                if abs(smooth_diff[elem]) <= max_diff and abs(diff_smooth[elem]) <= max_diff and \
 713                        merge[elem[0] * 2, elem[1] * 2] and merge[elem[0] * 2 + 1, elem[1] * 2] and \
 714                        merge[elem[0] * 2 + 1, elem[1] * 2 + 1] and merge[elem[0] * 2, elem[1] * 2 + 1]:
 715                    diff_smooth[elem] = 0
 716                    smooth_diff[elem] = 0
 717                    diff_diff[elem] = 0
 718                    weights_diff[elem] = True
 719                    allow_merge[elem] = True
 720                    end_ctrl = False
 721            except:
 722                pass
 723        # for elem in cand2:
 724        #     if smooth[elem] <= max_val:
 725        #         diff_smooth[elem] = 0
 726        #         smooth_diff[elem] = 0
 727        #         diff_diff[elem] = 0
 728
 729        image[x_dir:, :y_dir] = smooth_diff
 730        image[:x_dir, y_dir:] = diff_smooth
 731        image[x_dir:, y_dir:] = diff_diff
 732
 733        weights[x_dir:, :y_dir] = weights_diff
 734
 735        return image, weights, allow_merge, end_ctrl
 736
 737    def _runMako(self, folder, state):
 738        """
 739        Read the mako template (.mako) file from ../folder, and render the correct data file (.DATA) in folder.
 740
 741        Parameters
 742        ----------
 743        Folder : str, optional
 744            Folder for the ecl_100 run.
 745
 746        ensembleMember : int, optional
 747            Index of the ensemble member to run.
 748
 749        Changelog
 750        ---------
 751        - KF 14/9-2015
 752        """
 753
 754        # Check and add report time
 755        if hasattr(self, 'report'):
 756            for key in self.report:
 757                state[key] = self.report[key]
 758
 759        # Convert drilling order (float) to drilling queue (integer) - drilling order optimization
 760        # if "drillingorder" in en_dict:
 761        #     dorder = en_dict['drillingorder']
 762        #     en_dict['drillingqueue'] = np.argsort(dorder)[::-1][:len(dorder)]
 763
 764        # Add startdate
 765        if hasattr(self, 'startDate'):
 766            state['startdate'] = dt.datetime(
 767                self.startDate['year'], self.startDate['month'], self.startDate['day'])
 768
 769        # Add the coarsing values
 770        if hasattr(self, 'coarse'):
 771            state['coarse'] = self.coarse
 772
 773        # Look for the mako file
 774        lkup = TemplateLookup(directories=os.getcwd(),
 775                              input_encoding='utf-8')
 776
 777        # Get template
 778        # If we need the E300 run, define a E300 data file with _E300 added to the end.
 779        if hasattr(self, 'E300'):
 780            if self.E300:
 781                tmpl = lkup.get_template('%s.mako' % (self.file + '_E300'))
 782            else:
 783                tmpl = lkup.get_template('%s.mako' % self.file)
 784        else:
 785            tmpl = lkup.get_template('%s.mako' % self.file)
 786
 787        # use a context and render onto a file
 788        with open('{0}{1}'.format(folder + self.file, '.DATA'), 'w') as f:
 789            ctx = Context(f, **state)
 790            tmpl.render_context(ctx)
 791
 792    def get_sim_results(self, whichResponse, ext_data_info=None, member=None):
 793        """
 794        Read the output from simulator and store as an array. Optionally, if the DA method is based on an ensemble
 795        method, the output must be read inside a folder.
 796
 797        Parameters
 798        ----------
 799        file_rsm : str
 800            Summary file from ECL 100.
 801
 802        file_rst : str
 803            Restart file from ECL 100.
 804
 805        whichResponse : str
 806            Which of the responses is to be outputted (e.g., WBHP PRO-1, WOPR, PRESS, OILSAT, etc).
 807
 808        member : int, optional
 809            Ensemble member that is finished.
 810
 811        ext_data_info : tuple, optional
 812            Tuple containing the assimilation step information, including the place of assimilation (e.g., which TIME) and the
 813            index of this assimilation place.
 814
 815        Returns
 816        -------
 817        yFlow : array-like
 818            Array containing the response from ECL 100. The response type is chosen by the user in options['data_type'].
 819
 820        Notes
 821        -----
 822        - Modified the ecl package to allow reading the summary data directly, hence, we get cell, summary, and field data
 823        from the ecl package.
 824        - KF 29/10-2015
 825        - Modified the ecl package to read RFT files. To obtain, e.g,. RFT pressures form well 'PRO-1" whichResponse
 826        must be rft_PRESSURE PRO-1
 827        """
 828        # Check that we have no trailing spaces
 829        whichResponse = whichResponse.strip()
 830
 831        # if ensemble DA method
 832        if member is not None:
 833            # Get results
 834            if hasattr(self, 'ecl_case'):
 835                # En_XX/YYYY.DATA is the folder setup
 836                rt_mem = int(self.ecl_case.root.split('/')[0].split('_')[1])
 837                if rt_mem != member:  # wrong case
 838                    self.ecl_case = ecl.EclipseCase('En_' + str(member) + os.sep + self.file + '.DATA')
 839            else:
 840                self.ecl_case = ecl.EclipseCase('En_' + str(member) + os.sep + self.file + '.DATA')
 841            if ext_data_info[0] == 'days':
 842                time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \
 843                    dt.timedelta(days=ext_data_info[1])
 844                dates = self.ecl_case.by_date
 845                if time not in dates and 'date_slack' in self.input_dict:
 846                    slack = int(self.input_dict['date_slack'])
 847                    if slack > 0:
 848                        v = [el for el in dates if np.abs(
 849                            (el-time).total_seconds()) < slack]
 850                        if len(v) > 0:
 851                            time = v[0]
 852            else:
 853                time = ext_data_info[1]
 854
 855            # Check if the data is a field or well data, by checking if the well is defined
 856            if len(whichResponse.split(' ')) == 2:
 857                # if rft, search for rft_
 858                if 'rft_' in whichResponse:
 859                    # to speed up when performing the prediction step
 860                    if hasattr(self, 'rft_case'):
 861                        rt_mem = int(self.rft_case.root.split('/')[0].split('_')[1])
 862                        if rt_mem != member:
 863                            self.rft_case = ecl.EclipseRFT('En_' + str(member) + os.sep + self.file)
 864                    else:
 865                        self.rft_case = ecl.EclipseRFT(
 866                            'En_' + str(member) + os.sep + self.file)
 867                    # Get the data. Due to formating we can slice the property.
 868                    rft_prop = self.rft_case.rft_data(well=whichResponse.split(
 869                        ' ')[1], prop=whichResponse.split(' ')[0][4:])
 870                    # rft_data are collected for open connections. This may vary throughout the simulation, hence we
 871                    # must also collect the depth for the rft_data to check if all data is present
 872                    rft_depth = self.rft_case.rft_data(
 873                        well=whichResponse.split(' ')[1], prop='DEPTH')
 874                    # to check this we import the referance depth if this is available. If not we assume that the data
 875                    # is ok.
 876                    try:
 877                        ref_depth_f = np.load(whichResponse.split(
 878                            ' ')[1].upper() + '_rft_ref_depth.npz')
 879                        ref_depth = ref_depth_f[ref_depth_f.files[0]]
 880                        yFlow = np.array([])
 881                        interp = interpolate.interp1d(rft_depth, rft_prop, kind='linear', bounds_error=False,
 882                                                      fill_value=(rft_prop[0], rft_prop[-1]))
 883                        for d in ref_depth:
 884                            yFlow = np.append(yFlow, interp(d))
 885                    except:
 886                        yFlow = rft_prop
 887                else:
 888                    # If well, read the rsm file
 889                    if ext_data_info is not None:  # Get the data at a specific well and time
 890                        yFlow = self.ecl_case.summary_data(whichResponse, time)
 891            elif len(whichResponse.split(' ')) == 1:  # field data
 892                if whichResponse.upper() in ['FOPT', 'FWPT', 'FGPT', 'FWIT', 'FGIT']:
 893                    if ext_data_info is not None:
 894                        yFlow = self.ecl_case.summary_data(whichResponse, time)
 895                elif whichResponse.upper() in ['PERMX', 'PERMY', 'PERMZ', 'PORO', 'NTG', 'SATNUM',
 896                                               'MULTNUM', 'OPERNUM']:
 897                    yFlow = np.array([self.ecl_case.cell_data(whichResponse).flatten()[time]]) # assume that time is the index
 898                else:
 899                    yFlow = self.ecl_case.cell_data(whichResponse, time).flatten()
 900                    if yFlow is None:
 901                        yFlow = self.ecl_case.cell_data(whichResponse).flatten()
 902
 903            # store the run time. NB: elapsed must be defined in .DATA file for this to work
 904            if 'save_elapsed' in self.input_dict and len(self.run_time) <= member:
 905                self.run_time.extend(self.ecl_case.summary_data('ELAPSED', time))
 906
 907            # If we have performed coarsening, we store the number of active grid-cells
 908            if self.upscale is not None:
 909                # Get this number from INIT file
 910                with ecl.EclipseFile('En_' + str(member) + os.sep + self.file, 'INIT') as case:
 911                    intHead = case.get('INTEHEAD')
 912                # The active cell is element 12 of this vector, index 11 in python indexing...
 913                active_cells = intHead[11]
 914                if len(self.num_act) <= member:
 915                    self.num_act.extend([active_cells])
 916
 917        else:
 918            case = ecl.EclipseCase(self.file + '.DATA')
 919            if ext_data_info[0] == 'days':
 920                time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \
 921                    dt.timedelta(days=ext_data_info[1])
 922            else:
 923                time = ext_data_info[1]
 924
 925            # Check if the data is a field or well data, by checking if the well is defined
 926            if len(whichResponse.split(' ')) == 2:
 927                # if rft, search for rft_
 928                if 'rft_' in whichResponse:
 929                    rft_case = ecl.EclipseRFT(self.file)
 930                    # Get the data. Due to formating we can slice the property.
 931                    rft_prop = rft_case.rft_data(well=whichResponse.split(
 932                        ' ')[1], prop=whichResponse.split(' ')[0][4:])
 933                    # rft_data are collected for open connections. This may vary throughout the simulation, hence we
 934                    # must also collect the depth for the rft_data to check if all data is present
 935                    rft_depth = rft_case.rft_data(
 936                        well=whichResponse.split(' ')[1], prop='DEPTH')
 937                    try:
 938                        ref_depth_f = np.load(whichResponse.split(
 939                            ' ')[1].upper() + '_rft_ref_depth.npz')
 940                        ref_depth = ref_depth_f[ref_depth_f.files[0]]
 941                        yFlow = np.array([])
 942                        interp = interpolate.interp1d(rft_depth, rft_prop, kind='linear', bounds_error=False,
 943                                                      fill_value=(rft_prop[0], rft_prop[-1]))
 944                        for d in ref_depth:
 945                            yFlow = np.append(yFlow, interp(d))
 946                    except:
 947                        yFlow = rft_prop
 948                else:
 949                    # If well, read the rsm file
 950                    if ext_data_info is not None:  # Get the data at a specific well and time
 951                        yFlow = case.summary_data(whichResponse, time)
 952            elif len(whichResponse.split(' ')) == 1:
 953                if whichResponse in ['FOPT', 'FWPT', 'FGPT', 'FWIT', 'FGIT']:
 954                    if ext_data_info is not None:
 955                        yFlow = case.summary_data(whichResponse, time)
 956                else:
 957                    yFlow = case.cell_data(whichResponse, time).flatten()
 958                    if yFlow is None:
 959                        yFlow = case.cell_data(whichResponse).flatten()
 960
 961            # If we have performed coarsening, we store the number of active grid-cells
 962            if self.upscale is not None:
 963                # Get this number from INIT file
 964                with ecl.EclipseFile('En_' + str(member) + os.sep + self.file,'INIT') as case:
 965                    intHead = case.get('INTEHEAD')
 966                # The active cell is element 12 of this vector, index 11 in python indexing...
 967                active_cells = intHead[11]
 968                if len(self.num_act) <= member:
 969                    self.num_act.extend([active_cells])
 970
 971        return yFlow
 972
 973    def store_fwd_debug(self, assimstep):
 974        if 'fwddebug' in self.keys_fwd:
 975            # Init dict. of variables to save
 976            save_dict = {}
 977
 978            # Make sure "ANALYSISDEBUG" gives a list
 979            if isinstance(self.keys_fwd['fwddebug'], list):
 980                debug = self.keys_fwd['fwddebug']
 981            else:
 982                debug = [self.keys_fwd['fwddebug']]
 983
 984            # Loop over variables to store in save list
 985            for var in debug:
 986                # Save with key equal variable name and the actual variable
 987                if isinstance(eval('self.' + var), dict):
 988                    # save directly
 989                    np.savez('fwd_debug_%s_%i' % (var, assimstep), **eval('self.' + var))
 990                else:
 991                    save_dict[var] = eval('self.' + var)
 992
 993            np.savez('fwd_debug_%i' % (assimstep), **save_dict)
 994
 995    def write_to_grid(self, value, propname, path, dim, t_ind=None):
 996        if t_ind == None:
 997            trans_dict = {}
 998
 999            def _lookup(kw):
1000                return trans_dict[kw] if kw in trans_dict else kw
1001
1002            # Write a quantity to the grid as a grdecl file
1003            with open(path + propname + '.grdecl', 'wb') as fileobj:
1004                grdecl._write_kw(fileobj, propname, value, _lookup, dim)
1005        else:
1006            pass
1007            # some errors with rips
1008            # p = Process(target=_write_to_resinsight, args=(list(value[~value.mask]),propname, t_ind))
1009            # Find an open resinsight case
1010            # p.start()
1011            # time.sleep(1)
1012            # p.terminate()
1013
1014# def _write_to_resinsight(value, name,t_ind):
1015#     resinsight = rips.Instance.find()
1016#     case = resinsight.project.case(case_id=0)
1017#     case.set_active_cell_property(value, 'GENERATED', name,t_ind)
1018
1019
1020class ecl_100(eclipse):
1021    '''
1022    ecl_100 class
1023    '''
1024
1025    def call_sim(self, path=None, wait_for_proc=False):
1026        """
1027        Method for calling the ecl_100 simulator.
1028
1029        Parameters
1030        ----------
1031        Path : str
1032            Alternative folder for the ecl_100.data file.
1033
1034        Wait_for_proc : bool, optional
1035            Logical variable to wait for the simulator to finish. Default is False.
1036
1037        Returns
1038        -------
1039        .RSM : str
1040            Run summary file in the standard ECL format. Well data are collected from this file.
1041
1042        .RST : str
1043            Restart file in the standard ECL format. Pressure and saturation data are collected from this file.
1044
1045        .PRT : str
1046            Info file to be used for checking errors in the run.
1047
1048
1049        Changelog
1050        ---------
1051        - KF 14/9-2015
1052        """
1053
1054        # Filename
1055        if path is not None:
1056            filename = path + self.file
1057        else:
1058            filename = self.file
1059
1060        # Run the simulator:
1061        success = True
1062        try:
1063            with EclipseRunEnvironment(filename):
1064                com = ['eclrun', '--nocleanup', 'eclipse', filename + '.DATA']
1065                if 'sim_limit' in self.options:
1066                    call(com, stdout=DEVNULL, timeout=self.options['sim_limit'])
1067                else:
1068                    call(com, stdout=DEVNULL)
1069                raise ValueError
1070        except:
1071            print('\nError in the eclipse run.')  # add rerun?
1072            if not os.path.exists('Crashdump'):
1073                copytree(path, 'Crashdump')
1074            success = False
1075
1076        return success
1077
1078
1079class ecl_300(eclipse):
1080    '''
1081    eclipse 300 class
1082    '''
1083
1084    def call_sim(self, path=None, wait_for_proc=False):
1085        """
1086        Method for calling the ecl_300 simulator.
1087
1088        Parameters
1089        ----------
1090        Path : str
1091            Alternative folder for the ecl_100.data file.
1092
1093        Wait_for_proc : bool, optional
1094            Logical variable to wait for the simulator to finish. Default is False.
1095
1096        NOTE: For now, this option is only utilized in a single localization option.
1097
1098        Returns
1099        -------
1100        RSM : str
1101            Run summary file in the standard ECL format. Well data are collected from this file.
1102
1103        RST : str
1104            Restart file in the standard ECL format. Pressure and saturation data are collected from this file.
1105
1106        PRT : str
1107            Info file to be used for checking errors in the run.
1108
1109        Changelog
1110        ---------
1111        - KF 8/10-2015
1112
1113
1114        """
1115        # Filename
1116        if path is not None:
1117            filename = path + self.file
1118        else:
1119            filename = self.file
1120
1121        # Run the simulator:
1122        with EclipseRunEnvironment(filename):
1123            call(['eclrun', '--nocleanup', 'e300', filename + '.DATA'], stdout=DEVNULL)