simulator.flow_rock

  1from simulator.opm import flow
  2from importlib import import_module
  3import datetime as dt
  4import numpy as np
  5import os
  6from misc import ecl, grdecl
  7import shutil
  8import glob
  9from subprocess import Popen, PIPE
 10import mat73
 11from copy import deepcopy
 12from sklearn.cluster import KMeans
 13from sklearn.preprocessing import StandardScaler
 14
 15
 16class flow_sim2seis(flow):
 17    """
 18    Couple the OPM-flow simulator with a sim2seis simulator such that both reservoir quantities and petro-elastic
 19    quantities can be calculated. Inherit the flow class, and use super to call similar functions.
 20    """
 21
 22    def __init__(self, input_dict=None, filename=None, options=None):
 23        super().__init__(input_dict, filename, options)
 24        self._getpeminfo(input_dict)
 25
 26        self.dum_file_root = 'dummy.txt'
 27        self.dum_entry = str(0)
 28        self.date_slack = None
 29        if 'date_slack' in input_dict:
 30            self.date_slack = int(input_dict['date_slack'])
 31
 32        # If we want to extract, or evaluate, something uniquely from the ensemble specific run we can
 33        # run a user defined code to do this.
 34        self.saveinfo = None
 35        if 'savesiminfo' in input_dict:
 36            # Make sure "ANALYSISDEBUG" gives a list
 37            if isinstance(input_dict['savesiminfo'], list):
 38                self.saveinfo = input_dict['savesiminfo']
 39            else:
 40                self.saveinfo = [input_dict['savesiminfo']]
 41
 42        self.scale = []
 43
 44    def _getpeminfo(self, input_dict):
 45        """
 46        Get, and return, flow and PEM modules
 47        """
 48        if 'pem' in input_dict:
 49            self.pem_input = {}
 50            for elem in input_dict['pem']:
 51                if elem[0] == 'model':  # Set the petro-elastic model
 52                    self.pem_input['model'] = elem[1]
 53                if elem[0] == 'depth':  # provide the npz of depth values
 54                    self.pem_input['depth'] = elem[1]
 55                if elem[0] == 'actnum':  # the npz of actnum values
 56                    self.pem_input['actnum'] = elem[1]
 57                if elem[0] == 'baseline':  # the time for the baseline 4D measurment
 58                    self.pem_input['baseline'] = elem[1]
 59                if elem[0] == 'vintage':
 60                    self.pem_input['vintage'] = elem[1]
 61                    if not type(self.pem_input['vintage']) == list:
 62                        self.pem_input['vintage'] = [elem[1]]
 63                if elem[0] == 'ntg':
 64                    self.pem_input['ntg'] = elem[1]
 65                if elem[0] == 'press_conv':
 66                    self.pem_input['press_conv'] = elem[1]
 67                if elem[0] == 'compaction':
 68                    self.pem_input['compaction'] = True
 69                if elem[0] == 'overburden':  # the npz of overburden values
 70                    self.pem_input['overburden'] = elem[1]
 71                if elem[0] == 'percentile':  # use for scaling
 72                    self.pem_input['percentile'] = elem[1]
 73
 74            pem = getattr(import_module('simulator.rockphysics.' +
 75                          self.pem_input['model'].split()[0]), self.pem_input['model'].split()[1])
 76
 77            self.pem = pem(self.pem_input)
 78
 79        else:
 80            self.pem = None
 81
 82    def setup_fwd_run(self):
 83        super().setup_fwd_run()
 84
 85    def run_fwd_sim(self, state, member_i, del_folder=True):
 86        # The inherited simulator also has a run_fwd_sim. Call this.
 87        self.ensemble_member = member_i
 88        self.pred_data = super().run_fwd_sim(state, member_i, del_folder=True)
 89
 90        return self.pred_data
 91
 92    def call_sim(self, folder=None, wait_for_proc=False):
 93        # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first.
 94        # Then, get the pem.
 95        success = super().call_sim(folder, wait_for_proc)
 96
 97        if success:
 98            # need a if to check that we have correct sim2seis
 99            # copy relevant sim2seis files into folder.
