From c5477b9b39bcb5bb99f8ea3bafa53e892808156b Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Fri, 26 Jan 2018 08:50:25 -0500
Subject: [PATCH] Refactoring color simulator

---
 tests/lbpm_color_macro_simulator.cpp | 1610 ++++++++++----------------
 1 file changed, 627 insertions(+), 983 deletions(-)

diff --git a/tests/lbpm_color_macro_simulator.cpp b/tests/lbpm_color_macro_simulator.cpp
index bbc37c88..f42d1f1a 100644
--- a/tests/lbpm_color_macro_simulator.cpp
+++ b/tests/lbpm_color_macro_simulator.cpp
@@ -16,1041 +16,685 @@
 
 /*
  * Simulator for two-phase flow in porous media
- * James E. McClure 2013-2016
+ * James E. McClure 2013-2018
  */
 
 using namespace std;
 
 //*************************************************************************
-// Implementation of Two-Phase Immiscible LBM using CUDA
+// Implementation of Two-Phase Immiscible LBM 
 //*************************************************************************
-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)
 {
   // Initialize MPI
-  int provided_thread_support = -1;
-  MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
-  MPI_Comm comm;
-  MPI_Comm_dup(MPI_COMM_WORLD,&comm);
-  int rank = comm_rank(comm);
-  int nprocs = comm_size(comm);
-  if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE )
-    std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
-  { // Limit scope so variables that contain communicators will free before MPI_Finialize
+	int provided_thread_support = -1;
+	MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
+	MPI_Comm comm;
+	MPI_Comm_dup(MPI_COMM_WORLD,&comm);
+	int rank = comm_rank(comm);
+	int nprocs = comm_size(comm);
+	if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE )
+		std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
+	{ // Limit scope so variables that contain communicators will free before MPI_Finialize
 
-	// parallel domain size (# of sub-domains)
-	int nprocx,nprocy,nprocz;
-	int iproc,jproc,kproc;
+		// parallel domain size (# of sub-domains)
+		int nprocx,nprocy,nprocz;
+		int iproc,jproc,kproc;
 
-	MPI_Request req1[18],req2[18];
-	MPI_Status stat1[18],stat2[18];
+		MPI_Request req1[18],req2[18];
+		MPI_Status stat1[18],stat2[18];
 
-	if (rank == 0){
-		printf("********************************************************\n");
-		printf("Running Color LBM	\n");
-		printf("********************************************************\n");
-	}
+		if (rank == 0){
+			printf("********************************************************\n");
+			printf("Running Color LBM	\n");
+			printf("********************************************************\n");
+		}
 
-    PROFILE_ENABLE(1);
-    //PROFILE_ENABLE_TRACE();
-    //PROFILE_ENABLE_MEMORY();
-    PROFILE_SYNCHRONIZE();
-    PROFILE_START("Main");
-    Utilities::setErrorHandlers();
+		PROFILE_ENABLE(1);
+		//PROFILE_ENABLE_TRACE();
+		//PROFILE_ENABLE_MEMORY();
+		PROFILE_SYNCHRONIZE();
+		PROFILE_START("Main");
+		Utilities::setErrorHandlers();
 
-	// Variables that specify the computational domain  
-	string FILENAME;
-	unsigned int nBlocks, nthreads;
-	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;
-	double tau1, tau2, rho1,rho2;
-	double Fx,Fy,Fz,tol,err;
-	double alpha, beta;
-	double das, dbs, phi_s;
-	double din,dout;
-	double wp_saturation;
-	int BoundaryCondition;
-	int InitialCondition;
-//	bool pBC,Restart;
-	int i,j,k;
+		// Variables that specify the computational domain  
+		string FILENAME;
+		int Nx,Ny,Nz;		// local sub-domain size
+		double Lx,Ly,Lz;	// Domain length
+		double D = 1.0;		// reference length for non-dimensionalization
+		// Color Model parameters
+		int timestepMax;
+		double tau1, tau2, rho1,rho2;
+		double Fx,Fy,Fz,tol,err;
+		double alpha, beta;
+		int BoundaryCondition;
+		int InitialCondition;
+		//	bool pBC,Restart;
+		int i,j,k;
+		double din, dout, flux;
+		double inletA,inletB,outletA,outletB;
+		inletA=1.f;
+		inletB=0.f;
+		outletA=0.f;
+		outletB=1.f;
+		typeBC=4;
+		flux = 10.f;
+		dout=1.f;
 
-	// pmmc threshold values
-	//double fluid_isovalue,solid_isovalue;
-	//fluid_isovalue = 0.0;
-	//solid_isovalue = 0.0;
-	
-	int RESTART_INTERVAL=20000;
-	//int ANALYSIS_)INTERVAL=1000;	
-	int BLOB_ANALYSIS_INTERVAL=1000;
-	int timestep = -1;
+		int RESTART_INTERVAL=20000;
+		//int ANALYSIS_)INTERVAL=1000;	
+		int BLOB_ANALYSIS_INTERVAL=1000;
+		int timestep = -1;
 
-	if (rank==0){
-		//.............................................................
-		//		READ SIMULATION PARMAETERS FROM INPUT FILE
-		//.............................................................
-		ifstream input("Color.in");
-		if (input.is_open()){
-			// Line 1: model parameters (tau, alpha, beta, das, dbs)
-			input >> tau1;			// Viscosity non-wetting
-			input >> tau2;			// Viscosity wetting
-			input >> rho1;			// density non-wetting
-			input >> rho2;			// density wetting
-			input >> alpha;			// Surface Tension parameter
-			input >> beta;			// Width of the interface
-			input >> phi_s;			// value of phi at the solid surface
-			// Line 2: wetting phase saturation to initialize
-			input >> wp_saturation;
-			// Line 3: External force components (Fx,Fy, Fz)
-			input >> Fx;
-			input >> Fy;
-			input >> Fz;
-			// Line 4: Pressure Boundary conditions
-			input >> InitialCondition;
-			input >> BoundaryCondition;
-			input >> din;
-			input >> dout;
-			// Line 5: time-stepping criteria
-			input >> timestepMax;		// max no. of timesteps
-			input >> RESTART_INTERVAL;	// restart interval
-			input >> tol;		      	// error tolerance
-			// Line 6: Analysis options
-			input >> BLOB_ANALYSIS_INTERVAL; // interval to analyze blob states
+		if (rank==0){
 			//.............................................................
-		}
-		else{
-			// Set default values
-		    // Print warning
-			printf("WARNING: No input file provided (Color.in is missing)! Default parameters will be used. \n");
-			tau1 = tau2 = 1.0;
-			rho1 = rho2 = 1.0;
-			alpha=0.005;
-			beta= 0.9;
-			Fx = Fy = Fz = 0.0;
-			InitialCondition=0;
-			BoundaryCondition=0;
-			din=dout=1.0;
-			timestepMax=0;
-		}
+			//		READ SIMULATION PARMAETERS FROM INPUT FILE
+			//.............................................................
+			ifstream input("Color.in");
+			if (input.is_open()){
+				// Line 1: model parameters (tau, alpha, beta, das, dbs)
+				input >> tau1;			// Viscosity non-wetting
+				input >> tau2;			// Viscosity wetting
+				input >> rho1;			// density non-wetting
+				input >> rho2;			// density wetting
+				input >> alpha;			// Surface Tension parameter
+				input >> beta;			// Width of the interface
+				// Line 2:  External force components (Fx,Fy, Fz)
+				input >> Fx;
+				input >> Fy;
+				input >> Fz;
+				// Line 4: Pressure Boundary conditions
+				input >> InitialCondition;
+				input >> BoundaryCondition;
+				input >> din;
+				input >> dout;
+				// Line 5: time-stepping criteria
+				input >> timestepMax;		// max no. of timesteps
+				input >> RESTART_INTERVAL;	// restart interval
+				input >> tol;		      	// error tolerance
+				// Line 6: Analysis options
+				input >> BLOB_ANALYSIS_INTERVAL; // interval to analyze blob states
+				//.............................................................
+			}
+			else{
+				// Set default values
+				// Print warning
+				printf("WARNING: No input file provided (Color.in is missing)! Default parameters will be used. \n");
+				tau1 = tau2 = 1.0;
+				rho1 = rho2 = 1.0;
+				alpha=0.005;
+				beta= 0.9;
+				Fx = Fy = Fz = 0.0;
+				InitialCondition=0;
+				BoundaryCondition=0;
+				din=dout=1.0;
+				timestepMax=0;
+			}
 
-		//.......................................................................
-		// Reading the domain information file
-		//.......................................................................
-		ifstream domain("Domain.in");
-		if (input.is_open()){
-
-			domain >> nprocx;
-			domain >> nprocy;
-			domain >> nprocz;
-			domain >> Nx;
-			domain >> Ny;
-			domain >> Nz;
-			domain >> nspheres;
-			domain >> Lx;
-			domain >> Ly;
-			domain >> Lz;
 			//.......................................................................
-		}
-		else{
-			// Set default values
-		    // Print warning
-			printf("WARNING: No input file provided (Domain.in is missing)! Default parameters will be used. \n");
-			nprocx=nprocy=nprocz=1;
-			Nx=Ny=Nz=10;
-			nspheres=0;
-			Lx=Ly=Lz=1.0;
-		}
-	}
-	// **************************************************************
-	// Broadcast simulation parameters from rank 0 to all other procs
-	MPI_Barrier(comm);
-	//.................................................
-	MPI_Bcast(&tau1,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&tau2,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&rho1,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&rho2,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&alpha,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&beta,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&das,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&dbs,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&phi_s,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&wp_saturation,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&BoundaryCondition,1,MPI_INT,0,comm);
-	MPI_Bcast(&InitialCondition,1,MPI_INT,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(&RESTART_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);
-	//.................................................
-	
-	double flux = 0.f;
-	if (BoundaryCondition==4) flux = din*rho1; // mass flux must adjust for density (see formulation for details)
-
-    // Get the rank info
-    const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
-
-	MPI_Barrier(comm);
-	
-	// **************************************************************
-	// **************************************************************
-	double Ps = -(das-dbs)/(das+dbs);
-	//double xIntPos = log((1.0+phi_s)/(1.0-phi_s))/(2.0*beta); 	
-	
-	// Set the density values inside the solid based on the input value phi_s
- 	das = (phi_s+1.0)*0.5;
-	dbs = 1.0 - das;
-	
-	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 (non-wetting) = %f \n", tau1);
-		printf("tau (wetting) = %f \n", tau2);
-		printf("density (non-wetting) = %f \n", rho1);
-		printf("density (wetting) = %f \n", rho2);
-		printf("alpha = %f \n", alpha);		
-		printf("beta = %f \n", beta);
-		printf("das = %f \n", das);
-		printf("dbs = %f \n", dbs);
-		printf("gamma_{wn} = %f \n", 5.796*alpha);
-		printf("Force(x) = %f \n", Fx);
-		printf("Force(y) = %f \n", Fy);
-		printf("Force(z) = %f \n", Fz);
-		printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
-		printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
-		if (BoundaryCondition==0) printf("Periodic boundary conditions will applied \n");
-		if (BoundaryCondition==1) printf("Pressure boundary conditions will be applied \n");
-		if (BoundaryCondition==2) printf("Velocity boundary conditions will be applied \n");
-		if (BoundaryCondition==3) printf("Dynamic pressure boundary conditions will be applied \n");
-		if (BoundaryCondition==4) printf("Average flux boundary conditions will be applied \n");
-		if (InitialCondition==0) printf("Initial conditions assigned from phase ID file \n");
-		if (InitialCondition==1) printf("Initial conditions assigned from restart file \n");
-		printf("********************************************************\n");
-	}
-
-	// Initialized domain and averaging framework for Two-Phase Flow
-	bool pBC,velBC;
-	if (BoundaryCondition==1 || BoundaryCondition==3 || BoundaryCondition == 4)
-								pBC=true;
-	else						pBC=false;
-	if (BoundaryCondition==2)	velBC=true;
-	else						velBC=false;
-
-	bool Restart;
-	if (InitialCondition==1)    Restart=true;
-	else 						Restart=false;
-    NULL_USE(pBC);  NULL_USE(velBC);
-
-	// Full domain used for averaging (do not use mask for analysis)
-	Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
-	for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
-	std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
-    //   TwoPhase Averages(Dm);
-	Dm.CommInit(comm);
-
-	// Mask that excludes the solid phase
-	Domain Mask(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
-
-	MPI_Barrier(comm);
-
-	Nx+=2; Ny+=2; Nz += 2;
-	//Nx = Ny = Nz;	// Cubic domain
-
-	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];
-	double sum_local;
-	double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
-	if (BoundaryCondition > 0) 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);
-	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++){
-				int n = k*Nx*Ny+j*Nx+i;
-				id[n] = 0;
+			// Reading the domain information file
+			//.......................................................................
+			ifstream domain("Domain.in");
+			if (input.is_open()){
+				domain >> nprocx;
+				domain >> nprocy;
+				domain >> nprocz;
+				domain >> Nx;
+				domain >> Ny;
+				domain >> Nz;
+				domain >> nspheres;
+				domain >> Lx;
+				domain >> Ly;
+				domain >> Lz;
+				//.......................................................................
+			}
+			else{
+				// Set default values
+				// Print warning
+				printf("WARNING: No input file provided (Domain.in is missing)! Default parameters will be used. \n");
+				nprocx=nprocy=nprocz=1;
+				Nx=Ny=Nz=10;
+				nspheres=0;
+				Lx=Ly=Lz=1.0;
 			}
 		}
