import numpy as np import os import shutil import re import pandas as pd from tqdm.notebook import tqdm from pytictoc import TicToc import time # Ensemble Kalman Filter module from EnKF import EnKF t = TicToc() t.tic()

NofRealization = 200 dx = dy = dz = 30 nx=25 ny=25 nz=1 ngrid = nx * ny * nz Nofobservedwell = 8

Atime=np.linspace(20,600,30) Ptime=np.linspace(20,1800,90) APtime=np.linspace(620,1800,60) NofAtime=len(Atime)

Modelvariable = ["Permx"] Dynamicvariable = ["WOPR"] NofModelvariable = len(Modelvariable); NofDataType = len(Dynamicvariable) filename = 'SNESim_200En'

EnKF.MakeTstep(Atime[0],'TSTEP.DAT')

for i in range(1,NofRealization): # Assimilation folder if os.path.isdir('EnKF'+str(i)) == False: os.mkdir('EnKF'+str(i)) for j in range(1,NofAtime+1): shutil.copy('WOOGIE.DATA', 'EnKF' + str(i)+ '/' + str(j) + '_WOOGIE.DATA') shutil.copy('TSTEP.DAT', 'EnKF' + str(i) + '/TSTEP.DAT') shutil.copy('SOLUTION.DAT', 'EnKF' + str(i) + '/SOLUTION.DAT') # Initial Ensembels folder if os.path.isdir('ini'+str(i)) == False: os.mkdir('ini'+str(i)) shutil.copy('WOOGIE.DATA', 'ini' + str(i)+ '/WOOGIE.DATA') shutil.copy('refTSTEP.DAT', 'ini' + str(i) + '/TSTEP.DAT') shutil.copy('SOLUTION.DAT', 'ini' + str(i) + '/SOLUTION.DAT') # Reference folder if os.path.isdir('ref') == False: os.mkdir('ref') shutil.copy('WOOGIE.DATA','ref/WOOGIE.DATA') shutil.copy('refTSTEP.DAT','ref/TSTEP.DAT') shutil.copy('SOLUTION.DAT','ref/SOLUTION.DAT')

k_sand = 2000 k_clay = 20 RealizationName ,PermxAll = EnKF.GetParameters_(filename, True, k_sand, k_clay) refPermx = PermxAll['0'] # 일단 임의로 가정 enPermx = PermxAll.iloc[::,1::]

for Realization in RealizationName[1::]: EnKF.MakePermx(Realization ,enPermx, 'EnKF' + Realization +'/PERMX.DAT') EnKF.MakePermx(Realization ,enPermx, 'ini' + Realization +'/PERMX.DAT') EnKF.MakePermx('0' ,PermxAll, 'ref/PERMX.DAT')

EnKF.RunECL('ref','WOOGIE') refWOPR_A = EnKF.GetDynamicAssimilation(Atime, 'ref/WOOGIE') refWOPR_P, refWWCT_P, refFOPT_P, refFWPT_P = EnKF.GetDynamicPrediction(Ptime, 'ref/WOOGIE') refField=np.column_stack((refFOPT_P, refFWPT_P)) OBS = np.array([]) for k in range(1,Nofobservedwell+1): OBS = np.append(OBS,refWOPR_A[(NofDataType * NofAtime) * (k - 1)]) OBS = OBS.astype(np.float64)

WOPR_A =np.array([]) for i in tqdm(range(1,NofRealization), desc = 'Simulate Initial Ensembles'): time.sleep(0.5) EnKF.RunECL('ini' + str(i), 'WOOGIE') if i <= 1: WOPR_A = EnKF.GetDynamicAssimilation(Atime, 'ini1/WOOGIE') else: WOPR_A = np.column_stack((WOPR_A, EnKF.GetDynamicAssimilation(Atime, 'ini' + str(i) + '/WOOGIE')))

