OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

add a few handy python scripts for pre/post-processing purposes

Browse Source
This commit is contained in:
Rex Zhe Li 2017-09-24 14:13:22 +10:00
parent be59b70139
commit d50b965d81
18 changed files with 2231 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

213
example/Python/LBPM_Postprocess_ID_SignDist_Recon.py Normal file
View File

@ -0,0 +1,213 @@
#!/usr/bin/env python
import sys
import numpy as np
import matplotlib.pyplot as plt
from glob import glob
from netCDF4 import Dataset
def write_NetCDF_file_py(file_name,data_name,data_IN,data_IN_dtype,lx,ly,lz):
# Important !
# Note the dataOUT should have the shape (lz,ly,lx)
data_IN.shape=(lz,ly,lx) # make sure data_IN has the right shape
# open a new netCDF file for writing.
ncfile = Dataset(file_name,'w')
# create the output data.
# create the x, y and z dimensions.
ncfile.createDimension('x',lx)
ncfile.createDimension('y',ly)
ncfile.createDimension('z',lz)
# create the variable (4 byte integer in this case)
# first argument is name of variable, second is datatype, third is
# a tuple with the names of dimensions.
data = ncfile.createVariable(data_name,data_IN_dtype,('z','y','x'))
data[:] = data_IN
#data.voxel_unit = 'micrometer'
#data.voxel_size=5.7
ncfile.close()
print '**Info: The *.nc file is written successfully !'
#end def
# Check if there is a proper command line argument
if len(sys.argv) !=2:
sys.stderr.write('**Error: Usage: ' + sys.argv[0] + ' <Domain.in>\n')
sys.exit()
# end if
# Read 'Domain.in' to obtain the size of 'ID.xxxxxx' and 'SignDist.xxxxx'
f = open(sys.argv[1],'r')
lines = f.readlines()
# Read processors configuration
nproc_x,nproc_y,nproc_z=np.fromstring(lines[0].splitlines()[0],dtype = np.int32,sep=' ')
# Read subdomain dimensions
lx,ly,lz = np.fromstring(lines[1].splitlines()[0],dtype = np.int32,sep=' ')
# Read domain dimensions
LX,LY,LZ = np.fromstring(lines[3].splitlines()[0],dtype = np.float32,sep=' ')
LX=int(LX)
LY=int(LY)
LZ=int(LZ)
f.close()
## -------- Setup output data information -------------##
# 1. For *.nc file
# 1.1 for ID.* files
id_data_output_file='ID_All.nc'
id_data_output_variable='segmented'
# 1.2 for SignDist.* files
dist_data_output_file='SignDist_All.nc'
dist_data_output_variable='signed_distance'
# 1.3 for Phase field files
phase_data_output_file='Phase_All.nc'
phase_data_output_variable='phase_field'
# 2. For *.bin files
id_data_output_file_bin='ID_All.bin'
dist_data_output_file_bin='SignDist_All.bin'
phase_data_output_file_bin='Phase_All.bin'
# The size of the subdomain
lx+=2 # consider the halo layer
ly+=2
lz+=2
id_group = glob('ID.*')
id_group.sort()
dist_group = glob('SignDist.*')
dist_group.sort()
# Reconstruct ID.xxxxx data
for z_axis in np.arange(nproc_z):
for y_axis in np.arange(nproc_y):
for x_axis in np.arange(0,nproc_x,2):
if nproc_x%2==0: # nproc_x is an even number
temp_x1 = np.fromfile(id_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.int8)
temp_x2 = np.fromfile(id_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis+1],dtype=np.int8)
temp_x1.shape = (lz,ly,lx)
temp_x2.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]
temp_x2 = temp_x2[1:lz-1,1:ly-1,1:lx-1]
temp_x1 = np.concatenate((temp_x1,temp_x2),axis=2)
else: # nproc_x is an odd number
if nproc_x==1:
temp_x1 = np.fromfile(id_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.int8)
temp_x1.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]
elif x_axis < nproc_x-2: # i.e. nproc_x = 3 or 5 or 7 ...
temp_x1 = np.fromfile(id_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.int8)
temp_x2 = np.fromfile(id_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis+1],dtype=np.int8)
temp_x1.shape = (lz,ly,lx)
temp_x2.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]
temp_x2 = temp_x2[1:lz-1,1:ly-1,1:lx-1]
temp_x1 = np.concatenate((temp_x1,temp_x2),axis=2)
#end if
#end if
if x_axis == 0:
ID_x_temp = temp_x1.copy()
else:
ID_x_temp = np.concatenate((ID_x_temp,temp_x1),axis=2)
#end if
#end for
if nproc_x !=1 and nproc_x%2 != 0:
temp_x1 = np.fromfile(id_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+nproc_x-1],dtype=np.int8)
ID_x_temp = np.concatenate((ID_x_temp,temp_x1),axis=2)
#end if
if y_axis == 0:
ID_y_temp = ID_x_temp.copy()
else:
ID_y_temp = np.concatenate((ID_y_temp,ID_x_temp),axis=1)
#end if
#end for
if z_axis == 0:
ID_z_temp = ID_y_temp.copy()
else:
ID_z_temp = np.concatenate((ID_z_temp,ID_y_temp),axis=0)
#end if
#end for
## ---------------- ID decoding ------------------------- ##
# By convention, in LBPM-WIA, we have segmented image as:
# 0 -> solid phase
# 1 -> non-wetting phase
# 2 -> wetting phase
## ---------------- ID decoding ------------------------- ##
# By default, the inlet layer is non-wetting phase, and the outlet layer is wetting phase
# Also by default, the flow direction is along z-axis, and the inlet/outlet plane is in x-y plane
ID_z_temp[0,:] =1
ID_z_temp[-1,:]=2
ID_z_temp.tofile(id_data_output_file_bin)
write_NetCDF_file_py(id_data_output_file,id_data_output_variable,ID_z_temp,ID_z_temp.dtype.char,LX,LY,LZ)
# Create phase field based on the ID profiles
# By convention, in LBPM-WIA, we have segmented image as:
# 0 -> solid phase
# 1 -> non-wetting phase
# 2 -> wetting phase
phase = ID_z_temp.copy()
phase = phase.astype(np.float32)
#TODO: make this file read in Color.in file so phi_s can be assigned according to Color.in
phase[phase==0.0] = -0.55 # solid phase
phase[phase==1.0] = 1.0 # non-wetting phase
phase[phase==2.0] = -1.0 # wetting pahse
phase.tofile(phase_data_output_file_bin)
write_NetCDF_file_py(phase_data_output_file,phase_data_output_variable,phase,phase.dtype.char,LX,LY,LZ)
# Reconstruct SignDist.xxxxx data
for z_axis in np.arange(nproc_z):
for y_axis in np.arange(nproc_y):
for x_axis in np.arange(0,nproc_x,2):
if nproc_x%2==0: # nproc_x is an even number
temp_x1 = np.fromfile(dist_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.float64)
temp_x2 = np.fromfile(dist_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis+1],dtype=np.float64)
temp_x1.shape = (lz,ly,lx)
temp_x2.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]
temp_x2 = temp_x2[1:lz-1,1:ly-1,1:lx-1]
temp_x1 = np.concatenate((temp_x1,temp_x2),axis=2)
else: # nproc_x is an odd number
if nproc_x==1:
temp_x1 = np.fromfile(id_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.int8)
temp_x1.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]
elif x_axis < nproc_x-2: # i.e. nproc_x = 3 or 5 or 7 ...
temp_x1 = np.fromfile(dist_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.float64)
temp_x2 = np.fromfile(dist_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis+1],dtype=np.float64)
temp_x1.shape = (lz,ly,lx)
temp_x2.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]
temp_x2 = temp_x2[1:lz-1,1:ly-1,1:lx-1]
temp_x1 = np.concatenate((temp_x1,temp_x2),axis=2)
#end if
#end if
if x_axis == 0:
ID_x_temp = temp_x1.copy()
else:
ID_x_temp = np.concatenate((ID_x_temp,temp_x1),axis=2)
#end if
#end for
if nproc_x !=1 and nproc_x%2 !=0:
temp_x1 = np.fromfile(dist_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+nproc_x-1],dtype=np.float64)
ID_x_temp = np.concatenate((ID_x_temp,temp_x1),axis=2)
#end if
if y_axis == 0:
ID_y_temp = ID_x_temp.copy()
else:
ID_y_temp = np.concatenate((ID_y_temp,ID_x_temp),axis=1)
#end if
#end for
if z_axis == 0:
ID_z_temp = ID_y_temp.copy()
else:
ID_z_temp = np.concatenate((ID_z_temp,ID_y_temp),axis=0)
#end if
#end for
ID_z_temp.tofile(dist_data_output_file_bin)
write_NetCDF_file_py(dist_data_output_file,dist_data_output_variable,ID_z_temp,ID_z_temp.dtype.char,LX,LY,LZ)

196
example/Python/LBPM_Preprocessing_INPUT_writer_utils.py Normal file
View File

