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 Smoother module in Ensembles package
from Ensembles import ES
NofRealization = 101
dx = dy = dz = 50
nx=30
ny=30
nz=1
ngrid = nx * ny * nz
Nofobservedwell = 8
ASTEP = 1
source = os.getcwd().split('\\')[-1]
[(Method, Model)] = re.findall(r'(\w{2})\_(\w+)',source)
Atime=np.linspace(20,600,30)
Ptime=np.linspace(20,1800,90)
NofAtime=len(Atime)
NofPtime = len(Ptime)
Modelvariable = ["Permx"]
Dynamicvariable = ["WOPR"]
NofModelvariable = len(Modelvariable);
NofDataType = len(Dynamicvariable)
filename = 'SGSSim'
refNum = str(20)
ES.MakeTstep(NofAtime,Atime[0],'TSTEP.DAT') # Atime
if os.path.isdir('Log') == False:
os.mkdir('Log')
for i in range(1,NofRealization):
# Assimilation folder
if os.path.isdir(f'{Method}{i}') == False:
os.mkdir(f'{Method}{i}')
shutil.copy('WOOGIE.DATA', f'{Method}{i}/WOOGIE.DATA')
shutil.copy('TSTEP.DAT', f'{Method}{i}/TSTEP.DAT')
shutil.copy('SOLUTION.DAT', f'{Method}{i}/SOLUTION.DAT')
# Initial Ensembels folder
if os.path.isdir(f'ini{i}') == False:
os.mkdir(f'ini{i}')
shutil.copy('WOOGIE.DATA', f'ini{i}/WOOGIE.DATA')
shutil.copy('SOLUTION.DAT', f'ini{i}/SOLUTION.DAT')
ES.MakeTstep(NofPtime,Ptime[0],f'ini{i}/TSTEP.DAT')
# Reference folder
if os.path.isdir('ref') == False:
os.mkdir('ref')
shutil.copy('WOOGIE.DATA','ref/WOOGIE.DATA')
shutil.copy('SOLUTION.DAT','ref/SOLUTION.DAT')
ES.MakeTstep(NofPtime,Ptime[0],f'ref/TSTEP.DAT')
k_sand = 2000
k_clay = 20
RealizationName ,PermxAll = ES.GetParameters(filename, False, k_sand, k_clay)
refPermx = PermxAll[refNum]
encols = [str(i) for i in range(NofRealization)]
encols.remove(refNum)
enPermx = PermxAll[encols]
enPermx.columns = np.linspace(1,100,100, dtype= int).astype(str)
for Realization in enPermx.columns:
# ES.MakePermx(Realization ,enPermx, f'{Method}{Realization}/PERMX.DAT')
ES.MakePermx(Realization ,enPermx, f'ini{Realization}/PERMX.DAT')
ES.MakePermx(refNum ,PermxAll, 'ref/PERMX.DAT')
ES.RunECL(f'{source}/ref','WOOGIE')
refWOPR_A = ES.GetDynamicAssimilation(Atime, 'ref/WOOGIE', source)
refWOPR_P, refWWCT_P, refFOPT_P, refFWPT_P = ES.GetDynamicPrediction(Ptime, 'ref/WOOGIE', source)
OBS = refWOPR_A[np.newaxis,:]
for i in tqdm(range(1,NofRealization), desc = 'Simulate Initial Ensembles'):
time.sleep(0.5)
ES.RunECL(f'{source}/ini{i}', 'WOOGIE')
if i <= 1:
WOPR_A = ES.GetDynamicAssimilation(Atime, 'ini1/WOOGIE', source)
else:
WOPR_A = np.column_stack((WOPR_A, ES.GetDynamicAssimilation(Atime, f'ini{i}/WOOGIE', source)))
cols = [f'Real{jj}' for jj in range(1,NofRealization)]
for i in range(1, ASTEP+1):
'''
===================================================================
Assimilation Step
===================================================================
'''
EnMean = ES.GetMean(enPermx, source,'WOOGIE', Atime)
State_vector_Assimilation, Assimilated_Permx = ES.Ensemble_Smoother(OBS,enPermx.to_numpy(),WOPR_A ,EnMean, NofRealization -1, Nofobservedwell, NofAtime, 0.1)
iniEnMean = EnMean
static_vector = pd.DataFrame(data = Assimilated_Permx, columns = cols)
[ES.MakePermx(f'Real{j}', static_vector, f'{Method}{j}/PERMX.DAT')for j in range(1,NofRealization)]
'''
===================================================================
Prediction Step
===================================================================
'''
for j in tqdm(range(1,NofRealization), desc = f'Prediction No.{i} '):
ES.MakeTstep(NofPtime,Ptime[0],f'{Method}{j}/TSTEP.DAT')
ES.RunECL(f'{source}/{Method}{j}', 'WOOGIE')
ES.GetMean(enPermx, source,'WOOGIE', Ptime) # initial
WOPR_i, WWCT_i, FOPT_i, FWPT_i = ES.GetDynamicPrediction(Ptime, 'WOOGIE', source)
ES.GetMean(static_vector, source,'WOOGIE', Ptime) # Post
WOPR_M, WWCT_M, FOPT_M, FWPT_M = ES.GetDynamicPrediction(Ptime, 'WOOGIE', source)
# nparray를 데이터프레임으로
def cv2pd_W(raw,NofCol,time_):
df = pd.DataFrame(data = raw.reshape([NofCol,-1]).transpose(),
columns=[f'Well{i}' for i in range(1,NofCol+1)],
index = time_
)
return df
def cv2pd_F(raw, time_):
df = pd.DataFrame(data = raw, columns=['Field'], index = time_)
return df
WOPR_M = cv2pd_W(WOPR_M,Nofobservedwell, Ptime)
WWCT_M = cv2pd_W(WWCT_M,Nofobservedwell, Ptime)
FOPT_M = cv2pd_F(FOPT_M, Ptime)
FWPT_M = cv2pd_F(FWPT_M, Ptime)
iniWOPR_M = cv2pd_W(WOPR_i,Nofobservedwell, Ptime)
iniWWCT_M = cv2pd_W(WWCT_i,Nofobservedwell, Ptime)
iniFOPT_M = cv2pd_F(FOPT_i, Ptime)
iniFWPT_M = cv2pd_F(FWPT_i, Ptime)
from ImVisual import ImVisual as Img
ESWOPR, ESWWCT, ESFOPT, ESFWPT = ES.DataBox('ES', Nofobservedwell,NofRealization,Ptime, source)
iniWOPR, iniWWCT, iniFOPT, iniFWPT = ES.DataBox('ini', Nofobservedwell,NofRealization,Ptime,source)
Img.Heat_n_Hist(static_vector, refPermx, enPermx, Method, Model , nx, ny)
Img.Plot_by_Well(iniWOPR,iniWOPR_M,refWOPR_P, Nofobservedwell,NofRealization, Ptime,'ini', Model , 'WOPR')
Img.Plot_by_Well_ES(WOPR_M,refWOPR_P, Nofobservedwell,NofRealization, Ptime,Method, Model , 'WOPR')
Img.Plot_by_Well(ESWOPR,WOPR_M,refWOPR_P, Nofobservedwell,NofRealization, Ptime,f'{Method}_ALL', Model , 'WOPR')
Img.Plot_by_Well(iniWWCT,iniWWCT_M,refWWCT_P, Nofobservedwell,NofRealization, Ptime,'ini', Model , 'WWCT')
Img.Plot_by_Well_ES(WWCT_M,refWWCT_P, Nofobservedwell,NofRealization, Ptime,Method, Model , 'WWCT')
Img.Plot_by_Well(ESWWCT,WWCT_M,refWWCT_P, Nofobservedwell,NofRealization, Ptime,f'{Method}_ALL', Model , 'WWCT')
Img.Plot_Field(iniFOPT,iniFOPT_M,refFOPT_P, Ptime,'ini', Model , 'FOPT')
Img.Plot_Field_ES(FOPT_M,refFOPT_P, Ptime,Method, Model , 'FOPT')
Img.Plot_Field(ESFOPT,FOPT_M,refFOPT_P, Ptime,f'{Method}_ALL', Model , 'FOPT')
Img.Plot_Field(iniFWPT,iniFWPT_M,refFWPT_P, Ptime,'ini', Model , 'FWPT')
Img.Plot_Field_ES(FWPT_M,refFWPT_P, Ptime,Method, Model , 'FWPT')
Img.Plot_Field(ESFWPT,FWPT_M,refFWPT_P, Ptime,f'{Method}_ALL', Model , 'FWPT')