FOPT_STEP_Posteriori = np.array([]) FWPT_STEP_Posteriori = np.array([]) cols = ['Real' + str(jj) for jj in range(1,NofRealization)] for i in range(1,NofAtime+1): if i > 1: del WOPR_STEP_Posteriori, WWCT_STEP_Posteriori, FOPT_STEP_Posteriori, FWPT_STEP_Posteriori FOPT_STEP_Posteriori = np.array([]) FWPT_STEP_Posteriori = np.array([]) ''' =================================================================== Prediction Step =================================================================== ''' for j in tqdm(range(1,NofRealization), desc = 'Aissimilation No.{} '.format(i)): if i > 1: EnKF.MakePermx('Real' + str(j), static_vector, 'EnKF' + str(j) + '/PERMX.DAT') EnKF.RunECL('EnKF' + str(j), str(i) + '_WOOGIE') else: EnKF.RunECL('EnKF' + str(j), str(i) + '_WOOGIE') WOPR_A_Posteriori,WWCT_A_Posteriori, FOPT_A_Posteriori, FWPT_A_Posteriori = EnKF.GetDynamicPrediction(Atime[i-1], 'EnKF' + str(j) + '/' + str(i) + '_WOOGIE') if j == 1: WOPR_STEP_Posteriori = WOPR_A_Posteriori WWCT_STEP_Posteriori = WWCT_A_Posteriori else: WOPR_STEP_Posteriori = np.column_stack((WOPR_STEP_Posteriori, WOPR_A_Posteriori)) WWCT_STEP_Posteriori = np.column_stack((WWCT_STEP_Posteriori, WWCT_A_Posteriori)) FOPT_STEP_Posteriori = np.append(FOPT_STEP_Posteriori, FOPT_A_Posteriori) FWPT_STEP_Posteriori = np.append(FWPT_STEP_Posteriori, FWPT_A_Posteriori) WOPR_STEP = pd.DataFrame(data = WOPR_STEP_Posteriori, columns = cols) WWCT_STEP = pd.DataFrame(data = WWCT_STEP_Posteriori, columns = cols) FOPT_STEP = pd.DataFrame(data = FOPT_STEP_Posteriori[np.newaxis,:].reshape(1,-1), columns = cols) FWPT_STEP = pd.DataFrame(data = FWPT_STEP_Posteriori[np.newaxis,:].reshape(1,-1), columns = cols) WOPR_STEP['Label'] = i * np.ones([Nofobservedwell]).astype(np.int_) WWCT_STEP['Label'] = i * np.ones([Nofobservedwell]).astype(np.int_) FOPT_STEP['Label'] = i * np.ones(1).astype(np.int_) FWPT_STEP['Label'] = i * np.ones(1).astype(np.int_) if i == 1: WOPR_STEP_TOTAL = WOPR_STEP WWCT_STEP_TOTAL = WWCT_STEP FOPT_STEP_TOTAL = FOPT_STEP FWPT_STEP_TOTAL = FWPT_STEP else: WOPR_STEP_TOTAL = pd.concat([WOPR_STEP_TOTAL, WOPR_STEP]) WWCT_STEP_TOTAL = pd.concat([WWCT_STEP_TOTAL, WWCT_STEP]) FOPT_STEP_TOTAL = pd.concat([FOPT_STEP_TOTAL, FOPT_STEP]) FWPT_STEP_TOTAL = pd.concat([FWPT_STEP_TOTAL, FWPT_STEP]) ''' =================================================================== Assimilation Step =================================================================== ''' EnKF.MakeTstep(Atime[i-1], 'TSTEP.DAT') if i > 1: EnMean = EnKF.GetMean(static_vector, '../EnKF','WOOGIE') State_vector_Assimilation, En_Permx_Posteriori = EnKF.Ensemble_Kalman_Filter(OBS,En_Permx_Posteriori ,WOPR_STEP_Posteriori, EnMean, NofRealization -1, Nofobservedwell, 0.01) else: EnMean = EnKF.GetMean(enPermx, '../EnKF','WOOGIE') State_vector_Assimilation, En_Permx_Posteriori = EnKF.Ensemble_Kalman_Filter(OBS,enPermx.to_numpy(),WOPR_STEP_Posteriori ,EnMean, NofRealization -1, Nofobservedwell, 0.01) static_vector = pd.DataFrame(data = En_Permx_Posteriori, columns = cols) ''' =================================================================== Preparation Step for next step i + 1 =================================================================== ''' for kk in range(1,NofRealization): if i < NofAtime: EnKF.MakeTstep(Atime[i] - Atime[i-1], '../EnKF/EnKF{}/TSTEP.DATA'.format(kk)) EnKF.MakeSolution(i, 'WOOGIE', '../EnKF/EnKF{}/SOLUTION.DAT'.format(kk)) time_ = t.toc() min_ = int(time_ % 60) hr_ = int(min_ % 60) print('Total Elapsed Time = {} Hr {} min {:.2f} sec'.format(hr_ , min_ - hr_ * 60, time_ - min_ * 60 ))

