From 1eb4d9eeab3ea6f1cf19062fd81e47557480d7d9 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Sat, 3 Oct 2015 19:51:35 -0400
Subject: [PATCH 01/25] Added dynamic (time dependent) pressure boundary
 condition (3)

---
 tests/lbpm_color_simulator.cpp | 41 +++++++++++++++++++++++++++++++++-
 tests/lbpm_color_simulator.h   |  6 ++---
 2 files changed, 43 insertions(+), 4 deletions(-)

diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index 7685aaff..e4f61cff 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -290,6 +290,7 @@ int main(int argc, char **argv)
 		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 (InitialCondition==0) printf("Initial conditions assigned from phase ID file \n");
 		if (InitialCondition==1) printf("Initial conditions assigned from restart file \n");
 		printf("********************************************************\n");
@@ -297,10 +298,12 @@ int main(int argc, char **argv)
 
 	// Initialized domain and averaging framework for Two-Phase Flow
 	bool pBC,velBC;
-	if (BoundaryCondition==1)	pBC=true;
+	if (BoundaryCondition==1 || BoundaryCondition==3)
+								pBC=true;
 	else						pBC=false;
 	if (BoundaryCondition==2)	velBC=true;
 	else						velBC=false;
+
 	bool Restart;
 	if (InitialCondition==1)    Restart=true;
 	else 						Restart=false;
@@ -624,6 +627,25 @@ int main(int argc, char **argv)
 		SetPhiSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 	}
 