100            for file in glob.glob('sim2seis_config/*'):
101                shutil.copy(file, 'En_' + str(self.ensemble_member) + os.sep)
102
103            self.ecl_case = ecl.EclipseCase(
104                'En_' + str(self.ensemble_member) + os.sep + self.file + '.DATA')
105            grid = self.ecl_case.grid()
106
107            phases = self.ecl_case.init.phases
108            if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases:  # This should be extended
109                vintage = []
110                # loop over seismic vintages
111                for v, assim_time in enumerate(self.pem_input['vintage']):
112                    time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \
113                        dt.timedelta(days=assim_time)
114                    pem_input = {}
115                    # get active porosity
116                    tmp = self.ecl_case.cell_data('PORO')
117                    if 'compaction' in self.pem_input:
118                        multfactor = self.ecl_case.cell_data('PORV_RC', time)
119
120                        pem_input['PORO'] = np.array(
121                            multfactor[~tmp.mask]*tmp[~tmp.mask], dtype=float)
122                    else:
123                        pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float)
124                    # get active NTG if needed
125                    if 'ntg' in self.pem_input:
126                        if self.pem_input['ntg'] == 'no':
127                            pem_input['NTG'] = None
128                        else:
129                            tmp = self.ecl_case.cell_data('NTG')
130                            pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float)
131                    else:
132                        tmp = self.ecl_case.cell_data('NTG')
133                        pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float)
134
135                    for var in ['SWAT', 'SGAS', 'PRESSURE', 'RS']:
136                        tmp = self.ecl_case.cell_data(var, time)
137                        # only active, and conv. to float
138                        pem_input[var] = np.array(tmp[~tmp.mask], dtype=float)
139
140                    if 'press_conv' in self.pem_input:
141                        pem_input['PRESSURE'] = pem_input['PRESSURE'] * \
142                            self.pem_input['press_conv']
143
144                    tmp = self.ecl_case.cell_data('PRESSURE', 1)
145                    if hasattr(self.pem, 'p_init'):
146                        P_init = self.pem.p_init*np.ones(tmp.shape)[~tmp.mask]
147                    else:
148                        # initial pressure is first
149                        P_init = np.array(tmp[~tmp.mask], dtype=float)
150
151                    if 'press_conv' in self.pem_input:
152                        P_init = P_init*self.pem_input['press_conv']
153
154                    saturations = [1 - (pem_input['SWAT'] + pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)]
155                                   for ph in phases]
156                    # Get the pressure
157                    self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'],
158                                        ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init,
159                                        ensembleMember=self.ensemble_member)
160
161                    grdecl.write(f'En_{str(self.ensemble_member)}/Vs{v+1}.grdecl', {
162                                 'Vs': self.pem.getShearVel()*.1, 'DIMENS': grid['DIMENS']}, multi_file=False)
163                    grdecl.write(f'En_{str(self.ensemble_member)}/Vp{v+1}.grdecl', {
164                                 'Vp': self.pem.getBulkVel()*.1, 'DIMENS': grid['DIMENS']}, multi_file=False)
165                    grdecl.write(f'En_{str(self.ensemble_member)}/rho{v+1}.grdecl',
166                                 {'rho': self.pem.getDens(), 'DIMENS': grid['DIMENS']}, multi_file=False)
167
168            current_folder = os.getcwd()
169            run_folder = current_folder + os.sep + 'En_' + str(self.ensemble_member)
170            # The sim2seis is invoked via a shell script. The simulations provides outputs. Run, and get all output. Search
171            # for Done. If not finished in reasonable time -> kill
172            p = Popen(['./sim2seis.sh', run_folder], stdout=PIPE)
173            start = time
174            while b'done' not in p.stdout.readline():
175                pass
176
177            # Todo: handle sim2seis or pem error
178
179        return success
180
181    def extract_data(self, member):
182        # start by getting the data from the flow simulator
183        super().extract_data(member)
184
185        # get the sim2seis from file
186        for prim_ind in self.l_prim:
187            # Loop over all keys in pred_data (all data types)
188            for key in self.all_data_types:
189                if key in ['sim2seis']:
190                    if self.true_prim[1][prim_ind] in self.pem_input['vintage']:
191                        result = mat73.loadmat(f'En_{member}/Data_conv.mat')['data_conv']
192                        v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind])
193                        self.pred_data[prim_ind][key] = np.sum(
194                            np.abs(result[:, :, :, v]), axis=0).flatten()
195
196class flow_rock(flow):
197    """
198    Couple the OPM-flow simulator with a rock-physics simulator such that both reservoir quantities and petro-elastic
199    quantities can be calculated. Inherit the flow class, and use super to call similar functions.
200    """
201
202    def __init__(self, input_dict=None, filename=None, options=None):
203        super().__init__(input_dict, filename, options)
204        self._getpeminfo(input_dict)
205
206        self.dum_file_root = 'dummy.txt'
207        self.dum_entry = str(0)
208        self.date_slack = None
209        if 'date_slack' in input_dict:
210            self.date_slack = int(input_dict['date_slack'])
211
212        # If we want to extract, or evaluate, something uniquely from the ensemble specific run we can
213        # run a user defined code to do this.
214        self.saveinfo = None
215        if 'savesiminfo' in input_dict:
216            # Make sure "ANALYSISDEBUG" gives a list
217            if isinstance(input_dict['savesiminfo'], list):
218                self.saveinfo = input_dict['savesiminfo']
219            else:
220                self.saveinfo = [input_dict['savesiminfo']]
221
222        self.scale = []
223
224    def _getpeminfo(self, input_dict):
225        """
226        Get, and return, flow and PEM modules
227        """
228        if 'pem' in input_dict:
229            self.pem_input = {}
230            for elem in input_dict['pem']:
231                if elem[0] == 'model':  # Set the petro-elastic model
232                    self.pem_input['model'] = elem[1]
233                if elem[0] == 'depth':  # provide the npz of depth values
234                    self.pem_input['depth'] = elem[1]
235                if elem[0] == 'actnum':  # the npz of actnum values
236                    self.pem_input['actnum'] = elem[1]
237                if elem[0] == 'baseline':  # the time for the baseline 4D measurment
238                    self.pem_input['baseline'] = elem[1]
239                if elem[0] == 'vintage':
240                    self.pem_input['vintage'] = elem[1]
241                    if not type(self.pem_input['vintage']) == list:
242                        self.pem_input['vintage'] = [elem[1]]
243                if elem[0] == 'ntg':
244                    self.pem_input['ntg'] = elem[1]
245                if elem[0] == 'press_conv':
246                    self.pem_input['press_conv'] = elem[1]
247                if elem[0] == 'compaction':
248                    self.pem_input['compaction'] = True
249                if elem[0] == 'overburden':  # the npz of overburden values
250                    self.pem_input['overburden'] = elem[1]
251                if elem[0] == 'percentile':  # use for scaling
252                    self.pem_input['percentile'] = elem[1]
253
254            pem = getattr(import_module('simulator.rockphysics.' +
255                          self.pem_input['model'].split()[0]), self.pem_input['model'].split()[1])
256
257            self.pem = pem(self.pem_input)
258
259        else:
260            self.pem = None
261
262    def setup_fwd_run(self, redund_sim):
263        super().setup_fwd_run(redund_sim=redund_sim)
264
265    def run_fwd_sim(self, state, member_i, del_folder=True):
266        # The inherited simulator also has a run_fwd_sim. Call this.
267        self.ensemble_member = member_i
268        self.pred_data = super().run_fwd_sim(state, member_i, del_folder=True)
269
270        return self.pred_data
271
272    def call_sim(self, folder=None, wait_for_proc=False):
273        # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first.
274        # Then, get the pem.
275        success = super().call_sim(folder, wait_for_proc)
276
277        if success:
278            self.ecl_case = ecl.EclipseCase(
279                'En_' + str(self.ensemble_member) + os.sep + self.file + '.DATA')
280            phases = self.ecl_case.init.phases
281            #if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases:  # This should be extended
282            if 'WAT' in phases and 'GAS' in phases:
283                vintage = []
284                # loop over seismic vintages
285                for v, assim_time in enumerate(self.pem_input['vintage']):
286                    time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \
287                        dt.timedelta(days=assim_time)
288                    pem_input = {}
289                    # get active porosity
290                    tmp = self.ecl_case.cell_data('PORO')
291                    if 'compaction' in self.pem_input:
292                        multfactor = self.ecl_case.cell_data('PORV_RC', time)
293
294                        pem_input['PORO'] = np.array(
295                            multfactor[~tmp.mask]*tmp[~tmp.mask], dtype=float)
296                    else:
297                        pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float)
298                    # get active NTG if needed
299                    if 'ntg' in self.pem_input:
300                        if self.pem_input['ntg'] == 'no':
301                            pem_input['NTG'] = None
302                        else:
303                            tmp = self.ecl_case.cell_data('NTG')
304                            pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float)
305                    else:
306                        tmp = self.ecl_case.cell_data('NTG')
307                        pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float)
308
309                    pem_input['RS'] = None
310                    for var in ['SWAT', 'SGAS', 'PRESSURE', 'RS']:
311                        try:
312                            tmp = self.ecl_case.cell_data(var, time)
313                        except:
314                            pass
315                        # only active, and conv. to float
316                        pem_input[var] = np.array(tmp[~tmp.mask], dtype=float)
317
318                    if 'press_conv' in self.pem_input:
319                        pem_input['PRESSURE'] = pem_input['PRESSURE'] * \
320                            self.pem_input['press_conv']
321
322                    tmp = self.ecl_case.cell_data('PRESSURE', 1)
323                    if hasattr(self.pem, 'p_init'):
324                        P_init = self.pem.p_init*np.ones(tmp.shape)[~tmp.mask]
325                    else:
326                        # initial pressure is first
327                        P_init = np.array(tmp[~tmp.mask], dtype=float)
328
329                    if 'press_conv' in self.pem_input:
330                        P_init = P_init*self.pem_input['press_conv']
331
332                    saturations = [1 - (pem_input['SWAT'] + pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)]
333                                   for ph in phases]
334                    # Get the pressure
335                    self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'],
336                                        ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init,
337                                        ensembleMember=self.ensemble_member)
338                    # mask the bulkimp to get proper dimensions
339                    tmp_value = np.zeros(self.ecl_case.init.shape)
340                    tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp
341                    self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float,
342                                                   mask=deepcopy(self.ecl_case.init.mask))
343                    # run filter
344                    self.pem._filter()
345                    vintage.append(deepcopy(self.pem.bulkimp))
346
347            if hasattr(self.pem, 'baseline'):  # 4D measurement
348                base_time = dt.datetime(self.startDate['year'], self.startDate['month'],
349                                        self.startDate['day']) + dt.timedelta(days=self.pem.baseline)
350                # pem_input = {}
351                # get active porosity
352                tmp = self.ecl_case.cell_data('PORO')
353
354                if 'compaction' in self.pem_input:
355                    multfactor = self.ecl_case.cell_data('PORV_RC', base_time)
356
357                    pem_input['PORO'] = np.array(
358                        multfactor[~tmp.mask] * tmp[~tmp.mask], dtype=float)
359                else:
360                    pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float)
361
362                pem_input['RS'] = None
363                for var in ['SWAT', 'SGAS', 'PRESSURE', 'RS']:
364                    try:
365                        tmp = self.ecl_case.cell_data(var, base_time)
366                    except:
367                        pass
368                    # only active, and conv. to float
369                    pem_input[var] = np.array(tmp[~tmp.mask], dtype=float)
370
371                if 'press_conv' in self.pem_input:
372                    pem_input['PRESSURE'] = pem_input['PRESSURE'] * \
373                        self.pem_input['press_conv']
374
375                saturations = [1 - (pem_input['SWAT'] + pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)]
376                               for ph in phases]
377                # Get the pressure
378                self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'],
379                                    ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init,
380                                    ensembleMember=None)
381
382                # mask the bulkimp to get proper dimensions
383                tmp_value = np.zeros(self.ecl_case.init.shape)
384
385                tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp
386                # kill if values are inf or nan
387                assert not np.isnan(tmp_value).any()
388                assert not np.isinf(tmp_value).any()
389                self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float,
390                                               mask=deepcopy(self.ecl_case.init.mask))
391                self.pem._filter()
392
393                # 4D response
394                self.pem_result = []
395                for i, elem in enumerate(vintage):
396                    self.pem_result.append(elem - deepcopy(self.pem.bulkimp))
397            else:
398                for i, elem in enumerate(vintage):
399                    self.pem_result.append(elem)
400
401        return success
402
403    def extract_data(self, member):
404        # start by getting the data from the flow simulator
405        super().extract_data(member)
406
407        # get the sim2seis from file
408        for prim_ind in self.l_prim:
409            # Loop over all keys in pred_data (all data types)
410            for key in self.all_data_types:
411                if key in ['bulkimp']:
412                    if self.true_prim[1][prim_ind] in self.pem_input['vintage']:
413                        v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind])
414                        self.pred_data[prim_ind][key] = self.pem_result[v].data.flatten()
415
416class flow_barycenter(flow):
417    """
418    Couple the OPM-flow simulator with a rock-physics simulator such that both reservoir quantities and petro-elastic
419    quantities can be calculated. Inherit the flow class, and use super to call similar functions. In the end, the
420    barycenter and moment of interia for the bulkimpedance objects, are returned as observations. The objects are
421    identified using k-means clustering, and the number of objects are determined using and elbow strategy.
422    """
423
424    def __init__(self, input_dict=None, filename=None, options=None):
425        super().__init__(input_dict, filename, options)
426        self._getpeminfo(input_dict)
427
428        self.dum_file_root = 'dummy.txt'
429        self.dum_entry = str(0)
430        self.date_slack = None
431        if 'date_slack' in input_dict:
432            self.date_slack = int(input_dict['date_slack'])
433
434        # If we want to extract, or evaluate, something uniquely from the ensemble specific run we can
435        # run a user defined code to do this.
436        self.saveinfo = None
437        if 'savesiminfo' in input_dict:
438            # Make sure "ANALYSISDEBUG" gives a list
439            if isinstance(input_dict['savesiminfo'], list):
440                self.saveinfo = input_dict['savesiminfo']
441            else:
442                self.saveinfo = [input_dict['savesiminfo']]
443
444        self.scale = []
445        self.pem_result = []
446        self.bar_result = []
447
448    def _getpeminfo(self, input_dict):
449        """
450        Get, and return, flow and PEM modules
451        """
452        if 'pem' in input_dict:
453            self.pem_input = {}
454            for elem in input_dict['pem']:
455                if elem[0] == 'model':  # Set the petro-elastic model
456                    self.pem_input['model'] = elem[1]
457                if elem[0] == 'depth':  # provide the npz of depth values
458                    self.pem_input['depth'] = elem[1]
459                if elem[0] == 'actnum':  # the npz of actnum values
460                    self.pem_input['actnum'] = elem[1]
461                if elem[0] == 'baseline':  # the time for the baseline 4D measurment
462                    self.pem_input['baseline'] = elem[1]
463                if elem[0] == 'vintage':
464                    self.pem_input['vintage'] = elem[1]
465                    if not type(self.pem_input['vintage']) == list:
466                        self.pem_input['vintage'] = [elem[1]]
467                if elem[0] == 'ntg':
468                    self.pem_input['ntg'] = elem[1]
469                if elem[0] == 'press_conv':
470                    self.pem_input['press_conv'] = elem[1]
471                if elem[0] == 'compaction':
472                    self.pem_input['compaction'] = True
473                if elem[0] == 'overburden':  # the npz of overburden values
474                    self.pem_input['overburden'] = elem[1]
475                if elem[0] == 'percentile':  # use for scaling
476                    self.pem_input['percentile'] = elem[1]
477                if elem[0] == 'clusters':  # number of clusters for each barycenter
478                    self.pem_input['clusters'] = elem[1]
479
480            pem = getattr(import_module('simulator.rockphysics.' +
481                          self.pem_input['model'].split()[0]), self.pem_input['model'].split()[1])
482
483            self.pem = pem(self.pem_input)
484
485        else:
486            self.pem = None
487
488    def setup_fwd_run(self, redund_sim):
489        super().setup_fwd_run(redund_sim=redund_sim)
490
491    def run_fwd_sim(self, state, member_i, del_folder=True):
492        # The inherited simulator also has a run_fwd_sim. Call this.
493        self.ensemble_member = member_i
494        self.pred_data = super().run_fwd_sim(state, member_i, del_folder=True)
495
496        return self.pred_data
497
498    def call_sim(self, folder=None, wait_for_proc=False):
499        # the super run_fwd_sim will invoke call_sim. Modify this such that the fluid simulator is run first.
500        # Then, get the pem.
501        success = super().call_sim(folder, wait_for_proc)
502
503        if success:
504            self.ecl_case = ecl.EclipseCase(
505                'En_' + str(self.ensemble_member) + os.sep + self.file + '.DATA')
506            phases = self.ecl_case.init.phases
507            #if 'OIL' in phases and 'WAT' in phases and 'GAS' in phases:  # This should be extended
508            if 'WAT' in phases and 'GAS' in phases:
509                vintage = []
510                # loop over seismic vintages
511                for v, assim_time in enumerate(self.pem_input['vintage']):
512                    time = dt.datetime(self.startDate['year'], self.startDate['month'], self.startDate['day']) + \
513                        dt.timedelta(days=assim_time)
514                    pem_input = {}
515                    # get active porosity
516                    tmp = self.ecl_case.cell_data('PORO')
517                    if 'compaction' in self.pem_input:
518                        multfactor = self.ecl_case.cell_data('PORV_RC', time)
519
520                        pem_input['PORO'] = np.array(
521                            multfactor[~tmp.mask]*tmp[~tmp.mask], dtype=float)
522                    else:
523                        pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float)
524                    # get active NTG if needed
525                    if 'ntg' in self.pem_input:
526                        if self.pem_input['ntg'] == 'no':
527                            pem_input['NTG'] = None
528                        else:
529                            tmp = self.ecl_case.cell_data('NTG')
530                            pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float)
531                    else:
532                        tmp = self.ecl_case.cell_data('NTG')
533                        pem_input['NTG'] = np.array(tmp[~tmp.mask], dtype=float)
534
535                    pem_input['RS'] = None
536                    for var in ['SWAT', 'SGAS', 'PRESSURE', 'RS']:
537                        try:
538                            tmp = self.ecl_case.cell_data(var, time)
539                        except:
540                            pass
541                        # only active, and conv. to float
542                        pem_input[var] = np.array(tmp[~tmp.mask], dtype=float)
543
544                    if 'press_conv' in self.pem_input:
545                        pem_input['PRESSURE'] = pem_input['PRESSURE'] * \
546                            self.pem_input['press_conv']
547
548                    tmp = self.ecl_case.cell_data('PRESSURE', 1)
549                    if hasattr(self.pem, 'p_init'):
550                        P_init = self.pem.p_init*np.ones(tmp.shape)[~tmp.mask]
551                    else:
552                        # initial pressure is first
553                        P_init = np.array(tmp[~tmp.mask], dtype=float)
554
555                    if 'press_conv' in self.pem_input:
556                        P_init = P_init*self.pem_input['press_conv']
557
558                    saturations = [1 - (pem_input['SWAT'] + pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)]
559                                   for ph in phases]
560                    # Get the pressure
561                    self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'],
562                                        ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init,
563                                        ensembleMember=self.ensemble_member)
564                    # mask the bulkimp to get proper dimensions
565                    tmp_value = np.zeros(self.ecl_case.init.shape)
566                    tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp
567                    self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float,
568                                                   mask=deepcopy(self.ecl_case.init.mask))
569                    # run filter
570                    self.pem._filter()
571                    vintage.append(deepcopy(self.pem.bulkimp))
572
573            if hasattr(self.pem, 'baseline'):  # 4D measurement
574                base_time = dt.datetime(self.startDate['year'], self.startDate['month'],
575                                        self.startDate['day']) + dt.timedelta(days=self.pem.baseline)
576                # pem_input = {}
577                # get active porosity
578                tmp = self.ecl_case.cell_data('PORO')
579
580                if 'compaction' in self.pem_input:
581                    multfactor = self.ecl_case.cell_data('PORV_RC', base_time)
582
583                    pem_input['PORO'] = np.array(
584                        multfactor[~tmp.mask] * tmp[~tmp.mask], dtype=float)
585                else:
586                    pem_input['PORO'] = np.array(tmp[~tmp.mask], dtype=float)
587
588                pem_input['RS'] = None
589                for var in ['SWAT', 'SGAS', 'PRESSURE', 'RS']:
590                    try:
591                        tmp = self.ecl_case.cell_data(var, base_time)
592                    except:
593                        pass
594                    # only active, and conv. to float
595                    pem_input[var] = np.array(tmp[~tmp.mask], dtype=float)
596
597                if 'press_conv' in self.pem_input:
598                    pem_input['PRESSURE'] = pem_input['PRESSURE'] * \
599                        self.pem_input['press_conv']
600
601                saturations = [1 - (pem_input['SWAT'] + pem_input['SGAS']) if ph == 'OIL' else pem_input['S{}'.format(ph)]
602                               for ph in phases]
603                # Get the pressure
604                self.pem.calc_props(phases, saturations, pem_input['PRESSURE'], pem_input['PORO'],
605                                    ntg=pem_input['NTG'], Rs=pem_input['RS'], press_init=P_init,
606                                    ensembleMember=None)
607
608                # mask the bulkimp to get proper dimensions
609                tmp_value = np.zeros(self.ecl_case.init.shape)
610
611                tmp_value[self.ecl_case.init.actnum] = self.pem.bulkimp
612                # kill if values are inf or nan
613                assert not np.isnan(tmp_value).any()
614                assert not np.isinf(tmp_value).any()
615                self.pem.bulkimp = np.ma.array(data=tmp_value, dtype=float,
616                                               mask=deepcopy(self.ecl_case.init.mask))
617                self.pem._filter()
618
619                # 4D response
620                for i, elem in enumerate(vintage):
621                    self.pem_result.append(elem - deepcopy(self.pem.bulkimp))
622            else:
623                for i, elem in enumerate(vintage):
624                    self.pem_result.append(elem)
625
626            #  Extract k-means centers and interias for each element in pem_result
627            if 'clusters' in self.pem_input:
628                npzfile = np.load(self.pem_input['clusters'], allow_pickle=True)
629                n_clusters_list = npzfile['n_clusters_list']
630                npzfile.close()
631            else:
632                n_clusters_list = len(self.pem_result)*[2]
633            kmeans_kwargs = {"init": "random", "n_init": 10, "max_iter": 300, "random_state": 42}
634            for i, bulkimp in enumerate(self.pem_result):
635                std = np.std(bulkimp)
636                features = np.argwhere(np.squeeze(np.reshape(np.abs(bulkimp), self.ecl_case.init.shape,)) > 3 * std)
637                scaler = StandardScaler()
638                scaled_features = scaler.fit_transform(features)
639                kmeans = KMeans(n_clusters=n_clusters_list[i], **kmeans_kwargs)
640                kmeans.fit(scaled_features)
641                kmeans_center = np.squeeze(scaler.inverse_transform(kmeans.cluster_centers_))  # data / measurements
642                self.bar_result.append(np.append(kmeans_center, kmeans.inertia_))
643
644        return success
645
646    def extract_data(self, member):
647        # start by getting the data from the flow simulator
648        super().extract_data(member)
649
650        # get the barycenters and inertias
651        for prim_ind in self.l_prim:
652            # Loop over all keys in pred_data (all data types)
653            for key in self.all_data_types:
654                if key in ['barycenter']:
655                    if self.true_prim[1][prim_ind] in self.pem_input['vintage']:
656                        v = self.pem_input['vintage'].index(self.true_prim[1][prim_ind])
657                        self.pred_data[prim_ind][key] = self.bar_result[v].flatten()