diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index 8effe1fd..5d029456 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -9,8 +9,6 @@
 #include "common/Utilities.h"
 #include "models/ColorModel.h"
 
-//#define WRE_SURFACES
-
 /*
  * Simulator for two-phase flow in porous media
  * James E. McClure 2013-2014
@@ -71,29 +69,33 @@ int main( int argc, char **argv )
         if (SimulationMode == "development"){
         	double MLUPS=0.0;
         	int timestep = 0;
+		bool ContinueSimulation = true;
+        	int ANALYSIS_INTERVAL = ColorModel.timestepMax;
         	/* flow adaptor keys to control */
         	int SKIP_TIMESTEPS = 0;
         	int MAX_STEADY_TIME = 1000000;
+		double ENDPOINT_THRESHOLD = 0.1;
         	double FRACTIONAL_FLOW_INCREMENT = 0.05;
         	if (ColorModel.db->keyExists( "FlowAdaptor" )){
         		auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" );
         		MAX_STEADY_TIME = flow_db->getWithDefault<int>( "max_steady_timesteps", 1000000 );
         		SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 100000 );
         		FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
+        	        ENDPOINT_THRESHOLD = flow_db->getWithDefault<double>( "endpoint_threshold", 0.1);
         	}
-        	int ANALYSIS_INTERVAL = ColorModel.timestepMax;
         	if (ColorModel.analysis_db->keyExists( "analysis_interval" )){
         		ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar<int>( "analysis_interval" );
         	}
         	/* Launch the simulation */
         	FlowAdaptor Adapt(ColorModel);
         	runAnalysis analysis(ColorModel);
-        	while (ColorModel.timestep < ColorModel.timestepMax){
-            	/* this will run steady points */
-            	timestep += MAX_STEADY_TIME;
-            	MLUPS = ColorModel.Run(timestep);
-            	if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
-            	
+		if (ContinueSimulation){
+            	    /* this will run steady points */
+                    timestep += MAX_STEADY_TIME;
+		    MLUPS = ColorModel.Run(timestep);
+            	    if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
+            	    if (ColorModel.timestep > ColorModel.timestepMax) ContinueSimulation = false;
+		
             	/* update the fractional flow by adding mass */
             	int skip_time = 0;
             	timestep = ColorModel.timestep;
@@ -101,6 +103,17 @@ int main( int argc, char **argv )
     			double volB = ColorModel.Averages->gwb.V; 
     			double volA = ColorModel.Averages->gnb.V; 
     			double initialSaturation = volB/(volA + volB);
+
+			double vA_x = ColorModel.Averages->gnb.Px/ColorModel.Averages->gnb.M; 
+			double vA_y = ColorModel.Averages->gnb.Py/ColorModel.Averages->gnb.M; 
+			double vA_z = ColorModel.Averages->gnb.Pz/ColorModel.Averages->gnb.M; 
+			double vB_x = ColorModel.Averages->gwb.Px/ColorModel.Averages->gwb.M; 
+			double vB_y = ColorModel.Averages->gwb.Py/ColorModel.Averages->gwb.M; 
+			double vB_z = ColorModel.Averages->gwb.Pz/ColorModel.Averages->gwb.M; 			
+			double speedA = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z);
+			double speedB = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
+			if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false;
+		if (ContinueSimulation){
             	while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){
             		timestep += ANALYSIS_INTERVAL;
                 	Adapt.UpdateFractionalFlow(ColorModel);
@@ -109,11 +122,11 @@ int main( int argc, char **argv )
         			double volA = ColorModel.Averages->gnb.V;
         			SaturationChange = volB/(volA + volB) - initialSaturation;
             		skip_time += ANALYSIS_INTERVAL;
-            	}
+		}
             	if (rank==0) printf("  *********************************************************************  \n");
             	if (rank==0) printf("   Updated fractional flow with saturation change = %f  \n", SaturationChange);
             	if (rank==0) printf("  *********************************************************************  \n");
-
+		}
             	/* apply timestep skipping algorithm to accelerate steady-state */
             	/* skip_time = 0;
             	timestep = ColorModel.timestep;
@@ -125,7 +138,7 @@ int main( int argc, char **argv )
             	}
             	*/
             }
-            ColorModel.WriteDebug();
+		  //ColorModel.WriteDebug();
         }
 
         else