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

Fixing compile warnings

Browse Source
This commit is contained in:
Mark Berrill 2020-01-22 12:01:29 -05:00
parent f17bccb389
commit 78c2e710b9
37 changed files with 215 additions and 334 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
IO/netcdf.cpp
View File

@ -189,6 +189,7 @@ std::vector<size_t> getAttDim( int fid, const std::string& att )
{ {
std::vector<size_t> dim(1,0); std::vector<size_t> dim(1,0);
int err = nc_inq_attlen( fid, NC_GLOBAL, att.c_str(), dim.data() ); int err = nc_inq_attlen( fid, NC_GLOBAL, att.c_str(), dim.data() );
CHECK_NC_ERR( err );
return dim; return dim;
} }
std::vector<std::string> getVarNames( int fid ) std::vector<std::string> getVarNames( int fid )

10
analysis/SubPhase.cpp
View File

@ -169,7 +169,6 @@ void SubPhase::Basic(){
nb.reset(); wb.reset(); nb.reset(); wb.reset();
double nA,nB;
double count_w = 0.0; double count_w = 0.0;
double count_n = 0.0; double count_n = 0.0;
@ -297,8 +296,8 @@ void SubPhase::Basic(){
double saturation=gwb.V/(gwb.V + gnb.V); double saturation=gwb.V/(gwb.V + gnb.V);
double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M / Dm->Volume; double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M / Dm->Volume;
double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume; double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume;
double total_flow_rate = water_flow_rate + not_water_flow_rate; //double total_flow_rate = water_flow_rate + not_water_flow_rate;
double fractional_flow= water_flow_rate / total_flow_rate; //double fractional_flow = water_flow_rate / total_flow_rate;
double h = Dm->voxel_length; double h = Dm->voxel_length;
double krn = h*h*nu_n*not_water_flow_rate / force_mag ; double krn = h*h*nu_n*not_water_flow_rate / force_mag ;
@ -697,7 +696,8 @@ void SubPhase::Full(){
} }
void SubPhase::AggregateLabels(char *FILENAME){ void SubPhase::AggregateLabels( const std::string& filename )
{
int nx = Dm->Nx; int nx = Dm->Nx;
int ny = Dm->Ny; int ny = Dm->Ny;
@ -721,7 +721,7 @@ void SubPhase::AggregateLabels(char *FILENAME){
} }
MPI_Barrier(Dm->Comm); MPI_Barrier(Dm->Comm);
Dm->AggregateLabels(FILENAME); Dm->AggregateLabels( filename );
} }

2
analysis/SubPhase.h
View File

@ -101,7 +101,7 @@ public:
void Basic(); void Basic();
void Full(); void Full();
void Write(int time); void Write(int time);
void AggregateLabels(char *FILENAME); void AggregateLabels( const std::string& filename );
private: private:
FILE *TIMELOG; FILE *TIMELOG;

7
analysis/TwoPhase.cpp
View File