@ -0,0 +1,196 @@
#!/usr/bin/env python
import os
import numpy as np
def write_Domain_in_file(**kwargs):
# For most of the parameters in the Color.in file
# you wouldn't need to worry about them as by default
# they are all hard-coded here
Para={'nprocx':1}
Para['nprocy']=1
Para['nprocz']=1
Para['nx']=1
Para['ny']=1
Para['nz']=1
Para['nspheres']=0 # deprecated
Para['Lx']=100
Para['Ly']=100
Para['Lz']=100
for key in kwargs:
if key in Para.keys():
Para[key]=kwargs[key]
else:
print "Error: "+key+" is not a valid input !\n"
print "Error: Domain.in file is not writen."
return
#end if
#end for
# Check if the dompositoin requirement is satisfied
if (Para['nx']*Para['nprocx']>Para['Lx']) or \
(Para['ny']*Para['nprocy']>Para['Ly']) or \
(Para['nz']*Para['nprocz']>Para['Lz']):
print "Error: The decomposed size does not match the total domain size!"
return
#end if
# write the Domain.in file
ParamFile = open("Domain.in","w")
ParamFile.write("%i " % Para['nprocx'])
ParamFile.write("%i " % Para['nprocy'])
ParamFile.write("%i\n" % Para['nprocz'])
ParamFile.write("%i " % Para['nx'])
ParamFile.write("%i " % Para['ny'])
ParamFile.write("%i\n" % Para['nz'])
ParamFile.write("%i\n" % Para['nspheres'])
ParamFile.write("%.1f " % Para['Lx'])
ParamFile.write("%.1f " % Para['Ly'])
ParamFile.write("%.1f\n" % Para['Lz'])
ParamFile.close()
print "**Info: A Domain.in file is created."
#end def
def write_Segment_in_file(**kwargs):
# For most of the parameters in the Color.in file
# you wouldn't need to worry about them as by default
# they are all hard-coded here
Para={'Seg_data_name':'Micromodel_1_segmented.raw'}
Para['Nx']=100
Para['Ny']=100
Para['Nz']=100
Para['xStart']=0
Para['yStart']=0
Para['zStart']=0
for key in kwargs:
if key in Para.keys():
Para[key]=kwargs[key]
else:
print "Error: "+key+" is not a valid input !\n"
print "Error: Segmented.in file is not writen."
return
#end if
#end for
ParamFile = open("Segmented.in","w")
ParamFile.write("%s\n" % Para['Seg_data_name'])
ParamFile.write("%i " % Para['Nx'])
ParamFile.write("%i " % Para['Ny'])
ParamFile.write("%i\n" % Para['Nz'])
ParamFile.write("%i " % Para['xStart'])
ParamFile.write("%i " % Para['yStart'])
ParamFile.write("%i\n" % Para['zStart'])
ParamFile.close()
print "**Info: A Segmented.in file is created."
#end def
def write_Color_in_file(**kwargs):
# For most of the parameters in the Color.in file
# you wouldn't need to worry about them as by default
# they are all hard-coded here
Para={'tau':0.7}
Para['alpha']=0.001
Para['beta']=0.95
Para['phi_solid']=-0.98
Para['saturation']=0.7 # wetting phase saturation, depricated
Para['Fx']=0.0
Para['Fy']=0.0
Para['Fz']=0.0
Para['InitialCondition']=0
Para['BoundaryCondition']=0
Para['din']=0.0
Para['dout']=0.0
Para['TimeStepMax']=100005
Para['Restart_interval']=2000
Para['tolerance']=1e-5
Para['Blob_analysis_interval'] = 1000 #NOTE: not used
for key in kwargs:
if key in Para.keys():
Para[key]=kwargs[key]
else:
print "Error: "+key+" is not a valid input !\n"
print "Error: Color.in file is not writen."
return
#end if
#end for
# write the color.in file
ParamFile = open("Color.in","w")
ParamFile.write("%.3f\n" % Para['tau'])
ParamFile.write("%.3e " % Para['alpha'])
ParamFile.write("%.3f " % Para['beta'])
ParamFile.write("%.3f\n" % Para['phi_solid'])
ParamFile.write("%.2f\n" % Para['saturation'])
ParamFile.write("%.3e " % Para['Fx'])
ParamFile.write("%.3e " % Para['Fy'])
ParamFile.write("%.3e\n" % Para['Fz'])
ParamFile.write("%i " % Para['Restart'])
ParamFile.write("%i " % Para['pBC'])
ParamFile.write("%.3f " % Para['din'])
ParamFile.write("%.3f\n" % Para['dout'])
ParamFile.write("%i " % Para['maxtime'])
ParamFile.write("%i " % Para['interval'])
ParamFile.write("%.2e\n" % Para['tolerance'])
ParamFile.close()
print "**Info: A Color.in file is created."
#end def
def write_Color_Macro_in_file(**kwargs):
# For most of the parameters in the Color.in file
# you wouldn't need to worry about them as by default
# they are all hard-coded here
Para={'tau1':1.0}
Para['tau2'] = 1.0
Para['alpha']=0.001
Para['beta']=0.95
Para['phi_solid']=-0.98
Para['saturation']=0.7 # wetting phase saturation, depricated
Para['Fx']=0.0
Para['Fy']=0.0
Para['Fz']=0.0
Para['InitialCondition']=0
Para['BoundaryCondition']=0
Para['din']=0.0
Para['dout']=0.0
Para['TimeStepMax']=100005
Para['Restart_interval']=2000
Para['tolerance']=1e-5
Para['Blob_analysis_interval'] = 1000 #NOTE: not used
for key in kwargs:
if key in Para.keys():
Para[key]=kwargs[key]
else:
print "Error: "+key+" is not a valid input !\n"
print "Error: Color.in file is not writen."
return
#end if
#end for
# write the color.in file
ParamFile = open("Color.in","w")
ParamFile.write("%.3f " % Para['tau1'])
ParamFile.write("%.3f " % Para['tau2'])
ParamFile.write("%.3e " % Para['alpha'])
ParamFile.write("%.3f " % Para['beta'])
ParamFile.write("%.3f\n" % Para['phi_solid'])
ParamFile.write("%.2f\n" % Para['saturation'])
ParamFile.write("%.3e " % Para['Fx'])
ParamFile.write("%.3e " % Para['Fy'])
ParamFile.write("%.3e\n" % Para['Fz'])
ParamFile.write("%i " % Para['InitialCondition'])
ParamFile.write("%i " % Para['BoundaryCondition'])
ParamFile.write("%.3f " % Para['din'])
ParamFile.write("%.3f\n" % Para['dout'])
ParamFile.write("%i " % Para['TimeStepMax'])
ParamFile.write("%i " % Para['Restart_interval'])
ParamFile.write("%.2e\n" % Para['tolerance'])
#ParamFile.write("%i\n" % Para['Blob_analysis_interval'])
ParamFile.close()
print "**Info: A Color.in file is created."
#end def

87
example/Python/LBPM_Preprocessing_ReadFromProcessedSegImage.py Normal file
View File

@ -0,0 +1,87 @@
#!/usr/bin/env python
import numpy as np
from netCDF4 import Dataset
from LBPM_Preprocessing_utils import *
import matplotlib.pyplot as plt
# Load the *nc image
# NOTE: for the variable name, you might want to check the *.nc file with ncdump -h
image_red = read_NetCDF_file_py('Refined_subdomain_1_raw.nc','segmented')
# This *nc file is not a raw CT image data, which has been processed by me.
# It already has the format as :
# Whereas LBPM-WIA has the convention as:
# 0 -> solid phase
# 1 -> fluid phase
(lz_red,ly_red,lx_red) = image_red.shape
# Preprocess the image
# ---------------- By default: ---------------------------------------------
# 1. The flow direction is along positive z-axis
# z_min: inlet; z_max: outlet
# 2. The boundary condition, e.g. pressure boundary condition is applied
# on x-y planes
# --------------------------------------------------------------------------
# Steps of preprocessing:
# 1. wrap boundaries along z-axis with solid nodes
# 1.1 (optional) save the raw image for morphological drainge & opeing
# 2. change the original z_min and z_max layers to reservoir layers
# 3. set z_min reservoir layer to non-wetting phase reservoir
# 4. set z_max reservoir layer to wetting phase reservoir
# 5. add porous plate layer before the wetting phase reservoir
# 6. saturate the pore space with wetting phase
# 1. wrap the boundaries
# TODO: make a variable for wall_width
image_red[:,0,:] = 0
image_red[:,ly_red-1,:] = 0
image_red[:,:,0] = 0
image_red[:,:,lx_red-1] = 0
# 1.1 (optional) Save the image for morphological analysis
#print "**Info: Save data for morphological drainage and imbibition."
#image_red.tofile('subdomain_1_reduced_morph_2x.bin') # dont use this
# 2. Saturate the medium with wetting phase
image_red[image_red==1] = 2
# 3. Set up reservoir layers
# Input parameters for the reservoir layers
NWR_len = 1 # The thickness of the non-wetting phaer reservoir
WR_len = 1 # The thickness of the wetting phase reservoir
image_red[0:NWR_len,:] = 1
image_red[-WR_len:,:] = 2
# temp test
#plt.pcolormesh(image_red[:,1,:],cmap='hot')
#plt.axis('equal')
#plt.colorbar();
#plt.show()
# 4. setup up a porous plate saturated with **wetting phase**
# Input parameters for the porous plate
pp_size = 3 #NOTE: Although it is made as an variable, pp_size of 3 is recommended
pp_len = 3 # the thickness of the porous plate
pp = 2*np.ones((pp_size,pp_size),dtype=np.int8)
#pp = np.lib.pad(pp,((1,1),(1,1)),'constant',constant_values=0)
pp = np.lib.pad(pp,((1,1),(1,1)),'constant',constant_values=0)
(ly_pp,lx_pp) = pp.shape
pp_pad = np.tile(pp,(ly_red/ly_pp,lx_red/lx_pp))
l_residual=ly_red%ly_pp
pp_pad = np.lib.pad(pp_pad,((l_residual%2,0),(l_residual%2,0)),'constant',constant_values=0)
pp_pad = np.lib.pad(pp_pad,((l_residual/2,l_residual/2),(l_residual/2,l_residual/2)),'constant',constant_values=0)
# 5. Create the simulation domain
image_red[-(WR_len+pp_len):-WR_len,:,:]=pp_pad
# temp test
#plt.pcolormesh(image_red[-5,:],cmap='hot')
#plt.axis('equal')
#plt.grid(True)
#plt.colorbar()
#plt.show()
# write to *.nc file for futher check
#print "**Info: Save data for LBPM simulations."
#write_NetCDF_file_py('subdomain_1_reduced_lbpm_2x.nc','segmented',\
# image_red,image_red.dtype.char,lx_red,ly_red,lz_red)
#image_red.tofile('subdomain_1_reduced_lbpm_2x.bin')

103
example/Python/LBPM_Preprocessing_ReadFromRawCTimageData.py Normal file
View File

@ -0,0 +1,103 @@
#!/usr/bin/env python
import numpy as np
from netCDF4 import Dataset
from LBPM_Preprocessing_utils import *
import matplotlib.pyplot as plt
# Load the *nc image
# NOTE: for the data name, you might want to check the *.nc file with ncdump -h
image = read_NetCDF_file_py('subdomain_1_from_Anna.nc','segmented')
# The loaded data from the original CT image is a masked array
# with the pore space being masked.
image = np.ma.filled(image,fill_value=0)
(lz,ly,lx) = image.shape
# NOTE: the CT image has the convention as:
# 0 -> pore space
# 1 -> solid phase
# Whereas LBPM-WIA has the convention as:
# 0 -> solid phase
# 1 -> non-wetting phase
# 2 -> wetting phase
image = np.logical_not(image)
image = image.astype(np.int8)# Now: 0 -> solid phase
# 1 -> fluid phase
image_red = image[0:256,0:256,0:256]
(lz_red,ly_red,lx_red) = image_red.shape
## 0.0. Save the original chopped image, in simple binary format.
#image_red.tofile('subdomain_1.bin')
## 0.1. Save the original chopped image, in simple *nc format
write_NetCDF_file_py('subdomain_1.nc','segmented',\
image_red,image_red.dtype.char,lx_red,ly_red,lz_red)
image_red.tofile('subdomain_1.bin')
# Preprocess the image
# ---------------- By default: ---------------------------------------------
# 1. The flow direction is along positive z-axis
# z_min: inlet; z_max: outlet
# 2. The boundary condition, e.g. pressure boundary condition is applied
# on x-y planes
# --------------------------------------------------------------------------
# Steps of preprocessing:
# 1. wrap boundaries along z-axis with solid nodes
# 1.1 (optional) save the raw image for morphological drainge & opeing
# 2. change the original z_min and z_max layers to reservoir layers
# 3. set z_min reservoir layer to non-wetting phase reservoir
# 4. set z_max reservoir layer to wetting phase reservoir
# 5. add porous plate layer before the wetting phase reservoir
# 6. saturate the pore space with wetting phase
# 1. wrap the boundaries
# TODO: make a variable for wall_width
image_red[:,0,:] = 0
image_red[:,ly_red-1,:] = 0
image_red[:,:,0] = 0
image_red[:,:,lx_red-1] = 0
# 1.1 (optional) Save the image for morphological analysis
#print "**Info: Save data for morphological drainage and imbibition."
#image_red.tofile('subdomain_1_reduced_morph.bin') # dont use this
# 2. Saturate the medium with wetting phase
image_red[image_red==1] = 2
# 3. Set up reservoir layers
# Input parameters for the reservoir layers
NWR_len = 1 # The thickness of the non-wetting phaer reservoir
WR_len = 1 # The thickness of the wetting phase reservoir
image_red[0:NWR_len,:] = 1
image_red[-WR_len:,:] = 2
# temp test
#plt.pcolormesh(image_red[:,1,:],cmap='hot')
#plt.axis('equal')
#plt.colorbar();
#plt.show()
# 4. setup up a porous plate saturated with **wetting phase**
# Input parameters for the porous plate
pp_size = 3 #NOTE: Although it is made as an variable, pp_size of 3 is recommended
pp_len = 3 # the thickness of the porous plate
pp = 2*np.ones((pp_size,pp_size),dtype=np.int8)
#pp = np.lib.pad(pp,((1,1),(1,1)),'constant',constant_values=0)
pp = np.lib.pad(pp,((1,1),(1,1)),'constant',constant_values=0)
(ly_pp,lx_pp) = pp.shape
pp_pad = np.tile(pp,(ly_red/ly_pp,lx_red/lx_pp))
l_residual=ly_red%ly_pp
pp_pad = np.lib.pad(pp_pad,((l_residual%2,0),(l_residual%2,0)),'constant',constant_values=0)
pp_pad = np.lib.pad(pp_pad,((l_residual/2,l_residual/2),(l_residual/2,l_residual/2)),'constant',constant_values=0)
# 5. Create the simulation domain
image_red[-(WR_len+pp_len):-WR_len,:,:]=pp_pad
# temp test
#plt.pcolormesh(image_red[-5,:],cmap='hot')
#plt.axis('equal')
#plt.grid(True)
#plt.colorbar()
#plt.show()
# write to *.nc file for futher check
#print "**Info: Save data for LBPM simulations."
#write_NetCDF_file_py('subdomain_1_reduced_lbpm.nc','segmented',\
# image_red,image_red.dtype.char,lx_red,ly_red,lz_red)
#image_red.tofile('subdomain_1_reduced_lbpm.bin')

