From 7dc8f33d444bfc8af9a35a9f3dcc6a719be7f064 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Fri, 14 Apr 2017 14:06:13 -0400
Subject: [PATCH] Added support for multi-input domains in
 lbpm_segmented_decomp

---
 tests/lbpm_segmented_decomp.cpp | 540 ++++++++++++++++----------------
 1 file changed, 277 insertions(+), 263 deletions(-)

diff --git a/tests/lbpm_segmented_decomp.cpp b/tests/lbpm_segmented_decomp.cpp
index 99d6fe55..be0c5d65 100644
--- a/tests/lbpm_segmented_decomp.cpp
+++ b/tests/lbpm_segmented_decomp.cpp
@@ -18,294 +18,308 @@ int main(int argc, char **argv)
 	// Initialize MPI
 	int rank, nprocs;
 	MPI_Init(&argc,&argv);
-	
+
 	MPI_Comm comm = MPI_COMM_WORLD;
 	MPI_Comm_rank(comm,&rank);
 	MPI_Comm_size(comm,&nprocs);
 	{
-	
-	int SOLID=atoi(argv[1]);
-	int NWP=atoi(argv[2]);
-	//char NWP,SOLID;
-	//SOLID=argv[1][0];
-	//NWP=argv[2][0];
-	if (rank==0){
-	  printf("Solid Label: %i \n",SOLID);
-	  printf("NWP Label: %i \n",NWP);
-	}
-	
-    //.......................................................................
-    // Reading the domain information file
-    //.......................................................................
-    int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
-    double Lx, Ly, Lz;
-    int Nx,Ny,Nz;
-    int i,j,k,n;
-	int BC=0;
-    char Filename[40];
-    int xStart,yStart,zStart;
-  //  char fluidValue,solidValue;
 
-    std::vector<char> solidValues;
-    std::vector<char> nwpValues;
-    std::string line;
 
-    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;
+		bool MULTINPUT=false;
 
-    	ifstream image("Segmented.in");
-    	image >> Filename; 	// Name of data file containing segmented data
-    	image >> Nx;   		// size of the binary file
-    	image >> Ny;
-    	image >> Nz;
-    	image >> xStart;	// offset for the starting voxel
-    	image >> yStart;
-    	image >> zStart;
+		int NWP,SOLID,rank_offset;
+		SOLID=atoi(argv[1]);
+		NWP=atoi(argv[2]);
+		//char NWP,SOLID;
+		//SOLID=argv[1][0];
+		//NWP=argv[2][0];
+		if (rank==0){
+			printf("Solid Label: %i \n",SOLID);
+			printf("NWP Label: %i \n",NWP);
+		}
+		if (argc > 3){
+			rank_offset = atoi(argv[3]);
+		}
+		else{
+			MULTINPUT=true;
+			rank_offset=0;
+		}
 
-    }
-	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_Bcast(&Ny,1,MPI_INT,0,comm);
-	MPI_Bcast(&Ny,1,MPI_INT,0,comm);
-	MPI_Bcast(&Nz,1,MPI_INT,0,comm);
-	MPI_Bcast(&xStart,1,MPI_INT,0,comm);
-	MPI_Bcast(&yStart,1,MPI_INT,0,comm);
-	MPI_Bcast(&zStart,1,MPI_INT,0,comm);
-	//.................................................
-	MPI_Barrier(comm);
+		//.......................................................................
+		// Reading the domain information file
+		//.......................................................................
+		int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
+		double Lx, Ly, Lz;
+		int Nx,Ny,Nz;
+		int i,j,k,n;
+		int BC=0;
+		char Filename[40];
+		int xStart,yStart,zStart;
+		//  char fluidValue,solidValue;
 