@ -204,6 +204,7 @@ TwoPhase::~TwoPhase()
void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, DoubleArray &DistData) void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, DoubleArray &DistData)
{ {
NULL_USE( Beta );
/*double factor,temp,value; /*double factor,temp,value;
factor=0.5/Beta; factor=0.5/Beta;
// Initialize to -1,1 (segmentation) // Initialize to -1,1 (segmentation)
@ -627,8 +628,8 @@ void TwoPhase::ComputeLocal()
void TwoPhase::AssignComponentLabels() void TwoPhase::AssignComponentLabels()
{ {
int LabelNWP=1; //int LabelNWP=1;
int LabelWP=2; //int LabelWP=2;
// NOTE: labeling the wetting phase components is tricky! One sandstone media had over 800,000 components // NOTE: labeling the wetting phase components is tricky! One sandstone media had over 800,000 components
// NumberComponents_WP = ComputeGlobalPhaseComponent(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->rank_info,PhaseID,LabelWP,Label_WP); // NumberComponents_WP = ComputeGlobalPhaseComponent(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->rank_info,PhaseID,LabelWP,Label_WP);
// treat all wetting phase is connected // treat all wetting phase is connected
@ -1172,6 +1173,8 @@ void TwoPhase::Reduce()
void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT) void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT)
{ {
NULL_USE( viscosity );
NULL_USE( IFT );
awn_global *= D; awn_global *= D;
ans_global *= D; ans_global *= D;
ans_global *= D; ans_global *= D;

6
analysis/dcel.cpp
View File

@ -352,6 +352,8 @@ double DECL::EdgeAngle(int edge)
void Isosurface(DoubleArray &A, const double &v) void Isosurface(DoubleArray &A, const double &v)
{ {
NULL_USE( v );
Point P,Q; Point P,Q;
Point PlaceHolder; Point PlaceHolder;
Point C0,C1,C2,C3,C4,C5,C6,C7; Point C0,C1,C2,C3,C4,C5,C6,C7;
@ -562,7 +564,7 @@ void Isosurface(DoubleArray &A, const double &v)
if (P.z == 1.0 && Q.z == 1.0) HalfEdge[idx_edge][3] = -6; // ghost twin for z=1 face if (P.z == 1.0 && Q.z == 1.0) HalfEdge[idx_edge][3] = -6; // ghost twin for z=1 face
} }
// Find all the angles // Find all the angles
for (int idx=0; idx<EdgeCount; idx++){ /*for (int idx=0; idx<EdgeCount; idx++){
int V1=HalfEdge[idx][0]; int V1=HalfEdge[idx][0];
int V2=HalfEdge[idx][1]; int V2=HalfEdge[idx][1];
int T1= HalfEdge[idx_edge][2]; int T1= HalfEdge[idx_edge][2];
@ -570,7 +572,7 @@ void Isosurface(DoubleArray &A, const double &v)
if (twin == -1){ if (twin == -1){
} }
} }*/
// Map vertices to global coordinates // Map vertices to global coordinates
for (int idx=0;idx<VertexCount;idx++) { for (int idx=0;idx<VertexCount;idx++) {

18
analysis/morphology.cpp
View File

@ -34,9 +34,6 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
int nx = Dm->Nx; int nx = Dm->Nx;
int ny = Dm->Ny; int ny = Dm->Ny;
int nz = Dm->Nz; int nz = Dm->Nz;
int iproc = Dm->iproc();
int jproc = Dm->jproc();
int kproc = Dm->kproc();
int nprocx = Dm->nprocx(); int nprocx = Dm->nprocx();
int nprocy = Dm->nprocy(); int nprocy = Dm->nprocy();
int nprocz = Dm->nprocz(); int nprocz = Dm->nprocz();
@ -122,7 +119,6 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
int sendtag,recvtag; int sendtag,recvtag;
sendtag = recvtag = 7; sendtag = recvtag = 7;
int x,y,z;
int ii,jj,kk; int ii,jj,kk;
int Nx = nx; int Nx = nx;
int Ny = ny; int Ny = ny;
@ -336,9 +332,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
int nx = Dm->Nx; int nx = Dm->Nx;
int ny = Dm->Ny; int ny = Dm->Ny;
int nz = Dm->Nz; int nz = Dm->Nz;
int iproc = Dm->iproc();
int jproc = Dm->jproc();
int kproc = Dm->kproc();
int nprocx = Dm->nprocx(); int nprocx = Dm->nprocx();
int nprocy = Dm->nprocy(); int nprocy = Dm->nprocy();
int nprocz = Dm->nprocz(); int nprocz = Dm->nprocz();
@ -427,7 +420,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
int sendtag,recvtag; int sendtag,recvtag;
sendtag = recvtag = 7; sendtag = recvtag = 7;
int x,y,z;
int ii,jj,kk; int ii,jj,kk;
int Nx = nx; int Nx = nx;
int Ny = ny; int Ny = ny;
@ -693,17 +685,11 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
return final_void_fraction; return final_void_fraction;
} }
double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id, std::shared_ptr<Domain> Dm, double TargetGrowth){ double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id, std::shared_ptr<Domain> Dm, double TargetGrowth)
{
int Nx = Dm->Nx; int Nx = Dm->Nx;
int Ny = Dm->Ny; int Ny = Dm->Ny;
int Nz = Dm->Nz; int Nz = Dm->Nz;
int iproc = Dm->iproc();
int jproc = Dm->jproc();
int kproc = Dm->kproc();
int nprocx = Dm->nprocx();
int nprocy = Dm->nprocy();
int nprocz = Dm->nprocz();
int rank = Dm->rank(); int rank = Dm->rank();
double count=0.0; double count=0.0;

2
analysis/uCT.cpp
View File

@ -157,6 +157,7 @@ void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
// int depth = 5; // int depth = 5;
// float sigsq=0.1; // float sigsq=0.1;
int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq); int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq);
NULL_USE( nlm_count );
fillFloat.fill(NonLocalMean); fillFloat.fill(NonLocalMean);
} }
@ -201,6 +202,7 @@ void refine( const Array<float>& Dist_coarse,
// int depth = 3; // int depth = 3;
// float sigsq = 0.1; // float sigsq = 0.1;
int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq); int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq);
NULL_USE( nlm_count );
fillFloat.fill(NonLocalMean); fillFloat.fill(NonLocalMean);
segment( NonLocalMean, ID, 0.001 ); segment( NonLocalMean, ID, 0.001 );
for (size_t i=0; i<ID.length(); i++) { for (size_t i=0; i<ID.length(); i++) {

4
cmake/libraries.cmake
View File

@ -305,4 +305,8 @@ MACRO ( CONFIGURE_LBPM )
SET( CMAKE_BUILD_WITH_INSTALL_RPATH TRUE ) SET( CMAKE_BUILD_WITH_INSTALL_RPATH TRUE )
SET( CMAKE_INSTALL_RPATH ${CMAKE_INSTALL_RPATH} "${TIMER_DIRECTORY}" "${LBPM_INSTALL_DIR}/lib" ) SET( CMAKE_INSTALL_RPATH ${CMAKE_INSTALL_RPATH} "${TIMER_DIRECTORY}" "${LBPM_INSTALL_DIR}/lib" )
ENDIF() ENDIF()
# Suppress some common warnings
IF ( USING_GCC )
SET( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-reorder -Wno-unused-parameter")
ENDIF()
ENDMACRO () ENDMACRO ()

42
common/Domain.cpp
View File

@ -58,6 +58,9 @@ Domain::Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz,
recvData_xY(NULL), recvData_yZ(NULL), recvData_Xz(NULL), recvData_XY(NULL), recvData_YZ(NULL), recvData_XZ(NULL), recvData_xY(NULL), recvData_yZ(NULL), recvData_Xz(NULL), recvData_XY(NULL), recvData_YZ(NULL), recvData_XZ(NULL),
id(NULL) id(NULL)
{ {
NULL_USE( rnk );
NULL_USE( npy );
NULL_USE( npz );
// set up the neighbor ranks // set up the neighbor ranks
int myrank; int myrank;
MPI_Comm_rank( Comm, &myrank ); MPI_Comm_rank( Comm, &myrank );
@ -241,7 +244,7 @@ void Domain::initialize( std::shared_ptr<Database> db )
INSIST(nprocs == nproc[0]*nproc[1]*nproc[2],"Fatal error in processor count!"); INSIST(nprocs == nproc[0]*nproc[1]*nproc[2],"Fatal error in processor count!");
} }
void Domain::Decomp(std::string Filename) void Domain::Decomp( const std::string& Filename )
{ {
//....................................................................... //.......................................................................
// Reading the domain information file // Reading the domain information file
@ -251,7 +254,6 @@ void Domain::Decomp(std::string Filename)
int nprocs, nprocx, nprocy, nprocz, nx, ny, nz; int nprocs, nprocx, nprocy, nprocz, nx, ny, nz;
int64_t global_Nx,global_Ny,global_Nz; int64_t global_Nx,global_Ny,global_Nz;
int64_t i,j,k,n; int64_t i,j,k,n;
int BC=0;
int64_t xStart,yStart,zStart; int64_t xStart,yStart,zStart;
int checkerSize; int checkerSize;
//int inlet_layers_x, inlet_layers_y, inlet_layers_z; //int inlet_layers_x, inlet_layers_y, inlet_layers_z;
@ -331,7 +333,7 @@ void Domain::Decomp(std::string Filename)
if (RANK==0){ if (RANK==0){
printf("Input media: %s\n",Filename.c_str()); printf("Input media: %s\n",Filename.c_str());
printf("Relabeling %lu values\n",ReadValues.size()); printf("Relabeling %lu values\n",ReadValues.size());
for (int idx=0; idx<ReadValues.size(); idx++){ for (size_t idx=0; idx<ReadValues.size(); idx++){
int oldvalue=ReadValues[idx]; int oldvalue=ReadValues[idx];
int newvalue=WriteValues[idx]; int newvalue=WriteValues[idx];
printf("oldvalue=%d, newvalue =%d \n",oldvalue,newvalue); printf("oldvalue=%d, newvalue =%d \n",oldvalue,newvalue);
@ -374,7 +376,7 @@ void Domain::Decomp(std::string Filename)
n = k*global_Nx*global_Ny+j*global_Nx+i; n = k*global_Nx*global_Ny+j*global_Nx+i;
//char locval = loc_id[n]; //char locval = loc_id[n];
char locval = SegData[n]; char locval = SegData[n];
for (int idx=0; idx<ReadValues.size(); idx++){ for (size_t idx=0; idx<ReadValues.size(); idx++){
signed char oldvalue=ReadValues[idx]; signed char oldvalue=ReadValues[idx];
signed char newvalue=WriteValues[idx]; signed char newvalue=WriteValues[idx];
if (locval == oldvalue){ if (locval == oldvalue){
@ -387,10 +389,10 @@ void Domain::Decomp(std::string Filename)
} }
} }
if (RANK==0){ if (RANK==0){
for (int idx=0; idx<ReadValues.size(); idx++){ for (size_t idx=0; idx<ReadValues.size(); idx++){
long int label=ReadValues[idx]; long int label=ReadValues[idx];
long int count=LabelCount[idx]; long int count=LabelCount[idx];
printf("Label=%d, Count=%d \n",label,count); printf("Label=%ld, Count=%ld \n",label,count);
} }
} }
@ -475,7 +477,7 @@ void Domain::Decomp(std::string Filename)
printf("Checkerboard pattern at y outlet for %i layers \n",outlet_layers_y); printf("Checkerboard pattern at y outlet for %i layers \n",outlet_layers_y);
// use checkerboard pattern // use checkerboard pattern
for (int k = 0; k<global_Nz; k++){ for (int k = 0; k<global_Nz; k++){
for (int j = yStart + ny*nprocy - outlet_layers_y; i < yStart + ny*nprocy; j++){ for (int j = yStart + ny*nprocy - outlet_layers_y; j < yStart + ny*nprocy; j++){
for (int i = 0; i<global_Nx; i++){ for (int i = 0; i<global_Nx; i++){
if ( (i/checkerSize + k/checkerSize)%2 == 0){ if ( (i/checkerSize + k/checkerSize)%2 == 0){
// void checkers // void checkers
@ -587,7 +589,7 @@ void Domain::Decomp(std::string Filename)
MPI_Barrier(Comm); MPI_Barrier(Comm);
} }
void Domain::AggregateLabels(char *FILENAME){ void Domain::AggregateLabels( const std::string& filename ){
int nx = Nx; int nx = Nx;
int ny = Ny; int ny = Ny;
@ -664,8 +666,7 @@ void Domain::AggregateLabels(char *FILENAME){
} }
} }
// write the output // write the output
FILE *OUTFILE; FILE *OUTFILE = fopen(filename.c_str(),"wb");
OUTFILE = fopen(FILENAME,"wb");
fwrite(FullID,1,full_size,OUTFILE); fwrite(FullID,1,full_size,OUTFILE);
fclose(OUTFILE); fclose(OUTFILE);
} }
@ -1145,19 +1146,18 @@ void Domain::CommunicateMeshHalo(DoubleArray &Mesh)
} }
// Ideally stuff below here should be moved somewhere else -- doesn't really belong here // Ideally stuff below here should be moved somewhere else -- doesn't really belong here
void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, int Np) void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, size_t Np)
{ {
int q,n;
double value; double value;
ofstream File(FILENAME,ios::binary); ofstream File(FILENAME,ios::binary);
for (n=0; n<Np; n++){ for (size_t n=0; n<Np; n++){
// Write the two density values // Write the two density values
value = cDen[n]; value = cDen[n];
File.write((char*) &value, sizeof(value)); File.write((char*) &value, sizeof(value));
value = cDen[Np+n]; value = cDen[Np+n];
File.write((char*) &value, sizeof(value)); File.write((char*) &value, sizeof(value));
// Write the even distributions // Write the even distributions
for (q=0; q<19; q++){ for (size_t q=0; q<19; q++){
value = cfq[q*Np+n]; value = cfq[q*Np+n];
File.write((char*) &value, sizeof(value)); File.write((char*) &value, sizeof(value));
} }
@ -1166,16 +1166,15 @@ void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq
} }
void ReadCheckpoint(char *FILENAME, double *cPhi, double *cfq, int Np) void ReadCheckpoint(char *FILENAME, double *cPhi, double *cfq, size_t Np)
{ {
int q=0, n=0;
double value=0; double value=0;
ifstream File(FILENAME,ios::binary); ifstream File(FILENAME,ios::binary);
for (n=0; n<Np; n++){ for (size_t n=0; n<Np; n++){
File.read((char*) &value, sizeof(value)); File.read((char*) &value, sizeof(value));
cPhi[n] = value; cPhi[n] = value;
// Read the distributions // Read the distributions
for (q=0; q<19; q++){ for (size_t q=0; q<19; q++){
File.read((char*) &value, sizeof(value)); File.read((char*) &value, sizeof(value));
cfq[q*Np+n] = value; cfq[q*Np+n] = value;
} }
@ -1183,13 +1182,12 @@ void ReadCheckpoint(char *FILENAME, double *cPhi, double *cfq, int Np)
File.close(); File.close();
} }
void ReadBinaryFile(char *FILENAME, double *Data, int N) void ReadBinaryFile(char *FILENAME, double *Data, size_t N)
{ {
int n;
double value; double value;
ifstream File(FILENAME,ios::binary); ifstream File(FILENAME,ios::binary);
if (File.good()){ if (File.good()){
for (n=0; n<N; n++){ for (size_t n=0; n<N; n++){
// Write the two density values // Write the two density values
File.read((char*) &value, sizeof(value)); File.read((char*) &value, sizeof(value));
Data[n] = value; Data[n] = value;
@ -1197,7 +1195,7 @@ void ReadBinaryFile(char *FILENAME, double *Data, int N)
} }
} }
else { else {
for (n=0; n<N; n++) Data[n] = 1.2e-34; for (size_t n=0; n<N; n++) Data[n] = 1.2e-34;
} }
File.close(); File.close();
} }

10
common/Domain.h
View File

@ -177,11 +177,11 @@ public: // Public variables (need to create accessors instead)
signed char *id; signed char *id;
void ReadIDs(); void ReadIDs();
void Decomp(std::string Filename); void Decomp( const std::string& filename );
void CommunicateMeshHalo(DoubleArray &Mesh); void CommunicateMeshHalo(DoubleArray &Mesh);
void CommInit(); void CommInit();
int PoreCount(); int PoreCount();
void AggregateLabels(char *FILENAME); void AggregateLabels( const std::string& filename );
private: private:
@ -244,10 +244,10 @@ private:
}; };
void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, int Np); void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, size_t Np);
void ReadCheckpoint(char *FILENAME, double *cDen, double *cfq, int Np); void ReadCheckpoint(char *FILENAME, double *cDen, double *cfq, size_t Np);
void ReadBinaryFile(char *FILENAME, double *Data, int N); void ReadBinaryFile(char *FILENAME, double *Data, size_t N);
#endif #endif

6
common/MPI_Helpers.h
View File

@ -188,6 +188,12 @@ inline int sumReduce( MPI_Comm comm, int x )
MPI_Allreduce(&x,&y,1,MPI_INT,MPI_SUM,comm); MPI_Allreduce(&x,&y,1,MPI_INT,MPI_SUM,comm);
return y; return y;
} }
inline long long sumReduce( MPI_Comm comm, long long x )
{
long long y = 0;
MPI_Allreduce(&x,&y,1,MPI_LONG_LONG,MPI_SUM,comm);
return y;
}
inline bool sumReduce( MPI_Comm comm, bool x ) inline bool sumReduce( MPI_Comm comm, bool x )
{ {
int y = sumReduce( comm, x?1:0 ); int y = sumReduce( comm, x?1:0 );

25
common/ScaLBL.cpp
View File

@ -365,7 +365,7 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis
int idx,i,j,k,n; int idx,i,j,k,n;
// Check that Map has size matching sub-domain // Check that Map has size matching sub-domain
if (Map.size(0) != Nx) if ( (int) Map.size(0) != Nx)
ERROR("ScaLBL_Communicator::MemoryOptimizedLayout: Map array dimensions do not match! \n"); ERROR("ScaLBL_Communicator::MemoryOptimizedLayout: Map array dimensions do not match! \n");
// Initialize Map // Initialize Map
@ -1480,7 +1480,6 @@ void ScaLBL_Communicator::RecvHalo(double *data){
void ScaLBL_Communicator::RegularLayout(IntArray map, const double *data, DoubleArray &regdata){ void ScaLBL_Communicator::RegularLayout(IntArray map, const double *data, DoubleArray &regdata){
// Gets data from the device and stores in regular layout // Gets data from the device and stores in regular layout
int i,j,k,n,idx;
int Nx = map.size(0); int Nx = map.size(0);
int Ny = map.size(1); int Ny = map.size(1);
int Nz = map.size(2); int Nz = map.size(2);
@ -1492,11 +1491,10 @@ void ScaLBL_Communicator::RegularLayout(IntArray map, const double *data, Double
double value; double value;
TmpDat = new double [N]; TmpDat = new double [N];
ScaLBL_CopyToHost(&TmpDat[0],&data[0], N*sizeof(double)); ScaLBL_CopyToHost(&TmpDat[0],&data[0], N*sizeof(double));
for (k=0; k<Nz; k++){ for (int k=0; k<Nz; k++){
for (j=0; j<Ny; j++){ for (int j=0; j<Ny; j++){
for (i=0; i<Nx; i++){ for (int i=0; i<Nx; i++){
n=k*Nx*Ny+j*Nx+i; auto idx=map(i,j,k);
idx=map(i,j,k);
if (!(idx<0)){ if (!(idx<0)){
value=TmpDat[idx]; value=TmpDat[idx];
regdata(i,j,k)=value; regdata(i,j,k)=value;
@ -1510,8 +1508,9 @@ void ScaLBL_Communicator::RegularLayout(IntArray map, const double *data, Double
} }
void ScaLBL_Communicator::Color_BC_z(int *Map, double *Phi, double *Den, double vA, double vB){ void ScaLBL_Communicator::Color_BC_z(int *Map, double *Phi, double *Den, double vA, double vB)
double Value=(vA-vB)/(vA+vB); {
//double Value=(vA-vB)/(vA+vB);
if (kproc == 0) { if (kproc == 0) {
// Set the phase indicator field and density on the z inlet // Set the phase indicator field and density on the z inlet
ScaLBL_Color_BC_z(dvcSendList_z, Map, Phi, Den, vA, vB, sendCount_z, N); ScaLBL_Color_BC_z(dvcSendList_z, Map, Phi, Den, vA, vB, sendCount_z, N);
@ -1519,8 +1518,9 @@ void ScaLBL_Communicator::Color_BC_z(int *Map, double *Phi, double *Den, double
} }
} }
void ScaLBL_Communicator::Color_BC_Z(int *Map, double *Phi, double *Den, double vA, double vB){ void ScaLBL_Communicator::Color_BC_Z(int *Map, double *Phi, double *Den, double vA, double vB)
double Value=(vA-vB)/(vA+vB); {
//double Value=(vA-vB)/(vA+vB);
if (kproc == nprocz-1){ if (kproc == nprocz-1){
// Set the phase indicator field and density on the Z outlet // Set the phase indicator field and density on the Z outlet
ScaLBL_Color_BC_Z(dvcSendList_Z, Map, Phi, Den, vA, vB, sendCount_Z, N); ScaLBL_Color_BC_Z(dvcSendList_Z, Map, Phi, Den, vA, vB, sendCount_Z, N);
@ -1528,7 +1528,8 @@ void ScaLBL_Communicator::Color_BC_Z(int *Map, double *Phi, double *Den, double
} }
} }
void ScaLBL_Communicator::D3Q19_Pressure_BC_z(int *neighborList, double *fq, double din, int time){ void ScaLBL_Communicator::D3Q19_Pressure_BC_z(int *neighborList, double *fq, double din, int time)
{
//ScaLBL_D3Q19_Pressure_BC_z(int *LIST,fq,din,Nx,Ny,Nz); //ScaLBL_D3Q19_Pressure_BC_z(int *LIST,fq,din,Nx,Ny,Nz);
if (kproc == 0) { if (kproc == 0) {
if (time%2==0){ if (time%2==0){

53
models/ColorModel.cpp
View File

@ -186,7 +186,6 @@ void ScaLBL_ColorModel::ReadInput(){
if (color_db->keyExists( "image_sequence" )){ if (color_db->keyExists( "image_sequence" )){
auto ImageList = color_db->getVector<std::string>( "image_sequence"); auto ImageList = color_db->getVector<std::string>( "image_sequence");
int IMAGE_INDEX = color_db->getWithDefault<int>( "image_index", 0 ); int IMAGE_INDEX = color_db->getWithDefault<int>( "image_index", 0 );
int IMAGE_COUNT = ImageList.size();
std::string first_image = ImageList[IMAGE_INDEX]; std::string first_image = ImageList[IMAGE_INDEX];
Mask->Decomp(first_image); Mask->Decomp(first_image);
IMAGE_INDEX++; IMAGE_INDEX++;
@ -195,9 +194,9 @@ void ScaLBL_ColorModel::ReadInput(){
// Read the local domain data // Read the local domain data
auto input_id = readMicroCT( *domain_db, MPI_COMM_WORLD ); auto input_id = readMicroCT( *domain_db, MPI_COMM_WORLD );
// Fill the halo (assuming GCW of 1) // Fill the halo (assuming GCW of 1)
array<int,3> size0 = { input_id.size(0), input_id.size(1), input_id.size(2) }; array<int,3> size0 = { (int) input_id.size(0), (int) input_id.size(1), (int) input_id.size(2) };
ArraySize size1 = { Mask->Nx, Mask->Ny, Mask->Nz }; ArraySize size1 = { (size_t) Mask->Nx, (size_t) Mask->Ny, (size_t) Mask->Nz };
ASSERT( size1[0] == size0[0]+2 && size1[1] == size0[1]+2 && size1[2] == size0[2]+2 ); ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 );
fillHalo<signed char> fill( MPI_COMM_WORLD, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 ); fillHalo<signed char> fill( MPI_COMM_WORLD, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
Array<signed char> id_view; Array<signed char> id_view;
id_view.viewRaw( size1, Mask->id ); id_view.viewRaw( size1, Mask->id );
@ -216,7 +215,6 @@ void ScaLBL_ColorModel::ReadInput(){
// Generate the signed distance map // Generate the signed distance map
// Initialize the domain and communication // Initialize the domain and communication
Array<char> id_solid(Nx,Ny,Nz); Array<char> id_solid(Nx,Ny,Nz);
int count = 0;
// Solve for the position of the solid phase // Solve for the position of the solid phase
for (int k=0;k<Nz;k++){ for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){ for (int j=0;j<Ny;j++){
@ -233,7 +231,6 @@ void ScaLBL_ColorModel::ReadInput(){
for (int k=0;k<Nz;k++){ for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){ for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){ for (int i=0;i<Nx;i++){
int n=k*Nx*Ny+j*Nx+i;
// Initialize distance to +/- 1 // Initialize distance to +/- 1
Averages->SDs(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0; Averages->SDs(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
} }
@ -266,7 +263,7 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
double label_count_global[NLABELS]; double label_count_global[NLABELS];
// Assign the labels // Assign the labels
for (int idx=0; idx<NLABELS; idx++) label_count[idx]=0; for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
for (int k=0;k<Nz;k++){ for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){ for (int j=0;j<Ny;j++){
@ -294,7 +291,8 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
// Set Dm to match Mask // Set Dm to match Mask
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i]; for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
for (int idx=0; idx<NLABELS; idx++) label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]); for (size_t idx=0; idx<NLABELS; idx++)
label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
if (rank==0){ if (rank==0){
printf("Component labels: %lu \n",NLABELS); printf("Component labels: %lu \n",NLABELS);
@ -373,16 +371,16 @@ void ScaLBL_ColorModel::Create(){
} }
// check that TmpMap is valid // check that TmpMap is valid
for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){ for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
int n = TmpMap[idx]; auto n = TmpMap[idx];
if (n > Nx*Ny*Nz){ if (n > Nx*Ny*Nz){
printf("Bad value! idx=%i \n"); printf("Bad value! idx=%i \n", n);
TmpMap[idx] = Nx*Ny*Nz-1; TmpMap[idx] = Nx*Ny*Nz-1;
} }
} }
for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){ for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
int n = TmpMap[idx]; auto n = TmpMap[idx];
if (n > Nx*Ny*Nz){ if ( n > Nx*Ny*Nz ){
printf("Bad value! idx=%i \n"); printf("Bad value! idx=%i \n",n);
TmpMap[idx] = Nx*Ny*Nz-1; TmpMap[idx] = Nx*Ny*Nz-1;
} }
} }
@ -553,8 +551,9 @@ void ScaLBL_ColorModel::Run(){
} }
if (color_db->keyExists( "residual_endpoint_threshold" )){ if (color_db->keyExists( "residual_endpoint_threshold" )){
RESIDUAL_ENDPOINT_THRESHOLD = color_db->getScalar<double>( "residual_endpoint_threshold" ); RESIDUAL_ENDPOINT_THRESHOLD = color_db->getScalar<double>( "residual_endpoint_threshold" );
} }
NULL_USE( RESIDUAL_ENDPOINT_THRESHOLD );
if (color_db->keyExists( "noise_threshold" )){ if (color_db->keyExists( "noise_threshold" )){
NOISE_THRESHOLD = color_db->getScalar<double>( "noise_threshold" ); NOISE_THRESHOLD = color_db->getScalar<double>( "noise_threshold" );
USE_BUMP_RATE = true; USE_BUMP_RATE = true;
@ -874,7 +873,7 @@ void ScaLBL_ColorModel::Run(){
WriteHeader=true; WriteHeader=true;
kr_log_file = fopen("relperm.csv","a"); kr_log_file = fopen("relperm.csv","a");
if (WriteHeader) if (WriteHeader)
fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water eff.perm.oil.connected eff.perm.water.connected eff.perm.oil.disconnected eff.perm.water.disconnected cap.pressure cap.pressure.connected pressure.drop Ca M\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,pAB,viscous_pressure_drop,Ca,Mobility); fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water eff.perm.oil.connected eff.perm.water.connected eff.perm.oil.disconnected eff.perm.water.disconnected cap.pressure cap.pressure.connected pressure.drop Ca M\n");
fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,kAeff_disconnected,kBeff_disconnected,pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility); fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,kAeff_disconnected,kBeff_disconnected,pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility);
fclose(kr_log_file); fclose(kr_log_file);
@ -937,7 +936,7 @@ void ScaLBL_ColorModel::Run(){
else if (USE_SEED){ else if (USE_SEED){
delta_volume = volA*Dm->Volume - initial_volume; delta_volume = volA*Dm->Volume - initial_volume;
CURRENT_MORPH_TIMESTEPS += analysis_interval; CURRENT_MORPH_TIMESTEPS += analysis_interval;
double massChange = SeedPhaseField(seed_water); //double massChange = SeedPhaseField(seed_water);
if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", seed_water, delta_volume, delta_volume_target); if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", seed_water, delta_volume, delta_volume_target);
} }
else if (USE_MORPHOPEN_OIL){ else if (USE_MORPHOPEN_OIL){
@ -1010,7 +1009,6 @@ void ScaLBL_ColorModel::Run(){
double ScaLBL_ColorModel::ImageInit(std::string Filename){ double ScaLBL_ColorModel::ImageInit(std::string Filename){
bool suppress = false;
if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str()); if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str());
Mask->Decomp(Filename); Mask->Decomp(Filename);
for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
@ -1080,10 +1078,9 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){
ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,Dm->Comm); ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,Dm->Comm);
MPI_Barrier(Dm->Comm); MPI_Barrier(Dm->Comm);
int count_oil=0; long long count_connected=0;
int count_connected=0; long long count_porespace=0;
int count_porespace=0; long long count_water=0;
int count_water=0;
for (int k=1; k<nz-1; k++){ for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){ for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){ for (int i=1; i<nx-1; i++){
@ -1311,8 +1308,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
// 1. Copy phase field to CPU // 1. Copy phase field to CPU
ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double)); ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
double count,count_global,volume_initial,volume_final,volume_connected; double count = 0.f;
count = 0.f;
for (int k=1; k<Nz-1; k++){ for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){ for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){ for (int i=1; i<Nx-1; i++){
@ -1320,7 +1316,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
} }
} }
} }
volume_initial = sumReduce( Dm->Comm, count); double volume_initial = sumReduce( Dm->Comm, count);
/* /*
sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank); sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank);
FILE *INPUT = fopen(LocalRankFilename,"wb"); FILE *INPUT = fopen(LocalRankFilename,"wb");
@ -1352,16 +1348,16 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
} }
} }
} }
volume_connected = sumReduce( Dm->Comm, count); double volume_connected = sumReduce( Dm->Comm, count);
second_biggest = sumReduce( Dm->Comm, second_biggest); second_biggest = sumReduce( Dm->Comm, second_biggest);
int reach_x, reach_y, reach_z; /*int reach_x, reach_y, reach_z;
for (int k=0; k<Nz; k++){ for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){ for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){ for (int i=0; i<Nx; i++){
} }
} }
} }*/
// 3. Generate a distance map to the largest object -> phase_distance // 3. Generate a distance map to the largest object -> phase_distance
CalcDist(phase_distance,phase_id,*Dm); CalcDist(phase_distance,phase_id,*Dm);
@ -1417,7 +1413,6 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
for (int k=0; k<Nz; k++){ for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){ for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){ for (int i=0; i<Nx; i++){
int n = k*Nx*Ny + j*Nx + i;
double d = phase_distance(i,j,k); double d = phase_distance(i,j,k);
if (Averages->SDs(i,j,k) > 0.f){ if (Averages->SDs(i,j,k) > 0.f){
if (d < 3.f){ if (d < 3.f){
@ -1441,7 +1436,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
} }
} }
} }
volume_final= sumReduce( Dm->Comm, count); double volume_final= sumReduce( Dm->Comm, count);
delta_volume = (volume_final-volume_initial); delta_volume = (volume_final-volume_initial);
if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", delta_volume/volume_initial); if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", delta_volume/volume_initial);

1
models/DFHModel.cpp
View File

@ -114,7 +114,6 @@ void ScaLBL_DFHModel::SetDomain(){
} }
void ScaLBL_DFHModel::ReadInput(){ void ScaLBL_DFHModel::ReadInput(){
size_t readID;
//....................................................................... //.......................................................................
if (rank == 0) printf("Read input media... \n"); if (rank == 0) printf("Read input media... \n");
//....................................................................... //.......................................................................

4
models/MRTModel.cpp
View File

@ -94,7 +94,6 @@ void ScaLBL_MRTModel::SetDomain(){
void ScaLBL_MRTModel::ReadInput(){ void ScaLBL_MRTModel::ReadInput(){
int rank=Dm->rank(); int rank=Dm->rank();
size_t readID;
//....................................................................... //.......................................................................
//....................................................................... //.......................................................................
Mask->ReadIDs(); Mask->ReadIDs();
@ -106,7 +105,6 @@ void ScaLBL_MRTModel::ReadInput(){
// Generate the signed distance map // Generate the signed distance map
// Initialize the domain and communication // Initialize the domain and communication
Array<char> id_solid(Nx,Ny,Nz); Array<char> id_solid(Nx,Ny,Nz);
int count = 0;
// Solve for the position of the solid phase // Solve for the position of the solid phase
for (int k=0;k<Nz;k++){ for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){ for (int j=0;j<Ny;j++){
@ -122,7 +120,6 @@ void ScaLBL_MRTModel::ReadInput(){
for (int k=0;k<Nz;k++){ for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){ for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){ for (int i=0;i<Nx;i++){
int n=k*Nx*Ny+j*Nx+i;
// Initialize distance to +/- 1 // Initialize distance to +/- 1
Distance(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0; Distance(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
} }
@ -191,7 +188,6 @@ void ScaLBL_MRTModel::Run(){
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
Minkowski Morphology(Mask); Minkowski Morphology(Mask);
int SIZE=Np*sizeof(double);
if (rank==0){ if (rank==0){
bool WriteHeader=false; bool WriteHeader=false;

5
tests/GenerateSphereTest.cpp
View File

@ -363,7 +363,7 @@ int main(int argc, char **argv)
nspheres = domain_db->getScalar<int>( "nspheres"); nspheres = domain_db->getScalar<int>( "nspheres");
//printf("Set domain \n"); //printf("Set domain \n");
int BoundaryCondition=1; //int BoundaryCondition=1;
//Nz += 2; //Nz += 2;
//Nx = Ny = Nz; // Cubic domain //Nx = Ny = Nz; // Cubic domain
int N = Nx*Ny*Nz; int N = Nx*Ny*Nz;
@ -396,7 +396,7 @@ int main(int argc, char **argv)
int sum = 0; int sum = 0;
double sum_local; double sum_local;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs); double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
double porosity, pore_vol; double porosity;
//........................................................................... //...........................................................................
DoubleArray SignDist(Nx,Ny,Nz); DoubleArray SignDist(Nx,Ny,Nz);
//....................................................................... //.......................................................................
@ -450,7 +450,6 @@ int main(int argc, char **argv)
} }
} }
sum=0; sum=0;
pore_vol = 0.0;
for ( k=1;k<Nz-1;k++){ for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){ for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){ for ( i=1;i<Nx-1;i++){

4
tests/TestBubbleDFH.cpp
View File

@ -193,7 +193,7 @@ int main(int argc, char **argv)
// char value; // char value;
char *id; char *id;
id = new char[N]; id = new char[N];
double sum, sum_local; double sum;
//........................................................................... //...........................................................................
if (rank == 0) cout << "Setting up bubble..." << endl; if (rank == 0) cout << "Setting up bubble..." << endl;
double BubbleRadius = 15.5; // Radius of the capillary tube double BubbleRadius = 15.5; // Radius of the capillary tube
@ -516,7 +516,7 @@ int main(int argc, char **argv)
DoubleArray PhaseField(Nx,Ny,Nz); DoubleArray PhaseField(Nx,Ny,Nz);
ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField); ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField);
FILE *OUTFILE; FILE *OUTFILE;
sprintf(LocalRankFilename,"Phase.raw",rank); sprintf(LocalRankFilename,"Phase.raw");
OUTFILE = fopen(LocalRankFilename,"wb"); OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE); fwrite(PhaseField.data(),8,N,OUTFILE);
fclose(OUTFILE); fclose(OUTFILE);

6
tests/TestColorGradDFH.cpp
View File

@ -53,9 +53,6 @@ int main(int argc, char **argv)
int Nx = db->getVector<int>( "n" )[0]; int Nx = db->getVector<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1]; int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2]; int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0];
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "nproc" )[2];
if (rank==0){ if (rank==0){
printf("********************************************************\n"); printf("********************************************************\n");
@ -64,7 +61,7 @@ int main(int argc, char **argv)
} }
// Get the rank info // Get the rank info
std::shared_ptr<Domain> Dm(new Domain(db,comm)); auto Dm = std::make_shared<Domain>(db,comm);
Nx += 2; Nx += 2;
Ny += 2; Ny += 2;
Nz += 2; Nz += 2;
@ -111,7 +108,6 @@ int main(int argc, char **argv)
MPI_Barrier(comm); MPI_Barrier(comm);
//......................device distributions................................. //......................device distributions.................................
int dist_mem_size = Np*sizeof(double);
int neighborSize=18*Np*sizeof(int); int neighborSize=18*Np*sizeof(int);
if (rank==0) printf ("Allocating distributions \n"); if (rank==0) printf ("Allocating distributions \n");
int *NeighborList; int *NeighborList;

27
tests/TestCommD3Q19.cpp
View File

@ -49,7 +49,7 @@ extern void GlobalFlipScaLBL_D3Q19_Init(double *dist, IntArray Map, int Np, int
{1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0},{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1}, {1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0},{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1},
{0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}}; {0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}};
int q,i,j,k,n,N; int q,i,j,k;
int Cqx,Cqy,Cqz; // Discrete velocity int Cqx,Cqy,Cqz; // Discrete velocity
int x,y,z; // Global indices int x,y,z; // Global indices
int xn,yn,zn; // Global indices of neighbor int xn,yn,zn; // Global indices of neighbor
@ -59,8 +59,6 @@ extern void GlobalFlipScaLBL_D3Q19_Init(double *dist, IntArray Map, int Np, int
Y = Ny*nprocy; Y = Ny*nprocy;
Z = Nz*nprocz; Z = Nz*nprocz;
NULL_USE(Z); NULL_USE(Z);
N = (Nx+2)*(Ny+2)*(Nz+2); // size of the array including halo
for (k=0; k<Nz; k++){ for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){ for (j=0; j<Ny; j++){
@ -104,16 +102,13 @@ extern int GlobalCheckDebugDist(double *dist, IntArray Map, int Np, int Nx, int
{ {
int returnValue = 0; int returnValue = 0;
int q,i,j,k,n,N,idx; int q,i,j,k,idx;
int Cqx,Cqy,Cqz; // Discrete velocity
int x,y,z; // Global indices int x,y,z; // Global indices
int xn,yn,zn; // Global indices of neighbor
int X,Y,Z; // Global size int X,Y,Z; // Global size
X = Nx*nprocx; X = Nx*nprocx;
Y = Ny*nprocy; Y = Ny*nprocy;
Z = Nz*nprocz; Z = Nz*nprocz;
NULL_USE(Z); NULL_USE(Z);
N = (Nx+2)*(Ny+2)*(Nz+2); // size of the array including halo
for (k=0; k<Nz; k++){ for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){ for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){ for (i=0; i<Nx; i++){
@ -168,9 +163,6 @@ inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
//*************************************************************************************** //***************************************************************************************
int main(int argc, char **argv) int main(int argc, char **argv)
{ {
//*****************************************
// ***** MPI STUFF ****************
//*****************************************
// Initialize MPI // Initialize MPI
int rank,nprocs; int rank,nprocs;
MPI_Init(&argc,&argv); MPI_Init(&argc,&argv);
@ -178,10 +170,7 @@ int main(int argc, char **argv)
MPI_Comm_rank(comm,&rank); MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs); MPI_Comm_size(comm,&nprocs);
int check; int check;
{ {
MPI_Request req1[18],req2[18];
MPI_Status stat1[18],stat2[18];
if (rank == 0){ if (rank == 0){
printf("********************************************************\n"); printf("********************************************************\n");
@ -191,11 +180,8 @@ int main(int argc, char **argv)
// BGK Model parameters // BGK Model parameters
string FILENAME; string FILENAME;
unsigned int nBlocks, nthreads;
int timestepMax, interval;
double tau,Fx,Fy,Fz,tol;
// Domain variables // Domain variables
int i,j,k,n; int i,j,k;
// Load inputs // Load inputs
auto db = loadInputs( nprocs ); auto db = loadInputs( nprocs );
@ -223,8 +209,7 @@ int main(int argc, char **argv)
char LocalRankFilename[40]; char LocalRankFilename[40];
sprintf(LocalRankFilename,"ID.%05i",rank); sprintf(LocalRankFilename,"ID.%05i",rank);
char *id; auto id = new char[Nx*Ny*Nz];
id = new char[Nx*Ny*Nz];
/* if (rank==0) printf("Assigning phase ID from file \n"); /* if (rank==0) printf("Assigning phase ID from file \n");
if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n"); if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
@ -237,7 +222,7 @@ int main(int argc, char **argv)
for (k=0;k<Nz;k++){ for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){ for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){ for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i; int n = k*Nx*Ny+j*Nx+i;
id[n] = 1; id[n] = 1;
Dm->id[n] = id[n]; Dm->id[n] = id[n];
} }
@ -270,7 +255,7 @@ int main(int argc, char **argv)
for (k=1;k<Nz-1;k++){ for (k=1;k<Nz-1;k++){
for (j=1;j<Ny-1;j++){ for (j=1;j<Ny-1;j++){
for (i=1;i<Nx-1;i++){ for (i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i; int n = k*Nx*Ny+j*Nx+i;
if (id[n] == component){ if (id[n] == component){
sum_local+=1.0; sum_local+=1.0;
} }

20
tests/TestFluxBC.cpp
View File

@ -32,21 +32,14 @@ int main (int argc, char **argv)
} }
{ {
int i,j,k,n,Np; int i,j,k,n,Np;
bool pBC=true;
double Lx,Ly,Lz;
Lx = Ly = Lz = 1.f;
double din,dout; double din,dout;
int BC=1;
// Load inputs // Load inputs
auto db = loadInputs( nprocs ); auto db = loadInputs( nprocs );
int Nx = db->getVector<int>( "n" )[0]; int Nx = db->getVector<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1]; int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2]; int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0]; auto Dm = std::make_shared<Domain>(db,comm);
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "nproc" )[2];
std::shared_ptr<Domain> Dm(new Domain(db,comm));
Nx += 2; Ny+=2; Nz += 2; Nx += 2; Ny+=2; Nz += 2;
Nx = Ny = Nz; // Cubic domain Nx = Ny = Nz; // Cubic domain
@ -55,8 +48,7 @@ int main (int argc, char **argv)
//....................................................................... //.......................................................................
// Assign the phase ID // Assign the phase ID
//....................................................................... //.......................................................................
char *id; auto id = new char[N];
id = new char[N];
for (k=0;k<Nz;k++){ for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){ for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){ for (i=0;i<Nx;i++){
@ -160,9 +152,7 @@ int main (int argc, char **argv)
ScaLBL_DeviceBarrier(); MPI_Barrier(comm); ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_CopyToHost(&VEL[0],&dvc_vel[0],SIZE); ScaLBL_CopyToHost(&VEL[0],&dvc_vel[0],SIZE);
double err,value,Q; double Q = 0.f;
Q = 0.f;
k=1; k=1;
for (j=1;j<Ny-1;j++){ for (j=1;j<Ny-1;j++){
for (i=1;i<Nx-1;i++){ for (i=1;i<Nx-1;i++){
@ -176,7 +166,7 @@ int main (int argc, char **argv)
// respect backwards read / write!!! // respect backwards read / write!!!
printf("Inlet Flux: input=%f, output=%f \n",flux,Q); printf("Inlet Flux: input=%f, output=%f \n",flux,Q);
err = fabs(flux + Q); double err = fabs(flux + Q);
if (err > 1e-12){ if (err > 1e-12){
error = 1; error = 1;
printf(" Inlet error %f \n",err); printf(" Inlet error %f \n",err);
@ -185,7 +175,7 @@ int main (int argc, char **argv)
// Consider a larger number of timesteps and simulate flow // Consider a larger number of timesteps and simulate flow
double Fx, Fy, Fz; double Fx, Fy, Fz;
double tau = 1.0; double tau = 1.0;
double mu=(tau-0.5)/3.0; //double mu=(tau-0.5)/3.0;
double rlx_setA=1.0/tau; double rlx_setA=1.0/tau;
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
dout=1.f; dout=1.f;

14
tests/TestForceD3Q19.cpp
View File

@ -457,24 +457,16 @@ int main (int argc, char **argv)
double *x = new double[1]; double *x = new double[1];
ASSERT(x!=NULL); ASSERT(x!=NULL);
} }
// set the error code
// Note: the error code should be consistent across all processors
int error = 0;
int Np = 1; int Np = 1;
int Q = 9; //int Q = 9;
double Fx = 1.0; double Fx = 1.0;
double Fy = 1.0; double Fy = 1.0;
double Fz = 1.0; double Fz = 1.0;
double *dist; auto dist = new double [19*Np];
double * Velocity; //auto Velocity = new double [3*Np
dist = new double [19*Np];
Velocity = new double [3*Np];
for (int n=0; n<Np; n++){ for (int n=0; n<Np; n++){
dist[n] = 0.3333333333333333; dist[n] = 0.3333333333333333;

6
tests/TestInterfaceSpeed.cpp
View File

@ -24,7 +24,7 @@ int main (int argc, char *argv[])
MPI_Comm_rank(comm,&rank); MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs); MPI_Comm_size(comm,&nprocs);
int i,j,k,n; int i,j,k;
// Load inputs // Load inputs
string FILENAME = argv[1]; string FILENAME = argv[1];
@ -36,7 +36,7 @@ int main (int argc, char *argv[])
int Ny = domain_db->getVector<int>( "n" )[1]; int Ny = domain_db->getVector<int>( "n" )[1];
int Nz = domain_db->getVector<int>( "n" )[2]; int Nz = domain_db->getVector<int>( "n" )[2];
std::shared_ptr<Domain> Dm(new Domain(domain_db,comm)); auto Dm = std::make_shared<Domain>(domain_db,comm);
Nx+=2; Ny+=2; Nz+=2; Nx+=2; Ny+=2; Nz+=2;
@ -44,7 +44,7 @@ int main (int argc, char *argv[])
Dm->CommInit(); Dm->CommInit();
std::shared_ptr<TwoPhase> Averages(new TwoPhase(Dm)); auto Averages = std::make_shared<TwoPhase>(Dm);
int timestep=0; int timestep=0;
double Cx,Cy,Cz; double Cx,Cy,Cz;

6
tests/TestMap.cpp
View File

@ -56,11 +56,7 @@ int main(int argc, char **argv)
int Nx = db->getVector<int>( "n" )[0]; int Nx = db->getVector<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1]; int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2]; int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0]; auto Dm = std::make_shared<Domain>(db,comm);
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "nproc" )[2];
std::shared_ptr<Domain> Dm(new Domain(db,comm));
Nx += 2; Nx += 2;
Ny += 2; Ny += 2;

25
tests/TestMassConservationD3Q7.cpp
View File

@ -66,9 +66,6 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
int main(int argc, char **argv) int main(int argc, char **argv)
{ {
//*****************************************
// ***** MPI STUFF ****************
//*****************************************
// Initialize MPI // Initialize MPI
int rank,nprocs; int rank,nprocs;
MPI_Init(&argc,&argv); MPI_Init(&argc,&argv);
@ -76,19 +73,6 @@ int main(int argc, char **argv)
MPI_Comm_rank(comm,&rank); MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs); MPI_Comm_size(comm,&nprocs);
// parallel domain size (# of sub-domains) // parallel domain size (# of sub-domains)
int nprocx,nprocy,nprocz;
int iproc,jproc,kproc;
int sendtag,recvtag;
//*****************************************
// MPI ranks for all 18 neighbors
//**********************************
int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z;
int rank_xy,rank_XY,rank_xY,rank_Xy;
int rank_xz,rank_XZ,rank_xZ,rank_Xz;
int rank_yz,rank_YZ,rank_yZ,rank_Yz;
//**********************************
MPI_Request req1[18],req2[18];
MPI_Status stat1[18],stat2[18];
if (rank == 0){ if (rank == 0){
printf("********************************************************\n"); printf("********************************************************\n");
@ -110,7 +94,6 @@ int main(int argc, char **argv)
Ny = CM.Ny; Ny = CM.Ny;
Nz = CM.Nz; Nz = CM.Nz;
N = Nx*Ny*Nz; N = Nx*Ny*Nz;
int dist_mem_size = N*sizeof(double);
//CM.ReadInput(); //CM.ReadInput();
double radius=0.4*double(Nx); double radius=0.4*double(Nx);
@ -142,11 +125,9 @@ int main(int argc, char **argv)
CM.Run(); CM.Run();
int D3Q7[7][3]={{0,0,0},{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}}; int D3Q7[7][3]={{0,0,0},{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};
// Compare and make sure mass is conserved at every lattice site // Compare and make sure mass is conserved at every lattice site
double *Error; auto Error = new double[N];
Error = new double [N]; auto A_q = new double[7*Np];
double *A_q, *B_q; //auto B_q = new double[7*Np];
A_q = new double [7*Np];
B_q = new double [7*Np];
bool CleanCheck = true; bool CleanCheck = true;
double original,final, sum_q; double original,final, sum_q;
double total_mass_A_0 = 0.0; double total_mass_A_0 = 0.0;

1
tests/TestNetcdf.cpp
View File

@ -14,7 +14,6 @@ void load( const std::string& );
void test_NETCDF( UnitTest& ut ) void test_NETCDF( UnitTest& ut )
{ {
const int rank = comm_rank( MPI_COMM_WORLD ); const int rank = comm_rank( MPI_COMM_WORLD );
const int size = comm_size( MPI_COMM_WORLD );
int nprocx = 2; int nprocx = 2;
int nprocy = 2; int nprocy = 2;
int nprocz = 2; int nprocz = 2;

6
tests/TestSubphase.cpp
View File

@ -60,13 +60,11 @@ int main(int argc, char **argv)
} }
// Get the rank info // Get the rank info
std::shared_ptr<Domain> Dm(new Domain(db,comm)); auto Dm = std::make_shared<Domain>(db,comm);
// const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); auto Averages = std::make_shared<SubPhase>(Dm);
std::shared_ptr<SubPhase> Averages(new SubPhase(Dm));
Nx += 2; Nx += 2;
Ny += 2; Ny += 2;
Nz += 2; Nz += 2;
int N = Nx*Ny*Nz;
//....................................................................... //.......................................................................
for ( k=1;k<Nz-1;k++){ for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){ for ( j=1;j<Ny-1;j++){

3
tests/TestTopo3D.cpp
View File

@ -60,12 +60,11 @@ int main(int argc, char **argv)
} }
// Get the rank info // Get the rank info
std::shared_ptr<Domain> Dm(new Domain(db,comm)); auto Dm = std::make_shared<Domain>(db,comm);
Nx += 2; Nx += 2;
Ny += 2; Ny += 2;
Nz += 2; Nz += 2;
int N = Nx*Ny*Nz;
//....................................................................... //.......................................................................
for ( k=1;k<Nz-1;k++){ for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){ for ( j=1;j<Ny-1;j++){

7
tests/TestTorus.cpp
View File

@ -60,13 +60,12 @@ int main(int argc, char **argv)
} }
// Get the rank info // Get the rank info
std::shared_ptr<Domain> Dm(new Domain(db,comm)); auto Dm = std::make_shared<Domain>(db,comm);
// const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); // const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
std::shared_ptr<TwoPhase> Averages(new TwoPhase(Dm)); auto Averages = std::make_shared<TwoPhase>(Dm);
Nx += 2; Nx += 2;
Ny += 2; Ny += 2;
Nz += 2; Nz += 2;
int N = Nx*Ny*Nz;
//....................................................................... //.......................................................................
for ( k=1;k<Nz-1;k++){ for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){ for ( j=1;j<Ny-1;j++){
@ -142,7 +141,7 @@ int main(int argc, char **argv)
} }
} }
double beta = 0.95; //double beta = 0.95;
if (rank==0) printf("initializing the system \n"); if (rank==0) printf("initializing the system \n");
Averages->UpdateSolid(); Averages->UpdateSolid();

10
tests/TestTorusEvolve.cpp
View File

@ -60,12 +60,11 @@ int main(int argc, char **argv)
} }
// Get the rank info // Get the rank info
std::shared_ptr<Domain> Dm(new Domain(db,comm)); auto Dm = std::make_shared<Domain>(db,comm);
Nx += 2; Nx += 2;
Ny += 2; Ny += 2;
Nz += 2; Nz += 2;
int N = Nx*Ny*Nz;
//....................................................................... //.......................................................................
for ( k=1;k<Nz-1;k++){ for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){ for ( j=1;j<Ny-1;j++){
@ -98,14 +97,13 @@ int main(int argc, char **argv)
//....................................................................... //.......................................................................
// Assign the phase ID field based and the signed distance // Assign the phase ID field based and the signed distance
//....................................................................... //.......................................................................
double R1,R2,R;
double CX,CY,CZ; //CY1,CY2; double CX,CY,CZ; //CY1,CY2;
CX=Nx*nprocx*0.5; CX=Nx*nprocx*0.5;
CY=Ny*nprocy*0.5; CY=Ny*nprocy*0.5;
CZ=Nz*nprocz*0.5; CZ=Nz*nprocz*0.5;
R1 = (Nx-2)*nprocx*0.3; // middle radius auto R1 = (Nx-2)*nprocx*0.3; // middle radius
R2 = (Nx-2)*nprocx*0.1; // donut thickness auto R2 = (Nx-2)*nprocx*0.1; // donut thickness
R = 0.4*nprocx*(Nx-2); //auto R = 0.4*nprocx*(Nx-2);
Minkowski Object(Dm); Minkowski Object(Dm);

8
tests/lbpm_morph_pp.cpp
View File

@ -34,7 +34,6 @@ int main(int argc, char **argv)
//....................................................................... //.......................................................................
int n, nx, ny, nz; int n, nx, ny, nz;
char LocalRankFilename[40]; char LocalRankFilename[40];
char FILENAME[128];
string filename; string filename;
double SW; double SW;
@ -239,10 +238,9 @@ int main(int argc, char **argv)
} }
MPI_Barrier(comm); MPI_Barrier(comm);
sprintf(FILENAME,READFILE.c_str()); auto filename2 = READFILE + ".morph.raw";
sprintf(FILENAME+strlen(FILENAME),".morph.raw"); if (rank==0) printf("Writing file to: %s \n", filename2.c_str());
if (rank==0) printf("Writing file to: %s \n", FILENAME); Mask->AggregateLabels(filename2);
Mask->AggregateLabels(FILENAME);
} }
MPI_Barrier(comm); MPI_Barrier(comm);

31
tests/lbpm_morphdrain_pp.cpp
View File

@ -32,10 +32,7 @@ int main(int argc, char **argv)
//....................................................................... //.......................................................................
// Reading the domain information file // Reading the domain information file
//....................................................................... //.......................................................................
int n, nprocx, nprocy, nprocz, nx, ny, nz;
char LocalRankString[8];
char LocalRankFilename[40]; char LocalRankFilename[40];
char FILENAME[128];
string filename; string filename;
double SW,Rcrit_new; double SW,Rcrit_new;
@ -43,8 +40,10 @@ int main(int argc, char **argv)
filename=argv[1]; filename=argv[1];
Rcrit_new=0.f; Rcrit_new=0.f;
//SW=strtod(argv[2],NULL); //SW=strtod(argv[2],NULL);
} } else {
else ERROR("No input database provided\n"); ERROR("No input database provided\n");
}
NULL_USE( Rcrit_new );
// read the input database // read the input database
auto db = std::make_shared<Database>( filename ); auto db = std::make_shared<Database>( filename );
auto domain_db = db->getDatabase( "Domain" ); auto domain_db = db->getDatabase( "Domain" );
@ -62,19 +61,16 @@ int main(int argc, char **argv)
if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW); if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW);
// GenerateResidual(id,nx,ny,nz,Saturation); // GenerateResidual(id,nx,ny,nz,Saturation);
nx = size[0]; int nx = size[0];
ny = size[1]; int ny = size[1];
nz = size[2]; int nz = size[2];
nprocx = nproc[0];
nprocy = nproc[1];
nprocz = nproc[2];
int N = (nx+2)*(ny+2)*(nz+2); size_t N = (nx+2)*(ny+2)*(nz+2);
std::shared_ptr<Domain> Dm (new Domain(domain_db,comm)); std::shared_ptr<Domain> Dm (new Domain(domain_db,comm));
std::shared_ptr<Domain> Mask (new Domain(domain_db,comm)); std::shared_ptr<Domain> Mask (new Domain(domain_db,comm));
// std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC)); // std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
for (n=0; n<N; n++) Dm->id[n]=1; for (size_t n=0; n<N; n++) Dm->id[n]=1;
Dm->CommInit(); Dm->CommInit();
signed char *id; signed char *id;
@ -116,7 +112,6 @@ int main(int argc, char **argv)
for (int k=0;k<nz;k++){ for (int k=0;k<nz;k++){
for (int j=0;j<ny;j++){ for (int j=0;j<ny;j++){
for (int i=0;i<nx;i++){ for (int i=0;i<nx;i++){
int n = k*nx*ny+j*nx+i;
// Initialize distance to +/- 1 // Initialize distance to +/- 1
SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0; SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
} }
@ -158,7 +153,6 @@ int main(int argc, char **argv)
for (int k=0;k<nz;k++){ for (int k=0;k<nz;k++){
for (int j=0;j<ny;j++){ for (int j=0;j<ny;j++){
for (int i=0;i<nx;i++){ for (int i=0;i<nx;i++){
int n = k*nx*ny+j*nx+i;
// Initialize distance to +/- 1 // Initialize distance to +/- 1
SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0; SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
} }
@ -204,10 +198,9 @@ int main(int argc, char **argv)
} }
MPI_Barrier(comm); MPI_Barrier(comm);
sprintf(FILENAME,READFILE.c_str()); auto filename2 = READFILE + ".morphdrain.raw";
sprintf(FILENAME+strlen(FILENAME),".morphdrain.raw"); if (rank==0) printf("Writing file to: %s \n", filename2.data() );
if (rank==0) printf("Writing file to: %s \n", FILENAME); Mask->AggregateLabels( filename2 );
Mask->AggregateLabels(FILENAME);
} }
MPI_Barrier(comm); MPI_Barrier(comm);

26
tests/lbpm_morphopen_pp.cpp
View File

@ -32,10 +32,7 @@ int main(int argc, char **argv)
//....................................................................... //.......................................................................
// Reading the domain information file // Reading the domain information file
//....................................................................... //.......................................................................
int n, nprocx, nprocy, nprocz, nx, ny, nz;
char LocalRankString[8];
char LocalRankFilename[40]; char LocalRankFilename[40];
char FILENAME[128];
string filename; string filename;
double SW,Rcrit_new; double SW,Rcrit_new;
@ -45,6 +42,7 @@ int main(int argc, char **argv)
//SW=strtod(argv[2],NULL); //SW=strtod(argv[2],NULL);
} }
else ERROR("No input database provided\n"); else ERROR("No input database provided\n");
NULL_USE( Rcrit_new );
// read the input database // read the input database
auto db = std::make_shared<Database>( filename ); auto db = std::make_shared<Database>( filename );
auto domain_db = db->getDatabase( "Domain" ); auto domain_db = db->getDatabase( "Domain" );
@ -69,19 +67,16 @@ int main(int argc, char **argv)
if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW); if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW);
// GenerateResidual(id,nx,ny,nz,Saturation); // GenerateResidual(id,nx,ny,nz,Saturation);
nx = size[0]; int nx = size[0];
ny = size[1]; int ny = size[1];
nz = size[2]; int nz = size[2];
nprocx = nproc[0];
nprocy = nproc[1];
nprocz = nproc[2];
int N = (nx+2)*(ny+2)*(nz+2); size_t N = (nx+2)*(ny+2)*(nz+2);
std::shared_ptr<Domain> Dm (new Domain(domain_db,comm)); std::shared_ptr<Domain> Dm (new Domain(domain_db,comm));
std::shared_ptr<Domain> Mask (new Domain(domain_db,comm)); std::shared_ptr<Domain> Mask (new Domain(domain_db,comm));
// std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC)); // std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
for (n=0; n<N; n++) Dm->id[n]=1; for (size_t n=0; n<N; n++) Dm->id[n]=1;
Dm->CommInit(); Dm->CommInit();
signed char *id; signed char *id;
@ -119,7 +114,6 @@ int main(int argc, char **argv)
for (int k=0;k<nz;k++){ for (int k=0;k<nz;k++){
for (int j=0;j<ny;j++){ for (int j=0;j<ny;j++){
for (int i=0;i<nx;i++){ for (int i=0;i<nx;i++){
int n = k*nx*ny+j*nx+i;
// Initialize distance to +/- 1 // Initialize distance to +/- 1
SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0; SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
} }
@ -161,7 +155,6 @@ int main(int argc, char **argv)
for (int k=0;k<nz;k++){ for (int k=0;k<nz;k++){
for (int j=0;j<ny;j++){ for (int j=0;j<ny;j++){
for (int i=0;i<nx;i++){ for (int i=0;i<nx;i++){
int n = k*nx*ny+j*nx+i;
// Initialize distance to +/- 1 // Initialize distance to +/- 1
SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0; SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
} }
@ -207,10 +200,9 @@ int main(int argc, char **argv)
} }
MPI_Barrier(comm); MPI_Barrier(comm);
sprintf(FILENAME,READFILE.c_str()); auto filename2 = READFILE + ".morphopen.raw";
sprintf(FILENAME+strlen(FILENAME),".morphopen.raw"); if (rank==0) printf("Writing file to: %s \n", filename2.data());
if (rank==0) printf("Writing file to: %s \n", FILENAME); Mask->AggregateLabels(filename2);
Mask->AggregateLabels(FILENAME);
} }
MPI_Barrier(comm); MPI_Barrier(comm);

9
tests/lbpm_permeability_simulator.cpp
View File

@ -23,9 +23,6 @@ using namespace std;
int main(int argc, char **argv) int main(int argc, char **argv)
{ {
//*****************************************
// ***** MPI STUFF ****************
//*****************************************
// Initialize MPI // Initialize MPI
int rank,nprocs; int rank,nprocs;
MPI_Init(&argc,&argv); MPI_Init(&argc,&argv);
@ -33,10 +30,6 @@ int main(int argc, char **argv)
MPI_Comm_rank(comm,&rank); MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs); MPI_Comm_size(comm,&nprocs);
{ {
// parallel domain size (# of sub-domains)
int nprocx,nprocy,nprocz;
int iproc,jproc,kproc;
if (rank == 0){ if (rank == 0){
printf("********************************************************\n"); printf("********************************************************\n");
printf("Running Single Phase Permeability Calculation \n"); printf("Running Single Phase Permeability Calculation \n");
@ -44,10 +37,10 @@ int main(int argc, char **argv)
} }
// Initialize compute device // Initialize compute device
int device=ScaLBL_SetDevice(rank); int device=ScaLBL_SetDevice(rank);
NULL_USE( device );
ScaLBL_DeviceBarrier(); ScaLBL_DeviceBarrier();
MPI_Barrier(comm); MPI_Barrier(comm);
ScaLBL_MRTModel MRT(rank,nprocs,comm); ScaLBL_MRTModel MRT(rank,nprocs,comm);
auto filename = argv[1]; auto filename = argv[1];
MRT.ReadParams(filename); MRT.ReadParams(filename);

29
tests/lbpm_refine_pp.cpp
View File

@ -26,10 +26,9 @@ int main(int argc, char **argv)
//....................................................................... //.......................................................................
// Reading the domain information file // Reading the domain information file
//....................................................................... //.......................................................................
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz; double Lx, Ly, Lz;
Lx = Ly = Lz = 1.0;
int i,j,k,n; int i,j,k,n;
int BC=0;
string filename; string filename;
if (argc > 1){ if (argc > 1){
@ -47,12 +46,12 @@ int main(int argc, char **argv)
auto ReadValues = domain_db->getVector<char>( "ReadValues" ); auto ReadValues = domain_db->getVector<char>( "ReadValues" );
auto WriteValues = domain_db->getVector<char>( "WriteValues" ); auto WriteValues = domain_db->getVector<char>( "WriteValues" );
nx = size[0]; int nx = size[0];
ny = size[1]; int ny = size[1];
nz = size[2]; int nz = size[2];
nprocx = nproc[0]; int nprocx = nproc[0];
nprocy = nproc[1]; int nprocy = nproc[1];
nprocz = nproc[2]; int nprocz = nproc[2];
// Check that the number of processors >= the number of ranks // Check that the number of processors >= the number of ranks
if ( rank==0 ) { if ( rank==0 ) {
@ -66,10 +65,9 @@ int main(int argc, char **argv)
char LocalRankFilename[40]; char LocalRankFilename[40];
int rnx,rny,rnz; int rnx=2*nx;
rnx=2*nx; int rny=2*ny;
rny=2*ny; int rnz=2*nz;
rnz=2*nz;
if (rank==0) printf("Refining mesh to %i x %i x %i \n",rnx,rny,rnz); if (rank==0) printf("Refining mesh to %i x %i x %i \n",rnx,rny,rnz);
@ -128,13 +126,12 @@ int main(int argc, char **argv)
} }
} }
int ri,rj,rk,rn; //refined mesh indices
//char *RefineLabel; //char *RefineLabel;
//RefineLabel = new char [rnx*rny*rnz]; //RefineLabel = new char [rnx*rny*rnz];
Array <char> RefineLabel(rnx,rny,rnz); Array <char> RefineLabel(rnx,rny,rnz);
for (rk=1; rk<rnz-1; rk++){ for (int rk=1; rk<rnz-1; rk++){
for (rj=1; rj<rny-1; rj++){ for (int rj=1; rj<rny-1; rj++){
for (ri=1; ri<rnx-1; ri++){ for (int ri=1; ri<rnx-1; ri++){
n = rk*rnx*rny+rj*rnx+ri; n = rk*rnx*rny+rj*rnx+ri;
// starting node for each processor matches exactly // starting node for each processor matches exactly
i = (ri-1)/2+1; i = (ri-1)/2+1;

39
tests/lbpm_serial_decomp.cpp
View File

@ -47,11 +47,7 @@ int main(int argc, char **argv)
//....................................................................... //.......................................................................
// Reading the domain information file // Reading the domain information file
//....................................................................... //.......................................................................
int nprocs, nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
int64_t Nx,Ny,Nz;
int64_t i,j,k,n; int64_t i,j,k,n;
int BC=0;
int64_t xStart,yStart,zStart; int64_t xStart,yStart,zStart;
int checkerSize; int checkerSize;
int inlet_count_x, inlet_count_y, inlet_count_z; int inlet_count_x, inlet_count_y, inlet_count_z;
@ -112,25 +108,25 @@ int main(int argc, char **argv)
ReadType = "8bit"; ReadType = "8bit";
} }
nx = size[0]; int nx = size[0];
ny = size[1]; int ny = size[1];
nz = size[2]; int nz = size[2];
nprocx = nproc[0]; int nprocx = nproc[0];
nprocy = nproc[1]; int nprocy = nproc[1];
nprocz = nproc[2]; int nprocz = nproc[2];
Nx = SIZE[0]; long int Nx = SIZE[0];
Ny = SIZE[1]; long int Ny = SIZE[1];
Nz = SIZE[2]; long int Nz = SIZE[2];
printf("Input media: %s\n",Filename.c_str()); printf("Input media: %s\n",Filename.c_str());
printf("Relabeling %lu values\n",ReadValues.size()); printf("Relabeling %lu values\n",ReadValues.size());
for (int idx=0; idx<ReadValues.size(); idx++){ for (size_t idx=0; idx<ReadValues.size(); idx++){
int oldvalue=ReadValues[idx]; int oldvalue=ReadValues[idx];
int newvalue=WriteValues[idx]; int newvalue=WriteValues[idx];
printf("oldvalue=%d, newvalue =%d \n",oldvalue,newvalue); printf("oldvalue=%d, newvalue =%d \n",oldvalue,newvalue);
} }
nprocs=nprocx*nprocy*nprocz; int nprocs=nprocx*nprocy*nprocz;
char *SegData = NULL; char *SegData = NULL;
// Rank=0 reads the entire segmented data and distributes to worker processes // Rank=0 reads the entire segmented data and distributes to worker processes
@ -172,7 +168,7 @@ int main(int argc, char **argv)
n = k*Nx*Ny+j*Nx+i; n = k*Nx*Ny+j*Nx+i;
//char locval = loc_id[n]; //char locval = loc_id[n];
char locval = SegData[n]; char locval = SegData[n];
for (int idx=0; idx<ReadValues.size(); idx++){ for (size_t idx=0; idx<ReadValues.size(); idx++){
signed char oldvalue=ReadValues[idx]; signed char oldvalue=ReadValues[idx];
signed char newvalue=WriteValues[idx]; signed char newvalue=WriteValues[idx];
if (locval == oldvalue){ if (locval == oldvalue){
@ -185,10 +181,10 @@ int main(int argc, char **argv)
} }
} }
if (rank==0){ if (rank==0){
for (int idx=0; idx<ReadValues.size(); idx++){ for (size_t idx=0; idx<ReadValues.size(); idx++){
long int label=ReadValues[idx]; long int label=ReadValues[idx];
long int count=LabelCount[idx]; long int count=LabelCount[idx];
printf("Label=%d, Count=%d \n",label,count); printf("Label=%ld, Count=%ld \n",label,count);
} }
} }
@ -215,7 +211,7 @@ int main(int argc, char **argv)
printf("Checkerboard pattern at y inlet for %i layers \n",inlet_count_y); printf("Checkerboard pattern at y inlet for %i layers \n",inlet_count_y);
// use checkerboard pattern // use checkerboard pattern
for (int k = 0; k<Nz; k++){ for (int k = 0; k<Nz; k++){
for (int j = yStart; i < yStart+inlet_count_y; j++){ for (int j = yStart; j < yStart+inlet_count_y; j++){
for (int i = 0; i<Nx; i++){ for (int i = 0; i<Nx; i++){
if ( (i/checkerSize + k/checkerSize)%2 == 0){ if ( (i/checkerSize + k/checkerSize)%2 == 0){
// void checkers // void checkers
@ -272,7 +268,7 @@ int main(int argc, char **argv)
printf("Checkerboard pattern at y outlet for %i layers \n",outlet_count_y); printf("Checkerboard pattern at y outlet for %i layers \n",outlet_count_y);
// use checkerboard pattern // use checkerboard pattern
for (int k = 0; k<Nz; k++){ for (int k = 0; k<Nz; k++){
for (int j = yStart + ny*nprocy - outlet_count_y; i < yStart + ny*nprocy; j++){ for (int j = yStart + ny*nprocy - outlet_count_y; j < yStart + ny*nprocy; j++){
for (int i = 0; i<Nx; i++){ for (int i = 0; i<Nx; i++){
if ( (i/checkerSize + k/checkerSize)%2 == 0){ if ( (i/checkerSize + k/checkerSize)%2 == 0){
// void checkers // void checkers
@ -306,9 +302,6 @@ int main(int argc, char **argv)
} }
} }
// Get the rank info
int64_t N = (nx+2)*(ny+2)*(nz+2);
// number of sites to use for periodic boundary condition transition zone // number of sites to use for periodic boundary condition transition zone
int64_t z_transition_size = (nprocz*nz - (Nz - zStart))/2; int64_t z_transition_size = (nprocz*nz - (Nz - zStart))/2;
if (z_transition_size < 0) z_transition_size=0; if (z_transition_size < 0) z_transition_size=0;

66
tests/lbpm_uCT_pp.cpp
View File

@ -61,7 +61,7 @@ int main(int argc, char **argv)
auto L = domain_db->getVector<double>( "L" ); auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" ); auto size = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" ); auto nproc = domain_db->getVector<int>( "nproc" );
int BoundaryCondition = domain_db->getScalar<int>( "BC" ); //int BoundaryCondition = domain_db->getScalar<int>( "BC" );
int nx = size[0]; int nx = size[0];
int ny = size[1]; int ny = size[1];
int nz = size[2]; int nz = size[2];
@ -91,10 +91,10 @@ int main(int argc, char **argv)
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz); printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
printf("Number of MPI ranks used: %i \n", nprocs); printf("Number of MPI ranks used: %i \n", nprocs);
printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz); printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz);
printf("target value = %f \n",target); printf("target value = %f \n",target);
printf("background value = %f \n",background); printf("background value = %f \n",background);
printf("cylinder center = %i, %i, %i \n",center[0],center[1],center[2]); printf("cylinder center = %i, %i, %i \n",center[0],center[1],center[2]);
printf("cylinder radius = %f \n",CylRad); printf("cylinder radius = %f \n",CylRad);
} }
if ( nprocs < nprocx*nprocy*nprocz ){ if ( nprocs < nprocx*nprocy*nprocz ){
ERROR("Insufficient number of processors"); ERROR("Insufficient number of processors");
@ -196,19 +196,19 @@ int main(int argc, char **argv)
filter_src( *Dm[0], LOCVOL[0] ); filter_src( *Dm[0], LOCVOL[0] );
// Set up the mask to be distance to cylinder (crop outside cylinder) // Set up the mask to be distance to cylinder (crop outside cylinder)
if (rank==0) printf("Cropping with cylinder: %i, %i, %i, radius=%f \n",Dm[0]->nprocx()*Nx[0],Dm[0]->nprocy()*Ny[0],Dm[0]->nprocz()*Nz[0],CylRad); if (rank==0) printf("Cropping with cylinder: %i, %i, %i, radius=%f \n",Dm[0]->nprocx()*Nx[0],Dm[0]->nprocy()*Ny[0],Dm[0]->nprocz()*Nz[0],CylRad);
for (int k=0;k<Nz[0]+2;k++) { for (int k=0;k<Nz[0]+2;k++) {
for (int j=0;j<Ny[0]+2;j++) { for (int j=0;j<Ny[0]+2;j++) {
for (int i=0;i<Nx[0]+2;i++) { for (int i=0;i<Nx[0]+2;i++) {
float x= float(Dm[0]->iproc()*Nx[0]+i-1); //float x= float(Dm[0]->iproc()*Nx[0]+i-1);
float y= float (Dm[0]->jproc()*Ny[0]+j-1); float y= float (Dm[0]->jproc()*Ny[0]+j-1);
float z= float(Dm[0]->kproc()*Nz[0]+k-1); float z= float(Dm[0]->kproc()*Nz[0]+k-1);
float cx = float(center[0] - offset[0]); //float cx = float(center[0] - offset[0]);
float cy = float(center[1] - offset[1]); float cy = float(center[1] - offset[1]);
float cz = float(center[2] - offset[2]); float cz = float(center[2] - offset[2]);
// distance from the center line // distance from the center line
MASK(i,j,k) = sqrt((z-cz)*(z-cz) + (y-cy)*(y-cy)); MASK(i,j,k) = sqrt((z-cz)*(z-cz) + (y-cy)*(y-cy));
//if (sqrt(((z-cz)*(z-cz) + (y-cy)*(y-cy)) ) > CylRad) LOCVOL[0](i,j,k)=background; //if (sqrt(((z-cz)*(z-cz) + (y-cy)*(y-cy)) ) > CylRad) LOCVOL[0](i,j,k)=background;
} }
} }
} }
@ -219,18 +219,18 @@ int main(int argc, char **argv)
float THRESHOLD=0.5*(target+background); float THRESHOLD=0.5*(target+background);
float mean_plus=0; float mean_plus=0;
float mean_minus=0; float mean_minus=0;
float min_value = LOCVOL[0](0); float min_value = LOCVOL[0](0);
float max_value = LOCVOL[0](0); float max_value = LOCVOL[0](0);
int count_plus=0; int count_plus=0;
int count_minus=0; int count_minus=0;
for (int k=1;k<Nz[0]+1;k++) { for (int k=1;k<Nz[0]+1;k++) {
for (int j=1;j<Ny[0]+1;j++) { for (int j=1;j<Ny[0]+1;j++) {
for (int i=1;i<Nx[0]+1;i++) { for (int i=1;i<Nx[0]+1;i++) {
//LOCVOL[0](i,j,k) = MASK(i,j,k); //LOCVOL[0](i,j,k) = MASK(i,j,k);
if (MASK(i,j,k) < CylRad ){ if (MASK(i,j,k) < CylRad ){
auto tmp = LOCVOL[0](i,j,k); auto tmp = LOCVOL[0](i,j,k);
/* if ((tmp-background)*(tmp-target) > 0){ /* if ((tmp-background)*(tmp-target) > 0){
// direction to background / target is the same // direction to background / target is the same
if (fabs(tmp-target) > fabs(tmp-background)) tmp=background; // tmp closer to background if (fabs(tmp-target) > fabs(tmp-background)) tmp=background; // tmp closer to background
@ -241,20 +241,20 @@ int main(int argc, char **argv)
mean_plus += tmp; mean_plus += tmp;
count_plus++; count_plus++;
} }
else { else {
mean_minus += tmp; mean_minus += tmp;
count_minus++; count_minus++;
} }
if (tmp < min_value) min_value = tmp; if (tmp < min_value) min_value = tmp;
if (tmp > max_value) max_value = tmp; if (tmp > max_value) max_value = tmp;
} }
} }
} }
} }
count_plus=sumReduce( Dm[0]->Comm, count_plus); count_plus=sumReduce( Dm[0]->Comm, count_plus);
count_minus=sumReduce( Dm[0]->Comm, count_minus); count_minus=sumReduce( Dm[0]->Comm, count_minus);
if (rank==0) printf("minimum value=%f, max value=%f \n",min_value,max_value); if (rank==0) printf("minimum value=%f, max value=%f \n",min_value,max_value);
if (rank==0) printf("plus=%i, minus=%i \n",count_plus,count_minus); if (rank==0) printf("plus=%i, minus=%i \n",count_plus,count_minus);
ASSERT( count_plus > 0 && count_minus > 0 ); ASSERT( count_plus > 0 && count_minus > 0 );
MPI_Barrier(comm); MPI_Barrier(comm);
mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / count_plus; mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / count_plus;
@ -262,25 +262,25 @@ int main(int argc, char **argv)
MPI_Barrier(comm); MPI_Barrier(comm);
if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus); if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus);
//if (rank==0) printf("Scale the input data (size = %i) \n",LOCVOL[0].length()); //if (rank==0) printf("Scale the input data (size = %i) \n",LOCVOL[0].length());
for (size_t i=0; i<LOCVOL[0].length(); i++) { for (size_t i=0; i<LOCVOL[0].length(); i++) {
if ( MASK(i) > CylRad ){ if ( MASK(i) > CylRad ){
LOCVOL[0](i)=background; LOCVOL[0](i)=background;
} }
if ( LOCVOL[0](i) >= THRESHOLD ) { if ( LOCVOL[0](i) >= THRESHOLD ) {
auto tmp = LOCVOL[0](i)/ mean_plus; auto tmp = LOCVOL[0](i)/ mean_plus;
LOCVOL[0](i) = std::min( tmp, 1.0f ); LOCVOL[0](i) = std::min( tmp, 1.0f );
} }
else { else {
auto tmp = -LOCVOL[0](i)/mean_minus; auto tmp = -LOCVOL[0](i)/mean_minus;
LOCVOL[0](i) = std::max( tmp, -1.0f ); LOCVOL[0](i) = std::max( tmp, -1.0f );
} }
//LOCVOL[0](i) = MASK(i); //LOCVOL[0](i) = MASK(i);
} }
// Fill the source data for the coarse meshes // Fill the source data for the coarse meshes
if (rank==0) printf("Coarsen the mesh for N_levels=%i \n",N_levels); if (rank==0) printf("Coarsen the mesh for N_levels=%i \n",N_levels);
MPI_Barrier(comm); MPI_Barrier(comm);
PROFILE_START("CoarsenMesh"); PROFILE_START("CoarsenMesh");
for (int i=1; i<N_levels; i++) { for (int i=1; i<N_levels; i++) {
Array<float> filter(ratio[0],ratio[1],ratio[2]); Array<float> filter(ratio[0],ratio[1],ratio[2]);

10
tests/pmmc_cylinder.cpp
View File

@ -11,17 +11,11 @@
int main (int argc, char **argv) int main (int argc, char **argv)
{ {
// printf("Radius = %s \n,"RADIUS); // printf("Radius = %s \n,"RADIUS);
int SIZE = N*N*N;
int Nx,Ny,Nz; int Nx,Ny,Nz;
Nx = Ny = Nz = N; Nx = Ny = Nz = N;
int i,j,k,p,q,r; int i,j,k,p,q,r;
// double *Solid; // cylinder
// double *Phase; // region of the cylinder
// Solid = new double [SIZE];
// Phase = new double [SIZE];
DoubleArray SignDist(Nx,Ny,Nz); DoubleArray SignDist(Nx,Ny,Nz);
DoubleArray Phase(Nx,Ny,Nz); DoubleArray Phase(Nx,Ny,Nz);
double fluid_isovalue = 0.0; double fluid_isovalue = 0.0;
@ -36,9 +30,6 @@ int main (int argc, char **argv)
//........................................................................... //...........................................................................
double awn,ans,aws,lwns,nwp_volume; double awn,ans,aws,lwns,nwp_volume;
double As; double As;
double dEs,dAwn,dAns; // Global surface energy (calculated by rank=0)
double awn_global,ans_global,aws_global,lwns_global,nwp_volume_global;
double As_global;
// bool add=1; // Set to false if any corners contain nw-phase ( F > fluid_isovalue) // bool add=1; // Set to false if any corners contain nw-phase ( F > fluid_isovalue)
int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners
// int count_in=0,count_out=0; // int count_in=0,count_out=0;
@ -75,7 +66,6 @@ int main (int argc, char **argv)
int n_local_nws_pts; int n_local_nws_pts;
int c; int c;
int newton_steps = 0;
//........................................................................... //...........................................................................
int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo
IntArray cubeList(3,ncubes); IntArray cubeList(3,ncubes);