From 163f9a2d5d1aa364f2167e73f98c6b6dfb91d05f Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 21 Feb 2018 11:31:33 -0500
Subject: [PATCH] refactor lbpm_permeability_simulator

---
 tests/CMakeLists.txt                  |   2 +-
 tests/lbpm_permeability_simulator.cpp | 895 +++++++++++---------------
 2 files changed, 366 insertions(+), 531 deletions(-)

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index c9940957..9202f034 100755
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,8 +1,8 @@
 # Copy files for the tests
 #ADD_LBPM_EXECUTABLE( lbpm_nonnewtonian_simulator ) 
-#ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
 #ADD_LBPM_EXECUTABLE( lbpm_nondarcy_simulator )
 #ADD_LBPM_EXECUTABLE( lbpm_color_simulator )
+ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_sphere_pp )
 ADD_LBPM_EXECUTABLE( lbpm_random_pp )
diff --git a/tests/lbpm_permeability_simulator.cpp b/tests/lbpm_permeability_simulator.cpp
index cdb100d3..02201d44 100644
--- a/tests/lbpm_permeability_simulator.cpp
+++ b/tests/lbpm_permeability_simulator.cpp
@@ -20,78 +20,6 @@
 
 using namespace std;
 
-//*************************************************************************
-// Implementation of Steady State Single-Phase LBM for permeability measurement
-//*************************************************************************
-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<count; idx++){
-		n = list[idx];
-		sendbuf[idx] = ID[n];
-	}
-}
-//***************************************************************************************
-
-inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
-	// Fill in the phase ID values from neighboring processors
-	// This unpacks the values once they have been recieved from neighbors
-	int idx,n;
-
-	for (idx=0; idx<count; idx++){
-		n = list[idx];
-		ID[n] = recvbuf[idx];
-	}
-}
-
-//***************************************************************************************
-
-inline void ZeroHalo(double *Data, int Nx, int Ny, int Nz)
-{
-	int i,j,k,n;
-	for (k=0;k<Nz;k++){
-		for (j=0;j<Ny;j++){
-			i=0;
-			n = k*Nx*Ny+j*Nx+i;
-			Data[2*n] = 0.0;
-			Data[2*n+1] = 0.0;
-			i=Nx-1;
-			n = k*Nx*Ny+j*Nx+i;
-			Data[2*n] = 0.0;
-			Data[2*n+1] = 0.0;
-		}
-	}
-
-	for (k=0;k<Nz;k++){
-		for (i=0;i<Nx;i++){
-			j=0;
-			n = k*Nx*Ny+j*Nx+i;
-			Data[2*n] = 0.0;
-			Data[2*n+1] = 0.0;
-			j=Ny-1;
-			n = k*Nx*Ny+j*Nx+i;
-			Data[2*n] = 0.0;
-			Data[2*n+1] = 0.0;
-		}
-	}
-
-	for (j=0;j<Ny;j++){
-		for (i=0;i<Nx;i++){
-			k=0;
-			n = k*Nx*Ny+j*Nx+i;
-			Data[2*n] = 0.0;
-			Data[2*n+1] = 0.0;
-			k=Nz-1;
-			n = k*Nx*Ny+j*Nx+i;
-			Data[2*n] = 0.0;
-			Data[2*n+1] = 0.0;
-		}
-	}
-}
-//***************************************************************************************
-
 
 int main(int argc, char **argv)
 {
@@ -101,476 +29,383 @@ 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);
 	{
-	// parallel domain size (# of sub-domains)
-	int nprocx,nprocy,nprocz;
-	int iproc,jproc,kproc;
-	//*****************************************
-	// MPI ranks for all 18 neighbors
-	//**********************************
-	int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z;
-	int rank_xy,rank_XY,rank_xY,rank_Xy;
-	int rank_xz,rank_XZ,rank_xZ,rank_Xz;
-	int rank_yz,rank_YZ,rank_yZ,rank_Yz;
-	//**********************************
-	MPI_Request req1[18],req2[18];
-	MPI_Status stat1[18],stat2[18];
+		// parallel domain size (# of sub-domains)
+		int nprocx,nprocy,nprocz;
+		int iproc,jproc,kproc;
 
-	if (rank == 0){
-		printf("********************************************************\n");
-		printf("Running Single Phase Permeability Calculation \n");
-		printf("********************************************************\n");
-	}
-	
-	// Variables that specify the computational domain  
-	string FILENAME;
-	int Nx,Ny,Nz;		// local sub-domain size
-	int nspheres;		// number of spheres in the packing
-	double Lx,Ly,Lz;	// Domain length
-	double D = 1.0;		// reference length for non-dimensionalization
-	// Color Model parameters
-	int timestepMax, interval;
-	double tau,Fx,Fy,Fz,tol,err;
-	double din,dout;
-	bool pBC,Restart;
-	int i,j,k,n;
-
-	int RESTART_INTERVAL=20000;
-
-	if (rank==0){
-		//.............................................................
-		//		READ SIMULATION PARMAETERS FROM INPUT FILE
-		//.............................................................
-		ifstream input("Permeability.in");
-		// Line 1: model parameters (tau, alpha, beta, das, dbs)
-		input >> tau;			// Viscosity parameter
-		// Line 2: External force components (Fx,Fy, Fz)
-		input >> Fx;
-		input >> Fy;
-		input >> Fz;
-		// Line 3: Pressure Boundary conditions
-		input >> Restart;
-		input >> pBC;
-		input >> din;
-		input >> dout;
-		// Line 4: time-stepping criteria
-		input >> timestepMax;		// max no. of timesteps
-		input >> interval;			// restart interval
-		input >> tol;				// error tolerance
-		//.............................................................
-
-		//.......................................................................
-		// Reading the domain information file
-		//.......................................................................
-		ifstream domain("Domain.in");
-		domain >> nprocx;
-		domain >> nprocy;
-		domain >> nprocz;
-		domain >> Nx;
-		domain >> Ny;
-		domain >> Nz;
-		domain >> nspheres;
-		domain >> Lx;
-		domain >> Ly;
-		domain >> Lz;
-		//.......................................................................
-		
-	}
-	// **************************************************************
-	// Broadcast simulation parameters from rank 0 to all other procs
-	MPI_Barrier(comm);
-	//.................................................
-	MPI_Bcast(&tau,1,MPI_DOUBLE,0,comm);
-	//MPI_Bcast(&pBC,1,MPI_LOGICAL,0,comm);
-	//	MPI_Bcast(&Restart,1,MPI_LOGICAL,0,comm);
-	MPI_Bcast(&din,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&dout,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&Fx,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&Fy,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&Fz,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&timestepMax,1,MPI_INT,0,comm);
-	MPI_Bcast(&interval,1,MPI_INT,0,comm);
-	MPI_Bcast(&tol,1,MPI_DOUBLE,0,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);
-	
-	RESTART_INTERVAL=interval;
-	// **************************************************************
-	// **************************************************************
-	double rlxA = 1.f/tau;
-	double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA);
-
-
-	if (nprocs != nprocx*nprocy*nprocz){
-		printf("nprocx =  %i \n",nprocx);
-		printf("nprocy =  %i \n",nprocy);
-		printf("nprocz =  %i \n",nprocz);
-		INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!");
-	}
-
-	if (rank==0){
-		printf("********************************************************\n");
-		printf("tau = %f \n", tau);
-		printf("Force(x) = %.5g \n", Fx);
-		printf("Force(y) = %.5g \n", Fy);
-		printf("Force(z) = %.5g \n", Fz);
-		printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
-		printf("Process grid = %i x %i x %i\n",nprocx,nprocy,nprocz);
-		printf("********************************************************\n");
-	}
-
-	double viscosity=(tau-0.5)/3.0;
-	// Initialized domain and averaging framework for Two-Phase Flow
-	int BC=pBC;
-	Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
-	TwoPhase Averages(Dm);
-
-	InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc,
-			 	 	 rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z,
-			 	 	 rank_xy, rank_XY, rank_xY, rank_Xy, rank_xz, rank_XZ, rank_xZ, rank_Xz,
-			 	 	 rank_yz, rank_YZ, rank_yZ, rank_Yz );
-	 
-	MPI_Barrier(comm);
-
-	Nx += 2;	Ny += 2;	Nz += 2;
-
-	int N = Nx*Ny*Nz;
-	int dist_mem_size = N*sizeof(double);
-
-	//.......................................................................
-	if (rank == 0)	printf("Read input media... \n");
-	//.......................................................................
-	
-	//.......................................................................
-	// Filenames used
-	char LocalRankString[8];
-	char LocalRankFilename[40];
-	char LocalRestartFile[40];
-	char tmpstr[10];
-	sprintf(LocalRankString,"%05d",rank);
-	sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
-	sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
-	
-//	printf("Local File Name =  %s \n",LocalRankFilename);
-	// .......... READ THE INPUT FILE .......................................
-//	char value;
-	char *id;
-	id = new char[N];
-	int sum = 0;
-	double sum_local;
-	double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
-	if (pBC) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
-	double porosity, pore_vol;
-	//...........................................................................
-	if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
-
-	//.......................................................................
-	sprintf(LocalRankString,"%05d",rank);
-//	sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
-//	WriteLocalSolidID(LocalRankFilename, id, N);
-	sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
-	ReadBinaryFile(LocalRankFilename, Averages.SDs.data(), N);
-	MPI_Barrier(comm);
-	if (rank == 0) cout << "Domain set." << endl;
-	
-	//.......................................................................
-	// Assign the phase ID field based on the signed distance
-	//.......................................................................
-	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;
-				id[n] = 0;
-			}
+		if (rank == 0){
+			printf("********************************************************\n");
+			printf("Running Single Phase Permeability Calculation \n");
+			printf("********************************************************\n");
 		}