-    // 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 *SegData = NULL;
-    // Rank=0 reads the entire segmented data and distributes to worker processes
-    if (rank==0){
-    	printf("Dimensions of segmented image: %i x %i x %i \n",Nx,Ny,Nz);
-	SegData = new char[Nx*Ny*Nz];
-    	FILE *SEGDAT = fopen(Filename,"rb");
-    	if (SEGDAT==NULL) ERROR("Error reading segmented data");
-    	size_t ReadSeg;
-    	ReadSeg=fread(SegData,1,Nx*Ny*Nz,SEGDAT);
-        if (ReadSeg != size_t(Nx*Ny*Nz)) printf("lbpm_segmented_decomp: Error reading segmented data (rank=%i)\n",rank);
-    	fclose(SEGDAT);
-        printf("Read segmented data from %s \n",Filename);
-    }
-    MPI_Barrier(comm);
+		std::vector<char> solidValues;
+		std::vector<char> nwpValues;
+		std::string line;
 
-    // Get the rank info
-    int N = (nx+2)*(ny+2)*(nz+2);
-	Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
-	for (k=0;k<nz+2;k++){
-		for (j=0;j<ny+2;j++){
-			for (i=0;i<nx+2;i++){
-				n = k*(nx+2)*(ny+2)+j*(nx+2)+i;
-				Dm.id[n] = 1;
+		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;
+
+			ifstream image("Segmented.in");
+			image >> Filename; 	// Name of data file containing segmented data
+			image >> Nx;   		// size of the binary file
+			image >> Ny;
+			image >> Nz;
+			image >> xStart;	// offset for the starting voxel
+			image >> yStart;
+			image >> zStart;
+
+		}
+		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_Bcast(&Ny,1,MPI_INT,0,comm);
+		MPI_Bcast(&Ny,1,MPI_INT,0,comm);
+		MPI_Bcast(&Nz,1,MPI_INT,0,comm);
+		MPI_Bcast(&xStart,1,MPI_INT,0,comm);
+		MPI_Bcast(&yStart,1,MPI_INT,0,comm);
+		MPI_Bcast(&zStart,1,MPI_INT,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 *SegData = NULL;
+		// Rank=0 reads the entire segmented data and distributes to worker processes
+		if (rank==0){
+			printf("Dimensions of segmented image: %i x %i x %i \n",Nx,Ny,Nz);
+			SegData = new char[Nx*Ny*Nz];
+			FILE *SEGDAT = fopen(Filename,"rb");
+			if (SEGDAT==NULL) ERROR("Error reading segmented data");
+			size_t ReadSeg;
+			ReadSeg=fread(SegData,1,Nx*Ny*Nz,SEGDAT);
+			if (ReadSeg != size_t(Nx*Ny*Nz)) printf("lbpm_segmented_decomp: Error reading segmented data (rank=%i)\n",rank);
+			fclose(SEGDAT);
+			printf("Read segmented data from %s \n",Filename);
+		}
+		MPI_Barrier(comm);
+
+		// Get the rank info
+		int N = (nx+2)*(ny+2)*(nz+2);
+		Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
+		for (k=0;k<nz+2;k++){
+			for (j=0;j<ny+2;j++){
+				for (i=0;i<nx+2;i++){
+					n = k*(nx+2)*(ny+2)+j*(nx+2)+i;
+					Dm.id[n] = 1;
+				}
 			}
 		}
-	}
-	Dm.CommInit(comm);
+		Dm.CommInit(comm);
 
-	// number of sites to use for periodic boundary condition transition zone
-	int z_transition_size = (nprocz*nz - (Nz - zStart))/2;
-	if (z_transition_size < 0) z_transition_size=0;
+		// number of sites to use for periodic boundary condition transition zone
+		int z_transition_size = (nprocz*nz - (Nz - zStart))/2;
+		if (z_transition_size < 0) z_transition_size=0;
 
-	// Set up the sub-domains
-	if (rank==0){
-		printf("Distributing subdomains across %i processors \n",nprocs);
-		printf("Process grid: %i x %i x %i \n",Dm.nprocx,Dm.nprocy,Dm.nprocz);
-		printf("Subdomain size: %i \n",N);
-		printf("Size of transition region: %i \n", z_transition_size);
-		char *tmp;
-		tmp = new char[N];
-		for (int kp=0; kp<nprocz; kp++){
-			for (int jp=0; jp<nprocy; jp++){
-				for (int ip=0; ip<nprocx; ip++){
-					// rank of the process that gets this subdomain
-					int rnk = kp*Dm.nprocx*Dm.nprocy + jp*Dm.nprocx + ip;
-					// Pack and send the subdomain for rnk
-					for (k=0;k<nz+2;k++){
-						for (j=0;j<ny+2;j++){
-							for (i=0;i<nx+2;i++){
-								int x = xStart + ip*nx + i-1;
-								int y = yStart + jp*ny + j-1;
-						//		int z = zStart + kp*nz + k-1;
-								int z = zStart + kp*nz + k-1 - z_transition_size;
-								if (z<zStart) 	z=zStart;
-								if (!(z<Nz))	z=Nz-1;
-								int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
-								int nglobal = z*Nx*Ny+y*Nx+x;
-								tmp[nlocal] = SegData[nglobal];
-							}
-						}
-					}
-					if (rnk==0){
+		// Set up the sub-domains
+		if (rank==0){
+			printf("Distributing subdomains across %i processors \n",nprocs);
+			printf("Process grid: %i x %i x %i \n",Dm.nprocx,Dm.nprocy,Dm.nprocz);
+			printf("Subdomain size: %i \n",N);
+			printf("Size of transition region: %i \n", z_transition_size);
+			char *tmp;
+			tmp = new char[N];
+			for (int kp=0; kp<nprocz; kp++){
+				for (int jp=0; jp<nprocy; jp++){
+					for (int ip=0; ip<nprocx; ip++){
+						// rank of the process that gets this subdomain
+						int rnk = kp*Dm.nprocx*Dm.nprocy + jp*Dm.nprocx + ip;
+						// Pack and send the subdomain for rnk
 						for (k=0;k<nz+2;k++){
 							for (j=0;j<ny+2;j++){
 								for (i=0;i<nx+2;i++){
+									int x = xStart + ip*nx + i-1;
+									int y = yStart + jp*ny + j-1;
+									//		int z = zStart + kp*nz + k-1;
+									int z = zStart + kp*nz + k-1 - z_transition_size;
+									if (z<zStart) 	z=zStart;
+									if (!(z<Nz))	z=Nz-1;
 									int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
-									Dm.id[nlocal] = tmp[nlocal];
+									int nglobal = z*Nx*Ny+y*Nx+x;
+									tmp[nlocal] = SegData[nglobal];
 								}
 							}
 						}
-					}
-					else{
-						printf("Sending data to process %i \n", rnk);
-						MPI_Send(tmp,N,MPI_CHAR,rnk,15,comm);
+						if (rnk==0){
+							for (k=0;k<nz+2;k++){
+								for (j=0;j<ny+2;j++){
+									for (i=0;i<nx+2;i++){
+										int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
+										Dm.id[nlocal] = tmp[nlocal];
+									}
+								}
+							}
+						}
+						else{
+							printf("Sending data to process %i \n", rnk);
+							MPI_Send(tmp,N,MPI_CHAR,rnk,15,comm);
+						}
 					}
 				}
 			}
 		}
-	}
-	else{
-		// Recieve the subdomain from rank = 0
-		printf("Ready to recieve data %i at process %i \n", N,rank);
-		MPI_Recv(Dm.id,N,MPI_CHAR,0,15,comm,MPI_STATUS_IGNORE);
+		else{
+			// Recieve the subdomain from rank = 0
+			printf("Ready to recieve data %i at process %i \n", N,rank);
+			MPI_Recv(Dm.id,N,MPI_CHAR,0,15,comm,MPI_STATUS_IGNORE);
+		}
+		MPI_Barrier(comm);
+
+		nx+=2; ny+=2; nz+=2;
+		N=nx*ny*nz;
+
+		if (rank==0) printf("All sub-domains recieved \n");
+		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 (Dm.id[n]==char(SOLID))     Dm.id[n] = 0;
+					else if (Dm.id[n]==char(NWP))  Dm.id[n] = 1;
+					else                     Dm.id[n] = 2;
+
+				}
+			}
+		}
+		if (rank==0) printf("Domain set \n");
+
+		int count = 0;
+		int total = 0;
+		int countGlobal = 0;
+		int totalGlobal = 0;
+		for (k=1;k<nz-1;k++){
+			for (j=1;j<ny-1;j++){
+				for (i=1;i<nx-1;i++){
+					n=k*nx*ny+j*nx+i;
+					total++;
+					if (Dm.id[n] == 0){
+						count++;
+					}
+				}
+			}
+		}
+		MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
+		MPI_Allreduce(&total,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
+
+
+		float porosity = float(totalGlobal-countGlobal)/totalGlobal;
+		if (rank==0) printf("Porosity=%f\n",porosity);
+
+		if (rank==0){
+			int xstart = xStart;    // Is this correct?
+			int ystart = yStart;
+			int zstart = zStart;
+			//totalGlobal=(Nx-xstart)*(Ny-ystart)*(Nz-zstart);
+			countGlobal = 0;
+			for (k=zstart; k<zstart+nprocz*(nz-2); k++){
+				for (j=ystart; j<ystart+nprocy*(ny-2); j++){
+					for (i=xstart; i<xstart+nprocx*(nx-2); i++){
+
+						n=k*Nx*Ny+j*Nx+i;
+						if (n < Nx*Ny*Nz){
+							if (SegData[n] == char(SOLID)){
+								countGlobal++;
+							}
+						}
+					}
+				}
+			}
+			float porosity = float(totalGlobal-countGlobal)/totalGlobal;
+			printf("Original Porosity=%f\n",porosity);
+		}
+
+		count = 0;
+		total = 0;
+		countGlobal = 0;
+		totalGlobal = 0;
+		for (k=1;k<nz-1;k++){
+			for (j=1;j<ny-1;j++){
+				for (i=1;i<nx-1;i++){
+					n=k*nx*ny+j*nx+i;
+					if (Dm.id[n] != 0)      total++;
+					if (Dm.id[n] == 2)	count++;
+				}
+			}
+		}
+		MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
+		MPI_Allreduce(&total,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
+		float saturation = float(countGlobal)/totalGlobal;
+		if (rank==0) printf("wetting phase saturation=%f\n",saturation);
+
+
+		char LocalRankFilename[40];
+
+		sprintf(LocalRankFilename,"ID.%05i",rank+rank_offset);
+		FILE *ID = fopen(LocalRankFilename,"wb");
+		fwrite(Dm.id,1,N,ID);
+		//    fwrite(Distance.get(),8,Distance.length(),ID);
+		fclose(ID);
+
+		if (!MULTINPUT){
+
+			if (rank==0) printf("Writing symmetric domain reflection\n");
+			MPI_Barrier(comm);
+			int symrank,sympz;
+			sympz = 2*nprocz - Dm.kproc -1;
+			symrank = sympz*nprocx*nprocy + Dm.jproc*nprocx + Dm.iproc;
+
+			//    DoubleArray SymDist(nx,ny,nz);
+			char *symid;
+			symid = new char [N];
+			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;
+						int nsym=(nz-k-1)*nx*ny+j*nx+i;
+						symid[nsym] = Dm.id[n];
+						//SymDist(i,j,nz-k-1)=Distance(i,j,k);
+					}
+				}
+			}
+
+			sprintf(LocalRankFilename,"ID.%05i",symrank);
+			FILE *SYMID = fopen(LocalRankFilename,"wb");
+			//    fwrite(SymDist.get(),8,SymDist.length(),SYMDIST);
+			fwrite(symid,1,N,SYMID);
+			fclose(SYMID);
+		}
 	}
 	MPI_Barrier(comm);
-
-	nx+=2; ny+=2; nz+=2;
-	N=nx*ny*nz;
-
-	if (rank==0) printf("All sub-domains recieved \n");
-    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 (Dm.id[n]==char(SOLID))     Dm.id[n] = 0;
-				else if (Dm.id[n]==char(NWP))  Dm.id[n] = 1;
-				else                     Dm.id[n] = 2;
-				
-			}
-		}
-	}
-	if (rank==0) printf("Domain set \n");
-
-    int count = 0;
-    int total = 0;
-    int countGlobal = 0;
-    int totalGlobal = 0;
-    for (k=1;k<nz-1;k++){
-    	for (j=1;j<ny-1;j++){
-    		for (i=1;i<nx-1;i++){
-    			n=k*nx*ny+j*nx+i;
-    			total++;
-    			if (Dm.id[n] == 0){
-    				count++;
-    			}
-    		}
-    	}
-    }
-    MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
-    MPI_Allreduce(&total,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
-
-
-    float porosity = float(totalGlobal-countGlobal)/totalGlobal;
-    if (rank==0) printf("Porosity=%f\n",porosity);
-
-    if (rank==0){
-        int xstart = xStart;    // Is this correct?
-        int ystart = yStart;
-        int zstart = zStart;
-    	//totalGlobal=(Nx-xstart)*(Ny-ystart)*(Nz-zstart);
-    	countGlobal = 0;
-    	for (k=zstart; k<zstart+nprocz*(nz-2); k++){
-    		for (j=ystart; j<ystart+nprocy*(ny-2); j++){
-    			for (i=xstart; i<xstart+nprocx*(nx-2); i++){
-
-    				n=k*Nx*Ny+j*Nx+i;
-    				if (n < Nx*Ny*Nz){
-    					if (SegData[n] == char(SOLID)){
-    						countGlobal++;
-    					}
-    				}
-    			}
-    		}
-    	}
-    	float porosity = float(totalGlobal-countGlobal)/totalGlobal;
-    	printf("Original Porosity=%f\n",porosity);
-    }
-
-    count = 0;
-    total = 0;
-    countGlobal = 0;
-    totalGlobal = 0;
-    for (k=1;k<nz-1;k++){
-    	for (j=1;j<ny-1;j++){
-    		for (i=1;i<nx-1;i++){
-    			n=k*nx*ny+j*nx+i;
-    			if (Dm.id[n] != 0)      total++;
-    			if (Dm.id[n] == 2)	count++;
-    		}
-    	}
-    }
-    MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
-    MPI_Allreduce(&total,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
-    float saturation = float(countGlobal)/totalGlobal;
-    if (rank==0) printf("wetting phase saturation=%f\n",saturation);
-
-
-   char LocalRankFilename[40];
-
-    sprintf(LocalRankFilename,"ID.%05i",rank);
-    FILE *ID = fopen(LocalRankFilename,"wb");
-    fwrite(Dm.id,1,N,ID);
-//    fwrite(Distance.get(),8,Distance.length(),ID);
-    fclose(ID);
-
-    if (rank==0) printf("Writing symmetric domain reflection\n");
-    MPI_Barrier(comm);
-    int symrank,sympz;
-    sympz = 2*nprocz - Dm.kproc -1;
-    symrank = sympz*nprocx*nprocy + Dm.jproc*nprocx + Dm.iproc;
-
-//    DoubleArray SymDist(nx,ny,nz);
-    char *symid;
-    symid = new char [N];
-	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;
-				int nsym=(nz-k-1)*nx*ny+j*nx+i;
-				symid[nsym] = Dm.id[n];
-				//SymDist(i,j,nz-k-1)=Distance(i,j,k);
-			}
-		}
-    }
-
-    sprintf(LocalRankFilename,"ID.%05i",symrank);
-    FILE *SYMID = fopen(LocalRankFilename,"wb");
-//    fwrite(SymDist.get(),8,SymDist.length(),SYMDIST);
-    fwrite(symid,1,N,SYMID);
-    fclose(SYMID);
-    }
-    MPI_Barrier(comm);
-    MPI_Finalize();
+	MPI_Finalize();
 }