This commit is contained in:
James E McClure 2016-10-28 15:20:29 -04:00
commit 0c0c39dd81
13 changed files with 1645 additions and 40 deletions

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;

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;

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_simulator_basic )
ADD_LBPM_EXECUTABLE( lbpm_sphere_pp )
ADD_LBPM_EXECUTABLE( lbpm_random_pp )
ADD_LBPM_EXECUTABLE( lbpm_morphdrain_pp )

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;

View File

@ -0,0 +1,994 @@
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/TwoPhase.h"
#include "common/MPI_Helpers.h"
#include "ProfilerApp.h"
//#include "threadpool/thread_pool.h"
//#include "lbpm_color_simulator.h"
/*
* Simulator for two-phase flow in porous media
* James E. McClure 2013-2016
*/
using namespace std;
//*************************************************************************
// Implementation of Two-Phase Immiscible LBM using CUDA
//*************************************************************************
inline void PackID(int *list, int count, char *sendbuf, char *ID){
// Fill in the phase ID values from neighboring processors
// This packs up the values that need to be sent from one processor to another
int idx,n;
for (idx=0; idx<count; idx++){
n = list[idx];
sendbuf[idx] = ID[n];
}
}
//***************************************************************************************
inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
// Fill in the phase ID values from neighboring processors
// This unpacks the values once they have been recieved from neighbors
int idx,n;
for (idx=0; idx<count; idx++){
n = list[idx];
ID[n] = recvbuf[idx];
}
}
//***************************************************************************************
inline void ZeroHalo(double *Data, int Nx, int Ny, int Nz)
{
int i,j,k,n;
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
i=0;
n = k*Nx*Ny+j*Nx+i;
Data[2*n] = 0.0;
Data[2*n+1] = 0.0;
i=Nx-1;
n = k*Nx*Ny+j*Nx+i;
Data[2*n] = 0.0;
Data[2*n+1] = 0.0;
}
}
for (k=0;k<Nz;k++){
for (i=0;i<Nx;i++){
j=0;
n = k*Nx*Ny+j*Nx+i;
Data[2*n] = 0.0;
Data[2*n+1] = 0.0;
j=Ny-1;
n = k*Nx*Ny+j*Nx+i;
Data[2*n] = 0.0;
Data[2*n+1] = 0.0;
}
}
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
k=0;
n = k*Nx*Ny+j*Nx+i;
Data[2*n] = 0.0;
Data[2*n+1] = 0.0;
k=Nz-1;
n = k*Nx*Ny+j*Nx+i;
Data[2*n] = 0.0;
Data[2*n+1] = 0.0;
}
}
}
//***************************************************************************************
int main(int argc, char **argv)
{
// Initialize MPI
int provided_thread_support = -1;
MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
MPI_Comm comm;
MPI_Comm_dup(MPI_COMM_WORLD,&comm);
int rank = comm_rank(comm);
int nprocs = comm_size(comm);
if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE )
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
// parallel domain size (# of sub-domains)
int nprocx,nprocy,nprocz;
int iproc,jproc,kproc;
MPI_Request req1[18],req2[18];
MPI_Status stat1[18],stat2[18];
if (rank == 0){
printf("********************************************************\n");
printf("Running Color LBM \n");
printf("********************************************************\n");
}
PROFILE_ENABLE(1);
//PROFILE_ENABLE_TRACE();
//PROFILE_ENABLE_MEMORY();
PROFILE_SYNCHRONIZE();
PROFILE_START("Main");
Utilities::setErrorHandlers();
// Variables that specify the computational domain
string FILENAME;
unsigned int nBlocks, nthreads;
int Nx,Ny,Nz; // local sub-domain size
int nspheres; // number of spheres in the packing
double Lx,Ly,Lz; // Domain length
double D = 1.0; // reference length for non-dimensionalization
// Color Model parameters
int timestepMax;
double tau,Fx,Fy,Fz,tol,err;
double alpha, beta;
double das, dbs, phi_s;
double din,dout;
double wp_saturation;
int BoundaryCondition;
int InitialCondition;
// bool pBC,Restart;
int i,j,k;
// pmmc threshold values
//double fluid_isovalue,solid_isovalue;
//fluid_isovalue = 0.0;
//solid_isovalue = 0.0;
int RESTART_INTERVAL=20000;
//int ANALYSIS_)INTERVAL=1000;
int BLOB_ANALYSIS_INTERVAL=1000;
int timestep = -1;
if (rank==0){
//.............................................................
// READ SIMULATION PARMAETERS FROM INPUT FILE
//.............................................................
ifstream input("Color.in");
if (input.is_open()){
// Line 1: model parameters (tau, alpha, beta, das, dbs)
input >> tau; // Viscosity parameter
input >> alpha; // Surface Tension parameter
input >> beta; // Width of the interface
input >> phi_s; // value of phi at the solid surface
// Line 2: wetting phase saturation to initialize
input >> wp_saturation;
// Line 3: External force components (Fx,Fy, Fz)
input >> Fx;
input >> Fy;
input >> Fz;
// Line 4: Pressure Boundary conditions
input >> InitialCondition;
input >> BoundaryCondition;
input >> din;
input >> dout;
// Line 5: time-stepping criteria
input >> timestepMax; // max no. of timesteps
input >> RESTART_INTERVAL; // restart interval
input >> tol; // error tolerance
// Line 6: Analysis options
input >> BLOB_ANALYSIS_INTERVAL; // interval to analyze blob states
//.............................................................
}
else{
// Set default values
// Print warning
printf("WARNING: No input file provided (Color.in is missing)! Default parameters will be used. \n");
tau = 1.0;
alpha=0.005;
beta= 0.9;
Fx = Fy = Fz = 0.0;
InitialCondition=0;
BoundaryCondition=0;
din=dout=1.0;
timestepMax=0;
}
//.......................................................................
// Reading the domain information file
//.......................................................................
ifstream domain("Domain.in");
if (input.is_open()){
domain >> nprocx;
domain >> nprocy;
domain >> nprocz;
domain >> Nx;
domain >> Ny;
domain >> Nz;
domain >> nspheres;
domain >> Lx;
domain >> Ly;
domain >> Lz;
//.......................................................................
}
else{
// Set default values
// Print warning
printf("WARNING: No input file provided (Domain.in is missing)! Default parameters will be used. \n");
nprocx=nprocy=nprocz=1;
Nx=Ny=Nz=10;
nspheres=0;
Lx=Ly=Lz=1.0;
}
}
// **************************************************************
// Broadcast simulation parameters from rank 0 to all other procs
MPI_Barrier(comm);
//.................................................
MPI_Bcast(&tau,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&alpha,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&beta,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&das,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&dbs,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&phi_s,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&wp_saturation,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&BoundaryCondition,1,MPI_INT,0,comm);
MPI_Bcast(&InitialCondition,1,MPI_INT,0,comm);
MPI_Bcast(&din,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&dout,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Fx,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Fy,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Fz,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&timestepMax,1,MPI_INT,0,comm);
MPI_Bcast(&RESTART_INTERVAL,1,MPI_INT,0,comm);
MPI_Bcast(&tol,1,MPI_DOUBLE,0,comm);
// Computational domain
MPI_Bcast(&Nx,1,MPI_INT,0,comm);
MPI_Bcast(&Ny,1,MPI_INT,0,comm);
MPI_Bcast(&Nz,1,MPI_INT,0,comm);
MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm);
//.................................................
// Get the rank info
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
MPI_Barrier(comm);
// **************************************************************
// **************************************************************
double Ps = -(das-dbs)/(das+dbs);
double rlxA = 1.f/tau;
double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA);
//double xIntPos = log((1.0+phi_s)/(1.0-phi_s))/(2.0*beta);
// Set the density values inside the solid based on the input value phi_s
das = (phi_s+1.0)*0.5;
dbs = 1.0 - das;
if (nprocs != nprocx*nprocy*nprocz){
printf("nprocx = %i \n",nprocx);
printf("nprocy = %i \n",nprocy);
printf("nprocz = %i \n",nprocz);
INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!");
}
if (rank==0){
printf("********************************************************\n");
printf("tau = %f \n", tau);
printf("alpha = %f \n", alpha);
printf("beta = %f \n", beta);
printf("das = %f \n", das);
printf("dbs = %f \n", dbs);
printf("gamma_{wn} = %f \n", 5.796*alpha);
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",Nx,Ny,Nz);
printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
if (BoundaryCondition==0) printf("Periodic boundary conditions will applied \n");
if (BoundaryCondition==1) printf("Pressure boundary conditions will be applied \n");
if (BoundaryCondition==2) printf("Velocity boundary conditions will be applied \n");
if (BoundaryCondition==3) printf("Dynamic pressure boundary conditions will be applied \n");
if (InitialCondition==0) printf("Initial conditions assigned from phase ID file \n");
if (InitialCondition==1) printf("Initial conditions assigned from restart file \n");
printf("********************************************************\n");
}
// Initialized domain and averaging framework for Two-Phase Flow
bool pBC,velBC;
if (BoundaryCondition==1 || BoundaryCondition==3)
pBC=true;
else pBC=false;
if (BoundaryCondition==2) velBC=true;
else velBC=false;
bool Restart;
if (InitialCondition==1) Restart=true;
else Restart=false;
NULL_USE(pBC); NULL_USE(velBC);
// Full domain used for averaging (do not use mask for analysis)
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
// std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
TwoPhase Averages(Dm);
Dm.CommInit(comm);
// Mask that excludes the solid phase
Domain Mask(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
MPI_Barrier(comm);
Nx+=2; Ny+=2; Nz += 2;
//Nx = Ny = Nz; // Cubic domain
int N = Nx*Ny*Nz;
int dist_mem_size = N*sizeof(double);
//.......................................................................
if (rank == 0) printf("Read input media... \n");
//.......................................................................
//.......................................................................
// Filenames used
char LocalRankString[8];
char LocalRankFilename[40];
char LocalRestartFile[40];
char tmpstr[10];
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
// printf("Local File Name = %s \n",LocalRankFilename);
// .......... READ THE INPUT FILE .......................................
// char value;
char *id;
id = new char[N];
int sum = 0;
double sum_local;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
double porosity, pore_vol;
//...........................................................................
if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
//.......................................................................
sprintf(LocalRankString,"%05d",rank);
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
// WriteLocalSolidID(LocalRankFilename, id, N);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
ReadBinaryFile(LocalRankFilename, Averages->SDs.data(), N);
MPI_Barrier(comm);
if (rank == 0) cout << "Domain set." << endl;
//.......................................................................
// Assign the phase ID field based on the signed distance
//.......................................................................
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
id[n] = 0;
}
}
}
sum=0;
pore_vol = 0.0;
for ( k=0;k<Nz;k++){
for ( j=0;j<Ny;j++){
for ( i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
if (Averages->SDs(n) > 0.0){
id[n] = 2;
}
// compute the porosity (actual interface location used)
if (Averages->SDs(n) > 0.0){
sum++;
}
}
}
}
if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
sprintf(LocalRankFilename,"ID.%05i",rank);
size_t readID;
FILE *IDFILE = fopen(LocalRankFilename,"rb");
if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx");
readID=fread(id,1,N,IDFILE);
if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID (rank=%i) \n",rank);
fclose(IDFILE);
// Set up kstart, kfinish so that the reservoirs are excluded from averaging
int kstart,kfinish;
kstart = 1;
kfinish = Nz-1;
if (BoundaryCondition > 0 && Dm.kproc==0) kstart = 4;
if (BoundaryCondition > 0 && Dm.kproc==nprocz-1) kfinish = Nz-4;
// Compute the pore volume
sum_local = 0.0;
for ( k=kstart;k<kfinish;k++){
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
int n = k*Nx*Ny+j*Nx+i;
if (id[n] > 0){
sum_local += 1.0;
}
}
}
}
MPI_Allreduce(&sum_local,&pore_vol,1,MPI_DOUBLE,MPI_SUM,comm);
porosity = pore_vol*iVol_global;
if (rank==0) printf("Media porosity = %f \n",porosity);
//.........................................................
// If external boundary conditions are applied remove solid
if (BoundaryCondition > 0 && Dm.kproc == 0){
for (k=0; k<3; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
//id[n] = 1;
Averages->SDs(n) = max(Averages->SDs(n),1.0*(2.5-k));
}
}
}
}
if (BoundaryCondition > 0 && Dm.kproc == nprocz-1){
for (k=Nz-3; k<Nz; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
//id[n] = 2;
Averages->SDs(n) = max(Averages->SDs(n),1.0*(k-Nz+2.5));
}
}
}
}
//.........................................................
// don't perform computations at the eight corners
id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0;
id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
//.........................................................
// Initialize communication structures in averaging domain
for (i=0; i<Mask.Nx*Mask.Ny*Mask.Nz; i++) Mask.id[i] = id[i];
Mask.CommInit(comm);
//...........................................................................
if (rank==0) printf ("Create ScaLBL_Communicator \n");
// Create a communicator for the device
ScaLBL_Communicator ScaLBL_Comm(Mask);
// set reservoirs
if (BoundaryCondition > 0){
for ( k=0;k<Nz;k++){
for ( j=0;j<Ny;j++){
for ( i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
if (Dm.kproc==0 && k==0) id[n]=1;
if (Dm.kproc==0 && k==1) id[n]=1;
if (Dm.kproc==nprocz-1 && k==Nz-2) id[n]=2;
if (Dm.kproc==nprocz-1 && k==Nz-1) id[n]=2;
Mask.id[n] = id[n];
}
}
}
}
//...........device phase ID.................................................
if (rank==0) printf ("Copy phase ID to device \n");
char *ID;
AllocateDeviceMemory((void **) &ID, N); // Allocate device memory
// Don't compute in the halo
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
if (i==0 || i==Nx-1 || j==0 || j==Ny-1 || k==0 || k==Nz-1) id[n] = 0;
}
}
}
// Copy to the device
CopyToDevice(ID, id, N);
DeviceBarrier();
//...........................................................................
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................
// LBM variables
if (rank==0) printf ("Allocating distributions \n");
//......................device distributions.................................
double *f_even,*f_odd;
double *A_even,*A_odd,*B_even,*B_odd;
//...........................................................................
AllocateDeviceMemory((void **) &f_even, 10*dist_mem_size); // Allocate device memory
AllocateDeviceMemory((void **) &f_odd, 9*dist_mem_size); // Allocate device memory
AllocateDeviceMemory((void **) &A_even, 4*dist_mem_size); // Allocate device memory
AllocateDeviceMemory((void **) &A_odd, 3*dist_mem_size); // Allocate device memory
AllocateDeviceMemory((void **) &B_even, 4*dist_mem_size); // Allocate device memory
AllocateDeviceMemory((void **) &B_odd, 3*dist_mem_size); // Allocate device memory
//...........................................................................
double *Phi,*Den;
double *ColorGrad, *Velocity, *Pressure, *dvcSignDist;
//...........................................................................
AllocateDeviceMemory((void **) &Phi, dist_mem_size);
AllocateDeviceMemory((void **) &Pressure, dist_mem_size);
AllocateDeviceMemory((void **) &dvcSignDist, dist_mem_size);
AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size);
//...........................................................................
// Copy signed distance for device initialization
CopyToDevice(dvcSignDist, Averages->SDs.data(), dist_mem_size);
//...........................................................................
int logcount = 0; // number of surface write-outs
//...........................................................................
// MAIN VARIABLES INITIALIZED HERE
//...........................................................................
//...........................................................................
//...........................................................................
if (rank==0) printf("Setting the distributions, size = %i\n", N);
//...........................................................................
DeviceBarrier();
InitD3Q19(ID, f_even, f_odd, Nx, Ny, Nz);
InitDenColor(ID, Den, Phi, das, dbs, Nx, Ny, Nz);
DeviceBarrier();
//......................................................................
if (Restart == true){
if (rank==0){
printf("Reading restart file! \n");
ifstream restart("Restart.txt");
if (restart.is_open()){
restart >> timestep;
printf("Restarting from timestep =%i \n",timestep);
}
else{
printf("WARNING:No Restart.txt file, setting timestep=0 \n");
timestep=0;
}
}
MPI_Bcast(&timestep,1,MPI_INT,0,comm);
// Read in the restart file to CPU buffers
double *cDen = new double[2*N];
double *cDistEven = new double[10*N];
double *cDistOdd = new double[9*N];
ReadCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
// Copy the restart data to the GPU
CopyToDevice(f_even,cDistEven,10*N*sizeof(double));
CopyToDevice(f_odd,cDistOdd,9*N*sizeof(double));
CopyToDevice(Den,cDen,2*N*sizeof(double));
DeviceBarrier();
delete [] cDen;
delete [] cDistEven;
delete [] cDistOdd;
MPI_Barrier(comm);
}
//......................................................................
InitD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
InitD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
DeviceBarrier();
MPI_Barrier(comm);
//.......................................................................
// Once phase has been initialized, map solid to account for 'smeared' interface
for (i=0; i<N; i++) Averages->SDs(i) -= (1.0);
// Make sure the id match for the two domains
for (i=0; i<N; i++) Dm.id[i] = Mask.id[i];
//.......................................................................
// Finalize setup for averaging domain
//Averages->SetupCubes(Mask);
Averages->UpdateSolid();
//.......................................................................
//*************************************************************************
// Compute the phase indicator field and reset Copy, Den
//*************************************************************************
ComputePhi(ID, Phi, Den, N);
//*************************************************************************
DeviceBarrier();
ScaLBL_Comm.SendHalo(Phi);
ScaLBL_Comm.RecvHalo(Phi);
DeviceBarrier();
MPI_Barrier(comm);
//*************************************************************************
if (rank==0 && BoundaryCondition==1){
printf("Setting inlet pressure = %f \n", din);
printf("Setting outlet pressure = %f \n", dout);
}
if (BoundaryCondition==1 && Mask.kproc == 0) {
PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (rank==0 && BoundaryCondition==2){
printf("Setting inlet velocity = %f \n", din);
printf("Setting outlet velocity = %f \n", dout);
}
if (BoundaryCondition==2 && Mask.kproc == 0) {
ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
//ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
SetPhiSlice_z(Phi,1.0,Nx,Ny,Nz,0);
}
if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
//ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
SetPhiSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
}
// Set dynamic pressure boundary conditions
double dp, slope;
dp = slope = 0.0;
if (BoundaryCondition==3){
slope = (dout-din)/timestepMax;
dp = din + timestep*slope;
if (rank==0) printf("Change in pressure / time =%f \n",slope);
// set the initial value
din = 1.0+0.5*dp;
dout = 1.0-0.5*dp;
// set the initial boundary conditions
if (Mask.kproc == 0) {
PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (Mask.kproc == nprocz-1){
PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
}
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
ComputeVelocityD3Q19(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Averages->Phase_tplus.data(),Phi,N*sizeof(double));
//...........................................................................
//...........................................................................
// Copy the data for for the analysis timestep
//...........................................................................
// Copy the phase from the GPU -> CPU
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
CopyToHost(Averages->Press.data(),Pressure,N*sizeof(double));
CopyToHost(Averages->Vel_x.data(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages->Vel_y.data(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages->Vel_z.data(),&Velocity[2*N],N*sizeof(double));
//...........................................................................
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("No. of timesteps: %i \n", timestepMax);
//.......create and start timer............
double starttime,stoptime,cputime;
DeviceBarrier();
MPI_Barrier(comm);
starttime = MPI_Wtime();
//.........................................
err = 1.0;
double sat_w_previous = 1.01; // slightly impossible value!
if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol);
// Create the MeshDataStruct
fillHalo<double> fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
std::vector<IO::MeshDataStruct> meshData(1);
meshData[0].meshName = "domain";
meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );
std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
std::shared_ptr<IO::Variable> PressVar( new IO::Variable() );
std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
PhaseVar->name = "phase";
PhaseVar->type = IO::VolumeVariable;
PhaseVar->dim = 1;
PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(PhaseVar);
PressVar->name = "Pressure";
PressVar->type = IO::VolumeVariable;
PressVar->dim = 1;
PressVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(PressVar);
SignDistVar->name = "SignDist";
SignDistVar->type = IO::VolumeVariable;
SignDistVar->dim = 1;
SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(SignDistVar);
BlobIDVar->name = "BlobID";
BlobIDVar->type = IO::VolumeVariable;
BlobIDVar->dim = 1;
BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(BlobIDVar);
//************ MAIN ITERATION LOOP ***************************************/
PROFILE_START("Loop");
while (timestep < timestepMax && err > tol ) {
//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
PROFILE_START("Update");
//*************************************************************************
// Fused Color Gradient and Collision
//*************************************************************************
ColorCollideOpt( ID,f_even,f_odd,Phi,ColorGrad,
Velocity,Nx,Ny,Nz,rlxA,rlxB,alpha,beta,Fx,Fy,Fz);
//*************************************************************************
DeviceBarrier();
//*************************************************************************
// Pack and send the D3Q19 distributions
ScaLBL_Comm.SendD3Q19(f_even, f_odd);
//*************************************************************************
//*************************************************************************
// Carry out the density streaming step for mass transport
//*************************************************************************
MassColorCollideD3Q7(ID, A_even, A_odd, B_even, B_odd, Den, Phi,
ColorGrad, Velocity, beta, N, pBC);
//*************************************************************************
DeviceBarrier();
MPI_Barrier(comm);
//*************************************************************************
// Swap the distributions for momentum transport
//*************************************************************************
SwapD3Q19(ID, f_even, f_odd, Nx, Ny, Nz);
//*************************************************************************
DeviceBarrier();
MPI_Barrier(comm);
//*************************************************************************
// Wait for communications to complete and unpack the distributions
ScaLBL_Comm.RecvD3Q19(f_even, f_odd);
//*************************************************************************
DeviceBarrier();
//*************************************************************************
// Pack and send the D3Q7 distributions
ScaLBL_Comm.BiSendD3Q7(A_even, A_odd, B_even, B_odd);
//*************************************************************************
DeviceBarrier();
SwapD3Q7(ID, A_even, A_odd, Nx, Ny, Nz);
SwapD3Q7(ID, B_even, B_odd, Nx, Ny, Nz);
DeviceBarrier();
MPI_Barrier(comm);
//*************************************************************************
// Wait for communication and unpack the D3Q7 distributions
ScaLBL_Comm.BiRecvD3Q7(A_even, A_odd, B_even, B_odd);
//*************************************************************************
DeviceBarrier();
//..................................................................................
ComputeDensityD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
ComputeDensityD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
//*************************************************************************
// Compute the phase indicator field
//*************************************************************************
DeviceBarrier();
MPI_Barrier(comm);
ComputePhi(ID, Phi, Den, N);
//*************************************************************************
ScaLBL_Comm.SendHalo(Phi);
DeviceBarrier();
ScaLBL_Comm.RecvHalo(Phi);
//*************************************************************************
DeviceBarrier();
// Pressure boundary conditions
if (BoundaryCondition==1 && Mask.kproc == 0) {
PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
// Velocity boundary conditions
if (BoundaryCondition==2 && Mask.kproc == 0) {
ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
//ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
SetPhiSlice_z(Phi,1.0,Nx,Ny,Nz,0);
}
if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
//ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
SetPhiSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
}
if (BoundaryCondition==3){
// Increase the pressure difference
dp += slope;
din = 1.0+0.5*dp;
dout = 1.0-0.5*dp;
// set the initial boundary conditions
if (Mask.kproc == 0) {
PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (Mask.kproc == nprocz-1){
PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
}
//...................................................................................
MPI_Barrier(comm);
PROFILE_STOP("Update");
// Timestep completed!
timestep++;
// Run the analysis
//if (timestep > 5)
//...................................................................
if (timestep%1000 == 995){
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Averages.Phase_tplus,Phi,N*sizeof(double));
// Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
//...........................................................................
}
if (timestep%1000 == 0){
//...........................................................................
// Copy the data for for the analysis timestep
//...........................................................................
// Copy the phase from the GPU -> CPU
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Averages.Phase,Phi,N*sizeof(double));
CopyToHost(Averages.Press,Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x,&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y,&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z,&Velocity[2*N],N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
}
if (timestep%1000 == 5){
//...........................................................................
// Copy the phase indicator field for the later timestep
DeviceBarrier();
CopyToHost(Averages.Phase_tminus,Phi,N*sizeof(double));
// Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus);
//....................................................................
Averages.Initialize();
Averages.ComputeDelPhi();
Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn);
Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus);
Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus);
Averages.UpdateMeshValues();
Averages.ComputeLocal();
Averages.Reduce();
Averages.PrintAll(timestep);
*/ //....................................................................
}
if (timestep%RESTART_INTERVAL == 0){
if (pBC){
//err = fabs(sat_w - sat_w_previous);
//sat_w_previous = sat_w;
if (rank==0) printf("Timestep %i: change in saturation since last checkpoint is %f \n", timestep, err);
}
else{
// Not clear yet
}
// Copy the data to the CPU
CopyToHost(cDistEven,f_even,10*N*sizeof(double));
CopyToHost(cDistOdd,f_odd,9*N*sizeof(double));
CopyToHost(cDen,Den,2*N*sizeof(double));
// Read in the restart file to CPU buffers
WriteCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
}
// Save the timers
if ( timestep%50==0 )
PROFILE_SAVE("lbpm_color_simulator",1);
}
PROFILE_STOP("Loop");
//************************************************************************
DeviceBarrier();
MPI_Barrier(comm);
stoptime = MPI_Wtime();
if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep
cputime = (stoptime - starttime)/timestep;
// Performance obtained from each node
double MLUPS = double(Nx*Ny*Nz)/cputime/1000000;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("CPU time = %f \n", cputime);
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("********************************************************\n");
// ************************************************************************
/* // Perform component averaging and write tcat averages
Averages->Initialize();
Averages->ComponentAverages();
Averages->SortBlobs();
Averages->PrintComponents(timestep);
// ************************************************************************
int NumberComponents_NWP = ComputeGlobalPhaseComponent(Mask.Nx-2,Mask.Ny-2,Mask.Nz-2,Mask.rank_info,Averages->PhaseID,1,Averages->Label_NWP);
printf("Number of non-wetting phase components: %i \n ",NumberComponents_NWP);
DeviceBarrier();
CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
*/
/* Averages->WriteSurfaces(0);
sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
FILE *PHASE;
PHASE = fopen(LocalRankFilename,"wb");
fwrite(Averages->SDn.data(),8,N,PHASE);
fclose(PHASE);
*/
/* sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
FILE *PRESS;
PRESS = fopen(LocalRankFilename,"wb");
fwrite(Averages->Press.data(),8,N,PRESS);
fclose(PRESS);
CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
double * Grad;
Grad = new double [3*N];
CopyToHost(Grad,ColorGrad,3*N*sizeof(double));
sprintf(LocalRankFilename,"%s%s","ColorGrad.",LocalRankString);
FILE *GRAD;
GRAD = fopen(LocalRankFilename,"wb");
fwrite(Grad,8,3*N,GRAD);
fclose(GRAD);
*/
PROFILE_STOP("Main");
PROFILE_SAVE("lbpm_color_simulator",1);
// ****************************************************
MPI_Barrier(comm);
} // Limit scope so variables that contain communicators will free before MPI_Finialize
MPI_Comm_free(&comm);
MPI_Finalize();
}