41
example/Python/LBPM_Preprocessing_utils.py Normal file
View File

@ -0,0 +1,41 @@
#!/usr/bin/env python
import numpy as np
from netCDF4 import Dataset
def write_NetCDF_file_py(file_name,data_name,data_IN,data_IN_dtype,lx,ly,lz):
# NOTE: the dataIN should have the size of lz*ly*lx
data_IN.shape=(lz,ly,lx) # make sure data_IN has the right size
# open a new netCDF file for writing.
ncfile = Dataset(file_name,'w')
# create the output data.
# create the x, y and z dimensions.
ncfile.createDimension('x',lx)
ncfile.createDimension('y',ly)
ncfile.createDimension('z',lz)
# create the variable
# first argument is name of variable, second is datatype, third is
# a tuple with the names of dimensions.
data = ncfile.createVariable(data_name,data_IN_dtype,('z','y','x'))
data[:] = data_IN
# data.voxel_unit = 'micrometer'
# data.voxel_size=5.7
ncfile.close()
print '**Info: The *.nc file is written successfully !'
#end def
def read_NetCDF_file_py(file_name,data_name):
ncfile = Dataset(file_name,mode='r')
dataOUT = ncfile.variables[data_name][:]
ncfile.close()
return dataOUT
#end def

88
example/Python/LBPM_WIA_OutputData_Extraction_Preprocessing.py Normal file
View File

@ -0,0 +1,88 @@
#!/usr/bin/env python
import sys
import os
import numpy as np
from glob import glob
#from netCDF4 import Dataset
#import matplotlib.pyplot as plt
# Check if there is a proper command line argument
if len(sys.argv) !=2:
sys.stderr.write('**Error: Usage: ' + sys.argv[0] + ' <Domain.in>\n')
sys.exit()
# end if
# Read 'Domain.in' to obtain the size of 'ID.xxxxxx' and 'SignDist.xxxxx'
f = open(sys.argv[1],'r')
lines = f.readlines()
# Read processors configuration
nproc_x,nproc_y,nproc_z=np.fromstring(lines[0].splitlines()[0],dtype = np.int32,sep=' ')
# Read subdomain dimensions
lx,ly,lz = np.fromstring(lines[1].splitlines()[0],dtype = np.int32,sep=' ')
# Read domain dimensions
LX,LY,LZ = np.fromstring(lines[3].splitlines()[0],dtype = np.float32,sep=' ')
LX=int(LX)
LY=int(LY)
LZ=int(LZ)
f.close()
# The flag of writing data in binary format
# if True: write data to binary format (besides the normal *.nc format)
# if False: not write data to binary format
#flag_write_bin = int(sys.argv[2])
# read Domain.in
size_of_file=os.path.getsize('./00000') # the size of the semgented file in bytes
Size_of_Var = 8 #size_of(double) = 8 bytes
Num_of_Var = 4 # (Phase, Pressure, SignDist, BlobID)
#TODO: add input option so that user can choose what the output data they want,
# currently the output data is unnecessarily too big
# perhaps just live phase, pressure and sign_dist
offset_data_header_1 = len("Var: domain-00000-Phase: 1, 3, "+str(lx*ly*lz)+", "+str(Size_of_Var*lx*ly*lz)+", double\n")
offset_data_header_2 = len("Var: domain-00000-Pressure: 1, 3, "+str(lx*ly*lz)+", "+str(Size_of_Var*lx*ly*lz)+", double\n")
offset_data_header_3 = len("Var: domain-00000-SignDist: 1, 3, "+str(lx*ly*lz)+", "+str(Size_of_Var*lx*ly*lz)+", double\n")
offset_data_header_4 = len("Var: domain-00000-BlobID: 1, 3, "+str(lx*ly*lz)+", "+str(Size_of_Var*lx*ly*lz)+", double\n")
offset_Var_with_endian = Num_of_Var*lx*ly*lz*Size_of_Var + Num_of_Var*1
offset_pre_header = size_of_file-(offset_Var_with_endian+offset_data_header_1+offset_data_header_2+offset_data_header_3+offset_data_header_4)
data_group = glob('00*') #NOTE Might need to change this if your original image is segmented into too many (e.g. hundreds of) parts
data_group.sort()
if not data_group:
print 'Error: No data files: '+data_group
else:
for ii in range(len(data_group)):
print '**Info: Read data file: '+data_group[ii]
phase=np.memmap(data_group[ii],dtype=np.float64,\
offset=offset_pre_header+offset_data_header_1,\
shape=lx*ly*lz)
pressure=np.memmap(data_group[ii],dtype=np.float64,\
offset=offset_pre_header+offset_data_header_1+Size_of_Var*lx*ly*lz+len("\n")+offset_data_header_2,\
shape=lx*ly*lz)
blob_id=np.memmap(data_group[ii],dtype=np.float64,\
offset=offset_pre_header+offset_data_header_1+offset_data_header_2+offset_data_header_3+offset_data_header_4+3*Size_of_Var*lx*ly*lz+3*len("\n"),\
shape=lx*ly*lz)
print '**Info: Write data in normal binary format......'
phase.tofile('Phase.'+data_group[ii])
pressure.tofile('Pressure.'+data_group[ii])
blob_id.tofile('BlobID.'+data_group[ii])
print '**Info: Data extraction is completed.'
#end for
#end if

316
example/Python/LBPM_WIA_OutputData_Postprocessing_utils.py Normal file
View File

