OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Simpler morph processor

Browse Source
This commit is contained in:
James E McClure
2017-08-16 15:47:20 -04:00
parent 04ce3bbb68
commit 0342f0acd3
2 changed files with 202 additions and 247 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
tests/CMakeLists.txt
View File

@@ -8,6 +8,7 @@ ADD_LBPM_EXECUTABLE( lbpm_random_pp )
ADD_LBPM_EXECUTABLE( lbpm_refine_pp )
ADD_LBPM_EXECUTABLE( lbpm_morphdrain_pp )
ADD_LBPM_EXECUTABLE( lbpm_morphopen_pp )
ADD_LBPM_EXECUTABLE( lbpm_morph_pp )
ADD_LBPM_EXECUTABLE( lbpm_segmented_pp )
ADD_LBPM_EXECUTABLE( lbpm_segmented_decomp )
ADD_LBPM_EXECUTABLE( lbpm_serial_decomp )

254
tests/lbpm_morph_pp.cpp
View File

@@ -14,8 +14,8 @@
//*************************************************************************
// Morpohologica pre-processor
// Initialize phase distribution using morphological approach
// Morpohological drainage pre-processor
// Generate states on primary drainage using morphological approach
// Signed distance function is used to determine fluid configuration
//*************************************************************************
inline void PackID(int *list, int count, char *sendbuf, char *ID){
@@ -44,7 +44,6 @@ inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
//***************************************************************************************
int main(int argc, char **argv)
{
// Initialize MPI
@@ -54,10 +53,11 @@ int main(int argc, char **argv)
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
double Rcrit;
Rcrit=strtod(argv[1],NULL);
//double Rcrit;
double Rcrit=strtod(argv[1],NULL);
if (rank==0){
printf("Critical radius %f (voxels)\n",Rcrit);
//printf("Target saturation %f \n",SW);
}
// }
@@ -127,6 +127,26 @@ int main(int argc, char **argv)
}
}
}
/*// Exclude the maximum / minimum z rows
if (rank/(nprocx*nprocy)==0){
int k=0;
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n = k*nx*ny+j*nx+i;
Dm.id[n] = 0;
}
}
}
if (rank/(nprocx*nprocy)==nprocz-1){
int k=nz-1;
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n = k*nx*ny+j*nx+i;
Dm.id[n] = 0;
}
}
}
*/
Dm.CommInit(comm);
DoubleArray SignDist(nx,ny,nz);
@@ -138,137 +158,14 @@ int main(int argc, char **argv)
if (ReadSignDist != size_t(N)) printf("lbpm_morphdrain_pp: Error reading signed distance function (rank=%i)\n",rank);
fclose(DIST);
int count,countGlobal,totalGlobal;
count = 0;
double maxdist=0;
double maxdistGlobal;
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
n = k*nx*ny+j*nx+i;
if (SignDist(i,j,k) < 0.0) id[n] = 0;
else{
// initially saturated with wetting phase
id[n] = 2;
count++;
// extract maximum distance for critical radius
if ( SignDist(i,j,k) > maxdist) maxdist=SignDist(i,j,k);
}
}
}
}
// total Global is the number of nodes in the pore-space
MPI_Allreduce(&count,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
MPI_Allreduce(&maxdist,&maxdistGlobal,1,MPI_DOUBLE,MPI_MAX,comm);
float porosity=float(totalGlobal)/(nprocx*nprocy*nprocz*(nx-2)*(ny-2)*(nz-2));
if (rank==0) printf("Media Porosity: %f \n",porosity);
if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal);\
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_segmented_pp: Error reading ID (rank=%i) \n",rank);
fclose(IDFILE);
// 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<double> PoreSize;
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
n = k*nx*ny+j*nx+i;
if (SignDist(i,j,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));
}
}
}
}
}
// Compute min and max poresize
if (rank==0) printf(" computing local minimum and maximum... \n");
MinPoreSize=MaxPoreSize=PoreSize[0];
for (int idx=0; idx<PoreSize.size(); idx++){
if (PoreSize[idx] < MinPoreSize) MinPoreSize=PoreSize[idx];
if (PoreSize[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 (int idx=0; idx<PoreSize.size(); idx++){
double value = PoreSize[idx];
int myBin = 0;
while (MinPoreSize+myBin*BinWidth < value) myBin++;
//int((value-MinPoreSize)/double(BinWidth));
BinCounts[myBin]+=1;
}
if (rank==0) printf(" summing global bin counts... \n");
// Reduce the counts to generate the fhistogram at rank=0
MPI_Reduce(BinCounts,GlobalHistogram,NumBins,MPI_INT,MPI_SUM,0,comm);
FILE *PORESIZE;
if (rank==0){
PORESIZE=fopen("PoreSize.hist","w");
printf(" writing PoreSize.hist \n");
uint64_t PoreCount=0;
uint64_t Count;
double PoreVol=0.f;
for (int idx=0; idx<NumBins; idx++){
double BinCenter=MinPoreSize+idx*BinWidth;
Count=GlobalHistogram[idx];
PoreCount+=Count;
PoreVol+=Count*BinCenter*BinCenter*BinCenter;
fprintf(PORESIZE,"%lu %f\n",Count,BinCenter);
}
fclose(PORESIZE);
// Compute quartiles
//printf("Total pores: %lu\n",PoreCount);
double Q1,Q2,Q3,Q4;
//uint64_t Qval=PoreCount/4;
double Qval=PoreVol*0.25;
Q1=Q2=Q3=MinPoreSize;
Q4=MaxPoreSize;
//printf("Volume per quartile %f\n",Qval);
PoreVol=0.f;
for (int idx=0; idx<NumBins; idx++){
double BinCenter=MinPoreSize+idx*BinWidth;
Count=GlobalHistogram[idx];
PoreVol+=Count*BinCenter*BinCenter*BinCenter;
if (PoreVol<Qval) Q1+=BinWidth;
if (PoreVol<2*Qval) Q2+=BinWidth;
if (PoreVol<3*Qval) Q3+=BinWidth;
}
printf("Quartiles (volumetric) for pore size distribution \n");
printf("Q1 %f\n",Q1);
printf("Q2 %f\n",Q2);
printf("Q3 %f\n",Q3);
printf("Q4 %f\n",Q4);
}
MPI_Barrier(comm);
Dm.CommInit(comm);
int iproc = Dm.iproc;
@@ -276,7 +173,8 @@ int main(int argc, char **argv)
int kproc = Dm.kproc;
// Generate the NWP configuration
if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
//if (rank==0) printf("Performing morphological drainage with critical radius %f \n", Rcrit);
//if (rank==0) printf("Performing morphological drainage with target saturation %f \n", SW);
// GenerateResidual(id,nx,ny,nz,Saturation);
// Communication buffers
@@ -334,18 +232,46 @@ int main(int argc, char **argv)
int Nx = nx;
int Ny = ny;
int Nz = nz;
float sat = 0.f;
int GlobalNumber = 1;
int Window=int(Rcrit);
double GlobalNumber = 1.f;
double count,countGlobal,totalGlobal;
count = 0.f;
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n = k*nx*ny+j*nx+i;
if (SignDist(i,j,k) < 0.0) id[n] = 0;
else{
// initially saturated with wetting phase
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);
if (rank==0) printf("Starting morhpological drainage with critical radius = %f \n",Rcrit);
int imin,jmin,kmin,imax,jmax,kmax;
int LocalNumber=0;
int Window=round(Rcrit);
GlobalNumber = 1.0;
while (GlobalNumber > 0){
//if (rank==0) printf("GlobalNumber=%f \n",GlobalNumber);
double LocalNumber=GlobalNumber=0.f;
for(k=0; k<Nz; k++){
for(j=0; j<Ny; j++){
for(i=0; i<Nx; i++){
n = k*nx*ny + j*nx+i;
if (SignDist(i,j,k) > Rcrit){
if (id[n] == 1 && SignDist(i,j,k) > Rcrit){
// loop over the window and update
imin=max(1,i-Window);
jmin=max(1,j-Window);
@@ -359,7 +285,7 @@ int main(int argc, char **argv)
int nn = kk*nx*ny+jj*nx+ii;
double dsq = double((ii-i)*(ii-i)+(jj-j)*(jj-j)+(kk-k)*(kk-k));
if (id[nn] == 2 && dsq <= Rcrit*Rcrit){
LocalNumber++;
LocalNumber+=1.0;
id[nn]=1;
}
}
@@ -371,6 +297,7 @@ int main(int argc, char **argv)
}
}
// Pack and send the updated ID values
PackID(Dm.sendList_x, Dm.sendCount_x ,sendID_x, id);
PackID(Dm.sendList_X, Dm.sendCount_X ,sendID_X, id);
@@ -427,7 +354,7 @@ int main(int argc, char **argv)
recvID_yZ,Dm.recvCount_yZ,MPI_CHAR,Dm.rank_yZ,recvtag,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,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);
@@ -448,23 +375,50 @@ int main(int argc, char **argv)
UnpackID(Dm.recvList_YZ, Dm.recvCount_YZ ,recvID_YZ, id);
//......................................................................................
MPI_Allreduce(&LocalNumber,&GlobalNumber,1,MPI_INT,MPI_SUM,comm);
MPI_Allreduce(&LocalNumber,&GlobalNumber,1,MPI_DOUBLE,MPI_SUM,comm);
// Layer the inlet with NWP
if (kproc == 0){
for(j=0; j<Ny; j++){
for(i=0; i<Nx; i++){
n = j*nx+i;
// n = nx*ny + j*nx+i;
id[n]=1;
}
}
}
count = 0;
// Layer the outlet with WP
if (kproc == nprocz-1){
for(j=0; j<Ny; j++){
for(i=0; i<Nx; i++){
n = (nz-1)*nx*ny+j*nx+i;
// n = nx*ny + j*nx+i;
id[n]=2;
}
}
}
}
count = 1.f;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
n=k*Nx*Ny+j*Nx+i;
if (id[n] == 1){
count++;
if (id[n] == 2){
count+=1.0;
}
}
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
sat = float(countGlobal)/totalGlobal;
if (rank==0) printf("Final saturation=%f\n",sat);
MPI_Allreduce(&count,&countGlobal,1,MPI_DOUBLE,MPI_SUM,comm);
double sw = countGlobal/totalGlobal;
if (rank==0){
printf("Final saturation=%f\n",sw);
}
sprintf(LocalRankFilename,"ID.%05i",rank);
FILE *ID = fopen(LocalRankFilename,"wb");