From 19ada29a187bb44d27a6c39397329f74c84274f6 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Tue, 19 Mar 2019 20:24:28 -0400
Subject: [PATCH 01/53] summit script

---
 sample_scripts/configure_summit | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/sample_scripts/configure_summit b/sample_scripts/configure_summit
index 5f348f8a..a5b79ba4 100755
--- a/sample_scripts/configure_summit
+++ b/sample_scripts/configure_summit
@@ -11,9 +11,7 @@ cmake                                       \
     -D CMAKE_BUILD_TYPE:STRING=Release      \
     -D CMAKE_C_COMPILER:PATH=mpicc          \
     -D CMAKE_CXX_COMPILER:PATH=mpicxx       \
-    -D CMAKE_C_FLAGS="-std=c++11"                     \
-    -D CMAKE_CXX_FLAGS="-std=c++11"                   \
-    -D CMAKE_CXX_STANDARD=11                \
+    -D CMAKE_CXX_STANDARD=14                \
     -D USE_CUDA=1                           \
         -D CMAKE_CUDA_FLAGS="-arch sm_70 -Xptxas=-v -Xptxas -dlcm=cg -lineinfo" \
         -D CMAKE_CUDA_HOST_COMPILER="/sw/summit/gcc/6.4.0/bin/gcc" \
@@ -21,6 +19,12 @@ cmake                                       \
         -D MPI_COMPILER:BOOL=TRUE           \
         -D MPIEXEC=mpirun                   \
         -D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
+    -D USE_HDF5=1                        \
+        -D HDF5_DIRECTORY="/ccs/proj/bip176/TPLs/hdf5/"		   \
+    	-D HDF5_LIB="/ccs/proj/bip176/TPLs/hdf5/lib/libhdf5.a"	   \
+     -D USE_SILO=1                        \
+        -D SILO_LIB="/ccs/proj/bip176/TPLs/silo/lib/libsiloh5.a"    \
+   	-D SILO_DIRECTORY="/ccs/proj/bip176/TPLs/silo/"		     \
     -D USE_DOXYGEN:BOOL=false		        \
     -D USE_TIMER=0			                \
     ~/LBPM-WIA

From 72f7eb7c6f7c2e4dab2aa814cf38fe8e52549be5 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Sun, 24 Mar 2019 20:48:34 -0400
Subject: [PATCH 02/53] fix liekly bug in TestColorBubble

---
 models/ColorModel.cpp     | 2 +-
 tests/TestColorBubble.cpp | 1 +
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 208c1af5..419fe568 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -205,7 +205,7 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
 	for (int idx=0; idx<NLABELS; idx++)		label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
 
 	if (rank==0){
-		printf("Components labels: %lu \n",NLABELS);
+		printf("Component labels: %lu \n",NLABELS);
 		for (unsigned int idx=0; idx<NLABELS; idx++){
 			VALUE=LabelList[idx];
 			AFFINITY=AffinityList[idx];
diff --git a/tests/TestColorBubble.cpp b/tests/TestColorBubble.cpp
index 8025176b..10f343b1 100644
--- a/tests/TestColorBubble.cpp
+++ b/tests/TestColorBubble.cpp
@@ -49,6 +49,7 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
 					ColorModel.Mask->id[n]=1;
 					ColorModel.Mask->id[n]=1;
 				}
+				ColorModel.id[n] = ColorModel.Mask->id[n];
 			}
 		}
 	}

From ed48ffa264ce78fd6ead9e790f5db44eec44e7c9 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Mon, 25 Mar 2019 07:03:16 -0400
Subject: [PATCH 03/53] reconsidering blob removal rule

---
 models/ColorModel.cpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 419fe568..73b01302 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -760,11 +760,11 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 	}
 
 
-	if (volume_connected < 0.025*volume_initial){
+/*	if (volume_connected < 0.025*volume_initial){
 		// if connected volume is less than 2.5% just delete the whole thing
 		if (rank==0) printf("Connected region has shrunk to less than 2.5% of total fluid volume (remove the whole thing) \n");
 	}
-	else {
+	else { */
 		if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
 		double target_delta_volume_incremental = target_delta_volume;
 		if (fabs(target_delta_volume) > 0.01*volume_initial)  
@@ -799,7 +799,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 			}
 		}
 		fillDouble.fill(phase);
-	}
+	//}
 
 	count = 0.f;
 	for (int k=1; k<Nz-1; k++){

From 24f5fff3f37b6f1afc40dbcfa4768254e1b66dcc Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Mon, 25 Mar 2019 10:34:56 -0400
Subject: [PATCH 04/53] remove shared_ptr from subphase include

---
 analysis/SubPhase.h | 2 --
 1 file changed, 2 deletions(-)

diff --git a/analysis/SubPhase.h b/analysis/SubPhase.h
index ac5cbf76..b2efadf3 100644
--- a/analysis/SubPhase.h
+++ b/analysis/SubPhase.h
@@ -11,8 +11,6 @@
 #include "analysis/analysis.h"
 #include "analysis/distance.h"
 #include "analysis/Minkowski.h"
-
-#include "shared_ptr.h"
 #include "common/Utilities.h"
 #include "common/MPI_Helpers.h"
 #include "IO/MeshDatabase.h"

From d0a516ecf6a8ea5dc4d0620adda3760573247dcc Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Mon, 25 Mar 2019 10:50:58 -0400
Subject: [PATCH 05/53] reverse flow direction if connected pathway is removed

---
 models/ColorModel.cpp | 21 ++++++++++++++++-----
 models/ColorModel.h   |  1 +
 2 files changed, 17 insertions(+), 5 deletions(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 73b01302..173f632e 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -10,7 +10,7 @@ rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),rh
 Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0),
 Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
 {
-
+	REVERSE_FLOW_DIRECTION = false;
 }
 ScaLBL_ColorModel::~ScaLBL_ColorModel(){
 
@@ -621,11 +621,15 @@ void ScaLBL_ColorModel::Run(){
 					//  flow reversal criteria based on fractional flow rate
 					if (delta_volume_target < 0.0 &&
 							volA*flow_rate_A/(volA*flow_rate_A+volB*flow_rate_B) < RESIDUAL_ENDPOINT_THRESHOLD){
-						delta_volume_target *= (-1.0);
+						REVERSE_FLOW_DIRECTION = true;
 					}
 					else if (delta_volume_target > 0.0 &&
 							volB*flow_rate_B/(volA*flow_rate_A+volB*flow_rate_B) < RESIDUAL_ENDPOINT_THRESHOLD){
+						REVERSE_FLOW_DIRECTION = true;
+					}
+					if ( REVERSE_FLOW_DIRECTION ){
 						delta_volume_target *= (-1.0);
+						REVERSE_FLOW_DIRECTION = false;
 					}
 				}
 				else{
@@ -651,6 +655,12 @@ void ScaLBL_ColorModel::Run(){
 					MORPH_ADAPT = false;
 					CURRENT_STEADY_TIMESTEPS=0;
 				}
+				if ( REVERSE_FLOW_DIRECTION ){
+					// flow direction will reverse after next steady point
+					MORPH_ADAPT = false;
+					CURRENT_STEADY_TIMESTEPS=0;
+				}
+
 				MPI_Barrier(comm);
 			}
 			morph_timesteps += analysis_interval;
@@ -760,11 +770,12 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 	}
 
 
-/*	if (volume_connected < 0.025*volume_initial){
+	if (volume_connected < 0.025*volume_initial){
 		// if connected volume is less than 2.5% just delete the whole thing
 		if (rank==0) printf("Connected region has shrunk to less than 2.5% of total fluid volume (remove the whole thing) \n");
+		REVERSE_FLOW_DIRECTION = true;
 	}
-	else { */
+	else { 
 		if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
 		double target_delta_volume_incremental = target_delta_volume;
 		if (fabs(target_delta_volume) > 0.01*volume_initial)  
@@ -799,7 +810,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 			}
 		}
 		fillDouble.fill(phase);
-	//}
+	}
 
 	count = 0.f;
 	for (int k=1; k<Nz-1; k++){
diff --git a/models/ColorModel.h b/models/ColorModel.h
index 76241d80..805d8e22 100644
--- a/models/ColorModel.h
+++ b/models/ColorModel.h
@@ -32,6 +32,7 @@ public:
 	void WriteDebug();
 	
 	bool Restart,pBC;
+	bool REVERSE_FLOW_DIRECTION;
 	int timestep,timestepMax;
 	int BoundaryCondition;
 	double tauA,tauB,rhoA,rhoB,alpha,beta;

From d8f7d2274d6d0e7cafee9b0d42320667190f8a35 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 07:14:16 -0400
Subject: [PATCH 06/53] tweaking endpoint criterion

---
 models/ColorModel.cpp | 71 ++++++++++++++++++++++---------------------
 1 file changed, 37 insertions(+), 34 deletions(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 173f632e..a36166d5 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -652,10 +652,13 @@ void ScaLBL_ColorModel::Run(){
 					CURRENT_STEADY_TIMESTEPS=0;
 				}
 				else if (CURRENT_MORPH_TIMESTEPS > MAX_MORPH_TIMESTEPS) {
+					delta_volume = 0.0;
 					MORPH_ADAPT = false;
 					CURRENT_STEADY_TIMESTEPS=0;
 				}
 				if ( REVERSE_FLOW_DIRECTION ){
+					if (rank==0) printf("*****REVERSE FLOW DIRECTION***** \n");
+					delta_volume = 0.0;
 					// flow direction will reverse after next steady point
 					MORPH_ADAPT = false;
 					CURRENT_STEADY_TIMESTEPS=0;
@@ -770,47 +773,47 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 	}
 
 
-	if (volume_connected < 0.025*volume_initial){
-		// if connected volume is less than 2.5% just delete the whole thing
-		if (rank==0) printf("Connected region has shrunk to less than 2.5% of total fluid volume (remove the whole thing) \n");
+	if (volume_connected < 0.02*volume_initial){
+		// if connected volume is less than 2% just delete the whole thing
+		if (rank==0) printf("Connected region has shrunk to less than 2% of total fluid volume \n");
 		REVERSE_FLOW_DIRECTION = true;
 	}
-	else { 
-		if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
-		double target_delta_volume_incremental = target_delta_volume;
-		if (fabs(target_delta_volume) > 0.01*volume_initial)  
-			target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume);
-		delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental);
 
-		for (int k=0; k<Nz; k++){
-			for (int j=0; j<Ny; j++){
-				for (int i=0; i<Nx; i++){
-					if (phase_distance(i,j,k) < 0.0 ) phase_id(i,j,k) = 0;
-					else 		     				  phase_id(i,j,k) = 1;
-					//if (phase_distance(i,j,k) < 0.0 ) phase(i,j,k) = 1.0;
-				}
-			}
-		}	
+	if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
+	double target_delta_volume_incremental = target_delta_volume;
+	if (fabs(target_delta_volume) > 0.01*volume_initial)  
+		target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume);
+	delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental);
 
-		CalcDist(phase_distance,phase_id,*Dm); // re-calculate distance
-
-		// 5. Update phase indicator field based on new distnace
-		for (int k=0; k<Nz; k++){
-			for (int j=0; j<Ny; j++){
-				for (int i=0; i<Nx; i++){
-					int n = k*Nx*Ny + j*Nx + i;
-					double d = phase_distance(i,j,k);
-					if (Averages->SDs(i,j,k) > 0.f){
-						if (d < 3.f){
-							//phase(i,j,k) = -1.0;
-							phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);	
-						}
-					}
-				} 
+	for (int k=0; k<Nz; k++){
+		for (int j=0; j<Ny; j++){
+			for (int i=0; i<Nx; i++){
+				if (phase_distance(i,j,k) < 0.0 ) phase_id(i,j,k) = 0;
+				else 		     				  phase_id(i,j,k) = 1;
+				//if (phase_distance(i,j,k) < 0.0 ) phase(i,j,k) = 1.0;
 			}
 		}
-		fillDouble.fill(phase);
+	}	
+
+	CalcDist(phase_distance,phase_id,*Dm); // re-calculate distance
+
+	// 5. Update phase indicator field based on new distnace
+	for (int k=0; k<Nz; k++){
+		for (int j=0; j<Ny; j++){
+			for (int i=0; i<Nx; i++){
+				int n = k*Nx*Ny + j*Nx + i;
+				double d = phase_distance(i,j,k);
+				if (Averages->SDs(i,j,k) > 0.f){
+					if (d < 3.f){
+						//phase(i,j,k) = -1.0;
+						phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);	
+					}
+				}
+			} 
+		}
 	}
+	fillDouble.fill(phase);
+
 
 	count = 0.f;
 	for (int k=1; k<Nz-1; k++){

From 2378a1455393800cb90738bbaf892ec8e37d2a46 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 10:17:03 -0400
Subject: [PATCH 07/53] added analysis for connected interface

---
 analysis/SubPhase.cpp | 31 +++++++++++++++++++++++--------
 analysis/SubPhase.h   |  6 ++++--
 2 files changed, 27 insertions(+), 10 deletions(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index 9cc977e5..e59d91ef 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -41,10 +41,12 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 			fprintf(TIMELOG,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");			
 			fprintf(TIMELOG,"Kwc Kwd Kwi Knc Knd Kni ");				// kinetic energy
 			fprintf(TIMELOG,"Vwc Awc Hwc Xwc ");					 	// wc region 
-			fprintf(TIMELOG,"Vwd Awd Hwd Xwd Nwd ");					 	// wd region
+			fprintf(TIMELOG,"Vwd Awd Hwd Xwd Nwd ");					// wd region
 			fprintf(TIMELOG,"Vnc Anc Hnc Xnc ");					 	// nc region
-			fprintf(TIMELOG,"Vnd And Hnd Xnd Nnd ");					 	// nd region
-			fprintf(TIMELOG,"Vi Ai Hi Xi\n");					 		// interface region 
+			fprintf(TIMELOG,"Vnd And Hnd Xnd Nnd ");					// nd region
+			fprintf(TIMELOG,"Vi Ai Hi Xi ");					 		// interface region 
+			fprintf(TIMELOG,"Vic Aic Hic Xic Nic\n");					// interface region 
+
 			// stress tensor?
 		}
 
@@ -64,10 +66,11 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 		fprintf(TIMELOG,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");			
 		fprintf(TIMELOG,"Kwc Kwd Kwi Knc Knd Kni ");				// kinetic energy
 		fprintf(TIMELOG,"Vwc Awc Hwc Xwc ");					 	// wc region 
-		fprintf(TIMELOG,"Vwd Awd Hwd Xwd Nwd ");					 	// wd region
+		fprintf(TIMELOG,"Vwd Awd Hwd Xwd Nwd ");					// wd region
 		fprintf(TIMELOG,"Vnc Anc Hnc Xnc ");					 	// nc region
-		fprintf(TIMELOG,"Vnd And Hnd Xnd Nnd ");					 	// nd region
-		fprintf(TIMELOG,"Vi Ai Hi Xi\n");					 		// interface region 
+		fprintf(TIMELOG,"Vnd And Hnd Xnd Nnd ");					// nd region
+		fprintf(TIMELOG,"Vi Ai Hi Xi ");					 		// interface region 
+		fprintf(TIMELOG,"Vic Aic Hic Xic Nic\n");					// interface region 
 	}
 }
 
@@ -93,7 +96,8 @@ void SubPhase::Write(int timestep)
 		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
 		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",gnc.V, gnc.A, gnc.H, gnc.X);
 		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",gnd.V, gnd.A, gnd.H, gnd.X, gnd.Nc);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",giwn.V, giwn.A, giwn.H, giwn.X);