-	}
-	sum=0;
-	pore_vol = 0.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 (Averages.SDs(n) > 0.0){
-					id[n] = 2;	
-				}
-				// compute the porosity (actual interface location used)
-				if (Averages.SDs(n) > 0.0){
-					sum++;	
-				}
-			}
+
+		// Variables that specify the computational domain  
+		string FILENAME;
+		int Nx,Ny,Nz;		// local sub-domain size
+		int nspheres;		// number of spheres in the packing
+		double Lx,Ly,Lz;	// Domain length
+		double D = 1.0;		// reference length for non-dimensionalization
+		// Color Model parameters
+		int timestepMax, interval;
+		double tau,Fx,Fy,Fz,tol,err;
+		double din,dout;
+		bool pBC,Restart;
+		int i,j,k,n;
+
+		int RESTART_INTERVAL=20000;
+
+		if (rank==0){
+			//.............................................................
+			//		READ SIMULATION PARMAETERS FROM INPUT FILE
+			//.............................................................
+			ifstream input("Permeability.in");
+			// Line 1: model parameters (tau, alpha, beta, das, dbs)
+			input >> tau;			// Viscosity parameter
+			// Line 2: External force components (Fx,Fy, Fz)
+			input >> Fx;
+			input >> Fy;
+			input >> Fz;
+			// Line 3: Pressure Boundary conditions
+			input >> Restart;
+			input >> pBC;
+			input >> din;
+			input >> dout;
+			// Line 4: time-stepping criteria
+			input >> timestepMax;		// max no. of timesteps
+			input >> interval;			// restart interval
+			input >> tol;				// error tolerance
+			//.............................................................
+
+			//.......................................................................
+			// Reading the domain information file
+			//.......................................................................
+			ifstream domain("Domain.in");
+			domain >> nprocx;
+			domain >> nprocy;
+			domain >> nprocz;
+			domain >> Nx;
+			domain >> Ny;
+			domain >> Nz;
+			domain >> nspheres;
+			domain >> Lx;
+			domain >> Ly;
+			domain >> Lz;
+			//.......................................................................
+
 		}
