From 72a4e0462d59725ad76ad04cec56224714e2ae16 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 22 Feb 2018 13:54:46 -0500 Subject: [PATCH] Added sphere test to project --- tests/CMakeLists.txt | 2 +- tests/GenerateSphereTest.cpp | 568 +++++++++++++++++++++++++++++++++++ 2 files changed, 569 insertions(+), 1 deletion(-) create mode 100644 tests/GenerateSphereTest.cpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 9202f034..86dece80 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -19,7 +19,7 @@ ADD_LBPM_EXECUTABLE( lbpm_captube_pp ) ADD_LBPM_EXECUTABLE( lbpm_inkbottle_pp ) ADD_LBPM_EXECUTABLE( lbpm_plates_pp ) ADD_LBPM_EXECUTABLE( lbpm_squaretube_pp ) -#ADD_LBPM_EXECUTABLE( TestBubble ) +ADD_LBPM_EXECUTABLE( GenerateSphereTest ) ADD_LBPM_EXECUTABLE( ComponentLabel ) ADD_LBPM_EXECUTABLE( ColorToBinary ) ADD_LBPM_EXECUTABLE( BlobAnalysis ) diff --git a/tests/GenerateSphereTest.cpp b/tests/GenerateSphereTest.cpp new file mode 100644 index 00000000..add99341 --- /dev/null +++ b/tests/GenerateSphereTest.cpp @@ -0,0 +1,568 @@ +#include +#include +#include +#include +#include +#include +#include + +//#include "common/pmmc.h" +#include "common/Domain.h" +#include "common/MPI_Helpers.h" +#include "common/Communication.h" + +/* + * Pre-Processor to generate signed distance function from sphere packing + * to use as an input domain for lattice Boltzmann simulator + * James E. McClure 2014 + */ + +using namespace std; + +inline void PackID(int *list, int count, char *sendbuf, char *ID){ + // Fill in the phase ID values from neighboring processors + // This packs up the values that need to be sent from one processor to another + int idx,n; + + for (idx=0; idx maxdist) maxdist=SignDist(i,j,k); + } + } + } + } + // total Global is the number of nodes in the pore-space + MPI_Allreduce(&count,&totalGlobal,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&maxdist,&maxdistGlobal,1,MPI_DOUBLE,MPI_MAX,Dm.Comm); + double volume=double(nprocx*nprocy*nprocz)*double(nx-2)*double(ny-2)*double(nz-2); + double porosity=totalGlobal/volume; + if (rank==0) printf("Media Porosity: %f \n",porosity); + if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal); + + + // Generate the NWP configuration + //if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit); + if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW); + // GenerateResidual(id,nx,ny,nz,Saturation); + + // Communication buffers + char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z; + char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ; + char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ; + char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z; + char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ; + char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ; + // send buffers + sendID_x = new char [Dm.sendCount_x]; + sendID_y = new char [Dm.sendCount_y]; + sendID_z = new char [Dm.sendCount_z]; + sendID_X = new char [Dm.sendCount_X]; + sendID_Y = new char [Dm.sendCount_Y]; + sendID_Z = new char [Dm.sendCount_Z]; + sendID_xy = new char [Dm.sendCount_xy]; + sendID_yz = new char [Dm.sendCount_yz]; + sendID_xz = new char [Dm.sendCount_xz]; + sendID_Xy = new char [Dm.sendCount_Xy]; + sendID_Yz = new char [Dm.sendCount_Yz]; + sendID_xZ = new char [Dm.sendCount_xZ]; + sendID_xY = new char [Dm.sendCount_xY]; + sendID_yZ = new char [Dm.sendCount_yZ]; + sendID_Xz = new char [Dm.sendCount_Xz]; + sendID_XY = new char [Dm.sendCount_XY]; + sendID_YZ = new char [Dm.sendCount_YZ]; + sendID_XZ = new char [Dm.sendCount_XZ]; + //...................................................................................... + // recv buffers + recvID_x = new char [Dm.recvCount_x]; + recvID_y = new char [Dm.recvCount_y]; + recvID_z = new char [Dm.recvCount_z]; + recvID_X = new char [Dm.recvCount_X]; + recvID_Y = new char [Dm.recvCount_Y]; + recvID_Z = new char [Dm.recvCount_Z]; + recvID_xy = new char [Dm.recvCount_xy]; + recvID_yz = new char [Dm.recvCount_yz]; + recvID_xz = new char [Dm.recvCount_xz]; + recvID_Xy = new char [Dm.recvCount_Xy]; + recvID_xZ = new char [Dm.recvCount_xZ]; + recvID_xY = new char [Dm.recvCount_xY]; + recvID_yZ = new char [Dm.recvCount_yZ]; + recvID_Yz = new char [Dm.recvCount_Yz]; + recvID_Xz = new char [Dm.recvCount_Xz]; + recvID_XY = new char [Dm.recvCount_XY]; + recvID_YZ = new char [Dm.recvCount_YZ]; + recvID_XZ = new char [Dm.recvCount_XZ]; + //...................................................................................... + int sendtag,recvtag; + sendtag = recvtag = 7; + + int x,y,z; + int ii,jj,kk; + int Nx = nx; + int Ny = ny; + int Nz = nz; + + double sw_old=1.0; + double sw_new=1.0; + double sw_diff_old = 1.0; + double sw_diff_new = 1.0; + + // Increase the critical radius until the target saturation is met + double deltaR=0.05; // amount to change the radius in voxel units + double Rcrit_old; + double Rcrit_new; + + double GlobalNumber = 1.f; + int imin,jmin,kmin,imax,jmax,kmax; + + Rcrit_new = maxdistGlobal; + while (sw_new > SW) + { + sw_diff_old = sw_diff_new; + sw_old = sw_new; + Rcrit_old = Rcrit_new; + Rcrit_new -= deltaR*Rcrit_old; + int Window=round(Rcrit_new); + if (Window == 0) Window = 1; // If Window = 0 at the begining, after the following process will have sw=1.0 + // and sw Rcrit_new){ + // loop over the window and update + imin=max(1,i-Window); + jmin=max(1,j-Window); + kmin=max(1,k-Window); + imax=min(Nx-1,i+Window); + jmax=min(Ny-1,j+Window); + kmax=min(Nz-1,k+Window); + for (kk=kmin; kk 1) + nspheres=atoi(argv[1]); + else nspheres=0; + + if (rank==0){ + //....................................................................... + // Reading the domain information file + //....................................................................... + ifstream domain("Domain.in"); + domain >> nprocx; + domain >> nprocy; + domain >> nprocz; + domain >> Nx; + domain >> Ny; + domain >> Nz; + domain >> Lx; + domain >> Ly; + domain >> Lz; + //....................................................................... + } + // ************************************************************** + // Broadcast simulation parameters from rank 0 to all other procs + MPI_Barrier(comm); + //................................................. + // Computational domain + MPI_Bcast(&Nx,1,MPI_INT,0,comm); + MPI_Bcast(&Ny,1,MPI_INT,0,comm); + MPI_Bcast(&Nz,1,MPI_INT,0,comm); + MPI_Bcast(&nprocx,1,MPI_INT,0,comm); + MPI_Bcast(&nprocy,1,MPI_INT,0,comm); + MPI_Bcast(&nprocz,1,MPI_INT,0,comm); + MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); + //................................................. + MPI_Barrier(comm); + + // ************************************************************** + + if (nprocs != nprocx*nprocy*nprocz){ + printf("nprocx = %i \n",nprocx); + printf("nprocy = %i \n",nprocy); + printf("nprocz = %i \n",nprocz); + INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); + } + + InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc, + rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z, + rank_xy, rank_XY, rank_xY, rank_Xy, rank_xz, rank_XZ, rank_xZ, rank_Xz, + rank_yz, rank_YZ, rank_yZ, rank_Yz ); + + MPI_Barrier(comm); + + int BoundaryCondition=1; + Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + + Nz += 2; + Nx = Ny = Nz; // Cubic domain + int N = Nx*Ny*Nz; + + // Define Dm.Communication sub-domain -- everywhere + for (int k=0; k 0.0){ + id[n] = 2; + } + // compute the porosity (actual interface location used) + if (SignDist(n) > 0.0){ + sum++; + } + } + } + } + sum_local = 1.0*sum; + MPI_Allreduce(&sum_local,&porosity,1,MPI_DOUBLE,MPI_SUM,comm); + porosity = porosity*iVol_global; + if (rank==0) printf("Media porosity = %f \n",porosity); + + // Run Morphological opening to initialize 50% saturation + double SW=0.50; + if (rank==0) printf("MorphOpen: Initializing with saturation %f \n",SW); + MorphOpen(SignDist, id, Dm, Nx, Ny, Nz, rank, SW); + + //......................................................... + // don't perform computations at the eight corners + id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0; + id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0; + //......................................................... + + //....................................................................... + sprintf(LocalRankString,"%05d",rank); + sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString); + FILE *DIST = fopen(LocalRankFilename,"wb"); + if (DIST==NULL) ERROR("Error opening file: ID.xxxxx"); + fwrite(SignDist.data(),1,N,DIST); + fclose(DIST); + //...................................................................... + + //....................................................................... + sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString); + FILE *IDFILE = fopen(LocalRankFilename,"wb"); + if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx"); + fwrite(id,1,N,IDFILE); + fclose(IDFILE); + //...................................................................... + } + // **************************************************** + MPI_Barrier(comm); + MPI_Finalize(); + // **************************************************** +}