+		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",giwn.V, giwn.A, giwn.H, giwn.X);
+		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i\n",giwnc.V, giwnc.A, giwnc.H, giwnc.X, giwnc.Nc);
 		fflush(TIMELOG);
 	}
 	else{
@@ -108,7 +112,8 @@ void SubPhase::Write(int timestep)
 		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);
 		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",nc.V, nc.A, nc.H, nc.X);
 		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",nd.V, nd.A, nd.H, nd.X, nd.Nc);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",iwn.V, iwn.A, iwn.H, iwn.X);
+		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",iwn.V, iwn.A, iwn.H, iwn.X);
+		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",iwnc.V, iwnc.A, iwnc.H, iwnc.X);
 		fflush(TIMELOG);
 	}
 
@@ -429,6 +434,16 @@ void SubPhase::Full(){
 	giwn.A=sumReduce( Dm->Comm, iwn.A);
 	giwn.H=sumReduce( Dm->Comm, iwn.H);
 	giwn.X=sumReduce( Dm->Comm, iwn.X);
+	// measure only the connected part
+	iwnc.Nc = morph_i->MeasureConnectedPathway();
+	iwnc.V = morph_i->V(); 
+	iwnc.A = morph_i->A(); 
+	iwnc.H = morph_i->H(); 
+	iwnc.X = morph_i->X(); 
+	giwnc.V=sumReduce( Dm->Comm, iwnc.V);
+	giwnc.A=sumReduce( Dm->Comm, iwnc.A);
+	giwnc.H=sumReduce( Dm->Comm, iwnc.H);
+	giwnc.X=sumReduce( Dm->Comm, iwnc.X);
 	
 	double vol_nc_bulk = 0.0;
 	double vol_wc_bulk = 0.0;
diff --git a/analysis/SubPhase.h b/analysis/SubPhase.h
index b2efadf3..b0a68347 100644
--- a/analysis/SubPhase.h
+++ b/analysis/SubPhase.h
@@ -35,10 +35,12 @@ private:
 
 class interface{
 public:
+	int Nc;
 	double M,Px,Py,Pz,K;
 	double Mw,Mn,Pnx,Pny,Pnz,Pwx,Pwy,Pwz,Kw,Kn;
 	double V,A,H,X;
 	void reset(){
+		Nc = 0;
 		M=Px=Py=Pz=K=0.0;
 		V=A=H=X=0.0;
 		Mw=Mn=Pnx=Pny=Pnz=Pwx=Pwy=Pwz=Kw=Kn=0.0;
@@ -67,11 +69,11 @@ public:
 	 */
 	// local entities
 	phase wc,wd,wb,nc,nd,nb;
-	interface iwn;
+	interface iwn,iwnc;
 	
 	// global entities
 	phase gwc,gwd,gwb,gnc,gnd,gnb;
-	interface giwn;
+	interface giwn,giwnc;
 	//...........................................................................
     int Nx,Ny,Nz;
 	IntArray PhaseID;		// Phase ID array (solid=0, non-wetting=1, wetting=2)

From ebabdcaef0c47b1a2fc01bfe8e58cd499f663e78 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 10:26:02 -0400
Subject: [PATCH 08/53] fix interface averaging for mass /momentum

---
 analysis/SubPhase.cpp     | 15 +++++++++++++--
 tests/TestColorBubble.cpp | 34 ++++++++++++++++++++--------------
 2 files changed, 33 insertions(+), 16 deletions(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index e59d91ef..e24ffe4f 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -290,7 +290,7 @@ void SubPhase::Full(){
 	if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y;
 	if (Dm->inlet_layers_z > 0) kmin = Dm->inlet_layers_z;
 		
-	nd.reset();	nc.reset(); wd.reset();	wc.reset();	iwn.reset();
+	nd.reset();	nc.reset(); wd.reset();	wc.reset();	iwn.reset();	iwnc.reset();
 
  	Dm->CommunicateMeshHalo(Phi);
 	for (int k=1; k<Nz-1; k++){
@@ -543,7 +543,18 @@ void SubPhase::Full(){
 	gwc.Pz=sumReduce( Dm->Comm, wc.Pz);
 	gwc.K=sumReduce( Dm->Comm, wc.K);
 	gwc.p=sumReduce( Dm->Comm, wc.p);
-
+	
+	giwn.Mn=sumReduce( Dm->Comm, iwn.Mn);
+	giwn.Pnx=sumReduce( Dm->Comm, iwn.Pnx);
+	giwn.Pny=sumReduce( Dm->Comm, iwn.Pny);
+	giwn.Pnz=sumReduce( Dm->Comm, iwn.Pnz);
+	giwn.Kn=sumReduce( Dm->Comm, iwn.Kn);
+	giwn.Mw=sumReduce( Dm->Comm, iwn.Mw);
+	giwn.Pwx=sumReduce( Dm->Comm, iwn.Pwx);
+	giwn.Pwy=sumReduce( Dm->Comm, iwn.Pwy);
+	giwn.Pwz=sumReduce( Dm->Comm, iwn.Pwz);
+	giwn.Kw=sumReduce( Dm->Comm, iwn.Kw);
+	
 	// pressure averaging
 	if (vol_wc_bulk > 0.0)
 		wc.p = wc.p /vol_wc_bulk;
diff --git a/tests/TestColorBubble.cpp b/tests/TestColorBubble.cpp
index 10f343b1..0e6ea25a 100644
--- a/tests/TestColorBubble.cpp
+++ b/tests/TestColorBubble.cpp
@@ -14,13 +14,13 @@ using namespace std;
 inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius){
 	// initialize a bubble
 	int i,j,k,n;
-	int rank = ColorModel.Mask->rank();
-	int nprocx = ColorModel.Mask->nprocx();
-	int nprocy = ColorModel.Mask->nprocy();
-	int nprocz = ColorModel.Mask->nprocz();
-	int Nx = ColorModel.Mask->Nx;
-	int Ny = ColorModel.Mask->Ny;
-	int Nz = ColorModel.Mask->Nz;
+	int rank = ColorModel.Dm->rank();
+	int nprocx = ColorModel.Dm->nprocx();
+	int nprocy = ColorModel.Dm->nprocy();
+	int nprocz = ColorModel.Dm->nprocz();
+	int Nx = ColorModel.Dm->Nx;
+	int Ny = ColorModel.Dm->Ny;
+	int Nz = ColorModel.Dm->Nz;
 	if (rank == 0) cout << "Setting up bubble..." << endl;
 	for (k=0;k<Nz;k++){
 		for (j=0;j<Ny;j++){
@@ -34,14 +34,12 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
 	for (k=0;k<Nz;k++){
 		for (j=0;j<Ny;j++){
 			for (i=0;i<Nx;i++){
-				n = k*Nx*Ny + j*Nz + i;
-				int iglobal= i+(Nx-2)*ColorModel.Mask->iproc();
-				int jglobal= j+(Ny-2)*ColorModel.Mask->jproc();
-				int kglobal= k+(Nz-2)*ColorModel.Mask->kproc();
+				n = k*Nx*Ny + j*Nx + i;
+				double iglobal= double(i+(Nx-2)*ColorModel.Dm->iproc())-double((Nx-2)*nprocx)*0.5;
+				double jglobal= double(j+(Ny-2)*ColorModel.Dm->jproc())-double((Ny-2)*nprocy)*0.5;
+				double kglobal= double(k+(Nz-2)*ColorModel.Dm->kproc())-double((Nz-2)*nprocz)*0.5;
 				// Initialize phase position field for parallel bubble test
-				if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
-						+(jglobal-0.5*(Ny-2)*nprocy)*(jglobal-0.5*(Ny-2)*nprocy)
-						+(kglobal-0.5*(Nz-2)*nprocz)*(kglobal-0.5*(Nz-2)*nprocz) < BubbleRadius*BubbleRadius){
+				if ((iglobal*iglobal)+(jglobal*jglobal)+(kglobal*kglobal) < BubbleRadius*BubbleRadius){
 					ColorModel.Mask->id[n] = 2;
 					ColorModel.Mask->id[n] = 2;
 				}
@@ -50,9 +48,17 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
 					ColorModel.Mask->id[n]=1;
 				}
 				ColorModel.id[n] = ColorModel.Mask->id[n];
+				ColorModel.Dm->id[n] = ColorModel.Mask->id[n];
 			}
 		}
 	}
+	
+	FILE *OUTFILE;
+	char LocalRankFilename[40];
+	sprintf(LocalRankFilename,"Bubble.%05i.raw",rank);
+	OUTFILE = fopen(LocalRankFilename,"wb");
+	fwrite(ColorModel.id,1,Nx*Ny*Nz,OUTFILE);
+	fclose(OUTFILE);
 	// initialize the phase indicator field
 }
 //***************************************************************************************

From 2b84064d0c8eb88598c316b0b748129a5b849552 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 10:50:40 -0400
Subject: [PATCH 09/53] fixing flow rate to match force direction

---
 analysis/SubPhase.cpp | 9 +++++++--
 models/ColorModel.cpp | 9 ++++++---
 2 files changed, 13 insertions(+), 5 deletions(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index e24ffe4f..d5155237 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -208,9 +208,14 @@ void SubPhase::Basic(){
 	gnb.Pz=sumReduce( Dm->Comm, nb.Pz);
 
 	if (Dm->rank() == 0){
+		double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz)
+		double dir_x = Fx/force_mag;
+		double dir_y = Fy/force_mag;
+		double dir_z = Fz/force_mag;
+		
 		double saturation=gwb.V/(gwb.V + gnb.V);
-		double water_flow_rate=gwb.V*sqrt(gwb.Px*gwb.Px + gwb.Py*gwb.Py + gwb.Pz*gwb.Pz)/gwb.M;
-		double not_water_flow_rate=gnb.V*sqrt(gnb.Px*gnb.Px + gnb.Py*gnb.Py + gnb.Pz*gnb.Pz)/gnb.M;
+		double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M;
+		double not_water_flow_rate=gnb.V*sqrt(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M;
 		double total_flow_rate = water_flow_rate + not_water_flow_rate;
 		double fractional_flow= water_flow_rate / total_flow_rate;
 		printf("   water saturation = %f, fractional flow =%f \n",saturation,fractional_flow);
diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index a36166d5..a148538f 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -571,8 +571,12 @@ void ScaLBL_ColorModel::Run(){
 				double muA = rhoA*(tauA-0.5)/3.f; 
 				double muB = rhoB*(tauB-0.5)/3.f;				
 				
-				double flow_rate_A = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z);
-				double flow_rate_B = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
+				double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz)
+				double dir_x = Fx/force_mag;
+				double dir_y = Fy/force_mag;
+				double dir_z = Fz/force_mag;
+				double flow_rate_A = (vA_x*dir_x + vA_y*dir_y + vA_z*dir_z);
+				double flow_rate_B = (vB_x*dir_x + vB_y*dir_y + vB_z*dir_z);
 				double current_saturation = volB/(volA+volB);
 				double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double((Nx-2)*(Ny-2)*(Nz-2)*nprocs));
 
@@ -772,7 +776,6 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 		}
 	}
 
-
 	if (volume_connected < 0.02*volume_initial){
 		// if connected volume is less than 2% just delete the whole thing
 		if (rank==0) printf("Connected region has shrunk to less than 2% of total fluid volume \n");

From cb74f0e9abeb98b15f8a45be73447609da4278f1 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 11:17:56 -0400
Subject: [PATCH 10/53] added max steady timesteps for morph

---
 models/ColorModel.cpp | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index a148538f..9465ed0e 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -415,6 +415,7 @@ void ScaLBL_ColorModel::Run(){
 	bool USE_MORPH = false;
 	int analysis_interval = 1000; 	// number of timesteps in between in situ analysis 
 	int MAX_MORPH_TIMESTEPS = 50000; // maximum number of LBM timesteps to spend in morphological adaptation routine
+	int MAX_STEADY_TIMESTEPS = 200000;
 	int RAMP_TIMESTEPS = 0;//50000;		 // number of timesteps to run initially (to get a reasonable velocity field before other pieces kick in)
 	int morph_interval = 1000000;
 	int CURRENT_MORPH_TIMESTEPS=0;   // counter for number of timesteps spent in  morphological adaptation routine (reset each time)
@@ -452,7 +453,12 @@ void ScaLBL_ColorModel::Run(){
 	if (analysis_db->keyExists( "analysis_interval" )){
 		analysis_interval = analysis_db->getScalar<int>( "analysis_interval" );
 	}
-	
+	if (analysis_db->keyExists( "max_steady_timesteps" )){
+		MAX_STEADY_TIMESTEPS = analysis_db->getScalar<int>( "max_steady_timesteps" );
+	}
+	if (analysis_db->keyExists( "max_morph_timesteps" )){
+		MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" );
+	}
 	if (rank==0){
 		printf("********************************************************\n");
 		printf("No. of timesteps: %i \n", timestepMax);
@@ -582,7 +588,7 @@ void ScaLBL_ColorModel::Run(){
 
 				double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz);
 
-				if (fabs((Ca - Ca_previous)/Ca) < tolerance ){
+				if (fabs((Ca - Ca_previous)/Ca) < tolerance || CURRENT_STEADY_TIMESTEPS > MAX_STEADY_TIMESTEPS){
 					MORPH_ADAPT = true;
 					CURRENT_MORPH_TIMESTEPS=0;
 					delta_volume_target = (volA + volB)*morph_delta; // set target volume change

From 8a60057a2809eb08258b78f6aebc6c904d0fe2a2 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 11:18:43 -0400
Subject: [PATCH 11/53] fix small bug

---
 analysis/SubPhase.cpp | 2 +-
 models/ColorModel.cpp | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index d5155237..c1415ebc 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -208,7 +208,7 @@ void SubPhase::Basic(){
 	gnb.Pz=sumReduce( Dm->Comm, nb.Pz);
 
 	if (Dm->rank() == 0){
-		double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz)
+		double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
 		double dir_x = Fx/force_mag;
 		double dir_y = Fy/force_mag;
 		double dir_z = Fz/force_mag;
diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 9465ed0e..12dd19e1 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -577,7 +577,7 @@ void ScaLBL_ColorModel::Run(){
 				double muA = rhoA*(tauA-0.5)/3.f; 
 				double muB = rhoB*(tauB-0.5)/3.f;				
 				
-				double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz)
+				double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
 				double dir_x = Fx/force_mag;
 				double dir_y = Fy/force_mag;
 				double dir_z = Fz/force_mag;

From 9d60c0f949bf29251f35772118eb22e0ed262eb5 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 11:21:25 -0400
Subject: [PATCH 12/53] fix print bug

---
 models/ColorModel.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 12dd19e1..19d0f9d0 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -784,7 +784,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 
 	if (volume_connected < 0.02*volume_initial){
 		// if connected volume is less than 2% just delete the whole thing
-		if (rank==0) printf("Connected region has shrunk to less than 2% of total fluid volume \n");
+		if (rank==0) printf("Connected region has shrunk to less than 2 %% of total fluid volume \n");
 		REVERSE_FLOW_DIRECTION = true;
 	}
 

From aed56d558717a6b3d025d4e2d0c94415cfb87518 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 13:58:29 -0400
Subject: [PATCH 13/53] different limits on growth/ shrink

---
 analysis/SubPhase.cpp   | 120 ++++++++++++++++++++--------------------
 analysis/SubPhase.h     |   2 +-
 analysis/morphology.cpp |  20 +++++--
 3 files changed, 75 insertions(+), 67 deletions(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index c1415ebc..7b466a24 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -28,24 +28,24 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 
 	//.........................................
 	if (Dm->rank()==0){
-		TIMELOG = fopen("subphase.csv","a+");
-		if (fseek(TIMELOG,0,SEEK_SET) == fseek(TIMELOG,0,SEEK_CUR))
+		SUBPHASE = fopen("subphase.csv","a+");
+		if (ftell(SUBPHASE) == 0)
 		{
 			// If timelog is empty, write a short header to list the averages
-			//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
-			fprintf(TIMELOG,"time rn rw nun nuw Fx Fy Fz iftwn ");				
-			fprintf(TIMELOG,"pwc pwd pnc pnd ");						// pressures 
-			fprintf(TIMELOG,"Mwc Mwd Mwi Mnc Mnd Mni ");				// mass 
-			fprintf(TIMELOG,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x ");	// momentum 
-			fprintf(TIMELOG,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");			
-			fprintf(TIMELOG,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");			
-			fprintf(TIMELOG,"Kwc Kwd Kwi Knc Knd Kni ");				// kinetic energy
-			fprintf(TIMELOG,"Vwc Awc Hwc Xwc ");					 	// wc region 
-			fprintf(TIMELOG,"Vwd Awd Hwd Xwd Nwd ");					// wd region
-			fprintf(TIMELOG,"Vnc Anc Hnc Xnc ");					 	// nc region
-			fprintf(TIMELOG,"Vnd And Hnd Xnd Nnd ");					// nd region
-			fprintf(TIMELOG,"Vi Ai Hi Xi ");					 		// interface region 
-			fprintf(TIMELOG,"Vic Aic Hic Xic Nic\n");					// interface region 
+			//fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n");
+			fprintf(SUBPHASE,"time rn rw nun nuw Fx Fy Fz iftwn ");				
+			fprintf(SUBPHASE,"pwc pwd pnc pnd ");						// pressures 
+			fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni ");				// mass 
+			fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x ");	// momentum 
+			fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");			
+			fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");			
+			fprintf(SUBPHASE,"Kwc Kwd Kwi Knc Knd Kni ");				// kinetic energy
+			fprintf(SUBPHASE,"Vwc Awc Hwc Xwc ");					 	// wc region 
+			fprintf(SUBPHASE,"Vwd Awd Hwd Xwd Nwd ");					// wd region
+			fprintf(SUBPHASE,"Vnc Anc Hnc Xnc ");					 	// nc region
+			fprintf(SUBPHASE,"Vnd And Hnd Xnd Nnd ");					// nd region
+			fprintf(SUBPHASE,"Vi Ai Hi Xi ");					 		// interface region 
+			fprintf(SUBPHASE,"Vic Aic Hic Xic Nic\n");					// interface region 
 
 			// stress tensor?
 		}
@@ -56,21 +56,21 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 		sprintf(LocalRankString,"%05d",Dm->rank());
 		char LocalRankFilename[40];
 		sprintf(LocalRankFilename,"%s%s","subphase.csv.",LocalRankString);
-		TIMELOG = fopen(LocalRankFilename,"a+");
-		//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
-		fprintf(TIMELOG,"time rn rw nun nuw Fx Fy Fz iftwn ");				
-		fprintf(TIMELOG,"pwc pwd pnc pnd ");						// pressures 
-		fprintf(TIMELOG,"Mwc Mwd Mwi Mnc Mnd Mni ");				// mass 
-		fprintf(TIMELOG,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x ");	// momentum 
-		fprintf(TIMELOG,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");			
-		fprintf(TIMELOG,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");			
-		fprintf(TIMELOG,"Kwc Kwd Kwi Knc Knd Kni ");				// kinetic energy
-		fprintf(TIMELOG,"Vwc Awc Hwc Xwc ");					 	// wc region 
-		fprintf(TIMELOG,"Vwd Awd Hwd Xwd Nwd ");					// wd region
-		fprintf(TIMELOG,"Vnc Anc Hnc Xnc ");					 	// nc region
-		fprintf(TIMELOG,"Vnd And Hnd Xnd Nnd ");					// nd region
-		fprintf(TIMELOG,"Vi Ai Hi Xi ");					 		// interface region 
-		fprintf(TIMELOG,"Vic Aic Hic Xic Nic\n");					// interface region 
+		SUBPHASE = fopen(LocalRankFilename,"a+");
+		//fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n");
+		fprintf(SUBPHASE,"time rn rw nun nuw Fx Fy Fz iftwn ");				
+		fprintf(SUBPHASE,"pwc pwd pnc pnd ");						// pressures 
+		fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni ");				// mass 
+		fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x ");	// momentum 
+		fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");			
+		fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");			
+		fprintf(SUBPHASE,"Kwc Kwd Kwi Knc Knd Kni ");				// kinetic energy
+		fprintf(SUBPHASE,"Vwc Awc Hwc Xwc ");					 	// wc region 
+		fprintf(SUBPHASE,"Vwd Awd Hwd Xwd Nwd ");					// wd region
+		fprintf(SUBPHASE,"Vnc Anc Hnc Xnc ");					 	// nc region
+		fprintf(SUBPHASE,"Vnd And Hnd Xnd Nnd ");					// nd region
+		fprintf(SUBPHASE,"Vi Ai Hi Xi ");					 		// interface region 
+		fprintf(SUBPHASE,"Vic Aic Hic Xic Nic\n");					// interface region 
 	}
 }
 
@@ -78,43 +78,43 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 // Destructor
 SubPhase::~SubPhase()
 {
-    if ( TIMELOG!=NULL ) { fclose(TIMELOG); }
+    if ( SUBPHASE!=NULL ) { fclose(SUBPHASE); }
 
 }
 
 void SubPhase::Write(int timestep)
 {
 	if (Dm->rank()==0){
-		fprintf(TIMELOG,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn); 
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",gwc.p, gwd.p, gnc.p, gnd.p);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Px, gwd.Px, giwn.Px, gnc.Px, gnd.Px, giwn.Px);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Py, gwd.Py, giwn.Py, gnc.Py, gnd.Py, giwn.Py);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Pz, gwd.Pz, giwn.Pz, gnc.Pz, gnd.Pz, giwn.Pz);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",gwc.V, gwc.A, gwc.H, gwc.X);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",gnc.V, gnc.A, gnc.H, gnc.X);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",gnd.V, gnd.A, gnd.H, gnd.X, gnd.Nc);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",giwn.V, giwn.A, giwn.H, giwn.X);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i\n",giwnc.V, giwnc.A, giwnc.H, giwnc.X, giwnc.Nc);
-		fflush(TIMELOG);
+		fprintf(SUBPHASE,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn); 
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.p, gwd.p, gnc.p, gnd.p);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Px, gwd.Px, giwn.Px, gnc.Px, gnd.Px, giwn.Px);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Py, gwd.Py, giwn.Py, gnc.Py, gnd.Py, giwn.Py);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Pz, gwd.Pz, giwn.Pz, gnc.Pz, gnd.Pz, giwn.Pz);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.V, gwc.A, gwc.H, gwc.X);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gnc.V, gnc.A, gnc.H, gnc.X);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",gnd.V, gnd.A, gnd.H, gnd.X, gnd.Nc);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",giwn.V, giwn.A, giwn.H, giwn.X);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i\n",giwnc.V, giwnc.A, giwnc.H, giwnc.X, giwnc.Nc);
+		fflush(SUBPHASE);
 	}
 	else{
-		fprintf(TIMELOG,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",wc.p, wd.p, nc.p, nd.p);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Px, wd.Px, iwn.Px, nc.Px, nd.Px, iwn.Px);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Py, nc.Py, nd.Py, iwn.Py);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pz, nc.Pz, nd.Pz, iwn.Pz);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",wc.V, wc.A, wc.H, wc.X);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",nc.V, nc.A, nc.H, nc.X);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",nd.V, nd.A, nd.H, nd.X, nd.Nc);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",iwn.V, iwn.A, iwn.H, iwn.X);
-		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",iwnc.V, iwnc.A, iwnc.H, iwnc.X);
-		fflush(TIMELOG);
+		fprintf(SUBPHASE,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",wc.p, wd.p, nc.p, nd.p);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Px, wd.Px, iwn.Px, nc.Px, nd.Px, iwn.Px);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Py, nc.Py, nd.Py, iwn.Py);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pz, nc.Pz, nd.Pz, iwn.Pz);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",wc.V, wc.A, wc.H, wc.X);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",nc.V, nc.A, nc.H, nc.X);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",nd.V, nd.A, nd.H, nd.X, nd.Nc);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",iwn.V, iwn.A, iwn.H, iwn.X);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g\n",iwnc.V, iwnc.A, iwnc.H, iwnc.X);
+		fflush(SUBPHASE);
 	}
 
 }
diff --git a/analysis/SubPhase.h b/analysis/SubPhase.h
index b0a68347..738d2188 100644
--- a/analysis/SubPhase.h
+++ b/analysis/SubPhase.h
@@ -104,7 +104,7 @@ public:
 	
 private:
 	FILE *TIMELOG;
-
+	FILE *SUBPHASE;
 };
 
 #endif
diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp
index f86e2de8..07d3a13b 100644
--- a/analysis/morphology.cpp
+++ b/analysis/morphology.cpp
@@ -754,12 +754,20 @@ double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id,
 		if (morph_delta / morph_delta_previous > 2.0 ) morph_delta = morph_delta_previous*2.0;
 				
 		//MAX_DISPLACEMENT *= max(TargetGrowth/GrowthEstimate,1.25);
-		if (MAX_DISPLACEMENT > 1.0 ){
-			if (morph_delta > 0.0 ) morph_delta = 1.0;
-			else 					morph_delta = -1.0;
-			
-			//if (COUNT_FOR_LOOP > 2) COUNT_FOR_LOOP = 100;
-			COUNT_FOR_LOOP = 100; // exit loop if displacement is too large
+		
+		if (morph_delta > 0.0 ){
+			// object is growing
+			if (MAX_DISPLACEMENT > 3.0 ){
+				morph_delta = 3.0;
+				COUNT_FOR_LOOP = 100; // exit loop if displacement is too large
+			}
+		}
+		else{
+			// object is shrinking
+			if (MAX_DISPLACEMENT > 1.0 ){
+				morph_delta = -1.0;
+				COUNT_FOR_LOOP = 100; // exit loop if displacement is too large
+			}
 		}
 	}
 	if (rank == 0) printf("Final delta=%f \n",morph_delta);

From 7f9518e070e6c0ca6eba640f37a3b945f6179810 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 15:17:45 -0400
Subject: [PATCH 14/53] added check for header write

---
 analysis/SubPhase.cpp | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index 7b466a24..55db6e48 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -28,8 +28,15 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 
 	//.........................................
 	if (Dm->rank()==0){
+		bool WriteHeader=false;
+		SUBPHASE = fopen("subphase.csv","r");
+		if (SUBPHASE != NULL)
+			SUBPHASE.close();
+		else
+			WriteHeader=true;
+
 		SUBPHASE = fopen("subphase.csv","a+");
-		if (ftell(SUBPHASE) == 0)
+		if (WriteHeader)
 		{
 			// If timelog is empty, write a short header to list the averages
 			//fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n");

From 6b8ea6bdf0f78821520a0856bcd3585cd2fabddc Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 15:20:06 -0400
Subject: [PATCH 15/53] added header check

---
 analysis/SubPhase.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index 55db6e48..dbdcf33d 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -31,7 +31,7 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 		bool WriteHeader=false;
 		SUBPHASE = fopen("subphase.csv","r");
 		if (SUBPHASE != NULL)
-			SUBPHASE.close();
+			fclose(SUBPHASE);
 		else
 			WriteHeader=true;
 

From 2eddcfd8cd75d0d5afc6b6e920f7f26f9859afb7 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 15:30:34 -0400
Subject: [PATCH 16/53] increased default inlet/outlet layers to 10

---
 common/Domain.cpp | 74 +++++++++++++++++++++++++----------------------
 1 file changed, 39 insertions(+), 35 deletions(-)

diff --git a/common/Domain.cpp b/common/Domain.cpp
index 27a2359a..10dfbe14 100644
--- a/common/Domain.cpp
+++ b/common/Domain.cpp
@@ -514,43 +514,47 @@ void Domain::CommInit()
 
 void Domain::ReadIDs(){
 	// Read the IDs from input file
-  int nprocs=nprocx()*nprocy()*nprocz();
-    size_t readID;
-    char LocalRankString[8];
-    char LocalRankFilename[40];
-    //.......................................................................
-    if (rank() == 0)    printf("Read input media... \n");
-    //.......................................................................
-    sprintf(LocalRankString,"%05d",rank());
-    sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
-    // .......... READ THE INPUT FILE .......................................
-    if (rank()==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
-    sprintf(LocalRankFilename,"ID.%05i",rank());
-    FILE *IDFILE = fopen(LocalRankFilename,"rb");
-    if (IDFILE==NULL) ERROR("Domain::ReadIDs --  Error opening file: ID.xxxxx");
-    readID=fread(id,1,N,IDFILE);
-    if (readID != size_t(N)) printf("Domain::ReadIDs -- Error reading ID (rank=%i) \n",rank());
-    fclose(IDFILE);
-    // Compute the porosity
-    double sum;
-    double porosity;
-    double sum_local=0.0;
-    double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
-    if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
- 	//.........................................................
- 	// If external boundary conditions are applied remove solid
-    if (BoundaryCondition >  0  && kproc() == 0){
- 		for (int k=0; k<3; k++){
- 			for (int j=0;j<Ny;j++){
- 				for (int i=0;i<Nx;i++){
- 					int n = k*Nx*Ny+j*Nx+i;
- 					id[n] = 1;
- 				}                    
- 			}
+	int nprocs=nprocx()*nprocy()*nprocz();
+	size_t readID;
+	char LocalRankString[8];
+	char LocalRankFilename[40];
+	//.......................................................................
+	if (rank() == 0)    printf("Read input media... \n");
+	//.......................................................................
+	sprintf(LocalRankString,"%05d",rank());
+	sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
+	// .......... READ THE INPUT FILE .......................................
+	if (rank()==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
+	sprintf(LocalRankFilename,"ID.%05i",rank());
+	FILE *IDFILE = fopen(LocalRankFilename,"rb");
+	if (IDFILE==NULL) ERROR("Domain::ReadIDs --  Error opening file: ID.xxxxx");
+	readID=fread(id,1,N,IDFILE);
+	if (readID != size_t(N)) printf("Domain::ReadIDs -- Error reading ID (rank=%i) \n",rank());
+	fclose(IDFILE);
+	int inlet_layers=0;
+	int outlet_layers=0;
+	// Compute the porosity
+	double sum;
+	double porosity;
+	double sum_local=0.0;
+	double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
+	if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
+	//.........................................................
+	// If external boundary conditions are applied remove solid
+	if (BoundaryCondition >  0  && kproc() == 0){
+		inlet_layers=outlet_layers=10;
+		for (int k=0; k<inlet_layers; k++){
+			for (int j=0;j<Ny;j++){
+				for (int i=0;i<Nx;i++){
+					int n = k*Nx*Ny+j*Nx+i;
+					id[n] = 1;
+				}                    
+			}
  		}
  	}
     if (BoundaryCondition >  0  && kproc() == nprocz()-1){
- 		for (int k=Nz-3; k<Nz; k++){
+		outlet_layers=10;
+ 		for (int k=Nz-outlet_layers; k<Nz; k++){
  			for (int j=0;j<Ny;j++){
  				for (int i=0;i<Nx;i++){
  					int n = k*Nx*Ny+j*Nx+i;
@@ -559,7 +563,7 @@ void Domain::ReadIDs(){
  			}
  		}
  	}
-    for (int k=1;k<Nz-1;k++){
+    for (int k=inlet_layers+1; k<Nz-outlet_layers-1;k++){
         for (int j=1;j<Ny-1;j++){
             for (int i=1;i<Nx-1;i++){
                 int n = k*Nx*Ny+j*Nx+i;

From 79c3b0d1731153327c5c94ef30e1e0af0895ad27 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 15:42:26 -0400
Subject: [PATCH 17/53] added outlet layers to domain

---
 common/Domain.cpp | 30 ++++++++++++++++++++----------
 common/Domain.h   |  1 +
 2 files changed, 21 insertions(+), 10 deletions(-)

diff --git a/common/Domain.cpp b/common/Domain.cpp
index 10dfbe14..e88a6b25 100644
--- a/common/Domain.cpp
+++ b/common/Domain.cpp
@@ -77,6 +77,7 @@ Domain::Domain( std::shared_ptr<Database> db, MPI_Comm Communicator):
 	Lx(0), Ly(0), Lz(0), Volume(0), BoundaryCondition(0),
 	Comm(MPI_COMM_NULL),
 	inlet_layers_x(0), inlet_layers_y(0), inlet_layers_z(0),
+	outlet_layers_x(0), outlet_layers_y(0), outlet_layers_z(0),
 	sendCount_x(0), sendCount_y(0), sendCount_z(0), sendCount_X(0), sendCount_Y(0), sendCount_Z(0),
 	sendCount_xy(0), sendCount_yz(0), sendCount_xz(0), sendCount_Xy(0), sendCount_Yz(0), sendCount_xZ(0),
 	sendCount_xY(0), sendCount_yZ(0), sendCount_Xz(0), sendCount_XY(0), sendCount_YZ(0), sendCount_XZ(0),
@@ -126,6 +127,12 @@ void Domain::initialize( std::shared_ptr<Database> db )
 		inlet_layers_y = InletCount[1];
 		inlet_layers_z = InletCount[2];
 	}
+	if (d_db->keyExists( "OutletLayers" )){
+		auto OutletCount = d_db->getVector<int>( "OutletLayers" );
+		outlet_layers_x = OutletCount[0];
+		outlet_layers_y = OutletCount[1];
+		outlet_layers_z = OutletCount[2];
+	}
 	
     ASSERT( n.size() == 3u );
     ASSERT( L.size() == 3u );
@@ -144,9 +151,13 @@ void Domain::initialize( std::shared_ptr<Database> db )
     MPI_Comm_rank( Comm, &myrank );
 	rank_info = RankInfoStruct(myrank,nproc[0],nproc[1],nproc[2]);
 	// inlet layers only apply to lower part of domain
-	if (nproc[0] > 0) inlet_layers_x = 0;
-	if (nproc[1] > 0) inlet_layers_y = 0;
-	if (nproc[2] > 0) inlet_layers_z = 0;
+	if (rank_info.ix > 0) inlet_layers_x = 0;
+	if (rank_info.jy > 0) inlet_layers_y = 0;
+	if (rank_info.kz > 0) inlet_layers_z = 0;
+	// outlet layers only apply to top part of domain
+	if (rank_info.ix < nproc[0]-1 ) outlet_layers_x = 0;
+	if (rank_info.jy < nproc[1]-1) outlet_layers_y = 0;
+	if (rank_info.kz < nproc[2]-1) outlet_layers_z = 0;
     // Fill remaining variables
 	N = Nx*Ny*Nz;
 	Volume = nx*ny*nx*nproc[0]*nproc[1]*nproc[2]*1.0;
@@ -531,8 +542,7 @@ void Domain::ReadIDs(){
 	readID=fread(id,1,N,IDFILE);
 	if (readID != size_t(N)) printf("Domain::ReadIDs -- Error reading ID (rank=%i) \n",rank());
 	fclose(IDFILE);
-	int inlet_layers=0;
-	int outlet_layers=0;
+
 	// Compute the porosity
 	double sum;
 	double porosity;
@@ -542,8 +552,8 @@ void Domain::ReadIDs(){
 	//.........................................................
 	// If external boundary conditions are applied remove solid
 	if (BoundaryCondition >  0  && kproc() == 0){
-		inlet_layers=outlet_layers=10;
-		for (int k=0; k<inlet_layers; k++){
+    	if (inlet_layers_z < 4)	inlet_layers_z=4;
+		for (int k=0; k<inlet_layers_z; k++){
 			for (int j=0;j<Ny;j++){
 				for (int i=0;i<Nx;i++){
 					int n = k*Nx*Ny+j*Nx+i;
@@ -553,8 +563,8 @@ void Domain::ReadIDs(){
  		}
  	}
     if (BoundaryCondition >  0  && kproc() == nprocz()-1){
-		outlet_layers=10;
- 		for (int k=Nz-outlet_layers; k<Nz; k++){
+    	if (outlet_layers_z < 4)	outlet_layers_z=4;
+ 		for (int k=Nz-outlet_layers_z; k<Nz; k++){
  			for (int j=0;j<Ny;j++){
  				for (int i=0;i<Nx;i++){
  					int n = k*Nx*Ny+j*Nx+i;
@@ -563,7 +573,7 @@ void Domain::ReadIDs(){
  			}
  		}
  	}
-    for (int k=inlet_layers+1; k<Nz-outlet_layers-1;k++){
+    for (int k=inlet_layers_z+1; k<Nz-outlet_layers_z-1;k++){
         for (int j=1;j<Ny-1;j++){
             for (int i=1;i<Nx-1;i++){
                 int n = k*Nx*Ny+j*Nx+i;
diff --git a/common/Domain.h b/common/Domain.h
index 9fa87901..764ec957 100755
--- a/common/Domain.h
+++ b/common/Domain.h
@@ -110,6 +110,7 @@ public: // Public variables (need to create accessors instead)
     double Lx,Ly,Lz,Volume;
     int Nx,Ny,Nz,N;
     int inlet_layers_x, inlet_layers_y, inlet_layers_z;
+    int outlet_layers_x, outlet_layers_y, outlet_layers_z;
     RankInfoStruct rank_info;
 
     MPI_Comm Comm;        // MPI Communicator for this domain

From 9926cbc257c18dacdbdbb01fb20fa0819c7fd15a Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 15:53:11 -0400
Subject: [PATCH 18/53] added subphase analysis interval

---
 analysis/SubPhase.cpp    |  5 +++--
 analysis/runAnalysis.cpp | 42 +++++++++++++++++++++++++++++++++++++---
 analysis/runAnalysis.h   |  1 +
 3 files changed, 43 insertions(+), 5 deletions(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index dbdcf33d..f2597447 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -148,11 +148,12 @@ void SubPhase::Basic(){
 	if (Dm->BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
 
 	imin=jmin=1;
-	// If inlet layers exist use these as default
+	// If inlet/outlet layers exist use these as default
 	if (Dm->inlet_layers_x > 0) imin = Dm->inlet_layers_x;
 	if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y;
 	if (Dm->inlet_layers_z > 0) kmin = Dm->inlet_layers_z;
-		
+	if (Dm->outlet_layers_z > 0) kmax = Dm->outlet_layers_z;
+	
 	nb.reset(); wb.reset();
 /*
  	//Dm->CommunicateMeshHalo(Phi);
diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index 06f83171..abf87cb8 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -389,9 +389,9 @@ public:
             // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
         }
         if ( matches(type,AnalysisType::ComputeAverages) ) {
-            PROFILE_START("Compute subphase",1);
+            PROFILE_START("Compute basic averages",1);
             Averages.Basic();
-            PROFILE_STOP("Compute subphase",1);
+            PROFILE_STOP("Compute basic averages",1);
         }
     }
 private:
@@ -402,6 +402,29 @@ private:
     double beta;
 };
 
+class SubphaseWorkItem: public ThreadPool::WorkItemRet<void>
+{
+public:
+	SubphaseWorkItem( AnalysisType type_, int timestep_, SubPhase& Averages_ ):
+                type(type_), timestep(timestep_), Averages(Averages_){ }
+    ~SubphaseWorkItem() { }
+    virtual void run() {
+
+        if ( matches(type,AnalysisType::ComputeAverages) ) {
+            PROFILE_START("Compute subphase",1);
+            Averages.Full();
+            PROFILE_STOP("Compute subphase",1);
+        }
+    }
+private:
+    SubphaseWorkItem();
+    AnalysisType type;
+    int timestep;
+    SubPhase& Averages;
+    double beta;
+};
+
+
 
 /******************************************************************
  *  MPI comm wrapper for use with analysis                         *
@@ -479,6 +502,12 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
     d_blobid_interval = db->getScalar<int>( "blobid_interval" );
     d_visualization_interval = db->getScalar<int>( "visualization_interval" );
     auto restart_file = db->getScalar<std::string>( "restart_file" );
+	if (db->keyExists( "subphase_analysis_interval" ){
+		d_subphase_analysis_interval = db->getScalar<int>( "subphase_analysis_interval" );
+	}
+	else{
+		d_subphase_analysis_interval = INT_MAX;
+	}
     d_restartFile = restart_file + "." + rankString;
     d_rank = MPI_WORLD_RANK();
     writeIDMap(ID_map_struct(),0,id_map_filename);
@@ -892,9 +921,16 @@ void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, do
     // if ( matches(type,AnalysisType::ComputeAverages) ) {
     if ( timestep%d_analysis_interval == 0 ) {
         auto work = new BasicWorkItem(type,timestep,Averages);
-        work->add_dependency(d_wait_analysis);    // Make sure we are done using analysis before modifying
+        work->add_dependency(d_wait_subphase);    // Make sure we are done using analysis before modifying
+        work->add_dependency(d_wait_analysis);    
         d_wait_analysis = d_tpool.add_work(work);
     }
+    
+    if ( timestep%d_subphase_analysis_interval == 0 ) {
+        auto work = new BasicWorkItem(type,timestep,Averages);
+        work->add_dependency(d_wait_subphase);    // Make sure we are done using analysis before modifying
+        d_wait_subphase = d_tpool.add_work(work);
+    }
 
     if (timestep%d_restart_interval==0){
     	std::shared_ptr<double> cfq,cDen;
diff --git a/analysis/runAnalysis.h b/analysis/runAnalysis.h
index de2c1053..50e681eb 100644
--- a/analysis/runAnalysis.h
+++ b/analysis/runAnalysis.h
@@ -88,6 +88,7 @@ private:
     int d_Np;
     int d_rank;
     int d_restart_interval, d_analysis_interval, d_blobid_interval, d_visualization_interval;
+    int d_subphase_analysis_interval;
     double d_beta;
     bool d_regular;
     ThreadPool d_tpool;

From 3c664eadc3d603ceb7055431efb4a2e1fc78924a Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 15:54:00 -0400
Subject: [PATCH 19/53] added subphase analysis interval

---
 analysis/runAnalysis.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index abf87cb8..ae79df5d 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -491,6 +491,7 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
     ThreadPool::thread_id_t d_wait_analysis;
     ThreadPool::thread_id_t d_wait_vis;
     ThreadPool::thread_id_t d_wait_restart;
+    ThreadPool::thread_id_t d_wait_subphase;
 
     char rankString[20];
     sprintf(rankString,"%05d",Dm->rank());

From dcd79b4cebf53044d44271fa80f125e820f7d4ec Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 15:55:34 -0400
Subject: [PATCH 20/53] added subphase analysis interval

---
 analysis/runAnalysis.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/analysis/runAnalysis.h b/analysis/runAnalysis.h
index 50e681eb..788ce1b7 100644
--- a/analysis/runAnalysis.h
+++ b/analysis/runAnalysis.h
@@ -108,6 +108,7 @@ private:
     // Ids of work items to use for dependencies
     ThreadPool::thread_id_t d_wait_blobID;
     ThreadPool::thread_id_t d_wait_analysis;
+    ThreadPool::thread_id_t d_wait_subphase;
     ThreadPool::thread_id_t d_wait_vis;
     ThreadPool::thread_id_t d_wait_restart;
 

From 587f9bba62f79329dd8c2b4b476f46e9cf2f7021 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 15:56:22 -0400
Subject: [PATCH 21/53] added subphase analysis interval

---
 analysis/runAnalysis.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index ae79df5d..62802aa6 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -503,7 +503,7 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
     d_blobid_interval = db->getScalar<int>( "blobid_interval" );
     d_visualization_interval = db->getScalar<int>( "visualization_interval" );
     auto restart_file = db->getScalar<std::string>( "restart_file" );
-	if (db->keyExists( "subphase_analysis_interval" ){
+	if (db->keyExists( "subphase_analysis_interval" )){
 		d_subphase_analysis_interval = db->getScalar<int>( "subphase_analysis_interval" );
 	}
 	else{

From 281cd2c08c8e1cae58f42272abf2ffdea4d8e8b3 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 15:59:04 -0400
Subject: [PATCH 22/53] added c limits header

---
 analysis/runAnalysis.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/analysis/runAnalysis.h b/analysis/runAnalysis.h
index 788ce1b7..d9c96e6b 100644
--- a/analysis/runAnalysis.h
+++ b/analysis/runAnalysis.h
@@ -7,7 +7,7 @@
 #include "common/Communication.h"
 #include "common/ScaLBL.h"
 #include "threadpool/thread_pool.h"
-
+#include <limits.h>
 
 typedef std::shared_ptr<std::pair<int,IntArray>> BlobIDstruct;
 typedef std::shared_ptr<std::vector<BlobIDType>> BlobIDList;

From 790a18c59328cb643d45d9211340fcd4c1ab7339 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 16:05:11 -0400
Subject: [PATCH 23/53] update interface num components

---
 analysis/SubPhase.cpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index f2597447..0bf9fa51 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -457,7 +457,8 @@ void SubPhase::Full(){
 	giwnc.A=sumReduce( Dm->Comm, iwnc.A);
 	giwnc.H=sumReduce( Dm->Comm, iwnc.H);
 	giwnc.X=sumReduce( Dm->Comm, iwnc.X);
-	
+	giwnc.Nc = iwnc.Nc;
+
 	double vol_nc_bulk = 0.0;
 	double vol_wc_bulk = 0.0;
 	double vol_nd_bulk = 0.0;

From 28ef4fafc41124a23480429053b1b4639f1a2a35 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 16:57:18 -0400
Subject: [PATCH 24/53] fix bug in flow rate

---
 analysis/SubPhase.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index 0bf9fa51..468d15e4 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -223,7 +223,7 @@ void SubPhase::Basic(){
 		
 		double saturation=gwb.V/(gwb.V + gnb.V);
 		double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M;
-		double not_water_flow_rate=gnb.V*sqrt(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M;
+		double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M;
 		double total_flow_rate = water_flow_rate + not_water_flow_rate;
 		double fractional_flow= water_flow_rate / total_flow_rate;
 		printf("   water saturation = %f, fractional flow =%f \n",saturation,fractional_flow);

From c67fa8574fa33de47c5b9d9a85ce01468249ce9c Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 17:01:00 -0400
Subject: [PATCH 25/53] added steady interval

---
 models/ColorModel.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 19d0f9d0..c21275de 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -591,6 +591,7 @@ void ScaLBL_ColorModel::Run(){
 				if (fabs((Ca - Ca_previous)/Ca) < tolerance || CURRENT_STEADY_TIMESTEPS > MAX_STEADY_TIMESTEPS){
 					MORPH_ADAPT = true;
 					CURRENT_MORPH_TIMESTEPS=0;
+					CURRENT_STEADY_TIMESTEPS=0;
 					delta_volume_target = (volA + volB)*morph_delta; // set target volume change
 					Averages->Full();
 					Averages->Write(timestep);

From f51832cb5648918021fb89db79696d651ec7d85d Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 17:48:39 -0400
Subject: [PATCH 26/53] fix steady timestep print issue

---
 models/ColorModel.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index c21275de..6c480dee 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -591,7 +591,6 @@ void ScaLBL_ColorModel::Run(){
 				if (fabs((Ca - Ca_previous)/Ca) < tolerance || CURRENT_STEADY_TIMESTEPS > MAX_STEADY_TIMESTEPS){
 					MORPH_ADAPT = true;
 					CURRENT_MORPH_TIMESTEPS=0;
-					CURRENT_STEADY_TIMESTEPS=0;
 					delta_volume_target = (volA + volB)*morph_delta; // set target volume change
 					Averages->Full();
 					Averages->Write(timestep);
@@ -650,6 +649,7 @@ void ScaLBL_ColorModel::Run(){
 					}
 					morph_timesteps=0;
 				}
+				CURRENT_STEADY_TIMESTEPS=0;
 				Ca_previous = Ca;
 			}
 			if (MORPH_ADAPT ){

From 78f56bdb236c5d781c6c29e8e67b3cf14b7b4667 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 18:20:59 -0400
Subject: [PATCH 27/53] added vis to basic subphase

---
 analysis/runAnalysis.cpp | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index 62802aa6..11f2a615 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -952,6 +952,15 @@ void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, do
     	d_wait_restart = d_tpool.add_work(work1);
 
     }
+    
+    if (timestep%d_visualization_interval==0){
+        // Write the vis files
+        auto work = new WriteVisWorkItem( timestep, d_meshData, Averages, d_fillData, getComm() );
+        work->add_dependency(d_wait_analysis);
+        work->add_dependency(d_wait_subphase);
+        work->add_dependency(d_wait_vis);
+        d_wait_vis = d_tpool.add_work(work);
+    }
 
     PROFILE_STOP("run");
 }

From 3d038b232e54594cac3e1692784abe1bd7f379c8 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 18:24:31 -0400
Subject: [PATCH 28/53] added vis to basic subphase

---
 analysis/runAnalysis.cpp | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index 11f2a615..21dfea6f 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -498,17 +498,23 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
     d_N[0] = Dm->Nx;
     d_N[1] = Dm->Ny;
     d_N[2] = Dm->Nz;
+    
     d_restart_interval = db->getScalar<int>( "restart_interval" );
-    d_analysis_interval = db->getScalar<int>( "analysis_interval" );
-    d_blobid_interval = db->getScalar<int>( "blobid_interval" );
-    d_visualization_interval = db->getScalar<int>( "visualization_interval" );
-    auto restart_file = db->getScalar<std::string>( "restart_file" );
+    d_analysis_interval = db->getScalar<int>( "analysis_interval" );    
+    d_subphase_analysis_interval = INT_MAX;
+	d_visualization_interval = INT_MAX;
+	d_blobid_interval = INT_MAX;
+	if (db->keyExists( "blobid_interval" )){
+	    d_blobid_interval = db->getScalar<int>( "blobid_interval" );
+	}
+	if (db->keyExists( "visualization_interval" )){
+	    d_visualization_interval = db->getScalar<int>( "visualization_interval" );
+	}
 	if (db->keyExists( "subphase_analysis_interval" )){
 		d_subphase_analysis_interval = db->getScalar<int>( "subphase_analysis_interval" );
 	}
-	else{
-		d_subphase_analysis_interval = INT_MAX;
-	}
+	
+    auto restart_file = db->getScalar<std::string>( "restart_file" );
     d_restartFile = restart_file + "." + rankString;
     d_rank = MPI_WORLD_RANK();
     writeIDMap(ID_map_struct(),0,id_map_filename);

From 4f97bf78f242e44eb1225137fec06a3f0d02c247 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 18:26:49 -0400
Subject: [PATCH 29/53] added vis to basic subphase

---
 analysis/runAnalysis.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index 21dfea6f..36949ee0 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -961,7 +961,7 @@ void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, do
     
     if (timestep%d_visualization_interval==0){
         // Write the vis files
-        auto work = new WriteVisWorkItem( timestep, d_meshData, Averages, d_fillData, getComm() );
+        auto work = new IOWorkItem( timestep, d_meshData, Averages, d_fillData, getComm() );
         work->add_dependency(d_wait_analysis);
         work->add_dependency(d_wait_subphase);
         work->add_dependency(d_wait_vis);

From 6b16660e2d11e4db3f3972cb08efe221bed6e78b Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 21:35:14 -0400
Subject: [PATCH 30/53] switch checkers to allow void boundary

---
 tests/lbpm_serial_decomp.cpp | 27 +++++++++++++++------------
 1 file changed, 15 insertions(+), 12 deletions(-)

diff --git a/tests/lbpm_serial_decomp.cpp b/tests/lbpm_serial_decomp.cpp
index 738e387e..9c1c887e 100644
--- a/tests/lbpm_serial_decomp.cpp
+++ b/tests/lbpm_serial_decomp.cpp
@@ -87,6 +87,9 @@ int main(int argc, char **argv)
 	if (domain_db->keyExists( "checkerSize" )){
 		checkerSize = domain_db->getScalar<int>( "checkerSize" );
 	}
+	else {
+		checkerSize = SIZE[0];
+	}
 	auto ReadValues = domain_db->getVector<char>( "ReadValues" );
 	auto WriteValues = domain_db->getVector<char>( "WriteValues" );
 	auto ReadType = domain_db->getScalar<std::string>( "ReadType" );
@@ -158,13 +161,13 @@ int main(int argc, char **argv)
 			for (int j = 0; j<Ny; j++){
 				for (int i = xStart; i < xStart+inlet_count_x; i++){
 					if ( (j/checkerSize + k/checkerSize)%2 == 0){
-						// solid checkers
-						SegData[k*Nx*Ny+j*Nx+i] = 0;
-					}
-					else{
 						// void checkers
 						SegData[k*Nx*Ny+j*Nx+i] = 2;
 					}
+					else{
+						// solid checkers
+						SegData[k*Nx*Ny+j*Nx+i] = 0;
+					}
 				}
 			}
 		}
@@ -177,13 +180,13 @@ int main(int argc, char **argv)
 			for (int j = yStart; i < yStart+inlet_count_y; j++){
 				for (int i = 0; i<Nx; i++){
 					if ( (i/checkerSize + k/checkerSize)%2 == 0){
-						// solid checkers
-						SegData[k*Nx*Ny+j*Nx+i] = 0;
-					}
-					else{
 						// void checkers
 						SegData[k*Nx*Ny+j*Nx+i] = 2;
 					}
+					else{
+						// solid checkers
+						SegData[k*Nx*Ny+j*Nx+i] = 0;
+					}
 				}
 			}
 		}
@@ -196,13 +199,13 @@ int main(int argc, char **argv)
 			for (int j = 0; j<Ny; j++){
 				for (int i = 0; i<Nx; i++){
 					if ( (i/checkerSize+j/checkerSize)%2 == 0){
-						// solid checkers
-						SegData[k*Nx*Ny+j*Nx+i] = 0;
-					}
-					else{
 						// void checkers
 						SegData[k*Nx*Ny+j*Nx+i] = 2;
 					}
+					else{
+						// solid checkers
+						SegData[k*Nx*Ny+j*Nx+i] = 0;
+					}
 				}
 			}
 		}

From 5078343655f7fffadcadc8260bfcc372bbdbf6c5 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 27 Mar 2019 06:38:39 -0400
Subject: [PATCH 31/53] fix flow reversal

---
 models/ColorModel.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 6c480dee..a751bdeb 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -783,7 +783,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 		}
 	}
 
-	if (volume_connected < 0.02*volume_initial){
+	if (volume_connected < 0.02*volume_initial && target_delta_volume < 0.0){
 		// if connected volume is less than 2% just delete the whole thing
 		if (rank==0) printf("Connected region has shrunk to less than 2 %% of total fluid volume \n");
 		REVERSE_FLOW_DIRECTION = true;

From 21591f358975db170c7bf17a2ee7a13b1c6920f3 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 27 Mar 2019 07:03:56 -0400
Subject: [PATCH 32/53] basic logs to timelog.csv

---
 analysis/SubPhase.cpp | 58 +++++++++++++++++++++++++++++--------------
 1 file changed, 40 insertions(+), 18 deletions(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index 468d15e4..b69a043c 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -79,6 +79,23 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 		fprintf(SUBPHASE,"Vi Ai Hi Xi ");					 		// interface region 
 		fprintf(SUBPHASE,"Vic Aic Hic Xic Nic\n");					// interface region 
 	}
+	
+	if (Dm->rank()==0){
+		bool WriteHeader=false;
+		TIMELOG = fopen("timelog.csv","r");
+		if (TIMELOG != NULL)
+			fclose(TIMELOG);
+		else
+			WriteHeader=true;
+
+		TIMELOG = fopen("timelog.csv","a+");
+		if (WriteHeader)
+		{
+			// If timelog is empty, write a short header to list the averages
+			//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
+			fprintf(TIMELOG,"sw krw krn qw qn pw pn\n");				
+		}
+	}
 }
 
 
@@ -121,7 +138,6 @@ void SubPhase::Write(int timestep)
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",nd.V, nd.A, nd.H, nd.X, nd.Nc);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",iwn.V, iwn.A, iwn.H, iwn.X);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g\n",iwnc.V, iwnc.A, iwnc.H, iwnc.X);
-		fflush(SUBPHASE);
 	}
 
 }
@@ -155,22 +171,10 @@ void SubPhase::Basic(){
 	if (Dm->outlet_layers_z > 0) kmax = Dm->outlet_layers_z;
 	
 	nb.reset(); wb.reset();
-/*
- 	//Dm->CommunicateMeshHalo(Phi);
-	for (int k=1; k<Nz-1; k++){
-		for (int j=1; j<Ny-1; j++){
-			for (int i=1; i<Nx-1; i++){
-				// Compute all of the derivatives using finite differences
-				double fx = 0.5*(Phi(i+1,j,k) - Phi(i-1,j,k));
-				double fy = 0.5*(Phi(i,j+1,k) - Phi(i,j-1,k));
-				double fz = 0.5*(Phi(i,j,k+1) - Phi(i,j,k-1));
-				DelPhi(i,j,k) = sqrt(fx*fx+fy*fy+fz*fz);
-			}
-		}
-	}
-	*/
-	
+
 	double nA,nB;
+	double count_w = 0.0;
+	double count_n = 0.0;
 	for (k=kmin; k<kmax; k++){
 		for (j=jmin; j<Ny-1; j++){
 			for (i=imin; i<Nx-1; i++){
@@ -200,6 +204,14 @@ void SubPhase::Basic(){
 						wb.Py += rho_w*nB*Vel_y(n);
 						wb.Pz += rho_w*nB*Vel_z(n);
 					}
+					if ( phi > 0.99 ){
+						nb.P += Pressure(n);
+						count_n += 1.0;
+					}
+					else if ( phi < -0.99 ){
+						wb.P += Pressure(n);
+						count_w += 1.0;
+					}
 				}
 			}
 		}
@@ -214,19 +226,29 @@ void SubPhase::Basic(){
 	gnb.Px=sumReduce( Dm->Comm, nb.Px);
 	gnb.Py=sumReduce( Dm->Comm, nb.Py);
 	gnb.Pz=sumReduce( Dm->Comm, nb.Pz);
+	
+	count_w=sumReduce( Dm->Comm, count_w);
+	count_n=sumReduce( Dm->Comm, count_n);
+	gwb.p=sumReduce( Dm->Comm, wb.p) / count_w;
+	gnb.p=sumReduce( Dm->Comm, nb.p) / count_n;
 
 	if (Dm->rank() == 0){
 		double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
 		double dir_x = Fx/force_mag;
 		double dir_y = Fy/force_mag;
 		double dir_z = Fz/force_mag;
-		
+
 		double saturation=gwb.V/(gwb.V + gnb.V);
 		double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M;
 		double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M;
 		double total_flow_rate = water_flow_rate + not_water_flow_rate;
 		double fractional_flow= water_flow_rate / total_flow_rate;
-		printf("   water saturation = %f, fractional flow =%f \n",saturation,fractional_flow);
+		
+		double krn = nu_n*not_water_flow_rate / force_mag;
+		double krw = nu_w*water_flow_rate / force_mag;
+		//printf("   water saturation = %f, fractional flow =%f \n",saturation,fractional_flow);
+		fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,water_flow_rate,not_water_flow_rate, gwb.p, gnb.p); 
+		fflush(TIMELOG);
 	}
 
 }

From 19492ec4b4e6657bec0b6851b6db8aa8014c3aab Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 27 Mar 2019 07:10:43 -0400
Subject: [PATCH 33/53] subphase enabled in basic

---
 analysis/SubPhase.cpp    | 4 ++--
 analysis/runAnalysis.cpp | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index b69a043c..2cf2a206 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -205,11 +205,11 @@ void SubPhase::Basic(){
 						wb.Pz += rho_w*nB*Vel_z(n);
 					}
 					if ( phi > 0.99 ){
-						nb.P += Pressure(n);
+						nb.p += Pressure(n);
 						count_n += 1.0;
 					}
 					else if ( phi < -0.99 ){
-						wb.P += Pressure(n);
+						wb.p += Pressure(n);
 						count_w += 1.0;
 					}
 				}
diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index 36949ee0..2e357e11 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -934,7 +934,7 @@ void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, do
     }
     
     if ( timestep%d_subphase_analysis_interval == 0 ) {
-        auto work = new BasicWorkItem(type,timestep,Averages);
+        auto work = new SubphaseWorkItem(type,timestep,Averages);
         work->add_dependency(d_wait_subphase);    // Make sure we are done using analysis before modifying
         d_wait_subphase = d_tpool.add_work(work);
     }

From 231401aacacf5167cf1d98799a7f75d643fd7cd8 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 27 Mar 2019 07:25:57 -0400
Subject: [PATCH 34/53] default to z direction for timelog.csv

---
 analysis/SubPhase.cpp | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index 2cf2a206..2100a745 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -237,7 +237,11 @@ void SubPhase::Basic(){
 		double dir_x = Fx/force_mag;
 		double dir_y = Fy/force_mag;
 		double dir_z = Fz/force_mag;
-
+		if (force_mag == 0.0){
+			// default to z direction
+			dir_z = 1.0;
+			force_mag = 1.0;
+		}
 		double saturation=gwb.V/(gwb.V + gnb.V);
 		double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M;
 		double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M;

From 0b802f4be380e0fdaf8c8a786a004a8303755ed0 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 27 Mar 2019 07:33:42 -0400
Subject: [PATCH 35/53] default to z direction for timelog.csv

---
 analysis/SubPhase.cpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index 2100a745..2a7e61b7 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -239,6 +239,8 @@ void SubPhase::Basic(){
 		double dir_z = Fz/force_mag;
 		if (force_mag == 0.0){
 			// default to z direction
+			dir_x = 0.0;
+			dir_y = 0.0;
 			dir_z = 1.0;
 			force_mag = 1.0;
 		}

From e02eb967251391532657127752b80576adcf1138 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 27 Mar 2019 07:37:55 -0400
Subject: [PATCH 36/53] add dependencies for subphase analysis

---
 analysis/runAnalysis.cpp | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index 2e357e11..7e0e77d7 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -929,13 +929,16 @@ void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, do
     if ( timestep%d_analysis_interval == 0 ) {
         auto work = new BasicWorkItem(type,timestep,Averages);
         work->add_dependency(d_wait_subphase);    // Make sure we are done using analysis before modifying
-        work->add_dependency(d_wait_analysis);    
+        work->add_dependency(d_wait_analysis);  
+        work->add_dependency(d_wait_vis);
         d_wait_analysis = d_tpool.add_work(work);
     }
     
     if ( timestep%d_subphase_analysis_interval == 0 ) {
         auto work = new SubphaseWorkItem(type,timestep,Averages);
         work->add_dependency(d_wait_subphase);    // Make sure we are done using analysis before modifying
+        work->add_dependency(d_wait_analysis);  
+        work->add_dependency(d_wait_vis);
         d_wait_subphase = d_tpool.add_work(work);
     }
 

From c0b3b04172bd014ed03aa13f2f15fabf923932de Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 27 Mar 2019 07:39:10 -0400
Subject: [PATCH 37/53] add dependencies for subphase analysis

---
 analysis/runAnalysis.cpp | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index 7e0e77d7..a98a6448 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -410,11 +410,9 @@ public:
     ~SubphaseWorkItem() { }
     virtual void run() {
 
-        if ( matches(type,AnalysisType::ComputeAverages) ) {
-            PROFILE_START("Compute subphase",1);
-            Averages.Full();
-            PROFILE_STOP("Compute subphase",1);
-        }
+    	PROFILE_START("Compute subphase",1);
+    	Averages.Full();
+    	PROFILE_STOP("Compute subphase",1);
     }
 private:
     SubphaseWorkItem();

From bf55f928b854775f1970178ed77c9c0c47652712 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 27 Mar 2019 12:07:35 -0400
Subject: [PATCH 38/53] subphase not quite working with threadpool

---
 analysis/runAnalysis.cpp | 1 +
 analysis/runAnalysis.h   | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index a98a6448..eef31fd8 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -598,6 +598,7 @@ void runAnalysis::finish( )
     d_wait_blobID.reset();
     d_wait_analysis.reset();
     d_wait_vis.reset();
+    d_wait_subphase.reset();
     d_wait_restart.reset();
     // Syncronize
     MPI_Barrier( d_comm );
diff --git a/analysis/runAnalysis.h b/analysis/runAnalysis.h
index d9c96e6b..df73fafd 100644
--- a/analysis/runAnalysis.h
+++ b/analysis/runAnalysis.h
@@ -15,7 +15,7 @@ typedef std::shared_ptr<std::vector<BlobIDType>> BlobIDList;
 
 // Types of analysis
 enum class AnalysisType : uint64_t { AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02, 
-    CopySimState=0x04, ComputeAverages=0x08, CreateRestart=0x10, WriteVis=0x20 };
+    CopySimState=0x04, ComputeAverages=0x08, CreateRestart=0x10, WriteVis=0x20, ComputeSubphase=0x40 };
 
 
 //! Class to run the analysis in multiple threads

From a46c9104c3ed197671cb9a57c9549f2093758765 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 27 Mar 2019 12:59:56 -0400
Subject: [PATCH 39/53] write subphase -- dumb mistake

---
 analysis/runAnalysis.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index eef31fd8..a2b3801d 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -412,6 +412,7 @@ public:
 
     	PROFILE_START("Compute subphase",1);
     	Averages.Full();
+    	Averages.Write(timestep);
     	PROFILE_STOP("Compute subphase",1);
     }
 private:

From e712090f72afba8f66a9705f376eec498016520b Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Thu, 28 Mar 2019 15:46:52 -0400
Subject: [PATCH 40/53] read labels as signed char

---
 models/ColorModel.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index a751bdeb..5552e9bd 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -163,7 +163,7 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
 	char VALUE=0;
 	double AFFINITY=0.f;
 
-	auto LabelList = color_db->getVector<char>( "ComponentLabels" );
+	auto LabelList = color_db->getVector<signed char>( "ComponentLabels" );
 	auto AffinityList = color_db->getVector<double>( "ComponentAffinity" );
 
 	NLABELS=LabelList.size();

From d2cc5086d7c89cd6a05207ae44e668b1c9e69882 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Thu, 28 Mar 2019 17:52:00 -0400
Subject: [PATCH 41/53] revert signed char

---
 models/ColorModel.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 5552e9bd..a751bdeb 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -163,7 +163,7 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
 	char VALUE=0;
 	double AFFINITY=0.f;
 
-	auto LabelList = color_db->getVector<signed char>( "ComponentLabels" );
+	auto LabelList = color_db->getVector<char>( "ComponentLabels" );
 	auto AffinityList = color_db->getVector<double>( "ComponentAffinity" );
 
 	NLABELS=LabelList.size();

From a8bd5dca7e88dcb0cfce8291cf47994dfcf0e1c3 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Fri, 29 Mar 2019 05:23:49 -0400
Subject: [PATCH 42/53] trying to make label more robust

---
 models/ColorModel.cpp        |  6 +++---
 tests/lbpm_serial_decomp.cpp | 10 +++++-----
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 5552e9bd..a937f403 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -160,10 +160,10 @@ void ScaLBL_ColorModel::ReadInput(){
 void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
 {
 	size_t NLABELS=0;
-	char VALUE=0;
+	signed char VALUE=0;
 	double AFFINITY=0.f;
 
-	auto LabelList = color_db->getVector<signed char>( "ComponentLabels" );
+	auto LabelList = color_db->getVector<int>( "ComponentLabels" );
 	auto AffinityList = color_db->getVector<double>( "ComponentAffinity" );
 
 	NLABELS=LabelList.size();
@@ -210,7 +210,7 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
 			VALUE=LabelList[idx];
 			AFFINITY=AffinityList[idx];
 			double volume_fraction  = double(label_count_global[idx])/double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
-			printf("   label=%hhd, affinity=%f, volume fraction==%f\n",VALUE,AFFINITY,volume_fraction); 
+			printf("   label=%d, affinity=%f, volume fraction==%f\n",VALUE,AFFINITY,volume_fraction); 
 		}
 	}
 
diff --git a/tests/lbpm_serial_decomp.cpp b/tests/lbpm_serial_decomp.cpp
index 9c1c887e..9277797f 100644
--- a/tests/lbpm_serial_decomp.cpp
+++ b/tests/lbpm_serial_decomp.cpp
@@ -90,8 +90,8 @@ int main(int argc, char **argv)
 	else {
 		checkerSize = SIZE[0];
 	}
-	auto ReadValues = domain_db->getVector<char>( "ReadValues" );
-	auto WriteValues = domain_db->getVector<char>( "WriteValues" );
+	auto ReadValues = domain_db->getVector<int>( "ReadValues" );
+	auto WriteValues = domain_db->getVector<int>( "WriteValues" );
 	auto ReadType = domain_db->getScalar<std::string>( "ReadType" );
 	if (ReadType == "8bit"){
 	}
@@ -262,8 +262,8 @@ int main(int argc, char **argv)
 								n = k*(nx+2)*(ny+2) + j*(nx+2) + i;;
 								char locval = loc_id[n];
 								for (int idx=0; idx<ReadValues.size(); idx++){
-									char oldvalue=ReadValues[idx];
-									char newvalue=WriteValues[idx];
+									signed char oldvalue=ReadValues[idx];
+									signed char newvalue=WriteValues[idx];
 									if (locval == oldvalue){
 										loc_id[n] = newvalue;
 										LabelCount[idx]++;
@@ -288,7 +288,7 @@ int main(int argc, char **argv)
 		}
 	}
 	for (int idx=0; idx<ReadValues.size(); idx++){
-		char label=ReadValues[idx];
+		int label=ReadValues[idx];
 		int count=LabelCount[idx];
 		printf("Label=%d, Count=%d \n",label,count);
 	}

From 4e33d1c412f7b856d84cdad26eea5099e084d787 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Fri, 29 Mar 2019 06:04:00 -0400
Subject: [PATCH 43/53] initialize CPU phi with labels

---
 models/ColorModel.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index a937f403..ed5b355a 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -403,7 +403,7 @@ void ScaLBL_ColorModel::Initialize(){
 			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
 		}
 	}
-
+	ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double));
 }
 
 void ScaLBL_ColorModel::Run(){

From 2da6662e9971f1fab28850cc7312cd75235ec8a1 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Fri, 29 Mar 2019 06:07:31 -0400
Subject: [PATCH 44/53] initialize CPU phi with labels

---
 models/ColorModel.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index ed5b355a..93e44ae9 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -403,7 +403,7 @@ void ScaLBL_ColorModel::Initialize(){
 			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
 		}
 	}
-	ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double));
+	ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double));
 }
 
 void ScaLBL_ColorModel::Run(){

From a0fabd4b32ed795ab1a5049197f6005097f301d5 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Fri, 29 Mar 2019 06:46:48 -0400
Subject: [PATCH 45/53] make signed char explicit

---
 common/Domain.cpp            | 2 +-
 common/Domain.h              | 6 +++---
 common/ScaLBL.cpp            | 2 +-
 common/ScaLBL.h              | 2 +-
 models/ColorModel.h          | 2 +-
 tests/lbpm_serial_decomp.cpp | 4 ++--
 6 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/common/Domain.cpp b/common/Domain.cpp
index e88a6b25..931ac5dd 100644
--- a/common/Domain.cpp
+++ b/common/Domain.cpp
@@ -161,7 +161,7 @@ void Domain::initialize( std::shared_ptr<Database> db )
     // Fill remaining variables
 	N = Nx*Ny*Nz;
 	Volume = nx*ny*nx*nproc[0]*nproc[1]*nproc[2]*1.0;
-	id = new char[N];
+	id = new signed char[N];
 	memset(id,0,N);
 	BoundaryCondition = d_db->getScalar<int>("BC");
     int nprocs;
diff --git a/common/Domain.h b/common/Domain.h
index 764ec957..038627c2 100755
--- a/common/Domain.h
+++ b/common/Domain.h
@@ -170,7 +170,7 @@ public: // Public variables (need to create accessors instead)
     int *recvList_xY, *recvList_yZ, *recvList_Xz, *recvList_XY, *recvList_YZ, *recvList_XZ;
     //......................................................................................    
     // Solid indicator function
-    char *id;
+    signed char *id;
 
     void ReadIDs();
     void CommunicateMeshHalo(DoubleArray &Mesh);
@@ -179,8 +179,8 @@ public: // Public variables (need to create accessors instead)
 
 private:
 
-    void PackID(int *list, int count, char *sendbuf, char *ID);
-    void UnpackID(int *list, int count, char *recvbuf, char *ID);
+    void PackID(int *list, int count, signed char *sendbuf, signed char *ID);
+    void UnpackID(int *list, int count, signed char *recvbuf, signed char *ID);
     void CommHaloIDs();
     
 	//......................................................................................
diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp
index 0c4caa81..21656757 100644
--- a/common/ScaLBL.cpp
+++ b/common/ScaLBL.cpp
@@ -352,7 +352,7 @@ void ScaLBL_Communicator::D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, int *list,  i
 	delete [] ReturnDist;
 }
 
-int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, char *id, int Np){
+int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np){
 	/*
 	 * Generate a memory optimized layout
 	 *   id[n] == 0 implies that site n should be ignored (treat as a mask)
diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index f1eee8fe..a50ab7ed 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -158,7 +158,7 @@ public:
 	int FirstInterior();
 	int LastInterior();
 	
-	int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, char *id, int Np);
+	int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np);
 //	void MemoryOptimizedLayout(IntArray &Map, int *neighborList, char *id, int Np);
 //	void MemoryOptimizedLayoutFull(IntArray &Map, int *neighborList, char *id, int Np);
 //	void MemoryDenseLayout(IntArray &Map, int *neighborList, char *id, int Np);
diff --git a/models/ColorModel.h b/models/ColorModel.h
index 805d8e22..3913270f 100644
--- a/models/ColorModel.h
+++ b/models/ColorModel.h
@@ -57,7 +57,7 @@ public:
     std::shared_ptr<Database> analysis_db;
 
     IntArray Map;
-    char *id;    
+    signed char *id;    
 	int *NeighborList;
 	int *dvcMap;
 	double *fq, *Aq, *Bq;
diff --git a/tests/lbpm_serial_decomp.cpp b/tests/lbpm_serial_decomp.cpp
index 9277797f..41cb9686 100644
--- a/tests/lbpm_serial_decomp.cpp
+++ b/tests/lbpm_serial_decomp.cpp
@@ -115,8 +115,8 @@ int main(int argc, char **argv)
 	printf("Input media: %s\n",Filename.c_str());
 	printf("Relabeling %lu values\n",ReadValues.size());
 	for (int idx=0; idx<ReadValues.size(); idx++){
-		char oldvalue=ReadValues[idx];
-		char newvalue=WriteValues[idx];
+		int oldvalue=ReadValues[idx];
+		int newvalue=WriteValues[idx];
 		printf("oldvalue=%d, newvalue =%d \n",oldvalue,newvalue);
 	}
 

From 48af5bc1383ffe380ab4aa4533181a95d7621ab6 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Fri, 29 Mar 2019 06:54:17 -0400
Subject: [PATCH 46/53] deprecated redundant functionality

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

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index c6775e68..37561fbe 100755
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -12,9 +12,9 @@ ADD_LBPM_EXECUTABLE( lbpm_refine_pp )
 ADD_LBPM_EXECUTABLE( lbpm_morphdrain_pp )
 ADD_LBPM_EXECUTABLE( lbpm_morphopen_pp )
 #ADD_LBPM_EXECUTABLE( lbpm_morph_pp )
-ADD_LBPM_EXECUTABLE( lbpm_segmented_pp )
+#ADD_LBPM_EXECUTABLE( lbpm_segmented_pp )
 #ADD_LBPM_EXECUTABLE( lbpm_block_pp )
-ADD_LBPM_EXECUTABLE( lbpm_segmented_decomp )
+#ADD_LBPM_EXECUTABLE( lbpm_segmented_decomp )
 ADD_LBPM_EXECUTABLE( lbpm_serial_decomp )
 #ADD_LBPM_EXECUTABLE( lbpm_disc_pp )
 #ADD_LBPM_EXECUTABLE( lbpm_juanes_bench_disc_pp )

From 0405e285ac8fdaec771dccb5a7258ff81114a602 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Fri, 29 Mar 2019 06:56:24 -0400
Subject: [PATCH 47/53] explicit sign char

---
 models/ColorModel.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 93e44ae9..adadc89d 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -104,7 +104,7 @@ void ScaLBL_ColorModel::SetDomain(){
 	Mask  = std::shared_ptr<Domain>(new Domain(domain_db,comm));    // mask domain removes immobile phases
 	Nx+=2; Ny+=2; Nz += 2;
 	N = Nx*Ny*Nz;
-	id = new char [N];
+	id = new signed char [N];
 	for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1;               // initialize this way
 	//Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
 	Averages = std::shared_ptr<SubPhase> ( new SubPhase(Dm) ); // TwoPhase analysis object

From 2b36a4257d5c7f342a98f26aa10db428f283bee3 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Fri, 29 Mar 2019 07:26:07 -0400
Subject: [PATCH 48/53] explict signed char in morph tools

---
 analysis/morphology.cpp      | 176 +++++++++++++++++------------------
 tests/lbpm_morphdrain_pp.cpp |  22 ++---
 tests/lbpm_morphopen_pp.cpp  |  20 ++--
 3 files changed, 109 insertions(+), 109 deletions(-)

diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp
index 07d3a13b..57ead57e 100644
--- a/analysis/morphology.cpp
+++ b/analysis/morphology.cpp
@@ -1,7 +1,7 @@
 #include <analysis/morphology.h>
 // Implementation of morphological opening routine
 
-inline void PackID(int *list, int count, char *sendbuf, char *ID){
+inline void PackID(int *list, int count, signed char *sendbuf, signed char *ID){
 	// Fill in the phase ID values from neighboring processors
 	// This packs up the values that need to be sent from one processor to another
 	int idx,n;
@@ -13,7 +13,7 @@ inline void PackID(int *list, int count, char *sendbuf, char *ID){
 }
 //***************************************************************************************
 
-inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
+inline void UnpackID(int *list, int count, signed char *recvbuf, signed char *ID){
 	// Fill in the phase ID values from neighboring processors
 	// This unpacks the values once they have been recieved from neighbors
 	int idx,n;
@@ -25,7 +25,7 @@ inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
 }
 
 //***************************************************************************************
-double MorphOpen(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, double VoidFraction){
+double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction){
 	// SignDist is the distance to the object that you want to constaing the morphological opening
 	// VoidFraction is the the empty space where the object inst
 	// id is a labeled map
@@ -73,51 +73,51 @@ double MorphOpen(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, do
 	if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal);
 
 	// Communication buffers
-	char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
-	char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
-	char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
-	char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z;
-	char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ;
-	char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ;
+	signed char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
+	signed char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
+	signed char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
+	signed char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z;
+	signed char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ;
+	signed char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ;
 	// send buffers
-	sendID_x = new char [Dm->sendCount_x];
-	sendID_y = new char [Dm->sendCount_y];
-	sendID_z = new char [Dm->sendCount_z];
-	sendID_X = new char [Dm->sendCount_X];
-	sendID_Y = new char [Dm->sendCount_Y];
-	sendID_Z = new char [Dm->sendCount_Z];
-	sendID_xy = new char [Dm->sendCount_xy];
-	sendID_yz = new char [Dm->sendCount_yz];
-	sendID_xz = new char [Dm->sendCount_xz];
-	sendID_Xy = new char [Dm->sendCount_Xy];
-	sendID_Yz = new char [Dm->sendCount_Yz];
-	sendID_xZ = new char [Dm->sendCount_xZ];
-	sendID_xY = new char [Dm->sendCount_xY];
-	sendID_yZ = new char [Dm->sendCount_yZ];
-	sendID_Xz = new char [Dm->sendCount_Xz];
-	sendID_XY = new char [Dm->sendCount_XY];
-	sendID_YZ = new char [Dm->sendCount_YZ];
-	sendID_XZ = new char [Dm->sendCount_XZ];
+	sendID_x = new signed char [Dm->sendCount_x];
+	sendID_y = new signed char [Dm->sendCount_y];
+	sendID_z = new signed char [Dm->sendCount_z];
+	sendID_X = new signed char [Dm->sendCount_X];
+	sendID_Y = new signed char [Dm->sendCount_Y];
+	sendID_Z = new signed char [Dm->sendCount_Z];
+	sendID_xy = new signed char [Dm->sendCount_xy];
+	sendID_yz = new signed char [Dm->sendCount_yz];
+	sendID_xz = new signed char [Dm->sendCount_xz];
+	sendID_Xy = new signed char [Dm->sendCount_Xy];
+	sendID_Yz = new signed char [Dm->sendCount_Yz];
+	sendID_xZ = new signed char [Dm->sendCount_xZ];
+	sendID_xY = new signed char [Dm->sendCount_xY];
+	sendID_yZ = new signed char [Dm->sendCount_yZ];
+	sendID_Xz = new signed char [Dm->sendCount_Xz];
+	sendID_XY = new signed char [Dm->sendCount_XY];
+	sendID_YZ = new signed char [Dm->sendCount_YZ];
+	sendID_XZ = new signed char [Dm->sendCount_XZ];
 	//......................................................................................
 	// recv buffers
-	recvID_x = new char [Dm->recvCount_x];
-	recvID_y = new char [Dm->recvCount_y];
-	recvID_z = new char [Dm->recvCount_z];
-	recvID_X = new char [Dm->recvCount_X];
-	recvID_Y = new char [Dm->recvCount_Y];
-	recvID_Z = new char [Dm->recvCount_Z];
-	recvID_xy = new char [Dm->recvCount_xy];
-	recvID_yz = new char [Dm->recvCount_yz];
-	recvID_xz = new char [Dm->recvCount_xz];
-	recvID_Xy = new char [Dm->recvCount_Xy];
-	recvID_xZ = new char [Dm->recvCount_xZ];
-	recvID_xY = new char [Dm->recvCount_xY];
-	recvID_yZ = new char [Dm->recvCount_yZ];
-	recvID_Yz = new char [Dm->recvCount_Yz];
-	recvID_Xz = new char [Dm->recvCount_Xz];
-	recvID_XY = new char [Dm->recvCount_XY];
-	recvID_YZ = new char [Dm->recvCount_YZ];
-	recvID_XZ = new char [Dm->recvCount_XZ];
+	recvID_x = new signed char [Dm->recvCount_x];
+	recvID_y = new signed char [Dm->recvCount_y];
+	recvID_z = new signed char [Dm->recvCount_z];
+	recvID_X = new signed char [Dm->recvCount_X];
+	recvID_Y = new signed char [Dm->recvCount_Y];
+	recvID_Z = new signed char [Dm->recvCount_Z];
+	recvID_xy = new signed char [Dm->recvCount_xy];
+	recvID_yz = new signed char [Dm->recvCount_yz];
+	recvID_xz = new signed char [Dm->recvCount_xz];
+	recvID_Xy = new signed char [Dm->recvCount_Xy];
+	recvID_xZ = new signed char [Dm->recvCount_xZ];
+	recvID_xY = new signed char [Dm->recvCount_xY];
+	recvID_yZ = new signed char [Dm->recvCount_yZ];
+	recvID_Yz = new signed char [Dm->recvCount_Yz];
+	recvID_Xz = new signed char [Dm->recvCount_Xz];
+	recvID_XY = new signed char [Dm->recvCount_XY];
+	recvID_YZ = new signed char [Dm->recvCount_YZ];
+	recvID_XZ = new signed char [Dm->recvCount_XZ];
 	//......................................................................................
 	int sendtag,recvtag;
 	sendtag = recvtag = 7;
@@ -327,7 +327,7 @@ double morph_open()
 */
 
 //***************************************************************************************
-double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, double VoidFraction){
+double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction){
 	// SignDist is the distance to the object that you want to constaing the morphological opening
 	// VoidFraction is the the empty space where the object inst
 	// id is a labeled map
@@ -378,51 +378,51 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, d
 	if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal);
 
 	// Communication buffers
-	char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
-	char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
-	char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
-	char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z;
-	char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ;
-	char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ;
+	signed char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
+	signed char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
+	signed char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
+	signed char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z;
+	signed char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ;
+	signed char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ;
 	// send buffers
-	sendID_x = new char [Dm->sendCount_x];
-	sendID_y = new char [Dm->sendCount_y];
-	sendID_z = new char [Dm->sendCount_z];
-	sendID_X = new char [Dm->sendCount_X];
-	sendID_Y = new char [Dm->sendCount_Y];
-	sendID_Z = new char [Dm->sendCount_Z];
-	sendID_xy = new char [Dm->sendCount_xy];
-	sendID_yz = new char [Dm->sendCount_yz];
-	sendID_xz = new char [Dm->sendCount_xz];
-	sendID_Xy = new char [Dm->sendCount_Xy];
-	sendID_Yz = new char [Dm->sendCount_Yz];
-	sendID_xZ = new char [Dm->sendCount_xZ];
-	sendID_xY = new char [Dm->sendCount_xY];
-	sendID_yZ = new char [Dm->sendCount_yZ];
-	sendID_Xz = new char [Dm->sendCount_Xz];
-	sendID_XY = new char [Dm->sendCount_XY];
-	sendID_YZ = new char [Dm->sendCount_YZ];
-	sendID_XZ = new char [Dm->sendCount_XZ];
+	sendID_x = new signed char [Dm->sendCount_x];
+	sendID_y = new signed char [Dm->sendCount_y];
+	sendID_z = new signed char [Dm->sendCount_z];
+	sendID_X = new signed char [Dm->sendCount_X];
+	sendID_Y = new signed char [Dm->sendCount_Y];
+	sendID_Z = new signed char [Dm->sendCount_Z];
+	sendID_xy = new signed char [Dm->sendCount_xy];
+	sendID_yz = new signed char [Dm->sendCount_yz];
+	sendID_xz = new signed char [Dm->sendCount_xz];
+	sendID_Xy = new signed char [Dm->sendCount_Xy];
+	sendID_Yz = new signed char [Dm->sendCount_Yz];
+	sendID_xZ = new signed char [Dm->sendCount_xZ];
+	sendID_xY = new signed char [Dm->sendCount_xY];
+	sendID_yZ = new signed char [Dm->sendCount_yZ];
+	sendID_Xz = new signed char [Dm->sendCount_Xz];
+	sendID_XY = new signed char [Dm->sendCount_XY];
+	sendID_YZ = new signed char [Dm->sendCount_YZ];
+	sendID_XZ = new signed char [Dm->sendCount_XZ];
 	//......................................................................................
 	// recv buffers
-	recvID_x = new char [Dm->recvCount_x];
-	recvID_y = new char [Dm->recvCount_y];
-	recvID_z = new char [Dm->recvCount_z];
-	recvID_X = new char [Dm->recvCount_X];
-	recvID_Y = new char [Dm->recvCount_Y];
-	recvID_Z = new char [Dm->recvCount_Z];
-	recvID_xy = new char [Dm->recvCount_xy];
-	recvID_yz = new char [Dm->recvCount_yz];
-	recvID_xz = new char [Dm->recvCount_xz];
-	recvID_Xy = new char [Dm->recvCount_Xy];
-	recvID_xZ = new char [Dm->recvCount_xZ];
-	recvID_xY = new char [Dm->recvCount_xY];
-	recvID_yZ = new char [Dm->recvCount_yZ];
-	recvID_Yz = new char [Dm->recvCount_Yz];
-	recvID_Xz = new char [Dm->recvCount_Xz];
-	recvID_XY = new char [Dm->recvCount_XY];
-	recvID_YZ = new char [Dm->recvCount_YZ];
-	recvID_XZ = new char [Dm->recvCount_XZ];
+	recvID_x = new signed char [Dm->recvCount_x];
+	recvID_y = new signed char [Dm->recvCount_y];
+	recvID_z = new signed char [Dm->recvCount_z];
+	recvID_X = new signed char [Dm->recvCount_X];
+	recvID_Y = new signed char [Dm->recvCount_Y];
+	recvID_Z = new signed char [Dm->recvCount_Z];
+	recvID_xy = new signed char [Dm->recvCount_xy];
+	recvID_yz = new signed char [Dm->recvCount_yz];
+	recvID_xz = new signed char [Dm->recvCount_xz];
+	recvID_Xy = new signed char [Dm->recvCount_Xy];
+	recvID_xZ = new signed char [Dm->recvCount_xZ];
+	recvID_xY = new signed char [Dm->recvCount_xY];
+	recvID_yZ = new signed char [Dm->recvCount_yZ];
+	recvID_Yz = new signed char [Dm->recvCount_Yz];
+	recvID_Xz = new signed char [Dm->recvCount_Xz];
+	recvID_XY = new signed char [Dm->recvCount_XY];
+	recvID_YZ = new signed char [Dm->recvCount_YZ];
+	recvID_XZ = new signed char [Dm->recvCount_XZ];
 	//......................................................................................
 	int sendtag,recvtag;
 	sendtag = recvtag = 7;
diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp
index c9a58222..675a76de 100644
--- a/tests/lbpm_morphdrain_pp.cpp
+++ b/tests/lbpm_morphdrain_pp.cpp
@@ -59,8 +59,8 @@ int main(int argc, char **argv)
 		auto L = domain_db->getVector<double>( "L" );
 		auto size = domain_db->getVector<int>( "n" );
 		auto nproc = domain_db->getVector<int>( "nproc" );
-		auto ReadValues = domain_db->getVector<char>( "ReadValues" );
-		auto WriteValues = domain_db->getVector<char>( "WriteValues" );
+		auto ReadValues = domain_db->getVector<int>( "ReadValues" );
+		auto WriteValues = domain_db->getVector<int>( "WriteValues" );
 		SW = domain_db->getScalar<double>("Sw");
 
 		// Generate the NWP configuration
@@ -82,8 +82,8 @@ int main(int argc, char **argv)
 		for (n=0; n<N; n++) Dm->id[n]=1;
 		Dm->CommInit();
 
-		char *id;
-		id = new char [N];
+		signed char *id;
+		id = new signed char [N];
 		sprintf(LocalRankFilename,"ID.%05i",rank);
 		size_t readID;
 		FILE *IDFILE = fopen(LocalRankFilename,"rb");
@@ -131,14 +131,14 @@ int main(int argc, char **argv)
 		// calculate distance to non-wetting fluid
 		if (domain_db->keyExists( "HistoryLabels" )){
 			if (rank==0) printf("Relabel solid components that touch fluid 1 \n");
-			auto LabelList = domain_db->getVector<char>( "ComponentLabels" );
-			auto HistoryLabels = domain_db->getVector<char>( "HistoryLabels" );
+			auto LabelList = domain_db->getVector<int>( "ComponentLabels" );
+			auto HistoryLabels = domain_db->getVector<int>( "HistoryLabels" );
 			size_t NLABELS=LabelList.size();
 			if (rank==0){
 				for (unsigned int idx=0; idx < NLABELS; idx++){
-					char VALUE = LabelList[idx];
-					char NEWVAL = HistoryLabels[idx];
-					printf("    Relabel component %d as %d \n", VALUE, NEWVAL);
+					signed char VALUE = LabelList[idx];
+					signed char NEWVAL = HistoryLabels[idx];
+					printf("    Relabel component %hhd as %hhd \n", VALUE, NEWVAL);
 				}
 			}
 			for (int k=0;k<nz;k++){
@@ -169,8 +169,8 @@ int main(int argc, char **argv)
 						int n = k*nx*ny+j*nx+i;
 						signed char LOCVAL = id[n];
 						for (unsigned int idx=0; idx < NLABELS; idx++){
-							char VALUE=LabelList[idx];
-							char NEWVALUE=HistoryLabels[idx];
+							signed char VALUE=LabelList[idx];
+							signed char NEWVALUE=HistoryLabels[idx];
 							if (LOCVAL == VALUE){
 								idx = NLABELS;
 								if (SignDist(i,j,k) < 1.0){
diff --git a/tests/lbpm_morphopen_pp.cpp b/tests/lbpm_morphopen_pp.cpp
index 3efec5b5..d4809fd5 100644
--- a/tests/lbpm_morphopen_pp.cpp
+++ b/tests/lbpm_morphopen_pp.cpp
@@ -59,8 +59,8 @@ int main(int argc, char **argv)
 		auto L = domain_db->getVector<double>( "L" );
 		auto size = domain_db->getVector<int>( "n" );
 		auto nproc = domain_db->getVector<int>( "nproc" );
-		auto ReadValues = domain_db->getVector<char>( "ReadValues" );
-		auto WriteValues = domain_db->getVector<char>( "WriteValues" );
+		auto ReadValues = domain_db->getVector<int>( "ReadValues" );
+		auto WriteValues = domain_db->getVector<int>( "WriteValues" );
 		SW = domain_db->getScalar<double>("Sw");
 
 		// Generate the NWP configuration
@@ -82,8 +82,8 @@ int main(int argc, char **argv)
 		for (n=0; n<N; n++) Dm->id[n]=1;
 		Dm->CommInit();
 
-		char *id;
-		id = new char [N];
+		signed char *id;
+		id = new signed char [N];
 		sprintf(LocalRankFilename,"ID.%05i",rank);
 		size_t readID;
 		FILE *IDFILE = fopen(LocalRankFilename,"rb");
@@ -131,13 +131,13 @@ int main(int argc, char **argv)
 		// calculate distance to non-wetting fluid
 		if (domain_db->keyExists( "HistoryLabels" )){
 			if (rank==0) printf("Relabel solid components that touch fluid 1 \n");
-			auto LabelList = domain_db->getVector<char>( "ComponentLabels" );
-			auto HistoryLabels = domain_db->getVector<char>( "HistoryLabels" );
+			auto LabelList = domain_db->getVector<int>( "ComponentLabels" );
+			auto HistoryLabels = domain_db->getVector<int>( "HistoryLabels" );
 			size_t NLABELS=LabelList.size();
 			if (rank==0){
 				for (unsigned int idx=0; idx < NLABELS; idx++){
-					char VALUE = LabelList[idx];
-					char NEWVAL = HistoryLabels[idx];
+					signed char VALUE = LabelList[idx];
+					signed char NEWVAL = HistoryLabels[idx];
 					printf("    Relabel component %d as %d \n", VALUE, NEWVAL);
 				}
 			}
@@ -169,8 +169,8 @@ int main(int argc, char **argv)
 						int n = k*nx*ny+j*nx+i;
 						signed char LOCVAL = id[n];
 						for (unsigned int idx=0; idx < NLABELS; idx++){
-							char VALUE=LabelList[idx];
-							char NEWVALUE=HistoryLabels[idx];
+							signed char VALUE=LabelList[idx];
+							signed char NEWVALUE=HistoryLabels[idx];
 							if (LOCVAL == VALUE){
 								idx = NLABELS;
 								if (SignDist(i,j,k) < 1.0){

From c417559a7dd8a7f29b4d0166ffdb2ff335380140 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Fri, 29 Mar 2019 07:28:24 -0400
Subject: [PATCH 49/53] explicit signed char in morph header

---
 analysis/morphology.h | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/analysis/morphology.h b/analysis/morphology.h
index 4692475d..f69af5f2 100644
--- a/analysis/morphology.h
+++ b/analysis/morphology.h
@@ -3,6 +3,6 @@
 #include "common/Domain.h"
 #include "analysis/runAnalysis.h"
 
-double MorphOpen(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, double VoidFraction);
-double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, double VoidFraction);
-double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id, std::shared_ptr<Domain> Dm, double TargetVol);
\ No newline at end of file
+double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction);
+double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction);
+double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id, std::shared_ptr<Domain> Dm, double TargetVol);

From e9c41069d492267e3976f8b57a5fb50bedf2a2c7 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Sat, 30 Mar 2019 07:32:19 -0400
Subject: [PATCH 50/53] added wait before copying data

---
 analysis/runAnalysis.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index a2b3801d..319014c5 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -902,6 +902,7 @@ void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, do
 
     //if ( matches(type,AnalysisType::CopySimState) ) {
     if ( timestep%d_analysis_interval == 0 ) {
+        finish(); // can't copy if threads are still working on data
         // Copy the members of Averages to the cpu (phase was copied above)
         PROFILE_START("Copy-Pressure",1);
         ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np);

From 7c0233d79bcacc16243d8a48ab20f0fe1ae063b7 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 3 Apr 2019 11:27:29 -0400
Subject: [PATCH 51/53] print BC info / fix relperm.csv timestamp

---
 models/ColorModel.cpp | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index adadc89d..40f6e367 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -560,6 +560,10 @@ void ScaLBL_ColorModel::Run(){
 		//analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
 		analysis.basic( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
 
+		if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition > 0){
+			printf("....inlet pressure=%f \n",din);
+		}
+		
 		// allow initial ramp-up to get closer to steady state
 		if (timestep > RAMP_TIMESTEPS && timestep%analysis_interval == 0 && USE_MORPH){
 			analysis.finish();
@@ -649,7 +653,6 @@ void ScaLBL_ColorModel::Run(){
 					}
 					morph_timesteps=0;
 				}
-				CURRENT_STEADY_TIMESTEPS=0;
 				Ca_previous = Ca;
 			}
 			if (MORPH_ADAPT ){

From 754e13a5ce5b1473f6e88f0b7d26cd84b1290c05 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 3 Apr 2019 14:16:12 -0400
Subject: [PATCH 52/53] fixed interface momentum

---
 analysis/SubPhase.cpp | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index 2a7e61b7..4890a8d3 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -112,9 +112,9 @@ void SubPhase::Write(int timestep)
 		fprintf(SUBPHASE,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn); 
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.p, gwd.p, gnc.p, gnd.p);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
-		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Px, gwd.Px, giwn.Px, gnc.Px, gnd.Px, giwn.Px);
-		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Py, gwd.Py, giwn.Py, gnc.Py, gnd.Py, giwn.Py);
-		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Pz, gwd.Pz, giwn.Pz, gnc.Pz, gnd.Pz, giwn.Pz);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.V, gwc.A, gwc.H, gwc.X);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
@@ -128,9 +128,9 @@ void SubPhase::Write(int timestep)
 		fprintf(SUBPHASE,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",wc.p, wd.p, nc.p, nd.p);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn);
-		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Px, wd.Px, iwn.Px, nc.Px, nd.Px, iwn.Px);
-		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Py, nc.Py, nd.Py, iwn.Py);
-		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pz, nc.Pz, nd.Pz, iwn.Pz);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",wc.V, wc.A, wc.H, wc.X);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);

From ace0c78ff30b11d3b66e398f78a0233dd45a7cc1 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Thu, 4 Apr 2019 20:22:48 -0400
Subject: [PATCH 53/53] fixed bug in permeability calculation

---
 common/Domain.cpp   | 1 -
 common/Domain.h     | 2 ++
 models/MRTModel.cpp | 4 ++--
 3 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/common/Domain.cpp b/common/Domain.cpp
index 931ac5dd..f7427e0c 100644
--- a/common/Domain.cpp
+++ b/common/Domain.cpp
@@ -545,7 +545,6 @@ void Domain::ReadIDs(){
 
 	// Compute the porosity
 	double sum;
-	double porosity;
 	double sum_local=0.0;
 	double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
 	if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
diff --git a/common/Domain.h b/common/Domain.h
index 038627c2..e298836c 100755
--- a/common/Domain.h
+++ b/common/Domain.h
@@ -111,6 +111,7 @@ public: // Public variables (need to create accessors instead)
     int Nx,Ny,Nz,N;
     int inlet_layers_x, inlet_layers_y, inlet_layers_z;
     int outlet_layers_x, outlet_layers_y, outlet_layers_z;
+    double porosity;
     RankInfoStruct rank_info;
 
     MPI_Comm Comm;        // MPI Communicator for this domain
@@ -122,6 +123,7 @@ public: // Public variables (need to create accessors instead)
     //**********************************
     // MPI ranks for all 18 neighbors
     //**********************************
+    inline double Porosity() const { return porosity; }
     inline int iproc() const { return rank_info.ix; }
     inline int jproc() const { return rank_info.jy; }
     inline int kproc() const { return rank_info.kz; }
diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp
index fc1fa05f..8bb06fec 100644
--- a/models/MRTModel.cpp
+++ b/models/MRTModel.cpp
@@ -241,8 +241,8 @@ void ScaLBL_MRTModel::Run(){
 			Hs=sumReduce( Dm->Comm, Hs);
 			Xs=sumReduce( Dm->Comm, Xs);
 			if (rank==0) {
-				double h = Lz/double(Nz);
-				double absperm = h*h*mu*sqrt(vax*vax+vay*vay+vaz*vaz)/sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
+				double h = Lz/double((Nz-2)*nprocz);
+				double absperm = h*h*mu*Mask->Porosity()*sqrt(vax*vax+vay*vay+vaz*vaz)/sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
 				printf("     %f\n",absperm);
 				FILE * log_file = fopen("Permeability.csv","a");
 				fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu,