refactor morphopen

tests/lbpm_morphopen_pp.cpp

@@ -53,6 +53,7 @@ int main(int argc, char **argv)
 	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; 
@@ -62,9 +63,6 @@ int main(int argc, char **argv)
 			printf("Target saturation %f \n",SW);
 		}
 
-
-	//	}
-
 		//.......................................................................
 		// Reading the domain information file
 		//.......................................................................
@@ -72,66 +70,45 @@ int main(int argc, char **argv)
 		double Lx, Ly, Lz;
 		int i,j,k,n;
 		int BC=0;
+		//  char fluidValue,solidValue;
+		int MAXTIME=1000;
+		int READ_FROM_BLOCK=0;
 
-	if (rank==0){
+		char LocalRankString[8];
-		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;
+		string filename;
-	Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
+		double Rcrit_new, SW;
+		if (argc > 2){
+			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<Database>( filename );
+		auto domain_db = db->getDatabase( "Domain" );
 
-	nx+=2; ny+=2; nz+=2;
+		// Read domain parameters
-	int N = nx*ny*nz;
+		auto L = domain_db->getVector<double>( "L" );
-	char *id;
+		auto size = domain_db->getVector<int>( "n" );
-	id = new char[N];
+		auto nproc = domain_db->getVector<int>( "nproc" );
+		auto ReadValues = domain_db->getVector<char>( "ReadValues" );
+		auto WriteValues = domain_db->getVector<char>( "WriteValues" );
 
-	// Define communication sub-domain -- everywhere
+		nx = size[0];
-	for (int k=0; k<nz; k++){
+		ny = size[1];
-		for (int j=0; j<ny; j++){
+		nz = size[2];
-			for (int i=0; i<nx; i++){
+		nprocx = nproc[0];
-				n = k*nx*ny+j*nx+i;
+		nprocy = nproc[1];
-				Dm.id[n] = 1;
+		nprocz = nproc[2];
-			}
+
-		}
+		int N = (nx+2)*(ny+2)*(nz+2);
-	}
+
-	Dm.CommInit();
+		std::shared_ptr<Domain> Dm (new Domain(domain_db,comm));
+ //		std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
+		for (n=0; n<N; n++) Dm->id[n]=1;
+		Dm->CommInit();
 
 		DoubleArray SignDist(nx,ny,nz);
 		// Read the signed distance from file
@@ -169,106 +146,6 @@ int main(int argc, char **argv)
 		if (rank==0) printf("Media Porosity: %f \n",porosity);
 		if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal);\
 
-	/*	// 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 (size_t 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 (size_t 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");
-		int PoreCount=0;
-		int Count;
-		for (int idx=0; idx<NumBins; idx++){
-		  double BinCenter=MinPoreSize+idx*BinWidth;
-		  Count=GlobalHistogram[idx];
-		  PoreCount+=Count;
-			fprintf(PORESIZE,"%i %f\n",Count,BinCenter);
-		}
-		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<NumBins; idx++){
-		  double BinCenter=MinPoreSize+idx*BinWidth;
-		  Count+=GlobalHistogram[idx];
-		  if (Count<Qval) Q1+=BinWidth;
-		  if (Count<2*Qval) Q2+=BinWidth;
-		  if (Count<3*Qval) Q3+=BinWidth;
-		}
-
-		      printf("Quartiles 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();
 		int iproc = Dm.iproc();
@@ -517,6 +394,7 @@ int main(int argc, char **argv)
 		FILE *ID = fopen(LocalRankFilename,"wb");
 		fwrite(id,1,N,ID);
 		fclose(ID);
+	}
 
 	MPI_Barrier(comm);
 	MPI_Finalize();