''' =================================================================== Ensemble Kalman Filter Module Modified by woogie (original : K.B. Lee, 2014 with Matlab) data based on nummpy 원본은 perm 데이터가 vector로 사전 가공되어있음 수정한 본 파일은 SGeMS에서의 파일을 그대로 들고옴 =================================================================== ''' import numpy as np import pandas as pd import subprocess import shutil import os import re class EnKF: ''' =================================================================== [ Make Series ] Assimilation이 진행됨에 따라 순서대로 TESTP.DAT SOLTUTION.DAT PERMX.DAT 을 갱신하는 메소드들임 =================================================================== ''' def MakeTstep(TSTEP,filename): f = open(filename, 'w') f.write('TSTEP\n') f.write('{:s} /\n'.format(str(TSTEP))) f.close() def MakeSolution(ASTEP,filename, solution_filename): with open(solution_filename,'w') as f: f.write('RESTART\n') f.write(f'{ASTEP}_{filename} {ASTEP} /') def MakePermx(key, value, filename): value = np.array(value[key]) with open(filename, 'w') as f: f.write('PERMX\n') for i in range(len(value)): f.write(str(value[i]) + '\n') f.write('/') ''' =================================================================== [ GetParameters_ ] SGeMS Raw file에서 각 Ensemble 별로 permeability를 추출해내어 pandas의 DataFrame 형태로 데이터를 return함 =================================================================== ''' def GetParameters_(filename, BoolofChannel, k_sand, k_clay): with open(filename, 'r') as f: txt = f.read() lines = txt.split('\n') total_text = [line for line in lines] En_num = int(total_text[1]) En_name_ = total_text[2:En_num+2] En_name = [re.findall(r'\d+',En_name_[i]) for i in range(En_num)] En_name = np.array([ii[1] for ii in En_name]) facies = [] for grid in total_text[En_num+2:-1]: facies.append(list(map(int, grid.strip().split(' ')))) facies = np.array(facies) if BoolofChannel: facies[facies >= 1] = k_sand facies[facies < 1] = k_clay else: facies = np.exp(facies) parameters = pd.DataFrame(data=facies, columns=En_name) for kk in range(int(total_text[1])): if kk == 0: data = parameters[str(kk)].to_numpy() else: data = np.column_stack((data, parameters[str(kk)].to_numpy())) parameters = pd.DataFrame(data = data, columns= np.linspace(0,int(total_text[1])-1,int(total_text[1]), dtype=int).astype(np.str_)) return En_name, parameters ''' =================================================================== [ GetDynamicAssimilation ] Reference Model에서 Observation이 되는 Dynamic data를 읽어오기 위함 Initial Model에서의 사용은 시각화를 위함 본 코드에서 Dynamic parameter는 WOPR임 =================================================================== ''' def GetDynamicAssimilation(Atime, filename): with open('..//EnKF/'+filename + '.RSM','r') as f: file = f.read() lines = file.split('\n') lines = [line.strip() for line in lines] if lines[9][0] == str(0): dummy = 10 else: dummy = 9 if type(Atime) == type(np.float64(1)): timelen = 1 else: timelen = len(Atime) WOPR_lines = lines[dummy:dummy + timelen] Dynamic = [] for k in WOPR_lines: Dynamic.append(k.split()[2::]) WOPR = np.array(Dynamic).transpose().flatten().astype(np.float64) return WOPR ''' =================================================================== [ GetDynamicPrediction ] Assimilation에 사용되는 parameter는 WOPR이므로, 나머지는 사실상 Assimilation이 잘 되었는지(불확실성이 줄었는지) 확인하기 위함 =================================================================== ''' def GetDynamicPrediction(Ptime, filename): with open('..//EnKF/'+filename + '.RSM','r') as f: file = f.read() lines = file.split('\n') lines = [line.strip() for line in lines] if lines[9][0] == str(0): dummy = 10 else: dummy = 9 if type(Ptime) == type(np.float64(1)): timelen = 1 else: timelen = len(Ptime) WOPR_lines = lines[dummy:dummy + timelen] Dynamic = [] for k in WOPR_lines: Dynamic.append(k.split()[2::]) WOPR = np.array(Dynamic).transpose().flatten().astype(np.float64) Other_lines_1 = lines[dummy * 2 + timelen :(dummy + timelen) * 2] Dynamic = [] for k in Other_lines_1: Dynamic.append(k.split()[1::]) WWCT = np.array(Dynamic)[::,0:-1].transpose().flatten().astype(np.float64) FOPT = np.array(Dynamic)[::,-1].transpose().flatten().astype(np.float64) Other_lines_2 = lines[dummy * 3 + timelen * 2:(dummy + timelen) * 3] Dynamic = [] for k in Other_lines_2: Dynamic.append(k.split()[1::]) FWPT = np.array(Dynamic).transpose().flatten().astype(np.float64) return WOPR, WWCT, FOPT, FWPT ## 이번엔 static만 사용했으므로, dynamic은 추후에 확장할 예정 # mean vector 생성 ''' =================================================================== [ GetMean ] EnKF에서 최소자승법을 적용할 때, 참값을 ensemble의 평균으로 가정 따라서, State vector에서 Static parameter는 데이터별로 mean을 적용하고 그 값에 따라 시뮬레이션한 결과로 Observation을 새로 구한다. 이로써 [static mean vector ; Observation] 형태의 vector를 return한다. 💌 이번에는 static만 사용했으므로, dynamic parameter들은 추후에 사용할때 확장하여 사용한다. =================================================================== ''' def GetMean(static,dir,filename): enMean = np.mean(static.to_numpy(),1) enMean = pd.DataFrame(data = enMean, columns= ['mean']) EnKF.MakePermx('mean', enMean, 'PERMX.DAT') EnKF.RunECL(dir, filename) os.chdir('../Ensemble/EnKF') WOPR = EnKF.GetDynamicAssimilation([20],filename) enMean = np.append(enMean.to_numpy().reshape(1,-1), WOPR) return enMean ''' =================================================================== [ Run Series ] 어떠한 forward simulation을 이용할지 선정, format은 동일 RunECL : ECL100을 이용한 시뮬레이션(FDM based Sim.) RunTRACE3D : TRACE3D를 이용한 시뮬레이션(Streamline Sim.) 💌 본 코드는 ECL100기반으로 아웃풋을 받도록 설계되어 있으므로, TRACE3D를 이용할 경우 Dynamic data를 읽어들이는 방식을 수정하여야함 =================================================================== ''' def RunECL(directory, fileName): runCommand = 'C:\\ecl\\2022.1\\bin\\pc_x86_64\\eclipse.exe ' + fileName + '>NUL' os.chdir('../EnKF/' + directory) subprocess.run(runCommand, shell=True, check=True) os.chdir('..') return def RunTRACE3D(directory, fileName): runCommand = 'C:\TRACE3D\T3DMain.exe ' + fileName + '>NUL' os.chdir(directory) subprocess.run(runCommand, shell=True, check=True) os.chdir('..') return ''' =================================================================== [ Ensemble_Kalman_Filter ] Assimilation을 수행하는 핵심적인 코드 Static and Dynamic parameters, Ensemble mean vector, True Observation(=Result of Reference model)를 입력받아 Posteriori state vector 전체와 static vector를 return함. row size = ngrid | column size = Number of Ensemble인 matrix 반환 =================================================================== ''' def Ensemble_Kalman_Filter(Observation, Static_parameter, Dynamic_parameter,En_mean, Nofrealization, Nofobservedwell, StdErr_Dynamic): # Log conversion과 Static vector 구성 En_mean = En_mean[np.newaxis,:].reshape(-1,1) En_mean[:-Nofobservedwell,::] = np.log(En_mean[:-Nofobservedwell,::]) Static_parameter = np.log(Static_parameter) State_vector = np.vstack([Static_parameter, Dynamic_parameter]) # 측정행렬 H = np.zeros([Nofobservedwell ,State_vector.shape[0]]) for i in range(Nofobservedwell): H[-1 - i,-1 - i] = 1 # 추정오차공분산 L = State_vector - ( En_mean * np.ones([1,Nofrealization])) Cov_e = 1 / (Nofrealization-1) * L @ L.transpose() # 노이즈 추가 Observation = Observation.reshape(Nofobservedwell,-1) * np.ones([Nofobservedwell, Nofrealization]) + pow(StdErr_Dynamic,2) * np.random.rand(Nofobservedwell, Nofrealization) C_d = np.zeros([Nofobservedwell,Nofobservedwell]) for kk in range(Nofobservedwell): C_d[kk, kk] = pow(StdErr_Dynamic,2) # 칼만 게인 K = Cov_e @ H.transpose() @ np.linalg.inv(H @ Cov_e @ H.transpose() + C_d) # 교정 메인수식 State_vector_Assimilation = State_vector + K @ (Observation - H @ State_vector) static_Assimilation = State_vector_Assimilation[0:-Nofobservedwell ,::] # Exp conversion static_Assimilation = np.exp(static_Assimilation) State_vector_Assimilation[0:-Nofobservedwell ,::] = static_Assimilation return State_vector_Assimilation, static_Assimilation ''' =================================================================== [ DeleteFolder ] 모든 ini, EnKF 폴더 삭제 =================================================================== ''' def DeleteFolder(NofRealization): [shutil.rmtree('EnKF{}'.format(i)) for i in range(1,NofRealization)] [shutil.rmtree('ini{}'.format(i)) for i in range(1, NofRealization)]