Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

analysis/SubPhase.cpp

@@ -23,6 +23,7 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 	Vel_x.resize(Nx,Ny,Nz);         Vel_x.fill(0);	    // Gradient of the phase indicator field
 	Vel_y.resize(Nx,Ny,Nz);         Vel_y.fill(0);
 	Vel_z.resize(Nx,Ny,Nz);         Vel_z.fill(0);
+	Dissipation.resize(Nx,Ny,Nz);   Dissipation.fill(0);
 	SDs.resize(Nx,Ny,Nz);         	SDs.fill(0);
 	//.........................................
 
@@ -42,11 +43,12 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 			//fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n");
 			fprintf(SUBPHASE,"time rn rw nun nuw Fx Fy Fz iftwn wet ");				
 			fprintf(SUBPHASE,"pwc pwd pnc pnd ");						// pressures 
-			fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni ");				// mass 
+			fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni Msw Msn ");				// mass 
-			fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x ");	// momentum 
+			fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x Psw_x Psn_x ");	// momentum 
-			fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");			
+			fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y Psw_y Psn_y ");			
-			fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");			
+			fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z Psw_z Psn_z ");			
 			fprintf(SUBPHASE,"Kwc Kwd Kwi Knc Knd Kni ");				// kinetic energy
+			fprintf(SUBPHASE,"Dwc Dwd Dnc Dnd ");						// viscous dissipation
 			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
@@ -67,11 +69,12 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 		//fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n");
 		fprintf(SUBPHASE,"time rn rw nun nuw Fx Fy Fz iftwn wet ");				
 		fprintf(SUBPHASE,"pwc pwd pnc pnd ");						// pressures 
-		fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni ");				// mass 
+		fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni Msw Msn ");				// mass 
-		fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x ");	// momentum 
+		fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x Psw_x Psn_x ");	// momentum 
-		fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");			
+		fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y Psw_y Psn_y ");			
-		fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");			
+		fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z Psw_z Psn_z ");				
 		fprintf(SUBPHASE,"Kwc Kwd Kwi Knc Knd Kni ");				// kinetic energy
+		fprintf(SUBPHASE,"Dwc Dwd Dnc Dnd ");						// viscous dissipation		
 		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
@@ -111,11 +114,12 @@ void SubPhase::Write(int timestep)
 	if (Dm->rank()==0){
 		fprintf(SUBPHASE,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn,total_wetting_interaction_global); 
 		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.p, gwd.p, gnc.p, gnd.p);
-		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
+		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn, gifs.Mw, gifs.Mn);
-		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Px, gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx);
+		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Px, gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx, gifs.Pwx, gifs.Pnx);
-		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Py, gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny);
+		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Py, gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny, gifs.Pwy, gifs.Pny);
-		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Pz, gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz);
+		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Pz, gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz, gifs.Pwz, gifs.Pnz);
 		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
+		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.visc, gwd.visc, gnc.visc, gnd.visc);
 		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.V, gwc.A, gwc.H, gwc.X);
 		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
 		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gnc.V, gnc.A, gnc.H, gnc.X);
@@ -127,11 +131,12 @@ void SubPhase::Write(int timestep)
 	else{
 		fprintf(SUBPHASE,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn,total_wetting_interaction);
 		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.p, wd.p, nc.p, nd.p);
-		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn);
+		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn, ifs.Mw, ifs.Mn);
-		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx);
+		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx, ifs.Pwx, ifs.Pnx);
-		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny);
+		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny, ifs.Pwy, ifs.Pny);
-		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz);
+		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz, ifs.Pwz, ifs.Pnz);
 		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn);
+		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.visc, wd.visc, nc.visc, nd.visc);
 		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.V, wc.A, wc.H, wc.X);
 		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);
 		fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",nc.V, nc.A, nc.H, nc.X);
@@ -481,7 +486,7 @@ void SubPhase::Full(){
 	if (Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin += Dm->inlet_layers_z; 
 	if (Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax -= Dm->outlet_layers_z; 
 	*/
-	nd.reset();	nc.reset(); wd.reset();	wc.reset();	iwn.reset();	iwnc.reset();
+	nd.reset();	nc.reset(); wd.reset();	wc.reset();	iwn.reset(); iwnc.reset(); ifs.reset();
 
  	Dm->CommunicateMeshHalo(Phi);
 	for (int k=1; k<Nz-1; k++){
@@ -496,6 +501,33 @@ void SubPhase::Full(){
 		}
 	}
  	Dm->CommunicateMeshHalo(DelPhi);
+ 	
+ 	
+ 	Dm->CommunicateMeshHalo(Vel_x);
+ 	Dm->CommunicateMeshHalo(Vel_y);
+ 	Dm->CommunicateMeshHalo(Vel_z);
+	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 velocity gradients using finite differences
+				double phi = Phi(i,j,k);
+				double nu = nu_n + 0.5*(1.0-phi)*(nu_w-nu_n);
+				double rho = rho_n + 0.5*(1.0-phi)*(rho_w-rho_n);
+				double ux = 0.5*(Vel_x(i+1,j,k) - Vel_x(i-1,j,k));
+				double uy = 0.5*(Vel_x(i,j+1,k) - Vel_x(i,j-1,k));
+				double uz = 0.5*(Vel_x(i,j,k+1) - Vel_x(i,j,k-1));
+				double vx = 0.5*(Vel_y(i+1,j,k) - Vel_y(i-1,j,k));
+				double vy = 0.5*(Vel_y(i,j+1,k) - Vel_y(i,j-1,k));
+				double vz = 0.5*(Vel_y(i,j,k+1) - Vel_y(i,j,k-1));
+				double wx = 0.5*(Vel_z(i+1,j,k) - Vel_z(i-1,j,k));
+				double wy = 0.5*(Vel_z(i,j+1,k) - Vel_z(i,j-1,k));
+				double wz = 0.5*(Vel_z(i,j,k+1) - Vel_z(i,j,k-1));
+				if (SDs(i,j,k) > 2.0){
+					Dissipation(i,j,k) = 2*rho*nu*( ux*ux + vy*vy + wz*wz + 0.5*(vx + uy)*(vx + uy)+ 0.5*(vz + wy)*(vz + wy)+ 0.5*(uz + wx)*(uz + wx));					
+				}
+			}
+		}
+	}
 
  	/*  Set up geometric analysis of each region */
 	