-	}
-	
-	// Set up kstart, kfinish so that the reservoirs are excluded from averaging
-	int kstart,kfinish;
-	kstart = 1;
-	kfinish = Nz-1;
-	if (pBC && kproc==0)		kstart = 4;
-	if (pBC && kproc==nprocz-1)	kfinish = Nz-4;
-
-	// Compute the pore volume
-	sum_local = 0.0;
-	for ( k=kstart;k<kfinish;k++){
-		for ( j=1;j<Ny-1;j++){
-			for ( i=1;i<Nx-1;i++){
-				n = k*Nx*Ny+j*Nx+i;
-				if (id[n] > 0){
-					sum_local += 1.0;
-				}
-			}
-		}
-	}
-	MPI_Allreduce(&sum_local,&pore_vol,1,MPI_DOUBLE,MPI_SUM,comm);
-//	MPI_Allreduce(&sum_local,&porosity,1,MPI_DOUBLE,MPI_SUM,comm);
-	porosity = pore_vol*iVol_global;
-	if (rank==0) printf("Media porosity = %f \n",porosity);
-	//.........................................................
-	// If pressure boundary conditions are applied remove solid
-	if (pBC && kproc == 0){
-		for (k=0; k<3; k++){
-			for (j=0;j<Ny;j++){
-				for (i=0;i<Nx;i++){
-					n = k*Nx*Ny+j*Nx+i;
-					id[n] = 1;
-					Averages.SDs(n) = max(Averages.SDs(n),1.0*(2.5-k));
-				}					
-			}
-		}
-	}
-	if (pBC && kproc == nprocz-1){
-		for (k=Nz-3; k<Nz; k++){
-			for (j=0;j<Ny;j++){
-				for (i=0;i<Nx;i++){
-					n = k*Nx*Ny+j*Nx+i;
-					id[n] = 2;
-					Averages.SDs(n) = max(Averages.SDs(n),1.0*(k-Nz+2.5));
-				}					
-			}
-		}
-	}
-	//.........................................................
-	// don't perform computations at the eight corners
-	id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0;
-	id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
-	//.........................................................
-	// Initialize communication structures in averaging domain
-	for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = id[i];
-	Dm.CommInit(comm);
-
-	//...........................................................................
-	if (rank==0)	printf ("Create ScaLBL_Communicator \n");
-	// Create a communicator for the device
-	ScaLBL_Communicator ScaLBL_Comm(Dm);
-
-	//...........device phase ID.................................................
-	if (rank==0)	printf ("Copying phase ID to device \n");
-	char *ID;
-	ScaLBL_AllocateDeviceMemory((void **) &ID, N);						// Allocate device memory
-	// Copy to the device
-	ScaLBL_CopyToDevice(ID, id, N);
-	//...........................................................................
-
-	//...........................................................................
-	//				MAIN  VARIABLES ALLOCATED HERE
-	//...........................................................................
-	// LBM variables
-	if (rank==0)	printf ("Allocating distributions \n");
-	//......................device distributions.................................
-	double *f_even,*f_odd;
-	//...........................................................................
-	ScaLBL_AllocateDeviceMemory((void **) &f_even, 10*dist_mem_size);	// Allocate device memory
-	ScaLBL_AllocateDeviceMemory((void **) &f_odd, 9*dist_mem_size);	// Allocate device memory
-	//...........................................................................
-	double *Velocity, *Pressure, *dvcSignDist;
-	//...........................................................................
-	ScaLBL_AllocateDeviceMemory((void **) &Pressure, dist_mem_size);
-	ScaLBL_AllocateDeviceMemory((void **) &dvcSignDist, dist_mem_size);
-	ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
-	//...........................................................................
-
-	// Copy signed distance for device initialization
-	ScaLBL_CopyToDevice(dvcSignDist, Averages.SDs.data(), dist_mem_size);
-	//...........................................................................
-
-	int logcount = 0; // number of surface write-outs
-	
-	//...........................................................................
-	//				MAIN  VARIABLES INITIALIZED HERE
-	//...........................................................................
-	//...........................................................................
-	if (rank==0)	printf("Setting the distributions, size = %i\n", N);
-	//...........................................................................
-	ScaLBL_D3Q19_Init(ID, f_even, f_odd, Nx, Ny, Nz);
-	//......................................................................
-
-	//.......................................................................
-	// Finalize setup for averaging domain
-	//Averages.SetupCubes(Dm);
-	Averages.UpdateSolid();
-	// Initialize two phase flow variables (all wetting phase)
-	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;
-				Averages.Phase(i,j,k) = -1.0;
-				Averages.SDn(i,j,k) = Averages.Phase(i,j,k);
-				Averages.Phase_tplus(i,j,k) = Averages.SDn(i,j,k);
-				Averages.Phase_tminus(i,j,k) = Averages.SDn(i,j,k);
-				Averages.DelPhi(i,j,k) = 0.0;
-				Averages.Press(i,j,k) = 0.0;
-				Averages.Vel_x(i,j,k) = 0.0;
-				Averages.Vel_y(i,j,k) = 0.0;
-				Averages.Vel_z(i,j,k) = 0.0;
-			}
-		}
-	}
-
-	//.......................................................................
-
-	if (rank==0 && pBC){
-		printf("Setting inlet pressure = %f \n", din);
-		printf("Setting outlet pressure = %f \n", dout);
-	}
-	if (pBC && kproc == 0)	{
-		ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-	}
-		
-	if (pBC && kproc == nprocz-1){
-		ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-	}
-
-	int timestep = 0;
-	if (rank==0) printf("********************************************************\n");
-	if (rank==0)	printf("No. of timesteps: %i \n", timestepMax);
-
-	//.......create and start timer............
-	double starttime,stoptime,cputime;
-	MPI_Barrier(comm);
-	starttime = MPI_Wtime();
-	//.........................................
-	
-	double D32,Fo,Re,velocity,err1D,mag_force,vel_prev;
-	err = vel_prev = 1.0;
-	if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol);
-	//************ MAIN ITERATION LOOP ***************************************/
-	while (timestep < timestepMax && err > tol ){
-
-		//*************************************************************************
-		// Fused Color Gradient and Collision 
-		//*************************************************************************
-		ScaLBL_D3Q19_MRT( ID,f_even,f_odd,rlxA,rlxB,Fx,Fy,Fz,Nx,Ny,Nz);
-		//*************************************************************************
-
-		//*************************************************************************
-		// Pack and send the D3Q19 distributions
-		ScaLBL_Comm.SendD3Q19(f_even, f_odd);
-		//*************************************************************************
-		// 		Swap the distributions for momentum transport
-		//*************************************************************************
-		ScaLBL_D3Q19_Swap(ID, f_even, f_odd, Nx, Ny, Nz);
-		//*************************************************************************
-		// Wait for communications to complete and unpack the distributions
-		ScaLBL_Comm.RecvD3Q19(f_even, f_odd);
-		//*************************************************************************
-		
-		if (pBC && kproc == 0)	{
-			ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-		}
-			
-		if (pBC && kproc == nprocz-1){
-			ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-		}
-		//...................................................................................
-		ScaLBL_DeviceBarrier();
+		// **************************************************************
+		// Broadcast simulation parameters from rank 0 to all other procs
+		MPI_Barrier(comm);
+		//.................................................
+		MPI_Bcast(&tau,1,MPI_DOUBLE,0,comm);
+		//MPI_Bcast(&pBC,1,MPI_LOGICAL,0,comm);
+		//	MPI_Bcast(&Restart,1,MPI_LOGICAL,0,comm);
+		MPI_Bcast(&din,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&dout,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&Fx,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&Fy,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&Fz,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&timestepMax,1,MPI_INT,0,comm);
+		MPI_Bcast(&interval,1,MPI_INT,0,comm);
+		MPI_Bcast(&tol,1,MPI_DOUBLE,0,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);
 
-		// Timestep completed!
-		timestep++;
+		RESTART_INTERVAL=interval;
+		// **************************************************************
+		// **************************************************************
+		double rlxA = 1.f/tau;
+		double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA);
 