View File

@ -229,7 +229,7 @@ int main(int argc, char **argv)
printf("Force(y) = %f \n", Fy);
printf("Force(z) = %f \n", Fz);
printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
printf("Process grid = %i x %i x %i\n",nprocx,nprocy,nprocz);
printf("********************************************************\n");
}
@ -247,7 +247,6 @@ int main(int argc, char **argv)
MPI_Barrier(comm);
Nx += 2; Ny += 2; Nz += 2;
//Nx = Ny = Nz; // Cubic domain
int N = Nx*Ny*Nz;
int dist_mem_size = N*sizeof(double);

View File

@ -33,6 +33,7 @@ int main(int argc, char **argv)
printf("Solid Label: %i \n",SOLID);
printf("NWP Label: %i \n",NWP);
}
{
//.......................................................................
// Reading the domain information file
//.......................................................................

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
//.......................................................................
@ -442,6 +442,7 @@ int main(int argc, char **argv)
Averages.SortBlobs();
Averages.PrintComponents(timestep);
*/
}
MPI_Barrier(comm);
MPI_Finalize();
return 0;

View File

@ -38,7 +38,9 @@ int main(int argc, char **argv)
MPI_Request req1[18],req2[18];
MPI_Status stat1[18],stat2[18];
int ORIENTATION;
int ORIENTATION=2; //default: the tube is aligned with Z axis
//ORIENTATION = 0: tube is aligned with X axis
//ORIENTATION = 1: tube is aligned with Y axis
double TubeWidth =15.0;
int BC;
int BubbleTop,BubbleBottom;
@ -161,34 +163,66 @@ int main(int argc, char **argv)
if (ORIENTATION==0){
// square capillary tube aligned with the x direction
Averages.SDs(i,j,k) = TubeWidth/2 - fabs(i-0.5*Ny);
Averages.SDs(i,j,k) = min(Averages.SDs(i,j,k),TubeWidth/2-fabs(j-0.5*Nz));
Averages.SDs(i,j,k) = TubeWidth/2 - fabs(j-0.5*Ny);
Averages.SDs(i,j,k) = min(Averages.SDs(i,j,k),TubeWidth/2-fabs(k-0.5*Nz));
// Initialize phase positions
if (Averages.SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Dm.iproc*Nx+k<BubbleBottom){
id[n] = 2;
sum++;
}
else if (Dm.iproc*Nx+k<BubbleTop){
id[n] = 1;
sum++;
}
else{
id[n] = 2;
sum++;
}
}
else if (ORIENTATION==1){
// square capillary tube aligned with the y direction
Averages.SDs(i,j,k) = TubeWidth/2 - fabs(i-0.5*Nx);
Averages.SDs(i,j,k) = min(Averages.SDs(i,j,k),TubeWidth/2-fabs(j-0.5*Nz));
Averages.SDs(i,j,k) = min(Averages.SDs(i,j,k),TubeWidth/2-fabs(k-0.5*Nz));
// Initialize phase positions
if (Averages.SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Dm.jproc*Ny+k<BubbleBottom){
id[n] = 2;
sum++;
}
else if (Dm.jproc*Ny+k<BubbleTop){
id[n] = 1;
sum++;
}
else{
id[n] = 2;
sum++;
}
}
else {
else { //
// square capillary tube aligned with the z direction
Averages.SDs(i,j,k) = TubeWidth/2 - fabs(i-0.5*Nx);
Averages.SDs(i,j,k) = min(Averages.SDs(i,j,k),TubeWidth/2-fabs(j-0.5*Ny));
}
// Initialize phase positions
if (Averages.SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Dm.kproc*Nz+k<BubbleBottom){
id[n] = 2;
sum++;
}
else if (Dm.kproc*Nz+k<BubbleTop){
id[n] = 1;
sum++;
}
else{
id[n] = 2;
sum++;
// Initialize phase positions
if (Averages.SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Dm.kproc*Nz+k<BubbleBottom){
id[n] = 2;
sum++;
}
else if (Dm.kproc*Nz+k<BubbleTop){
id[n] = 1;
sum++;
}
else{
id[n] = 2;
sum++;
}
}
}
}

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

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