-	}
-	double sum=0;
-	pore_vol = 0.0;
-	for ( k=0;k<Nz;k++){
-		for ( j=0;j<Ny;j++){
-			for ( i=0;i<Nx;i++){
-				int 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++;	
-				}
-			}
+		// **************************************************************
+		// Broadcast simulation parameters from rank 0 to all other procs
+		MPI_Barrier(comm);
+		//.................................................
+		MPI_Bcast(&tau1,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&tau2,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&rho1,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&rho2,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&alpha,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&beta,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&BoundaryCondition,1,MPI_INT,0,comm);
+		MPI_Bcast(&InitialCondition,1,MPI_INT,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(&RESTART_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);
+		//.................................................
+
+		double flux = 0.f;
+		if (BoundaryCondition==4) flux = din*rho1; // mass flux must adjust for density (see formulation for details)
+
+		// Get the rank info
+		const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
+
+		MPI_Barrier(comm);
+
+		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("Initialize from segmented data: solid=0, NWP=1, WP=2 \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_segmented_pp: Error reading ID (rank=%i) \n",rank);
-	fclose(IDFILE);
-
-    // Calculate the actual inlet cross-sectional area
-	double inlet_area_local=0.f;
-	double InletArea;
-    if (Dm.kproc==0)
-    {
-        for ( k=1;k<2;k++){
-            for ( j=1;j<Ny-1;j++){
-                for ( i=1;i<Nx-1;i++){
-                    int n = k*Nx*Ny+j*Nx+i;
-                    if (id[n] > 0){
-                        inlet_area_local+=1.f;
-                    }
-                }
-            }
-        }
-    
-    }
-	MPI_Allreduce(&inlet_area_local,&InletArea,1,MPI_DOUBLE,MPI_SUM,comm);
-
-	// Set up kstart, kfinish so that the reservoirs are excluded from averaging
-	int kstart,kfinish;
-	kstart = 1;
-	kfinish = Nz-1;
-	if (BoundaryCondition >  0 && Dm.kproc==0)		kstart = 4;
-	if (BoundaryCondition >  0 && Dm.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++){
-				int n = k*Nx*Ny+j*Nx+i;
-				if (id[n] > 0){
-					sum_local += 1.0;
-				}
-			}
+		if (rank==0){
+			printf("********************************************************\n");
+			printf("tau (non-wetting) = %f \n", tau1);
+			printf("tau (wetting) = %f \n", tau2);
+			printf("density (non-wetting) = %f \n", rho1);
+			printf("density (wetting) = %f \n", rho2);
+			printf("alpha = %f \n", alpha);		
+			printf("beta = %f \n", beta);
+			printf("gamma_{wn} = %f \n", 5.796*alpha);
+			printf("Force(x) = %f \n", Fx);
+			printf("Force(y) = %f \n", Fy);
+			printf("Force(z) = %f \n", Fz);
+			printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
+			printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
+			if (BoundaryCondition==0) printf("Periodic boundary conditions will applied \n");
+			if (BoundaryCondition==1) printf("Pressure boundary conditions will be applied \n");
+			if (BoundaryCondition==2) printf("Velocity boundary conditions will be applied \n");
+			if (BoundaryCondition==3) printf("Dynamic pressure boundary conditions will be applied \n");
+			if (BoundaryCondition==4) printf("Average flux boundary conditions will be applied \n");
+			if (InitialCondition==0) printf("Initial conditions assigned from phase ID file \n");
+			if (InitialCondition==1) printf("Initial conditions assigned from restart file \n");
+			printf("********************************************************\n");
 		}
-	}
 
-	MPI_Allreduce(&sum_local,&pore_vol,1,MPI_DOUBLE,MPI_SUM,comm);
-	porosity = pore_vol*iVol_global;
-	if (rank==0) printf("Media porosity = %f \n",porosity);
-	//.........................................................
-	/*	// If external boundary conditions are applied remove solid
-	if (BoundaryCondition >  0  && Dm.kproc == 0){
-	  // zero out the first layer so that no comms are performed on that boundary
-	  for (j=0; j<Ny; j++){
-	    for (i=0; i<Nx; i++){
-	      int n = j*Nx+i;
-	      id[n] = 0;
-	    }
-	  }
-		for (k=0; k<3; k++){
+		// Initialized domain and averaging framework for Two-Phase Flow
+		bool pBC,velBC;
+		if (BoundaryCondition==1 || BoundaryCondition==3 || BoundaryCondition == 4)
+			pBC=true;
+		else						pBC=false;
+		if (BoundaryCondition==2)	velBC=true;
+		else						velBC=false;
+
+		bool Restart;
+		if (InitialCondition==1)    Restart=true;
+		else 						Restart=false;
+		NULL_USE(pBC);  NULL_USE(velBC);
+
+		// Full domain used for averaging (do not use mask for analysis)
+		Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
+		for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
+		std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
+		//   TwoPhase Averages(Dm);
+		Dm.CommInit(comm);
+
+		// Mask that excludes the solid phase
+		Domain Mask(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
+		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];
+		double sum_local;
+		double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
+		if (BoundaryCondition > 0) 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;
+
+		//.......................................................................
+		// Read the signed distance
+		sprintf(LocalRankString,"%05d",rank);
+		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++){
 					int n = k*Nx*Ny+j*Nx+i;
-					//id[n] = 1;
-					//Averages.SDs(n) = max(Averages.SDs(n),1.0*(2.5-k));
-					Averages->SDs(n) = max(Averages->SDs(n),1.0*(2.5-k));
-
-				}					
+					id[n] = 0;
+				}
 			}
 		}
-	}
-	if (BoundaryCondition >  0  && Dm.kproc == nprocz-1){
-	  // zero out the last layer so that no comms are performed on that boundary
-	  for (j=0; j<Ny; j++){
-	    for (i=0; i<Nx; i++){
-	      int n = (Nz-1)*Nx*Ny+j*Nx+i;
-	      id[n] = 0;
-	    }
-	  }
-		for (k=Nz-3; k<Nz; k++){
-			for (j=0;j<Ny;j++){
-				for (i=0;i<Nx;i++){
-					int n = k*Nx*Ny+j*Nx+i;
-					//id[n] = 2;
-					//Averages.SDs(n) = max(Averages.SDs(n),1.0*(k-Nz+2.5));
-					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;
-	//.........................................................
-
-	// set reservoirs
-	if (BoundaryCondition > 0){
+		double sum=0;
+		pore_vol = 0.0;
 		for ( k=0;k<Nz;k++){
 			for ( j=0;j<Ny;j++){
 				for ( i=0;i<Nx;i++){
 					int n = k*Nx*Ny+j*Nx+i;
-					if (id[n] > 0){
-						if (Dm.kproc==0 && k==0)			id[n]=1;
-						if (Dm.kproc==0 && k==1)			id[n]=1;
-						if (Dm.kproc==nprocz-1 && k==Nz-2)	id[n]=2;
-						if (Dm.kproc==nprocz-1 && k==Nz-1)	id[n]=2;
+					if (Averages->SDs(n) > 0.0){
+						id[n] = 2;	
+					}
+					// compute the porosity (actual interface location used)
+					if (Averages->SDs(n) > 0.0){
+						sum++;	
 					}
-					Mask.id[n] = id[n];
 				}
 			}
 		}
-	}
 
-	// Initialize communication structures in averaging domain
-	for (i=0; i<Mask.Nx*Mask.Ny*Mask.Nz; i++) Mask.id[i] = id[i];
-	Mask.CommInit(comm);
-
-	//...........................................................................
-	if (rank==0)	printf ("Create ScaLBL_Communicator \n");
-	// Create a communicator for the device
-	ScaLBL_Communicator ScaLBL_Comm(Mask);
-
-	//...........device phase ID.................................................
-	if (rank==0)	printf ("Copy phase ID to device \n");
-	char *ID;
-	ScaLBL_AllocateDeviceMemory((void **) &ID, N);						// Allocate device memory
-	// Don't compute in the halo
-	for (k=0;k<Nz;k++){
-		for (j=0;j<Ny;j++){
-			for (i=0;i<Nx;i++){
-				int n = k*Nx*Ny+j*Nx+i;
-				if (i==0 || i==Nx-1 || j==0 || j==Ny-1 || k==0 || k==Nz-1)	id[n] = 0;
-			}
-		}
-	}
-	// Copy to the device
-	ScaLBL_CopyToDevice(ID, id, N);
-	ScaLBL_DeviceBarrier();
-	//...........................................................................
-
-	//...........................................................................
-	//				MAIN  VARIABLES ALLOCATED HERE
-	//...........................................................................
-	// LBM variables
-	if (rank==0)	printf ("Allocating distributions \n");
-	//......................device distributions.................................
-	double *f_even,*f_odd;
-	double *A_even,*A_odd,*B_even,*B_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
-	ScaLBL_AllocateDeviceMemory((void **) &A_even, 4*dist_mem_size);	// Allocate device memory
-	ScaLBL_AllocateDeviceMemory((void **) &A_odd, 3*dist_mem_size);	// Allocate device memory
-	ScaLBL_AllocateDeviceMemory((void **) &B_even, 4*dist_mem_size);	// Allocate device memory
-	ScaLBL_AllocateDeviceMemory((void **) &B_odd, 3*dist_mem_size);	// Allocate device memory
-	//...........................................................................
-	double *Phi,*Den;
-	double *ColorGrad, *Velocity, *Pressure, *dvcSignDist;
-	//...........................................................................
-	ScaLBL_AllocateDeviceMemory((void **) &Phi, dist_mem_size);
-	ScaLBL_AllocateDeviceMemory((void **) &Pressure, dist_mem_size);
-	ScaLBL_AllocateDeviceMemory((void **) &dvcSignDist, dist_mem_size);
-	ScaLBL_AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
-	ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
-	ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size);
-	//...........................................................................
-
-	// Copy signed distance for device initialization
-	//ScaLBL_CopyToDevice(dvcSignDist, Averages.SDs.data(), dist_mem_size);
-	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_DeviceBarrier();
-	ScaLBL_D3Q19_Init(ID, f_even, f_odd, Nx, Ny, Nz);
-	ScaLBL_Color_Init(ID, Den, Phi, das, dbs, Nx, Ny, Nz);
-	ScaLBL_DeviceBarrier();
-	//......................................................................
-
-	if (Restart == true){
-
-	  
-	  if (rank==0){
-	    printf("Reading restart file! \n");
-	    ifstream restart("Restart.txt");
-	    if (restart.is_open()){
-	      restart  >> timestep;
-	      printf("Restarting from timestep =%i \n",timestep);
-	    }
-	    else{
-	      printf("WARNING:No Restart.txt file, setting timestep=0 \n");
-	      timestep=0;
-	    }
-	  }
-	  MPI_Bcast(&timestep,1,MPI_INT,0,comm);
-
-		// Read in the restart file to CPU buffers
-	    double *cDen = new double[2*N];
-	    double *cDistEven = new double[10*N];
-	    double *cDistOdd = new double[9*N];
-		ReadCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
-		// Copy the restart data to the GPU
-		ScaLBL_CopyToDevice(f_even,cDistEven,10*N*sizeof(double));
-		ScaLBL_CopyToDevice(f_odd,cDistOdd,9*N*sizeof(double));
-		ScaLBL_CopyToDevice(Den,cDen,2*N*sizeof(double));
-		ScaLBL_DeviceBarrier();
-	    delete [] cDen;
-	    delete [] cDistEven;
-	    delete [] cDistOdd;
-		MPI_Barrier(comm);
-	}
-
-	//......................................................................
-	ScaLBL_D3Q7_Init(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
-	ScaLBL_D3Q7_Init(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
-	ScaLBL_DeviceBarrier();
-	MPI_Barrier(comm);
-	//.......................................................................
-	// Once phase has been initialized, map solid to account for 'smeared' interface
-	//for (i=0; i<N; i++)	Averages.SDs(i) -= (1.0);
-	// Make sure the id match for the two domains
-	for (i=0; i<N; i++)	Dm.id[i] = Mask.id[i];
-	//.......................................................................
-	// Finalize setup for averaging domain
-	//Averages.SetupCubes(Mask);
-//	Averages.UpdateSolid();
-	Averages->UpdateSolid();
-
-	//.......................................................................
-	
-	//*************************************************************************
-	// 		Compute the phase indicator field and reset Copy, Den
-	//*************************************************************************
-	ScaLBL_ComputePhaseField(ID, Phi, Den, N);
-	//*************************************************************************
-	ScaLBL_DeviceBarrier();
-	ScaLBL_Comm.SendHalo(Phi);
-	ScaLBL_Comm.RecvHalo(Phi);
-	ScaLBL_DeviceBarrier();
-	MPI_Barrier(comm);
-	//*************************************************************************
-
-	if (rank==0 && BoundaryCondition==1){
-		printf("Setting inlet pressure = %f \n", din);
-		printf("Setting outlet pressure = %f \n", dout);
-	}
-	if (BoundaryCondition==1 && Mask.kproc == 0)	{
-		ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-    	ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-		ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
-	}
+		if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \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_segmented_pp: Error reading ID (rank=%i) \n",rank);
+		fclose(IDFILE);
 		
-	if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
-		ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-    	ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-		ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
-	}
-
-	if (rank==0 && BoundaryCondition==2){
-		printf("Setting inlet velocity = %f \n", din);
-		printf("Setting outlet velocity = %f \n", dout);
-	}
-	if (BoundaryCondition==2 && Mask.kproc == 0)	{
-		ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-    	ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-		ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
-	}
-
-	if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
-		ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-    	ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-		ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
-	}
-
-	// Set dynamic pressure boundary conditions
-	double dp, slope;
-	dp = slope = 0.0;
-	if (BoundaryCondition==3){
-		slope = abs(dout-din)/timestepMax;
-		dp = abs(timestep)*slope;
-		if (rank==0) printf("Change in pressure / time =%.3e \n",slope);
-		// set the initial value
-		din  = 1.0+0.5*dp;
-		dout = 1.0-0.5*dp;
-		// set the initial boundary conditions
-		if (Mask.kproc == 0)	{
-			ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-        	ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-			ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
-
+		//.......................................................................
+		// Compute the media porosity, assign phase labels and solid composition
+		//.......................................................................
+		double sum;
+		double sum_local=0.0, porosity;
+		int Np=0;  // number of local pore nodes
+		double *PhaseLabel;
+		PhaseLabel = new double[N];
+		Dm.AssignComponentLabels(PhaseLabel);
+		//.......................................................................
+		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 (id[n] > 0){
+						sum_local+=1.0;
+						Np++;
+					}
+				}
+			}
 		}
-		if (Mask.kproc == nprocz-1){
-			ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-        	ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-			ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
-		}
-	}
+		MPI_Allreduce(&sum_local,&sum,1,MPI_DOUBLE,MPI_SUM,comm);
+		porosity = sum*iVol_global;
+		if (rank==0) printf("Media porosity = %f \n",porosity);
 
-	// Set flux boundary condition
-	double SumLocal,SumGlobal;
-	if (BoundaryCondition==4){
-	  double SumLocal=0.f;
-	  if (rank==0){
-	    printf("Using flux boundary condition \n");
-	    printf("   Volumetric flux: Q = %f \n",flux);
-	    printf("   Inlet area  = %f \n",InletArea);
-	    printf("   outlet pressure: %f \n",dout);
+		//.........................................................
+		// 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;
+		//.........................................................
 
-	  }
-	  SumLocal=0.0;
-	  if (pBC && Dm.kproc == 0){
-	     SumLocal = ScaLBL_D3Q19_Flux_BC_z(ID,f_even,f_odd,flux,InletArea,Nx,Ny,Nz);
-	  }
-	  MPI_Allreduce(&SumLocal,&SumGlobal,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
-	  din = flux / InletArea + SumGlobal;
-	  
-	  if (pBC && Dm.kproc == 0){
-	    if (rank==0) printf("Flux = %.3e, Computed inlet pressure: %f \n",flux,din);
-	    ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-	    ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-		ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
-		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
-		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
-	  }
+		// Initialize communication structures in averaging domain
+		for (i=0; i<Mask.Nx*Mask.Ny*Mask.Nz; i++) Mask.id[i] = id[i];
+		Mask.CommInit(comm);
 
-	  if (pBC && Dm.kproc == nprocz-1){
-	   // if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout);
-	    ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-	    ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-		ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
-		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
-		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
-	  }
-	}
-
-	ScaLBL_D3Q19_Pressure(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
-	ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-
-	//...........................................................................
-	// Copy the phase indicator field for the earlier timestep
-	ScaLBL_DeviceBarrier();
-	//ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double));
-	ScaLBL_CopyToHost(Averages->Phase_tplus.data(),Phi,N*sizeof(double));
-
-	//...........................................................................
-	//...........................................................................
-	// 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_CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
-	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));
-	//...........................................................................
-	
-	if (rank==0) printf("********************************************************\n");
-	if (rank==0)	printf("No. of timesteps: %i \n", timestepMax);
-
-	//.......create and start timer............
-	double starttime,stoptime,cputime;
-	ScaLBL_DeviceBarrier();
-	MPI_Barrier(comm);
-	starttime = MPI_Wtime();
-	//.........................................
-	
-	err = 1.0; 	
-	double sat_w_previous = 1.01; // slightly impossible value!
-	if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol);
-    // Create the thread pool
-    int N_threads = 4;
-    if ( provided_thread_support < MPI_THREAD_MULTIPLE )
-        N_threads = 0;
-    if ( N_threads > 0 ) {
-        // Set the affinity
-        int N_procs = ThreadPool::getNumberOfProcessors();
-        std::vector<int> procs(N_procs);
-        for (int i=0; i<N_procs; i++)
-            procs[i] = i;
-        ThreadPool::setProcessAffinity(procs);
-    }
-    ThreadPool tpool(N_threads);
-    
-    // Create the MeshDataStruct
-    fillHalo<double> fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
-    std::vector<IO::MeshDataStruct> meshData(1);
-    meshData[0].meshName = "domain";
-    meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );
-    std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
-    std::shared_ptr<IO::Variable> PressVar( new IO::Variable() );
-    std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
-    std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
-    PhaseVar->name = "phase";
-    PhaseVar->type = IO::VariableType::VolumeVariable;
-    PhaseVar->dim = 1;
-    PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
-    meshData[0].vars.push_back(PhaseVar);
-    PressVar->name = "Pressure";
-    PressVar->type = IO::VariableType::VolumeVariable;
-    PressVar->dim = 1;
-    PressVar->data.resize(Nx-2,Ny-2,Nz-2);
-    meshData[0].vars.push_back(PressVar);
-    SignDistVar->name = "SignDist";
-    SignDistVar->type = IO::VariableType::VolumeVariable;
-    SignDistVar->dim = 1;
-    SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
-    meshData[0].vars.push_back(SignDistVar);
-    BlobIDVar->name = "BlobID";
-    BlobIDVar->type = IO::VariableType::VolumeVariable;
-    BlobIDVar->dim = 1;
-    BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
-    meshData[0].vars.push_back(BlobIDVar);
-
-	//************ MAIN ITERATION LOOP ***************************************/
-    PROFILE_START("Loop");
-    BlobIDstruct last_ids, last_index;
-    BlobIDList last_id_map;
-    writeIDMap(ID_map_struct(),0,id_map_filename);
-    AnalysisWaitIdStruct work_ids;
-	while (timestep < timestepMax && err > tol ) {
-        //if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
-        PROFILE_START("Update");
-
-		//*************************************************************************
-		// Fused Color Gradient and Collision 
-		//*************************************************************************
-		ScaLBL_D3Q19_ColorCollide_gen( ID,f_even,f_odd,Phi,ColorGrad,
-							 Velocity,Nx,Ny,Nz,tau1,tau2,rho1,rho2,alpha,beta,Fx,Fy,Fz);
-		//*************************************************************************
-
-		ScaLBL_DeviceBarrier();
-		//*************************************************************************
-		// Pack and send the D3Q19 distributions
-		ScaLBL_Comm.SendD3Q19(f_even, f_odd);
-		//*************************************************************************
-
-		//*************************************************************************
-		// 		Carry out the density streaming step for mass transport
-		//*************************************************************************
-		ScaLBL_D3Q7_ColorCollideMass(ID, A_even, A_odd, B_even, B_odd, Den, Phi,
-								ColorGrad, Velocity, beta, N, pBC);
-		//*************************************************************************
-
-		ScaLBL_DeviceBarrier();
-		MPI_Barrier(comm);
-		//*************************************************************************
-		// 		Swap the distributions for momentum transport
-		//*************************************************************************
-  		ScaLBL_D3Q19_Swap(ID, f_even, f_odd, Nx, Ny, Nz);
-		//*************************************************************************
-
-		ScaLBL_DeviceBarrier();
-		MPI_Barrier(comm);
-		//*************************************************************************
-		// Wait for communications to complete and unpack the distributions
-		ScaLBL_Comm.RecvD3Q19(f_even, f_odd);
-		//*************************************************************************
-
-		ScaLBL_DeviceBarrier();
-		//*************************************************************************
-		// Pack and send the D3Q7 distributions
-		ScaLBL_Comm.BiSendD3Q7(A_even, A_odd, B_even, B_odd);
-		//*************************************************************************
-
-		ScaLBL_DeviceBarrier();
-		ScaLBL_D3Q7_Swap(ID, A_even, A_odd, Nx, Ny, Nz);
-		ScaLBL_D3Q7_Swap(ID, B_even, B_odd, Nx, Ny, Nz);
-
-		ScaLBL_DeviceBarrier();
-		MPI_Barrier(comm);
-
-		//*************************************************************************
-		// Wait for communication and unpack the D3Q7 distributions
-		ScaLBL_Comm.BiRecvD3Q7(A_even, A_odd, B_even, B_odd);
-		//*************************************************************************
-
-		ScaLBL_DeviceBarrier();
-		//..................................................................................
-		ScaLBL_D3Q7_Density(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
-		ScaLBL_D3Q7_Density(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
-		//*************************************************************************
-		// 		Compute the phase indicator field 
-		//*************************************************************************
-		ScaLBL_DeviceBarrier();
-		MPI_Barrier(comm);
-
-		ScaLBL_ComputePhaseField(ID, Phi, Den, N);
-		//*************************************************************************
-		ScaLBL_Comm.SendHalo(Phi);
-		ScaLBL_DeviceBarrier();
-		ScaLBL_Comm.RecvHalo(Phi);
-		//*************************************************************************
-
-		ScaLBL_DeviceBarrier();
+		//...........................................................................
+		if (rank==0)	printf ("Create ScaLBL_Communicator \n");
+		// Create a communicator for the device (will use optimized layout)
+		ScaLBL_Communicator ScaLBL_Comm(Mask);
+		//Create a second communicator based on the regular data layout
+		ScaLBL_Communicator ScaLBL_Comm_Regular(Mask);
 		
-		// Pressure boundary conditions
-		if (BoundaryCondition==1 && Mask.kproc == 0)	{
-			ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-	        ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-			ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
-		}
-		if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
-			ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-	        ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-			ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
-		}
-
-		// Velocity boundary conditions
-		if (BoundaryCondition==2 && Mask.kproc == 0)	{
-			ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-	        ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-			ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
-		}
-		if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
-			ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-	        ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-			ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
-		}
-
-		if (BoundaryCondition==3){
-			// Increase the pressure difference
-			dp += slope;
-			din  = 1.0+0.5*dp;
-			dout = 1.0-0.5*dp;
-			// set the initial boundary conditions
-			if (Mask.kproc == 0)	{
-				ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-	            ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-				ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			    ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
-			}
-			if (Mask.kproc == nprocz-1){
-				ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-	            ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-				ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			    ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
-			}
-		}
-
-		// Set flux boundary condition
-		if (BoundaryCondition==4){
-			SumLocal=0.0;
-			if (pBC && Dm.kproc == 0){
-				SumLocal = ScaLBL_D3Q19_Flux_BC_z(ID,f_even,f_odd,flux,InletArea,Nx,Ny,Nz);
-			}
-			MPI_Allreduce(&SumLocal,&SumGlobal,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
-			din = flux / InletArea + SumGlobal;
-
-			if (pBC && Dm.kproc == 0){
-				ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-				ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-				ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-				ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
-			}
-
-			if (pBC && Dm.kproc == nprocz-1){
-				// if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout);
-				ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-				ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
-				ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-				ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
-			}
-		}
-
+		if (rank==0)	printf ("Set up memory efficient layout \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);
 		MPI_Barrier(comm);
-		PROFILE_STOP("Update");
 
-		// Timestep completed!
-		timestep++;
-  
-        // Run the analysis
-		  run_analysis(timestep,RESTART_INTERVAL,rank_info,*Averages,last_ids,last_index,last_id_map,
-	            Nx,Ny,Nz,pBC,beta,err,Phi,Pressure,Velocity,ID,f_even,f_odd,Den,
-	            LocalRestartFile,meshData,fillData,tpool,work_ids);
+		//...........................................................................
+		//				MAIN  VARIABLES ALLOCATED HERE
+		//...........................................................................
+		// LBM variables
+		if (rank==0)	printf ("Allocating distributions \n");
+		//......................device distributions.................................
+		int dist_mem_size = Np*sizeof(double);
+		if (rank==0)	printf ("Allocating distributions \n");
 
+		int *NeighborList;
+		int *dvcMap;
+		//		double *f_even,*f_odd;
+		double *fq, *Aq, *Bq;
+		double *Den, *Phi;
+		double *ColorGrad;
+		double *Vel;
+		double *Pressure;
+		
+		//...........................................................................
+		ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
+		ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np);
 
-        // Save the timers
-        if ( timestep%50==0 )
-            PROFILE_SAVE("lbpm_color_simulator",1);
-	}
-    tpool.wait_pool_finished();
-    PROFILE_STOP("Loop");
-	//************************************************************************
-	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;
+		ScaLBL_AllocateDeviceMemory((void **) &fq, 19*dist_mem_size);
+		ScaLBL_AllocateDeviceMemory((void **) &Aq, 7*dist_mem_size);
+		ScaLBL_AllocateDeviceMemory((void **) &Bq, 7*dist_mem_size);
+		ScaLBL_AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
+		ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Nx*Ny*Nz);		
+		ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np);
+		ScaLBL_AllocateDeviceMemory((void **) &Vel, 3*sizeof(double)*Np);
+		ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np);
+		
+		//...........................................................................
+		// Update GPU data structures
+		if (rank==0)	printf ("Setting up device map and neighbor list \n");
+		int *TmpMap;
+		TmpMap=new int[Np];
+		for (k=1; k<Nz-1; k++){
+			for (j=1; j<Ny-1; j++){
+				for (i=1; i<Nx-1; i++){
+					int idx=Map(i,j,k);
+					if (!(idx < 0))
+						TmpMap[idx] = k*Nx*Ny+j*Nx+i;
+				}
+			}
+		}
+		//for (int idx=0; idx<Np; idx++) printf("Map=%i\n",TmpMap[idx]);
+
+		ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
+		ScaLBL_DeviceBarrier();
+		delete [] TmpMap;
+		
+		// copy the neighbor list 
+		ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
+		// initialize phi based on PhaseLabel (include solid component labels)
+		ScaLBL_CopyToDevice(Phi, PhaseLabel, N*sizeof(double));
+		//...........................................................................
+
+		if (rank==0)	printf ("Initializing distributions \n");
+		// Initialize the phase field and variables
+		ScaLBL_D3Q19_Init(fq, Np);
+		if (rank==0)	printf ("Initializing phase field \n");
+		ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, Np);
+		if (Dm.kproc==0){
+			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
+			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
+			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
+		}
+		if (Dm.kproc == nprocz-1){
+			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
+			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
+			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
+		}
+		
+
+		if (Restart == true){
+			if (rank==0){
+				printf("Reading restart file! \n");
+				ifstream restart("Restart.txt");
+				if (restart.is_open()){
+					restart  >> timestep;
+					printf("Restarting from timestep =%i \n",timestep);
+				}
+				else{
+					printf("WARNING:No Restart.txt file, setting timestep=0 \n");
+					timestep=0;
+				}
+			}
+			MPI_Bcast(&timestep,1,MPI_INT,0,comm);
+
+			// Read in the restart file to CPU buffers
+			double *cDen = new double[2*N];
+			double *cDistEven = new double[10*N];
+			double *cDistOdd = new double[9*N];
+			ReadCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
+			// Copy the restart data to the GPU
+			ScaLBL_CopyToDevice(f_even,cDistEven,10*N*sizeof(double));
+			ScaLBL_CopyToDevice(f_odd,cDistOdd,9*N*sizeof(double));
+			ScaLBL_CopyToDevice(Den,cDen,2*N*sizeof(double));
+			ScaLBL_DeviceBarrier();
+			delete [] cDen;
+			delete [] cDistEven;
+			delete [] cDistOdd;
+			MPI_Barrier(comm);
+		}
+
+		//.......................................................................
+		// Once phase has been initialized, map solid to account for 'smeared' interface
+		//for (i=0; i<N; i++)	Averages.SDs(i) -= (1.0);
+		// Make sure the id match for the two domains
+		for (i=0; i<N; i++)	Dm.id[i] = Mask.id[i];
+		//.......................................................................
+		// Finalize setup for averaging domain
+		Averages->UpdateSolid();
+		//.......................................................................		
+		ScaLBL_D3Q19_Pressure(fq,Pressure,Np);
+    	ScaLBL_D3Q19_Momentum(fq,Vel,Np);
+		//...........................................................................
+		// Copy the phase indicator field for the earlier timestep
+		ScaLBL_DeviceBarrier();
+		ScaLBL_CopyToHost(Averages->Phase_tplus.data(),Phi,N*sizeof(double));
+		//...........................................................................
+		// Copy the data for for the analysis timestep
+		//...........................................................................
+		// Copy the phase from the GPU -> CPU
+		//...........................................................................
+		ScaLBL_DeviceBarrier();
+		ScaLBL_CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
+		ScaLBL_CopyToHost(Averages->Press.data(),Pressure,Np*sizeof(double));
+		ScaLBL_CopyToHost(Averages->Vel_x.data(),&Velocity[0],Np*sizeof(double));
+		ScaLBL_CopyToHost(Averages->Vel_y.data(),&Velocity[N],Np*sizeof(double));
+		ScaLBL_CopyToHost(Averages->Vel_z.data(),&Velocity[2*N],Np*sizeof(double));
 	
