diff --git a/ExampleData/Cart2D.DATA b/ExampleData/Cart2D.DATA new file mode 100644 index 0000000..5dedeb5 --- /dev/null +++ b/ExampleData/Cart2D.DATA @@ -0,0 +1,76 @@ + +RUNSPEC ===================================================================== + +TITLE +validation + +START +1 'JAN' 2000 / + +FIELD + +FMTIN + +UNIFIN + +FMTOUT + +UNIFOUT + + +WATER + + +DIMENS +50 50 1 / + +EQLDIMS + 1 100 10 1 1 / + +REGDIMS + 1 1 0 10 / + +TABDIMS + 1 1 16 12 1 12 20 1 / + +WELLDIMS + 3 2 2 3 / + +TRACERS +0 1 0 0/ + +GRID =============================================================== + +GRIDFILE + 1 1 / + +INIT + + +DX +2500*20 +/ + +DY +2500*20 +/ + +DZ +2500*6 +/ + +TOPS +2500*1000 +/ + +PORO +2500*0.3 / + +PERMX +2500*300 / + +PERMY +2500*300 / + +PERMZ +2500*50 / diff --git a/ExampleData/Cart2D_Fault.GRDECL b/ExampleData/Cart2D_Fault.GRDECL new file mode 100644 index 0000000..6a4fdf1 --- /dev/null +++ b/ExampleData/Cart2D_Fault.GRDECL @@ -0,0 +1,84 @@ +GRID + +DIMENS +15 8 1 / + +TOPS +9310 9352 9342 12*9342 +9340 9318 9318 9338 9322 4*9322 9315 9309 9308 9308 2*9308 +9340 9321 9317 9325 9325 9315 9302 9300 9300 9298 9298 9295 9298 9302 9308 +9342 9320 9320 9310 9315 9317 9310 9280 9298 9290 9290 9294 9292 9300 9300 +9340 9340 9335 9320 9295 9295 9298 9290 9294 9297 9295 9291 3*9291 +5*9290 9290 9290 9286 9295 9288 5*9288 +7*9296 9296 9298 9282 9295 4*9295 +7*9310 9310 9288 9279 5*9279 / + +DX +294.12 382.35 500 794.12 470.58 323.53 352.94 588.24 617.65 441.17 426.47 426.47 441.17 352.94 411.76 +294.12 382.35 500 794.12 470.58 323.53 352.94 588.24 617.65 441.17 426.47 426.47 441.17 352.94 411.76 +294.12 382.35 500 794.12 470.58 323.53 352.94 588.24 617.65 441.17 426.47 426.47 441.17 352.94 411.76 +294.12 382.35 500 794.12 470.58 323.53 352.94 588.24 617.65 441.17 426.47 426.47 441.17 352.94 411.76 +294.12 382.35 500 794.12 470.58 323.53 352.94 588.24 617.65 441.17 426.47 426.47 441.17 352.94 411.76 +294.12 382.35 500 794.12 470.58 323.53 352.94 588.24 617.65 441.17 426.47 426.47 441.17 352.94 411.76 +294.12 382.35 500 794.12 470.58 323.53 352.94 588.24 617.65 441.17 426.47 426.47 441.17 352.94 411.76 +294.12 382.35 500 794.12 470.58 323.53 352.94 588.24 617.65 441.17 426.47 426.47 441.17 352.94 411.76 / + +DY +15*500 +15*205.88 +15*382.35 +15*529.41 +15*441.17 +15*441.17 +15*441.17 +15*352.94 / + +DZ +5 12 11 12*11 +10 10 32 14 1 4*1 4 5.5 2 3*2 +10 42 40 28 18 15 13 12 12 15 13 8 10 5 3 +12 36 33 40 34 34 33 32 27 24 22 17 17 5 5 +7 7 6 24 42 42 40 30 18 10 8 5.5 3*5.5 +5*5 5 7 12 10 3 5*3 +8*6 13 6 24 4*24 +8*0 5 4 5*4 / + +PORO +19.2 19.4 19.7 12*0 +19.2 19.51 20.0 20.5 21.0 4*0 21.2 20.5 20.0 3*0 +19.1 19.5 20.2 20.75 21.2 21.5 21.55 22.2 22.2 21.4 21.0 20.5 20.4 22.25 20.0 +18.8 19.5 20.4 21.25 21.54 22.23 22.50 22.7 23.0 22.0 21.6 21.3 20.5 20.4 19.9 +0 19.0 20.0 21.0 22.0 22.5 23.0 24.0 23.5 22.5 22.1 21.6 3*0 +5*0 22.3 22.7 24.0 23.3 21.8 5*0 +7*0 23.8 23.1 22.0 21.7 4*0 +7*0 23.0 22.5 21.8 5*0 / + +PERMX +272 275 261 12*261 +273 276 273 260 248 4*248 260 270 278 3*278 +270 279 284 270 265 267 268 271 270 272 275 285 285 277 273 +261 275 298 288 280 278 275 285 275 285 300 288 270 265 260 +255 261 275 286 290 282 275 290 280 276 279 270 3*260 +6*270 267 281 296 275 5*265 +8*279 283 275 5*265 +8*257 270 268 5*260 / + +PERMY +217.6 220 208.8 12*245 +218.4 220.8 220 208 5*198.4 208 216 4*222.4 +216 223.2 168.4 220 265 213.6 214.4 216.8 216 217.6 220.8 228 228 221.8 218.4 +208.8 220 238.4 230.4 224 222.4 220 228 220 228 240 230.4 216 212 208 +204 208.8 220 228.8 232 225.6 220 232 224 220.8 223.2 216 3*208 +6*216 213.6 224.8 236.8 6*220 +8*223.2 226.4 220 5*212 +8*205.6 216 214.4 5*208 / + +PERMZ +0.1 0.1 0.1 12*0.1 +15*0.1 +0.1 10*1 4*0.1 +0.1 11*1 3*0.1 +0.1 10*1 4*0.1 +5*0.1 10*1 +7*0.1 8*1 +15*0.1 / diff --git a/Example_Builder.ipynb b/Example_Builder.ipynb new file mode 100644 index 0000000..a19e598 --- /dev/null +++ b/Example_Builder.ipynb @@ -0,0 +1,131 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "from GRDECL2VTK import *" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Grid Type=Cartesian Grid\n", + " Grid Dimension(NX,NY,NZ): (50 x 50 x 15)\n", + " NumOfGrids=37500\n", + "[Geometry] Converting GRDECL to Paraview Hexahedron mesh data....\n", + " NumOfPoints 300000\n", + " NumOfCells 37500\n", + " .....Done!\n", + "[Output] Writing \"Results\\Cart3D.vtu\" Paraview file....Done!\n", + "NPSL file [Results\\Cart3D_permx.txt] successfully genetrated, pelase use NPSL to load it!\n", + "NPSL file [Results\\Cart3D_permy.txt] successfully genetrated, pelase use NPSL to load it!\n", + "NPSL file [Results\\Cart3D_permz.txt] successfully genetrated, pelase use NPSL to load it!\n", + "NPSL file [Results\\Cart3D_poro.txt] successfully genetrated, pelase use NPSL to load it!\n" + ] + } + ], + "source": [ + "#Model=GeologyModel(filename='./ExampleData/Cart2D.DATA')\n", + "#Model=GeologyModel(filename='./ExampleData/Cart2D_Fault.GRDECL')\n", + "Model=GeologyModel()\n", + "Model.fname='Cart3D'\n", + "Model.buildCartGrid(physDims=[50.0,50.0,10.0],gridDims=[50,50,15])\n", + "\n", + "Model.UpdateCellData(varname=\"PERMX\",val=200)\n", + "Model.UpdateCellData(varname=\"PERMY\",val=200)\n", + "Model.UpdateCellData(varname=\"PERMZ\",val=10)\n", + "\n", + "Model.UpdateCellData(varname=\"PERMX\",val=50,nx_range=(10,15),ny_range=(10,15))\n", + "\n", + "\n", + "Model.GRDECL2VTK()\n", + "Model.Write2VTU()\n", + "Model.WriteNPSL()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([1., 1., 1., 1., 6., 6., 6., 6.])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Model.GRDECL_Data.TOPS" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "a,b,c=[1,2,3]\n", + "a,b,c=np.zeros((3,10))" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([1., 1., 1.])" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.ones(3)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/GRDECL2VTK.py b/GRDECL2VTK.py index 8abed8b..af618a0 100644 --- a/GRDECL2VTK.py +++ b/GRDECL2VTK.py @@ -27,7 +27,7 @@ from GRDECL_Parser import * from GRDECL_FaultProcess import * #from GRDECL_CADExporter import * - +from utils import * ############################################# # @@ -65,10 +65,57 @@ def readGRDECL(self,filename): self.GRDECL_Data.fname=filename self.GRDECL_Data.read_GRDECL() + def buildCartGrid(self,physDims=[100.0,100.0,10.0],gridDims=[10,10,1]): + #* Create simple cartesian grid + self.GRDECL_Data.buildCartGrid(physDims,gridDims) + def GRDECL2VTK(self): #* Convert corner point grid/cartesian grid into VTK unstructure grid print('[Geometry] Converting GRDECL to Paraview Hexahedron mesh data....') NX,NY,NZ=self.GRDECL_Data.NX,self.GRDECL_Data.NY,self.GRDECL_Data.NZ + + if(self.GRDECL_Data.GRID_type=='Cartesian'): + #1. Collect points from Cartesian data + debug=0 + DX,DY,DZ,TOPS=self.GRDECL_Data.DX,self.GRDECL_Data.DY,self.GRDECL_Data.DZ,self.GRDECL_Data.TOPS + coordX,coordY,coordZ=Cartesian2UnstructGrid(DX,DY,DZ,TOPS,NX,NY,NZ) + + Points = vtk.vtkPoints() + Points.SetNumberOfPoints(2*NX*2*NY*2*NZ) #=2*NX*2*NY*2*NZ + + ptsid=0 + for k in range(2*NZ): + for j in range(2*NY): + for i in range(2*NX): + x,y,z=coordX[i][j][k],coordY[i][j][k],coordZ[i][j][k] + Points.SetPoint(ptsid,[x,y,z]) + ptsid+=1 + self.VTK_Grids.SetPoints(Points) + + #2. Recover Cells which follows the cartesian data + cellArray = vtk.vtkCellArray() + Cell=vtk.vtkHexahedron() + + cellid=0 + for k in range(NZ): + for j in range(NY): + for i in range(NX): + idx_GB=[getIJK(2*i,2*j,2*k,2*NX,2*NY,2*NZ), #index-0 + getIJK(2*i+1,2*j,2*k,2*NX,2*NY,2*NZ), #index-1 + getIJK(2*i+1,2*j+1,2*k,2*NX,2*NY,2*NZ), #index-2 + getIJK(2*i,2*j+1,2*k,2*NX,2*NY,2*NZ), #index-3 + getIJK(2*i,2*j,2*k+1,2*NX,2*NY,2*NZ), #index-4 + getIJK(2*i+1,2*j,2*k+1,2*NX,2*NY,2*NZ), #index-5 + getIJK(2*i+1,2*j+1,2*k+1,2*NX,2*NY,2*NZ),#index-6 + getIJK(2*i,2*j+1,2*k+1,2*NX,2*NY,2*NZ)] #index-7 + for pi in range(8): + Cell.GetPointIds().SetId(pi,idx_GB[pi]) + cellArray.InsertNextCell(Cell) + cellid+=1 + + self.VTK_Grids.SetCells(Cell.GetCellType(), cellArray) + + if(self.GRDECL_Data.GRID_type=='CornerPoint'): #1.Collect Points from the raw CornerPoint data [ZCORN]&[COORD] @@ -107,14 +154,14 @@ def GRDECL2VTK(self): self.VTK_Grids.SetCells(Cell.GetCellType(), cellArray) - print(" NumOfPoints",self.VTK_Grids.GetNumberOfPoints()) - print(" NumOfCells",self.VTK_Grids.GetNumberOfCells()) - - #3. Load grid properties data if applicable - for keyword,data in self.GRDECL_Data.SpatialDatas.items(): - self.AppendScalarData(keyword,data) + + print(" NumOfPoints",self.VTK_Grids.GetNumberOfPoints()) + print(" NumOfCells",self.VTK_Grids.GetNumberOfCells()) + #3. Load grid properties data if applicable + for keyword,data in self.GRDECL_Data.SpatialDatas.items(): + self.AppendScalarData(keyword,data) - print('.....Done!') + print(' .....Done!') def decomposeModel(self): '''#* Identify and extract boundary/falut faces @@ -166,13 +213,90 @@ def decomposeModel(self): DomainMarker3D=np.tile(DomainMarker2D,self.GRDECL_Data.NZ) self.AppendScalarData('SubVolumes',DomainMarker3D) - + def AppendScalarData(self,name,numpy_array): #* Append scalar cell data (numpy array) into vtk object data = ns.numpy_to_vtk(numpy_array.ravel(order='F'),deep=True, array_type=vtk.VTK_FLOAT) data.SetName(str(name)) data.SetNumberOfComponents(1) self.VTK_Grids.GetCellData().AddArray(data) + + def UpdateCellData(self,varname,val=100,var="PERMX",nx_range=(1,-1),ny_range=(1,-1),nz_range=(1,-1)): + """Update/modify Permeability/Porosity field with given grid block range + + Arguments + --------- + var -- The varable name you want to update, e.g PERMX, PERMY or PERMZ + val -- The varable value you want to update + nx_range -- The specifc grid range in x for updating, 1-based index + nx_range -- The specifc grid range in y for updating, 1-based index + nz_range -- The specifc grid range in z for updating, 1-based index + + Author:Bin Wang(binwang.0213@gmail.com) + Date: Feb. 2018 + """ + + assert varname in self.GRDECL_Data.SpatialDatas, "[Error] Variable [%s] is not existed!" + + nx_range=np.array(nx_range) + ny_range=np.array(ny_range) + nz_range=np.array(nz_range) + + #If no nx,ny,nz range are defined, all perm will be updated + if(nx_range[1] == -1): + nx_range[1] = self.GRDECL_Data.NX + if(ny_range[1] == -1): + ny_range[1] = self.GRDECL_Data.NY + if(nz_range[1] == -1): + nz_range[1] = self.GRDECL_Data.NZ + + #Convert 1-based index to 0-based index + nx_range=nx_range-1 + ny_range=ny_range-1 + nz_range=nz_range-1 + + #Update perm field with specific grid range + for k in range(self.GRDECL_Data.NZ): + for j in range(self.GRDECL_Data.NY): + for i in range(self.GRDECL_Data.NX): + if(i>=nx_range[0] and i<=nx_range[1]): + if(j>=ny_range[0] and j<=ny_range[1]): + if(k>=nz_range[0] and k<=nz_range[1]): + ijk = getIJK(i, j, k, self.GRDECL_Data.NX, self.GRDECL_Data.NY, self.GRDECL_Data.NZ) + self.GRDECL_Data.SpatialDatas[varname][ijk]=val + + + def WriteNPSL(self): + """Write the permeability/porosity field for NPSL + + filename_permx.txt + filename_permy.txt + filename_permz.txt + filename_poro.txt + + Arguments + --------- + filename -- The surname of permeability and porosity data files + + Programmer: Bin Wang (yin.feng@louisiana.edu) + Creation: Feb, 2018 + """ + basename=os.path.splitext(os.path.basename(self.fname))[0] + if not os.path.exists("Results"): + os.makedirs('Results') + fnames=[os.path.join('Results',basename + '_permx.txt'), + os.path.join('Results',basename + '_permy.txt'), + os.path.join('Results',basename + '_permz.txt'), + os.path.join('Results',basename + '_poro.txt')] + + np.savetxt(fnames[0], self.GRDECL_Data.SpatialDatas['PERMX'], delimiter="\n",fmt='%1.4f') + np.savetxt(fnames[1], self.GRDECL_Data.SpatialDatas['PERMY'], delimiter="\n",fmt='%1.4f') + np.savetxt(fnames[2], self.GRDECL_Data.SpatialDatas['PERMZ'], delimiter="\n",fmt='%1.4f') + np.savetxt(fnames[3], self.GRDECL_Data.SpatialDatas['PORO'], delimiter="\n",fmt='%1.4f') + + for name in fnames: + print('NPSL file [%s] successfully genetrated, pelase use NPSL to load it!' % (name)) + def Write2VTU(self): basename=os.path.splitext(os.path.basename(self.fname))[0] diff --git a/GRDECL_Parser.py b/GRDECL_Parser.py index dd50e91..c513ed0 100644 --- a/GRDECL_Parser.py +++ b/GRDECL_Parser.py @@ -35,6 +35,7 @@ float,float,float ] + class GRDECL_Parser: def __init__(self,filename='',nx=0,ny=0,nz=0): """Eclipse Input file(GRDECL) Parser @@ -69,6 +70,7 @@ def __init__(self,filename='',nx=0,ny=0,nz=0): #Petrophysics data Keywords self.SpatialDatas={} + self.SkipedKeywords=0 #Read GRDECL file when initializing the class if(len(filename)>0): @@ -104,12 +106,15 @@ def read_GRDECL(self): GoodFlag=0 for i,block in enumerate(contents_in_block):#Keyword, Block-wise + #Clean the data where no spliter \ provided + block=scanKeyword(block) + blockData_raw=block.strip().split() Keyword='' DataArray=[] if(len(blockData_raw)>1): if(blockData_raw[0]=='ECHO'): #This keyword may next to real keyword - Keyword,DataArray=blockData_raw[1],blockData_raw[2:] + Keyword,DataArray=blockData_raw[1],blockData_raw[2:] else: Keyword,DataArray=blockData_raw[0],blockData_raw[1:] @@ -119,6 +124,7 @@ def read_GRDECL(self): self.GRID_type='Cartesian' self.NX,self.NY,self.NZ=DataArray[0],DataArray[1],DataArray[2] self.N=self.NX*self.NY*self.NZ + print(" Grid Type=%s Grid" %(self.GRID_type)) print(" Grid Dimension(NX,NY,NZ): (%s x %s x %s)"%(self.NX,self.NY,self.NZ)) print(" NumOfGrids=%s"%(self.N)) print(' NumOfKeywords=%s'%(NumKeywords)) @@ -130,6 +136,7 @@ def read_GRDECL(self): self.GRID_type='CornerPoint' self.NX,self.NY,self.NZ=DataArray[0],DataArray[1],DataArray[2] self.N=self.NX*self.NY*self.NZ + print(" Grid Type=%s" %(self.GRID_type)) print(" Grid Dimension(NX,NY,NZ): (%s x %s x %s)"%(self.NX,self.NY,self.NZ)) print(" NumOfGrids=%s"%(self.N)) print(' NumOfKeywords=%s'%(NumKeywords)) @@ -139,14 +146,35 @@ def read_GRDECL(self): if(self.GRID_type=='NaN'):#Skip unnecessary keywords continue + + if(Keyword in SupportKeyWords): #We need parse the special format in + DataArray=self.parseDataArray(DataArray) + #print(Keyword,DataArray) + #Read Grid spatial information, x,y,z ordering + #Corner point cell if(Keyword=='COORD'):# Pillar coords assert len(DataArray)==6*(self.NX+1)*(self.NY+1),'[Error] Incompatible COORD data size!' self.COORD=np.array(DataArray,dtype=float) elif(Keyword=='ZCORN'):# Depth coords assert len(DataArray)==8*self.N, '[Error] Incompatible ZCORN data size!' self.ZCORN=np.array(DataArray,dtype=float) + + #Cartesian cell + elif(Keyword=='DX'):# Grid size in X dir + assert len(DataArray)==self.N, '[Error] Incompatible DX data size!' + self.DX=np.array(DataArray,dtype=float) + elif(Keyword=='DY'):# Grid size in Y dir + assert len(DataArray)==self.N, '[Error] Incompatible DY data size!' + self.DY=np.array(DataArray,dtype=float) + elif(Keyword=='DZ'):# Grid size in Z dir + assert len(DataArray)==self.N, '[Error] Incompatible DZ data size!' + self.DZ=np.array(DataArray,dtype=float) + elif(Keyword=='TOPS'):# TOP position + assert len(DataArray)==self.N, '[Error] Incompatible TOPS data size!' + self.TOPS=np.array(DataArray,dtype=float) + #Read Grid Properties information else: self.LoadVar(Keyword,DataArray,DataSize=self.N) @@ -154,7 +182,58 @@ def read_GRDECL(self): f.close() assert GoodFlag==1,'Can not find grid dimension info, [SPECGRID] or [DIMENS]!' print('.....Done!') + + + #Genetrate TOPS for cartesian grid if TOPS if not given + if(self.GRID_type=='Cartesian' and len(self.TOPS)==0): + self.TOPS=np.zeros(self.N) + for k in range(self.NZ-1): + for j in range(self.NY): + for i in range(self.NX): + ijk=getIJK(i,j,k,self.NX,self.NY,self.NZ) + ijk_next=getIJK(i,j,k+1,self.NX,self.NY,self.NZ) + self.TOPS[ijk_next] = self.TOPS[ijk] + self.DZ[ijk] + + def buildCartGrid(self,physDim=[100.0,100.0,10.0],gridDims=[10,10,1]): + """Build simple cartesian grid + + Arguments + --------- + physDim -- physical dimensions of system + gridDims -- grid dimension of system + + Author:Bin Wang(binwang.0213@gmail.com) + Date: Feb. 2019 + """ + self.NX,self.NY,self.NZ=gridDims + self.N=self.NX*self.NY*self.NZ + self.GRID_type='Cartesian' + + + print(" Grid Type=%s Grid" %(self.GRID_type)) + print(" Grid Dimension(NX,NY,NZ): (%s x %s x %s)"%(self.NX,self.NY,self.NZ)) + print(" NumOfGrids=%s"%(self.N)) + + #Assign value to cart grid + self.DX=np.ones(self.N)*physDim[0]/self.NX + self.DY=np.ones(self.N)*physDim[1]/self.NY + self.DZ=np.ones(self.N)*physDim[2]/self.NZ + self.TOPS=np.zeros(self.N) + for k in range(self.NZ-1): + for j in range(self.NY): + for i in range(self.NX): + ijk=getIJK(i,j,k,self.NX,self.NY,self.NZ) + ijk_next=getIJK(i,j,k+1,self.NX,self.NY,self.NZ) + self.TOPS[ijk_next] = self.TOPS[ijk] + self.DZ[ijk] + + #Build up basic spatial propertis + self.SpatialDatas["PERMX"]=np.ones(self.N)*10.0 + self.SpatialDatas["PERMY"]=np.ones(self.N)*10.0 + self.SpatialDatas["PERMZ"]=np.ones(self.N)*10.0 + self.SpatialDatas["PORO"]=np.ones(self.N)*0.3 + + def LoadVar(self,Keyword,DataArray,DataSize): """Load varables into class example: @@ -163,12 +242,45 @@ def LoadVar(self,Keyword,DataArray,DataSize): Date: Sep. 2018 """ if(Keyword in SupportKeyWords):#KeyWords Check - assert len(DataArray)==DataSize,'\n [Error] Incompatible data size!' + assert len(DataArray)==DataSize,'\n [Error-%s] Incompatible data size! %d-%d' %(Keyword,len(DataArray),DataSize) KeywordID=SupportKeyWords.index(Keyword) print(' [%s] '%(Keyword),end='') self.SpatialDatas[Keyword]=np.array(DataArray,dtype=KeyWordsDatatypes[KeywordID]) else: - print('\n [Warnning] Unsupport keywords[%s]' % (Keyword)) + if(self.SkipedKeywords==0):print() + print(' [Warnning] Unsupport keywords[%s]' % (Keyword)) + self.SkipedKeywords+=1 + + def parseDataArray(self,DataArray): + """Parse special dataArray format in GRDECL + example: + 5*3.0=[3.0 3.0 3.0 3.0 3.0] + 1.0 2*3.0 5.0=[1.0 3.0 3.0 5.0] + + Author:Bin Wang(binwang.0213@gmail.com) + Date: Sep. 2018 + """ + + data=[] + error_count=0 + for value in DataArray: + if(is_number(value)==2): + num,val=value.split('*') + for i in range(int(num)): data.append(val) + elif(is_number(value)==1): + data.append(value) + else: + error_count+=1 + + if(error_count>0): + print(DataArray) + + assert error_count==0, '[Error] Can not find any numeric value!' + + + + return data + def read_IncludeFile(self,filename_include,NumData): """Read Include data file @@ -528,6 +640,14 @@ def RemoveCommentLines(data,commenter='--'): newdata.append(line) return '\n'.join(newdata) +def scanKeyword(data): + #scan and find the keyword + #e.g. ['INIT','DX','2500*30.0'] -> ['DX','2500*30.0'] + for key in SupportKeyWords: + if (key in data) and (data.find(key)!=0): + return data[data.find(key):-1] + return data + def is_number(s): #Determine a string is a number or not #Used in [read_GRDECL] [getBlkdata] @@ -544,7 +664,7 @@ def is_number(s): except (TypeError, ValueError): pass - try: + try: #Special format N*val= [val val val ....] num, val = s.split('*') return 2 except ValueError: diff --git a/utils.py b/utils.py new file mode 100644 index 0000000..affa190 --- /dev/null +++ b/utils.py @@ -0,0 +1,74 @@ +import numpy as np + + +def getIJK(i,j,k,NX,NY,NZ): + #Convert index [i,j,k] to a flat 3D matrix index [ijk] + return i + (NX)*(j + k*(NY)) + +def Cartesian2UnstructGrid(DX,DY,DZ,TOPS,NX,NY,NZ): + #Convert self.DX to coordiantes-X of unstructure gridblocks + #Used in [write_VTU] + debug=0 + + coordX,coordY,coordZ=np.zeros((2*NX,2*NY,2*NZ)),np.zeros((2*NX,2*NY,2*NZ)),np.zeros((2*NX,2*NY,2*NZ)) + for k in range(2*NZ): + for j in range(2*NY): + for i in range(2*NX): + #print(DX[i][j]) + i_GB,j_GB,k_GB=(int)(i/2),(int)(j/2),(int)(k/2) + + #CoordX + if(i==0): + coordX[0][j][k]=0 + if(debug):print('First Node',(0),'-') + if(i>0 and i%2==0): + ijk_GB=getIJK(i_GB-1,j_GB,k_GB,NX,NY,NZ) + coordX[i-1][j][k]=DX[ijk_GB] #odd + coordX[i][j][k]=DX[ijk_GB] #even + if(debug):print('pair',(i-1,i),(i_GB-1,j_GB)) + if(i_GB>1): + coordX[i-1][j][k]+=coordX[i-1-1][j][k] + coordX[i][j][k]=coordX[i-1][j][k] + if(i==2*NX-1): + ijk_GB=getIJK(NX-1,j_GB,k_GB,NX,NY,NZ) + coordX[2*NX-1][j][k]=DX[ijk_GB]+coordX[2*NX-1-1][j][k] + if(debug):print('Last Node',(2*NX-1),(NX-1,j_GB)) + + #CoordY + if(j==0): + coordY[i][0][k]=0 + if(debug and i==0 and k==0):print('First Node',(0),'-') + if(j>0 and j%2==0): + ijk_GB=getIJK(i_GB,j_GB-1,k_GB,NX,NY,NZ) + coordY[i][j-1][k]=DY[ijk_GB] #odd + coordY[i][j][k]=DY[ijk_GB] #even + if(debug and i==0 and k==0):print('pair',(j-1,j),(i_GB,j_GB-1)) + if(j_GB>1): + coordY[i][j-1][k]+=coordY[i][j-1-1][k] + coordY[i][j][k]=coordY[i][j-1][k] + if(j==2*NY-1): + ijk_GB=getIJK(i_GB,NY-1,k_GB,NX,NY,NZ) + coordY[i][2*NY-1][k]=DY[ijk_GB]+coordY[i][2*NY-1-1][k] + if(debug and i==0 and k==0):print('Last Node',(2*NY-1),(i_GB,NY-1)) + + #CoordZ + ijk_GB=getIJK(i_GB,j_GB,k_GB,NX,NY,NZ) + + if(k==0): + coordZ[i][j][0]=TOPS[ijk_GB] + if(debug and i==0 and j==0):print('First Node',(0),'-') + if(k>0 and k%2==0): + ijk_GB=getIJK(i_GB,j_GB,k_GB-1,NX,NY,NZ) + coordZ[i][j][k-1]=DZ[ijk_GB] #odd + coordZ[i][j][k]=DZ[ijk_GB] #even + if(debug and i==0 and j==0):print('pair',(k-1,k),(i_GB,k_GB-1)) + if(k_GB>1): + coordZ[i][j][k-1]+=coordZ[i][j][k-1-1] + coordZ[i][j][k]=coordZ[i][j][k-1] + if(k==2*NZ-1): + ijk_GB=getIJK(i_GB,j_GB,NZ-1,NX,NY,NZ) + coordZ[i][j][2*NZ-1]=DZ[ijk_GB]+coordZ[i][j][2*NZ-1-1] + if(debug and i==0 and j==0):print('Last Node',(2*NZ-1),(i_GB,NZ-1)) + + return coordX,coordY,coordZ +