-		if (timestep%500 == 0){
-			//...........................................................................
-			// Copy the data for for the analysis timestep
-			//...........................................................................
-			// Copy the phase from the GPU -> CPU
-			//...........................................................................
-			ScaLBL_DeviceBarrier();
-			ScaLBL_D3Q19_Pressure(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
-			ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-			ScaLBL_CopyToHost(Averages.Press.data(),Pressure,N*sizeof(double));
-			ScaLBL_CopyToHost(Averages.Vel_x.data(),&Velocity[0],N*sizeof(double));
-			ScaLBL_CopyToHost(Averages.Vel_y.data(),&Velocity[N],N*sizeof(double));
-			ScaLBL_CopyToHost(Averages.Vel_z.data(),&Velocity[2*N],N*sizeof(double));
 
-			// Way more work than necessary -- this is just to get the solid interfacial area!!
-			Averages.Initialize();
-			Averages.UpdateMeshValues();
-			Averages.ComputeLocal();
-			Averages.Reduce();
-			
-			double vawx = -Averages.vaw_global(0);
-			double vawy = -Averages.vaw_global(1);
-			double vawz = -Averages.vaw_global(2);
-		       if (rank==0){
-				// ************* DIMENSIONLESS FORCHEIMER EQUATION *************************
-				//  Dye, A.L., McClure, J.E., Gray, W.G. and C.T. Miller
-				//  Description of Non-Darcy Flows in Porous Medium Systems
-				//  Physical Review E 87 (3), 033012
-				//  Fo := density*D32^3*(density*force) / (viscosity^2)
-				//	Re := density*D32*velocity / viscosity
-				//  Fo = a*Re + b*Re^2
-				// *************************************************************************
-				//viscosity = (tau-0.5)*0.333333333333333333;
-				D32 = 6.0*(Dm.Volume-Averages.vol_w_global)/Averages.As_global;
-				printf("Sauter Mean Diameter = %f \n",D32);
-				mag_force = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
-				Fo = D32*D32*D32*mag_force/viscosity/viscosity;
-				// .... 1-D flow should be aligned with force ...
-				velocity = vawx*Fx/mag_force + vawy*Fy/mag_force + vawz*Fz/mag_force;
-				err1D = fabs(velocity-sqrt(vawx*vawx+vawy*vawy+vawz*vawz))/velocity;
-				//.......... Computation of the Reynolds number Re ..............
-				Re = D32*velocity/viscosity;
-				printf("Force: %.5g,%.5g,%.5g \n",Fx,Fy,Fz);
-				printf("Velocity: %.5g,%.5g,%.5g \n",vawx,vawy,vawz);
-				printf("Relative error for 1D representation: %.5g \n",err1D);
-				printf("Dimensionless force: %5g \n", Fo);
-				printf("Reynolds number: %.5g \n", Re);
-				printf("Dimensionless Permeability (k/D^2): %.5g \n", Re/Fo);
-			}
-			
+		if (nprocs != nprocx*nprocy*nprocz){
+			printf("nprocx =  %i \n",nprocx);
+			printf("nprocy =  %i \n",nprocy);
+			printf("nprocz =  %i \n",nprocz);
+			INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!");
 		}
-	}
-	//************************************************************************/
-	ScaLBL_DeviceBarrier();
-	MPI_Barrier(comm);
-	stoptime = MPI_Wtime();
-	if (rank==0) printf("-------------------------------------------------------------------\n");
-	// Compute the walltime per timestep
-	cputime = (stoptime - starttime)/timestep;
-	// Performance obtained from each node
-	double MLUPS = double(Nx*Ny*Nz)/cputime/1000000;
-	
-	if (rank==0) printf("********************************************************\n");
-	if (rank==0) printf("CPU time = %f \n", cputime);
-	if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
-	MLUPS *= nprocs;
-	if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
-	if (rank==0) printf("********************************************************\n");
 
-	NULL_USE(RESTART_INTERVAL);
+		if (rank==0){
+			printf("********************************************************\n");
+			printf("tau = %f \n", tau);
+			printf("Force(x) = %.5g \n", Fx);
+			printf("Force(y) = %.5g \n", Fy);
+			printf("Force(z) = %.5g \n", Fz);
+			printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
+			printf("Process grid = %i x %i x %i\n",nprocx,nprocy,nprocz);
+			printf("********************************************************\n");
+		}
+
+		double viscosity=(tau-0.5)/3.0;
+		// Initialized domain and averaging framework for Two-Phase Flow
+		int BC=pBC;
+		Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
+		TwoPhase Averages(Dm);
+
+		MPI_Barrier(comm);
+
+		Nx += 2;	Ny += 2;	Nz += 2;
+
+		int N = Nx*Ny*Nz;
+		int dist_mem_size = N*sizeof(double);
+
+		//.......................................................................
+		if (rank == 0)	printf("Read input media... \n");
+		//.......................................................................
+
+		//.......................................................................
+		// Filenames used
+		char LocalRankString[8];
+		char LocalRankFilename[40];
+		char LocalRestartFile[40];
+		char tmpstr[10];
+		sprintf(LocalRankString,"%05d",rank);
+		sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
+		sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
+
+		//	printf("Local File Name =  %s \n",LocalRankFilename);
+		// .......... READ THE INPUT FILE .......................................
+		//	char value;
+		char *id;
+		id = new char[N];
+		int sum = 0;
+		double sum_local;
+		double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
+		if (pBC) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
+		double porosity, pore_vol;
+		//...........................................................................
+		if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
+
+		//.......................................................................
+		sprintf(LocalRankString,"%05d",rank);
+		//	sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
+		//	WriteLocalSolidID(LocalRankFilename, id, N);
+		sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
+		ReadBinaryFile(LocalRankFilename, Averages.SDs.data(), N);
+		MPI_Barrier(comm);
+		if (rank == 0) cout << "Domain set." << endl;
+
+		//.......................................................................
+		// Assign the phase ID field based on the signed distance
+		//.......................................................................
+		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;
+					Dm.id[n] = 0;
+				}
+			}
+		}
+
+		//.......................................................................
+		// Compute the media porosity
+		//.......................................................................
+		double sum,porosity;
+		double sum_local=0.0;	pore_vol = 0.0;
+		int Np=0;  // number of local pore nodes
+		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 (Averages.SDs(n) > 0.0){
+						Dm.id[n] = 2;	
+					}
+					// compute the porosity (actual interface location used)
+					if (Averages.SDs(n) > 0.0){
+						sum_local+=1.0;
+						Np++;
+					}
+				}
+			}
+		}
+		MPI_Barrier(comm);
+		MPI_Allreduce(&sum_local,&sum,1,MPI_DOUBLE,MPI_SUM,comm);
+		porosity = sum*iVol_global;
+		if (rank==0) printf("Media porosity = %f \n",porosity);
+
+		MPI_Barrier(comm);
+		if (rank == 0) cout << "Domain set." << endl;
+		if (rank==0)	printf ("Create ScaLBL_Communicator \n");
+		//...........................................................................
+		if (rank==0)	printf ("Create ScaLBL_Communicator \n");
+		Dm.CommInit(comm);
+		// Create a communicator for the device
+		ScaLBL_Communicator ScaLBL_Comm(Dm);
+
+		// LBM variables
+		if (rank==0)	printf ("Allocating distributions \n");
+
+		int neighborSize=18*Np*sizeof(int);
+		int *neighborList;
+		IntArray Map(Nx,Ny,Nz);
+
+		neighborList= new int[18*Np];
+		ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np);
+		// ScaLBL_Comm.MemoryDenseLayoutFull(Map,neighborList,Dm.id,Np);    // this was how I tested for correctness
+
+		MPI_Barrier(comm);
+
+		//......................device distributions.................................
+		int dist_mem_size = Np*sizeof(double);
+
+		int *NeighborList;
+		//		double *f_even,*f_odd;
+		double * dist;
+		double * Velocity;
+		//...........................................................................
+		ScaLBL_AllocateDeviceMemory((void **) &dist, 19*dist_mem_size);
+		ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
+		ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
+		ScaLBL_AllocateDeviceMemory((void **) &Pressure, 3*sizeof(double)*Np);
+		ScaLBL_CopyToDevice(NeighborList,     neighborList, neighborSize);
+		//...........................................................................
+
+		//...........................................................................
+		if (rank==0)	printf("Setting the distributions, size = %i\n", N);
+		//...........................................................................
+
+		// Finalize setup for averaging domain
+		//Averages.SetupCubes(Dm);
+		Averages.UpdateSolid();
+		// Initialize two phase flow variables (all wetting phase)
+		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;
+					Averages.Phase(i,j,k) = -1.0;
+					Averages.SDn(i,j,k) = Averages.Phase(i,j,k);
+					Averages.Phase_tplus(i,j,k) = Averages.SDn(i,j,k);
+					Averages.Phase_tminus(i,j,k) = Averages.SDn(i,j,k);
+					Averages.DelPhi(i,j,k) = 0.0;
+					Averages.Press(i,j,k) = 0.0;
+					Averages.Vel_x(i,j,k) = 0.0;
+					Averages.Vel_y(i,j,k) = 0.0;
+					Averages.Vel_z(i,j,k) = 0.0;
+				}
+			}
+		}
+
+		//.......................................................................
+
+		ScaLBL_D3Q19_Init(dist, Np);
+
+		int timestep = 0;
+		if (rank==0) printf("********************************************************\n");
+		if (rank==0)	printf("No. of timesteps: %i \n", timestepMax);
+
+		//.......create and start timer............
+		double starttime,stoptime,cputime;
+		MPI_Barrier(comm);
+		starttime = MPI_Wtime();
+		//.........................................
+
+		double D32,Fo,Re,velocity,err1D,mag_force,vel_prev;
+		err = vel_prev = 1.0;
+		if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol);
+		//************ MAIN ITERATION LOOP ***************************************/
+		while (timestep < timestepMax && err > tol ){
+
+			timestep++;
+			ScaLBL_Comm.SendD3Q19AA(dist); //READ FROM NORMAL
+			ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, ScaLBL_Comm.next, Np, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
+			ScaLBL_Comm.RecvD3Q19AA(dist); //WRITE INTO OPPOSITE
+			ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, 0, ScaLBL_Comm.next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
+			ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+
+			timestep++;
+			ScaLBL_Comm.SendD3Q19AA(dist); //READ FORM NORMAL
+			ScaLBL_D3Q19_AAeven_MRT(dist, ScaLBL_Comm.next, Np, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
+			ScaLBL_Comm.RecvD3Q19AA(dist); //WRITE INTO OPPOSITE
+			ScaLBL_D3Q19_AAeven_MRT(dist, 0, ScaLBL_Comm.next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
+			ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+			//************************************************************************/
+
+			if (timestep%500 == 0){
+				//...........................................................................
+				// Copy the data for for the analysis timestep
+				//...........................................................................
+				// Copy the phase from the GPU -> CPU
+				//...........................................................................
+				ScaLBL_DeviceBarrier();
+		        ScaLBL_D3Q19_Pressure(fq,Pressure,Np);
+		        ScaLBL_D3Q19_Momentum(fq,Velocity,Np);
+		        
+		        ScaLBL_Comm.RegularLayout(Map,Pressure,Averages.Press);
+		        ScaLBL_Comm.RegularLayout(Map,&Velocity[0],Averages.Vel_x);
+		        ScaLBL_Comm.RegularLayout(Map,&Velocity[Np],Averages.Vel_y);
+		        ScaLBL_Comm.RegularLayout(Map,&Velocity[2*Np],Averages.Vel_z);
+
+				// Way more work than necessary -- this is just to get the solid interfacial area!!
+				Averages.Initialize();
+				Averages.UpdateMeshValues();
+				Averages.ComputeLocal();
+				Averages.Reduce();
+
+				double vawx = -Averages.vaw_global(0);
+				double vawy = -Averages.vaw_global(1);
+				double vawz = -Averages.vaw_global(2);
+				if (rank==0){
+					// ************* DIMENSIONLESS FORCHEIMER EQUATION *************************
+					//  Dye, A.L., McClure, J.E., Gray, W.G. and C.T. Miller
+					//  Description of Non-Darcy Flows in Porous Medium Systems
+					//  Physical Review E 87 (3), 033012
+					//  Fo := density*D32^3*(density*force) / (viscosity^2)
+					//	Re := density*D32*velocity / viscosity
+					//  Fo = a*Re + b*Re^2
+					// *************************************************************************
+					//viscosity = (tau-0.5)*0.333333333333333333;
+					D32 = 6.0*(Dm.Volume-Averages.vol_w_global)/Averages.As_global;
+					printf("Sauter Mean Diameter = %f \n",D32);
+					mag_force = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
+					Fo = D32*D32*D32*mag_force/viscosity/viscosity;
+					// .... 1-D flow should be aligned with force ...
+					velocity = vawx*Fx/mag_force + vawy*Fy/mag_force + vawz*Fz/mag_force;
+					err1D = fabs(velocity-sqrt(vawx*vawx+vawy*vawy+vawz*vawz))/velocity;
+					//.......... Computation of the Reynolds number Re ..............
+					Re = D32*velocity/viscosity;
+					printf("Force: %.5g,%.5g,%.5g \n",Fx,Fy,Fz);
+					printf("Velocity: %.5g,%.5g,%.5g \n",vawx,vawy,vawz);
+					printf("Relative error for 1D representation: %.5g \n",err1D);
+					printf("Dimensionless force: %5g \n", Fo);
+					printf("Reynolds number: %.5g \n", Re);
+					printf("Dimensionless Permeability (k/D^2): %.5g \n", Re/Fo);
+				}
+
+			}
+		}
+		//************************************************************************/
+		ScaLBL_DeviceBarrier();
+		MPI_Barrier(comm);
+		stoptime = MPI_Wtime();
+		if (rank==0) printf("-------------------------------------------------------------------\n");
+		// Compute the walltime per timestep
+		cputime = (stoptime - starttime)/timestep;
+		// Performance obtained from each node
+		double MLUPS = double(Nx*Ny*Nz)/cputime/1000000;
+
+		if (rank==0) printf("********************************************************\n");
+		if (rank==0) printf("CPU time = %f \n", cputime);
+		if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
+		MLUPS *= nprocs;
+		if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
+		if (rank==0) printf("********************************************************\n");
+
+		NULL_USE(RESTART_INTERVAL);
 	}
 	// ****************************************************
 	MPI_Barrier(comm);