@@ -654,13 +686,33 @@ void SubPhase::Full(){
 					double ux = Vel_x(n);
 					double uy = Vel_y(n);
 					double uz = Vel_z(n);
+					double visc = Dissipation(n);
 					
-					if (DelPhi(n) > 1e-3 && SDs(n) < 3.0 ){
+					if (DelPhi(n) > 1e-3 ){
-						// film region
+						// get the normal vector 
 						double nx = 0.5*(Phi(i+1,j,k)-Phi(i-1,j,k));
 						double ny = 0.5*(Phi(i,j+1,k)-Phi(i,j-1,k));
 						double nz = 0.5*(Phi(i,j,k+1)-Phi(i,j,k-1));
-						InterfaceTransportMeasures(  beta,  rho_w,  rho_n,  nA, nB, nx, ny, nz, ux, uy, uz, iwn);
+						if (SDs(n) > 2.5){
+							// not a film region
+							InterfaceTransportMeasures(  beta,  rho_w,  rho_n,  nA, nB, nx, ny, nz, ux, uy, uz, iwn);
+						}
+						else{
+							// films that are close to the wetting fluid
+							if ( morph_w->distance(i,j,k) < 2.5 && phi > 0.0){
+								ifs.Mw += rho_w;
+								ifs.Pwx += rho_w*ux;
+								ifs.Pwy += rho_w*uy;
+								ifs.Pwz += rho_w*uz;
+							}
+							// films that are close to the NWP 
+							if ( morph_n->distance(i,j,k) < 2.5 && phi < 0.0){
+								ifs.Mn += rho_n;
+								ifs.Pnx += rho_n*ux;
+								ifs.Pny += rho_n*uy;
+								ifs.Pnz += rho_n*uz;
+							}
+						}
 					}
 					else if ( phi > 0.0){
 						if (morph_n->label(i,j,k) > 0 ){
@@ -691,6 +743,7 @@ void SubPhase::Full(){
 							nd.Py += nA*rho_n*uy;
 							nd.Pz += nA*rho_n*uz;
 							nd.K += nA*rho_n*(ux*ux + uy*uy + uz*uz);
+							nd.visc += visc;
 						}
 						else{
 							nA = 1.0;
@@ -699,6 +752,7 @@ void SubPhase::Full(){
 							nc.Py += nA*rho_n*uy;
 							nc.Pz += nA*rho_n*uz;
 							nc.K += nA*rho_n*(ux*ux + uy*uy + uz*uz);
+							nc.visc += visc;
 						}
 					}
 					else{
@@ -710,6 +764,7 @@ void SubPhase::Full(){
 							wd.Py += nB*rho_w*uy;
 							wd.Pz += nB*rho_w*uz;
 							wd.K += nB*rho_w*(ux*ux + uy*uy + uz*uz);
+							wd.visc += visc;
 						}
 						else{
 							nB = 1.0;
@@ -718,6 +773,7 @@ void SubPhase::Full(){
 							wc.Py += nB*rho_w*uy;
 							wc.Pz += nB*rho_w*uz;
 							wc.K += nB*rho_w*(ux*ux + uy*uy + uz*uz);
+							wc.visc += visc;
 						}
 					}
 				}
@@ -730,25 +786,29 @@ void SubPhase::Full(){
 	gnd.Py=Dm->Comm.sumReduce(  nd.Py);
 	gnd.Pz=Dm->Comm.sumReduce(  nd.Pz);
 	gnd.K=Dm->Comm.sumReduce(  nd.K);
+	gnd.visc=Dm->Comm.sumReduce(  nd.visc);
 
 	gwd.M=Dm->Comm.sumReduce(  wd.M);
 	gwd.Px=Dm->Comm.sumReduce(  wd.Px);
 	gwd.Py=Dm->Comm.sumReduce(  wd.Py);
 	gwd.Pz=Dm->Comm.sumReduce(  wd.Pz);
 	gwd.K=Dm->Comm.sumReduce(  wd.K);
+	gwd.visc=Dm->Comm.sumReduce(  wd.visc);
 	
 	gnc.M=Dm->Comm.sumReduce(  nc.M);
 	gnc.Px=Dm->Comm.sumReduce(  nc.Px);
 	gnc.Py=Dm->Comm.sumReduce(  nc.Py);
 	gnc.Pz=Dm->Comm.sumReduce(  nc.Pz);
 	gnc.K=Dm->Comm.sumReduce(  nc.K);
+	gnc.visc=Dm->Comm.sumReduce(  nc.visc);
 
 	gwc.M=Dm->Comm.sumReduce(  wc.M);
 	gwc.Px=Dm->Comm.sumReduce(  wc.Px);
 	gwc.Py=Dm->Comm.sumReduce(  wc.Py);
 	gwc.Pz=Dm->Comm.sumReduce(  wc.Pz);
 	gwc.K=Dm->Comm.sumReduce(  wc.K);
-	
+	gwc.visc=Dm->Comm.sumReduce(  wc.visc);
+
 	giwn.Mn=Dm->Comm.sumReduce(  iwn.Mn);
 	giwn.Pnx=Dm->Comm.sumReduce(  iwn.Pnx);
 	giwn.Pny=Dm->Comm.sumReduce(  iwn.Pny);
@@ -760,6 +820,15 @@ void SubPhase::Full(){
 	giwn.Pwz=Dm->Comm.sumReduce(  iwn.Pwz);
 	giwn.Kw=Dm->Comm.sumReduce(  iwn.Kw);
 	
+	gifs.Mn= Dm->Comm.sumReduce(  ifs.Mn);
+	gifs.Pnx=Dm->Comm.sumReduce(  ifs.Pnx);
+	gifs.Pny=Dm->Comm.sumReduce(  ifs.Pny);
+	gifs.Pnz=Dm->Comm.sumReduce(  ifs.Pnz);
+	gifs.Mw= Dm->Comm.sumReduce(  ifs.Mw);
+	gifs.Pwx=Dm->Comm.sumReduce(  ifs.Pwx);
+	gifs.Pwy=Dm->Comm.sumReduce(  ifs.Pwy);
+	gifs.Pwz=Dm->Comm.sumReduce(  ifs.Pwz);	
+	
 	// pressure averaging
 	gnc.p=Dm->Comm.sumReduce(  nc.p);
 	gnd.p=Dm->Comm.sumReduce(  nd.p);

analysis/SubPhase.h

@@ -22,10 +22,11 @@ class phase{
 public:
 	int Nc;
 	double p;
-	double M,Px,Py,Pz,K;
+	double M,Px,Py,Pz,K,visc;
 	double V,A,H,X;
 	void reset(){
 		p=M=Px=Py=Pz=K=0.0;
+		visc=0.0;
 		V=A=H=X=0.0;
 		Nc=1;
 	}
@@ -70,10 +71,12 @@ public:
 	// local entities
 	phase wc,wd,wb,nc,nd,nb,solid;
 	interface iwn,iwnc;
+	interface ifs;
 	
 	// global entities
 	phase gwc,gwd,gwb,gnc,gnd,gnb,gsolid;
 	interface giwn,giwnc;
+	interface gifs;
 	/* fluid-solid wetting interaction */
 	double total_wetting_interaction, count_wetting_interaction;
 	double total_wetting_interaction_global, count_wetting_interaction_global;
@@ -92,6 +95,7 @@ public:
 	DoubleArray Vel_x;		// velocity field
 	DoubleArray Vel_y;
 	DoubleArray Vel_z;
+	DoubleArray Dissipation;
 	DoubleArray SDs;
 
 	std::shared_ptr<Minkowski> morph_w;

analysis/runAnalysis.cpp

@@ -312,15 +312,21 @@ public:
             fillData.copy( Averages.Vel_z, VelzData );
         }
 
+        if ( vis_db->getWithDefault<bool>( "save_dissipation", false ) ) {
+            ASSERT( visData[0].vars[5]->name == "ViscousDissipation" );
+            Array<double> &ViscousDissipation = visData[0].vars[5]->data;
+            fillData.copy( Averages.Dissipation, ViscousDissipation );
+        }
+
         if ( vis_db->getWithDefault<bool>( "save_distance", false ) ) {
-            ASSERT( visData[0].vars[5]->name == "SignDist" );
+            ASSERT( visData[0].vars[6]->name == "SignDist" );
-            Array<double> &SignData = visData[0].vars[5]->data;
+            Array<double> &SignData = visData[0].vars[6]->data;
             fillData.copy( Averages.SDs, SignData );
         }
 
         if ( vis_db->getWithDefault<bool>( "save_connected_components", false ) ) {
-            ASSERT( visData[0].vars[6]->name == "BlobID" );
+            ASSERT( visData[0].vars[7]->name == "BlobID" );
-            Array<double> &BlobData = visData[0].vars[6]->data;
+            Array<double> &BlobData = visData[0].vars[7]->data;
             fillData.copy( Averages.morph_n->label, BlobData );
         }
 
@@ -668,6 +674,7 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> input_db, const RankInfoStru
     auto VxVar       = std::make_shared<IO::Variable>();
     auto VyVar       = std::make_shared<IO::Variable>();
     auto VzVar       = std::make_shared<IO::Variable>();
+    auto ViscousDissipationVar = std::make_shared<IO::Variable>();
     auto SignDistVar = std::make_shared<IO::Variable>();
     auto BlobIDVar   = std::make_shared<IO::Variable>();
 
@@ -704,6 +711,14 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> input_db, const RankInfoStru
         VzVar->data.resize( d_n[0], d_n[1], d_n[2] );
         d_meshData[0].vars.push_back( VzVar );
     }
+    
+    if ( vis_db->getWithDefault<bool>( "save_dissipation", false ) ) {
+        ViscousDissipationVar->name = "ViscousDissipation";
+    	ViscousDissipationVar->type = IO::VariableType::VolumeVariable;
+    	ViscousDissipationVar->dim  = 1;
+    	ViscousDissipationVar->data.resize( d_n[0], d_n[1], d_n[2] );
+        d_meshData[0].vars.push_back( ViscousDissipationVar );
+    }
 
     if ( vis_db->getWithDefault<bool>( "save_distance", false ) ) {
         SignDistVar->name = "SignDist";

common/Domain.cpp

@@ -399,7 +399,7 @@ void Domain::Decomp( const std::string& Filename )
 				for (int i = 0; i<global_Nx; i++){
 					n = k*global_Nx*global_Ny+j*global_Nx+i;
 					//char locval = loc_id[n];
-					char locval = SegData[n];
+					signed char locval = SegData[n];
 					for (size_t idx=0; idx<ReadValues.size(); idx++){
 						signed char oldvalue=ReadValues[idx];
 						signed char newvalue=WriteValues[idx];

common/ScaLBL.h

@@ -104,12 +104,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map,
 		double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
 
 extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, 
-        double *Phi, double *GreySolidW, double *Poros,double *Perm,double *Vel, double *Pressure,
+        double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure,
         double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
 		double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np);
 
 extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, 
-		double *Phi, double *GreySolidW, double *Poros,double *Perm,double *Vel,double *Pressure, 
+		double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure, 
         double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
 		double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np);
 

cpu/Color.cpp

@@ -2613,7 +2613,6 @@ extern "C" void ScaLBL_D3Q7_AAeven_Color(int *Map, double *Aq, double *Bq, doubl
 	double a1,b1,a2,b2,nAB,delta;
 	double C,nx,ny,nz; //color gradient magnitude and direction
 	double ux,uy,uz;
-	double phi;
 	// Instantiate mass transport distributions
 	// Stationary value - distribution 0
 	for (int n=start; n<finish; n++){
@@ -2780,7 +2779,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_PhaseField(int *neighborList, int *Map, double
 
 extern "C" void ScaLBL_D3Q7_AAeven_PhaseField(int *Map, double *Aq, double *Bq, double *Den, double *Phi, 
 			int start, int finish, int Np){
-	int idx,nread;
+	int idx;
 	double fq,nA,nB;
 	for (int n=start; n<finish; n++){
 		

cpu/Greyscale.cpp

@@ -17,7 +17,6 @@
 
 extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
                                               double *Poros,double *Perm, double *Velocity, double *Pressure){
-	int n;
 	// conserved momemnts
 	double rho,vx,vy,vz,v_mag;
     double ux,uy,uz,u_mag;
@@ -262,7 +261,6 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finis
 
 extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, int start, int finish, int Np, double rlx,  double rlx_eff, double Gx, double Gy, double Gz, 
                                              double *Poros,double *Perm, double *Velocity,double *Pressure){
-	int n;
 	// conserved momemnts
 	double rho,vx,vy,vz,v_mag;
     double ux,uy,uz,u_mag;
@@ -278,7 +276,6 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, in
     double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
     double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
 
-	int nread;
 	for (int n=start; n<finish; n++){
 		
 		// q=0
@@ -568,7 +565,6 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, in
 
 extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int finish, int Np, double rlx,  double rlx_eff, double Gx, double Gy, double Gz,
                                               double *Poros,double *Perm, double *Velocity, double Den,double *Pressure){
-	int n;
 	double vx,vy,vz,v_mag;
     double ux,uy,uz,u_mag;
     double pressure;//defined for this incompressible model
@@ -1057,7 +1053,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int
 
 extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx,  double rlx_eff, double Gx, double Gy, double Gz, 
                                              double *Poros,double *Perm, double *Velocity, double Den,double *Pressure){
-	int n, nread;
+	int nread;
 	double vx,vy,vz,v_mag;
     double ux,uy,uz,u_mag;
     double pressure;//defined for this incompressible model
@@ -1212,6 +1208,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dis
 
 		// q=9
 		nread = neighborList[n+8*Np];
+		
 		fq = dist[nread];
 		pressure += fq;
 		m1 += 8.0*fq;
@@ -1583,7 +1580,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dis
 
 extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,double *Poros,double *Perm, double *Velocity,double rho0,double *Pressure){
 
-	int n, nread;
+	int nread;
 	int nr1,nr2,nr3,nr4,nr5,nr6;
 	int nr7,nr8,nr9,nr10;
 	int nr11,nr12,nr13,nr14;
@@ -1720,6 +1717,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist
 		//nread = neighborList[n+6*Np];
 		//fq = dist[nread];
 		nr7 = neighborList[n+6*Np];
+		
 		fq = dist[nr7];
 		rho += fq;
 		m1 += 8.0*fq;
@@ -1969,6 +1967,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist
         Fz = rho0*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
         if (porosity==1.0){
             Fx=rho0*Gx;
+            
             Fy=rho0*Gy;
             Fz=rho0*Gz;
         }
@@ -2135,7 +2134,6 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist
 
 extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_MRT(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,double *Poros,double *Perm, double *Velocity,double rho0,double *Pressure){
 
-	int n;
 	double vx,vy,vz,v_mag;
     double ux,uy,uz,u_mag;
     double pressure;//defined for this incompressible model

cpu/GreyscaleColor.cpp

@@ -1340,10 +1340,10 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, doubl
 
 //CP: capillary penalty
 // also turn off recoloring for grey nodes
-extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, 
+extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
-		double *Phi, double *Psi, double *GreySolidGrad, double *Poros,double *Perm,double *Velocity,double *Pressure,
+		 double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure,
-        double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
+        double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff,double alpha, double beta,
-		double Gx, double Gy, double Gz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){
+		double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
 
 	int n,nn,ijk,nread;
 	int nr1,nr2,nr3,nr4,nr5,nr6;
@@ -1370,12 +1370,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
     //double c0, c1; //Guo's model parameters
     double tau_eff;
     double mu_eff;//kinematic viscosity
-    double nx_gs,ny_gs,nz_gs;//grey-solid color gradient
-    double nx_phase,ny_phase,nz_phase,C_phase;
     double Fx,Fy,Fz;
-	double gp1,gp2,gp4,gp6,gp8,gp9,gp10,gp11,gp12,gp13,gp14,gp15,gp16,gp17,gp18;
-	double gp3,gp5,gp7;
     double Fcpx,Fcpy,Fcpz;//capillary penalty force
+    double W;//greyscale wetting strength
+    double Sn_grey,Sw_grey;
 
 	const double mrt_V1=0.05263157894736842;
 	const double mrt_V2=0.012531328320802;
@@ -1394,11 +1392,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
 		// read the component number densities
 		nA = Den[n];
 		nB = Den[Np + n];
+
         porosity = Poros[n];
         perm = Perm[n];
-        nx_gs = GreySolidGrad[n+0*Np];
+        W = GreySolidW[n];
-        ny_gs = GreySolidGrad[n+1*Np];
+        Sn_grey = GreySn[n];
-        nz_gs = GreySolidGrad[n+2*Np];
+        Sw_grey = GreySw[n];
 
 		// compute phase indicator field
 		phi=(nA-nB)/(nA+nB);
@@ -1420,100 +1419,63 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
 		//........................................................................
 		nn = ijk-1;							// neighbor index (get convention)
 		m1 = Phi[nn];						// get neighbor for phi - 1
-        gp1 = Psi[nn];
 		//........................................................................
 		nn = ijk+1;							// neighbor index (get convention)
 		m2 = Phi[nn];						// get neighbor for phi - 2
-        gp2 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideY;							// neighbor index (get convention)
 		m3 = Phi[nn];					// get neighbor for phi - 3
-        gp3 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideY;							// neighbor index (get convention)
 		m4 = Phi[nn];					// get neighbor for phi - 4
-        gp4 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideZ;						// neighbor index (get convention)
 		m5 = Phi[nn];					// get neighbor for phi - 5
-        gp5 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideZ;						// neighbor index (get convention)
 		m6 = Phi[nn];					// get neighbor for phi - 6
-        gp6 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideY-1;						// neighbor index (get convention)
 		m7 = Phi[nn];					// get neighbor for phi - 7
-        gp7 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideY+1;						// neighbor index (get convention)
 		m8 = Phi[nn];					// get neighbor for phi - 8
-        gp8 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideY-1;						// neighbor index (get convention)
 		m9 = Phi[nn];					// get neighbor for phi - 9
-        gp9 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideY+1;						// neighbor index (get convention)
 		m10 = Phi[nn];					// get neighbor for phi - 10
-        gp10 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideZ-1;						// neighbor index (get convention)
 		m11 = Phi[nn];					// get neighbor for phi - 11
-        gp11 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideZ+1;						// neighbor index (get convention)
 		m12 = Phi[nn];					// get neighbor for phi - 12
-        gp12 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideZ-1;						// neighbor index (get convention)
 		m13 = Phi[nn];					// get neighbor for phi - 13
-        gp13 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideZ+1;						// neighbor index (get convention)
 		m14 = Phi[nn];					// get neighbor for phi - 14
-        gp14 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideZ-strideY;					// neighbor index (get convention)
 		m15 = Phi[nn];					// get neighbor for phi - 15
-        gp15 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideZ+strideY;					// neighbor index (get convention)
 		m16 = Phi[nn];					// get neighbor for phi - 16
-        gp16 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideZ-strideY;					// neighbor index (get convention)
 		m17 = Phi[nn];					// get neighbor for phi - 17
-        gp17 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideZ+strideY;					// neighbor index (get convention)
 		m18 = Phi[nn];					// get neighbor for phi - 18
-        gp18 = Psi[nn];
 		//............Compute the Color Gradient...................................
-		nx_phase = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+		nx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
-		ny_phase = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+		ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
-		nz_phase = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+		nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
-		C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase);
 
-        //correct the normal color gradient by considering the effect of grey solid
+		//............Compute the Greyscale Potential Gradient.....................
-        nx = nx_phase + (1.0-porosity)*nx_gs; 
-        ny = ny_phase + (1.0-porosity)*ny_gs; 
-        nz = nz_phase + (1.0-porosity)*nz_gs; 
-        if (C_phase==0.0){
-            nx = nx_phase; 
-            ny = ny_phase;
-            nz = nz_phase;
-        }
-
-		//...........Normalize the Color Gradient.................................
-		C = sqrt(nx*nx+ny*ny+nz*nz);
-		double ColorMag = C;
-		if (C==0.0) ColorMag=1.0;
-		nx = nx/ColorMag;
-		ny = ny/ColorMag;
-		nz = nz/ColorMag;		
-
-			//............Compute the Greyscale Potential Gradient.....................
 //			Fcpx = 0.0;
 //			Fcpy = 0.0;
 //			Fcpz = 0.0;
@@ -1534,14 +1496,24 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
 //                //ny = Fcpy/Fcp_mag;
 //                //nz = Fcpz/Fcp_mag;
 //            }
-        Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+        Fcpx = nx;
-        Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+        Fcpy = ny;
-        Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+        Fcpz = nz;
+        double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
+        if (Fcp_mag==0.0) Fcpx=Fcpy=Fcpz=0.0;
         //NOTE for open node (porosity=1.0),Fcp=0.0
         Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm); 
         Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm); 
         Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm); 
 
+		//...........Normalize the Color Gradient.................................
+		C = sqrt(nx*nx+ny*ny+nz*nz);
+		double ColorMag = C;
+		if (C==0.0) ColorMag=1.0;
+		nx = nx/ColorMag;
+		ny = ny/ColorMag;
+		nz = nz/ColorMag;		
+
 		// q=0
 		fq = dist[n];
 		rho = fq;
@@ -1893,6 +1865,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
 		//........................................................................
 		//..............carry out relaxation process..............................
 		//..........Toelke, Fruediger et. al. 2006................................
+        //---------------- NO higher-order force -------------------------------//
 		if (C == 0.0)	nx = ny = nz = 0.0;
 		m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1);
 		m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2);
@@ -1917,6 +1890,43 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
 		m16 = m16 + rlx_setB*( - m16);
 		m17 = m17 + rlx_setB*( - m17);
 		m18 = m18 + rlx_setB*( - m18);
+        //----------------------------------------------------------------------//
+
+        //----------------With higher-order force ------------------------------//
+		//if (C == 0.0)	nx = ny = nz = 0.0;
+		//m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1)
+        //        + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity;
+		//m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2)
+        //        + (1-0.5*rlx_setA)*11*(-Fx*ux-Fy*uy-Fz*uz)/porosity;
+        //jx = jx + Fx;
+		//m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
+        //        + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
+        //jy = jy + Fy;
+		//m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6)
+        //        + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
+        //jz = jz + Fz;
+		//m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
+        //        + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
+		//m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9)
+        //        + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity;
+		////m10 = m10 + rlx_setA*( - m10);
+        //m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10)
+        //          + (1-0.5*rlx_setA)*(-2*Fx*ux+Fy*uy+Fz*uz)/porosity;
+		//m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11)
+        //          + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity;
+		////m12 = m12 + rlx_setA*( - m12);
+        //m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12)
+        //          + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity;
+		//m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13);
+        //          + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity;
+		//m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14);
+        //          + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity;
+		//m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15);
+        //          + (1-0.5*rlx_setA)*(Fz*ux+Fx*uz)/porosity;
+		//m16 = m16 + rlx_setB*( - m16);
+		//m17 = m17 + rlx_setB*( - m17);
+		//m18 = m18 + rlx_setB*( - m18);
+        //----------------------------------------------------------------------//
 
 		//.................inverse transformation......................................................
 		// q=0
@@ -2049,6 +2059,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
 		// Cq = {1,0,0}, {0,1,0}, {0,0,1}
 		delta = beta*nA*nB*nAB*0.1111111111111111*nx;
 		if (!(nA*nB*nAB>0)) delta=0;
+        //----------------newly added for better control of recoloring---------------//
+        if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+        if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+        //---------------------------------------------------------------------------//
         if (RecoloringOff==true && porosity !=1.0) delta=0;
 		a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
 		b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@@ -2068,6 +2082,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
 		// Cq = {0,1,0}
 		delta = beta*nA*nB*nAB*0.1111111111111111*ny;
 		if (!(nA*nB*nAB>0)) delta=0;
+        //----------------newly added for better control of recoloring---------------//
+        if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+        if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+        //---------------------------------------------------------------------------//
         if (RecoloringOff==true && porosity !=1.0) delta=0;
 		a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
 		b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@@ -2088,6 +2106,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
 		// Cq = {0,0,1}
 		delta = beta*nA*nB*nAB*0.1111111111111111*nz;
 		if (!(nA*nB*nAB>0)) delta=0;
+        //----------------newly added for better control of recoloring---------------//
+        if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+        if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+        //---------------------------------------------------------------------------//
         if (RecoloringOff==true && porosity !=1.0) delta=0;
 		a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
 		b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@@ -2109,9 +2131,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
 //CP: capillary penalty
 // also turn off recoloring for grey nodes
 extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, 
-        double *Phi,double *Psi, double *GreySolidGrad, double *Poros,double *Perm,double *Velocity,double *Pressure,
+        double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure, 
         double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
-		double Gx, double Gy, double Gz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){
+		double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
 
 	int ijk,nn,n;
 	double fq;
@@ -2127,18 +2149,16 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 	double a1,b1,a2,b2,nAB,delta;
 	double C,nx,ny,nz; //color gradient magnitude and direction
 	double phi,tau,rho0,rlx_setA,rlx_setB;
-
+    double W;//greyscale wetting strength
+    double Sn_grey,Sw_grey;
+    
     //double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002)
     double porosity;
     double perm;//voxel permeability
     //double c0, c1; //Guo's model parameters
     double tau_eff;
     double mu_eff;//kinematic viscosity
-    double nx_gs,ny_gs,nz_gs;//grey-solid color gradient
-    double nx_phase,ny_phase,nz_phase,C_phase;
     double Fx,Fy,Fz;
-	double gp1,gp2,gp4,gp6,gp8,gp9,gp10,gp11,gp12,gp13,gp14,gp15,gp16,gp17,gp18;
-	double gp3,gp5,gp7;
     double Fcpx,Fcpy,Fcpz;//capillary penalty force
 
 	const double mrt_V1=0.05263157894736842;
@@ -2155,15 +2175,15 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 	const double mrt_V12=0.04166666666666666;
 
 	for (n=start; n<finish; n++){
-
 		// read the component number densities
 		nA = Den[n];
 		nB = Den[Np + n];
+
         porosity = Poros[n];
         perm = Perm[n];
-        nx_gs = GreySolidGrad[n+0*Np];
+        W = GreySolidW[n];
-        ny_gs = GreySolidGrad[n+1*Np];
+        Sn_grey = GreySn[n];
-        nz_gs = GreySolidGrad[n+2*Np];
+        Sw_grey = GreySw[n];
 
 		// compute phase indicator field
 		phi=(nA-nB)/(nA+nB);
@@ -2185,100 +2205,63 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 		//........................................................................
 		nn = ijk-1;							// neighbor index (get convention)
 		m1 = Phi[nn];						// get neighbor for phi - 1
-        gp1 = Psi[nn];
 		//........................................................................
 		nn = ijk+1;							// neighbor index (get convention)
 		m2 = Phi[nn];						// get neighbor for phi - 2
-        gp2 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideY;							// neighbor index (get convention)
 		m3 = Phi[nn];					// get neighbor for phi - 3
-        gp3 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideY;							// neighbor index (get convention)
 		m4 = Phi[nn];					// get neighbor for phi - 4
-        gp4 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideZ;						// neighbor index (get convention)
 		m5 = Phi[nn];					// get neighbor for phi - 5
-        gp5 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideZ;						// neighbor index (get convention)
 		m6 = Phi[nn];					// get neighbor for phi - 6
-        gp6 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideY-1;						// neighbor index (get convention)
 		m7 = Phi[nn];					// get neighbor for phi - 7
-        gp7 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideY+1;						// neighbor index (get convention)
 		m8 = Phi[nn];					// get neighbor for phi - 8
-        gp8 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideY-1;						// neighbor index (get convention)
 		m9 = Phi[nn];					// get neighbor for phi - 9
-        gp9 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideY+1;						// neighbor index (get convention)
 		m10 = Phi[nn];					// get neighbor for phi - 10
-        gp10 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideZ-1;						// neighbor index (get convention)
 		m11 = Phi[nn];					// get neighbor for phi - 11
-        gp11 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideZ+1;						// neighbor index (get convention)
 		m12 = Phi[nn];					// get neighbor for phi - 12
-        gp12 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideZ-1;						// neighbor index (get convention)
 		m13 = Phi[nn];					// get neighbor for phi - 13
-        gp13 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideZ+1;						// neighbor index (get convention)
 		m14 = Phi[nn];					// get neighbor for phi - 14
-        gp14 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideZ-strideY;					// neighbor index (get convention)
 		m15 = Phi[nn];					// get neighbor for phi - 15
-        gp15 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideZ+strideY;					// neighbor index (get convention)
 		m16 = Phi[nn];					// get neighbor for phi - 16
-        gp16 = Psi[nn];
 		//........................................................................
 		nn = ijk+strideZ-strideY;					// neighbor index (get convention)
 		m17 = Phi[nn];					// get neighbor for phi - 17
-        gp17 = Psi[nn];
 		//........................................................................
 		nn = ijk-strideZ+strideY;					// neighbor index (get convention)
 		m18 = Phi[nn];					// get neighbor for phi - 18
-        gp18 = Psi[nn];
 		//............Compute the Color Gradient...................................
-		nx_phase = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+		nx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
-		ny_phase = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+		ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
-		nz_phase = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+		nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
-		C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase);
 
-        //correct the normal color gradient by considering the effect of grey solid
+		//............Compute the Greyscale Potential Gradient.....................
-        nx = nx_phase + (1.0-porosity)*nx_gs; 
-        ny = ny_phase + (1.0-porosity)*ny_gs; 
-        nz = nz_phase + (1.0-porosity)*nz_gs; 
-        if (C_phase==0.0){
-            nx = nx_phase; 
-            ny = ny_phase;
-            nz = nz_phase;
-        }
-
-		//...........Normalize the Color Gradient.................................
-		C = sqrt(nx*nx+ny*ny+nz*nz);
-		double ColorMag = C;
-		if (C==0.0) ColorMag=1.0;
-		nx = nx/ColorMag;
-		ny = ny/ColorMag;
-		nz = nz/ColorMag;		
-
-			//............Compute the Greyscale Potential Gradient.....................
 //			Fcpx = 0.0;
 //			Fcpy = 0.0;
 //			Fcpz = 0.0;
@@ -2299,14 +2282,23 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 //                ny = Fcpy/Fcp_mag;
 //                nz = Fcpz/Fcp_mag;
 //            }
-        Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+        Fcpx = nx; 
-        Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+        Fcpy = ny; 
-        Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+        Fcpz = nz; 
+        double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
+        if (Fcp_mag==0.0) Fcpx=Fcpy=Fcpz=0.0;
         //NOTE for open node (porosity=1.0),Fcp=0.0
         Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm); 
         Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm); 
         Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm); 
 
+		//...........Normalize the Color Gradient.................................
+		C = sqrt(nx*nx+ny*ny+nz*nz);
+		double ColorMag = C;
+		if (C==0.0) ColorMag=1.0;
+		nx = nx/ColorMag;
+		ny = ny/ColorMag;
+		nz = nz/ColorMag;		
 
 		// q=0
 		fq = dist[n];
@@ -2608,6 +2600,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 		//........................................................................
 		//..............carry out relaxation process..............................
 		//..........Toelke, Fruediger et. al. 2006................................
+        //---------------- NO higher-order force -------------------------------//
 		if (C == 0.0)	nx = ny = nz = 0.0;
 		m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1);
 		m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2);
@@ -2632,6 +2625,43 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 		m16 = m16 + rlx_setB*( - m16);
 		m17 = m17 + rlx_setB*( - m17);
 		m18 = m18 + rlx_setB*( - m18);
+        //----------------------------------------------------------------------//
+
+        //----------------With higher-order force ------------------------------//
+		//if (C == 0.0)	nx = ny = nz = 0.0;
+		//m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1)
+        //        + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity;
+		//m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2)
+        //        + (1-0.5*rlx_setA)*11*(-Fx*ux-Fy*uy-Fz*uz)/porosity;
+        //jx = jx + Fx;
+		//m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
+        //        + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
+        //jy = jy + Fy;
+		//m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6)
+        //        + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
+        //jz = jz + Fz;
+		//m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
+        //        + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
+		//m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9)
+        //        + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity;
+		////m10 = m10 + rlx_setA*( - m10);
+        //m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10)
+        //          + (1-0.5*rlx_setA)*(-2*Fx*ux+Fy*uy+Fz*uz)/porosity;
+		//m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11)
+        //          + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity;
+		////m12 = m12 + rlx_setA*( - m12);
+        //m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12)
+        //          + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity;
+		//m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13);
+        //          + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity;
+		//m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14);
+        //          + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity;
+		//m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15);
+        //          + (1-0.5*rlx_setA)*(Fz*ux+Fx*uz)/porosity;
+		//m16 = m16 + rlx_setB*( - m16);
+		//m17 = m17 + rlx_setB*( - m17);
+		//m18 = m18 + rlx_setB*( - m18);
+        //----------------------------------------------------------------------//
 
 		//.................inverse transformation......................................................
 		// q=0
@@ -2748,6 +2778,10 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 		// Cq = {1,0,0}, {0,1,0}, {0,0,1}
 		delta = beta*nA*nB*nAB*0.1111111111111111*nx;
 		if (!(nA*nB*nAB>0)) delta=0;
+        //----------------newly added for better control of recoloring---------------//
+        if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+        if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+        //---------------------------------------------------------------------------//
         if (RecoloringOff==true && porosity !=1.0) delta=0;
 		a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
 		b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@@ -2764,6 +2798,10 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 		// Cq = {0,1,0}
 		delta = beta*nA*nB*nAB*0.1111111111111111*ny;
 		if (!(nA*nB*nAB>0)) delta=0;
+        //----------------newly added for better control of recoloring---------------//
+        if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+        if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+        //---------------------------------------------------------------------------//
         if (RecoloringOff==true && porosity !=1.0) delta=0;
 		a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
 		b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@@ -2779,6 +2817,10 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 		// Cq = {0,0,1}
 		delta = beta*nA*nB*nAB*0.1111111111111111*nz;
 		if (!(nA*nB*nAB>0)) delta=0;
+        //----------------newly added for better control of recoloring---------------//
+        if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+        if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+        //---------------------------------------------------------------------------//
         if (RecoloringOff==true && porosity !=1.0) delta=0;
 		a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
 		b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@@ -2790,6 +2832,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 		Aq[6*Np+n] = a2;
 		Bq[6*Np+n] = b2;
 		//...............................................
+
 	}
 }
 