@ -0,0 +1,316 @@
#!/usr/bin/env python
#import os
#import sys
import numpy as np
#import matplotlib.pyplot as plt
from glob import glob
from netCDF4 import Dataset
def Writer_NetCDF_multiVar_LBPMWIA(file_name,data_IN,LX,LY,LZ):
#TODO: add some extra **kwargs so that users can pick a paricular variable name
# and write out what data they want to write out
#NOTE (LX,LY,LZ) is the size of a single variable
#NOTE data_IN contains 1D data with the format (Phase, Pressure, BlobID)
# open a new netCDF file for writing.
ncfile = Dataset(file_name,'w')
# create the output data.
# create the x, y and z dimensions.
ncfile.createDimension('x',LX)
ncfile.createDimension('y',LY)
ncfile.createDimension('z',LZ)
# create the variable (4 byte integer in this case)
# first argument is name of variable, second is datatype, third is
# a tuple with the names of dimensions.
Phase = ncfile.createVariable('Phase',np.float32,('z','y','x'))
Pressure = ncfile.createVariable('Pressure',np.float32,('z','y','x'))
BlobID = ncfile.createVariable('BlobID',np.float32,('z','y','x'))
Phase[:] = data_IN[:LX*LY*LZ].reshape((LZ,LY,LX))
Pressure[:] = data_IN[LX*LY*LZ:2*LX*LY*LZ].reshape((LZ,LY,LX))
BlobID[:] = data_IN[2*LX*LY*LZ:].reshape((LZ,LY,LX))
#data.voxel_unit = 'micrometer'
#data.voxel_size=5.7
ncfile.close()
print '**Info: The *.nc file is written successfully !'
#end def
def Writer_NetCDF_multiVar_Restart_LBPMWIA(file_name,data_IN,LX,LY,LZ):
#NOTE (LX,LY,LZ) is the size of a single variable
#NOTE data_IN contains 1D data with the format (Phase, Pressure, BlobID)
# open a new netCDF file for writing.
ncfile = Dataset(file_name,'w')
# Total grid size
N=LX*LY*LZ
# create the x, y and z dimensions.
ncfile.createDimension('x',LX)
ncfile.createDimension('y',LY)
ncfile.createDimension('z',LZ)
# create the variable (4 byte integer in this case)
# first argument is name of variable, second is datatype, third is
# a tuple with the names of dimensions.
den_nw=ncfile.createVariable('DensityNW',np.float32,('z','y','x'))
den_w=ncfile.createVariable('DensityW',np.float32,('z','y','x'))
f0=ncfile.createVariable('f0',np.float32,('z','y','x'))
f1=ncfile.createVariable('f1',np.float32,('z','y','x'))
f2=ncfile.createVariable('f2',np.float32,('z','y','x'))
f3=ncfile.createVariable('f3',np.float32,('z','y','x'))
f4=ncfile.createVariable('f4',np.float32,('z','y','x'))
f5=ncfile.createVariable('f5',np.float32,('z','y','x'))
f6=ncfile.createVariable('f6',np.float32,('z','y','x'))
f7=ncfile.createVariable('f7',np.float32,('z','y','x'))
f8=ncfile.createVariable('f8',np.float32,('z','y','x'))
f9=ncfile.createVariable('f9',np.float32,('z','y','x'))
f10=ncfile.createVariable('f10',np.float32,('z','y','x'))
f11=ncfile.createVariable('f11',np.float32,('z','y','x'))
f12=ncfile.createVariable('f12',np.float32,('z','y','x'))
f13=ncfile.createVariable('f13',np.float32,('z','y','x'))
f14=ncfile.createVariable('f14',np.float32,('z','y','x'))
f15=ncfile.createVariable('f15',np.float32,('z','y','x'))
f16=ncfile.createVariable('f16',np.float32,('z','y','x'))
f17=ncfile.createVariable('f17',np.float32,('z','y','x'))
f18=ncfile.createVariable('f18',np.float32,('z','y','x'))
den_nw[:] = data_IN[ :1*N].reshape(LZ,LY,LX)
den_w[:] = data_IN[1*N:2*N].reshape(LZ,LY,LX)
f0[:] = data_IN[2*N:3*N].reshape(LZ,LY,LX)
f1[:] = data_IN[3*N:4*N].reshape(LZ,LY,LX)
f2[:] = data_IN[4*N:5*N].reshape(LZ,LY,LX)
f3[:] = data_IN[5*N:6*N].reshape(LZ,LY,LX)
f4[:] = data_IN[6*N:7*N].reshape(LZ,LY,LX)
f5[:] = data_IN[7*N:8*N].reshape(LZ,LY,LX)
f6[:] = data_IN[8*N:9*N].reshape(LZ,LY,LX)
f7[:] = data_IN[9*N:10*N].reshape(LZ,LY,LX)
f8[:] = data_IN[10*N:11*N].reshape(LZ,LY,LX)
f9[:] = data_IN[11*N:12*N].reshape(LZ,LY,LX)
f10[:] = data_IN[12*N:13*N].reshape(LZ,LY,LX)
f11[:] = data_IN[13*N:14*N].reshape(LZ,LY,LX)
f12[:] = data_IN[14*N:15*N].reshape(LZ,LY,LX)
f13[:] = data_IN[15*N:16*N].reshape(LZ,LY,LX)
f14[:] = data_IN[16*N:17*N].reshape(LZ,LY,LX)
f15[:] = data_IN[17*N:18*N].reshape(LZ,LY,LX)
f16[:] = data_IN[18*N:19*N].reshape(LZ,LY,LX)
f17[:] = data_IN[19*N:20*N].reshape(LZ,LY,LX)
f18[:] = data_IN[20*N:21*N].reshape(LZ,LY,LX)
ncfile.close()
print '**Info: The *.nc file is written successfully !'
#end def
def Writer_NetCDF_singleVar_LBPMWIA(file_name,data_name,data_IN,LX,LY,LZ):
#NOTE (LX,LY,LZ) is the size of a single variable
#NOTE data_IN contains 1D data with the format (Phase, Pressure, BlobID)
# open a new netCDF file for writing.
ncfile = Dataset(file_name,'w')
# create the output data.
# create the x, y and z dimensions.
ncfile.createDimension('x',LX)
ncfile.createDimension('y',LY)
ncfile.createDimension('z',LZ)
# create the variable (4 byte integer in this case)
# first argument is name of variable, second is datatype, third is
# a tuple with the names of dimensions.
data = ncfile.createVariable(data_name,np.float32,('z','y','x'))
data[:] = data_IN.reshape((LZ,LY,LX))
#data.voxel_unit = 'micrometer'
#data.voxel_size=5.7
ncfile.close()
print '**Info: The *.nc file is written successfully !'
#end def
def Reconstruct_single_data(data_name_group,proc_config,grid_config):
#NOTE: This function is not used for reconstructing ID.* and SignDist.*
#data_name_group is a string, e.g. 'Phase.00*'
#proc_config is a tuple with (nproc_x,nproc_y,nproc_z), the processor configuration
#NOTE: grid_config is a tuple with (lx,ly,lz) is the size of a segmented domain
data_group=glob(data_name_group)
data_group.sort()
nproc_x, nproc_y, nproc_z = proc_config
lx,ly,lz = grid_config
print '**Info: data reconstruction: '+data_group[0]
for z_axis in np.arange(nproc_z):
for y_axis in np.arange(nproc_y):
for x_axis in np.arange(0,nproc_x,2):
if nproc_x%2==0: # nproc_x is an even number
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.float64)
temp_x2 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis+1],dtype=np.float64)
temp_x1.shape = (lz,ly,lx)
temp_x2.shape = (lz,ly,lx)
#temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]# for ID.* and SignDist.* which have one extra buffer layer around the data
#temp_x2 = temp_x2[1:lz-1,1:ly-1,1:lx-1]
temp_x1 = np.concatenate((temp_x1,temp_x2),axis=2)
else: # nproc_x is an odd number
if nproc_x==1:
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.float64)
temp_x1.shape = (lz,ly,lx)
elif x_axis < nproc_x-2:
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.float64)
temp_x2 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis+1],dtype=np.float64)
temp_x1.shape = (lz,ly,lx)
temp_x2.shape = (lz,ly,lx)
#temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]
#temp_x2 = temp_x2[1:lz-1,1:ly-1,1:lx-1]
temp_x1 = np.concatenate((temp_x1,temp_x2),axis=2)
#end if
#end if
if x_axis == 0:
ID_x_temp = temp_x1.copy()
else:
ID_x_temp = np.concatenate((ID_x_temp,temp_x1),axis=2)
#end if
#end for
if nproc_x !=1 and nproc_x%2 != 0:
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+nproc_x-1],dtype=np.float64)
ID_x_temp = np.concatenate((ID_x_temp,temp_x1),axis=2)
#end if
if y_axis == 0:
ID_y_temp = ID_x_temp.copy()
else:
ID_y_temp = np.concatenate((ID_y_temp,ID_x_temp),axis=1)
#end if
#end for
if z_axis == 0:
ID_z_temp = ID_y_temp.copy()
else:
ID_z_temp = np.concatenate((ID_z_temp,ID_y_temp),axis=0)
#end if
#end for
print '**Info: data reconstruction completed for '+data_group[0]
return ID_z_temp.flatten()
#end def
def Reconstruct_SignDist_data(data_name_group,proc_config,grid_config):
#NOTE: This function is only used for reconstructing ID.* and SignDist.*
#data_name_group is a string, e.g. 'ID.00*'
#proc_config is a tuple with (nproc_x,nproc_y,nproc_z), the processor configuration
#NOTE: grid_config is a tuple with (lx,ly,lz) is the size of a segmented domain
data_group=glob(data_name_group)
data_group.sort()
nproc_x, nproc_y, nproc_z = proc_config
lx,ly,lz = grid_config
print '**Info: data reconstruction: '+data_group[0]
for z_axis in np.arange(nproc_z):
for y_axis in np.arange(nproc_y):
for x_axis in np.arange(0,nproc_x,2):
if nproc_x%2==0: # nproc_x is an even number
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.float64)
temp_x2 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis+1],dtype=np.float64)
temp_x1.shape = (lz,ly,lx)
temp_x2.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]# for ID.* and SignDist.* which have one extra buffer layer around the data
temp_x2 = temp_x2[1:lz-1,1:ly-1,1:lx-1]
temp_x1 = np.concatenate((temp_x1,temp_x2),axis=2)
else: # nproc_x is an odd number
if nproc_x==1:
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.float64)
temp_x1.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]
elif x_axis < nproc_x-2: # i.e. nproc_x = 3 or 5 or 7 ...
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.float64)
temp_x2 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis+1],dtype=np.float64)
temp_x1.shape = (lz,ly,lx)
temp_x2.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]
temp_x2 = temp_x2[1:lz-1,1:ly-1,1:lx-1]
temp_x1 = np.concatenate((temp_x1,temp_x2),axis=2)
#end if
#end if
if x_axis == 0:
ID_x_temp = temp_x1.copy()
else:
ID_x_temp = np.concatenate((ID_x_temp,temp_x1),axis=2)
#end if
#end for
if nproc_x !=1 and nproc_x%2 != 0:
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+nproc_x-1],dtype=np.float64)
ID_x_temp = np.concatenate((ID_x_temp,temp_x1),axis=2)
#end if
if y_axis == 0:
ID_y_temp = ID_x_temp.copy()
else:
ID_y_temp = np.concatenate((ID_y_temp,ID_x_temp),axis=1)
#end if
#end for
if z_axis == 0:
ID_z_temp = ID_y_temp.copy()
else:
ID_z_temp = np.concatenate((ID_z_temp,ID_y_temp),axis=0)
#end if
#end for
print '**Info: Signed distance data reconstruction completed for '+data_group[0]
return ID_z_temp.flatten()
#end def
def Reconstruct_ID_data(data_name_group,proc_config,grid_config):
#NOTE: This function is only used for reconstructing ID.* and SignDist.*
#data_name_group is a string, e.g. 'ID.00*'
#proc_config is a tuple with (nproc_x,nproc_y,nproc_z), the processor configuration
#NOTE: grid_config is a tuple with (lx,ly,lz) is the size of a segmented domain
data_group=glob(data_name_group)
data_group.sort()
nproc_x, nproc_y, nproc_z = proc_config
lx,ly,lz = grid_config
print '**Info: data reconstruction: '+data_group[0]
for z_axis in np.arange(nproc_z):
for y_axis in np.arange(nproc_y):
for x_axis in np.arange(0,nproc_x,2):
if nproc_x%2==0: # nproc_x is an even number
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.int8)
temp_x2 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis+1],dtype=np.int8)
temp_x1.shape = (lz,ly,lx)
temp_x2.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]# for ID.* and SignDist.* which have one extra buffer layer around the data
temp_x2 = temp_x2[1:lz-1,1:ly-1,1:lx-1]
temp_x1 = np.concatenate((temp_x1,temp_x2),axis=2)
else: # nproc_x is an odd number
if nproc_x==1:
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.int8)
temp_x1.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]
elif x_axis < nproc_x-2: # i.e. nproc_x = 3 or 5 or 7 ...
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis],dtype=np.int8)
temp_x2 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+x_axis+1],dtype=np.int8)
temp_x1.shape = (lz,ly,lx)
temp_x2.shape = (lz,ly,lx)
temp_x1 = temp_x1[1:lz-1,1:ly-1,1:lx-1]
temp_x2 = temp_x2[1:lz-1,1:ly-1,1:lx-1]
temp_x1 = np.concatenate((temp_x1,temp_x2),axis=2)
#end if
#end if
if x_axis == 0:
ID_x_temp = temp_x1.copy()
else:
ID_x_temp = np.concatenate((ID_x_temp,temp_x1),axis=2)
#end if
#end for
if nproc_x !=1 and nproc_x%2 != 0:
temp_x1 = np.fromfile(data_group[z_axis*nproc_y*nproc_x+y_axis*nproc_x+nproc_x-1],dtype=np.int8)
ID_x_temp = np.concatenate((ID_x_temp,temp_x1),axis=2)
#end if
if y_axis == 0:
ID_y_temp = ID_x_temp.copy()
else:
ID_y_temp = np.concatenate((ID_y_temp,ID_x_temp),axis=1)
#end if
#end for
if z_axis == 0:
ID_z_temp = ID_y_temp.copy()
else:
ID_z_temp = np.concatenate((ID_z_temp,ID_y_temp),axis=0)
#end if
#end for
print '**Info: Phase ID data reconstruction completed for '+data_group[0]
return ID_z_temp.flatten()
#end def

65
example/Python/LBPM_WIA_OutputData_Reconstruction.py Normal file
View File

