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

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
James E McClure committed 2016-11-01 08:33:31 -04:00
2 parents 2c0c00c114 1be9311d8a
commit 2b470d0f43
14 files changed
+1569 -99
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

IO/Writer.cpp
+2 -1
View File
@@ -71,7 +71,8 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
ERROR("Unknown mesh");
}
fclose(fid);
std::unique(mesh_entry.variables.begin(),mesh_entry.variables.end());
std::sort( mesh_entry.variables.begin(), mesh_entry.variables.end() );
mesh_entry.variables.erase( std::unique( mesh_entry.variables.begin(), mesh_entry.variables.end() ), mesh_entry.variables.end() );
meshes_written.push_back(mesh_entry);
}
return meshes_written;
analysis/analysis.cpp
+2 -3
View File
@@ -339,9 +339,8 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_
neighbors.push_back( rank_info.rank[1][2][1] );
neighbors.push_back( rank_info.rank[1][1][0] );
neighbors.push_back( rank_info.rank[1][1][2] );
std::unique(neighbors.begin(),neighbors.end());
if ( std::find(neighbors.begin(),neighbors.end(),rank) != neighbors.end() )
neighbors.erase(std::find(neighbors.begin(),neighbors.end(),rank));
std::sort( neighbors.begin(), neighbors.end() );
neighbors.erase( std::unique( neighbors.begin(), neighbors.end() ), neighbors.end() );
// Create a map of all local ids to the neighbor ids
std::map<int64_t,global_id_info_struct> map;
std::set<int64_t> local;
tests/CMakeLists.txt
+1
View File
@@ -1,6 +1,7 @@
# Copy files for the tests
ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
ADD_LBPM_EXECUTABLE( lbpm_color_simulator )
ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator )
ADD_LBPM_EXECUTABLE( lbpm_sphere_pp )
ADD_LBPM_EXECUTABLE( lbpm_random_pp )
ADD_LBPM_EXECUTABLE( lbpm_morphdrain_pp )
tests/lbpm_color_macro_simulator.cpp
+1003
View File
File diff suppressed because it is too large. Load diff
tests/lbpm_color_simulator.h
+23 -15
View File
@@ -61,18 +61,19 @@ public:
std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ ):
timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_) { }
phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_)
{
MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
}
~BlobIdentificationWorkItem1() { MPI_Comm_free(&newcomm); }
virtual void run() {
ThreadPool::WorkItem::d_state = 1; // Change state to in progress
// Compute the global blob id and compare to the previous version
PROFILE_START("Identify blobs",1);
MPI_Comm newcomm;
MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
double vF = 0.0;
double vS = -1.0; // one voxel buffer region around solid
IntArray& ids = new_index->second;
new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,newcomm);
MPI_Comm_free(&newcomm);
PROFILE_STOP("Identify blobs",1);
ThreadPool::WorkItem::d_state = 2; // Change state to finished
}
@@ -85,6 +86,7 @@ private:
const DoubleArray& dist;
BlobIDstruct last_id, new_index, new_id;
BlobIDList new_list;
MPI_Comm newcomm;
};
class BlobIdentificationWorkItem2: public ThreadPool::WorkItem
{
@@ -93,13 +95,15 @@ public:
std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ ):
timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_) { }
phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_)
{
MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
}
~BlobIdentificationWorkItem2() { MPI_Comm_free(&newcomm); }
virtual void run() {
ThreadPool::WorkItem::d_state = 1; // Change state to in progress
// Compute the global blob id and compare to the previous version
PROFILE_START("Identify blobs maps",1);
MPI_Comm newcomm;
MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
const IntArray& ids = new_index->second;
static int max_id = -1;
new_id->first = new_index->first;
@@ -118,7 +122,6 @@ public:
getNewIDs(map,max_id,*new_list);
writeIDMap(map,timestep,id_map_filename);
}
MPI_Comm_free(&newcomm);
PROFILE_STOP("Identify blobs maps",1);
ThreadPool::WorkItem::d_state = 2; // Change state to finished
}
@@ -131,6 +134,7 @@ private:
const DoubleArray& dist;
BlobIDstruct last_id, new_index, new_id;
BlobIDList new_list;
MPI_Comm newcomm;
};
@@ -140,7 +144,11 @@ class WriteVisWorkItem: public ThreadPool::WorkItem
public:
WriteVisWorkItem( int timestep_, std::vector<IO::MeshDataStruct>& visData_,
TwoPhase& Avgerages_, fillHalo<double>& fillData_ ):
timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_) {}
timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_)
{
MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
}
~WriteVisWorkItem() { MPI_Comm_free(&newcomm); }
virtual void run() {
ThreadPool::WorkItem::d_state = 1; // Change state to in progress
PROFILE_START("Save Vis",1);
@@ -156,10 +164,7 @@ public:
fillData.copy(Averages.Press,PressData);
fillData.copy(Averages.SDs,SignData);
fillData.copy(Averages.Label_NWP,BlobData);
MPI_Comm newcomm;
MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
IO::writeData( timestep, visData, 2, newcomm );
MPI_Comm_free(&newcomm);
PROFILE_STOP("Save Vis",1);
ThreadPool::WorkItem::d_state = 2; // Change state to finished
};
@@ -169,8 +174,10 @@ private:
std::vector<IO::MeshDataStruct>& visData;
TwoPhase& Averages;
fillHalo<double>& fillData;
MPI_Comm newcomm;
};
// Helper class to run the analysis from within a thread
// Note: Averages will be modified after the constructor is called
class AnalysisWorkItem: public ThreadPool::WorkItem
@@ -180,6 +187,7 @@ public:
BlobIDstruct ids, BlobIDList id_list_, double beta_ ):
type(type_), timestep(timestep_), Averages(Averages_),
blob_ids(ids), id_list(id_list_), beta(beta_) { }
~AnalysisWorkItem() { }
virtual void run() {
ThreadPool::WorkItem::d_state = 1; // Change state to in progress
Averages.NumberComponents_NWP = blob_ids->first;
@@ -243,7 +251,7 @@ void run_analysis( int timestep, int restart_interval,
// Identify blobs and update global ids in time
type = static_cast<AnalysisType>( type | IdentifyBlobs );
}
#ifdef USE_CUDA
/* #ifdef USE_CUDA
if ( tpool.getQueueSize()<=3 && tpool.getNumThreads()>0 && timestep%50==0 ) {
// Keep a few blob identifications queued up to keep the processors busy,
// allowing us to track the blobs as fast as possible
@@ -251,7 +259,7 @@ void run_analysis( int timestep, int restart_interval,
type = static_cast<AnalysisType>( type | IdentifyBlobs );
}
#endif
*/
if ( timestep%ANALYSIS_INTERVAL == 0 ) {
// Copy the averages to the CPU (and identify blobs)
type = static_cast<AnalysisType>( type | CopySimState );
@@ -369,7 +377,7 @@ void run_analysis( int timestep, int restart_interval,
// Retain the timestep associated with the restart files
if (rank==0){
FILE *Rst = fopen("Restart.txt","w");
fprintf(Rst,"%i\n",timestep);
fprintf(Rst,"%i\n",timestep+5);
fclose(Rst);
}
// Write the restart file (using a seperate thread)
tests/lbpm_permeability_simulator.cpp
+3 -4
View File
@@ -228,8 +228,8 @@ int main(int argc, char **argv)
printf("Force(x) = %f \n", Fx);
printf("Force(y) = %f \n", Fy);
printf("Force(z) = %f \n", Fz);
printf("Sub-domain size = %i x %i x %i\n",Nz,Nz,Nz);
printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
printf("Process grid = %i x %i x %i\n",nprocx,nprocy,nprocz);
printf("********************************************************\n");
}
@@ -246,8 +246,7 @@ int main(int argc, char **argv)
MPI_Barrier(comm);
Nz += 2;
//Nx = Ny = Nz; // Cubic domain
Nx += 2; Ny += 2; Nz += 2;
int N = Nx*Ny*Nz;
int dist_mem_size = N*sizeof(double);
tests/lbpm_segmented_decomp.cpp
+32 -16
View File
@@ -18,17 +18,22 @@ int main(int argc, char **argv)
// Initialize MPI
int rank, nprocs;
MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
int SOLID=atoi(argv[1]);
{
int SOLID=atoi(argv[1]);
int NWP=atoi(argv[2]);
//char NWP,SOLID;
//SOLID=argv[1][0];
//NWP=argv[2][0];
if (rank==0){
printf("Solid Label %i \n",SOLID);
printf("NWP Label %i \n",NWP);
printf("Solid Label: %i \n",SOLID);
printf("NWP Label: %i \n",NWP);
}
//.......................................................................
// Reading the domain information file
//.......................................................................
@@ -189,22 +194,15 @@ int main(int argc, char **argv)
nx+=2; ny+=2; nz+=2;
N=nx*ny*nz;
//if (rank==0) printf("WARNING: assumed that INPUT: WP(1),NWP(2) OUTPUT will be NWP(1),WP(2) (lbpm_segmented_decomp) \n");
if (rank==0) printf("All sub-domains recieved \n");
// Really need a better way to do this -- this is to flip convention for a particular data set
for (k=0;k<nz;k++){
for (j=0;j<ny;j++){
for (i=0;i<nx;i++){
n = k*nx*ny+j*nx+i;
if (Dm.id[n]==char(SOLID)) Dm.id[n] = 0;
else if (Dm.id[n]==char(NWP)) Dm.id[n] = 1;
else Dm.id[n] = 2;
else if (Dm.id[n]==char(NWP)) Dm.id[n] = 1;
else Dm.id[n] = 2;
// if (Dm.id[n] == 1) Dm.id[n]=2;
//else if (Dm.id[n] == 2) Dm.id[n]=1;
// Initialize the solid phase
// if (Dm.id[n] == 0) id[n] = 0;
// else id[n] = 1;
}
}
}
@@ -230,6 +228,24 @@ int main(int argc, char **argv)
float porosity = float(totalGlobal-countGlobal)/totalGlobal;
if (rank==0) printf("Porosity=%f\n",porosity);
count = 0;
total = 0;
countGlobal = 0;
totalGlobal = 0;
for (k=1;k<nz-1;k++){
for (j=1;j<ny-1;j++){
for (i=1;i<nx-1;i++){
n=k*nx*ny+j*nx+i;
if (Dm.id[n] != 0) total++;
if (Dm.id[n] == 2) count++;
}
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
MPI_Allreduce(&total,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
float saturation = float(countGlobal)/totalGlobal;
if (rank==0) printf("wetting phase saturation=%f\n",saturation);
char LocalRankFilename[40];
@@ -264,7 +280,7 @@ int main(int argc, char **argv)
// fwrite(SymDist.get(),8,SymDist.length(),SYMDIST);
fwrite(symid,1,N,SYMID);
fclose(SYMID);
}
MPI_Barrier(comm);
MPI_Finalize();
}
tests/lbpm_segmented_pp.cpp
+2 -1
View File
@@ -185,7 +185,7 @@ int main(int argc, char **argv)
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
{
//.......................................................................
// Reading the domain information file
//.......................................................................
@@ -443,6 +443,7 @@ int main(int argc, char **argv)
Averages.SortBlobs();
Averages.PrintComponents(timestep);
*/
}
MPI_Barrier(comm);
MPI_Finalize();
return 0;
workflows/Preprocess/preprocess.py
+123
View File
@@ -0,0 +1,123 @@
#!/usr/bin/env python
import numpy as np
import sys
import os
import shutil
from preprocess_utils import *
# 1. this script take one input argument: experiment.csv
# 2. Users need to put two things into this script
# 2.1 surface tension in physical units
# 2.2 physical depth
# 3. Read csv + convert to LBM information
# 4. write Color.in,-> this needs info from LBM pressure
# 5. write Segmented.in -> this needs names of segmented image data
# 6. write Domain.in -> this needs info on how users want to decompose the domain
# Check if there is a proper command line argument
if len(sys.argv) !=2:
sys.stderr.write('Usage: ' + sys.argv[0] + ' <Input experiment file>\n')
sys.exit()
# end if
#experiment_file = 'experiment.csv' # Give the name of the experiment file (e.g. *.csv)
experiment_file = sys.argv[1] # Give the name of the experiment file (e.g. *.csv)
seg_image_data_suffix = '_segmented.raw'# suffix of the segmented image data
# (should be consistent with what's in the segmenting script)
# TODO: It'd be good to find a better way to do this
image_format = '.tiff'
process_name = 'drainage'
#### Users need to put information here ####
ift = 24.0 # dyne/cm
Depth = 8.8 # micron
# A list of all default values in 'Color.in'
# Para['tau'] = 0.7
# Para['alpha'] = 0.005
# Para['beta'] = 0.95
# Para['phi_solid'] = -1.0
# Para['saturation'] = 0.0
# Para['Fx'] = 0.0
# Para['Fy'] = 0.0
# Para['Fz'] = 0.0
# Para['Restart'] = 0
# Para['pBC'] = 1
# Para['din'] = 1.001
# Para['dout'] = 0.999
# Para['maxtime'] = 100005
# Para['interval'] = 2000
# Para['tolerance'] = 1e-5
# ***** Update any variables in 'Color.in', using names given in the key ***** #
alpha = 0.01
# **************************************************************************** #
# A list of all default values in 'Domain.in'
# Para['nprocx'] = 1
# Para['nprocy'] = 2
# Para['nprocz'] = 2
# Para['nx'] = 1
# Para['ny'] = 2
# Para['nz'] = 2
# Para['nspheres'] = 0 # deprecated
# Para['Lx'] = 10
# Para['Ly'] = 500
# Para['Lz'] = 500
# ***** Update any variables in 'Domain.in', using names given in the key ***** #
# ***************************************************************************** #
# A list of all default values in 'Segmented.in'
# Para['file_name'] = 'Micromodel_1_segmented.raw'
# Para['Nx'] = 10
# Para['Ny'] = 500
# Para['Nz'] = 500
# Para['xStart'] = 0
# Para['yStart'] = 0
# Para['zStart'] = 0
# ***** Update any variables in 'Segmented.in', using names given in the key ***** #
# ******************************************************************************** #
# Extract key parameters for LBM simulation from the experimental input *.csv file
(Seg_data_name,din,dout)=get_LBM_parameters(experiment_file,seg_image_data_suffix,image_format,\
ift=ift,Depth=Depth)
# Now 'name_for_Segmented_in' should match the name of segmented data files that are already generated
# Write out 'Color.in', 'Domain.in' and 'Segmented.in' files
cwd = os.getcwd()
for k in range(Seg_data_name.size):
tag = k+1 # tag number for different folders to be created
dir_name = process_name+'_'+str(tag)
print "Creating folder : "+dir_name
if not os.path.exists(dir_name):
os.mkdir(dir_name)
#end if
# Either move the corresponding '*.raw' data file to the folder 'dir_name'
#os.rename('./'+Seg_data_name[k],'./'+dir_name+'/'+Seg_data_name[k])
# Or copy the corresponding '*.raw' data file to the folder 'dir_name'
shutil.copy('./'+Seg_data_name[k],'./'+dir_name)
# Change to the corresponding folder and write all input files
os.chdir(dir_name)
write_Color_in_file(din=din[k],dout=dout[k],alpha=alpha)
write_Segment_in_file(Seg_data_name=Seg_data_name[k])
write_Domain_in_file()
os.chdir(cwd)
#end for
workflows/Preprocess/preprocess_utils.py
+272
View File
@@ -0,0 +1,272 @@
#!/usr/bin/env python
import os
import numpy as np
import csv
def get_LBM_parameters(csv_file_name,base_name_suffix,image_format,ift=24,Depth=8.8,**kwargs):
# 'ift': dyne/cm
# 'Depth': micron
# Users need to provide the following information
# 1. surface tension in physical unit
# 2. physical depth in physical unit
# 3. experimental file e.g. *.csv
# 4. Other optional information, which is summarised in a dictionary
# 'Para', including:
# Para['D']: characteristic length
# Para['alpha']: LBM parameters, controlling the surface tension
# Para['fitting']: LBM parameters, extracted from the bubble test
# Experimental definitions - units are converted in terms of N, m
micron=1e-6
cm=1e-2
dyne=1e-5
kPa=1e3
Para={'D':30.0} # characteristic length
#length=500*micron
# LBM related parameters
# TODO: for parameters like 'alpha', you should find a way to
# communicate with other functions
# It cannot be hard-coded here
Para['alpha'] = 0.01
Para['fitting'] = 5.796
# Check users' input arguments
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: LBM pressure boundaries are not set successfully !"
return
#end if
#end for
IFT = Para['alpha'] * Para['fitting']
# 'ift' is the surface tension from the experiment in physical unit
ift=ift*dyne/cm # This converts the unit of 'ift' to SI unit
# Process the experimental data
# 'EXP_data' : original experimental data
# It is a recorded numpy array, which means that its component
# can be accessed by 'EXP_data.sw' or 'EXP_data['sw']'
# About 'EXP_data' see more information from the function 'read_csv'
EXP_data = read_csv(csv_file_name)
# A few notes for the experimental data
# 'Jwn': mean curvature in physical unit (e.g. 1/[micron])
# 'pwn': pressure difference in physical (e.g. kPa)
# Overall we need to map the measured physical pressures to get the principle radius of curvature R1
# and associated curvature J1, and similarly R2 and J2 in the model's depth direction
# J1: principal curvature 1
# J2: principal curvature 2 along the model's depth direction
# 'pc' is the dimensionless mean curvature scaled by the length scale 'D32'
# 'pc' is extracted from experimentally measured mean curvature 'Jwn'
pc = Para['D']*micron*EXP_data.Jwn*(1.0/micron)
# Alternatively, the dimensionless mean curvature can also be extracted from the
# experimentally measured pressure difference (i.e. capillary pressure)
# 'pwn' is the dimensionless mean curvature scaled by the length scale 'D32'
# pwn = Para['D']*micron*EXP_data.pwn*kPa/ift
# Curvature is fixed in the micromodel "depth" direction
J2 = Para['D']*micron/(Depth*micron/2.0)
# infer the curvature in the other direction
J1 = pc-J2
# determine the LBM pressure difference
dp=(J1+J2)*IFT/Para['D']
# determine the boundary pressure values for Color.in
din = 1.0+0.5*dp
dout = 1.0-0.5*dp
# Generate names of segmented image data for Segmented.in file
# Need the input 'base_name_suffix' (e.g. '_segmented.raw')
data_name_for_Segmented_in = EXP_data.image_name.copy()
data_name_for_Segmented_in = data_name_for_Segmented_in.replace(image_format,base_name_suffix)
return (data_name_for_Segmented_in,din,dout)
#end def
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']=2
Para['nprocz']=2
Para['nx']=1
Para['ny']=2
Para['nz']=2
Para['nspheres']=0 # deprecated
Para['Lx']=10
Para['Ly']=500
Para['Lz']=500
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 "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']=10
Para['Ny']=500
Para['Nz']=500
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 "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.01
Para['beta']=0.95
Para['phi_solid']=-1.0
Para['saturation']=0.0
Para['Fx']=0.0
Para['Fy']=0.0
Para['Fz']=0.0
Para['Restart']=0
Para['pBC']=1
Para['din']=1.001
Para['dout']=0.999
Para['maxtime']=100005
Para['interval']=2000
Para['tolerance']=1e-5
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("%f\n" % Para['tau'])
ParamFile.write("%f " % Para['alpha'])
ParamFile.write("%f " % Para['beta'])
ParamFile.write("%f\n" % Para['phi_solid'])
ParamFile.write("%f\n" % Para['saturation'])
ParamFile.write("%f " % Para['Fx'])
ParamFile.write("%f " % Para['Fy'])
ParamFile.write("%f\n" % Para['Fz'])
ParamFile.write("%i " % Para['Restart'])
ParamFile.write("%i " % Para['pBC'])
ParamFile.write("%f " % Para['din'])
ParamFile.write("%f\n" % Para['dout'])
ParamFile.write("%i " % Para['maxtime'])
ParamFile.write("%i " % Para['interval'])
ParamFile.write("%e\n" % Para['tolerance'])
ParamFile.close()
print "A Color.in file is created."
#end def
def read_csv(csv_file_name):
#TODO Haven't thought about a better way of doing this
# Right now I just hard-code the possible variables from the experiment
# which means users are forced to prepare their *.csv file as listed here
image_name = []
sw = []
Jwn = []
pwn = []
with open(csv_file_name,"r") as f:
for line in f:
reader=csv.reader(f,delimiter=' ') #Note the header is skipped by default
for row in reader:
image_name.append(row[0])
sw.append(float(row[1]))
Jwn.append(float(row[2]))
pwn.append(float(row[3]))
#end for
#end for
#end with
# Convter the list to numpy array
image_name_array = np.asarray(image_name,dtype=str)
sw_array = np.asarray(sw,dtype=np.float32)
Jwn_array = np.asarray(Jwn,dtype=np.float32)
pwn_array = np.asarray(pwn,dtype=np.float32)
# Restore the original shape of the experimental data
experiment_data = np.column_stack((image_name_array,sw_array,Jwn_array,pwn_array))
# Unfortunately the column stack will cast all float data to strings
experiment_data = np.core.records.fromrecords(experiment_data,names='image_name,sw,Jwn,pwn')
dt=experiment_data.dtype.descr
dt[1] = (dt[1][0],'float32')
dt[2] = (dt[2][0],'float32')
dt[3] = (dt[3][0],'float32')
experiment_data = experiment_data.astype(dt)
return experiment_data
#end def
workflows/SegmentingImage/Check_SignDist_and_ID.py
+35
View File
@@ -0,0 +1,35 @@
#!/usr/bin/env python
import numpy as np
import matplotlib.pylab as plt
lx = 3
ly = 500
lz = 500
# Load the data
SignDist = np.fromfile('SignDist.00000',dtype = np.float64)
SignDist.shape = (lz+2,ly+2,lx+2)
ID = np.fromfile('ID.00000',dtype = np.uint8)
ID.shape = (lz+2,ly+2,lx+2)
# Plot
plt.figure(1)
plt.title('SignDist Map')
plt.pcolormesh(SignDist[:,:,2])
plt.grid(True)
plt.axis('equal')
plt.colorbar()
plt.figure(2)
plt.title('ID.xxxxx')
plt.pcolormesh(ID[:,:,2],cmap='hot')
plt.grid(True)
plt.axis('equal')
plt.show()
workflows/SegmentingImage/SegmentMicromodelTiff_Group.py
+27 -38
View File
@@ -5,36 +5,32 @@ import matplotlib.pylab as plt
from glob import glob
from PIL import Image # import Python Imaging Library (PIL)
import sys
from SegmentMicromodelTiff_utils import convert_image_to_array,read_input_parameters
from SegmentMicromodelTiff_utils import *
# Check if there is a proper command line argument
if len(sys.argv) !=2:
sys.stderr.write('Usage: ' + sys.argv[0] + ' <Input parameter file>\n')
print 'Input parameter file: '
print 'line 0: the base name of the experimental images\n',\
'line 1: imin imax\n',\
'line 2: imin_b imax_b\n',\
'line 3: top_layer bottom_layer DEPTH\n',\
'line 4: threshold_nw threshold_s'
sys.exit()
# end if
input_parameter_file = sys.argv[1] # Give the name of the input parameter file
base_name_output = '_segmented.raw'# The string will be appended at the end of the
# name of the input image as an output file name
BaseName_output_init_config = '_InitConfig_segmented.raw'# The string will be appended at the end of the
# name of the input image as an output file name
BaseName_output_domain = '_segmented.raw'# The string will be appended at the end of the
# name of the input image as an output file name
image_format = '.tiff' # format of experimental images
# Read other input parameters
# Read other input parameters to get NX, NY and NZ
Para = read_input_parameters(input_parameter_file)
imin = Para['imin']
imax = Para['imax']
imin_b = Para['imin_b'] # greater than 'imin'
imax_b = Para['imax_b'] # smaller than 'imax'
top_layer = Para['top_layer']
bot_layer = Para['bot_layer']
DEPTH = Para['DEPTH']
NY = imax - imin
NY = Para['imax'] - Para['imin']
NZ = NY
NX = DEPTH+top_layer+bot_layer
# parameters for segmentation
threshold_nw = Para['threshold_nw'] # NW phase: RGB values > threshold_nw
threshold_s = Para['threshold_s'] # solid: RGB values < threshold_s
# W phase: threshold_s <= RGB values <= threshold_nw
NX = Para['DEPTH']+Para['top_layer']+Para['bot_layer']
# Load a group of image files with 'base_name' in the names of files
# e.g. 'Micromodel_1.tiff' etc.
@@ -50,31 +46,24 @@ else:
print "------ Micromodel dimensions (NX, NY, NZ) are (%i, %i, %i) ------" % (NX,NY,NZ)
# Get an array from the input .tiff image
im_array = convert_image_to_array(file_input_single,imin,imax)
# Initialize the segmented image 'ID'
# NOTE: 'ID' has the dimension (height, width, DEPTH)
ID = np.zeros((NZ,NY,NX),dtype=np.uint8)
# Map the picels to the 3D geometry
# 1. Identify the non-wetting phase
ID[im_array>threshold_nw,top_layer:top_layer+DEPTH] = 1
# 2. Identify the wetting phase
ID[np.logical_and(im_array>=threshold_s,im_array<=threshold_nw),\
top_layer:top_layer+DEPTH] = 2
# 3. Post boundary retreatment along y-axis. Note z-axis is always the flow direction
ID[:,:imin_b-imin,top_layer:top_layer+DEPTH]=0
ID[:,ID.shape[1]-(imax-imax_b):,top_layer:top_layer+DEPTH]=0
im_array = convert_image_to_array(file_input_single,Para['imin'],Para['imax'])
# Get initial fluid configuration shown in the *.Tiff image
ID = get_initial_config(im_array,Para,NX,NY,NZ)
# Get only the domain (0: solid, 1: fluid)
PorousMedium = get_porous_medium(im_array,Para,NX,NY,NZ)
# Calculate the porosity
POROSITY = 1.0 - 1.0*ID[ID==0].size/ID.size
print "Porosity of the micromodel is: "+str(POROSITY)
# Write the segmeted data to the output
# set the output file names (need to chop off say the '.tiff' part)
output_file_name = file_input_single[:file_input_single.find('.')]+base_name_output
output_file_name_domain = file_input_single[:file_input_single.find('.')]+BaseName_output_domain
output_file_name_init_config = file_input_single[:file_input_single.find('.')]+BaseName_output_init_config
# Write out the segmented data (in binary format)
print "The segmented image is written to "+output_file_name+" now..."
ID.tofile(output_file_name)
print "The segmented fluid configuration is written to "+output_file_name_init_config+" now..."
print "The segmented porous domain is written to "+output_file_name_domain+" now..."
ID.tofile(output_file_name_init_config)
PorousMedium.tofile(output_file_name_domain)
# NOTE when you want to load the binary file
# Do as follows:
# ID_reload = np.fromfile(file_name,dtype=np.uint8)
@@ -83,7 +72,7 @@ else:
#end if
# In case you want to test the segmentation
### In case you want to test the segmentation
#plt.figure(1)
#plt.subplot(1,2,1)
#plt.title('original RGB figure')
@@ -92,7 +81,7 @@ else:
#plt.colorbar()
#plt.subplot(1,2,2)
#plt.title('Segmented image')
# This will show the last-processed segmented image
### This will show the last-processed segmented image
#cmap = plt.cm.get_cmap('hot',3) #Plot 3 discrete colors for NW, W and Solids
#cax=plt.pcolormesh(ID[:,:,NX/2],cmap=cmap,vmin=-0.5,vmax=2.5);
#cbar = plt.colorbar(cax,ticks=[0,1,2])
workflows/SegmentingImage/SegmentMicromodelTiff_utils.py
+41 -18
View File
@@ -48,21 +48,44 @@ def read_input_parameters(input_file_name):
return output_file
#end def
#def get_segmented_array(image_array,parameters,NX,NY,NZ):
# # 'image_array' is a numpy array of RGB values, with the same size as micromodel
# # 'parameters' is a dictionary of input parameters
#
# # Initialize the segmented image 'ID'
# # NOTE: 'ID' has the dimension (DEPTH, height, width)
# ID = np.zeros((NX,NZ,NY),dtype=np.uint8)
# # Map the picels to the 3D geometry
# # 1. Identify the non-wetting phase
# ID[top_layer:top_layer+DEPTH,im_array>threshold_nw] = 1
# # 2. Identify the wetting phase
# ID[top_layer:top_layer+DEPTH,\
# np.logical_and(im_array>=threshold_s,im_array<=threshold_nw)] = 2
# # 3. Post boundary retreatment along y-axis
# ID[top_layer:top_layer+DEPTH,:,:imin_b-imin]=0
# ID[top_layer:top_layer+DEPTH,:,ID.shape[2]-(imax-imax_b):]=0
# return ID
##end def
def get_initial_config(image_array,Para,NX,NY,NZ):
# This function will give the same fluid configuration shown in the *.Tiff image
# Function input:
# 'image_array' is a numpy array of RGB values, with the same size as micromodel
# 'Para' is a dictionary of input parameters
# Initialize the segmented image 'ID'
# NOTE: 'ID' has the dimension (height, width, DEPTH) or (NZ,NY,NX)
ID = np.zeros((NZ,NY,NX),dtype=np.uint8)
# Map the pixels to the 3D geometry
# Rules: NW phase: RGB values > threshold_nw
# Solid: RGB values < threshold_s
# 1. Identify the non-wetting phase
ID[image_array>Para['threshold_nw'],Para['top_layer']:Para['top_layer']+Para['DEPTH']] = 1
# 2. Identify the wetting phase
ID[np.logical_and(image_array>=Para['threshold_s'],image_array<=Para['threshold_nw']),\
Para['top_layer']:Para['top_layer']+Para['DEPTH']] = 2
# 3. Post boundary retreatment along y-axis. Note z-axis is always the flow direction
ID[:,:Para['imin_b']-Para['imin'],Para['top_layer']:Para['top_layer']+Para['DEPTH']]=0
ID[:,ID.shape[1]-(Para['imax']-Para['imax_b']):,Para['top_layer']:Para['top_layer']+Para['DEPTH']]=0
return ID
#end def
def get_porous_medium(image_array,Para,NX,NY,NZ):
# This function only gives the domain (0->solid,1->fluid) of the *.Tiff image
# Function input:
# 'image_array' is a numpy array of RGB values, with the same size as micromodel
# 'Para' is a dictionary of input parameters
domain = np.zeros((NZ,NY,NX),dtype=np.uint8)
# Map the pixels to the 3D geometry
# Rules: NW phase: RGB values > threshold_nw
# Solid: RGB values < threshold_s
# 1. Identify fluid phase
domain[image_array>=Para['threshold_s'],Para['top_layer']:Para['top_layer']+Para['DEPTH']]=1
# 2. Post boundary retreatment along y-axis. Note z-axis is always the flow direction
domain[:,:Para['imin_b']-Para['imin'],Para['top_layer']:Para['top_layer']+Para['DEPTH']]=0
domain[:,domain.shape[1]-(Para['imax']-Para['imax_b']):,Para['top_layer']:Para['top_layer']+Para['DEPTH']]=0
return domain
#end def
workflows/imbibition/setup-imbibition.py
+3 -3
View File
@@ -4,7 +4,7 @@ import csv
import glob
name=sys.argv[1]
numpts=int(sys.argv[2])
#numpts=int(sys.argv[2])
process="drain"
cwd=os.getcwd()
@@ -13,7 +13,7 @@ print("Drainage should initialize imbibition!")
print("Base name is "+name)
# Parameters for color LBM
tau=0.7
tau=0.8
alpha=0.005
beta=0.95
phisol=-1.0
@@ -25,7 +25,7 @@ Restart=1
pBC=3
din=1.001
dout=0.999
maxtime=100005
maxtime=2000005
interval=50000
tolerance=1e-5