cuda/GreyscaleColor.cu

@@ -1450,7 +1450,7 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist,
 //CP: capillary penalty
 // also turn off recoloring for grey nodes
 __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
-		 double *Phi, double *GreySolidW, double *Poros,double *Perm, double *Velocity, double *Pressure,
+		 double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure,
          double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff,double alpha, double beta,
 		double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
 
@@ -1478,6 +1478,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
     double Fx,Fy,Fz;
     double Fcpx,Fcpy,Fcpz;//capillary penalty force
     double W;//greyscale wetting strength
+    double Sn_grey,Sw_grey;
 
 	const double mrt_V1=0.05263157894736842;
 	const double mrt_V2=0.012531328320802;
@@ -1504,6 +1505,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
             porosity = Poros[n];
             perm = Perm[n];
             W = GreySolidW[n];
+            Sn_grey = GreySn[n];
+            Sw_grey = GreySw[n];
 
 			// compute phase indicator field
 			phi=(nA-nB)/(nA+nB);
@@ -2165,6 +2168,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 			// Cq = {1,0,0}, {0,1,0}, {0,0,1}
 			delta = beta*nA*nB*nAB*0.1111111111111111*nx;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@@ -2184,6 +2191,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 			// Cq = {0,1,0}
 			delta = beta*nA*nB*nAB*0.1111111111111111*ny;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@@ -2204,6 +2215,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 			// Cq = {0,0,1}
 			delta = beta*nA*nB*nAB*0.1111111111111111*nz;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@@ -2226,7 +2241,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 //CP: capillary penalty
 // also turn off recoloring for grey nodes
 __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, 
