From 9098f725ec8eda235183cb6274a9c46e000b68d4 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 20 Feb 2019 17:59:28 -0500 Subject: [PATCH 01/13] initialized variables for inlet rule --- tests/lbpm_serial_decomp.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/lbpm_serial_decomp.cpp b/tests/lbpm_serial_decomp.cpp index bcdd9fc8..738e387e 100644 --- a/tests/lbpm_serial_decomp.cpp +++ b/tests/lbpm_serial_decomp.cpp @@ -58,6 +58,10 @@ int main(int argc, char **argv) // char fluidValue,solidValue; xStart=yStart=zStart=0; + inlet_count_x = 0; + inlet_count_y = 0; + inlet_count_z = 0; + checkerSize = 32; // read the input database auto db = std::make_shared( filename ); auto domain_db = db->getDatabase( "Domain" ); From 5ccb06646d39106c02dc40018a80a0e33cade52b Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 20 Feb 2019 18:33:41 -0500 Subject: [PATCH 02/13] refactor morphdrain pp --- tests/lbpm_morphdrain_pp.cpp | 368 +++++++---------------------------- 1 file changed, 73 insertions(+), 295 deletions(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index 70798e88..3ca4c7a9 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -134,6 +134,8 @@ int main(int argc, char **argv) // Initialize the domain and communication Array id_solid(nx,ny,nz); DoubleArray SignDist(nx,ny,nz); + DoubleArray phase(nx,ny,nz); + IntArray phase_label(nx,ny,nz); // Solve for the position of the solid phase for (int k=0;k 1){ id[n] = 2; - count+=1.0; } } } } - // total Global is the number of nodes in the pore-space - MPI_Allreduce(&count,&totalGlobal,1,MPI_DOUBLE,MPI_SUM,comm); - double porosity=totalGlobal/(double(nprocx*nprocy*nprocz)*double(nx-2)*double(ny-2)*double(nz-2)); - if (rank==0) printf("Media Porosity: %f \n",porosity); - - double radius,Rcrit_new; - radius = 0.0; - // Layer the inlet with NWP - if (kproc == 0){ - for(j=0; j radius){ - radius=SignDist(i,j,0); - } + // calculate distance to non-wetting fluid + if (domain_db->keyExists( "HistoryLabels" )){ + if (rank==0) printf("Relabel solid components that touch fluid 1 \n"); + auto LabelList = domain_db->getVector( "ComponentLabels" ); + auto HistoryLabels = domain_db->getVector( "HistoryLabels" ); + size_t NLABELS=LabelList.size(); + if (rank==0){ + for (unsigned int idx=0; idx < NLABELS; idx++){ + char VALUE = LabelList[idx]; + char NEWVAL = HistoryLabels[idx]; + printf(" Relabel component %d as %d \n", VALUE, NEWVAL); } } - } - - MPI_Allreduce(&radius,&Rcrit_new,1,MPI_DOUBLE,MPI_MAX,comm); - - if (rank==0) printf("Starting morhpological drainage with critical radius = %f \n",Rcrit_new); - - int imin,jmin,kmin,imax,jmax,kmax; - - // Decrease the critical radius until the target saturation is met - double deltaR=0.01; // amount to change the radius in voxel units - double Rcrit_old; - double sw_old=1.0; // initial WP saturation set to one - double sw_new=1.0; // initial WP saturation set to one - double sw_diff_old = 1.0; - double sw_diff_new = 1.0; - - while (sw_new > SW && Rcrit_new > 0.99){ - - Rcrit_old = Rcrit_new; - Rcrit_new -= deltaR*Rcrit_old;// decrease critical radius - sw_old = sw_new; - sw_diff_old = sw_diff_new; - - int Window=round(Rcrit_new); - GlobalNumber = 1.0; - - while (GlobalNumber > 0){ - - //if (rank==0) printf("GlobalNumber=%f \n",GlobalNumber); - double LocalNumber=GlobalNumber=0.f; - for(k=0; k 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 Date: Wed, 20 Feb 2019 18:37:06 -0500 Subject: [PATCH 03/13] refactor morphdrain pp --- tests/lbpm_morphdrain_pp.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index 3ca4c7a9..b08e4d7c 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -12,6 +12,8 @@ #include "common/Array.h" #include "common/Domain.h" #include "analysis/distance.h" +#include "analysis/morphology.h" +#include "analysis/runAnalysis.h" //************************************************************************* // Morpohological drainage pre-processor @@ -181,7 +183,7 @@ int main(int argc, char **argv) // Extract only the connected part BlobIDstruct new_index; double vF=0.0; double vS=0.0; - ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm.rank_info,phase,SignDist,vF,vS,phase_label,Dm.comm); + ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm.rank_info,phase,SignDist,vF,vS,phase_label,Dm.Comm); MPI_Barrier(comm); for (int k=0;k Date: Wed, 20 Feb 2019 18:37:49 -0500 Subject: [PATCH 04/13] refactor morphdrain pp --- tests/lbpm_morphdrain_pp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index b08e4d7c..c65f2697 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -165,7 +165,7 @@ int main(int argc, char **argv) CalcDist(SignDist,id_solid,Dm); // Run the morphological opening - MorphOpen(SignDist, id, Dm, SW); + MorphOpen(SignDist, id, *Dm, SW); for (int k=0;k Date: Wed, 20 Feb 2019 18:40:44 -0500 Subject: [PATCH 05/13] refactor morphdrain pp --- tests/lbpm_morphdrain_pp.cpp | 350 ++++++++++++++++------------------- 1 file changed, 159 insertions(+), 191 deletions(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index c65f2697..49e04b7b 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -16,207 +16,99 @@ #include "analysis/runAnalysis.h" //************************************************************************* -// Morpohological drainage pre-processor -// Generate states on primary drainage using morphological approach +// Morpohologica pre-processor +// Initialize phase distribution using morphological approach // Signed distance function is used to determine fluid configuration //************************************************************************* -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 1){ - filename=argv[1]; - } - else ERROR("No input database provided\n"); - // read the input database - auto db = std::make_shared( filename ); - auto domain_db = db->getDatabase( "Domain" ); - - // Read domain parameters - auto L = domain_db->getVector( "L" ); - auto size = domain_db->getVector( "n" ); - auto nproc = domain_db->getVector( "nproc" ); - double SW = domain_db->getScalar( "Sw" ); - - nx = size[0]; - ny = size[1]; - nz = size[2]; - nprocx = nproc[0]; - nprocy = nproc[1]; - nprocz = nproc[2]; - - // 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(domain_db,comm); - - 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 1){ + filename=argv[1]; + Rcrit_new=0.f; + //SW=strtod(argv[2],NULL); } - } + else ERROR("No input database provided\n"); + // read the input database + auto db = std::make_shared( filename ); + auto domain_db = db->getDatabase( "Domain" ); - Dm.CommInit(); + // 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" ); + SW = domain_db->getScalar("Sw"); - sprintf(LocalRankFilename,"ID.%05i",rank); - size_t readID; - FILE *IDFILE = fopen(LocalRankFilename,"rb"); - if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx"); - readID=fread(id,1,N,IDFILE); - if (readID != size_t(N)) printf("lbpm_morphdrain_pp: Error reading ID (rank=%i) \n",rank); - fclose(IDFILE); + // 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); + + nx = size[0]; + ny = size[1]; + nz = size[2]; + nprocx = nproc[0]; + nprocy = nproc[1]; + nprocz = nproc[2]; - int xdim,ydim,zdim; - xdim=Dm.Nx-2; - ydim=Dm.Ny-2; - zdim=Dm.Nz-2; - fillHalo fillData(Dm.Comm, Dm.rank_info,{xdim,ydim,zdim},{1,1,1},0,1); - - // Generate the signed distance map - // Initialize the domain and communication - Array id_solid(nx,ny,nz); - DoubleArray SignDist(nx,ny,nz); - DoubleArray phase(nx,ny,nz); - IntArray phase_label(nx,ny,nz); + int N = (nx+2)*(ny+2)*(nz+2); - // Solve for the position of the solid phase - for (int k=0;k 0) id_solid(i,j,k) = 1; - else id_solid(i,j,k) = 0; - } - } - } - // Initialize the signed distance function - for (int k=0;k 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(); - if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); - CalcDist(SignDist,id_solid,Dm); + char *id; + id = new char [N]; + sprintf(LocalRankFilename,"ID.%05i",rank); + size_t readID; + FILE *IDFILE = fopen(LocalRankFilename,"rb"); + if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx"); + readID=fread(id,1,N,IDFILE); + if (readID != size_t(N)) printf("lbpm_morphopen_pp: Error reading ID (rank=%i) \n",rank); + fclose(IDFILE); - // Run the morphological opening - MorphOpen(SignDist, id, *Dm, SW); + nx+=2; ny+=2; nz+=2; + // Generate the signed distance map + // Initialize the domain and communication + Array id_solid(nx,ny,nz); + DoubleArray SignDist(nx,ny,nz); + DoubleArray phase(nx,ny,nz); + IntArray phase_label(nx,ny,nz); - for (int k=0;k 1){ - id[n] = 2; - } - } - } - } - - // calculate distance to non-wetting fluid - if (domain_db->keyExists( "HistoryLabels" )){ - if (rank==0) printf("Relabel solid components that touch fluid 1 \n"); - auto LabelList = domain_db->getVector( "ComponentLabels" ); - auto HistoryLabels = domain_db->getVector( "HistoryLabels" ); - size_t NLABELS=LabelList.size(); - if (rank==0){ - for (unsigned int idx=0; idx < NLABELS; idx++){ - char VALUE = LabelList[idx]; - char NEWVAL = HistoryLabels[idx]; - printf(" Relabel component %d as %d \n", VALUE, NEWVAL); - } - } + // Solve for the position of the solid phase for (int k=0;k 0) id_solid(i,j,k) = 1; + else id_solid(i,j,k) = 0; } } } @@ -230,32 +122,108 @@ int main(int argc, char **argv) } } } - CalcDist(SignDist,id_solid,Dm); - // re-label IDs near the non-wetting fluid + + if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); + CalcDist(SignDist,id_solid,*Dm); + + MPI_Barrier(comm); + + // Run the morphological opening + MorphOpen(SignDist, id, Dm, SW); + for (int k=0;k 1){ + id[n] = 2; + } + } + } + } + + + // calculate distance to non-wetting fluid + if (domain_db->keyExists( "HistoryLabels" )){ + if (rank==0) printf("Relabel solid components that touch fluid 1 \n"); + auto LabelList = domain_db->getVector( "ComponentLabels" ); + auto HistoryLabels = domain_db->getVector( "HistoryLabels" ); + size_t NLABELS=LabelList.size(); + if (rank==0){ + for (unsigned int idx=0; idx < NLABELS; idx++){ + char VALUE = LabelList[idx]; + char NEWVAL = HistoryLabels[idx]; + printf(" Relabel component %d as %d \n", VALUE, NEWVAL); + } + } + for (int k=0;k Date: Wed, 20 Feb 2019 18:41:36 -0500 Subject: [PATCH 06/13] refactor morphdrain pp --- tests/lbpm_morphdrain_pp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index 49e04b7b..4184e4ba 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -147,7 +147,7 @@ int main(int argc, char **argv) // Extract only the connected part BlobIDstruct new_index; double vF=0.0; double vS=0.0; - ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm.rank_info,phase,SignDist,vF,vS,phase_label,Dm.Comm); + ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm); MPI_Barrier(comm); for (int k=0;k Date: Thu, 21 Feb 2019 10:06:29 -0500 Subject: [PATCH 07/13] refactor morphdrain --- analysis/morphology.cpp | 311 +++++++++++++++++++++++++++++++++++ analysis/morphology.h | 2 + tests/lbpm_morphdrain_pp.cpp | 37 +---- 3 files changed, 315 insertions(+), 35 deletions(-) diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index c15620c2..64d6ff62 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -326,6 +326,317 @@ double morph_open() } */ +//*************************************************************************************** +double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, double VoidFraction){ + // SignDist is the distance to the object that you want to constaing the morphological opening + // VoidFraction is the the empty space where the object inst + // id is a labeled map + // Dm contains information about the domain structure + + int nx = Dm->Nx; + int ny = Dm->Ny; + 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(); + + DoubleArray phase(Nx,Ny,Nz); + IntArray phase_label(Nx,Ny,Nz); + + int n; + double final_void_fraction; + double count,countGlobal,totalGlobal; + count = 0.f; + double maxdist=-200.f; + double maxdistGlobal; + for (int k=1; k maxdist) maxdist=SignDist(i,j,k); + if ( SignDist(i,j,k) > 0.0 ){ + count += 1.0; + id[n] = 2; + } + } + } + } + + MPI_Barrier(Dm->Comm); + + // 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 volume_fraction=totalGlobal/volume; + if (rank==0) printf("Volume fraction for morphological opening: %f \n",volume_fraction); + if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal); + + // 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 void_fraction_old=1.0; + double void_fraction_new=1.0; + double void_fraction_diff_old = 1.0; + double void_fraction_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 GlobalNumber = 1.f; + int imin,jmin,kmin,imax,jmax,kmax; + + double Rcrit_new = maxdistGlobal; + //if (argc>2){ + // Rcrit_new = strtod(argv[2],NULL); + // if (rank==0) printf("Max. distance =%f, Initial critical radius = %f \n",maxdistGlobal,Rcrit_new); + //} + + while (void_fraction_new > VoidFraction) + { + void_fraction_diff_old = void_fraction_diff_new; + void_fraction_old = void_fraction_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; kksendList_x, Dm->sendCount_x ,sendID_x, id); + PackID(Dm->sendList_X, Dm->sendCount_X ,sendID_X, id); + PackID(Dm->sendList_y, Dm->sendCount_y ,sendID_y, id); + PackID(Dm->sendList_Y, Dm->sendCount_Y ,sendID_Y, id); + PackID(Dm->sendList_z, Dm->sendCount_z ,sendID_z, id); + PackID(Dm->sendList_Z, Dm->sendCount_Z ,sendID_Z, id); + PackID(Dm->sendList_xy, Dm->sendCount_xy ,sendID_xy, id); + PackID(Dm->sendList_Xy, Dm->sendCount_Xy ,sendID_Xy, id); + PackID(Dm->sendList_xY, Dm->sendCount_xY ,sendID_xY, id); + PackID(Dm->sendList_XY, Dm->sendCount_XY ,sendID_XY, id); + PackID(Dm->sendList_xz, Dm->sendCount_xz ,sendID_xz, id); + PackID(Dm->sendList_Xz, Dm->sendCount_Xz ,sendID_Xz, id); + PackID(Dm->sendList_xZ, Dm->sendCount_xZ ,sendID_xZ, id); + PackID(Dm->sendList_XZ, Dm->sendCount_XZ ,sendID_XZ, id); + PackID(Dm->sendList_yz, Dm->sendCount_yz ,sendID_yz, id); + PackID(Dm->sendList_Yz, Dm->sendCount_Yz ,sendID_Yz, id); + PackID(Dm->sendList_yZ, Dm->sendCount_yZ ,sendID_yZ, id); + PackID(Dm->sendList_YZ, Dm->sendCount_YZ ,sendID_YZ, id); + //...................................................................................... + MPI_Sendrecv(sendID_x,Dm->sendCount_x,MPI_CHAR,Dm->rank_x(),sendtag, + recvID_X,Dm->recvCount_X,MPI_CHAR,Dm->rank_X(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_X,Dm->sendCount_X,MPI_CHAR,Dm->rank_X(),sendtag, + recvID_x,Dm->recvCount_x,MPI_CHAR,Dm->rank_x(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_y,Dm->sendCount_y,MPI_CHAR,Dm->rank_y(),sendtag, + recvID_Y,Dm->recvCount_Y,MPI_CHAR,Dm->rank_Y(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_Y,Dm->sendCount_Y,MPI_CHAR,Dm->rank_Y(),sendtag, + recvID_y,Dm->recvCount_y,MPI_CHAR,Dm->rank_y(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_z,Dm->sendCount_z,MPI_CHAR,Dm->rank_z(),sendtag, + recvID_Z,Dm->recvCount_Z,MPI_CHAR,Dm->rank_Z(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_Z,Dm->sendCount_Z,MPI_CHAR,Dm->rank_Z(),sendtag, + recvID_z,Dm->recvCount_z,MPI_CHAR,Dm->rank_z(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_xy,Dm->sendCount_xy,MPI_CHAR,Dm->rank_xy(),sendtag, + recvID_XY,Dm->recvCount_XY,MPI_CHAR,Dm->rank_XY(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_XY,Dm->sendCount_XY,MPI_CHAR,Dm->rank_XY(),sendtag, + recvID_xy,Dm->recvCount_xy,MPI_CHAR,Dm->rank_xy(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_Xy,Dm->sendCount_Xy,MPI_CHAR,Dm->rank_Xy(),sendtag, + recvID_xY,Dm->recvCount_xY,MPI_CHAR,Dm->rank_xY(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_xY,Dm->sendCount_xY,MPI_CHAR,Dm->rank_xY(),sendtag, + recvID_Xy,Dm->recvCount_Xy,MPI_CHAR,Dm->rank_Xy(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_xz,Dm->sendCount_xz,MPI_CHAR,Dm->rank_xz(),sendtag, + recvID_XZ,Dm->recvCount_XZ,MPI_CHAR,Dm->rank_XZ(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_XZ,Dm->sendCount_XZ,MPI_CHAR,Dm->rank_XZ(),sendtag, + recvID_xz,Dm->recvCount_xz,MPI_CHAR,Dm->rank_xz(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_Xz,Dm->sendCount_Xz,MPI_CHAR,Dm->rank_Xz(),sendtag, + recvID_xZ,Dm->recvCount_xZ,MPI_CHAR,Dm->rank_xZ(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_xZ,Dm->sendCount_xZ,MPI_CHAR,Dm->rank_xZ(),sendtag, + recvID_Xz,Dm->recvCount_Xz,MPI_CHAR,Dm->rank_Xz(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_yz,Dm->sendCount_yz,MPI_CHAR,Dm->rank_yz(),sendtag, + recvID_YZ,Dm->recvCount_YZ,MPI_CHAR,Dm->rank_YZ(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_YZ,Dm->sendCount_YZ,MPI_CHAR,Dm->rank_YZ(),sendtag, + recvID_yz,Dm->recvCount_yz,MPI_CHAR,Dm->rank_yz(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_Yz,Dm->sendCount_Yz,MPI_CHAR,Dm->rank_Yz(),sendtag, + recvID_yZ,Dm->recvCount_yZ,MPI_CHAR,Dm->rank_yZ(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_yZ,Dm->sendCount_yZ,MPI_CHAR,Dm->rank_yZ(),sendtag, + recvID_Yz,Dm->recvCount_Yz,MPI_CHAR,Dm->rank_Yz(),recvtag,Dm->Comm,MPI_STATUS_IGNORE); + //...................................................................................... + UnpackID(Dm->recvList_x, Dm->recvCount_x ,recvID_x, id); + UnpackID(Dm->recvList_X, Dm->recvCount_X ,recvID_X, id); + UnpackID(Dm->recvList_y, Dm->recvCount_y ,recvID_y, id); + UnpackID(Dm->recvList_Y, Dm->recvCount_Y ,recvID_Y, id); + UnpackID(Dm->recvList_z, Dm->recvCount_z ,recvID_z, id); + UnpackID(Dm->recvList_Z, Dm->recvCount_Z ,recvID_Z, id); + UnpackID(Dm->recvList_xy, Dm->recvCount_xy ,recvID_xy, id); + UnpackID(Dm->recvList_Xy, Dm->recvCount_Xy ,recvID_Xy, id); + UnpackID(Dm->recvList_xY, Dm->recvCount_xY ,recvID_xY, id); + UnpackID(Dm->recvList_XY, Dm->recvCount_XY ,recvID_XY, id); + UnpackID(Dm->recvList_xz, Dm->recvCount_xz ,recvID_xz, id); + UnpackID(Dm->recvList_Xz, Dm->recvCount_Xz ,recvID_Xz, id); + UnpackID(Dm->recvList_xZ, Dm->recvCount_xZ ,recvID_xZ, id); + UnpackID(Dm->recvList_XZ, Dm->recvCount_XZ ,recvID_XZ, id); + UnpackID(Dm->recvList_yz, Dm->recvCount_yz ,recvID_yz, id); + UnpackID(Dm->recvList_Yz, Dm->recvCount_Yz ,recvID_Yz, id); + UnpackID(Dm->recvList_yZ, Dm->recvCount_yZ ,recvID_yZ, id); + UnpackID(Dm->recvList_YZ, Dm->recvCount_YZ ,recvID_YZ, id); + //...................................................................................... + MPI_Allreduce(&LocalNumber,&GlobalNumber,1,MPI_DOUBLE,MPI_SUM,Dm->Comm); + + + for (int k=0;krank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm); + MPI_Barrier(comm); + + for (int k=0;k 1){ + id[n] = 2; + } + } + } + } + + count = 0.f; + for (int k=1; kComm); + void_fraction_new = countGlobal/totalGlobal; + void_fraction_diff_new = abs(void_fraction_new-VoidFraction); + if (rank==0){ + printf(" %f ",void_fraction_new); + printf(" %f\n",Rcrit_new); + } + } + + if (void_fraction_diff_new &id, std::shared_ptr Dm, double TargetGrowth){ int Nx = Dm->Nx; diff --git a/analysis/morphology.h b/analysis/morphology.h index a543d6e3..4692475d 100644 --- a/analysis/morphology.h +++ b/analysis/morphology.h @@ -1,6 +1,8 @@ // Morphological opening routine #include "common/Array.h" #include "common/Domain.h" +#include "analysis/runAnalysis.h" double MorphOpen(DoubleArray &SignDist, char *id, std::shared_ptr Dm, double VoidFraction); +double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, double VoidFraction); double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array &id, std::shared_ptr Dm, double TargetVol); \ No newline at end of file diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index 4184e4ba..9f655f37 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -98,8 +98,6 @@ int main(int argc, char **argv) // Initialize the domain and communication Array id_solid(nx,ny,nz); DoubleArray SignDist(nx,ny,nz); - DoubleArray phase(nx,ny,nz); - IntArray phase_label(nx,ny,nz); // Solve for the position of the solid phase for (int k=0;krank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm); - MPI_Barrier(comm); - - for (int k=0;k 1){ - id[n] = 2; - } - } - } - } - - + MorphDrain(SignDist, id, Dm, SW); + // calculate distance to non-wetting fluid if (domain_db->keyExists( "HistoryLabels" )){ if (rank==0) printf("Relabel solid components that touch fluid 1 \n"); From 3f4aa7f4c0d6550a6c1ca3b9d7b6fe521889b91f Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 21 Feb 2019 10:08:35 -0500 Subject: [PATCH 08/13] refactor morphdrain --- analysis/morphology.cpp | 31 +++++++++++++++---------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 64d6ff62..414fa14f 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -344,8 +344,8 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, d int nprocz = Dm->nprocz(); int rank = Dm->rank(); - DoubleArray phase(Nx,Ny,Nz); - IntArray phase_label(Nx,Ny,Nz); + DoubleArray phase(nx,ny,nz); + IntArray phase_label(nx,ny,nz); int n; double final_void_fraction; @@ -569,11 +569,10 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, d //...................................................................................... MPI_Allreduce(&LocalNumber,&GlobalNumber,1,MPI_DOUBLE,MPI_SUM,Dm->Comm); - - for (int k=0;k Dm, d // Extract only the connected part BlobIDstruct new_index; double vF=0.0; double vS=0.0; - ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,Dm->rank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm); + ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm); MPI_Barrier(comm); - for (int k=0;k 1){ id[n] = 2; } @@ -601,10 +600,10 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, d } count = 0.f; - for (int k=1; k Date: Thu, 21 Feb 2019 10:09:31 -0500 Subject: [PATCH 09/13] refactor morphdrain --- analysis/morphology.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 414fa14f..3deb6b97 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -586,7 +586,7 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, d BlobIDstruct new_index; double vF=0.0; double vS=0.0; ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm); - MPI_Barrier(comm); + MPI_Barrier(Dm->Comm); for (int k=1; k Date: Thu, 21 Feb 2019 13:52:47 -0500 Subject: [PATCH 10/13] analyze WP connectivity in morphdrain --- analysis/morphology.cpp | 48 ++++++++++++++++++++++++++++++++++-- tests/lbpm_morphdrain_pp.cpp | 1 - 2 files changed, 46 insertions(+), 3 deletions(-) diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 3deb6b97..33c193ae 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -582,7 +582,7 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, d } } - // Extract only the connected part + // Extract only the connected part of NWP BlobIDstruct new_index; double vF=0.0; double vS=0.0; ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm); @@ -598,13 +598,45 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, d } } } + + + /* + * Extract only the connected part of NWP + */ + for (int k=1; krank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm); + MPI_Barrier(Dm->Comm); + + for (int k=1; k 1){ + id[n] = 20; + } + } + } + } + // done count = 0.f; for (int k=1; k 1){ count+=1.0; } } @@ -633,6 +665,18 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, d printf("Final critical radius=%f\n",Rcrit_old); } } + // label all WP components as 2 + for (int k=1; k 1){ + id[n] = 2; + } + } + } + } + return final_void_fraction; } diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index 9f655f37..c9a58222 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -13,7 +13,6 @@ #include "common/Domain.h" #include "analysis/distance.h" #include "analysis/morphology.h" -#include "analysis/runAnalysis.h" //************************************************************************* // Morpohologica pre-processor From bc1a152a1a7a0a248234c41b8e819c1c8c901970 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 21 Feb 2019 13:56:19 -0500 Subject: [PATCH 11/13] comment analyze WP connectivity in morphdrain --- analysis/morphology.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 33c193ae..158e53eb 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -602,7 +602,6 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, d /* * Extract only the connected part of NWP - */ for (int k=1; k Dm, d } } // done + */ count = 0.f; for (int k=1; k Date: Thu, 21 Feb 2019 18:15:40 -0500 Subject: [PATCH 12/13] combine solid +WP when analyzing WP connectivity in morphdrain --- analysis/morphology.cpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 158e53eb..8a0ffed5 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -602,6 +602,7 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, d /* * Extract only the connected part of NWP + */ for (int k=1; k Dm, d if (id[n] == 2){ phase(i,j,k) = 1.0; } - else + else if (id[n] == 1){ + // nwp phase(i,j,k) = -1.0; + } + else{ + // treat solid as WP since films can connect + phase(i,j,k) = 1.0; + } } } } @@ -629,7 +636,6 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, d } } // done - */ count = 0.f; for (int k=1; k Date: Thu, 21 Feb 2019 18:18:31 -0500 Subject: [PATCH 13/13] comment analyze WP connectivity in morphdrain --- analysis/morphology.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 8a0ffed5..e7bac352 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -602,7 +602,6 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, d /* * Extract only the connected part of NWP - */ for (int k=1; k Dm, d } } // done + */ count = 0.f; for (int k=1; k