diff --git a/tests/lbpm_morphopen_pp.cpp b/tests/lbpm_morphopen_pp.cpp index 6f054d9f..d38c3e82 100644 --- a/tests/lbpm_morphopen_pp.cpp +++ b/tests/lbpm_morphopen_pp.cpp @@ -50,473 +50,351 @@ int main(int argc, char **argv) // Initialize MPI int rank, nprocs; MPI_Init(&argc,&argv); - MPI_Comm comm = MPI_COMM_WORLD; + MPI_Comm comm = MPI_COMM_WORLD; MPI_Comm_rank(comm,&rank); MPI_Comm_size(comm,&nprocs); + { - //double Rcrit_new=1.f; // Hard-coded 'Rcrit' to avoid any calculations under resolutions. - double Rcrit_new=0.f; - double SW=strtod(argv[1],NULL); - if (rank==0){ - //printf("Critical radius %f (voxels)\n",Rcrit); - printf("Target saturation %f \n",SW); - } - - - // } - - //....................................................................... - // Reading the domain information file - //....................................................................... - int nprocx, nprocy, nprocz, nx, ny, nz, nspheres; - double Lx, Ly, Lz; - int i,j,k,n; - int BC=0; - - if (rank==0){ - ifstream domain("Domain.in"); - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> nx; - domain >> ny; - domain >> nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - - } - 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(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - //................................................. - MPI_Barrier(comm); - - // Check that the number of processors >= the number of ranks - if ( rank==0 ) { - printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz); - 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); - } - if ( nprocs < nprocx*nprocy*nprocz ){ - ERROR("Insufficient number of processors"); - } - - char LocalRankFilename[40]; - - int BoundaryCondition=1; - Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); - - nx+=2; ny+=2; nz+=2; - int N = nx*ny*nz; - char *id; - id = new char[N]; - - // Define communication sub-domain -- everywhere - for (int k=0; k maxdist) maxdist=SignDist(i,j,k); - } - } + char LocalRankString[8]; + char LocalRankFilename[40]; + + string filename; + double Rcrit_new, SW; + if (argc > 2){ + filename=argv[1]; + Rcrit_new=0.f; + SW=strtod(argv[2],NULL); } - } - // total Global is the number of nodes in the pore-space - MPI_Allreduce(&count,&totalGlobal,1,MPI_DOUBLE,MPI_SUM,comm); - MPI_Allreduce(&maxdist,&maxdistGlobal,1,MPI_DOUBLE,MPI_MAX,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);\ + else ERROR("No input database provided\n"); + // read the input database + auto db = std::make_shared( filename ); + auto domain_db = db->getDatabase( "Domain" ); - /* // Generate a histogram of pore size distribution - // Get all local pore sizes (local maxima) - if (rank==0) printf("Generating a histogram of pore sizes \n"); - int NumBins=100; - int *BinCounts; - BinCounts = new int [NumBins]; - int *GlobalHistogram; - GlobalHistogram = new int [NumBins]; - double GlobalValue; - double BinWidth,MinPoreSize,MaxPoreSize; - std::vector PoreSize; - for (int k=1; k 0.0){ - // Generate a list of all local maxima (each processor -- aggregate these later) - if ( SignDist(i,j,k) > SignDist(i+1,j,k) && SignDist(i,j,k) > SignDist(i-1,j,k) && - SignDist(i,j,k) > SignDist(i,j+1,k) && SignDist(i,j,k) > SignDist(i,j-1,k) && - SignDist(i,j,k) > SignDist(i,j,k+1) && SignDist(i,j,k) > SignDist(i,j,k-1) && - SignDist(i,j,k) > SignDist(i+1,j+1,k) && SignDist(i,j,k) > SignDist(i-1,j+1,k) && - SignDist(i,j,k) > SignDist(i+1,j-1,k) && SignDist(i,j,k) > SignDist(i-1,j-1,k) && - SignDist(i,j,k) > SignDist(i+1,j,k+1) && SignDist(i,j,k) > SignDist(i-1,j,k+1) && - SignDist(i,j,k) > SignDist(i+1,j,k-1) && SignDist(i,j,k) > SignDist(i-1,j,k-1) && - SignDist(i,j,k) > SignDist(i,j+1,k+1) && SignDist(i,j,k) > SignDist(i,j-1,k+1) && - SignDist(i,j,k) > SignDist(i,j+1,k-1) && SignDist(i,j,k) > SignDist(i,j-1,k-1) && - SignDist(i,j,k) > SignDist(i+1,j+1,k+1) && SignDist(i,j,k) > SignDist(i-1,j-1,k-1) && - SignDist(i,j,k) > SignDist(i+1,j-1,k+1) && SignDist(i,j,k) > SignDist(i-1,j+1,k-1) && - SignDist(i,j,k) > SignDist(i-1,j+1,k+1) && SignDist(i,j,k) > SignDist(i+1,j-1,k-1) && - SignDist(i,j,k) > SignDist(i+1,j+1,k-1) && SignDist(i,j,k) > SignDist(i-1,j-1,k+1)){ - // save the size of each pore - PoreSize.push_back(SignDist(i,j,k)); + // Read domain parameters + auto L = domain_db->getVector( "L" ); + auto size = domain_db->getVector( "n" ); + auto nproc = domain_db->getVector( "nproc" ); + auto ReadValues = domain_db->getVector( "ReadValues" ); + auto WriteValues = domain_db->getVector( "WriteValues" ); + + nx = size[0]; + ny = size[1]; + nz = size[2]; + nprocx = nproc[0]; + nprocy = nproc[1]; + nprocz = nproc[2]; + + int N = (nx+2)*(ny+2)*(nz+2); + + std::shared_ptr Dm (new Domain(domain_db,comm)); + // std::shared_ptr Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC)); + for (n=0; nid[n]=1; + Dm->CommInit(); + + DoubleArray SignDist(nx,ny,nz); + // Read the signed distance from file + sprintf(LocalRankFilename,"SignDist.%05i",rank); + FILE *DIST = fopen(LocalRankFilename,"rb"); + size_t ReadSignDist; + ReadSignDist=fread(SignDist.data(),8,N,DIST); + if (ReadSignDist != size_t(N)) printf("lbpm_morphopen_pp: Error reading signed distance function (rank=%i)\n",rank); + fclose(DIST); + + double count,countGlobal,totalGlobal; + count = 0.f; + double maxdist=0.f; + double maxdistGlobal; + for (int k=0; k maxdist) maxdist=SignDist(i,j,k); } } } } - } - // Compute min and max poresize - if (rank==0) printf(" computing local minimum and maximum... \n"); - MinPoreSize=MaxPoreSize=PoreSize[0]; - for (size_t idx=0; idx MaxPoreSize) MaxPoreSize=PoreSize[idx]; - } - // reduce to get global minimum and maximum - MPI_Allreduce(&MinPoreSize,&GlobalValue,1,MPI_DOUBLE,MPI_MIN,comm); - MinPoreSize=GlobalValue; - MPI_Allreduce(&MaxPoreSize,&GlobalValue,1,MPI_DOUBLE,MPI_MAX,comm); - MaxPoreSize=GlobalValue; - //if (rank==0) printf(" MaxPoreSize %f\n", MaxPoreSize); - //if (rank==0) printf(" MinPoreSize %f\n", MinPoreSize); - // Generate histogram counts - if (rank==0) printf(" generating local bin counts... \n"); - BinWidth=(MaxPoreSize-MinPoreSize)/double(NumBins); - for (size_t idx=0; idx2){ + Rcrit_new = strtod(argv[2],NULL); + if (rank==0) printf("Max. distance =%f, Initial critical radius = %f \n",maxdistGlobal,Rcrit_new); } - fclose(PORESIZE); - // Compute quartiles - double Q1,Q2,Q3,Q4; - int Qval=PoreCount/4; - Q1=Q2=Q3=MinPoreSize; - Q4=MaxPoreSize; - Count=0; - for (int idx=0; idx 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; kk2){ - Rcrit_new = strtod(argv[2],NULL); - if (rank==0) printf("Max. distance =%f, Initial critical radius = %f \n",maxdistGlobal,Rcrit_new); - } - 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