diff --git a/common/Domain.h b/common/Domain.h index 6f9bd84f..e31c3459 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -34,6 +34,7 @@ struct Domain{ BlobLabel.resize(Nx,Ny,Nz); BlobGraph.resize(18,MAX_BLOB_COUNT,MAX_BLOB_COUNT); BoundaryCondition = BC; + rank_info=RankInfoStruct(rank,nprocx,nprocy,nprocz); } ~Domain(); @@ -41,9 +42,10 @@ struct Domain{ int Nx,Ny,Nz,N; int iproc,jproc,kproc; int nprocx,nprocy,nprocz; - double Lx,Ly,Lz,Volume; + double Lx,Ly,Lz,Volume; int rank; int BoundaryCondition; + RankInfoStruct rank_info; MPI_Group Group; // Group of processors associated with this domain MPI_Comm Comm; // MPI Communicator for this domain @@ -775,6 +777,142 @@ void Domain::BlobComm(MPI_Comm Communicator){ UnpackBlobData(recvList_YZ, recvCount_YZ ,recvBuf_YZ, BlobLabelData); //...................................................................................... } +inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ + /* + * This routine converts the data in the Distance array to a signed distance + * by solving the equation df/dt = sign(1-|grad f|), where Distance provides + * the values of f on the mesh associated with domain Dm + * It has been tested with segmented data initialized to values [-1,1] + * and will converge toward the signed distance to the surface bounding the associated phases + */ + + int Q=26; + int q,i,j,k,n; + double dt=0.1; + int in,jn,kn,nn; + double Dqx,Dqy,Dqz,Dx,Dy,Dz,W; + double nx,ny,nz,Cqx,Cqy,Cqz,sign,norm; + + const static int D3Q27[26][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,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},{1,1,1},{-1,-1,-1},{1,1,-1},{-1,-1,1}, + {-1,1,-1},{1,-1,1},{1,-1,-1},{-1,1,1}}; + + double weights[26]; + // Compute the weights from the finite differences + for (q=0; q fillData(Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1); + + int count = 0; + while (count < timesteps){ + + // Communicate the halo of values + fillData.fill(Distance); + + // Execute the next timestep + for (k=1;k 0.0 && !(Distance(i,j,k)*Distance(i,j,k) < 1.0) ){ + for (q=0; q<26; q++){ + Cqx = 1.0*D3Q27[q][0]; + Cqy = 1.0*D3Q27[q][1]; + Cqz = 1.0*D3Q27[q][2]; + // get the associated neighbor + in = i + D3Q27[q][0]; + jn = j + D3Q27[q][1]; + kn = k + D3Q27[q][2]; + + // make sure the neighbor is in the domain (periodic BC) + /* if (in < 0 ) in +=Nx; + * don't need this in parallel since MPI handles the halos + if (jn < 0 ) jn +=Ny; + if (kn < 0 ) kn +=Nz; + if (!(in < Nx) ) in -=Nx; + if (!(jn < Ny) ) jn -=Ny; + if (!(kn < Nz) ) kn -=Nz; + // symmetric boundary + if (in < 0 ) in = i; + if (jn < 0 ) jn = j; + if (kn < 0 ) kn = k; + if (!(in < Nx) ) in = i; + if (!(jn < Ny) ) jn = k; + if (!(kn < Nz) ) kn = k; + */ + + // Compute the gradient using upwind finite differences + Dqx = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqx; + Dqy = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqy; + Dqz = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqz; + + // Only include upwind derivatives + if (sign*(nx*Cqx + ny*Cqy + nz*Cqz) < 0.0 ){ + + Dx += Dqx; + Dy += Dqy; + Dz += Dqz; + W += weights[q]; + } + } + // Normalize by the weight to get the approximation to the gradient + Dx /= W; + Dy /= W; + Dz /= W; + + norm = sqrt(Dx*Dx+Dy*Dy+Dz*Dz); + } + else{ + norm = 0.0; + } + Distance(i,j,k) += dt*sign*(1.0 - norm); + + // Disallow any change in phase + if (Distance(i,j,k)*2.0*(ID[n]-1.0) < 0) Distance(i,j,k) = -Distance(i,j,k); + } + } + } + + count++; + } +} + inline void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad) { diff --git a/common/TwoPhase.h b/common/TwoPhase.h index afb70172..9b356892 100644 --- a/common/TwoPhase.h +++ b/common/TwoPhase.h @@ -550,7 +550,7 @@ void TwoPhase::ComputeLocalBlob(){ int i,j,k,n,label; double vF,vS; vF = 0.0; vS= -1.0; - const RankInfoStruct rank_info(Dm.rank,Dm.nprocx,Dm.nprocy,Dm.nprocz); +// const RankInfoStruct rank_info(Dm.rank,Dm.nprocx,Dm.nprocy,Dm.nprocz); 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}}; // get the maximum label locally -- then compute number of global blobs @@ -563,7 +563,7 @@ void TwoPhase::ComputeLocalBlob(){ nblobs_global+=1; if (Dm.rank==0) printf("Number of blobs is %i \n",nblobs_global); - nblobs_global = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info, + nblobs_global = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,Dm.rank_info, Phase,SDs,vF,vS,BlobLabel); if (Dm.rank==0) printf("Number of blobs is %i \n",nblobs_global); diff --git a/tests/BlobAnalyzeParallel.cpp b/tests/BlobAnalyzeParallel.cpp index ac0d9093..cc86a15b 100644 --- a/tests/BlobAnalyzeParallel.cpp +++ b/tests/BlobAnalyzeParallel.cpp @@ -7,7 +7,9 @@ #include #include "common/Communication.h" #include "analysis/analysis.h" -#include "ProfilerApp.h" +#ifdef PROFILE + #include "ProfilerApp.h" +#endif #include "TwoPhase.h" //#include "Domain.h" @@ -98,10 +100,12 @@ int main(int argc, char **argv) MPI_Init(&argc,&argv); MPI_Comm_rank(MPI_COMM_WORLD,&rank); MPI_Comm_size(MPI_COMM_WORLD,&nprocs); - PROFILE_ENABLE(0); +#ifdef PROFILE + PROFILE_ENABLE(0); PROFILE_DISABLE_TRACE(); PROFILE_SYNCHRONIZE(); PROFILE_START("main"); +#endif if ( rank==0 ) { printf("-----------------------------------------------------------\n"); @@ -147,7 +151,7 @@ int main(int argc, char **argv) int BC=0; // Get the rank info Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); - const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); + // const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); TwoPhase Averages(Dm); int N = (nx+2)*(ny+2)*(nz+2); @@ -157,6 +161,7 @@ int main(int argc, char **argv) readRankData( rank, nx+2, ny+2, nz+2, Phase, SignDist ); // Communication the halos + const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); fillHalo fillData(rank_info,nx,ny,nz,1,1,1,0,1); fillData.fill(Phase); fillData.fill(SignDist); @@ -166,7 +171,7 @@ int main(int argc, char **argv) double vF=0.0; double vS=0.0; IntArray GlobalBlobID; - int nblobs = ComputeGlobalBlobIDs(nx,ny,nz,rank_info, + int nblobs = ComputeGlobalBlobIDs(nx,ny,nz,Dm.rank_info, Phase,SignDist,vF,vS,GlobalBlobID); if ( rank==0 ) { printf("Identified %i blobs\n",nblobs); } @@ -191,9 +196,6 @@ int main(int argc, char **argv) //......................................................................... // Populate the arrays needed to perform averaging if (rank==0) printf("Populate arrays \n"); - // Compute porosity - double porosity,sum,sum_global; - sum=0.0; for (int k=0; k 0.0){ - Dm.id[n]=1; - sum += 1.0; - } - else Dm.id[n]=0; + } } } delete [] DistEven; delete [] DistOdd; + // Compute porosity + double porosity,sum,sum_global; + sum=0.0; + for (int k=1; k 0.0){ + Dm.id[n]=1; + sum += 1.0; + } + else Dm.id[n]=0; + } + } + } MPI_Allreduce(&sum,&sum_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); porosity = sum_global/Dm.Volume; if (rank==0) printf("Porosity = %f \n",porosity); @@ -366,9 +379,10 @@ int main(int argc, char **argv) /*FILE *BLOBS = fopen("Blobs.dat","wb"); fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS); fclose(BLOBS);*/ - +#ifdef PROFILE PROFILE_STOP("main"); PROFILE_SAVE("BlobIdentifyParallel",false); +#endif MPI_Barrier(MPI_COMM_WORLD); MPI_Finalize(); return 0; diff --git a/tests/BlobIdentifyParallel.cpp b/tests/BlobIdentifyParallel.cpp index e772f5da..7ceae67c 100644 --- a/tests/BlobIdentifyParallel.cpp +++ b/tests/BlobIdentifyParallel.cpp @@ -8,8 +8,9 @@ #include "common/pmmc.h" #include "common/Communication.h" #include "analysis/analysis.h" -#include "ProfilerApp.h" - +#ifdef PROFILE + #include "ProfilerApp.h" +#endif //#include "Domain.h" @@ -50,10 +51,12 @@ int main(int argc, char **argv) MPI_Init(&argc,&argv); MPI_Comm_rank(MPI_COMM_WORLD,&rank); MPI_Comm_size(MPI_COMM_WORLD,&nprocs); - PROFILE_ENABLE(0); +#ifdef PROFILE + PROFILE_ENABLE(0); PROFILE_DISABLE_TRACE(); PROFILE_SYNCHRONIZE(); PROFILE_START("main"); +#endif if ( rank==0 ) { printf("-----------------------------------------------------------\n"); @@ -121,10 +124,11 @@ int main(int argc, char **argv) /*FILE *BLOBS = fopen("Blobs.dat","wb"); fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS); fclose(BLOBS);*/ - +#ifdef PROFILE PROFILE_STOP("main"); PROFILE_SAVE("BlobIdentifyParallel",false); - MPI_Barrier(MPI_COMM_WORLD); +#endif + MPI_Barrier(MPI_COMM_WORLD); MPI_Finalize(); return 0; } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 71547f4f..de72bc88 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -22,6 +22,7 @@ ADD_LBPM_TEST( TestSphereCurvature ) ADD_LBPM_TEST( TestContactAngle ) ADD_LBPM_TEST_1_2_4( TestTwoPhase ) ADD_LBPM_TEST_PARALLEL( TestTwoPhase 8 ) +ADD_LBPM_TEST_PARALLEL( TestSegDist 8 ) ADD_LBPM_TEST_1_2_4( testCommunication ) ADD_LBPM_TEST_1_2_4( testUtilities ) ADD_LBPM_TEST( TestWriter ) diff --git a/tests/TestComputeSignDist.cpp b/tests/TestComputeSignDist.cpp index e444190f..62491b1d 100644 --- a/tests/TestComputeSignDist.cpp +++ b/tests/TestComputeSignDist.cpp @@ -11,13 +11,13 @@ #include #include -inline void SSO(DoubleArray &Distance, IntArray &ID, int timesteps){ +inline void SSO(DoubleArray &Distance, char *ID, int timesteps){ int Q=26; - int q,i,j,k; - int Nx = ID.m; - int Ny = ID.n; - int Nz = ID.o; + int q,i,j,k,n; + int Nx = Distance.m; + int Ny = Distance.n; + int Nz = Distance.o; const static int D3Q27[26][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,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},{1,1,1},{-1,-1,-1},{1,1,-1},{-1,-1,1}, @@ -30,40 +30,58 @@ inline void SSO(DoubleArray &Distance, IntArray &ID, int timesteps){ } // Initialize the Distance from ID - for (i=0; i 0.0){ - for (q=0; q<27; q++){ + for (q=0; q<26; q++){ Cqx = 1.0*D3Q27[q][0]; Cqy = 1.0*D3Q27[q][1]; Cqz = 1.0*D3Q27[q][2]; @@ -72,12 +90,19 @@ inline void SSO(DoubleArray &Distance, IntArray &ID, int timesteps){ jn = j + D3Q27[q][1]; kn = k + D3Q27[q][2]; // make sure the neighbor is in the domain (periodic BC) - if (in < 0 ) in +=Nx; + /* if (in < 0 ) in +=Nx; if (jn < 0 ) jn +=Ny; if (kn < 0 ) kn +=Nz; if (!(in < Nx) ) in -=Nx; if (!(jn < Ny) ) jn -=Ny; if (!(kn < Nz) ) kn -=Nz; + */ // symmetric boundary + if (in < 0 ) in = i; + if (jn < 0 ) jn = j; + if (kn < 0 ) kn = k; + if (!(in < Nx) ) in = i; + if (!(jn < Ny) ) jn = k; + if (!(kn < Nz) ) kn = k; // 1-D index nn = kn*Nx*Ny + jn*Nx + in; @@ -103,20 +128,18 @@ inline void SSO(DoubleArray &Distance, IntArray &ID, int timesteps){ norm = sqrt(Dx*Dx+Dy*Dy+Dz*Dz); } else{ - norm = 0.7; + norm = 0.0; } Distance.data[n] += dt*sign*(1.0 - norm); - // Disallow any change in phase position - if (Distance.data[n]*ID.data[n] < 0) Distance.data[n] = -Distance.data[n]; + // Disallow any change in phase + if (Distance.data[n]*2.0*(ID[n]-1.0) < 0) Distance.data[n] = -Distance.data[n]; } } } count++; } - - } int main(){ @@ -136,28 +159,31 @@ int main(){ double fxm,fym,fzm,fxp,fyp,fzp; double fxy,fXy,fxY,fXY,fxz,fXz,fxZ,fXZ,fyz,fYz,fyZ,fYZ; double nx,ny,nz; - int ip,im,jp,jm,kp,km; int count = 0; double sign; - double dt=0.01; + double dt=1.0; printf("Nx=%i, Ny=%i, Nz= %i, \n",Nx,Ny,Nz); - short int *id; + char *id; #ifdef READMEDIA Nx = 347; Ny = 347; Nz = 235; + Nx = 512; + Ny = 512; + Nz = 512; N = Nx*Ny*Nz; - id = new short int [N]; - FILE *INPUT = fopen("Solid.dat",'rb'); - fread(id,4,N*,Ny*Nz,INPUT) + id = new char [N]; + + FILE *INPUT = fopen("Solid.dat","rb"); + fread(id,1,Nx*Ny*Nz,INPUT); fclose(INPUT); #else - id = new short int [N]; + id = new char [N]; double BubbleRadius = 5; // Initialize the bubble for (k=0;k 1.0e-6){ + dt=1.0; + while (count < 10 && dt > 1.0e-6){ err = 0.0; for (k=0;k 0){ + if (id[n] < 0){ if ( nx > 0.0) f_x = fxp; else f_x = fxm; if ( ny > 0.0) f_y = fyp; @@ -348,8 +391,8 @@ int main(){ count++; } +*/ - */ FILE *DIST; DIST = fopen("SignDist","wb"); fwrite(Distance.data,8,N,DIST); diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp new file mode 100644 index 00000000..841c4f42 --- /dev/null +++ b/tests/TestSegDist.cpp @@ -0,0 +1,106 @@ +// Compute the signed distance from a digitized image +// Two phases are present +// Phase 1 has value -1 +// Phase 2 has value 1 +// this code uses the segmented image to generate the signed distance + +#include +#include +#include +#include +#include +#include +#include + +int main(int argc, char **argv) +{ + // Initialize MPI + int rank, nprocs; + MPI_Init(&argc,&argv); + MPI_Comm_rank(MPI_COMM_WORLD,&rank); + MPI_Comm_size(MPI_COMM_WORLD,&nprocs); + + int i,j,k,n,nn; + int iproc,jproc,kproc; + int nx,ny,nz; + int Nx, Ny, Nz, N; + int nprocx, nprocy, nprocz, nspheres; + double Lx, Ly, Lz; + Nx = Ny = Nz = 50; + nx = ny = nz = 50; + N = Nx*Ny*Nz; + nprocx=nprocy=nprocz=2; + Lx = Ly = Lz = 1.0; + int BC=0; + + + if (nprocs != 8){ + ERROR("TestSegDist: Number of MPI processes must be equal to 8"); + } + + if (nprocx !=2 || nprocz !=2 || nprocy !=2 ){ + ERROR("TestSegDist: MPI process grid must be 2x2x2"); + } + + // Get the rank info + Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); + Dm.CommInit(MPI_COMM_WORLD); + + nx+=2; ny+=2; nz+=2; + N = nx*ny*nz; + int count = 0; + + char *id; + id = new char [N]; + double BubbleRadius = 25; + // Initialize the bubble + int x,y,z; + for (k=1;k