View File

@ -0,0 +1,104 @@
#!/usr/bin/env python
import numpy as np
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
# 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')
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
image_format = '.tiff' # format of experimental images
# Read other input parameters
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
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
# Load a group of image files with 'base_name' in the names of files
# e.g. 'Micromodel_1.tiff' etc.
file_input_group = glob('*'+Para['base_name']+'*'+image_format) # need to match the image format
# 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 "Processing image "+file_input_single+" now..."
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
# 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
# Write out the segmented data (in binary format)
print "The segmented image is written to "+output_file_name+" now..."
ID.tofile(output_file_name)
# NOTE when you want to load the binary file
# Do as follows:
# ID_reload = np.fromfile(file_name,dtype=np.uint8)
#end for
#end if
# In case you want to test the segmentation
#plt.figure(1)
#plt.subplot(1,2,1)
#plt.title('original RGB figure')
#plt.pcolormesh(im_array);
#plt.axis('equal')
#plt.colorbar()
#plt.subplot(1,2,2)
#plt.title('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])
#plt.axis('equal')
#plt.show()

View File

@ -0,0 +1,68 @@
#!/usr/bin/env python
import numpy as np
from PIL import Image # import Python Imaging Library (PIL)
def convert_image_to_array(file_name,imin,imax):
# file_name: the name of the input image (a string)
# imin: (left) boundary that defines the output array
# imax: (right) boundary that defines the output array
# Note: assume a grayscale image (whose R, G, and B values are the same)
# Return an numpy array with RGB values that has the same size as the wanted micromodel
im = Image.open(file_name).convert('RGB')
im_array = np.array(im) # 'im_array' has dimension (height, width, 3)
# Grayscale input image is assumed, thus its R, G, and B values are the same
im_array = im_array[:,:,0] # 2D array with dimension (height, width)
im_array = im_array[:,imin:imax] # Chop the original image to a desirable size
#im_array = np.flipud(im_array[:,imin:imax]) # You might want to flip the Z-axis
return im_array
#end def
def read_input_parameters(input_file_name):
# input_file_name is a string
# The *.in file has the following lines
# line 0: the base name of the experimental images
# line 1: imin imax
# line 2: imin_b imax_b
# line 3: top_layer bottom_layer DEPTH
# line 4: threshold_nw threshold_s
# Note: 'threshold_nw' means: NW phase: RGB values > threshold_nw
# 'threshold_s' means: solid: RGB values < threshold_s
f = open(input_file_name,'r') # read-only
lines = f.readlines()
output_file = {'base_name':lines[0].splitlines()[0]}
line1_array = np.fromstring(lines[1].splitlines()[0],dtype=np.int32,sep=' ')
line2_array = np.fromstring(lines[2].splitlines()[0],dtype=np.int32,sep=' ')
line3_array = np.fromstring(lines[3].splitlines()[0],dtype=np.int32,sep=' ')
line4_array = np.fromstring(lines[4].splitlines()[0],dtype=np.int32,sep=' ')
output_file['imin']=line1_array[0]
output_file['imax']=line1_array[1]
output_file['imin_b']=line2_array[0]
output_file['imax_b']=line2_array[1]
output_file['top_layer']=line3_array[0]
output_file['bot_layer']=line3_array[1]
output_file['DEPTH']=line3_array[2]
output_file['threshold_nw']=line4_array[0]
output_file['threshold_s'] =line4_array[1]
f.close()
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