-        double *Phi, double *GreySolidW, double *Poros,double *Perm, double *Velocity, double *Pressure, 
+        double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure, 
         double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
 		double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
 	int ijk,nn,n;
@@ -2249,6 +2264,7 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
     double Fx,Fy,Fz;
     double Fcpx,Fcpy,Fcpz;//capillary penalty force
     double W;//greyscale wetting strength
+    double Sn_grey,Sw_grey;
 
 	const double mrt_V1=0.05263157894736842;
 	const double mrt_V2=0.012531328320802;
@@ -2276,6 +2292,8 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
             porosity = Poros[n];
             perm = Perm[n];
             W = GreySolidW[n];
+            Sn_grey = GreySn[n];
+            Sw_grey = GreySw[n];
 
 			// compute phase indicator field
 			phi=(nA-nB)/(nA+nB);
@@ -2870,6 +2888,10 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
 			// Cq = {1,0,0}, {0,1,0}, {0,0,1}
 			delta = beta*nA*nB*nAB*0.1111111111111111*nx;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@@ -2886,6 +2908,10 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
 			// Cq = {0,1,0}
 			delta = beta*nA*nB*nAB*0.1111111111111111*ny;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@@ -2901,6 +2927,10 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
 			// Cq = {0,0,1}
 			delta = beta*nA*nB*nAB*0.1111111111111111*nz;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@@ -4500,12 +4530,13 @@ extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, doubl
 
 //Model-1 & 4 with capillary pressure penalty
 extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, 
-        double *Phi, double *GreySolidW, double *Poros,double *Perm,double *Vel, double *Pressure,
+        double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure,
         double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
 		double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
 
-	dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, Poros, Perm, Vel, Pressure,
+	dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm, Vel, Pressure,
             rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np);
+	
 	cudaError_t err = cudaGetLastError();
 	if (cudaSuccess != err){
 		printf("CUDA error in ScaLBL_D3Q19_AAeven_GreyscaleColor_CP: %s \n",cudaGetErrorString(err));
@@ -4515,11 +4546,11 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 
 //Model-1 & 4 with capillary pressure penalty
 extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, 
-		double *Phi, double *GreySolidW, double *Poros,double *Perm,double *Vel,double *Pressure, 
+		double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure, 
         double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
 		double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
 
-	dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, Poros, Perm,Vel,Pressure,
+	dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm,Vel,Pressure,
 			rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np);
 
 	cudaError_t err = cudaGetLastError();

docs/source/install.rst

@@ -2,7 +2,161 @@
 Installation
 ============
 
-LBPM requires CMake, MPI and HDF5 as required dependencies. 
+Sample scripts to configure and build LBPM are included in the `sample_scripts` directory. To build this package, MPI and cmake are required, along with a compiler that supports the C++-14 standard. Building in source is not supported.
+
+The essential dependencies needed to build LBPM are:
+
+1.  cmake (version 3.9 or higher)
+2.  MPI
+3.  HDF5
+4.  silo
+
+*************************
+Building Dependencies 
+*************************
+(skip if they are already installed on your system)
+
+1. Download third-party library source files
+
+.. code-block:: bash
+
+   zlib-1.2.11.tar.gz
+   hdf5-1.8.12.tar.gz 
+   silo-4.10.2.tar.gz 
+
+
+2. Set up path where you want to install 
+
+.. code-block:: bash
+
+   export MPI_DIR=/path/to/mpi/
+   export LBPM_ZLIB_DIR=/path/to/zlib
+   export LBPM_HDF5_DIR=/path/to/hdf5
+   export LBPM_SILO_DIR=/path/to/silo
+
+3. Build third-party library dependencies 
+
+.. code-block:: bash
+
+   tar -xzvf zlib-1.2.11.tar.gz 
+   tar -xzvf hdf5-1.8.12.tar.gz 
+   tar -xzvf silo-4.10.2.tar.gz 
+   cd zlib-1.2.11
+
+   ./configure --prefix=$LBPM_ZLIB_DIR && make && make install
+
+   cd ../hdf5-1.8.12
+
+   CC=$MPI_DIR/bin/mpicc  CXX=$MPI_DIR/bin/mpicxx CXXFLAGS="-fPIC -O3 -std=c++14" \
+   ./configure --prefix=$LBPM_HDF5_DIR --enable-parallel --enable-shared --with-zlib=$LBPM_ZLIB_DIR 
+    make && make install
+
+
+   cd ../silo-4.10.2
+
+   CC=$MPI_DIR/bin/mpicc  CXX=$MPI_DIR/bin/mpicxx CXXFLAGS="-fPIC -O3 -std=c++14" \
+   ./configure --prefix=$LBPM_SILO_DIR -with-hdf5=$LBPM_HDF5_DIR/include,$LBPM_HDF5_DIR/lib --enable-static 
+    make && make install
+ 
+*************************
+Building LBPM 
+*************************
+
+Many HPC systems will include all of these dependencies, and LBPM can be built simply by setting the paths to the associated libraries in the cmake configure line. 
+
+The steps typically used to build LBPM are as follows:
+
+1. Set environment variables for the source directory `$LBPM_WIA_SOURCE` and build directory `$LBPM_WIA_DIR`  
+
+.. code-block:: bash
+
+   export LBPM_SOURCE=/path/to/source/LBPM  
+   export LBPM_DIR=/path/to/build/LBPM
+
+2. Set environment variables for the path to HDF5 and SILO (required),  and optionally for TimerUtility and NetCDF (optional)  
+
+.. code-block:: bash
+
+   export LBPM_HDF5_DIR=/path/to/hdf5  
+   export LBPM_SILO_DIR=/path/to/silo 
+   export LBPM_TIMER_DIR=/path/to/timer 
+   export LBPM_NETCDF_DIR=/path/to/netcdf
+
+3. Create the build directory and navigate to it  
+
+.. code-block:: bash
+
+   mkdir $LBPM_WIA_DIR  
+   cd $LBPM_WIA_DIR
+
+4. Configure the project. Numerous scripts exist to build LBPM on different HPC clusters, which are available in the  `$LBPM_SOURCE/sample_scripts/` directory. It is often possible to run these scripts directly if one exists for a system similar to the one you are building on. For a standard CPU build:  
+
+.. code-block:: bash
+
+    cmake                                           \  
+        -D CMAKE_BUILD_TYPE:STRING=Release          \   
+        -D CMAKE_C_COMPILER:PATH=mpicc              \   
+        -D CMAKE_CXX_COMPILER:PATH=mpicxx           \  
+        -D CMAKE_C_FLAGS="-fPIC"                    \  
+        -D CMAKE_CXX_FLAGS="-fPIC"                  \  
+        -D CMAKE_CXX_STD=14                         \  
+        -D USE_TIMER=0                              \  
+            -D TIMER_DIRECTORY=$LBPM_TIMER_DIR     \  
+        -D USE_NETCDF=0                             \  
+            -D NETCDF_DIRECTORY=$LBPM_NETCDF_DIR   \  
+        -D USE_SILO=1                               \  
+           -D HDF5_DIRECTORY=$LBPM_HDF5_DIR         \  
+           -D SILO_DIRECTORY=$LBPM_SILO_DIR         \  
+        -D USE_CUDA=0                               \
+        $LBPM_SOURCE  
+
+For GPU support, it is necessary to have CUDA along with a GPU-aware MPI implementation. Otherwise, the LBPM routines should behave identically irrespective of the underlying hardware.  
+
+.. code-block:: bash
+
+    cmake                                           \  
+        -D CMAKE_BUILD_TYPE:STRING=Release          \   
+        -D CMAKE_C_COMPILER:PATH=mpicc              \   
+        -D CMAKE_CXX_COMPILER:PATH=mpicxx           \  
+        -D CMAKE_C_FLAGS="-fPIC"                    \  
+        -D CMAKE_CXX_FLAGS="-fPIC"                  \  
+        -D CMAKE_CXX_STD=14                         \  
+        -D USE_TIMER=0                              \  
+            -D TIMER_DIRECTORY=$LBPM_TIMER_DIR     \  
+        -D USE_NETCDF=0                             \  
+            -D NETCDF_DIRECTORY=$LBPM_NETCDF_DIR   \  
+        -D USE_SILO=1                               \  
+           -D HDF5_DIRECTORY=$LBPM_HDF5_DIR         \  
+           -D SILO_DIRECTORY=$LBPM_SILO_DIR         \  
+        -D USE_CUDA=1                               \  
+        -D CMAKE_CUDA_FLAGS="-arch sm_70"           \  
+        $LBPM_SOURCE   
+
+5. Build the project (using four cores to build)   
+
+.. code-block:: bash
+  
+   make -j4
+
+6. Install the project  
+
+.. code-block:: bash
+
+   make install
+
+7. Run the tests to make sure they execute correctly (on a cluster, it is recommended to run these using the batch system rather than on the head node)  
+
+.. code-block:: bash
+
+   ctest
+
+
+*************************
+Sample Scripts
+*************************
+
+The LBPM repository contains sample scripts showing successful CMake configuration, build and
+install steps for a range of systems. Refer to the project sub-directory below for these examples.
 
 .. code-block:: bash
 

docs/source/userGuide/models/color/protocols/centrifuge.rst

@@ -3,10 +3,41 @@ Color model -- Centrifuge Protocol
 ======================================
 
 The centrifuge protocol is designed to mimic SCAL centrifuge experiments that
-are used to infer the capillary pressure. The centrifuge protocol 
+are used to infer the capillary pressure. The LBPM centrifuge protocol is
+constructed as an unsteady simulation with constant pressure boundary conditions
+and zero pressure drop across the sample. This will enforce the following key values
+
+* ``BC = 3`` -- constant pressure boundary condition
+* ``din = 1.0`` -- inlet pressure value 
+* ``dout = 1.0`` -- outlet pressure value
+
+By default LBPM will populate the inlet reservoir with fluid A (usually the non-wetting fluid)
+and the outlet reservoir with fluid B (usually water). Flow is induced by setting an external
+body force to generate displacement in the ``z`` direction. If the body force is set to
+zero, e.g. ``F = 0, 0, 0``, the simulation will produce spontaneous imbibition, with the
+balancing point being determined based on zero pressure drop across the sample. Setting
+an external body force will shift the capillary pressure. Setting a positive force will
+cause fluid A to be forced into the sample. Once steady conditions are achieved,
+the pressure of fluid A will be larger than fluid B. Alternatively, if the driving force is
+negative then fluid B will be forced into the sample, and the steady-state configuration
+will stabilize to a configuration where fluid B has a larger pressure compared to fluid A.
+The capillary pressure is thereby inferred based on the body force.
+
+In a conventional SCAL experiment the centrifugal forces are proportional to the density
+difference between fluids. While this is still true for LBPM simulation, the body force will
+still be effective even if there is no difference in density between the fluids.
+This is because a positive body force will favor a larger saturation of fluid A
+(positive capillary pressure ) whereas a negative body force will favor a lower
+saturation of fluid A (negative capillary pressure). 
+
+
+To enable the ``centrifuge`` protocol such that the pressure of fluid B is higher than
+fluid A, the following keys can be set. Increasing the body force will lead to a larger
+capillary pressure
 
 .. code-block:: bash
-
+		
-    image_sequence = "image1.raw", "image2.raw"
+    protocol = "centrifuge"
+    F = 0, 0, -1.0e-5
     
 

hip/GreyscaleColor.cu

@@ -1450,9 +1450,9 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist,
 //CP: capillary penalty
 // also turn off recoloring for grey nodes
 __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
-		 double *Phi, double *Psi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity, double *Pressure,
+		 double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure,
          double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff,double alpha, double beta,
-		double Gx, double Gy, double Gz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){
+		double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
 
 	int n,nn,ijk,nread;
 	int nr1,nr2,nr3,nr4,nr5,nr6;
@@ -1462,8 +1462,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 	double fq;
 	// conserved momemnts
 	double rho,jx,jy,jz;
-	//double vx,vy,vz,v_mag;
-    //double ux,uy,uz,u_mag;
     double ux,uy,uz;
 	// non-conserved moments
 	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
@@ -1473,18 +1471,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 	double C,nx,ny,nz; //color gradient magnitude and direction
 	double phi,tau,rho0,rlx_setA,rlx_setB;
 
-    //double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002)
     double porosity;
     double perm;//voxel permeability
-    //double c0, c1; //Guo's model parameters
     double tau_eff;
     double mu_eff;//kinematic viscosity
-    double nx_gs,ny_gs,nz_gs;//grey-solid color gradient
-    double nx_phase,ny_phase,nz_phase,C_phase;
     double Fx,Fy,Fz;
-	double gp1,gp2,gp4,gp6,gp8,gp9,gp10,gp11,gp12,gp13,gp14,gp15,gp16,gp17,gp18;
-	double gp3,gp5,gp7;
     double Fcpx,Fcpy,Fcpz;//capillary penalty force
+    double W;//greyscale wetting strength
+    double Sn_grey,Sw_grey;
 
 	const double mrt_V1=0.05263157894736842;
 	const double mrt_V2=0.012531328320802;
@@ -1510,9 +1504,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 
             porosity = Poros[n];
             perm = Perm[n];
-            nx_gs = GreySolidGrad[n+0*Np];
+            W = GreySolidW[n];
-            ny_gs = GreySolidGrad[n+1*Np];
+            Sn_grey = GreySn[n];
-            nz_gs = GreySolidGrad[n+2*Np];
+            Sw_grey = GreySw[n];
 
 			// compute phase indicator field
 			phi=(nA-nB)/(nA+nB);
@@ -1534,98 +1528,61 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 			//........................................................................
 			nn = ijk-1;							// neighbor index (get convention)
 			m1 = Phi[nn];						// get neighbor for phi - 1
-            gp1 = Psi[nn];
 			//........................................................................
 			nn = ijk+1;							// neighbor index (get convention)
 			m2 = Phi[nn];						// get neighbor for phi - 2
-            gp2 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideY;							// neighbor index (get convention)
 			m3 = Phi[nn];					// get neighbor for phi - 3
-            gp3 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideY;							// neighbor index (get convention)
 			m4 = Phi[nn];					// get neighbor for phi - 4
-            gp4 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideZ;						// neighbor index (get convention)
 			m5 = Phi[nn];					// get neighbor for phi - 5
-            gp5 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideZ;						// neighbor index (get convention)
 			m6 = Phi[nn];					// get neighbor for phi - 6
-            gp6 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideY-1;						// neighbor index (get convention)
 			m7 = Phi[nn];					// get neighbor for phi - 7
-            gp7 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideY+1;						// neighbor index (get convention)
 			m8 = Phi[nn];					// get neighbor for phi - 8
-            gp8 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideY-1;						// neighbor index (get convention)
 			m9 = Phi[nn];					// get neighbor for phi - 9
-            gp9 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideY+1;						// neighbor index (get convention)
 			m10 = Phi[nn];					// get neighbor for phi - 10
-            gp10 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideZ-1;						// neighbor index (get convention)
 			m11 = Phi[nn];					// get neighbor for phi - 11
-            gp11 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideZ+1;						// neighbor index (get convention)
 			m12 = Phi[nn];					// get neighbor for phi - 12
-            gp12 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideZ-1;						// neighbor index (get convention)
 			m13 = Phi[nn];					// get neighbor for phi - 13
-            gp13 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideZ+1;						// neighbor index (get convention)
 			m14 = Phi[nn];					// get neighbor for phi - 14
-            gp14 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideZ-strideY;					// neighbor index (get convention)
 			m15 = Phi[nn];					// get neighbor for phi - 15
-            gp15 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideZ+strideY;					// neighbor index (get convention)
 			m16 = Phi[nn];					// get neighbor for phi - 16
-            gp16 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideZ-strideY;					// neighbor index (get convention)
 			m17 = Phi[nn];					// get neighbor for phi - 17
-            gp17 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideZ+strideY;					// neighbor index (get convention)
 			m18 = Phi[nn];					// get neighbor for phi - 18
-            gp18 = Psi[nn];
 			//............Compute the Color Gradient...................................
-			nx_phase = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+			nx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
-			ny_phase = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+			ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
-			nz_phase = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+			nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
-			C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase);
-
-            //correct the normal color gradient by considering the effect of grey solid
-            nx = nx_phase + (1.0-porosity)*nx_gs; 
-            ny = ny_phase + (1.0-porosity)*ny_gs; 
-            nz = nz_phase + (1.0-porosity)*nz_gs; 
-            if (C_phase==0.0){
-                nx = nx_phase; 
-                ny = ny_phase;
-                nz = nz_phase;
-            }
-
-			//...........Normalize the Color Gradient.................................
-			C = sqrt(nx*nx+ny*ny+nz*nz);
-			double ColorMag = C;
-			if (C==0.0) ColorMag=1.0;
-			nx = nx/ColorMag;
-			ny = ny/ColorMag;
-			nz = nz/ColorMag;		
 
 			//............Compute the Greyscale Potential Gradient.....................
 //			Fcpx = 0.0;
@@ -1648,14 +1605,24 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 //                //ny = Fcpy/Fcp_mag;
 //                //nz = Fcpz/Fcp_mag;
 //            }
-            Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+            Fcpx = nx;
-            Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+            Fcpy = ny;
-            Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+            Fcpz = nz;
+            double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
+            if (Fcp_mag==0.0); Fcpx=Fcpy=Fcpz=0.0;
             //NOTE for open node (porosity=1.0),Fcp=0.0
             Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm); 
             Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm); 
             Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm); 
 
+			//...........Normalize the Color Gradient.................................
+			C = sqrt(nx*nx+ny*ny+nz*nz);
+			double ColorMag = C;
+			if (C==0.0) ColorMag=1.0;
+			nx = nx/ColorMag;
+			ny = ny/ColorMag;
+			nz = nz/ColorMag;		
+
 			// q=0
 			fq = dist[n];
 			rho = fq;
@@ -2201,6 +2168,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 			// Cq = {1,0,0}, {0,1,0}, {0,0,1}
 			delta = beta*nA*nB*nAB*0.1111111111111111*nx;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@@ -2220,6 +2191,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 			// Cq = {0,1,0}
 			delta = beta*nA*nB*nAB*0.1111111111111111*ny;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@@ -2240,6 +2215,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 			// Cq = {0,0,1}
 			delta = beta*nA*nB*nAB*0.1111111111111111*nz;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@@ -2262,15 +2241,13 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
 //CP: capillary penalty
 // also turn off recoloring for grey nodes
 __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, 
-        double *Phi, double *Psi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity, double *Pressure, 
+        double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure, 
         double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
-		double Gx, double Gy, double Gz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){
+		double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
 	int ijk,nn,n;
 	double fq;
 	// conserved momemnts
 	double rho,jx,jy,jz;
-	//double vx,vy,vz,v_mag;
-    //double ux,uy,uz,u_mag;
     double ux,uy,uz;
 	// non-conserved moments
 	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
@@ -2280,18 +2257,14 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
 	double C,nx,ny,nz; //color gradient magnitude and direction
 	double phi,tau,rho0,rlx_setA,rlx_setB;
 
-    //double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002)
     double porosity;
     double perm;//voxel permeability
-    //double c0, c1; //Guo's model parameters
     double tau_eff;
     double mu_eff;//kinematic viscosity
-    double nx_gs,ny_gs,nz_gs;//grey-solid color gradient
-    double nx_phase,ny_phase,nz_phase,C_phase;
     double Fx,Fy,Fz;
-	double gp1,gp2,gp4,gp6,gp8,gp9,gp10,gp11,gp12,gp13,gp14,gp15,gp16,gp17,gp18;
-	double gp3,gp5,gp7;
     double Fcpx,Fcpy,Fcpz;//capillary penalty force
+    double W;//greyscale wetting strength
+    double Sn_grey,Sw_grey;
 
 	const double mrt_V1=0.05263157894736842;
 	const double mrt_V2=0.012531328320802;
@@ -2315,11 +2288,12 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
 			// read the component number densities
 			nA = Den[n];
 			nB = Den[Np + n];
+
             porosity = Poros[n];
             perm = Perm[n];
-            nx_gs = GreySolidGrad[n+0*Np];
+            W = GreySolidW[n];
-            ny_gs = GreySolidGrad[n+1*Np];
+            Sn_grey = GreySn[n];
-            nz_gs = GreySolidGrad[n+2*Np];
+            Sw_grey = GreySw[n];
 
 			// compute phase indicator field
 			phi=(nA-nB)/(nA+nB);
@@ -2341,98 +2315,61 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
 			//........................................................................
 			nn = ijk-1;							// neighbor index (get convention)
 			m1 = Phi[nn];						// get neighbor for phi - 1
-            gp1 = Psi[nn];
 			//........................................................................
 			nn = ijk+1;							// neighbor index (get convention)
 			m2 = Phi[nn];						// get neighbor for phi - 2
-            gp2 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideY;							// neighbor index (get convention)
 			m3 = Phi[nn];					// get neighbor for phi - 3
-            gp3 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideY;							// neighbor index (get convention)
 			m4 = Phi[nn];					// get neighbor for phi - 4
-            gp4 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideZ;						// neighbor index (get convention)
 			m5 = Phi[nn];					// get neighbor for phi - 5
-            gp5 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideZ;						// neighbor index (get convention)
 			m6 = Phi[nn];					// get neighbor for phi - 6
-            gp6 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideY-1;						// neighbor index (get convention)
 			m7 = Phi[nn];					// get neighbor for phi - 7
-            gp7 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideY+1;						// neighbor index (get convention)
 			m8 = Phi[nn];					// get neighbor for phi - 8
-            gp8 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideY-1;						// neighbor index (get convention)
 			m9 = Phi[nn];					// get neighbor for phi - 9
-            gp9 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideY+1;						// neighbor index (get convention)
 			m10 = Phi[nn];					// get neighbor for phi - 10
-            gp10 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideZ-1;						// neighbor index (get convention)
 			m11 = Phi[nn];					// get neighbor for phi - 11
-            gp11 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideZ+1;						// neighbor index (get convention)
 			m12 = Phi[nn];					// get neighbor for phi - 12
-            gp12 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideZ-1;						// neighbor index (get convention)
 			m13 = Phi[nn];					// get neighbor for phi - 13
-            gp13 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideZ+1;						// neighbor index (get convention)
 			m14 = Phi[nn];					// get neighbor for phi - 14
-            gp14 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideZ-strideY;					// neighbor index (get convention)
 			m15 = Phi[nn];					// get neighbor for phi - 15
-            gp15 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideZ+strideY;					// neighbor index (get convention)
 			m16 = Phi[nn];					// get neighbor for phi - 16
-            gp16 = Psi[nn];
 			//........................................................................
 			nn = ijk+strideZ-strideY;					// neighbor index (get convention)
 			m17 = Phi[nn];					// get neighbor for phi - 17
-            gp17 = Psi[nn];
 			//........................................................................
 			nn = ijk-strideZ+strideY;					// neighbor index (get convention)
 			m18 = Phi[nn];					// get neighbor for phi - 18
-            gp18 = Psi[nn];
 			//............Compute the Color Gradient...................................
-			nx_phase = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+			nx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
-			ny_phase = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+			ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
-			nz_phase = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+			nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
-			C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase);
-
-            //correct the normal color gradient by considering the effect of grey solid
-            nx = nx_phase + (1.0-porosity)*nx_gs; 
-            ny = ny_phase + (1.0-porosity)*ny_gs; 
-            nz = nz_phase + (1.0-porosity)*nz_gs; 
-            if (C_phase==0.0){
-                nx = nx_phase; 
-                ny = ny_phase;
-                nz = nz_phase;
-            }
-
-			//...........Normalize the Color Gradient.................................
-			C = sqrt(nx*nx+ny*ny+nz*nz);
-			double ColorMag = C;
-			if (C==0.0) ColorMag=1.0;
-			nx = nx/ColorMag;
-			ny = ny/ColorMag;
-			nz = nz/ColorMag;		
 
 			//............Compute the Greyscale Potential Gradient.....................
 //			Fcpx = 0.0;
@@ -2455,14 +2392,23 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
 //                ny = Fcpy/Fcp_mag;
 //                nz = Fcpz/Fcp_mag;
 //            }
-            Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+            Fcpx = nx; 
-            Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+            Fcpy = ny; 
-            Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+            Fcpz = nz; 
+            double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
+            if (Fcp_mag==0.0); Fcpx=Fcpy=Fcpz=0.0;
             //NOTE for open node (porosity=1.0),Fcp=0.0
             Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm); 
             Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm); 
             Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm); 
 
+			//...........Normalize the Color Gradient.................................
+			C = sqrt(nx*nx+ny*ny+nz*nz);
+			double ColorMag = C;
+			if (C==0.0) ColorMag=1.0;
+			nx = nx/ColorMag;
+			ny = ny/ColorMag;
+			nz = nz/ColorMag;		
 
 			// q=0
 			fq = dist[n];
@@ -2942,6 +2888,10 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
 			// Cq = {1,0,0}, {0,1,0}, {0,0,1}
 			delta = beta*nA*nB*nAB*0.1111111111111111*nx;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@@ -2958,6 +2908,10 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
 			// Cq = {0,1,0}
 			delta = beta*nA*nB*nAB*0.1111111111111111*ny;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@@ -2973,6 +2927,10 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
 			// Cq = {0,0,1}
 			delta = beta*nA*nB*nAB*0.1111111111111111*nz;
 			if (!(nA*nB*nAB>0)) delta=0;
+            //----------------newly added for better control of recoloring---------------//
+            if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; 
+            if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; 
+            //---------------------------------------------------------------------------//
             if (RecoloringOff==true && porosity !=1.0) delta=0;
 			a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
 			b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@@ -2989,6 +2947,7 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
 	}
 }
 
+
 __global__ void dvc_ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, double *Bq, int start, int finish, int Np){
 	int idx;
 	double nA,nB;
@@ -4572,12 +4531,12 @@ extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, doubl
 
 //Model-1 & 4 with capillary pressure penalty
 extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, 
-        double *Phi,double *Psi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel, double *Pressure,
+        double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure,
         double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
-		double Fx, double Fy, double Fz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){
+		double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
 
-	dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, Psi, GreySolidGrad, Poros, Perm, Vel, Pressure,
+	dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm, Vel, Pressure,
-            rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, W, strideY, strideZ, start, finish, Np);
+            rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np);
 	hipError_t err = hipGetLastError();
 	if (hipSuccess != err){
 		printf("hip error in ScaLBL_D3Q19_AAeven_GreyscaleColor_CP: %s \n",hipGetErrorString(err));
@@ -4587,12 +4546,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
 
 //Model-1 & 4 with capillary pressure penalty
 extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, 
-		double *Phi, double *Psi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel,double *Pressure, 
+		double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure, 
         double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
-		double Fx, double Fy, double Fz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){
+		double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
 
-	dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, Psi, GreySolidGrad, Poros, Perm,Vel,Pressure,
+	dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm,Vel,Pressure,
-			rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, W, strideY, strideZ, start, finish, Np);
+			rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np);
 
 	hipError_t err = hipGetLastError();
 	if (hipSuccess != err){
@@ -4600,51 +4559,3 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *M
 	}
 }
 
-//extern "C" void ScaLBL_Update_GreyscalePotential(int *Map, double *Phi, double *Psi, double *Poro, double *Perm, double alpha, double W, 
-//		int start, int finish, int Np){
-//
-//	dvc_ScaLBL_Update_GreyscalePotential<<<NBLOCKS,NTHREADS >>>(Map, Phi, Psi, Poro, Perm, alpha, W, start, finish, Np);
-//
-//	hipError_t err = hipGetLastError();
-//	if (hipSuccess != err){
-//		printf("hip error in ScaLBL_Update_GreyscalePotential: %s \n",hipGetErrorString(err));
-//	}
-//}
-
-////Model-2&3
-//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den, 
-//        double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel, 
-//        double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
-//		double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){
-//
-//	//cudaProfilerStart();
-//	//cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor, cudaFuncCachePreferL1);
-//
-//	dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidGrad, Poros, Perm, Vel, 
-//            rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np);
-//	hipError_t err = hipGetLastError();
-//	if (hipSuccess != err){
-//		printf("hip error in ScaLBL_D3Q19_AAeven_GreyscaleColor: %s \n",hipGetErrorString(err));
-//	}
-//	//cudaProfilerStop();
-//
-//}
-//
-////Model-2&3
-//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, 
-//		double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel, 
-//        double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
-//		double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){
-//
-//	//cudaProfilerStart();
-//	//cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor, cudaFuncCachePreferL1);
-//	
-//	dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi,  GreySolidGrad, Poros, Perm,Vel, 
-//			rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np);
-//
-//	hipError_t err = hipGetLastError();
-//	if (hipSuccess != err){
-//		printf("hip error in ScaLBL_D3Q19_AAodd_GreyscaleColor: %s \n",hipGetErrorString(err));
-//	}
-//	//cudaProfilerStop();
-//}

models/ColorModel.cpp

@@ -219,8 +219,16 @@ void ScaLBL_ColorModel::ReadParams(string filename){
 			if (rank==0) printf("WARNING: protocol (core flooding) supports only volumetric flux boundary condition \n");
 		}
 		domain_db->putScalar<int>( "BC", BoundaryCondition );
+		if (color_db->keyExists( "capillary_number" )){
+			double capillary_number = color_db->getScalar<double>( "capillary_number" );
+			double MuB = rhoB*(tauB - 0.5)/3.0;
+			double IFT = 6.0*alpha;
+			double CrossSectionalArea = (double) (nprocx*(Nx-2)*nprocy*(Ny-2));
+			flux = Dm->Porosity()*CrossSectionalArea*IFT*capillary_number/MuB;
+			if (rank==0) printf("  protocol (core flooding): set flux=%f to achieve Ca=%f \n",flux, capillary_number);
+		}
+		color_db->putScalar<double>( "flux", flux );
 	} 
-	
 }
 
 void ScaLBL_ColorModel::SetDomain(){
@@ -739,7 +747,6 @@ double ScaLBL_ColorModel::Run(int returntime){
 			double flow_rate_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_z*dir_z);
 			double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(5.796*alpha);
 
-
 			bool isSteady = false;
 			if ( (fabs((Ca - Ca_previous)/Ca) < tolerance &&  CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS))
 				isSteady = true;
@@ -800,7 +807,16 @@ double ScaLBL_ColorModel::Run(int returntime){
 					double vBd_x = Averages->gwd.Px/Mass_w; 
 					double vBd_y = Averages->gwd.Py/Mass_w; 
 					double vBd_z = Averages->gwd.Pz/Mass_w;
-
+					// film contribution
+					double Mfn = Averages->gifs.Mn; 
+					double Mfw = Averages->gifs.Mw; 
+					double vfn_x = Averages->gifs.Pnx; 
+					double vfn_y = Averages->gifs.Pny; 
+					double vfn_z = Averages->gifs.Pnz; 
+					double vfw_x = Averages->gifs.Pwx; 
+					double vfw_y = Averages->gifs.Pwy; 
+					double vfw_z = Averages->gifs.Pwz;
+					
 					double flow_rate_A_connected = Vol_nc*(vAc_x*dir_x + vAc_y*dir_y + vAc_z*dir_z);
 					double flow_rate_B_connected = Vol_wc*(vBc_x*dir_x + vBc_y*dir_y + vBc_z*dir_z);
 					double flow_rate_A_disconnected = (Vol_nd)*(vAd_x*dir_x + vAd_y*dir_y + vAd_z*dir_z);
@@ -815,6 +831,19 @@ double ScaLBL_ColorModel::Run(int returntime){
 					double kAeff = h*h*muA*(flow_rate_A)/(force_mag);
 					double kBeff = h*h*muB*(flow_rate_B)/(force_mag);
 
+					/* flow rate = [volume of fluid] x [momentum of fluid] / [mass of fluid] */
+					/* fluid A eats the films */
+					double flow_rate_A_filmA = (flow_rate_A*Averages->gnb.M + volA*(vfn_x*dir_x + vfn_y*dir_y + vfn_z*dir_z))/(Averages->gnb.M + Mfn);
+					double flow_rate_B_filmA = (flow_rate_B*Averages->gwb.M - volB*(vfn_x*dir_x + vfn_y*dir_y + vfn_z*dir_z))/(Averages->gwb.M - (rhoB/rhoA)*Mfn);
+					/* fluid B eats the films */
+					double flow_rate_A_filmB = (flow_rate_A*Averages->gnb.M - volA*(vfw_x*dir_x + vfw_y*dir_y + vfw_z*dir_z))/(Averages->gnb.M - (rhoA/rhoB)*Mfw);
+					double flow_rate_B_filmB = (flow_rate_B*Averages->gwb.M + volB*(vfw_x*dir_x + vfw_y*dir_y + vfw_z*dir_z))/(Averages->gwb.M + Mfw);
+					/* effective permeability uncertainty limits */
+					double kAeff_filmA = h*h*muA*(flow_rate_A_filmA)/(force_mag);
+					double kBeff_filmA = h*h*muB*(flow_rate_B_filmA)/(force_mag);
+					double kAeff_filmB = h*h*muA*(flow_rate_A_filmB)/(force_mag);
+					double kBeff_filmB = h*h*muB*(flow_rate_B_filmB)/(force_mag);
+					
 					double viscous_pressure_drop = (rhoA*volA + rhoB*volB)*force_mag;
 					double Mobility = muA/muB;
 
@@ -825,10 +854,14 @@ double ScaLBL_ColorModel::Run(int returntime){
 					else
 						WriteHeader=true;
 					kr_log_file = fopen("relperm.csv","a");
-					if (WriteHeader)
+					if (WriteHeader){
-						fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water eff.perm.oil.connected eff.perm.water.connected eff.perm.oil.disconnected eff.perm.water.disconnected cap.pressure cap.pressure.connected pressure.drop Ca M\n");
+						fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water eff.perm.oil.connected eff.perm.water.connected eff.perm.oil.disconnected eff.perm.water.disconnected ");
-
+						fprintf(kr_log_file,"eff.perm.oil.upper.bound eff.perm.water.lower.bound eff.perm.oil.lower.bound eff.perm.water.upper.bound ");
-					fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_TIMESTEP,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,kAeff_disconnected,kBeff_disconnected,pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility);
+						fprintf(kr_log_file,"cap.pressure cap.pressure.connected pressure.drop Ca M\n");
+					}
+					fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",CURRENT_TIMESTEP,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,kAeff_disconnected,kBeff_disconnected);
+					fprintf(kr_log_file,"%.5g %.5g %.5g %.5g ",kAeff_filmA, kBeff_filmA, kAeff_filmB,kBeff_filmB);
+					fprintf(kr_log_file,"%.5g %.5g %.5g %.5g %.5g\n",pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility);
 					fclose(kr_log_file);
 
 					printf("  Measured capillary number %f \n ",Ca);
@@ -2086,12 +2119,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 	ScaLBL_CopyToHost(Vel_y, &M.Velocity[Np], Np*sizeof(double));
 	ScaLBL_CopyToHost(Vel_z, &M.Velocity[2*Np], Np*sizeof(double));
 
-	/* DEBUG STRUCTURES */
 	int Nx = M.Nx;  int Ny = M.Ny;  int Nz = M.Nz;
-	//int N = Nx*Ny*Nz;
-	//double * DebugMassA, *DebugMassB;
-	//DebugMassA = new double[Np];
-	//DebugMassB = new double[Np];
 
 	/* compute the total momentum */
 	vax = vay = vaz = 0.0;
@@ -2187,111 +2215,9 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 			}
 		}
 	}
-	/*	   for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){	    
-		   phi = Phase[n];
-		   vx = Vel_x[n];
-		   vy = Vel_y[n];
-		   vz = Vel_z[n];
-		   double local_momentum = sqrt(vx*vx+vy*vy+vz*vz);
-		   if (phi > 0.0){
-			   LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_A;
-			   Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+3*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE;
-		   }
-		   else{
-			   LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B;
-			   Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+3*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   DebugMassB[n] = LOCAL_MASS_CHANGE;
-		   }
-	   }
 
-	  for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){
-		   phi = Phase[n];
-		   vx = Vel_x[n];
-		   vy = Vel_y[n];
-		   vz = Vel_z[n];
-		   double local_momentum = sqrt(vx*vx+vy*vy+vz*vz);
-		   if (phi > 0.0){
-			   LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_A;
-			   Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+3*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE;
-		   }
-		   else{
-			   LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B;
-			   Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+3*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
-			   DebugMassB[n] = LOCAL_MASS_CHANGE;
-		   }
-	  }
-	 */
 	if (M.rank == 0) printf("Update Fractional Flow: change mass of fluid B by %f \n",TOTAL_MASS_CHANGE/mass_b_global);
 
-	/* Print out debugging info with mass update  
-	  // initialize the array
-	  double value;
-	  char LocalRankFilename[40];
-	  DoubleArray regdata(Nx,Ny,Nz);
-	  regdata.fill(0.f);
-	  for (int k=0; k<Nz; k++){
-		  for (int j=0; j<Ny; j++){
-			  for (int i=0; i<Nx; i++){
-				  int idx=M.Map(i,j,k);
-				  if (!(idx<0)){
-					  value=DebugMassA[idx];
-					  regdata(i,j,k)=value;
-				  }
-			  }
-		  }
-	  }
-
-	  FILE *AFILE;
-	  sprintf(LocalRankFilename,"dA.%05i.raw",M.rank);
-	  AFILE = fopen(LocalRankFilename,"wb");
-	  fwrite(regdata.data(),8,N,AFILE);
-	  fclose(AFILE);
-
-	  regdata.fill(0.f);
-	  for (int k=0; k<Nz; k++){
-		  for (int j=0; j<Ny; j++){
-			  for (int i=0; i<Nx; i++){
-				  int idx=M.Map(i,j,k);
-				  if (!(idx<0)){
-					  value=DebugMassB[idx];
-					  regdata(i,j,k)=value;
-				  }
-			  }
-		  }
-	  }
-	  FILE *BFILE;
-	  sprintf(LocalRankFilename,"dB.%05i.raw",M.rank);
-	  BFILE = fopen(LocalRankFilename,"wb");
-	  fwrite(regdata.data(),8,N,BFILE);
-	  fclose(BFILE);
-	 */
-
 	// Need to initialize Aq, Bq, Den, Phi directly
 	//ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
 	ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
@@ -2301,101 +2227,368 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 }
 
 void FlowAdaptor::Flatten(ScaLBL_ColorModel &M){	
-	
-	  int Np = M.Np;
-	  double dA, dB;
-	  
-	  double *Aq_tmp, *Bq_tmp;
-	  
-	  Aq_tmp = new double [7*Np];
-	  Bq_tmp = new double [7*Np];
 
-	  ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double));
+	int Np = M.Np;
-	  ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
+	double dA, dB;
 
-	  for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){
+	double *Aq_tmp, *Bq_tmp;
-		  dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
-		  dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
-		  if (dA > 1.0){
-			  double mass_change = dA - 1.0;
-			  Aq_tmp[n] -= 0.333333333333333*mass_change;
-			  Aq_tmp[n+Np] -= 0.111111111111111*mass_change;
-			  Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
-			  Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
-			  Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
-			  Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
-			  Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
-		  }
-		  if (dB > 1.0){
-			  double mass_change = dB - 1.0;
-			  Bq_tmp[n] -= 0.333333333333333*mass_change;
-			  Bq_tmp[n+Np] -= 0.111111111111111*mass_change;
-			  Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
-			  Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
-			  Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
-			  Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
-			  Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
-		  }
-	  }
-	  for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){
-		  dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
-		  dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
-		  if (dA > 1.0){
-			  double mass_change = dA - 1.0;
-			  Aq_tmp[n] -= 0.333333333333333*mass_change;
-			  Aq_tmp[n+Np] -= 0.111111111111111*mass_change;
-			  Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
-			  Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
-			  Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
-			  Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
-			  Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
-		  }
-		  if (dB > 1.0){
-			  double mass_change = dB - 1.0;
-			  Bq_tmp[n] -= 0.333333333333333*mass_change;
-			  Bq_tmp[n+Np] -= 0.111111111111111*mass_change;
-			  Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
-			  Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
-			  Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
-			  Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
-			  Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
-		  }
-	   }
 
-	  ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
+	Aq_tmp = new double [7*Np];
-	  ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
+	Bq_tmp = new double [7*Np];
+
+	ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double));
+	ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
+
+	for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){
+		dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
+		dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
+		if (dA > 1.0){
+			double mass_change = dA - 1.0;
+			Aq_tmp[n] -= 0.333333333333333*mass_change;
+			Aq_tmp[n+Np] -= 0.111111111111111*mass_change;
+			Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
+			Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
+			Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
+			Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
+			Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
+		}
+		if (dB > 1.0){
+			double mass_change = dB - 1.0;
+			Bq_tmp[n] -= 0.333333333333333*mass_change;
+			Bq_tmp[n+Np] -= 0.111111111111111*mass_change;
+			Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
+			Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
+			Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
+			Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
+			Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
+		}
+	}
+	for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){
+		dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
+		dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
+		if (dA > 1.0){
+			double mass_change = dA - 1.0;
+			Aq_tmp[n] -= 0.333333333333333*mass_change;
+			Aq_tmp[n+Np] -= 0.111111111111111*mass_change;
+			Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
+			Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
+			Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
+			Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
+			Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
+		}
+		if (dB > 1.0){
+			double mass_change = dB - 1.0;
+			Bq_tmp[n] -= 0.333333333333333*mass_change;
+			Bq_tmp[n+Np] -= 0.111111111111111*mass_change;
+			Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
+			Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
+			Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
+			Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
+			Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
+		}
+	}
+
+	ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
+	ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
 }
 
 double FlowAdaptor::MoveInterface(ScaLBL_ColorModel &M){	
-	
-    double INTERFACE_CUTOFF = M.color_db->getWithDefault<double>( "move_interface_cutoff", 0.1 );
-    double MOVE_INTERFACE_FACTOR = M.color_db->getWithDefault<double>( "move_interface_factor", 10.0 );
 
-    ScaLBL_CopyToHost( phi.data(), M.Phi, Nx*Ny*Nz* sizeof( double ) );
+	double INTERFACE_CUTOFF = M.color_db->getWithDefault<double>( "move_interface_cutoff", 0.1 );
-    /* compute the local derivative of phase indicator field */
+	double MOVE_INTERFACE_FACTOR = M.color_db->getWithDefault<double>( "move_interface_factor", 10.0 );
-    double beta = M.beta;
+
-    double factor = 0.5/beta;
+	ScaLBL_CopyToHost( phi.data(), M.Phi, Nx*Ny*Nz* sizeof( double ) );
-    double total_interface_displacement = 0.0;
+	/* compute the local derivative of phase indicator field */
-    double total_interface_sites = 0.0;
+	double beta = M.beta;
-    for (int n=0; n<Nx*Ny*Nz; n++){
+	double factor = 0.5/beta;
-    	/* compute the distance to the interface */
+	double total_interface_displacement = 0.0;
-    	double value1 = M.Averages->Phi(n);
+	double total_interface_sites = 0.0;
-    	double dist1 = factor*log((1.0+value1)/(1.0-value1));
+	for (int n=0; n<Nx*Ny*Nz; n++){
-    	double value2 = phi(n);
+		/* compute the distance to the interface */
-    	double dist2 = factor*log((1.0+value2)/(1.0-value2));
+		double value1 = M.Averages->Phi(n);
-    	phi_t(n) = value2;
+		double dist1 = factor*log((1.0+value1)/(1.0-value1));
-    	if (value1 < INTERFACE_CUTOFF && value1 > -1*INTERFACE_CUTOFF && value2 < INTERFACE_CUTOFF && value2 > -1*INTERFACE_CUTOFF ){
+		double value2 = phi(n);
-    		/* time derivative of distance */
+		double dist2 = factor*log((1.0+value2)/(1.0-value2));
-    		double dxdt = 0.125*(dist2-dist1);
+		phi_t(n) = value2;
-    		/* extrapolate to move the distance further */
+		if (value1 < INTERFACE_CUTOFF && value1 > -1*INTERFACE_CUTOFF && value2 < INTERFACE_CUTOFF && value2 > -1*INTERFACE_CUTOFF ){
-    		double dist3 = dist2 + MOVE_INTERFACE_FACTOR*dxdt;
+			/* time derivative of distance */
-    		/* compute the new phase interface */
+			double dxdt = 0.125*(dist2-dist1);
-    		phi_t(n) = (2.f*(exp(-2.f*beta*(dist3)))/(1.f+exp(-2.f*beta*(dist3))) - 1.f);
+			/* extrapolate to move the distance further */
-    		total_interface_displacement += fabs(MOVE_INTERFACE_FACTOR*dxdt);
+			double dist3 = dist2 + MOVE_INTERFACE_FACTOR*dxdt;
-    		total_interface_sites += 1.0;
+			/* compute the new phase interface */
-    	}
+			phi_t(n) = (2.f*(exp(-2.f*beta*(dist3)))/(1.f+exp(-2.f*beta*(dist3))) - 1.f);
+			total_interface_displacement += fabs(MOVE_INTERFACE_FACTOR*dxdt);
+			total_interface_sites += 1.0;
+		}
 	}
-    ScaLBL_CopyToDevice( M.Phi, phi_t.data(), Nx*Ny*Nz* sizeof( double ) );	
+	ScaLBL_CopyToDevice( M.Phi, phi_t.data(), Nx*Ny*Nz* sizeof( double ) );	
-    return total_interface_sites;
+	return total_interface_sites;
 }
 
+double FlowAdaptor::ShellAggregation(ScaLBL_ColorModel &M, const double target_delta_volume){
+
+	const RankInfoStruct rank_info(M.rank,M.nprocx,M.nprocy,M.nprocz);
+	auto rank = M.rank;
+	auto Nx = M.Nx;  auto Ny = M.Ny; auto Nz = M.Nz;
+	auto N = Nx*Ny*Nz;
+	double vF = 0.f;
+	double vS = 0.f;
+	double delta_volume;
+	double WallFactor = 1.0;
+	bool USE_CONNECTED_NWP = false;
+
+	DoubleArray phase(Nx,Ny,Nz);
+	IntArray phase_label(Nx,Ny,Nz);;
+	DoubleArray phase_distance(Nx,Ny,Nz);
+	Array<char> phase_id(Nx,Ny,Nz);
+	fillHalo<double> fillDouble(M.Dm->Comm,M.Dm->rank_info,{Nx-2,Ny-2,Nz-2},{1,1,1},0,1);
+
+	// Basic algorithm to 
+	// 1. Copy phase field to CPU
+	ScaLBL_CopyToHost(phase.data(), M.Phi, N*sizeof(double));
+
+	double count = 0.f;
+	for (int k=1; k<Nz-1; k++){
+		for (int j=1; j<Ny-1; j++){
+			for (int i=1; i<Nx-1; i++){
+				if (phase(i,j,k) > 0.f && M.Averages->SDs(i,j,k) > 0.f) count+=1.f;
+			}
+		}
+	}
+	double volume_initial = M.Dm->Comm.sumReduce(  count);
+	double PoreVolume = M.Dm->Volume*M.Dm->Porosity();
+	/*ensure target isn't an absurdly small fraction of pore volume */
+	if (volume_initial < target_delta_volume*PoreVolume){
+		volume_initial = target_delta_volume*PoreVolume;
+	}
+
+	// 2. Identify connected components of phase field -> phase_label
+
+	double volume_connected = 0.0;
+	double second_biggest = 0.0;
+	if (USE_CONNECTED_NWP){
+		ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,M.Averages->SDs,vF,vS,phase_label,M.Dm->Comm);
+		M.Dm->Comm.barrier();
+
+		// only operate on component "0"ScaLBL_ColorModel &M, 
+		count = 0.0;
+
+		for (int k=0; k<Nz; k++){
+			for (int j=0; j<Ny; j++){
+				for (int i=0; i<Nx; i++){
+					int label = phase_label(i,j,k);
+					if (label == 0 ){
+						phase_id(i,j,k) = 0;
+						count += 1.0;
+					}
+					else 		
+						phase_id(i,j,k) = 1;
+					if (label == 1 ){
+						second_biggest += 1.0;
+					}
+				}
+			}
+		}	
+		volume_connected = M.Dm->Comm.sumReduce(  count);
+		second_biggest = M.Dm->Comm.sumReduce(  second_biggest);
+	}
+	else {
+		// use the whole NWP 
+		for (int k=0; k<Nz; k++){
+			for (int j=0; j<Ny; j++){
+				for (int i=0; i<Nx; i++){
+					if (M.Averages->SDs(i,j,k) > 0.f){
+						if (phase(i,j,k) > 0.f ){
+							phase_id(i,j,k) = 0;
+						}
+						else {
+							phase_id(i,j,k) = 1;
+						}
+					}
+					else {
+						phase_id(i,j,k) = 1;
+					}
+				}
+			}
+		}
+	}
+
+	// 3. Generate a distance map to the largest object -> phase_distance
+	CalcDist(phase_distance,phase_id,*M.Dm);
+
+	double temp,value;
+	double factor=0.5/M.beta;
+	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) < 3.f ){
+					value = phase(i,j,k);
+					if (value > 1.f)   value=1.f;
+					if (value < -1.f)  value=-1.f;
+					// temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm.
+					temp = -factor*log((1.0+value)/(1.0-value));
+					/// use this approximation close to the object
+					if (fabs(value) < 0.8 && M.Averages->SDs(i,j,k) > 1.f ){
+						phase_distance(i,j,k) = temp;
+					}
+					// erase the original object
+					phase(i,j,k) = -1.0;
+				}
+			}
+		}
+	}
+
+
+	if (rank==0) printf("Pathway volume / next largest ganglion %f \n",volume_connected/second_biggest );
+
+	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(M.Averages->SDs,phase_distance,phase_id,M.Averages->Dm, target_delta_volume_incremental, WallFactor);
+
+	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;
+			}
+		}
+	}	
+
+	CalcDist(phase_distance,phase_id,*M.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++){
+				double d = phase_distance(i,j,k);
+				if (M.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*M.beta*d))/(1.f+exp(-2.f*M.beta*d))-1.f);	
+					}
+				}
+			} 
+		}
+	}
+	fillDouble.fill(phase);
+
+	count = 0.f;
+	for (int k=1; k<Nz-1; k++){
+		for (int j=1; j<Ny-1; j++){
+			for (int i=1; i<Nx-1; i++){
+				if (phase(i,j,k) > 0.f && M.Averages->SDs(i,j,k) > 0.f){
+					count+=1.f;
+				}
+			}
+		}
+	}
+	double volume_final= M.Dm->Comm.sumReduce(  count);
+
+	delta_volume = (volume_final-volume_initial);
+	if (rank == 0)  printf("Shell Aggregation: change fluid volume fraction by %f \n", delta_volume/volume_initial);
+	if (rank == 0)  printf("   new saturation =  %f \n", volume_final/(M.Mask->Porosity()*double((Nx-2)*(Ny-2)*(Nz-2)*M.nprocs)));
+
+	// 6. copy back to the device
+	//if (rank==0)  printf("MorphInit: copy data  back to device\n");
+	ScaLBL_CopyToDevice(M.Phi,phase.data(),N*sizeof(double));
+
+	// 7. Re-initialize phase field and density
+	ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, 0, M.ScaLBL_Comm->LastExterior(), M.Np);
+	ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, M.ScaLBL_Comm->FirstInterior(), M.ScaLBL_Comm->LastInterior(), M.Np);
+	auto BoundaryCondition = M.BoundaryCondition;
+	if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
+		if (M.Dm->kproc()==0){
+			ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,0);
+			ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,1);
+			ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,2);
+		}
+		if (M.Dm->kproc() == M.nprocz-1){
+			ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-1);
+			ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-2);
+			ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-3);
+		}
+	}
+	return delta_volume;
+}
+
+
+double FlowAdaptor::SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil){
+  srand(time(NULL));
+  auto rank = M.rank;
+  auto Np = M.Np; 
+  double mass_loss =0.f;
+  double count =0.f;
+  double *Aq_tmp, *Bq_tmp;
+  
+  Aq_tmp = new double [7*Np];
+  Bq_tmp = new double [7*Np];
+
+  ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double));
+  ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
+  
+ 
+   for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){
+    double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
+    double dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
+    double dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
+    double phase_id = (dA - dB) / (dA + dB);
+    if (phase_id > 0.0){
+      Aq_tmp[n] -= 0.3333333333333333*random_value;
+      Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
+      
+      Bq_tmp[n] += 0.3333333333333333*random_value;
+      Bq_tmp[n+Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
+    }
+    mass_loss += random_value*seed_water_in_oil;
+  }
+
+  for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){
+    double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
+    double dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
+    double dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
+    double phase_id = (dA - dB) / (dA + dB);
+    if (phase_id > 0.0){
+      Aq_tmp[n] -= 0.3333333333333333*random_value;
+      Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
+      
+      Bq_tmp[n] += 0.3333333333333333*random_value;
+      Bq_tmp[n+Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
+    }
+    mass_loss += random_value*seed_water_in_oil;
+  }
+
+  count= M.Dm->Comm.sumReduce(  count);
+  mass_loss= M.Dm->Comm.sumReduce(  mass_loss);
+  if (rank == 0) printf("Remove mass %f from %f voxels \n",mass_loss,count);
+
+  // Need to initialize Aq, Bq, Den, Phi directly
+  //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
+  ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
+  ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
+
+  return(mass_loss);
+}

models/ColorModel.h

@@ -115,7 +115,9 @@ public:
 	~FlowAdaptor();
 	double MoveInterface(ScaLBL_ColorModel &M);
 	double ImageInit(ScaLBL_ColorModel &M, std::string Filename);
+	double ShellAggregation(ScaLBL_ColorModel &M, const double delta_volume);
 	double UpdateFractionalFlow(ScaLBL_ColorModel &M);
+	double SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil);
 	void Flatten(ScaLBL_ColorModel &M);
 	DoubleArray phi;
 	DoubleArray phi_t;

models/GreyscaleColorModel.cpp

@@ -309,18 +309,26 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
     //    oil-wet   < 0
     //    neutral   = 0 (i.e. no penalty) 
 	double *GreySolidW_host  = new double [Np];
+	double *GreySn_host      = new double [Np];
+	double *GreySw_host      = new double [Np];
 
 	size_t NLABELS=0;
 	signed char VALUE=0;
 	double AFFINITY=0.f;
+    double Sn,Sw;//end-point saturation of greynodes set by users
 
 	auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
 	auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
+	auto SnList = greyscaleColor_db->getVector<double>( "GreySnList" );
+	auto SwList = greyscaleColor_db->getVector<double>( "GreySwList" );
 
 	NLABELS=LabelList.size();
 	if (NLABELS != AffinityList.size()){
 		ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
 	}
+	if (NLABELS != SnList.size() || NLABELS != SwList.size()){
+		ERROR("Error: GreySolidLabels, GreySnList, and GreySwList must be the same length! \n");
+	}
 
 	for (int k=0;k<Nz;k++){
 		for (int j=0;j<Ny;j++){
@@ -328,16 +336,22 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
 				int n = k*Nx*Ny+j*Nx+i;
 				VALUE=id[n];
 	            AFFINITY=0.f;//all nodes except the specified grey nodes have grey-solid affinity = 0.0
+                Sn=99.0;
+                Sw=-99.0;
 				// Assign the affinity from the paired list
 				for (unsigned int idx=0; idx < NLABELS; idx++){
 					if (VALUE == LabelList[idx]){
 						AFFINITY=AffinityList[idx];
+                        Sn = SnList[idx];
+                        Sw = SwList[idx];
 						idx = NLABELS;
 					}
 				}
 				int idx = Map(i,j,k);
 				if (!(idx < 0)){
                     GreySolidW_host[idx] = AFFINITY;
+                    GreySn_host[idx]     = Sn;
+                    GreySw_host[idx]     = Sw;
                 }
 			}
 		}
@@ -348,15 +362,22 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
 		for (unsigned int idx=0; idx<NLABELS; idx++){
 			VALUE=LabelList[idx];
 			AFFINITY=AffinityList[idx];
-			printf("   grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY); 
+			Sn=SnList[idx];
+			Sw=SwList[idx];
+			//printf("   grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY); 
+			printf("   grey-solid label=%d, grey-solid affinity=%.3g, grey-solid Sn=%.3g, grey-solid Sw=%.3g\n",VALUE,AFFINITY,Sn,Sw); 
 		}
 		printf("NOTE: grey-solid affinity>0: water-wet || grey-solid affinity<0: oil-wet \n");
 	}
     
     
 	ScaLBL_CopyToDevice(GreySolidW, GreySolidW_host, Np*sizeof(double));
+	ScaLBL_CopyToDevice(GreySn, GreySn_host, Np*sizeof(double));
+	ScaLBL_CopyToDevice(GreySw, GreySw_host, Np*sizeof(double));
 	ScaLBL_Comm->Barrier();
     delete [] GreySolidW_host;
+    delete [] GreySn_host;
+    delete [] GreySw_host;
 }
 ////----------------------------------------------------------------------------------------------------------//
 
@@ -598,6 +619,8 @@ void ScaLBL_GreyscaleColorModel::Create(){
     //ScaLBL_AllocateDeviceMemory((void **) &GreySolidPhi, sizeof(double)*Nx*Ny*Nz);		
     //ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np);		
     ScaLBL_AllocateDeviceMemory((void **) &GreySolidW, sizeof(double)*Np);		
+    ScaLBL_AllocateDeviceMemory((void **) &GreySn, sizeof(double)*Np);		
+    ScaLBL_AllocateDeviceMemory((void **) &GreySw, sizeof(double)*Np);		
     ScaLBL_AllocateDeviceMemory((void **) &Porosity_dvc, sizeof(double)*Np);
     ScaLBL_AllocateDeviceMemory((void **) &Permeability_dvc, sizeof(double)*Np);
 	//...........................................................................
@@ -921,7 +944,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
 		// Halo exchange for phase field
 		ScaLBL_Comm_Regular->SendHalo(Phi);
         //Model-1&4 with capillary pressure penalty for grey nodes
-        ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
+        ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
                 rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, 
                 alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
         //Model-1&4
@@ -950,7 +973,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
 		}
 
         //Model-1&4 with capillary pressure penalty for grey nodes
-        ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
+        ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
                 rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
                 alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
         //Model-1&4
@@ -983,7 +1006,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
 		}
 		ScaLBL_Comm_Regular->SendHalo(Phi);
         //Model-1&4 with capillary pressure penalty for grey nodes
-        ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure, 
+        ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, 
                 rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
                 alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
         //Model-1&4
@@ -1012,7 +1035,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
 		}
 
         //Model-1&4 with capillary pressure penalty for grey nodes
-        ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
+        ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
                 rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
                 alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
         //Model-1&4

models/GreyscaleColorModel.h

@@ -68,6 +68,8 @@ public:
     //double *GreySolidPhi; //Model 2 & 3
     //double *GreySolidGrad;//Model 1 & 4
     double *GreySolidW;
+    double *GreySn;
+    double *GreySw;
 	//double *ColorGrad;
 	double *Velocity;
 	double *Pressure;

models/StokesModel.cpp

@@ -59,6 +59,7 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
     BoundaryCondition = 0;
 	if (stokes_db->keyExists( "BC" )){
 		BoundaryCondition = stokes_db->getScalar<int>( "BC" );
+
 	}
 	if (stokes_db->keyExists( "tolerance" )){
 		tolerance = stokes_db->getScalar<double>( "tolerance" );
@@ -110,6 +111,7 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
 void ScaLBL_StokesModel::ReadParams(string filename){
     //NOTE the max time step is left unspecified
 
+
 	// read the input database 
 	db = std::make_shared<Database>( filename );
 	domain_db = db->getDatabase( "Domain" );
@@ -298,8 +300,10 @@ void ScaLBL_StokesModel::AssignZetaPotentialSolid(double *zeta_potential_solid)
 		ERROR("Error: LB Stokes Solver: SolidLabels and ZetaPotentialSolidList must be the same length! \n");
 	}
 
-	double label_count[NLABELS];
+	double *label_count;
-	double label_count_global[NLABELS];
+	double *label_count_global;
+	label_count = new double [NLABELS];
+	label_count_global = new double [NLABELS];
 
 	for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
 

tests/lbpm_color_simulator.cpp

@@ -66,7 +66,6 @@ int main( int argc, char **argv )
 		// structure and allocate variables
 		ColorModel.Initialize(); // initializing the model will set initial conditions for variables
 
-
 		if (SimulationMode == "development"){
 			double MLUPS=0.0;
 			int timestep = 0;
@@ -82,13 +81,18 @@ int main( int argc, char **argv )
 			int SKIP_TIMESTEPS = 0;
 			int MAX_STEADY_TIME = 1000000;
 			double ENDPOINT_THRESHOLD = 0.1;
-			double FRACTIONAL_FLOW_INCREMENT = 0.05;
+			double FRACTIONAL_FLOW_INCREMENT = 0.0; // this will skip the flow adaptor if not enabled
+			double SEED_WATER = 0.0;
 			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 );
+				SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 50000 );
-				FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
 				ENDPOINT_THRESHOLD = flow_db->getWithDefault<double>( "endpoint_threshold", 0.1);
+				/* protocol specific key values */
+				if (PROTOCOL == "fractional flow")
+					FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
+				if (PROTOCOL == "seed water")
+					SEED_WATER = flow_db->getWithDefault<double>( "seed_water", 0.01);
 			}
 			/* analysis keys*/
 			int ANALYSIS_INTERVAL = ColorModel.timestepMax;
@@ -106,9 +110,26 @@ int main( int argc, char **argv )
 				if (ColorModel.timestep > ColorModel.timestepMax){
 					ContinueSimulation = false;
 				}
+				
+				/* Load a new image if image sequence is specified */
+				if (PROTOCOL == "image sequence"){
+					IMAGE_INDEX++;
+					if (IMAGE_INDEX < IMAGE_COUNT){
+						std::string next_image = ImageList[IMAGE_INDEX];
+						if (rank==0) printf("***Loading next image in sequence (%i) ***\n",IMAGE_INDEX);
+						ColorModel.color_db->putScalar<int>("image_index",IMAGE_INDEX);
+						Adapt.ImageInit(ColorModel, next_image);
+					}
+					else{
+						if (rank==0) printf("Finished simulating image sequence \n");
+						ColorModel.timestep =  ColorModel.timestepMax;
+					}
+				}
+				/*********************************************************/
 				/* update the fluid configuration with the flow adapter */
 				int skip_time = 0;
 				timestep = ColorModel.timestep;
+				/* get the averaged flow measures computed internally for the last simulation point*/
 				double SaturationChange = 0.0;
 				double volB = ColorModel.Averages->gwb.V; 
 				double volA = ColorModel.Averages->gnb.V; 
@@ -121,6 +142,7 @@ int main( int argc, char **argv )
 				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);
+				/* stop simulation if previous point was sufficiently close to the endpoint*/
 				if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false;
 				if (ContinueSimulation){
 					while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){
@@ -128,21 +150,16 @@ int main( int argc, char **argv )
 						if (PROTOCOL == "fractional flow")	{							
 							Adapt.UpdateFractionalFlow(ColorModel);
 						}
-						else if (PROTOCOL == "image sequence"){
+						else if (PROTOCOL == "shell aggregation"){
-							// Use image sequence
+							double target_volume_change = FRACTIONAL_FLOW_INCREMENT*initialSaturation - SaturationChange;
-							IMAGE_INDEX++;
+							Adapt.ShellAggregation(ColorModel,target_volume_change);
-							if (IMAGE_INDEX < IMAGE_COUNT){
-								std::string next_image = ImageList[IMAGE_INDEX];
-								if (rank==0) printf("***Loading next image in sequence (%i) ***\n",IMAGE_INDEX);
-								ColorModel.color_db->putScalar<int>("image_index",IMAGE_INDEX);
-								Adapt.ImageInit(ColorModel, next_image);
-							}
-							else{
-								if (rank==0) printf("Finished simulating image sequence \n");
-								ColorModel.timestep =  ColorModel.timestepMax;
-							}
 						}
+						else if (PROTOCOL == "seed water"){
+							Adapt.SeedPhaseField(ColorModel,SEED_WATER);
+						}
+						/* Run some LBM timesteps to let the system relax a bit */
 						MLUPS = ColorModel.Run(timestep);
+						/* Recompute the volume fraction now that the system has adjusted */
 						double volB = ColorModel.Averages->gwb.V; 
 						double volA = ColorModel.Averages->gnb.V;
 						SaturationChange = volB/(volA + volB) - initialSaturation;
@@ -150,22 +167,12 @@ int main( int argc, char **argv )
 					}
 					if (rank==0) printf("  *********************************************************************  \n");
 					if (rank==0) printf("   Updated fractional flow with saturation change = %f  \n", SaturationChange);
+					if (rank==0) printf("   Used protocol = %s  \n", PROTOCOL.c_str());
 					if (rank==0) printf("  *********************************************************************  \n");
 				}
-				/* apply timestep skipping algorithm to accelerate steady-state */
+				/*********************************************************/
-				/* skip_time = 0;
-            	timestep = ColorModel.timestep;
-            	while (skip_time < SKIP_TIMESTEPS){
-            		timestep += ANALYSIS_INTERVAL;
-                	MLUPS = ColorModel.Run(timestep);
-            		Adapt.MoveInterface(ColorModel);
-            		skip_time += ANALYSIS_INTERVAL;
-            	}
-				 */
 			}
-			//ColorModel.WriteDebug();
 		}
-
 		else
 			ColorModel.Run();