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/TwoPhase.cpp b/common/TwoPhase.cpp
index df6abae7..c4437363 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+");
@@ -175,38 +174,41 @@ 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];
-		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);
-*/
-	for (int k=0; k<Nz; k++){
+	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++){
 				DistData(i,j,k) = ColorData(i,j,k);
 			}
 		}
 	}
+	*/
 //	for (int n=0; n<Nx*Ny*Nz; n++)	DistData[n] = ColorData[n];
 }
 
@@ -255,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;
 }
 
 
@@ -303,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);
 	//...........................................................................
@@ -437,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 p=0; p <n_nw_tris; p++) wwndnw += InterfaceSpeed(p);
+
 				}
 				// wns common curve averages
 				if (n_local_nws_pts > 0){
@@ -446,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);
@@ -1070,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);
@@ -1115,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;
 	}
@@ -1174,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;
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/common/pmmc.h b/common/pmmc.h
index 1b0720c3..73144127 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,23 +4342,27 @@ 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);
 		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);
 			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;
 
-			AvgVel(0) += temp*zeta*x;
-			AvgVel(1) += temp*zeta*y;
-			AvgVel(2) += temp*zeta*z;
+			//temp = sqrt(temp)/norm; <--- what was I thinking with this? (James)
+
+			// Compute the average
+			AvgVel(0) += area*zeta*x;
+			AvgVel(1) += area*zeta*y;
+			AvgVel(2) += area*zeta*z;
+			AvgSpeed(r) = zeta*area;
 		}
 	}
 	//.............................................................................
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..9d21578d
--- /dev/null
+++ b/tests/ColorToBinary.cpp
@@ -0,0 +1,261 @@
+// 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;
+							//........................................................................
+							SignDist(iglobal,jglobal,kglobal) = Temp[n];
+							//........................................................................
+						}
+					}
+				}
+				
+				sprintf(LocalRankFilename,"%s%s","Restart.",LocalRankString);
+				ReadFromRank(LocalRankFilename,Phase,nx,ny,nz,iproc,jproc,kproc);
+				
+			}
+		}
+	}
+	printf("Read %i ranks of %s \n",nprocs,BaseFilename);
+	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 
+					PhaseID[n] = 0;
+				}
+				else if (Phase(i,j,k) < 0.0){
+					// wetting phase
+					PhaseID[n] = 2;
+				}
+				else {
+					// non-wetting phase
+					PhaseID[n] = 1;
+				}
+			}
+		}
+	}
+
+	FILE *OUTFILE;
+	OUTFILE = fopen("ID.dat","wb");
+	fwrite(Dm.id,1,Nx*Ny*Nz,OUTFILE);
+	fclose(OUTFILE);
+
+	OUTFILE = fopen("Phase.dat","wb");
+	fwrite(Phase.get(),8,Nx*Ny*Nz,OUTFILE);
+	fclose(OUTFILE);
+
+	OUTFILE = fopen("SignDist.dat","wb");
+	fwrite(SignDist.get(),8,Nx*Ny*Nz,OUTFILE);
+	fclose(OUTFILE);
+
+	
+	// ****************************************************
+	MPI_Barrier(comm);
+	MPI_Finalize();
+	// ****************************************************
+}
+
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index 7685aaff..d0f6ad36 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;
@@ -575,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);
@@ -624,6 +627,26 @@ 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){
+		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;
+		// 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);
 
@@ -684,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;
@@ -692,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");
 
 		//*************************************************************************
@@ -787,6 +842,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 += 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);
@@ -794,12 +866,15 @@ 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,
-            LocalRestartFile,tpool,work_ids);
+            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");
@@ -833,40 +908,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 91d24ee2..43351c4b 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -1,17 +1,22 @@
 // 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 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;
 };
 
@@ -28,7 +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:
@@ -123,6 +127,44 @@ private:
 };
 
 
+// Helper class to write the vis file from a thread
+class WriteVisWorkItem: public ThreadPool::WorkItem
+{
+public:
+    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);
+        IO::writeData( timestep, visData, 2, newcomm );
+        MPI_Comm_free(&newcomm);
+        PROFILE_STOP("Save Vis",1);
+        ThreadPool::WorkItem::d_state = 2;  // Change state to finished
+    };
+private:
+    WriteVisWorkItem();
+    int timestep;
+    std::vector<IO::MeshDataStruct>& visData;
+    TwoPhase& Averages;
+    fillHalo<double>& fillData;
+};
+
+
 // Helper class to run the analysis from within a thread
 // Note: Averages will be modified after the constructor is called
 class AnalysisWorkItem: public ThreadPool::WorkItem
@@ -170,6 +212,8 @@ private:
     double beta;
 };
 
+
+
 // Function to start the analysis
 void run_analysis( int timestep, int restart_interval, 
     const RankInfoStruct& rank_info, TwoPhase& Averages,
@@ -177,7 +221,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;
 
@@ -191,7 +236,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 +244,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 );
@@ -213,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 )
@@ -238,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-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-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-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));
         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 ) {
@@ -282,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);
     }
 
@@ -295,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");
 }
 
 
+