+	// Set dynamic pressure boundary conditions
+	double dp, slope;
+	if (BoundaryCondition==3){
+		dp = din;
+		// set the initial value
+		din  = 1.0+0.5*dp;
+		dout = 1.0-0.5*dp;
+		slope = (dout-din)/timestepMax;
+		// set the initial boundary conditions
+		if (Dm.kproc == 0)	{
+			PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
+			ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+		}
+		if (Dm.kproc == nprocz-1){
+			PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
+			ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+		}
+	}
+
 	ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
 	ComputeVelocityD3Q19(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
 
@@ -787,6 +809,23 @@ int main(int argc, char **argv)
 			//ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 			SetPhiSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 		}
+
+		if (BoundaryCondition==3){
+			// Increase the pressure difference
+			dp += timestep*slope;
+			din  = 1.0+0.5*dp;
+			dout = 1.0-0.5*dp;
+			// set the initial boundary conditions
+			if (Dm.kproc == 0)	{
+				PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
+				ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			}
+			if (Dm.kproc == nprocz-1){
+				PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
+				ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			}
+		}
+
 		//...................................................................................
 
 		MPI_Barrier(comm);
diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index 45cd1507..ede8e377 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -2,7 +2,7 @@
 #include "common/Array.h"
 
 #define ANALYSIS_INTERVAL 1000
-#define BLOBID_INTERVAL 1000
+#define BLOBID_INTERVAL 250
 
 enum AnalysisType{ AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02, 
     CopyAverages=0x04, CalcDist=0x08, CreateRestart=0x10 };
@@ -191,7 +191,7 @@ void run_analysis( int timestep, int restart_interval,
         // Identify blobs and update global ids in time
         type = static_cast<AnalysisType>( type | IdentifyBlobs );
     }
-/*    #ifdef USE_CUDA
+    #ifdef USE_CUDA
         if ( tpool.getQueueSize()<=3 && tpool.getNumThreads()>0 && timestep%50==0 ) {
             // Keep a few blob identifications queued up to keep the processors busy,
             // allowing us to track the blobs as fast as possible
@@ -199,7 +199,7 @@ void run_analysis( int timestep, int restart_interval,
             type = static_cast<AnalysisType>( type | IdentifyBlobs );
         }
     #endif
- */
+
     if ( timestep%ANALYSIS_INTERVAL == 0 ) {
         // Copy the averages to the CPU (and identify blobs)
         type = static_cast<AnalysisType>( type | CopyAverages );

From a13b8a313d4042ea3e3b08bf77894721ad6398da Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Sun, 4 Oct 2015 17:08:22 -0400
Subject: [PATCH 02/25] Dynamic pressure BC implemented

---
 tests/lbpm_color_simulator.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index e4f61cff..d6f3e213 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -630,11 +630,12 @@ int main(int argc, char **argv)
 	// Set dynamic pressure boundary conditions
 	double dp, slope;
 	if (BoundaryCondition==3){
+		slope = (dout-din)/timestepMax;
 		dp = din;
+		if (rank==0) printf("Change in pressure / time =%f \n",slope);
 		// set the initial value
 		din  = 1.0+0.5*dp;
 		dout = 1.0-0.5*dp;
-		slope = (dout-din)/timestepMax;
 		// set the initial boundary conditions
 		if (Dm.kproc == 0)	{
 			PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
@@ -812,7 +813,7 @@ int main(int argc, char **argv)
 
 		if (BoundaryCondition==3){
 			// Increase the pressure difference
-			dp += timestep*slope;
+			dp += slope;
 			din  = 1.0+0.5*dp;
 			dout = 1.0-0.5*dp;
 			// set the initial boundary conditions
@@ -833,7 +834,7 @@ int main(int argc, char **argv)
 
 		// Timestep completed!
 		timestep++;
-
+  
         // Run the analysis, blob identification, and write restart files
         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,

From d9996ed1b0c8ee26e027482a166a8dd6269928b9 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Tue, 6 Oct 2015 16:55:18 -0400
Subject: [PATCH 03/25] Performance testing signed distance conversion on CPU

---
 common/TwoPhase.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/common/TwoPhase.cpp b/common/TwoPhase.cpp
index df6abae7..0e4e8fbe 100644
--- a/common/TwoPhase.cpp
+++ b/common/TwoPhase.cpp
@@ -175,7 +175,7 @@ TwoPhase::~TwoPhase()
 
 void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, DoubleArray &DistData)
 {
-/*	double factor,temp,value;
+	double factor,temp,value;
 	factor=0.5/Beta;
 	for (int n=0; n<Nx*Ny*Nz; n++){
 	 	value = ColorData[n];
@@ -199,14 +199,15 @@ void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, Double
 	}
 
 	SSO(DistData,Dm.id,Dm,10);
-*/
-	for (int k=0; k<Nz; k++){
+
+	/*	for (int k=0; k<Nz; k++){
 		for (int j=0; j<Ny; j++){
 			for (int i=0; i<Nx; i++){
 				DistData(i,j,k) = ColorData(i,j,k);
 			}
 		}
 	}
+	*/
 //	for (int n=0; n<Nx*Ny*Nz; n++)	DistData[n] = ColorData[n];
 }
 

From 2a1d97239e399a3c4a87550f7ba377d408e82bf4 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Wed, 4 Nov 2015 20:26:17 -0500
Subject: [PATCH 04/25] Added signed distance

---
 common/TwoPhase.cpp            | 28 ++++++++++++++++------------
 tests/lbpm_color_simulator.cpp |  2 +-
 2 files changed, 17 insertions(+), 13 deletions(-)

diff --git a/common/TwoPhase.cpp b/common/TwoPhase.cpp
index 0e4e8fbe..32397457 100644
--- a/common/TwoPhase.cpp
+++ b/common/TwoPhase.cpp
@@ -177,29 +177,31 @@ void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, Double
 {
 	double factor,temp,value;
 	factor=0.5/Beta;
-	for (int n=0; n<Nx*Ny*Nz; n++){
-	 	value = ColorData[n];
-		if (value > 0.999 ) DistData[n] = 4.0;
-		else if (value < -0.999 ) DistData[n] = -4.0;
-		else 	DistData[n] = factor*log((1.0+value)/(1.0-value));
-		if (DistData[n] > 1.0)  DistData[n] = 1.0;
-		if (DistData[n] < -1.0) DistData[n] = -1.0;
-	}
 	// Initialize to -1,1 (segmentation)
 	for (int k=0; k<Nz; k++){
 		for (int j=0; j<Ny; j++){
 			for (int i=0; i<Nx; i++){
 				value = ColorData(i,j,k);
-				temp = factor*log((1.0+value)/(1.0-value));
-				if (temp > 1.0) DistData(i,j,k) = 1.0;
-				else if (temp < -1.0) DistData(i,j,k) = -1.0;
+			 	temp = factor*log((1.0+value)/(1.0-value));
+				if (value > 0.8) DistData(i,j,k) = 2.94*factor;
+				else if (value < -0.8) DistData(i,j,k) = -2.94*factor;
 				else DistData(i,j,k) = temp;
+				// Basic threshold
+				//if (value > 0) DistData(i,j,k) = 1.0;
+				//else DistData(i,j,k) = -1.0;
 			}
 		}
 	}
 
-	SSO(DistData,Dm.id,Dm,10);
+	SSO(DistData,Dm.id,Dm,40);
 
+        for (int k=0; k<Nz; k++){
+	  for (int j=0; j<Ny; j++){
+	    for (int i=0; i<Nx; i++){
+	      DistData(i,j,k) += 1.0;
+	    }
+	  }
+	}	
 	/*	for (int k=0; k<Nz; k++){
 		for (int j=0; j<Ny; j++){
 			for (int i=0; i<Nx; i++){
@@ -304,6 +306,8 @@ void TwoPhase::UpdateMeshValues()
 {
 	int i,j,k,n;
 	//...........................................................................
+	Dm.CommunicateMeshHalo(SDn);
+	//...........................................................................
 	// Compute the gradients of the phase indicator and signed distance fields
 	pmmc_MeshGradient(SDn,SDn_x,SDn_y,SDn_z,Nx,Ny,Nz);
 	//...........................................................................
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index d6f3e213..7c829e91 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -578,7 +578,7 @@ int main(int argc, char **argv)
 	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); //
+	// for (i=0; i<N; i++)	Averages->SDs(i) -= (1.0); //
 	//.......................................................................
 	// Finalize setup for averaging domain
 	//Averages->SetupCubes(Dm);

From fe5caee295fa32d7bcae583b1c40b3064c6b3686 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 09:24:30 -0500
Subject: [PATCH 05/25] Added visualization to threadpool work item Restart

---
 tests/lbpm_color_simulator.cpp | 34 -----------------------------
 tests/lbpm_color_simulator.h   | 39 ++++++++++++++++++++++++++++++++++
 2 files changed, 39 insertions(+), 34 deletions(-)

diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index 7c829e91..2a084655 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -873,40 +873,6 @@ int main(int argc, char **argv)
 	DeviceBarrier();
 	CopyToHost(Averages->Phase.get(),Phi,N*sizeof(double));
 */
-    // 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::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::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::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::VolumeVariable;
-    BlobIDVar->dim = 1;
-    BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
-    meshData[0].vars.push_back(BlobIDVar);
-    
-    fillData.copy(Averages->SDn,PhaseVar->data);
-    fillData.copy(Averages->SDs,SignDistVar->data);
-    fillData.copy(Averages->Label_NWP,BlobIDVar->data);
-    IO::writeData( 0, meshData, 2, comm );
     
 /*	Averages->WriteSurfaces(0);
 
diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index ede8e377..3f31145b 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -28,6 +28,45 @@ public:
         PROFILE_START("Save Checkpoint",1);
         WriteCheckpoint(filename,cDen.get(),cDistEven.get(),cDistOdd.get(),N);
         PROFILE_STOP("Save Checkpoint",1);
+
+        // Write  VisIT files
+        PROFILE_START("Save Vis",1);
+        // 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::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::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::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::VolumeVariable;
+        BlobIDVar->dim = 1;
+        BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
+        meshData[0].vars.push_back(BlobIDVar);
+
+        fillData.copy(Averages->SDn,PhaseVar->data);
+        fillData.copy(Averages->SDs,SignDistVar->data);
+        fillData.copy(Averages->Label_NWP,BlobIDVar->data);
+        IO::writeData( 0, meshData, 2, comm );
+        PROFILE_STOP("Save Vis",1);
+
         PROFILE_SAVE("lbpm_color_simulator",1);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     };

From 4431239de8e2332016ee20bde637965c5d5fffcb Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 10:13:05 -0500
Subject: [PATCH 06/25] Working on adding vis to workpool (doesn't compile)

---
 tests/lbpm_color_simulator.h | 95 +++++++++++++++++++++---------------
 1 file changed, 57 insertions(+), 38 deletions(-)

diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index 3f31145b..96a7dfad 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -29,44 +29,6 @@ public:
         WriteCheckpoint(filename,cDen.get(),cDistEven.get(),cDistOdd.get(),N);
         PROFILE_STOP("Save Checkpoint",1);
 
-        // Write  VisIT files
-        PROFILE_START("Save Vis",1);
-        // 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::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::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::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::VolumeVariable;
-        BlobIDVar->dim = 1;
-        BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
-        meshData[0].vars.push_back(BlobIDVar);
-
-        fillData.copy(Averages->SDn,PhaseVar->data);
-        fillData.copy(Averages->SDs,SignDistVar->data);
-        fillData.copy(Averages->Label_NWP,BlobIDVar->data);
-        IO::writeData( 0, meshData, 2, comm );
-        PROFILE_STOP("Save Vis",1);
-
         PROFILE_SAVE("lbpm_color_simulator",1);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     };
@@ -161,6 +123,63 @@ private:
     BlobIDList new_list;
 };
 
+class WriteVisWorkItem: public ThreadPool::WorkItem
+{
+public:
+	WriteVisWorkItem(AnalysisType type_, Domain &Dm_, TwoPhase& Averages_):
+	type(type_), Domain(Dm_), Averages(Averages_) { }
+
+    virtual void run() {
+        ThreadPool::WorkItem::d_state = 1;  // Change state to in progress
+
+        // Write  VisIT files
+        PROFILE_START("Save Vis",1);
+        // 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::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::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::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::VolumeVariable;
+        BlobIDVar->dim = 1;
+        BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
+        meshData[0].vars.push_back(BlobIDVar);
+
+        fillData.copy(Averages->SDn,PhaseVar->data);
+        fillData.copy(Averages->SDs,SignDistVar->data);
+        fillData.copy(Averages->Label_NWP,BlobIDVar->data);
+        IO::writeData( 0, meshData, 2, comm );
+        PROFILE_STOP("Save Vis",1);
+
+        PROFILE_SAVE("lbpm_color_simulator",1);
+        ThreadPool::WorkItem::d_state = 2;  // Change state to finished
+    };
+private:
+    WriteRestartWorkItem();
+    AnalysisType type;
+    Domain& Dm;
+    TwoPhase& Averages;
+};
+
 
 // Helper class to run the analysis from within a thread
 // Note: Averages will be modified after the constructor is called

From 192ab6dd53ca7afba1d1a6456e833cba1bb932df Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 10:18:59 -0500
Subject: [PATCH 07/25] Working on adding vis to workpool (doesn't compile)

---
 tests/lbpm_color_simulator.h | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index 96a7dfad..29f2aa2a 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -126,8 +126,8 @@ private:
 class WriteVisWorkItem: public ThreadPool::WorkItem
 {
 public:
-	WriteVisWorkItem(AnalysisType type_, Domain &Dm_, TwoPhase& Averages_):
-	type(type_), Domain(Dm_), Averages(Averages_) { }
+	WriteVisWorkItem(AnalysisType type_, TwoPhase& Averages_):
+	type(type_), Averages(Averages_) { }
 
     virtual void run() {
         ThreadPool::WorkItem::d_state = 1;  // Change state to in progress
@@ -135,10 +135,13 @@ public:
         // Write  VisIT files
         PROFILE_START("Save Vis",1);
         // Create the MeshDataStruct
-        fillHalo<double> fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
+        Nx=Averages.Dm.Nx;
+        Ny=Averages.Dm.Ny;
+        Nz=Averages.Dm.Nz;
+        fillHalo<double> fillData(Averages.Dm.Comm,Averages.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) );
+        meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Averages.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() );
@@ -176,8 +179,8 @@ public:
 private:
     WriteRestartWorkItem();
     AnalysisType type;
-    Domain& Dm;
     TwoPhase& Averages;
+    int Nx,Ny,Nz;
 };
 
 

From 27e4e81c997d2cc9f539aed9d5e819dd69e95f8b Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 10:21:41 -0500
Subject: [PATCH 08/25] Working on adding vis to workpool (doesn't compile)

---
 tests/lbpm_color_simulator.h | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index 29f2aa2a..957d9d49 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -138,6 +138,9 @@ public:
         Nx=Averages.Dm.Nx;
         Ny=Averages.Dm.Ny;
         Nz=Averages.Dm.Nz;
+        Lx=Averages.Dm.Lx;
+        Ly=Averages.Dm.Ly;
+        Lz=Averages.Dm.Lz;
         fillHalo<double> fillData(Averages.Dm.Comm,Averages.Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
         std::vector<IO::MeshDataStruct> meshData(1);
         meshData[0].meshName = "domain";
@@ -177,10 +180,11 @@ public:
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     };
 private:
-    WriteRestartWorkItem();
+    WriteVisWorkItem();
     AnalysisType type;
     TwoPhase& Averages;
     int Nx,Ny,Nz;
+    double Lx,Ly,Lz;
 };
 
 

From 7e529a905787918249faf92f837883f741343d4e Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 10:26:55 -0500
Subject: [PATCH 09/25] Working on adding vis to workpool (doesn't compile)

---
 tests/lbpm_color_simulator.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index 957d9d49..d052d996 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -141,7 +141,7 @@ public:
         Lx=Averages.Dm.Lx;
         Ly=Averages.Dm.Ly;
         Lz=Averages.Dm.Lz;
-        fillHalo<double> fillData(Averages.Dm.Comm,Averages.Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
+//        fillHalo<double> fillData(Averages.Dm.Comm,Averages.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(Averages.Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );

From 64b94dd929a86e9a40aa44125672df635642cd16 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 10:28:58 -0500
Subject: [PATCH 10/25] Working on adding vis to workpool (doesn't compile)

---
 tests/lbpm_color_simulator.h | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index d052d996..d33e6a1b 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -131,7 +131,8 @@ public:
 
     virtual void run() {
         ThreadPool::WorkItem::d_state = 1;  // Change state to in progress
-
+        MPI_Comm newcomm;
+        MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
         // Write  VisIT files
         PROFILE_START("Save Vis",1);
         // Create the MeshDataStruct
@@ -141,7 +142,7 @@ public:
         Lx=Averages.Dm.Lx;
         Ly=Averages.Dm.Ly;
         Lz=Averages.Dm.Lz;
-//        fillHalo<double> fillData(Averages.Dm.Comm,Averages.Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
+        fillHalo<double> fillData(newcomm,Averages.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(Averages.Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );

From d5592c14e6fd30a8a934c1fde7cd5a2972b34717 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 10:32:54 -0500
Subject: [PATCH 11/25] Working on adding vis to workpool (doesn't compile)

---
 tests/lbpm_color_simulator.cpp | 1 -
 tests/lbpm_color_simulator.h   | 5 +++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index 2a084655..645ae670 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -7,7 +7,6 @@
 #include <fstream>
 
 #include "common/ScaLBL.h"
-#include "common/Communication.h"
 #include "common/TwoPhase.h"
 #include "common/MPI_Helpers.h"
 
diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index d33e6a1b..e19160af 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -1,5 +1,6 @@
 // Run the analysis, blob identification, and write restart files
 #include "common/Array.h"
+#include "common/Communication.h"
 
 #define ANALYSIS_INTERVAL 1000
 #define BLOBID_INTERVAL 250
@@ -15,7 +16,6 @@ struct AnalysisWaitIdStruct {
     ThreadPool::thread_id_t restart;
 };
 
-
 // Helper class to write the restart file from a seperate thread
 class WriteRestartWorkItem: public ThreadPool::WorkItem
 {
@@ -39,7 +39,6 @@ private:
     const int N;
 };
 
-
 // Helper class to compute the blob ids
 static const std::string id_map_filename = "lbpm_id_map.txt";
 typedef std::shared_ptr<std::pair<int,IntArray> > BlobIDstruct;
@@ -175,6 +174,8 @@ public:
         fillData.copy(Averages->SDs,SignDistVar->data);
         fillData.copy(Averages->Label_NWP,BlobIDVar->data);
         IO::writeData( 0, meshData, 2, comm );
+
+        MPI_Comm_free(&newcomm);
         PROFILE_STOP("Save Vis",1);
 
         PROFILE_SAVE("lbpm_color_simulator",1);

From c2b0cdc012e4ba665385a3599b9b957ef0c69836 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 18:33:52 -0500
Subject: [PATCH 12/25] Debugging vis writer

---
 tests/lbpm_color_simulator.h | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index e19160af..73b36b1b 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -170,10 +170,10 @@ public:
         BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
         meshData[0].vars.push_back(BlobIDVar);
 
-        fillData.copy(Averages->SDn,PhaseVar->data);
-        fillData.copy(Averages->SDs,SignDistVar->data);
-        fillData.copy(Averages->Label_NWP,BlobIDVar->data);
-        IO::writeData( 0, meshData, 2, comm );
+        fillData.copy(Averages.SDn,PhaseVar->data);
+        fillData.copy(Averages.SDs,SignDistVar->data);
+        fillData.copy(Averages.Label_NWP,BlobIDVar->data);
+        IO::writeData( 0, meshData, 2, newcomm );
 
         MPI_Comm_free(&newcomm);
         PROFILE_STOP("Save Vis",1);
@@ -187,6 +187,10 @@ private:
     TwoPhase& Averages;
     int Nx,Ny,Nz;
     double 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() );
 };
 
 

From b9ba20246a01dbe2e18fba716d745e25a9502d9a Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 18:35:31 -0500
Subject: [PATCH 13/25] Debugging vis writer

---
 tests/lbpm_color_simulator.h | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index 73b36b1b..29e58c03 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -187,10 +187,6 @@ private:
     TwoPhase& Averages;
     int Nx,Ny,Nz;
     double 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() );
 };
 
 

From 625c63bbb253c2585867566f5479884286db545e Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 19:11:41 -0500
Subject: [PATCH 14/25] Retaining the interface speed locally from
 pmmc_InterfaceSpeed

---
 common/pmmc.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/common/pmmc.h b/common/pmmc.h
index 1b0720c3..a413820b 100644
--- a/common/pmmc.h
+++ b/common/pmmc.h
@@ -4359,6 +4359,7 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
 			AvgVel(0) += temp*zeta*x;
 			AvgVel(1) += temp*zeta*y;
 			AvgVel(2) += temp*zeta*z;
+			SurfaceValues(r) = zeta;
 		}
 	}
 	//.............................................................................

From 88954a86ba33862e810f0dda27a414368cd50416 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 19:51:10 -0500
Subject: [PATCH 15/25] Fixing possible bug in interface speed

---
 common/pmmc.h | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/common/pmmc.h b/common/pmmc.h
index a413820b..bb593baa 100644
--- a/common/pmmc.h
+++ b/common/pmmc.h
@@ -4317,12 +4317,12 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
 //--------------------------------------------------------------------------------------------------------
 inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray &P_y, DoubleArray &P_z,
 									DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles,
-									  DoubleArray &SurfaceVector, DoubleArray &SurfaceValues, DoubleArray &AvgVel, 
+									  DoubleArray &SurfaceVector, DoubleArray &AvgSpeed, DoubleArray &AvgVel,
 									  int i, int j, int k, int npts, int ntris)
 {
 	Point A,B,C,P;
 	double x,y,z;
-	double s,s1,s2,s3,temp;
+	double s,s1,s2,s3,area;
 	double norm, zeta;
 
 	TriLinPoly Px,Py,Pz,Pt;
@@ -4342,7 +4342,7 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
 		s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
 		s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
 		s = 0.5*(s1+s2+s3);
-		temp = s*(s-s1)*(s-s2)*(s-s3);
+		area = sqrt(s*(s-s1)*(s-s2)*(s-s3));
 		// Compute the centroid P
 		P.x = 0.33333333333333333*(A.x+B.x+C.x);
 		P.y = 0.33333333333333333*(A.y+B.y+C.y);
@@ -4353,13 +4353,15 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
 			z = Pz.eval(P);
 			norm = sqrt(x*x+y*y+z*z);
 			if (norm==0.0) norm=1.0;
+			// Compute the interface speed from time derivative and gradient (Level Set Equation)
 			zeta = -Pt.eval(P) / norm;
-			temp = sqrt(temp)/norm;
+			//temp = sqrt(temp)/norm; <--- what was I thinking with this? (James)
 
-			AvgVel(0) += temp*zeta*x;
-			AvgVel(1) += temp*zeta*y;
-			AvgVel(2) += temp*zeta*z;
-			SurfaceValues(r) = zeta;
+			// Compute the average
+			AvgVel(0) += area*zeta*x;
+			AvgVel(1) += area*zeta*y;
+			AvgVel(2) += area*zeta*z;
+			AvgSpeed(r) = zeta*area;
 		}
 	}
 	//.............................................................................

From b788ba4dc48cb443585f55dd5eb756330e35a6ee Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 20:00:49 -0500
Subject: [PATCH 16/25] Added new quantity wwndnw to averaging for two-phase

---
 common/TwoPhase.cpp | 14 ++++++++++----
 common/TwoPhase.h   |  1 +
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/common/TwoPhase.cpp b/common/TwoPhase.cpp
index 32397457..e356d6aa 100644
--- a/common/TwoPhase.cpp
+++ b/common/TwoPhase.cpp
@@ -67,7 +67,7 @@ TwoPhase::TwoPhase(Domain &dm):
     Jwn(0), Jwn_global(0), Kwn(0), Kwn_global(0), KNwns(0), KNwns_global(0),
     KGwns(0), KGwns_global(0), trawn(0), trawn_global(0), trJwn(0), trJwn_global(0),
     trRwn(0), trRwn_global(0), nwp_volume_global(0), wp_volume_global(0),
-    As_global(0), dEs(0), dAwn(0), dAns(0)
+    As_global(0), dEs(0), dAwn(0), dAns(0), wwndnw(0), wwndnw_global(0)
 {
 	Nx=dm.Nx; Ny=dm.Ny; Nz=dm.Nz;
 	Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz*1.0;
@@ -147,8 +147,7 @@ TwoPhase::TwoPhase(Domain &dm):
 			fprintf(TIMELOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz ");				// Orientation tensors
 			fprintf(TIMELOG,"Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz ");
 			fprintf(TIMELOG,"Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz ");
-			fprintf(TIMELOG,"trawn trJwn trRwn Euler Kn Jn An\n");					// trimmed curvature & minkowski measures
-			//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
+			fprintf(TIMELOG,"trawn trJwn trRwn wwndnw Euler Kn Jn An\n");					// trimmed curvature & minkowski measures
 		}
 
 		NWPLOG = fopen("components.NWP.tcat","a+");
@@ -258,6 +257,7 @@ void TwoPhase::Initialize()
 	Jwn = Kwn = efawns = 0.0;
 	trJwn = trawn = trRwn = 0.0;
 	euler = Jn = An = Kn = 0.0;
+	wwndnw = 0.0;
 }
 
 
@@ -442,6 +442,9 @@ void TwoPhase::ComputeLocal()
 					// Compute the normal speed of the interface
 					pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris,
 							NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
+
+					for (int pt=0; pt <n_nw_tris; pt++) wwndnw += InterfaceSpeed(r);
+
 				}
 				// wns common curve averages
 				if (n_local_nws_pts > 0){
@@ -451,6 +454,7 @@ void TwoPhase::ComputeLocal()
 					pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z,SDs_x,SDs_y,SDs_z,
 							local_nws_pts,i,j,k,n_local_nws_pts);
 
+
 					pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y,
 							SDs_z, KNwns_values, KGwns_values, KNwns, KGwns,
 							nws_pts, n_nws_pts, i, j, k);
@@ -1075,6 +1079,7 @@ void TwoPhase::Reduce()
 	MPI_Allreduce(&KGwns,&KGwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
 	MPI_Allreduce(&KNwns,&KNwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
 	MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
+	MPI_Allreduce(&wwndnw,&wwndnw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
 	// Phase averages
 	MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
 	MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
@@ -1120,6 +1125,7 @@ void TwoPhase::Reduce()
 	if (awn_global > 0.0){
 		Jwn_global /= awn_global;
 		Kwn_global /= awn_global;
+		wwndnw_global /= awn_global;
 		for (i=0; i<3; i++) vawn_global(i) /= awn_global;
 		for (i=0; i<6; i++)	Gwn_global(i) /= awn_global;
 	}
@@ -1179,7 +1185,7 @@ void TwoPhase::PrintAll(int timestep)
 				Gns_global(0),Gns_global(1),Gns_global(2),Gns_global(3),Gns_global(4),Gns_global(5));	// orientation of ns interface
 		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",
 				Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5));	// orientation of ws interface
-		fprintf(TIMELOG,"%.5g %.5g %.5g ",trawn_global, trJwn_global, trRwn_global);						// Trimmed curvature
+		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",trawn_global, trJwn_global, trRwn_global, wwndnw_global);		// Trimmed curvature
 		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",euler_global, Kn_global, Jn_global, An_global);			// minkowski measures
 		fflush(TIMELOG);
 	}
diff --git a/common/TwoPhase.h b/common/TwoPhase.h
index 82a1c264..898c0d21 100644
--- a/common/TwoPhase.h
+++ b/common/TwoPhase.h
@@ -95,6 +95,7 @@ public:
 	double nwp_volume_global;					// volume for the non-wetting phase
 	double wp_volume_global;					// volume for the wetting phase
 	double As_global;
+	double wwndnw, wwndnw_global;
 	double dEs,dAwn,dAns;						// Global surface energy (calculated by rank=0)
 	DoubleArray van;
 	DoubleArray vaw;

From c44c5e8811bbe58f2f5e8beb3e42d96d5617474c Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 20:02:06 -0500
Subject: [PATCH 17/25] Added new quantity wwndnw to averaging for two-phase

---
 common/TwoPhase.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/common/TwoPhase.cpp b/common/TwoPhase.cpp
index e356d6aa..c4437363 100644
--- a/common/TwoPhase.cpp
+++ b/common/TwoPhase.cpp
@@ -443,7 +443,7 @@ void TwoPhase::ComputeLocal()
 					pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris,
 							NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
 
-					for (int pt=0; pt <n_nw_tris; pt++) wwndnw += InterfaceSpeed(r);
+					for (int p=0; p <n_nw_tris; p++) wwndnw += InterfaceSpeed(p);
 
 				}
 				// wns common curve averages

From a81f9e57767b9b57349e2e10f26863c8e8d17617 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 5 Nov 2015 20:03:28 -0500
Subject: [PATCH 18/25] Added new quantity wwndnw to averaging for two-phase

---
 common/pmmc.h | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/common/pmmc.h b/common/pmmc.h
index bb593baa..73144127 100644
--- a/common/pmmc.h
+++ b/common/pmmc.h
@@ -4347,7 +4347,7 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
 		P.x = 0.33333333333333333*(A.x+B.x+C.x);
 		P.y = 0.33333333333333333*(A.y+B.y+C.y);
 		P.z = 0.33333333333333333*(A.z+B.z+C.z);
-		if (temp > 0.0){
+		if (area > 0.0){
 			x = Px.eval(P);
 			y = Py.eval(P);
 			z = Pz.eval(P);
@@ -4355,6 +4355,7 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
 			if (norm==0.0) norm=1.0;
 			// Compute the interface speed from time derivative and gradient (Level Set Equation)
 			zeta = -Pt.eval(P) / norm;
+
 			//temp = sqrt(temp)/norm; <--- what was I thinking with this? (James)
 
 			// Compute the average

From f53f53adf31d9b29712e9e7eb0dfc55cf1d38054 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Fri, 6 Nov 2015 23:40:36 -0500
Subject: [PATCH 19/25] Added binary file writer

---
 tests/CMakeLists.txt    |   1 +
 tests/ColorToBinary.cpp | 251 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 252 insertions(+)
 create mode 100644 tests/ColorToBinary.cpp

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 4a4459a5..1560278a 100755
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -12,6 +12,7 @@ ADD_LBPM_EXECUTABLE( lbpm_BlobAnalysis )
 ADD_LBPM_EXECUTABLE( TestBubble )
 ADD_LBPM_EXECUTABLE( BasicSimulator )
 ADD_LBPM_EXECUTABLE( ComponentLabel )
+ADD_LBPM_EXECUTABLE( ColorToBinary )
 ADD_LBPM_EXECUTABLE( BlobAnalysis )
 ADD_LBPM_EXECUTABLE( BlobIdentify )
 ADD_LBPM_EXECUTABLE( BlobIdentifyParallel )
diff --git a/tests/ColorToBinary.cpp b/tests/ColorToBinary.cpp
new file mode 100644
index 00000000..e3049a67
--- /dev/null
+++ b/tests/ColorToBinary.cpp
@@ -0,0 +1,251 @@
+// Sequential component labeling for two phase systems
+// Reads parallel simulation data and performs connectivity analysis
+// and averaging on a blob-by-blob basis
+// James E. McClure 2015
+
+#include <iostream>
+#include <math.h>
+#include "analysis/analysis.h"
+#include "common/TwoPhase.h"
+
+
+using namespace std;
+
+
+inline void ReadFromRank(char *FILENAME, DoubleArray &Phase, int nx, int ny, int nz, int iproc, int jproc, int kproc)
+{
+	int i,j,k,q,n,N;
+	int iglobal,jglobal,kglobal;
+	double value;
+	double denA,denB;
+	double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
+	double f10,f11,f12,f13,f14,f15,f16,f17,f18;
+	double vx,vy,vz;
+	
+	N = nx*ny*nz;
+	
+	double *Den, *DistEven, *DistOdd;
+	
+	Den = new double[2*N];
+	DistEven = new double[10*N];
+	DistOdd = new double[9*N];
+
+	ifstream File(FILENAME,ios::binary);
+	for (n=0; n<N; n++){
+		// Write the two density values
+		File.read((char*) &value, sizeof(value));
+		Den[n] = value;
+		//	if (n== 66276)	printf("Density a  = %f \n",value);
+		File.read((char*) &value, sizeof(value));
+		Den[N+n] = value;
+
+		//	if (n== 66276)	printf("Density b  = %f \n",value);
+		// Read the even distributions
+		for (q=0; q<10; q++){
+			File.read((char*) &value, sizeof(value));
+			DistEven[q*N+n] = value;
+		}
+		// Read the odd distributions
+		for (q=0; q<9; q++){
+			File.read((char*) &value, sizeof(value));
+			DistOdd[q*N+n] = value;
+		}
+	}
+	File.close();
+	
+	// Compute the phase field, pressure and velocity
+	for (k=1; k<nz-1; k++){
+		for (j=1; j<ny-1; j++){
+			for (i=1; i<nz-1; i++){
+				//........................................................................
+				n = k*nx*ny+j*nx+i;
+				//........................................................................
+				denA = Den[n];
+				denB = Den[N+n];
+				//........................................................................
+				f0 = DistEven[n];
+				f2 = DistEven[N+n];
+				f4 = DistEven[2*N+n];
+				f6 = DistEven[3*N+n];
+				f8 = DistEven[4*N+n];
+				f10 = DistEven[5*N+n];
+				f12 = DistEven[6*N+n];
+				f14 = DistEven[7*N+n];
+				f16 = DistEven[8*N+n];
+				f18 = DistEven[9*N+n];
+				//........................................................................
+				f1 = DistOdd[n];
+				f3 = DistOdd[1*N+n];
+				f5 = DistOdd[2*N+n];
+				f7 = DistOdd[3*N+n];
+				f9 = DistOdd[4*N+n];
+				f11 = DistOdd[5*N+n];
+				f13 = DistOdd[6*N+n];
+				f15 = DistOdd[7*N+n];
+				f17 = DistOdd[8*N+n];
+				//........................................................................
+				//.................Compute the pressure....................................
+				value = 0.3333333333333333*(f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17);
+				//........................................................................
+				//.................Compute the velocity...................................
+				vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
+				vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
+				vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
+				//........................................................................
+				// save values in global arrays
+				//........................................................................
+				iglobal = iproc*(nx-2)+i;
+				jglobal = jproc*(ny-2)+j;
+				kglobal = kproc*(nz-2)+k;
+				//........................................................................				
+				Phase(iglobal,jglobal,kglobal) = (denA-denB)/(denA+denB);
+				//Pressure(iglobal,jglobal,kglobal) = value;
+				//Vel_x(iglobal,jglobal,kglobal) = vx;
+				//Vel_y(iglobal,jglobal,kglobal) = vy;
+				//Vel_z(iglobal,jglobal,kglobal) = vz;
+				//........................................................................
+			}
+		}
+	}
+	
+	delete Den;
+	delete DistEven;
+	delete DistOdd;
+}
+
+int main(int argc, char **argv)
+{
+	// Initialize MPI
+	int rank,nprocs;
+	MPI_Init(&argc,&argv);
+    MPI_Comm comm = MPI_COMM_WORLD;
+	MPI_Comm_rank(comm,&rank);
+	MPI_Comm_size(comm,&nprocs);
+
+	printf("----------------------------------------------------------\n");
+	printf("Creating single Binary file from restart (8-bit integer)\n");
+	printf("ID=0 (solid) ID=1 (non-wetting), ID=2 (wetting) \n");
+	printf("----------------------------------------------------------\n");
+
+	if (nprocs != 1) INSIST(nprocs == 1,"Error: ComponentLabel --serial case!");
+
+	//.......................................................................
+	int nprocx,nprocy,nprocz;
+	int Nx, Ny, Nz;
+	int nx,ny,nz;
+	int nspheres;
+	double Lx,Ly,Lz;
+	//.......................................................................
+	int i,j,k,n;
+	int iproc,jproc,kproc;
+	//.......................................................................
+	// 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;
+	//.......................................................................
+
+	nx+=2;
+	ny+=2;
+	nz+=2;
+	
+	nprocs = nprocx*nprocy*nprocz;
+	printf("Number of MPI ranks: %i \n", nprocs);
+	int BoundaryCondition=0;
+	Nx = (nx-2)*nprocx;
+	Ny = (ny-2)*nprocy;
+	Nz = (nz-2)*nprocz;
+	Domain Dm(Nx,Ny,Nz,rank,1,1,1,Lx,Ly,Lz,BoundaryCondition);
+	Nx+=2; Ny+=2; Nz+=2;
+	printf("Full domain size: %i x %i x %i  \n", Nx,Ny,Nz);
+	
+	DoubleArray Phase(Nx,Ny,Nz);
+	DoubleArray SignDist(Nx,Ny,Nz);
+	
+	// Filenames used
+	char LocalRankString[8];
+	char LocalRankFilename[40];
+	char BaseFilename[20];
+	                 
+	int proc,iglobal,kglobal,jglobal;
+
+	double * Temp;
+	Temp = new double[nx*ny*nz];
+
+	// read the files and populate main arrays
+	for ( kproc=0; kproc<nprocz; kproc++){
+		for ( jproc=0; jproc<nprocy; jproc++){
+			for ( iproc=0; iproc<nprocx; iproc++){
+				
+				proc = kproc*nprocx*nprocy + jproc*nprocx + iproc;
+
+				sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
+				ReadBinaryFile(LocalRankFilename, Temp, nx*ny*nz);	
+				for (k=1; k<nz-1; k++){
+					for (j=1; j<ny-1; j++){
+						for (i=1; i<nz-1; i++){
+
+							//........................................................................
+							n = k*nx*ny+j*nx+i;
+							//........................................................................
+							iglobal = iproc*(nx-2)+i;
+							jglobal = jproc*(ny-2)+j;
+							kglobal = kproc*(nz-2)+k;
+							//........................................................................
+							SDs(iglobal,jglobal,kglobal) = Temp[n];
+							//........................................................................
+						}
+					}
+				}
+				
+				sprintf(LocalRankFilename,"%s%s","Restart.",LocalRankString);
+				ReadFromRank(LocalRankFilename,Phase,Press,Vel_x,Vel_y,Vel_z,nx,ny,nz,iproc,jproc,kproc);
+				
+			}
+		}
+	}
+	printf("Read %i ranks of %s \n",nprocs,BaseFilename);
+	delete Temp;
+	
+	// Initializing the blob ID
+	for (k=0; k<Nz; k++){
+		for (j=0; j<Ny; j++){
+			for (i=0; i<Nx; i++){
+				n = k*Nx*Ny+j*Nx+i;
+				if (Averages.SDs(i,j,k) < 0.0){
+					// Solid phase 
+					Dm.id[n] = 0;
+				}
+				else if (Averages.Phase(i,j,k) < 0.0){
+					// wetting phase
+					Dm.id[n] = 2;
+				}
+				else {
+					// non-wetting phase
+					Dm.id[n] = 1;
+				}
+			}
+		}
+	}
+
+	FILE *OUTFILE;
+	OUTFILE = fopen("ID.dat","wb");
+	fwrite(&Dm.id,1,Nx*Ny*Nz,OUTFILE);
+	fclose(OUTFILE);
+
+	
+	// ****************************************************
+	MPI_Barrier(comm);
+	MPI_Finalize();
+	// ****************************************************
+}
+

From 725d1957e54e27f64cb76aca0a5dee42082fc72b Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Fri, 6 Nov 2015 23:44:15 -0500
Subject: [PATCH 20/25] Debugging ColorToBinary

---
 tests/ColorToBinary.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/ColorToBinary.cpp b/tests/ColorToBinary.cpp
index e3049a67..8ed57e60 100644
--- a/tests/ColorToBinary.cpp
+++ b/tests/ColorToBinary.cpp
@@ -208,7 +208,7 @@ int main(int argc, char **argv)
 				}
 				
 				sprintf(LocalRankFilename,"%s%s","Restart.",LocalRankString);
-				ReadFromRank(LocalRankFilename,Phase,Press,Vel_x,Vel_y,Vel_z,nx,ny,nz,iproc,jproc,kproc);
+				ReadFromRank(LocalRankFilename,Phase,nx,ny,nz,iproc,jproc,kproc);
 				
 			}
 		}

From 6695d874ff139935d3078db84e7c986dc9d27204 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Fri, 6 Nov 2015 23:46:00 -0500
Subject: [PATCH 21/25] Debugging ColorToBinary

---
 tests/ColorToBinary.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/ColorToBinary.cpp b/tests/ColorToBinary.cpp
index 8ed57e60..8d554ad0 100644
--- a/tests/ColorToBinary.cpp
+++ b/tests/ColorToBinary.cpp
@@ -221,11 +221,11 @@ int main(int argc, char **argv)
 		for (j=0; j<Ny; j++){
 			for (i=0; i<Nx; i++){
 				n = k*Nx*Ny+j*Nx+i;
-				if (Averages.SDs(i,j,k) < 0.0){
+				if (SignDist(i,j,k) < 0.0){
 					// Solid phase 
 					Dm.id[n] = 0;
 				}
-				else if (Averages.Phase(i,j,k) < 0.0){
+				else if (Phase(i,j,k) < 0.0){
 					// wetting phase
 					Dm.id[n] = 2;
 				}

From f1d6bc371da40aac5a7153f9dda32463acc40c4c Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Fri, 6 Nov 2015 23:47:32 -0500
Subject: [PATCH 22/25] Debugging ColorToBinary

---
 tests/ColorToBinary.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/ColorToBinary.cpp b/tests/ColorToBinary.cpp
index 8d554ad0..7e391adb 100644
--- a/tests/ColorToBinary.cpp
+++ b/tests/ColorToBinary.cpp
@@ -201,7 +201,7 @@ int main(int argc, char **argv)
 							jglobal = jproc*(ny-2)+j;
 							kglobal = kproc*(nz-2)+k;
 							//........................................................................
-							SDs(iglobal,jglobal,kglobal) = Temp[n];
+							SignDist(iglobal,jglobal,kglobal) = Temp[n];
 							//........................................................................
 						}
 					}

From 498830e81d357cafa5fde896f6e0f36c75d60416 Mon Sep 17 00:00:00 2001
From: Mark Berrill <berrillma@ornl.gov>
Date: Tue, 10 Nov 2015 14:47:52 -0500
Subject: [PATCH 23/25] Moving writing of vis data to thread

---
 common/Array.hpp               |   2 +-
 common/Communication.hpp       |  44 ++++++------
 common/MPI_Helpers.hpp         |   4 +-
 common/Utilities.cpp           |   6 +-
 tests/lbpm_color_simulator.cpp |  37 +++++++++-
 tests/lbpm_color_simulator.h   | 122 ++++++++++++++++++---------------
 6 files changed, 128 insertions(+), 87 deletions(-)

diff --git a/common/Array.hpp b/common/Array.hpp
index fe85b41c..894382c4 100644
--- a/common/Array.hpp
+++ b/common/Array.hpp
@@ -141,7 +141,7 @@ void Array<TYPE>::resize( const std::vector<size_t>& N )
     // Store the old data
     const size_t ndim_max = sizeof(d_N)/sizeof(size_t);
     std::vector<size_t> N1(ndim_max,1), N2(ndim_max,1);
-    for (size_t d=0; d<d_ndim; d++)
+    for (int d=0; d<d_ndim; d++)
         N1[d] = d_N[d];
     for (size_t d=0; d<N.size(); d++)
         N2[d] = N[d];
diff --git a/common/Communication.hpp b/common/Communication.hpp
index 965dfbc7..56c205cf 100644
--- a/common/Communication.hpp
+++ b/common/Communication.hpp
@@ -207,27 +207,27 @@ template<class TYPE1, class TYPE2>
 void fillHalo<TYPE>::copy( const Array<TYPE1>& src, Array<TYPE2>& dst )
 {
     PROFILE_START("fillHalo::copy",1);
-    ASSERT(src.size(0)==nx||src.size(0)==nx+2*ngx);
-    ASSERT(dst.size(0)==nx||dst.size(0)==nx+2*ngx);
-    bool src_halo = src.size(0)==nx+2*ngx;
-    bool dst_halo = dst.size(0)==nx+2*ngx;
+    ASSERT( (int)src.size(0)==nx || (int)src.size(0)==nx+2*ngx );
+    ASSERT( (int)dst.size(0)==nx || (int)dst.size(0)==nx+2*ngx );
+    bool src_halo = (int)src.size(0)==nx+2*ngx;
+    bool dst_halo = (int)dst.size(0)==nx+2*ngx;
     if ( src_halo ) {
-        ASSERT(src.size(0)==nx+2*ngx);
-        ASSERT(src.size(1)==ny+2*ngy);
-        ASSERT(src.size(2)==nz+2*ngz);
+        ASSERT((int)src.size(0)==nx+2*ngx);
+        ASSERT((int)src.size(1)==ny+2*ngy);
+        ASSERT((int)src.size(2)==nz+2*ngz);
     } else {
-        ASSERT(src.size(0)==nx);
-        ASSERT(src.size(1)==ny);
-        ASSERT(src.size(2)==nz);
+        ASSERT((int)src.size(0)==nx);
+        ASSERT((int)src.size(1)==ny);
+        ASSERT((int)src.size(2)==nz);
     }
     if ( dst_halo ) {
-        ASSERT(dst.size(0)==nx+2*ngx);
-        ASSERT(dst.size(1)==ny+2*ngy);
-        ASSERT(dst.size(2)==nz+2*ngz);
+        ASSERT((int)dst.size(0)==nx+2*ngx);
+        ASSERT((int)dst.size(1)==ny+2*ngy);
+        ASSERT((int)dst.size(2)==nz+2*ngz);
     } else {
-        ASSERT(dst.size(0)==nx);
-        ASSERT(dst.size(1)==ny);
-        ASSERT(dst.size(2)==nz);
+        ASSERT((int)dst.size(0)==nx);
+        ASSERT((int)dst.size(1)==ny);
+        ASSERT((int)dst.size(2)==nz);
     }
     if ( src_halo == dst_halo ) {
         // Src and dst halos match
@@ -235,18 +235,18 @@ void fillHalo<TYPE>::copy( const Array<TYPE1>& src, Array<TYPE2>& dst )
             dst(i) = src(i);
     } else if ( src_halo && !dst_halo ) {
         // Src has halos
-        for (size_t k=0; k<nz; k++) {
-            for (size_t j=0; j<ny; j++) {
-                for (size_t i=0; i<nx; i++) {
+        for (int k=0; k<nz; k++) {
+            for (int j=0; j<ny; j++) {
+                for (int i=0; i<nx; i++) {
                     dst(i,j,k) = src(i+ngx,j+ngy,k+ngz);
                 }
             }
         }
     } else if ( !src_halo && dst_halo ) {
         // Dst has halos
-        for (size_t k=0; k<nz; k++) {
-            for (size_t j=0; j<ny; j++) {
-                for (size_t i=0; i<nx; i++) {
+        for (int k=0; k<nz; k++) {
+            for (int j=0; j<ny; j++) {
+                for (int i=0; i<nx; i++) {
                     dst(i+ngx,j+ngy,k+ngz) = src(i,j,k);
                 }
             }
diff --git a/common/MPI_Helpers.hpp b/common/MPI_Helpers.hpp
index 53708368..85261cf1 100644
--- a/common/MPI_Helpers.hpp
+++ b/common/MPI_Helpers.hpp
@@ -27,7 +27,7 @@ void pack( const std::vector<TYPE>& rhs, char *buffer )
     size_t size = rhs.size();
     memcpy(buffer,&size,sizeof(size_t));
     size_t pos = sizeof(size_t);
-    for (int i=0; i<rhs.size(); i++) {
+    for (size_t i=0; i<rhs.size(); i++) {
         pack(rhs[i],&buffer[pos]);
         pos += packsize(rhs[i]);
     }
@@ -40,7 +40,7 @@ void unpack( std::vector<TYPE>& data, const char *buffer )
     data.clear();
     data.resize(size);
     size_t pos = sizeof(size_t);
-    for (int i=0; i<data.size(); i++) {
+    for (size_t i=0; i<data.size(); i++) {
         unpack(data[i],&buffer[pos]);
         pos += packsize(data[i]);
     }
diff --git a/common/Utilities.cpp b/common/Utilities.cpp
index 55369141..3e8eca8b 100644
--- a/common/Utilities.cpp
+++ b/common/Utilities.cpp
@@ -236,9 +236,9 @@ size_t Utilities::getMemoryUsage()
     size_t N_bytes = 0;
     #if defined(USE_LINUX)
         struct mallinfo meminfo = mallinfo();
-        size_t size_hblkhd = static_cast<size_t>( meminfo.hblkhd );
-        size_t size_uordblks = static_cast<size_t>( meminfo.uordblks );
-        N_bytes = static_cast<size_t>( size_hblkhd + size_uordblks );
+        size_t size_hblkhd = static_cast<unsigned int>( meminfo.hblkhd );
+        size_t size_uordblks = static_cast<unsigned int>( meminfo.uordblks );
+        N_bytes = size_hblkhd + size_uordblks;
     #elif defined(USE_MAC)
         struct task_basic_info t_info;
         mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT;
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index 645ae670..c84b6259 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -7,6 +7,7 @@
 #include <fstream>
 
 #include "common/ScaLBL.h"
+#include "common/Communication.h"
 #include "common/TwoPhase.h"
 #include "common/MPI_Helpers.h"
 
@@ -706,6 +707,37 @@ int main(int argc, char **argv)
         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::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::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::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::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");
 	int timestep = -1;
@@ -714,6 +746,7 @@ int main(int argc, char **argv)
     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");
 
 		//*************************************************************************
@@ -837,8 +870,8 @@ int main(int argc, char **argv)
         // Run the analysis, blob identification, and write restart files
         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,tpool,work_ids);
-
+            LocalRestartFile,meshData,fillData,tpool,work_ids);
+        PROFILE_SAVE("lbpm_color_simulator",false);
 	}
     tpool.wait_pool_finished();
     PROFILE_STOP("Loop");
diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index 29e58c03..80c97515 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -1,21 +1,26 @@
 // Run the analysis, blob identification, and write restart files
 #include "common/Array.h"
 #include "common/Communication.h"
+#include "common/MPI_Helpers.h"
+#include "IO/MeshDatabase.h"
 
+//#define ANALYSIS_INTERVAL 6
 #define ANALYSIS_INTERVAL 1000
 #define BLOBID_INTERVAL 250
 
 enum AnalysisType{ AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02, 
-    CopyAverages=0x04, CalcDist=0x08, CreateRestart=0x10 };
+    CopyAverages=0x04, CalcDist=0x08, CreateRestart=0x10, WriteVis=0x20 };
 
 
 // Structure used to store ids
 struct AnalysisWaitIdStruct {
     ThreadPool::thread_id_t blobID;
     ThreadPool::thread_id_t analysis;
+    ThreadPool::thread_id_t vis;
     ThreadPool::thread_id_t restart;
 };
 
+
 // Helper class to write the restart file from a seperate thread
 class WriteRestartWorkItem: public ThreadPool::WorkItem
 {
@@ -28,8 +33,6 @@ public:
         PROFILE_START("Save Checkpoint",1);
         WriteCheckpoint(filename,cDen.get(),cDistEven.get(),cDistOdd.get(),N);
         PROFILE_STOP("Save Checkpoint",1);
-
-        PROFILE_SAVE("lbpm_color_simulator",1);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     };
 private:
@@ -39,6 +42,7 @@ private:
     const int N;
 };
 
+
 // Helper class to compute the blob ids
 static const std::string id_map_filename = "lbpm_id_map.txt";
 typedef std::shared_ptr<std::pair<int,IntArray> > BlobIDstruct;
@@ -122,71 +126,42 @@ private:
     BlobIDList new_list;
 };
 
+
+// Helper class to write the vis file from a thread
 class WriteVisWorkItem: public ThreadPool::WorkItem
 {
 public:
-	WriteVisWorkItem(AnalysisType type_, TwoPhase& Averages_):
-	type(type_), Averages(Averages_) { }
-
+    WriteVisWorkItem( int timestep_, std::vector<IO::MeshDataStruct>& visData_,
+        TwoPhase& Avgerages_, fillHalo<double>& fillData_ ):
+        timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_) {}
     virtual void run() {
         ThreadPool::WorkItem::d_state = 1;  // Change state to in progress
+        PROFILE_START("Save Vis",1);
+        ASSERT(visData[0].vars[0]->name=="phase");
+        ASSERT(visData[0].vars[1]->name=="Pressure");
+        ASSERT(visData[0].vars[2]->name=="SignDist");
+        ASSERT(visData[0].vars[3]->name=="BlobID");
+        Array<double>& PhaseData = visData[0].vars[0]->data;
+        Array<double>& PressData = visData[0].vars[1]->data;
+        Array<double>& SignData  = visData[0].vars[2]->data;
+        Array<double>& BlobData  = visData[0].vars[3]->data;
+        fillData.copy(Averages.SDn,PhaseData);
+        fillData.copy(Averages.Press,PressData);
+        fillData.copy(Averages.SDs,SignData);
+        fillData.copy(Averages.Label_NWP,BlobData);
         MPI_Comm newcomm;
         MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
-        // Write  VisIT files
-        PROFILE_START("Save Vis",1);
-        // Create the MeshDataStruct
-        Nx=Averages.Dm.Nx;
-        Ny=Averages.Dm.Ny;
-        Nz=Averages.Dm.Nz;
-        Lx=Averages.Dm.Lx;
-        Ly=Averages.Dm.Ly;
-        Lz=Averages.Dm.Lz;
-        fillHalo<double> fillData(newcomm,Averages.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(Averages.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::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::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::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::VolumeVariable;
-        BlobIDVar->dim = 1;
-        BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
-        meshData[0].vars.push_back(BlobIDVar);
-
-        fillData.copy(Averages.SDn,PhaseVar->data);
-        fillData.copy(Averages.SDs,SignDistVar->data);
-        fillData.copy(Averages.Label_NWP,BlobIDVar->data);
-        IO::writeData( 0, meshData, 2, newcomm );
-
+        IO::writeData( timestep, visData, 2, newcomm );
         MPI_Comm_free(&newcomm);
         PROFILE_STOP("Save Vis",1);
-
-        PROFILE_SAVE("lbpm_color_simulator",1);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     };
 private:
     WriteVisWorkItem();
-    AnalysisType type;
+    int timestep;
+    std::vector<IO::MeshDataStruct>& visData;
     TwoPhase& Averages;
-    int Nx,Ny,Nz;
-    double Lx,Ly,Lz;
+    fillHalo<double>& fillData;
 };
 
 
@@ -225,6 +200,7 @@ public:
             Averages.PrintComponents(timestep);
             PROFILE_STOP("Compute dist",1);
         }
+        PROFILE_SAVE("lbpm_color_simulator",false);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     }
 private:
@@ -237,6 +213,8 @@ private:
     double beta;
 };
 
+
+
 // Function to start the analysis
 void run_analysis( int timestep, int restart_interval, 
     const RankInfoStruct& rank_info, TwoPhase& Averages,
@@ -244,7 +222,8 @@ void run_analysis( int timestep, int restart_interval,
     int Nx, int Ny, int Nz, bool pBC, double beta, double err,
     const double *Phi, double *Pressure, const double *Velocity, 
     const char *ID, const double *f_even, const double *f_odd, const double *Den, 
-    const char *LocalRestartFile, ThreadPool& tpool, AnalysisWaitIdStruct& wait )
+    const char *LocalRestartFile, std::vector<IO::MeshDataStruct>& visData, fillHalo<double>& fillData,
+    ThreadPool& tpool, AnalysisWaitIdStruct& wait )
 {
     int N = Nx*Ny*Nz;
 
@@ -266,7 +245,6 @@ void run_analysis( int timestep, int restart_interval,
             type = static_cast<AnalysisType>( type | IdentifyBlobs );
         }
     #endif
-
     if ( timestep%ANALYSIS_INTERVAL == 0 ) {
         // Copy the averages to the CPU (and identify blobs)
         type = static_cast<AnalysisType>( type | CopyAverages );
@@ -280,6 +258,12 @@ void run_analysis( int timestep, int restart_interval,
         // Write the restart file
         type = static_cast<AnalysisType>( type | CreateRestart );
     }
+    if (timestep%restart_interval == 0) {
+        // Write the visualization data
+        type = static_cast<AnalysisType>( type | WriteVis );
+        type = static_cast<AnalysisType>( type | CopyAverages );
+        type = static_cast<AnalysisType>( type | IdentifyBlobs );
+    }
     
     // Return if we are not doing anything
     if ( type == AnalyzeNone )
@@ -305,14 +289,22 @@ void run_analysis( int timestep, int restart_interval,
     }
     if ( (type&CopyAverages) != 0 ) {
         // Copy the members of Averages to the cpu (phase was copied above)
+        // Wait 
+        PROFILE_START("Copy-Wait",1);
+        tpool.wait(wait.analysis);
+        tpool.wait(wait.vis);   // Make sure we are done using analysis before modifying
+        PROFILE_STOP("Copy-Wait",1);
+        PROFILE_START("Copy-Pressure",1);
         ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
         memcpy(Averages.Phase.get(),phase->get(),N*sizeof(double));
         DeviceBarrier();
+        PROFILE_STOP("Copy-Pressure",1);
+        PROFILE_START("Copy-Averages",1);
         CopyToHost(Averages.Press.get(),Pressure,N*sizeof(double));
         CopyToHost(Averages.Vel_x.get(),&Velocity[0],N*sizeof(double));
         CopyToHost(Averages.Vel_y.get(),&Velocity[N],N*sizeof(double));
         CopyToHost(Averages.Vel_z.get(),&Velocity[2*N],N*sizeof(double));
-
+        PROFILE_STOP("Copy-Averages",1);
     }
     std::shared_ptr<double> cDen, cDistEven, cDistOdd;
     if ( (type&CreateRestart) != 0 ) {
@@ -349,6 +341,7 @@ void run_analysis( int timestep, int restart_interval,
             type,timestep,Averages,last_index,last_id_map,beta);
         work->add_dependency(wait.blobID);
         work->add_dependency(wait.analysis);
+        work->add_dependency(wait.vis);     // Make sure we are done using analysis before modifying
         wait.analysis = tpool.add_work(work);
     }
 
@@ -362,12 +355,27 @@ void run_analysis( int timestep, int restart_interval,
         } else {
             // Not clear yet
         }
+        // Wait for previous restart files to finish writing (not necessary, but helps to ensure memory usage is limited)
+        tpool.wait(wait.restart);
         // Write the restart file (using a seperate thread)
         WriteRestartWorkItem *work = new WriteRestartWorkItem(LocalRestartFile,cDen,cDistEven,cDistOdd,N);
         work->add_dependency(wait.restart);
         wait.restart = tpool.add_work(work);
     }
+
+    // Save the results for visualization
+    if ( (type&CreateRestart) != 0 ) {
+        // Wait for previous restart files to finish writing (not necessary, but helps to ensure memory usage is limited)
+        tpool.wait(wait.vis);
+        // Write the vis files
+        ThreadPool::WorkItem *work = new WriteVisWorkItem( timestep, visData, Averages, fillData );
+        work->add_dependency(wait.blobID);
+        work->add_dependency(wait.analysis);
+        work->add_dependency(wait.vis);
+        wait.vis = tpool.add_work(work);
+    }
     PROFILE_STOP("start_analysis");
 }
 
 
+

From 09c5d2ac21cf42c77b57deb09a5db507ae8d11ab Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Tue, 10 Nov 2015 16:04:44 -0500
Subject: [PATCH 24/25] Binary file aggregation

---
 tests/ColorToBinary.cpp | 20 ++++++++++++++++----
 1 file changed, 16 insertions(+), 4 deletions(-)

diff --git a/tests/ColorToBinary.cpp b/tests/ColorToBinary.cpp
index 7e391adb..cedd6b84 100644
--- a/tests/ColorToBinary.cpp
+++ b/tests/ColorToBinary.cpp
@@ -217,21 +217,23 @@ int main(int argc, char **argv)
 	delete Temp;
 	
 	// Initializing the blob ID
+	char *PhaseID;
+	PhaseID = new char (Nx*Ny*Nz);
 	for (k=0; k<Nz; k++){
 		for (j=0; j<Ny; j++){
 			for (i=0; i<Nx; i++){
 				n = k*Nx*Ny+j*Nx+i;
 				if (SignDist(i,j,k) < 0.0){
 					// Solid phase 
-					Dm.id[n] = 0;
+					PhaseID[n] = 0;
 				}
 				else if (Phase(i,j,k) < 0.0){
 					// wetting phase
-					Dm.id[n] = 2;
+					PhaseID[n] = 2;
 				}
 				else {
 					// non-wetting phase
-					Dm.id[n] = 1;
+					PhaseID[n] = 1;
 				}
 			}
 		}
@@ -239,7 +241,17 @@ int main(int argc, char **argv)
 
 	FILE *OUTFILE;
 	OUTFILE = fopen("ID.dat","wb");
-	fwrite(&Dm.id,1,Nx*Ny*Nz,OUTFILE);
+	fwrite(Dm.id,1,Nx*Ny*Nz,OUTFILE);
+	fclose(OUTFILE);
+
+	FILE *OUTFILE;
+	OUTFILE = fopen("Phase.dat","wb");
+	fwrite(Phase,8,Nx*Ny*Nz,OUTFILE);
+	fclose(OUTFILE);
+
+	FILE *OUTFILE;
+	OUTFILE = fopen("SignDist.dat","wb");
+	fwrite(SignDist,8,Nx*Ny*Nz,OUTFILE);
 	fclose(OUTFILE);
 
 	

From 6dd8708204a70cac2e13e8f366064252feef29e7 Mon Sep 17 00:00:00 2001
From: Mark Berrill <berrillma@ornl.gov>
Date: Tue, 10 Nov 2015 16:42:55 -0500
Subject: [PATCH 25/25] Minor performance optimizations for
 lbpm_color_simulator

---
 tests/ColorToBinary.cpp        |  6 ++----
 tests/lbpm_color_simulator.cpp |  5 ++++-
 tests/lbpm_color_simulator.h   | 11 +++++------
 3 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/tests/ColorToBinary.cpp b/tests/ColorToBinary.cpp
index cedd6b84..9d21578d 100644
--- a/tests/ColorToBinary.cpp
+++ b/tests/ColorToBinary.cpp
@@ -244,14 +244,12 @@ int main(int argc, char **argv)
 	fwrite(Dm.id,1,Nx*Ny*Nz,OUTFILE);
 	fclose(OUTFILE);
 
-	FILE *OUTFILE;
 	OUTFILE = fopen("Phase.dat","wb");
-	fwrite(Phase,8,Nx*Ny*Nz,OUTFILE);
+	fwrite(Phase.get(),8,Nx*Ny*Nz,OUTFILE);
 	fclose(OUTFILE);
 
-	FILE *OUTFILE;
 	OUTFILE = fopen("SignDist.dat","wb");
-	fwrite(SignDist,8,Nx*Ny*Nz,OUTFILE);
+	fwrite(SignDist.get(),8,Nx*Ny*Nz,OUTFILE);
 	fclose(OUTFILE);
 
 	
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index c84b6259..d0f6ad36 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -871,7 +871,10 @@ int main(int argc, char **argv)
         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);
-        PROFILE_SAVE("lbpm_color_simulator",false);
+
+        // Save the timers
+        if ( timestep%50==0 )
+            PROFILE_SAVE("lbpm_color_simulator",1);
 	}
     tpool.wait_pool_finished();
     PROFILE_STOP("Loop");
diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index 80c97515..2cc5bec6 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -200,7 +200,6 @@ public:
             Averages.PrintComponents(timestep);
             PROFILE_STOP("Compute dist",1);
         }
-        PROFILE_SAVE("lbpm_color_simulator",false);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     }
 private:
@@ -290,16 +289,16 @@ void run_analysis( int timestep, int restart_interval,
     if ( (type&CopyAverages) != 0 ) {
         // Copy the members of Averages to the cpu (phase was copied above)
         // Wait 
+        PROFILE_START("Copy-Pressure",1);
+        ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
+        DeviceBarrier();
+        PROFILE_STOP("Copy-Pressure",1);
         PROFILE_START("Copy-Wait",1);
         tpool.wait(wait.analysis);
         tpool.wait(wait.vis);   // Make sure we are done using analysis before modifying
         PROFILE_STOP("Copy-Wait",1);
-        PROFILE_START("Copy-Pressure",1);
-        ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
-        memcpy(Averages.Phase.get(),phase->get(),N*sizeof(double));
-        DeviceBarrier();
-        PROFILE_STOP("Copy-Pressure",1);
         PROFILE_START("Copy-Averages",1);
+        memcpy(Averages.Phase.get(),phase->get(),N*sizeof(double));
         CopyToHost(Averages.Press.get(),Pressure,N*sizeof(double));
         CopyToHost(Averages.Vel_x.get(),&Velocity[0],N*sizeof(double));
         CopyToHost(Averages.Vel_y.get(),&Velocity[N],N*sizeof(double));