@ -0,0 +1,65 @@
#!/usr/bin/env python
import sys
import os
import numpy as np
#import matplotlib.pyplot as plt
from glob import glob
#from netCDF4 import Dataset
from LBPM_WIA_OutputData_Postprocessing_utils import *
# Check if there is a proper command line argument
if len(sys.argv) !=2:
sys.stderr.write('**Error: Usage: ' + sys.argv[0] + ' <Domain.in>\n')
sys.exit()
# end if
# Read 'Domain.in' --------------------------------------------------------------------------
f = open(sys.argv[1],'r')
lines = f.readlines()
# Read processors configuration
nproc_x,nproc_y,nproc_z=np.fromstring(lines[0].splitlines()[0],dtype = np.int32,sep=' ')
# Read subdomain dimensions
lx,ly,lz = np.fromstring(lines[1].splitlines()[0],dtype = np.int32,sep=' ')
# Read domain dimensions
LX,LY,LZ = np.fromstring(lines[3].splitlines()[0],dtype = np.float32,sep=' ')
LX=int(LX)
LY=int(LY)
LZ=int(LZ)
f.close()
# -------------------------------------------------------------------------------------------
phase=Reconstruct_single_data('Phase.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))
pressure=Reconstruct_single_data('Pressure.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))
blob_id=Reconstruct_single_data('BlobID.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))
output_data=np.hstack((phase,pressure,blob_id))
#output_data_name_bin='Output_data.bin'
folder_path=os.getcwd()
output_data_name_nc ='OutputData_'+os.path.basename(folder_path)+'.nc'
# Write to binary file
#output_data.tofile(output_data_name_bin)
#print '**Info: Output data \''+output_data_name_bin+'\' is written.'
# Write to *nc file
Writer_NetCDF_multiVar_LBPMWIA(output_data_name_nc,output_data.astype(np.float32),LX,LY,LZ)
# Clean up temporary files
files=glob('Phase.0*')
for file_tmp in files:
print "**Warning: Remove file: "+file_tmp
os.remove(file_tmp)
#end for
files=glob('Pressure.0*')
for file_tmp in files:
print "**Warning: Remove file: "+file_tmp
os.remove(file_tmp)
#end for
files=glob('BlobID.0*')
for file_tmp in files:
print "**Warning: Remove file: "+file_tmp
os.remove(file_tmp)
#end for

73
example/Python/LBPM_WIA_Restart_Extraction_Preprocessing.py Normal file
View File

@ -0,0 +1,73 @@
#!/usr/bin/env python
import sys
import os
import numpy as np
from glob import glob
#from netCDF4 import Dataset
#import matplotlib.pyplot as plt
# Check if there is a proper command line argument
if len(sys.argv) !=2:
sys.stderr.write('**Error: Usage: ' + sys.argv[0] + ' <Domain.in>\n')
sys.exit()
# end if
# Read 'Domain.in' to obtain the size of 'ID.xxxxxx' and 'SignDist.xxxxx'
f = open(sys.argv[1],'r')
lines = f.readlines()
# Read processors configuration
nproc_x,nproc_y,nproc_z=np.fromstring(lines[0].splitlines()[0],dtype = np.int32,sep=' ')
# Read subdomain dimensions
lx,ly,lz = np.fromstring(lines[1].splitlines()[0],dtype = np.int32,sep=' ')
# Read domain dimensions
LX,LY,LZ = np.fromstring(lines[3].splitlines()[0],dtype = np.float32,sep=' ')
LX=int(LX)
LY=int(LY)
LZ=int(LZ)
f.close()
# The total size of the subdomain including the halo layer
N=(lx+2)*(ly+2)*(lz+2)
# Read 'Restart.00*' binary files
# NOTE: Currently there are in total 21 sets of data in the Restart files
# They are: { rho_nw, rho_w, f_Even, f_odd }
# f_even -> { f0, f2, f4, f6, f8, f10, f12, f14, f16, f18 }
# f_odd -> { f1, f3, f5, f7, f9, f11, f13, f15, f17 }
data_group = glob('Restart.00*')
data_group.sort()
if not data_group:
print 'Error: No data files: '+data_group
else:
for ii in range(len(data_group)):
print '**Info: Read data file: '+data_group[ii]
Restart=np.fromfile(data_group[ii],dtype=np.float64)
Restart.shape = (N,21)
print '**Info: Write data in normal binary format......'
Restart[:,0].tofile('DenNW.'+data_group[ii][len('Restart.'):])
Restart[:,1].tofile('DenW.'+data_group[ii][len('Restart.'):])
Restart[:,2].tofile('f0.'+data_group[ii][len('Restart.'):])
Restart[:,3].tofile('f2.'+data_group[ii][len('Restart.'):])
Restart[:,4].tofile('f4.'+data_group[ii][len('Restart.'):])
Restart[:,5].tofile('f6.'+data_group[ii][len('Restart.'):])
Restart[:,6].tofile('f8.'+data_group[ii][len('Restart.'):])
Restart[:,7].tofile('f10.'+data_group[ii][len('Restart.'):])
Restart[:,8].tofile('f12.'+data_group[ii][len('Restart.'):])
Restart[:,9].tofile('f14.'+data_group[ii][len('Restart.'):])
Restart[:,10].tofile('f16.'+data_group[ii][len('Restart.'):])
Restart[:,11].tofile('f18.'+data_group[ii][len('Restart.'):])
Restart[:,12].tofile('f1.'+data_group[ii][len('Restart.'):])
Restart[:,13].tofile('f3.'+data_group[ii][len('Restart.'):])
Restart[:,14].tofile('f5.'+data_group[ii][len('Restart.'):])
Restart[:,15].tofile('f7.'+data_group[ii][len('Restart.'):])
Restart[:,16].tofile('f9.'+data_group[ii][len('Restart.'):])
Restart[:,17].tofile('f11.'+data_group[ii][len('Restart.'):])
Restart[:,18].tofile('f13.'+data_group[ii][len('Restart.'):])
Restart[:,19].tofile('f15.'+data_group[ii][len('Restart.'):])
Restart[:,20].tofile('f17.'+data_group[ii][len('Restart.'):])
print '**Info: Data extraction is completed.'
#end for
#end if

81
example/Python/LBPM_WIA_Restart_Reconstruction.py Normal file
View File

@ -0,0 +1,81 @@
#!/usr/bin/env python
import sys
import os
import numpy as np
#import matplotlib.pyplot as plt
from glob import glob
#from netCDF4 import Dataset
from LBPM_WIA_OutputData_Postprocessing_utils import *
# Check if there is a proper command line argument
if len(sys.argv) !=2:
sys.stderr.write('**Error: Usage: ' + sys.argv[0] + ' <Domain.in>\n')
sys.exit()
# end if
# Read 'Domain.in' --------------------------------------------------------------------------
f = open(sys.argv[1],'r')
lines = f.readlines()
# Read processors configuration
nproc_x,nproc_y,nproc_z=np.fromstring(lines[0].splitlines()[0],dtype = np.int32,sep=' ')
# Read subdomain dimensions
lx,ly,lz = np.fromstring(lines[1].splitlines()[0],dtype = np.int32,sep=' ')
# Read domain dimensions
LX,LY,LZ = np.fromstring(lines[3].splitlines()[0],dtype = np.float32,sep=' ')
LX=int(LX)
LY=int(LY)
LZ=int(LZ)
f.close()
# NOTE: Do not forget the halo layer attached on each subdomain !
lx = lx + 2
ly = ly + 2
lz = lz + 2
# -------------------------------------------------------------------------------------------
print "**Info: Reconstructing individual data sets according to the processor grid..."
dataOUT = Reconstruct_SignDist_data('DenNW.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('DenW.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f0.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f1.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f2.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f3.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f4.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f5.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f6.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f7.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f8.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f9.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f10.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f11.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f12.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f13.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f14.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f15.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f16.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f17.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
dataOUT = np.hstack((dataOUT,Reconstruct_SignDist_data('f18.0*',(nproc_x,nproc_y,nproc_z),(lx,ly,lz))))
#output_data_name_bin='Output_data.bin'
output_data_name_nc ='Restart_dbg_TimeStep_.nc'
# Write to binary file
#output_data.tofile(output_data_name_bin)
#print '**Info: Output data \''+output_data_name_bin+'\' is written.'
# Write to *nc file
Writer_NetCDF_multiVar_Restart_LBPMWIA(output_data_name_nc,dataOUT.astype(np.float32),LX,LY,LZ)
# Clean up temporary files
files=glob('Den*')
for file_tmp in files:
print "**Warning: Remove file: "+file_tmp
os.remove(file_tmp)
#end for
files=glob('f*.0*')
for file_tmp in files:
print "**Warning: Remove file: "+file_tmp
os.remove(file_tmp)
#end for

119
example/Python/lbpm_morph_generate_configs.py Normal file
View File

@ -0,0 +1,119 @@
#!/usr/bin/env python
import numpy as np
import matplotlib.pylab as plt
from glob import glob
import sys
#from dist_func_utils import *
import scipy.stats as stats
import scipy.ndimage.morphology as morphology
import scipy.ndimage.measurements as measurements
from morphological_analysis_utils import *
# Check if there is a proper command line argument
#if len(sys.argv) !=2:
# sys.stderr.write('Usage: ' + sys.argv[0] + ' <Base name of the segmented image data>\n')
# sys.exit()
# end if
base_name = 'MicroModel' # Give the name of the input parameter file
#base_name = 'SignDist' # Give the name of the input parameter file
#base_name = sys.argv[1] # Give the name of the input parameter file
SegImage_format = '.raw'
# Load a group of segmented image data files with 'base_name' in the names of files
# e.g. 'Micromodel_1_segmented.raw' etc.
file_input_group = glob('*'+base_name+'*'+SegImage_format) # need to match the image data format
# Dimensions for segmented image
lz = 500
ly = 500
lx = 12
#lz = 200
#ly = 200
#lx = 200
R_critical = 4.5
# Process all imported experimental images
if not file_input_group:
print 'Error: Input files cannot be found ! '
else:
for ii in range(len(file_input_group)):
file_input_single = file_input_group[ii]
print "Analyzing segmented image data "+file_input_single+" now..."
# The imported data is a double-precision signed distance map
image_raw = np.fromfile(file_input_single,dtype=np.uint8)
#image_raw = np.fromfile(file_input_single,dtype=np.float64)
image_raw.shape=(lz,ly,lx)
edt_dist = morphology.distance_transform_edt(image_raw)
#edt_dist = edt_dist[:,:,5]
# wrap the medium with one layer of solid phase
#image_raw = np.lib.pad(image_raw,((0,0),(1,1),(1,1)),'constant',constant_values=0.0)
drain_config,Sw_drain = generate_morph_drain_config(edt_dist,R_critical)
imbib_config,Sw_imbib = generate_morph_imbib_config(edt_dist,R_critical)
#drain_config,Sw_drain = generate_morph_drain_config(image_raw,R_critical)
#imbib_config,Sw_imbib = generate_morph_imbib_config(image_raw,R_critical)
#Sw_drain = generate_morph_drain_curv_3D(image_raw,R_critical)
#imbib_config,Sw_imbib = generate_morph_imbib_config2_3D(image_raw,R_critical)
#Sw_imbib = generate_morph_imbib_curv2_3D(image_raw,R_critical)
print 'Rcrit: '+str(R_critical)+', Drain_Sw = '+str(Sw_drain)
print 'Rcrit: '+str(R_critical)+', Imbib_Sw = '+str(Sw_imbib)
# Save the configuration files
# drain_config.tofile('Drain_config_Rcrit_'+str(R_critical)+'.raw')
# imbib_config.tofile('Imbib_config_Rcrit_'+str(R_critical)+'.raw')
# Load the saved data
# drain_config = np.fromfile('Drain_config_Rcrit_'+str(R_critical)+'.raw',dtype=np.uint8)
# imbib_config = np.fromfile('Imbib_config_Rcrit_'+str(R_critical)+'.raw',dtype=np.uint8)
# drain_config.shape = (lz,ly,lx)
# imbib_config.shape = (lz,ly,lx)
plt.figure(1)
plt.subplot(1,2,1)
plt.title('Drainage: Rcrit='+str(R_critical)+' Sw='+str(Sw_drain))
plt.pcolormesh(drain_config[:,:,lx/2],cmap='hot')
#plt.pcolormesh(drain_config[20,:,:],cmap='hot')
plt.axis('equal')
plt.grid(True)
plt.subplot(1,2,2)
plt.title('Imbibition: Rcrit='+str(R_critical)+' Sw='+str(Sw_imbib))
plt.pcolormesh(imbib_config[:,:,lx/2],cmap='hot')
#plt.pcolormesh(imbib_config[20,:,:],cmap='hot')
plt.axis('equal')
plt.grid(True)
# plt.figure(1)
# plt.subplot(1,2,1)
# plt.title('Drainage: Rcrit='+str(R_critical)+' Sw='+str(Sw_drain))
# plt.pcolormesh(drain_config,cmap='hot')
# plt.axis('equal')
# plt.grid(True)
#
# plt.subplot(1,2,2)
# plt.title('Imbibition: Rcrit='+str(R_critical)+' Sw='+str(Sw_imbib))
# plt.pcolormesh(imbib_config,cmap='hot')
# plt.axis('equal')
# plt.grid(True)
plt.figure(2)
plt.plot(Sw_drain,1.0/R_critical,'ro',markersize=6,label='Drainage')
plt.plot(Sw_imbib,1.0/R_critical,'b^',markersize=6,label='Imbibition')
plt.legend(loc='best')
plt.grid(True)
plt.xlim([0,1.0])
plt.ylim([0,2.0])
plt.show()
#end if
#end for

399
example/Python/lbpm_morph_utils.py Normal file
View File

@ -0,0 +1,399 @@
#!/usr/bin/env python
# TODO: double check the issue of the "view of the array" for all the following functions
# make sure it is the copy of the array, not the view of the array is used for any
# subprocesses.
import numpy as np
import scipy.stats as stats
import scipy.ndimage.morphology as morphology
import scipy.ndimage.measurements as measurements
#from dist_func_utils import *
def generate_morph_imbib_config(seg_image_input,R_critical):
# This funciton works well at high-Pc-low-Sw
# There is no retarded Sw when the curvature is high
# The method of this functions:
# 1. Perform erosion with the radius of R_critical
# 2. Perform dilation with the radius of R_critical
if seg_image_input.ndim ==2 :
seg_image = seg_image_input>0.0
pore_vol = 1.0*seg_image.sum()
radius = R_critical
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion with radius of R_critical
seg_image_ero = morphology.binary_erosion(seg_image,structure=circle,border_value=1)
# NOTE: Be careful with the 'border_value' of the erosion - should be 'border_value=1'
# Step 2: Perform dilation with radius of R_critical
seg_image_dil = morphology.binary_dilation(seg_image_ero,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of the array after dilation is 'bool'
seg_image_dil = seg_image_dil.astype(np.uint8)
seg_image_dil[np.logical_not(seg_image_dil)]=2
seg_image_dil[np.logical_not(seg_image)] = 0
Sw = (seg_image_dil==2).sum()/pore_vol
else: # 3D
seg_image = seg_image_input>0.0
pore_vol = 1.0*seg_image.sum()
radius = R_critical
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 + (grid[2]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion with radius of R_critical
seg_image_ero = morphology.binary_erosion(seg_image,structure=circle,border_value=1)
# NOTE: Be careful with the 'border_value' of the erosion - should be 'border_value=1'
# Step 2: Perform dilation with radius of R_critical
seg_image_dil = morphology.binary_dilation(seg_image_ero,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of the array after dilation is 'bool'
seg_image_dil = seg_image_dil.astype(np.uint8)
seg_image_dil[np.logical_not(seg_image_dil)]=2
seg_image_dil[np.logical_not(seg_image)] = 0
Sw = (seg_image_dil==2).sum()/pore_vol
#end if
return (seg_image_dil,Sw)
#end def
def generate_morph_imbib_curv(seg_image_input,R_critical):
# This funciton works well at high-Pc-low-Sw
# There is no retarded Sw when the curvature is high
# The method of this functions:
# 1. Perform erosion with the radius of R_critical
# 2. Perform dilation with the radius of R_critical
# *********************************************************************
# Input: seg_image: a well shaped segmented image with size (lz,ly,lx)
# seg_image has values as : NW phase -> 1
# W phase -> 2
# solid phase -> 0
# *********************************************************************
if seg_image_input.ndim == 2:
pore_vol = 1.0*(seg_image_input>0.0).sum()
radius = R_critical
print 'Morphological Imbibition: processing critical radius: '+str(radius)+' now......'
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion with radius of R_critical
seg_image_ero = morphology.binary_erosion(seg_image_input>0.0,structure=circle,border_value=1)
# NOTE: Be careful with the 'border_value' of the erosion - should be 'border_value=1'
# Step 2: Perform dilation with radius of R_critical
seg_image_ero_dil = morphology.binary_dilation(seg_image_ero,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of the array after dilation is 'bool'
seg_image_ero_dil[seg_image_input<=0.0]=False
Sw = 1.0 - seg_image_ero_dil.sum()/pore_vol
else:
pore_vol = 1.0*(seg_image_input>0.0).sum()
radius = R_critical
print 'Morphological Imbibition: processing critical radius: '+str(radius)+' now......'
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 + (grid[2]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion with radius of R_critical
seg_image_ero = morphology.binary_erosion(seg_image_input>0.0,structure=circle,border_value=1)
# NOTE: Be careful with the 'border_value' of the erosion - should be 'border_value=1'
# Step 2: Perform dilation with radius of R_critical
seg_image_ero_dil = morphology.binary_dilation(seg_image_ero,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of the array after dilation is 'bool'
seg_image_ero_dil[seg_image_input<=0.0]=False
Sw = 1.0 - seg_image_ero_dil.sum()/pore_vol
#end if
return Sw
#end def
def generate_morph_drain_config(seg_image_input,R_critical):
# The method for this function follows Hilper & Miller AWR(2001)
# 1. Perform erosion for the pore space with radius of R_critical
# 2. Label the eroded pore space, and leave only the pore space that is still
# connected with the non-wetting phase reservoir
# 3. Perform the dilation for the labelled pore space with radius of R_critical
# **************************************************************************
# Input: seg_image: a well shaped segmented image with size (lz,ly,lx)
# seg_image has values as : NW phase -> 1
# W phase -> 2
# solid phase -> 0
# **************************************************************************
if seg_image_input.ndim ==2:
seg_image = seg_image_input>0.0
pore_vol = 1.0*seg_image.sum()
radius = R_critical
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion on the pore space
# NOTE: the dtype of 'seg_im_ero' is 'bool'
seg_im_ero = morphology.binary_erosion(seg_image,structure=circle,border_value=1)
# NOTE: 'border_value' for erosion should be 'True'
# Step 2: Label the eroded pore space
# NOTE: Assume the NW phase reservoir is at the first layer of the domain
# i.e. at seg_image[0,:] - adjust it if this does not suit your need
# For erosion, assume that diagonals are not considered
# For erosion, assume that diagonals are not considered
seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(2,1))
#seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(2,2))
# NOTE: Here I assume the inlet is at the first layer of the array's axis=2 (i.e. domain[0,:,:])\
# You can always change to any other layers as the inlet for this drainage.
label_check = seg_im_ero_label_temp[0,seg_im_ero_label_temp[0,:]!=0]
label_check = np.unique(label_check)
# NOTE the following lines are only for your to check things
# ******************** For check *******************************#
# It assign the labelled array as: NW -> 1, W -> 2, Solid -> 0
#seg_im_ero_label_show = seg_im_ero_label.copy()
#seg_im_ero_label_show[seg_im_ero_label_show !=1] = 2
#seg_im_ero_label_show[np.logical_not(seg_image_2d)]=0
# ******************** End: for check **************************#
seg_im_ero_label = np.zeros_like(seg_im_ero_label_temp,dtype=bool)
for labels in label_check:
seg_im_ero_label = np.logical_or(seg_im_ero_label,seg_im_ero_label_temp==labels)
seg_im_ero_label = seg_im_ero_label.astype(np.uint8)
# Step 3: perform dilation on the labelled pore space
seg_im_ero_label_dil = morphology.binary_dilation(seg_im_ero_label,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of 'seg_im_ero_label_dil' is 'bool'
seg_im_ero_label_dil = seg_im_ero_label_dil.astype(np.uint8)
seg_im_ero_label_dil[np.logical_not(seg_im_ero_label_dil)]=2
seg_im_ero_label_dil[np.logical_not(seg_image)]=0
Sw = (seg_im_ero_label_dil==2).sum()/pore_vol
else: # 3D porous medium
seg_image = seg_image_input>0.0
pore_vol = 1.0*seg_image.sum()
radius = R_critical
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 + (grid[2]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion on the pore space
# NOTE: the dtype of 'seg_im_ero' is 'bool'
seg_im_ero = morphology.binary_erosion(seg_image,structure=circle,border_value=1)
# NOTE: 'border_value' for erosion should be 'True'
# Step 2: Label the eroded pore space
# NOTE: Assume the NW phase reservoir is at the first layer of the domain
# i.e. at seg_image[0,:] - adjust it if this does not suit your need
# For erosion, assume that diagonals are not considered
seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(3,1))
#seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(3,3))
# NOTE: Here I assume the inlet is at the first layer of the array's axis=2 (i.e. domain[0,:,:])\
# You can always change to any other layers as the inlet for this drainage.
label_check = seg_im_ero_label_temp[0,seg_im_ero_label_temp[0,:]!=0]
label_check = np.unique(label_check)
# NOTE the following lines are only for your to check things
# ******************** For check *******************************#
# It assign the labelled array as: NW -> 1, W -> 2, Solid -> 0
#seg_im_ero_label_show = seg_im_ero_label.copy()
#seg_im_ero_label_show[seg_im_ero_label_show !=1] = 2
#seg_im_ero_label_show[np.logical_not(seg_image_2d)]=0
# ******************** End: for check **************************#
seg_im_ero_label = np.zeros_like(seg_im_ero_label_temp,dtype=bool)
for labels in label_check:
seg_im_ero_label = np.logical_or(seg_im_ero_label,seg_im_ero_label_temp==labels)
seg_im_ero_label = seg_im_ero_label.astype(np.uint8)
# Step 3: perform dilation on the labelled pore space
seg_im_ero_label_dil = morphology.binary_dilation(seg_im_ero_label,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of 'seg_im_ero_label_dil' is 'bool'
seg_im_ero_label_dil = seg_im_ero_label_dil.astype(np.uint8)
seg_im_ero_label_dil[np.logical_not(seg_im_ero_label_dil)]=2
seg_im_ero_label_dil[np.logical_not(seg_image)]=0
Sw = (seg_im_ero_label_dil==2).sum()/pore_vol
#end if
return (seg_im_ero_label_dil,Sw)
#end def
def generate_morph_drain_curv(seg_image_input,R_critical):
# The method for this function follows Hilper & Miller AWR(2001)
# 1. Perform erosion for the pore space with radius of R_critical
# 2. Label the eroded pore space, and leave only the pore space that is still
# connected with the non-wetting phase reservoir
# 3. Perform the dilation for the labelled pore space with radius of R_critical
# ****************************************************************************
# Currently I am provided with a 3D SignDist image which has positive values
# in the pore space and 0.0 in the solid phase.
# ****************************************************************************
if seg_image_input.ndim == 2:
pore_vol = 1.0*(seg_image_input>0.0).sum()
radius = R_critical
print 'Morphological Drainage: processing critical radius: '+str(radius)+' now......'
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion on the pore space
# NOTE: the dtype of 'seg_im_ero' is 'bool'
seg_im_ero = morphology.binary_erosion(seg_image_input>0.0,structure=circle,border_value=1)
# NOTE: 'border_value' for erosion should be 'True'
# Step 2: Label the eroded pore space
# NOTE: Assume the NW phase reservoir is at the first layer of the domain
# i.e. at seg_image[0,:] - adjust it if this does not suit your need
# For erosion, assume that diagonals are not considered
seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(2,1))
#seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(2,2))
# NOTE: Here I assume the inlet is at the first layer of the array's axis=2 (i.e. domain[0,:,:])\
# You can always change to any other layers as the inlet for this drainage.
label_check = seg_im_ero_label_temp[0,seg_im_ero_label_temp[0,:]!=0]
label_check = np.unique(label_check)
# NOTE the following lines are only for your to check things
# ******************** For check *******************************#
# It assign the labelled array as: NW -> 1, W -> 2, Solid -> 0
#seg_im_ero_label_show = seg_im_ero_label.copy()
#seg_im_ero_label_show[seg_im_ero_label_show !=1] = 2
#seg_im_ero_label_show[np.logical_not(seg_image_2d)]=0
# ******************** End: for check **************************#
seg_im_ero_label = np.zeros_like(seg_im_ero_label_temp,dtype=bool)
for labels in label_check:
seg_im_ero_label = np.logical_or(seg_im_ero_label,seg_im_ero_label_temp==labels)
#seg_im_ero_label = seg_im_ero_label.astype(np.uint8)
# Step 3: perform dilation on the labelled pore space
seg_im_ero_label_dil = morphology.binary_dilation(seg_im_ero_label,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of 'seg_im_ero_label_dil' is 'bool'
seg_im_ero_label_dil[seg_image_input<=0.0]=False
Sw = 1.0 - seg_im_ero_label_dil.sum()/pore_vol
else:
pore_vol = 1.0*(seg_image_input>0.0).sum()
radius = R_critical
print 'Morphological Drainage: processing critical radius: '+str(radius)+' now......'
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 + (grid[2]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion on the pore space
# NOTE: the dtype of 'seg_im_ero' is 'bool'
seg_im_ero = morphology.binary_erosion(seg_image_input>0.0,structure=circle,border_value=1)
# NOTE: 'border_value' for erosion should be 'True'
# Step 2: Label the eroded pore space
# NOTE: Assume the NW phase reservoir is at the first layer of the domain
# i.e. at seg_image[0,:] - adjust it if this does not suit your need
# For erosion, assume that diagonals are not considered
seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(3,1))
#seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(3,3))
# NOTE: Here I assume the inlet is at the first layer of the array's axis=2 (i.e. domain[0,:,:])\
# You can always change to any other layers as the inlet for this drainage.
label_check = seg_im_ero_label_temp[0,seg_im_ero_label_temp[0,:]!=0]
label_check = np.unique(label_check)
# NOTE the following lines are only for your to check things
# ******************** For check *******************************#
# It assign the labelled array as: NW -> 1, W -> 2, Solid -> 0
#seg_im_ero_label_show = seg_im_ero_label.copy()
#seg_im_ero_label_show[seg_im_ero_label_show !=1] = 2
#seg_im_ero_label_show[np.logical_not(seg_image_2d)]=0
# ******************** End: for check **************************#
seg_im_ero_label = np.zeros_like(seg_im_ero_label_temp,dtype=bool)
for labels in label_check:
seg_im_ero_label = np.logical_or(seg_im_ero_label,seg_im_ero_label_temp==labels)
#seg_im_ero_label = seg_im_ero_label.astype(np.uint8)
# Step 3: perform dilation on the labelled pore space
seg_im_ero_label_dil = morphology.binary_dilation(seg_im_ero_label,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of 'seg_im_ero_label_dil' is 'bool'
seg_im_ero_label_dil[seg_image_input<=0.0]=False
Sw = 1.0 - seg_im_ero_label_dil.sum()/pore_vol
#end if
return Sw
#end def
def extract_shape(domain):
if domain.ndim == 3:
where_tube = np.where(domain)
zStart = where_tube[0].min()
zEnd = where_tube[0].max()
yStart = where_tube[1].min()
yEnd = where_tube[1].max()
xStart = where_tube[2].min()
xEnd = where_tube[2].max()
domain_seg = domain[zStart:zEnd+1,yStart:yEnd+1,xStart:xEnd+1].copy()
# IMPORTANT: if you have "domain_seg = domain[yStart:yEnd+1,xStart:xEnd+1]"
# then "domain_seg" is only a view of domain, and later on you have
# any changes on your "domain_seg", the "domain" will also be changed
# correspondingly, which might introduce unexpected conflicts and errors
else: # domain.ndim == 2
where_tube = np.where(domain)
yStart = where_tube[0].min()
yEnd = where_tube[0].max()
xStart = where_tube[1].min()
xEnd = where_tube[1].max()
domain_seg = domain[yStart:yEnd+1,xStart:xEnd+1].copy()
return domain_seg
#end def

49
example/Python/lbpm_process_1D_Snw_profiles.py Normal file
View File

@ -0,0 +1,49 @@
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from glob import glob
def read_NetCDF_file_py(file_name,data_name):
ncfile = Dataset(file_name,mode='r')
dataOUT = ncfile.variables[data_name][:]
ncfile.close()
return dataOUT
#end def
# Load the domain
domain = read_NetCDF_file_py('domain1_256_reduced_lbpm.nc','segmented')
vol_pore_1D = np.sum(np.sum(domain,axis=2),axis=1)*1.0
# Set the files & names of the phase data
data_group_prefix="OutputData_vis*.nc"
data_group = glob(data_group_prefix)
#data_group.sort()
data_name="Phase"
output_data = np.empty((0,),dtype=np.float64)
if not data_group:
print "**Error: No input file group: "+data_group_prefix
else:
for ii in np.arange(len(data_group)):
data_single = data_group[ii]
print "**Info: Start processing file "+data_single+" now..."
time_step = int(data_single[data_single.find("_vis")+len("_vis"):data_single.find(".nc")])
phase = read_NetCDF_file_py(data_group[ii],data_name)
domain_nw = (phase>0.0).astype(np.int8)
vol_nw_1D = np.sum(np.sum(domain_nw,axis=2),axis=1)*1.0
Snw_1D = vol_nw_1D/vol_pore_1D
Snw_1D[Snw_1D==0.0] = 1.0e-10
output_data = np.hstack((output_data,np.hstack((time_step,Snw_1D))))
#end for
#end if
output_data.shape = (len(data_group),output_data.size/len(data_group))
output_data = output_data[output_data[:,0].argsort(),:] #sort output data according to time_step
print "**Info: Save the 1D Snw data now."
np.savetxt('1D_Snw_data.txt',output_data)

54
example/Python/lbpm_segment_phase_field.py Normal file
View File

@ -0,0 +1,54 @@
import numpy as np
#import matplotlib.pyplot as plt
from netCDF4 import Dataset
#from glob import glob
def read_NetCDF_file_py(file_name,data_name):
ncfile = Dataset(file_name,mode='r')
dataOUT = ncfile.variables[data_name][:]
ncfile.close()
return dataOUT
#end def
def write_NetCDF_file_py(file_name,data_name,data_IN,data_IN_dtype,lx,ly,lz):
# Important !
# Note the dataOUT should have the shape (lz,ly,lx)
data_IN.shape=(lz,ly,lx) # make sure data_IN has the right shape
# open a new netCDF file for writing.
ncfile = Dataset(file_name,'w')
# create the output data.
# create the x, y and z dimensions.
ncfile.createDimension('x',lx)
ncfile.createDimension('y',ly)
ncfile.createDimension('z',lz)
# create the variable (4 byte integer in this case)
# first argument is name of variable, second is datatype, third is
# a tuple with the names of dimensions.
data = ncfile.createVariable(data_name,data_IN_dtype,('z','y','x'))
data[:] = data_IN
#data.voxel_unit = 'micrometer'
#data.voxel_size=5.7
ncfile.close()
print '**Info: The *.nc file is written successfully !'
#end def
input_nc_file_name = 'OutputData_vis200000.nc'
input_nc_file_var_name = 'Phase'
output_nc_file_name = input_nc_file_name[:-3]+'_PhaseSeg.nc'
output_nc_file_var_name = 'Phase_seg'
phase=read_NetCDF_file_py(input_nc_file_name,input_nc_file_var_name)
(lz,ly,lx) = phase.shape
phase[phase>0.0]=1
phase[phase<0.0]=0
phase = phase.astype(np.int32)
write_NetCDF_file_py(output_nc_file_name,output_nc_file_var_name,phase,phase.dtype.char,lx,ly,lz)

96
example/Python/lbpm_segmented_decomp_and_pp.py Normal file
View File

@ -0,0 +1,96 @@
#!/usr/bin/env python
import sys
import numpy as np
import skfmm
#NOTE: need python package 'array_split'
# pip install array-split
# Ref: https://pypi.python.org/pypi/array-split/0.1.3
# For more information on using array_split function, please see
# http://array-split.readthedocs.io/en/latest/reference/generated/array_split.array_split.html#array_split.array_split
from array_split import array_split
# Check if there is a proper command line argument
if len(sys.argv) !=3:
sys.stderr.write('**Error: Usage: ' + sys.argv[0] + '<Segmented.in> <Domain.in>\n')
sys.exit()
# end if
## ***************** Read 'Segmented.in' ******************** ##
print "**Info: Reading Segmented.in file."
f = open(sys.argv[1],'r')
lines = f.readlines()
# Read the file name of the segmented image
seg_data_name=lines[0]
seg_data_name=seg_data_name[:-1]
# Read the (original) size of the segmeted image
LX,LY,LZ = np.fromstring(lines[1].splitlines()[0],dtype = np.int32,sep=' ')
f.close()
## ******************************* Read 'Domain.in' ************************************* ##
print "**Info: Reading Domain.in file."
f = open(sys.argv[2],'r')
lines = f.readlines()
# Read processors configuration
nproc_x,nproc_y,nproc_z=np.fromstring(lines[0].splitlines()[0],dtype = np.int32,sep=' ')
# Read subdomain dimensions
lx,ly,lz = np.fromstring(lines[1].splitlines()[0],dtype = np.int32,sep=' ')
# Read domain dimensions
LX_temp,LY_temp,LZ_temp = np.fromstring(lines[3].splitlines()[0],dtype = np.float32,sep=' ')
LX_temp=int(LX_temp)
LY_temp=int(LY_temp)
LZ_temp=int(LZ_temp)
f.close()
# Double check if the input of 'Domain.in' is consistent with that in the 'Segmented.in'
if LX != LX_temp or LY != LY_temp or LZ != LZ_temp:
print "**Error: Inconsistent image size input between Domain.in and Segmented.in !"
print " The following domain decomposition may fail due to this inconsistency."
#end if
print "**Info: Start preparing the domain decomposition and ID field labelling..."
print "**Info: The node type convention for LBPM-WIA is as follows:"
print "*********************************"
print "** 0 -> solid phase **"
print "** 1 -> non-wetting phase **"
print "** 2 -> wetting phase **"
print "*********************************"
# Load the segmented image
# NOTE: Assume the input segmented image has ID field as follows:
# 0 -> solid nodes
# 1 or 2 -> fluid nodes
seg_data = np.fromfile(seg_data_name,dtype=np.int8)
#TODO: Make a function for calculating the signed distance, so to save the memory
seg_data_domain = seg_data.copy() # seg_data_domain is used for calculating the signed distance
seg_data_domain[seg_data_domain == 0] = -1
seg_data_domain[seg_data_domain > 0] = 1
seg_data.shape = (LZ,LY,LX)
seg_data_dist = skfmm.distance(seg_data_domain.reshape((LZ,LY,LX)))
# Decomposition
seg_data_decomp_list = array_split(seg_data,axis=[nproc_z,nproc_y,nproc_x])
seg_data_dist_decomp_list = array_split(seg_data_dist,axis=[nproc_z,nproc_y,nproc_x])
#TODO: check the size of the decomposed data, comparing to the loaded subdomain size (lz,ly,lx)
# Write the decomposed data
ID_prefix="ID."
SignDist_prefix="SignDist."
halo_layer = 1 #The length of the halo layer
for ii in np.arange(nproc_z*nproc_y*nproc_x):
ID_name=ID_prefix+"%05i" % ii
SignDist_name=SignDist_prefix+"%05i" % ii
temp = seg_data_decomp_list[ii]
temp = np.lib.pad(temp,((halo_layer,halo_layer),(halo_layer,halo_layer),(halo_layer,halo_layer)),'edge')
print "**Info: Write to file: "+ID_name
temp.tofile(ID_name)
temp = seg_data_dist_decomp_list[ii]
temp = np.lib.pad(temp,((halo_layer,halo_layer),(halo_layer,halo_layer),(halo_layer,halo_layer)),'edge')
temp.tofile(SignDist_name)
print "**Info: Write to file: "+SignDist_name
#end for

44
example/Python/lbpm_solid_coordinate_number.py Normal file
View File

@ -0,0 +1,44 @@
import numpy as np
from glob import glob
from lbpm_solid_coordinate_number_utils import *
import matplotlib.pyplot as plt
#NOTE: I could have read the 'Domain.in' files to read the information about the subdomain
# but let's assume for now all the subdomains are strictly cubic
# set the name of the full domain in *nc format
domain_file_name = "domain1_256.nc"
stats_OUT = solid_coord_fulldomain(domain_file_name)
Aws = cal_Aws_fulldomain(domain_file_name)*1.0 # Convert Aws to a float number
stats_OUT = np.vstack((stats_OUT,np.array([99,Aws]))) #The code 99 is dummy - I just want to attach the Aws to the stats_OUT data
np.savetxt(domain_file_name[:-len(".nc")]+"_stats.txt",stats_OUT)
# Trial plot
plt.figure(1)
plt.semilogy(stats_OUT[1:-1,0],stats_OUT[1:-1,1]/stats_OUT[-1,-1],'ro-')
plt.ylabel('Partial Aws / Total Aws')
plt.xlabel('Number of solid neighbours')
plt.grid(True)
plt.show()
#TODO: make the routine that can analyse individual subdomains and agglomerate the itemfreq data from all subdomains
## Load the ID field "ID.00*"
#ID_prefix="ID."
#halo_layer = 1 # the length of the halo layer
#id_group = glob(ID_prefix+'*')
#id_group.sort()
#
#
#if not id_group:
# print 'Error: No data files: '+id_group
#else:
# for ii in range(len(id_group)):
# print '**Info: Read data file: '+id_group[ii]
# print "**Info: Start analysing the solid coordinate number......"
# # call function here

112
example/Python/lbpm_solid_coordinate_number_utils.py Normal file
View File

@ -0,0 +1,112 @@
import numpy as np
from scipy import stats
import scipy.ndimage.morphology as morphology
from netCDF4 import Dataset
def read_NetCDF_file_py(file_name,data_name):
ncfile = Dataset(file_name,mode='r')
dataOUT = ncfile.variables[data_name][:]
ncfile.close()
return dataOUT
#end def
def solid_coord_fulldomain(id_field_file):
# 'id_field_file' is the file name of the full ID field in *nc format
# NOTE: This input *nc file is not the raw CT file (which has 1->solid; 0->fluid)
# Assume the fulldomain is raw, i.e. no postprocessed reservoir layers or porous plate
# Assume the node type convention for LBPM-WIA simulations, which says
# id = 0 -> solid phase
# id = 1 -> non-wetting phase
# id = 2 -> wetting phase
# -------------------------------------
print "**Info: Load the image file: "+id_field_file
domain = read_NetCDF_file_py(id_field_file,'segmented')
print "**Info: Start analysing the solid coordinate number ......"
domain = np.logical_not(domain) # Now 1 -> solid nodes; 0 -> fluid nodes
domain = domain.astype(np.int8)
# Define the D3Q19 lattice
cx=np.array([0, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1, 0, 0, 0, 0])
cy=np.array([0, 0, 0, 1, -1, 0, 0, 1, -1, -1, 1, 0, 0, 0, 0, 1, -1, 1,-1])
cz=np.array([0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1])
z_axis = 2
y_axis = 1
x_axis = 0
domain_temp = np.zeros_like(domain)
for idx in np.arange(1,19):
domain_temp += np.roll(np.roll(np.roll(domain,cx[idx],axis=x_axis),cy[idx],axis=y_axis),cz[idx],axis=z_axis)
#end for
domain_temp = domain_temp[domain==0] # only extract the coordinate number for pore space nodes
# NOTE that we have 0 -> fluid nodes and 1 -> solid nodes in domain
return stats.itemfreq(domain_temp)
#end def
def cal_Aws_fulldomain(id_field_file):
# Calculate the (total) fluid-solid interfacial area for an input porous medium 'Aws'
# 'id_field_file' is the file name of the full ID field in *nc format
# NOTE: This input *nc file is not the raw CT file (which has 1->solid; 0->fluid)
# Assume the fulldomain is raw, i.e. no postprocessed reservoir layers or porous plate
# Assume the node type convention for LBPM-WIA simulations, which says
# id = 0 -> solid phase
# id = 1 -> non-wetting phase
# id = 2 -> wetting phase
# -------------------------------------
print "**Info: Load the image file: "+id_field_file
domain = read_NetCDF_file_py(id_field_file,'segmented')
domain[domain>0]=1 # in case there are some nodes with id=2
print "**Info: Start calculating the fluid-solid interfacial area ......"
# Generate a D3Q19 structure unit cell
D3Q19=morphology.generate_binary_structure(3,2)
domain_dil=morphology.binary_dilation(domain,structure=D3Q19,border_value=0).astype(np.int8)
#NOTE: It is important to have 'border_value' set as 0 for morphological dilation !
domain_dil = domain_dil - domain
if (domain_dil<0).sum() > 0:
print "**Error: The domain for calculating the fluid-solid interfacial area is not set properly !"
#end if
Aws=domain_dil.sum()
return Aws
#end def
#TODO: implement a similar function to analyse only the subdomain
# and find a way to concatenate the itemfeq data from all subdomains
def solid_coord_subdomain(id_field_file):
#NOTE: This function is not done yet *******************************************
# 'id_field_file' is the file name of the ID field, i.e. "ID.00*"
# Assume the input doamin is a subdomain segmented by lbpm_segment_decomp.cpp
# thus there is an extra halo layer wrapping the subdomain
# Assume the subdomain is strictly cubic
# Assume the subdomain is raw, i.e. originally segmented CT image with no reservoir layers or porous plate
# Assume the node type convention for LBPM-WIA simulations, which says
# id = 0 -> solid phase
# id = 1 -> non-wetting phase
# id = 3 -> wetting phase
# -------------------------------------
domain = np.fromfile(id_field_file,dtype=np.int8)
l_domain = int(pow(domain.size,1.0/3.0))
domain.shape=(l_domain,l_domain,l_domain)
domain = np.logical_not(domain) # Now 1 -> solid nodes; 0 -> fluid nodes
domain = domain.astype(np.int8)
# Define the D3Q19 lattice
cx=np.array([0, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1, 0, 0, 0, 0])
cy=np.array([0, 0, 0, 1, -1, 0, 0, 1, -1, -1, 1, 0, 0, 0, 0, 1, -1, 1,-1])
cz=np.array([0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1])
z_axis = 2
y_axis = 1
x_axis = 0
domain_temp = np.zeros_like(domain)
for idx in np.arange(1,19):
domain_temp += np.roll(np.roll(np.roll(domain,cx[idx],axis=x_axis),cy[idx],axis=y_axis),cz[idx],axis=z_axis)
#end for
#end def

95
example/Python/lbpm_timelog_analysis_v1.py Normal file
View File

@ -0,0 +1,95 @@
#!/usr/bin/env python
import sys
import numpy as np
import matplotlib.pyplot as plt
# Check if there is a proper command line argument
if len(sys.argv) !=2:
sys.stderr.write('**Error: Usage: ' + sys.argv[0] + ' <Color.in>\n')
sys.exit()
# end if
## *********************** Read 'Color.in' ***************************** ##
print "**Info: Reading Color.in file."
f = open(sys.argv[1],'r')
lines = f.readlines()
tau=float(lines[0])
alpha, beta, phi_s = np.fromstring(lines[1].splitlines()[0],dtype=np.float64, sep=' ')
# Read the 3-D body forces
Fx, Fy, Fz = np.fromstring(lines[3].splitlines()[0],dtype=np.float64, sep=' ')
# There can be more to be read from Color.in file
f.close()
## ******************** END: Read 'Color.in' *************************** ##
# Load the data file
print "**Info: Reading timelog.tcat file."
data_all = np.genfromtxt('timelog.tcat',names=True)
# Load individual data sets
pw_all=data_all['pw']
#pw=pw_all[np.isfinite(pw_all)]
sw=data_all['sw']
sw=sw[np.isfinite(pw_all)]
vanz=data_all['vanz']
vanz=vanz[np.isfinite(pw_all)]
time_step=data_all['time']
time_step=time_step[np.isfinite(pw_all)]
Jwn = data_all['Jwn']
Jwn=Jwn[np.isfinite(pw_all)]
## ********************************** LBPM-WIA related parameters *********************************** ##
# (lattice) interfacial tension
#TODO: need to load Color.in file
IFT_conv = 5.796
IFT = alpha*IFT_conv
viscosity = (tau-0.5)/3.0 # lattice dynamic viscosity
cos_theta=phi_s #NOTE: this is only an approximation
# Plot Capillary number
#Ca=vanz*viscosity/IFT/cos_theta
Ca=vanz*viscosity/IFT
plt.figure(1)
plt.subplot(1,2,1)
plt.plot(time_step,Ca,'ks--',label='Total Ca')
plt.plot(np.array([0,time_step.max()]),np.array([0,0]),'r-',linewidth=2)
#plt.plot(time_step,vanz,'ro--')
plt.grid(True)
plt.legend(loc='best')
plt.xlabel('Time step (LBM)')
plt.ylabel('Capillary number (w.r.t non-wetting phase)')
plt.subplot(1,2,2)
plt.semilogy(time_step[Ca>0],Ca[Ca>0],'b>--',label='Positive Ca')
plt.semilogy(time_step[Ca<0],np.abs(Ca[Ca<0]),'ro--',label='Negative Ca')
#plt.plot(time_step,vanz,'ro--')
plt.grid(True)
plt.legend(loc='best')
plt.xlabel('Time step (LBM)')
plt.ylabel('Capillary number (w.r.t non-wetting phase)')
plt.figure(2)
plt.plot(time_step,sw,'ro--')
plt.grid(True)
plt.xlabel('Time step (LBM)')
plt.ylabel('Wetting phase saturation Sw')
plt.figure(3)
plt.plot(time_step,-1*Jwn,'b>--')
plt.grid(True)
plt.xlabel('Time step (LBM)')
plt.ylabel('Interfacial mean curvature (-1.0*Jwn)')
plt.show()