From bf98d9aa232d49e6d70e53a202a6c11cfce04cf8 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 15:37:39 -0400 Subject: [PATCH 01/28] Adding support for signed distance from segmented images in parallel --- common/Domain.h | 4 +- common/TwoPhase.h | 118 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 121 insertions(+), 1 deletion(-) diff --git a/common/Domain.h b/common/Domain.h index 6f9bd84f..738f99a2 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; + const RankInfoStruct; MPI_Group Group; // Group of processors associated with this domain MPI_Comm Comm; // MPI Communicator for this domain diff --git a/common/TwoPhase.h b/common/TwoPhase.h index afb70172..1fffd915 100644 --- a/common/TwoPhase.h +++ b/common/TwoPhase.h @@ -289,6 +289,7 @@ public: void UpdateMeshValues(); void UpdateSolid(); void ComputeDelPhi(); + void SSO(DoubleArray &Distance, char *ID, int timesteps); void ColorToSignedDistance(double Beta, double *ColorData, double *DistData); void ComputeLocal(); void ComputeLocalBlob(); @@ -300,6 +301,123 @@ public: }; +inline void TwoPhase::SSO(DoubleArray &Distance, char *ID, int timesteps){ + + int Q=26; + int q,i,j,k,n; + 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 0.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; + 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; + + // 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++; + } +} + + void TwoPhase::ColorToSignedDistance(double Beta, double *ColorData, double *DistData){ double factor,temp,value; From f56a4410d9d056e4ad1fc2721252f07fd2984faa Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 15:41:31 -0400 Subject: [PATCH 02/28] Adding support for signed distance from segmented images in parallel --- tests/BlobAnalyzeParallel.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/BlobAnalyzeParallel.cpp b/tests/BlobAnalyzeParallel.cpp index ac0d9093..693dff04 100644 --- a/tests/BlobAnalyzeParallel.cpp +++ b/tests/BlobAnalyzeParallel.cpp @@ -7,7 +7,7 @@ #include #include "common/Communication.h" #include "analysis/analysis.h" -#include "ProfilerApp.h" +//#include "ProfilerApp.h" #include "TwoPhase.h" //#include "Domain.h" @@ -98,10 +98,10 @@ 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); - PROFILE_DISABLE_TRACE(); - PROFILE_SYNCHRONIZE(); - PROFILE_START("main"); + // PROFILE_ENABLE(0); + // PROFILE_DISABLE_TRACE(); + // PROFILE_SYNCHRONIZE(); + // PROFILE_START("main"); if ( rank==0 ) { printf("-----------------------------------------------------------\n"); @@ -367,8 +367,8 @@ int main(int argc, char **argv) fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS); fclose(BLOBS);*/ - PROFILE_STOP("main"); - PROFILE_SAVE("BlobIdentifyParallel",false); +// PROFILE_STOP("main"); +// PROFILE_SAVE("BlobIdentifyParallel",false); MPI_Barrier(MPI_COMM_WORLD); MPI_Finalize(); return 0; From 5c933c15bbfc4f79452f6a54c84dde8f18946c99 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 15:42:40 -0400 Subject: [PATCH 03/28] Adding support for signed distance from segmented images in parallel --- common/Domain.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common/Domain.h b/common/Domain.h index 738f99a2..42de4272 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -45,7 +45,7 @@ struct Domain{ double Lx,Ly,Lz,Volume; int rank; int BoundaryCondition; - const RankInfoStruct; + const RankInfoStruct rank_info; MPI_Group Group; // Group of processors associated with this domain MPI_Comm Comm; // MPI Communicator for this domain From e5f14f99a6407fad68c5e6a10224bc31562b513f Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 15:44:49 -0400 Subject: [PATCH 04/28] Adding support for signed distance from segmented images in parallel -- debugging --- common/Domain.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common/Domain.h b/common/Domain.h index 42de4272..991a8ad8 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -45,7 +45,7 @@ struct Domain{ double Lx,Ly,Lz,Volume; int rank; int BoundaryCondition; - const RankInfoStruct rank_info; + RankInfoStruct rank_info; MPI_Group Group; // Group of processors associated with this domain MPI_Comm Comm; // MPI Communicator for this domain From c2d84c4ff1fd2d1b39d1c30514ce48ea44e3389f Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 15:49:43 -0400 Subject: [PATCH 05/28] Added compiler directives for PROFILING in tests/BlobAnalyzieParallel --- tests/BlobAnalyzeParallel.cpp | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/tests/BlobAnalyzeParallel.cpp b/tests/BlobAnalyzeParallel.cpp index 693dff04..c9298b87 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); - // PROFILE_DISABLE_TRACE(); - // PROFILE_SYNCHRONIZE(); - // PROFILE_START("main"); +#ifdef PROFILE + PROFILE_ENABLE(0); + PROFILE_DISABLE_TRACE(); + PROFILE_SYNCHRONIZE(); + PROFILE_START("main"); +#endif if ( rank==0 ) { printf("-----------------------------------------------------------\n"); @@ -366,9 +370,10 @@ int main(int argc, char **argv) /*FILE *BLOBS = fopen("Blobs.dat","wb"); fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS); fclose(BLOBS);*/ - -// PROFILE_STOP("main"); -// PROFILE_SAVE("BlobIdentifyParallel",false); +#ifdef PROFILE + PROFILE_STOP("main"); + PROFILE_SAVE("BlobIdentifyParallel",false); +#endif MPI_Barrier(MPI_COMM_WORLD); MPI_Finalize(); return 0; From 16061f12950d7ea5a8140b6afc6d68ba8681be8c Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 15:52:53 -0400 Subject: [PATCH 06/28] Added compiler directives for PROFILING in tests/BlobIdentifyParallel --- tests/BlobIdentifyParallel.cpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) 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; } From 57817dc94aff6402b964ff29dd1a68512326fb7e Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 16:02:47 -0400 Subject: [PATCH 07/28] Addings support for parallel signed distance froms segmented -- pieces in place --- common/TwoPhase.h | 37 +++++++++++++++++++++++------------ tests/BlobAnalyzeParallel.cpp | 6 +++--- 2 files changed, 28 insertions(+), 15 deletions(-) diff --git a/common/TwoPhase.h b/common/TwoPhase.h index 1fffd915..1dec4a6f 100644 --- a/common/TwoPhase.h +++ b/common/TwoPhase.h @@ -289,7 +289,7 @@ public: void UpdateMeshValues(); void UpdateSolid(); void ComputeDelPhi(); - void SSO(DoubleArray &Distance, char *ID, int timesteps); + void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps); void ColorToSignedDistance(double Beta, double *ColorData, double *DistData); void ComputeLocal(); void ComputeLocalBlob(); @@ -301,7 +301,7 @@ public: }; -inline void TwoPhase::SSO(DoubleArray &Distance, char *ID, int timesteps){ +inline void TwoPhase::SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ int Q=26; int q,i,j,k,n; @@ -323,20 +323,27 @@ inline void TwoPhase::SSO(DoubleArray &Distance, char *ID, int timesteps){ double nx,ny,nz,Cqx,Cqy,Cqz,sign,norm; // double f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; + fillHalo fillData(Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1); + + printf("Number of timesteps is %i \n",timesteps); printf("Mesh is %i,%i,%i \n",Nx,Ny,Nz); while (count < timesteps){ printf("count=%i \n",count); - for (k=0;k 0.0){ @@ -364,20 +373,24 @@ inline void TwoPhase::SSO(DoubleArray &Distance, char *ID, int timesteps){ 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 + // 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; @@ -668,7 +681,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 @@ -681,7 +694,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 c9298b87..af013a5a 100644 --- a/tests/BlobAnalyzeParallel.cpp +++ b/tests/BlobAnalyzeParallel.cpp @@ -151,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); @@ -161,7 +161,7 @@ int main(int argc, char **argv) readRankData( rank, nx+2, ny+2, nz+2, Phase, SignDist ); // Communication the halos - fillHalo fillData(rank_info,nx,ny,nz,1,1,1,0,1); + fillHalo fillData(Dm.rank_info,nx,ny,nz,1,1,1,0,1); fillData.fill(Phase); fillData.fill(SignDist); @@ -170,7 +170,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); } From e7eb11f2282b3b6a8310e84f614aa3dfecbebaf1 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 17:54:16 -0400 Subject: [PATCH 08/28] Moved signed distance computation SSO from TwoPhase.h to Domain.h --- common/Domain.h | 137 ++++++++++++++++++++++++++++++++++++++++++++++ common/TwoPhase.h | 131 -------------------------------------------- 2 files changed, 137 insertions(+), 131 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index 991a8ad8..70a8217b 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -777,6 +777,143 @@ 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; + 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,Nx-2,Ny-2,Nz-2,1,1,1,0,1); + + while (count < timesteps){ + + printf("count=%i \n",count); + + // 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; + */ + + // 1-D index + nn = kn*Nx*Ny + jn*Nx + in; + + // 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 1dec4a6f..9b356892 100644 --- a/common/TwoPhase.h +++ b/common/TwoPhase.h @@ -289,7 +289,6 @@ public: void UpdateMeshValues(); void UpdateSolid(); void ComputeDelPhi(); - void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps); void ColorToSignedDistance(double Beta, double *ColorData, double *DistData); void ComputeLocal(); void ComputeLocalBlob(); @@ -301,136 +300,6 @@ public: }; -inline void TwoPhase::SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ - - int Q=26; - int q,i,j,k,n; - 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,Nx-2,Ny-2,Nz-2,1,1,1,0,1); - - - printf("Number of timesteps is %i \n",timesteps); - printf("Mesh is %i,%i,%i \n",Nx,Ny,Nz); - - while (count < timesteps){ - - printf("count=%i \n",count); - - // Communicate the halo of values - fillData.fill(Distance); - - for (k=1;k 0.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; - */ - - // 1-D index - nn = kn*Nx*Ny + jn*Nx + in; - - // 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++; - } -} - - void TwoPhase::ColorToSignedDistance(double Beta, double *ColorData, double *DistData){ double factor,temp,value; From ae9b8482965172c45094601a1fddc783627deee8 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 17:55:25 -0400 Subject: [PATCH 09/28] Added tests/TestSegDist.cpp to test parallel converseion for segemented data to signed distance function --- tests/TestComputeSignDist.cpp | 113 +++++++++++++++++++++++----------- tests/TestSegDist.cpp | 108 ++++++++++++++++++++++++++++++++ 2 files changed, 186 insertions(+), 35 deletions(-) create mode 100644 tests/TestSegDist.cpp 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..c659abaf --- /dev/null +++ b/tests/TestSegDist.cpp @@ -0,0 +1,108 @@ +// 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 + +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 nx,ny,nz; + int Nx, Ny, Nz, N; + Nx = Ny = Nz = 50; + N = Nx*Ny*Nz; + + if (nprocs != 8){ + ERROR("TestSegDist: Number of MPI processes must be equal to 8"); + } + + //....................................................................... + // Reading the domain information file + //....................................................................... + int nprocx, nprocy, nprocz, nx, ny, nz, nspheres; + double Lx, Ly, Lz; + ifstream domain("Domain.in"); + domain >> nprocx; + domain >> nprocy; + domain >> nprocz; + domain >> nx; + domain >> ny; + domain >> nz; + domain >> nspheres; + domain >> Lx; + domain >> Ly; + domain >> Lz; + + nprocx=nprocy=nprocz=2; + if (nprocx !=2 || nprocz !=2 || nprocy !=2 ){ + ERROR("TestSegDist: MPI process grid must be 2x2x2"); + } + + int BC=0; + // Get the rank info + Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); + nx+=2; ny+=2; nz+=2; + int count = 0; + + char *id; + id = new char [N]; + double BubbleRadius = 5; + // Initialize the bubble + int x,y,z; + for (k=1;k Date: Wed, 3 Jun 2015 17:57:59 -0400 Subject: [PATCH 10/28] Added unit test for segmented -> signed distance --- tests/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) 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 ) From ce199b4d3d5fb76774cfd55e5dc72bd4f6e18c61 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 18:04:43 -0400 Subject: [PATCH 11/28] Fixed bug in Domain.h: SSo --- common/Domain.h | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index 70a8217b..b5c39be6 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -788,6 +788,11 @@ inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ int Q=26; int q,i,j,k,n; + double dt=0.25; + 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}, @@ -799,14 +804,9 @@ inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ weights[q] = sqrt(1.0*(D3Q27[q][0]*D3Q27[q][0]) + 1.0*(D3Q27[q][1]*D3Q27[q][1]) + 1.0*(D3Q27[q][2]*D3Q27[q][2])); } + fillHalo fillData(Dm.rank_info,Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,1,1,1,0,1); + int count = 0; - double dt=0.25; - int in,jn,kn,nn; - double Dqx,Dqy,Dqz,Dx,Dy,Dz,W; - double nx,ny,nz,Cqx,Cqy,Cqz,sign,norm; - - fillHalo fillData(Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1); - while (count < timesteps){ printf("count=%i \n",count); @@ -815,11 +815,11 @@ inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ fillData.fill(Distance); // Execute the next timestep - for (k=1;k Date: Wed, 3 Jun 2015 21:09:41 -0400 Subject: [PATCH 12/28] Added headers to tests/TestSegDist --- tests/TestSegDist.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index c659abaf..8af0f242 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -10,6 +10,7 @@ #include #include #include +#include int main(int argc, char **argv) { From f2846c0d7b98771506dd913a6aa86a52710f9521 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 21:11:50 -0400 Subject: [PATCH 13/28] debugging tests/TestSegDist --- tests/TestSegDist.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index 8af0f242..a26b431a 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -21,8 +21,11 @@ int main(int argc, char **argv) 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; N = Nx*Ny*Nz; @@ -33,8 +36,6 @@ int main(int argc, char **argv) //....................................................................... // Reading the domain information file //....................................................................... - int nprocx, nprocy, nprocz, nx, ny, nz, nspheres; - double Lx, Ly, Lz; ifstream domain("Domain.in"); domain >> nprocx; domain >> nprocy; From 34765d92c28c8d2f338a9ddfa1d347fa868fa9df Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 3 Jun 2015 21:50:24 -0400 Subject: [PATCH 14/28] debugging tests/TestSegDist --- tests/TestSegDist.cpp | 22 +++++----------------- 1 file changed, 5 insertions(+), 17 deletions(-) diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index a26b431a..579545c5 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -27,33 +27,21 @@ int main(int argc, char **argv) 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"); } - //....................................................................... - // Reading the domain information file - //....................................................................... - ifstream domain("Domain.in"); - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> nx; - domain >> ny; - domain >> nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - - nprocx=nprocy=nprocz=2; if (nprocx !=2 || nprocz !=2 || nprocy !=2 ){ ERROR("TestSegDist: MPI process grid must be 2x2x2"); } - int BC=0; // Get the rank info Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); nx+=2; ny+=2; nz+=2; From 24a11aa45ad2875fa184b2861dbd89866d280d1c Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 4 Jun 2015 06:58:28 -0400 Subject: [PATCH 15/28] Fixing bugs in SSO --- common/Domain.h | 6 ++---- tests/TestSegDist.cpp | 1 + tests/TestTwoPhase.cpp | 1 - 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index b5c39be6..65c07da3 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -34,7 +34,6 @@ 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(); @@ -45,7 +44,6 @@ struct Domain{ 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 @@ -803,8 +801,8 @@ inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ for (q=0; q fillData(Dm.rank_info,Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,1,1,1,0,1); + const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); + fillHalo fillData(rank_info,Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,1,1,1,0,1); int count = 0; while (count < timesteps){ diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index 579545c5..094c0b9c 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -44,6 +44,7 @@ int main(int argc, char **argv) // Get the rank info Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); + nx+=2; ny+=2; nz+=2; int count = 0; diff --git a/tests/TestTwoPhase.cpp b/tests/TestTwoPhase.cpp index 387e4b1e..da673ec5 100644 --- a/tests/TestTwoPhase.cpp +++ b/tests/TestTwoPhase.cpp @@ -106,7 +106,6 @@ int main(int argc, char **argv) Averages.PrintAll(timestep); //.................................................................... - if (rank==0){ FILE *PHASE; PHASE = fopen("Phase.00000","wb"); From 017a40842225baeb61fbfdd60c954c8f22dcf89e Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 4 Jun 2015 07:00:31 -0400 Subject: [PATCH 16/28] debugging SSO --- common/Domain.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common/Domain.h b/common/Domain.h index 65c07da3..010747f9 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -801,7 +801,7 @@ inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ for (q=0; q fillData(rank_info,Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,1,1,1,0,1); int count = 0; From 965257e9cdc021d0294823f1fa4fcb0f964a9ba5 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 4 Jun 2015 07:02:40 -0400 Subject: [PATCH 17/28] debugging SSO --- common/Domain.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/common/Domain.h b/common/Domain.h index 010747f9..51e3dc71 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(); @@ -44,6 +45,7 @@ struct Domain{ 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 From 57d3e282cb7f6785ca3ad68585c6668592c85f98 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 4 Jun 2015 07:26:31 -0400 Subject: [PATCH 18/28] debugging SSO --- common/Domain.h | 8 ++++++-- tests/BlobAnalyzeParallel.cpp | 3 ++- tests/TestSegDist.cpp | 4 ++++ 3 files changed, 12 insertions(+), 3 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index 51e3dc71..c25ae3a5 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -803,8 +803,12 @@ inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ for (q=0; q fillData(rank_info,Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,1,1,1,0,1); + + int xdim,ydim,zdim; + xdim=Dm.Nx-2; + ydim=Dm.Ny-2; + zdim=Dm.Nz-2; + fillHalo fillData(Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1); int count = 0; while (count < timesteps){ diff --git a/tests/BlobAnalyzeParallel.cpp b/tests/BlobAnalyzeParallel.cpp index af013a5a..8f11293a 100644 --- a/tests/BlobAnalyzeParallel.cpp +++ b/tests/BlobAnalyzeParallel.cpp @@ -161,7 +161,8 @@ int main(int argc, char **argv) readRankData( rank, nx+2, ny+2, nz+2, Phase, SignDist ); // Communication the halos - fillHalo fillData(Dm.rank_info,nx,ny,nz,1,1,1,0,1); + 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); diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index 094c0b9c..c1458daa 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -45,6 +45,10 @@ int main(int argc, char **argv) // Get the rank info Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); + if (rank==0){ + printf() + } + nx+=2; ny+=2; nz+=2; int count = 0; From 57f0699bf5279694a8e252d7634199b9bec4eaba Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 4 Jun 2015 07:29:49 -0400 Subject: [PATCH 19/28] debugging SSO --- tests/TestSegDist.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index c1458daa..2946c7f4 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -45,9 +45,6 @@ int main(int argc, char **argv) // Get the rank info Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); - if (rank==0){ - printf() - } nx+=2; ny+=2; nz+=2; int count = 0; From 698da7798b1bfb8fe045a7f6859f615dc3ca0a93 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 4 Jun 2015 07:42:08 -0400 Subject: [PATCH 20/28] debugging SSO --- tests/TestSegDist.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index 2946c7f4..f1e55ddd 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -45,10 +45,9 @@ int main(int argc, char **argv) // Get the rank info Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); - nx+=2; ny+=2; nz+=2; int count = 0; - + char *id; id = new char [N]; double BubbleRadius = 5; @@ -84,6 +83,10 @@ int main(int argc, char **argv) } } + if (rank==0) printf("Nx = %i \n",(int)Distance.size(0)); + if (rank==0) printf("Ny = %i \n",(int)Distance.size(1)); + if (rank==0) printf("Nz = %i \n",(int)Distance.size(2)); + printf("Initialized! Converting to Signed Distance function \n"); SSO(Distance,id,Dm,10); From 3ca1cf649b99d8ac08498d5b43979dcea085fbeb Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 4 Jun 2015 07:44:38 -0400 Subject: [PATCH 21/28] debugging SSO --- tests/TestSegDist.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index f1e55ddd..cf002d15 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -44,6 +44,7 @@ int main(int argc, char **argv) // 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; int count = 0; From 772857aef15dc8416a0480a008289869ae2af502 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 4 Jun 2015 07:47:11 -0400 Subject: [PATCH 22/28] debugging SSO --- tests/TestSegDist.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index cf002d15..aa2773e4 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -57,9 +57,9 @@ int main(int argc, char **argv) for (k=1;k Date: Thu, 4 Jun 2015 07:49:03 -0400 Subject: [PATCH 23/28] debugging SSO --- tests/TestSegDist.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index aa2773e4..9b71f983 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -72,7 +72,7 @@ int main(int argc, char **argv) } } - DoubleArray Distance(Nx,Ny,Nz); + DoubleArray Distance(nx,ny,nz); // Initialize the signed distance function for (k=0;k Date: Thu, 4 Jun 2015 07:52:13 -0400 Subject: [PATCH 24/28] tests/TestSegDist seems to be working now --- common/Domain.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index c25ae3a5..c6245c43 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -813,8 +813,6 @@ inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ int count = 0; while (count < timesteps){ - printf("count=%i \n",count); - // Communicate the halo of values fillData.fill(Distance); From 4cc6bdba2aff8d608d601900892680e350ecfe5f Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 4 Jun 2015 08:11:46 -0400 Subject: [PATCH 25/28] Debugging unit test TestSegDist --- tests/TestSegDist.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index 9b71f983..ce3b883d 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -60,6 +60,7 @@ int main(int argc, char **argv) x = (nx-2)*Dm.iproc+i; y = (ny-2)*Dm.jproc+j; z = (nz-2)*Dm.kproc+k; + n = k*nx*ny+j*nx+i; // Initialize phase positions if ((x-nx+1)*(x-nx+1)+(y-ny+1)*(y-ny+1)+(z-nz+1)*(z-nz+1) < BubbleRadius*BubbleRadius){ @@ -89,7 +90,7 @@ int main(int argc, char **argv) if (rank==0) printf("Nz = %i \n",(int)Distance.size(2)); printf("Initialized! Converting to Signed Distance function \n"); - SSO(Distance,id,Dm,10); + SSO(Distance,id,Dm,1); char LocalRankFilename[40]; sprintf(LocalRankFilename,"Dist.%05i",rank); From 89353b8159c74b5d0949e5c59d9438ec75fc933b Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 4 Jun 2015 08:15:34 -0400 Subject: [PATCH 26/28] Debugging unit test TestSegDist --- tests/TestSegDist.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index ce3b883d..68bece76 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -47,6 +47,7 @@ int main(int argc, char **argv) Dm.CommInit(MPI_COMM_WORLD); nx+=2; ny+=2; nz+=2; + N = nx*ny*nz; int count = 0; char *id; From ba2043429232338432f0bfe99733005033d07b46 Mon Sep 17 00:00:00 2001 From: James McClure Date: Thu, 4 Jun 2015 08:28:26 -0400 Subject: [PATCH 27/28] Unit test for tests/TestSegDist works --- common/Domain.h | 2 +- tests/TestSegDist.cpp | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index c6245c43..e31c3459 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -788,7 +788,7 @@ inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ int Q=26; int q,i,j,k,n; - double dt=0.25; + 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; diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index 68bece76..841c4f42 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -52,7 +52,7 @@ int main(int argc, char **argv) char *id; id = new char [N]; - double BubbleRadius = 5; + double BubbleRadius = 25; // Initialize the bubble int x,y,z; for (k=1;k Date: Thu, 4 Jun 2015 09:46:49 -0400 Subject: [PATCH 28/28] Trying to figure out weird porosity bug in tests/BlobAnalyzeParallel --- tests/BlobAnalyzeParallel.cpp | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/tests/BlobAnalyzeParallel.cpp b/tests/BlobAnalyzeParallel.cpp index 8f11293a..cc86a15b 100644 --- a/tests/BlobAnalyzeParallel.cpp +++ b/tests/BlobAnalyzeParallel.cpp @@ -196,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);