-	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");
-	
-	// ************************************************************************
+		double TmpDat;
+		TmpDat = new double [Np];
+    	ScaLBL_CopyToHost(&TmpDat[0],&Pressure[0],SIZE);
+    	ScaLBL_Comm.RegularLayout(Map,TmpDat,Averages->Press.data());
+		//...........................................................................
 
-    PROFILE_STOP("Main");
-    PROFILE_SAVE("lbpm_color_simulator",1);
-	// ****************************************************
-	MPI_Barrier(comm);
-  } // Limit scope so variables that contain communicators will free before MPI_Finialize
-  MPI_Comm_free(&comm);
-  MPI_Finalize();
+		if (rank==0) printf("********************************************************\n");
+		if (rank==0)	printf("No. of timesteps: %i \n", timestepMax);
+
+		//.......create and start timer............
+		double starttime,stoptime,cputime;
+		ScaLBL_DeviceBarrier();
+		MPI_Barrier(comm);
+		starttime = MPI_Wtime();
+		//.........................................
+
+		err = 1.0; 	
+		double sat_w_previous = 1.01; // slightly impossible value!
+		if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol);
+		// Create the thread pool
+		int N_threads = 4;
+		if ( provided_thread_support < MPI_THREAD_MULTIPLE )
+			N_threads = 0;
+		if ( N_threads > 0 ) {
+			// Set the affinity
+			int N_procs = ThreadPool::getNumberOfProcessors();
+			std::vector<int> procs(N_procs);
+			for (int i=0; i<N_procs; i++)
+				procs[i] = i;
+			ThreadPool::setProcessAffinity(procs);
+		}
+		ThreadPool tpool(N_threads);
+
+		// Create the MeshDataStruct
+		fillHalo<double> fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
+		std::vector<IO::MeshDataStruct> meshData(1);
+		meshData[0].meshName = "domain";
+		meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );
+		std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
+		std::shared_ptr<IO::Variable> PressVar( new IO::Variable() );
+		std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
+		std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
+		PhaseVar->name = "phase";
+		PhaseVar->type = IO::VariableType::VolumeVariable;
+		PhaseVar->dim = 1;
+		PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
+		meshData[0].vars.push_back(PhaseVar);
+		PressVar->name = "Pressure";
+		PressVar->type = IO::VariableType::VolumeVariable;
+		PressVar->dim = 1;
+		PressVar->data.resize(Nx-2,Ny-2,Nz-2);
+		meshData[0].vars.push_back(PressVar);
+		SignDistVar->name = "SignDist";
+		SignDistVar->type = IO::VariableType::VolumeVariable;
+		SignDistVar->dim = 1;
+		SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
+		meshData[0].vars.push_back(SignDistVar);
+		BlobIDVar->name = "BlobID";
+		BlobIDVar->type = IO::VariableType::VolumeVariable;
+		BlobIDVar->dim = 1;
+		BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
+		meshData[0].vars.push_back(BlobIDVar);
+
+		//************ MAIN ITERATION LOOP ***************************************/
+		PROFILE_START("Loop");
+		BlobIDstruct last_ids, last_index;
+		BlobIDList last_id_map;
+		writeIDMap(ID_map_struct(),0,id_map_filename);
+		AnalysisWaitIdStruct work_ids;
+		while (timestep < timestepMax && err > tol ) {
+			//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
+			PROFILE_START("Update");
+			// *************ODD TIMESTEP*************
+			timestep++;
+			// Compute the Phase indicator field
+			// Read for Aq, Bq happens in this routine (requires communication)
+			ScaLBL_Comm.BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
+			ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm.next, Np, Np);
+			ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
+			ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm.next, Np);
+
+			// Halo exchange for phase field
+			ScaLBL_Comm_Regular.SendHalo(Phi);
+			ScaLBL_Comm_Regular.RecvHalo(Phi);
+			
+			// Perform the collision operation
+			ScaLBL_Comm.SendD3Q19AA(fq); //READ FROM NORMAL
+			ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB,
+					alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm.next, Np, Np);
+			ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
+			// Set BCs
+			if (typeBC > 0){
+				ScaLBL_Comm.Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
+				ScaLBL_Comm.Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
+			}
+			if (typeBC == 3){
+				ScaLBL_Comm.D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
+				ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+			}
+			if (typeBC == 4){
+				din = ScaLBL_Comm.D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
+				ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+			}
+			ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB,
+					alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm.next, Np);
+			ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+
+			// *************EVEN TIMESTEP*************
+			timestep++;
+			// Compute the Phase indicator field
+			ScaLBL_Comm.BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
+			ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm.next, Np, Np);
+			ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
+			ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm.next, Np);
+
+			// Halo exchange for phase field
+			ScaLBL_Comm_Regular.SendHalo(Phi);
+			ScaLBL_Comm_Regular.RecvHalo(Phi);
+
+			// Perform the collision operation
+			ScaLBL_Comm.SendD3Q19AA(fq); //READ FORM NORMAL
+			ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB,
+					alpha, beta, Fx, Fy, Fz,  Nx, Nx*Ny, ScaLBL_Comm.next, Np, Np);
+			ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
+			// Set boundary conditions
+			if (typeBC > 0){
+				ScaLBL_Comm.Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
+				ScaLBL_Comm.Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
+			}
+			if (typeBC == 3){
+				ScaLBL_Comm.D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
+				ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+			}
+			else if (typeBC == 4){
+				din = ScaLBL_Comm.D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
+				ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+			}
+			ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB,
+					alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm.next, Np);
+			ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+			//************************************************************************
+			
+			MPI_Barrier(comm);
+			PROFILE_STOP("Update");
+
+			// Run the analysis
+			run_analysis(timestep,RESTART_INTERVAL,rank_info,*Averages,last_ids,last_index,last_id_map,
+					Nx,Ny,Nz,pBC,beta,err,Phi,Pressure,Velocity,ID,f_even,f_odd,Den,
+					LocalRestartFile,meshData,fillData,tpool,work_ids);
+
+			// Save the timers
+			if ( timestep%50==0 )
+				PROFILE_SAVE("lbpm_color_simulator",1);
+		}
+		tpool.wait_pool_finished();
+		PROFILE_STOP("Loop");
+		//************************************************************************
+		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");
+
+		// ************************************************************************
+
+		PROFILE_STOP("Main");
+		PROFILE_SAVE("lbpm_color_simulator",1);
+		// ****************************************************
+		MPI_Barrier(comm);
+	} // Limit scope so variables that contain communicators will free before MPI_Finialize
+	MPI_Comm_free(&comm);
+	MPI_Finalize();
 }