diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index d3833bcb..df87fcc4 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -72,6 +72,27 @@ 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 Fx, double Fy, double Fz, 
                                              double *Poros,double *Perm, double *Velocity,double Den,double *Pressure);
 
+// GREYSCALE COLOR MODEL
+
+extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(double *dist, double *Aq, double *Bq, double *Den,
+                double *DenGradA, double *DenGradB, double *SolidForce, int start, int finish, int Np,
+                double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double Gsc, double Gx, double Gy, double Gz,
+                double *Poros,double *Perm, double *Velocity,double *Pressure);
+
+extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, double *dist, double *Aq, double *Bq, double *Den,
+                double *DenGradA, double *DenGradB, double *SolidForce, int start, int finish, int Np,
+                double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double Gsc, double Gx, double Gy, double Gz,
+                double *Poros,double *Perm, double *Velocity,double *Pressure);
+
+extern "C" void ScaLBL_D3Q7_GreyColorIMRT_Init(double *Den, double *Aq, double *Bq, int start, int finish, int Np);
+
+extern "C" void ScaLBL_D3Q19_GreyColorIMRT_Init(double *dist, double *Den, double rhoA, double rhoB, int Np);
+
+extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleColorDensity(int *NeighborList, double *Aq, double *Bq, double *Den, int start, int finish, int Np);
+
+extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleColorDensity(double *Aq, double *Bq, double *Den, int start, int finish, int Np);
+
+extern "C" void ScaLBL_D3Q19_GreyscaleColor_Gradient(int *neighborList, double *Den, double *DenGrad, int start, int finish, int Np);
 
 // MRT MODEL
 extern "C" void ScaLBL_D3Q19_AAeven_MRT(double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx,
diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp
index 2c0e686d..564eb96d 100644
--- a/cpu/D3Q19.cpp
+++ b/cpu/D3Q19.cpp
@@ -84,33 +84,6 @@ extern "C" void ScaLBL_D3Q19_Init(double *dist, int Np)
 	}
 }
 
-
-extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den)
-{
-	int n;
-	for (n=0; n<Np; n++){
-		dist[n] = Den - 0.6666666666666667;
-		dist[Np+n] = 0.055555555555555555;		//double(100*n)+1.f;
-		dist[2*Np+n] = 0.055555555555555555;	//double(100*n)+2.f;
-		dist[3*Np+n] = 0.055555555555555555;	//double(100*n)+3.f;
-		dist[4*Np+n] = 0.055555555555555555;	//double(100*n)+4.f;
-		dist[5*Np+n] = 0.055555555555555555;	//double(100*n)+5.f;
-		dist[6*Np+n] = 0.055555555555555555;	//double(100*n)+6.f;
-		dist[7*Np+n] = 0.0277777777777778;   //double(100*n)+7.f;
-		dist[8*Np+n] = 0.0277777777777778;   //double(100*n)+8.f;
-		dist[9*Np+n] = 0.0277777777777778;   //double(100*n)+9.f;
-		dist[10*Np+n] = 0.0277777777777778;  //double(100*n)+10.f;
-		dist[11*Np+n] = 0.0277777777777778;  //double(100*n)+11.f;
-		dist[12*Np+n] = 0.0277777777777778;  //double(100*n)+12.f;
-		dist[13*Np+n] = 0.0277777777777778;  //double(100*n)+13.f;
-		dist[14*Np+n] = 0.0277777777777778;  //double(100*n)+14.f;
-		dist[15*Np+n] = 0.0277777777777778;  //double(100*n)+15.f;
-		dist[16*Np+n] = 0.0277777777777778;  //double(100*n)+16.f;
-		dist[17*Np+n] = 0.0277777777777778;  //double(100*n)+17.f;
-		dist[18*Np+n] = 0.0277777777777778;  //double(100*n)+18.f;
-	}
-}
-
 //*************************************************************************
 extern "C" void ScaLBL_D3Q19_Swap(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz)
 {
diff --git a/cpu/Greyscale.cpp b/cpu/Greyscale.cpp
index b4b017c8..8ed95a6a 100644
--- a/cpu/Greyscale.cpp
+++ b/cpu/Greyscale.cpp
@@ -1565,3 +1565,30 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dis
 	}
 }
 
+
+extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den)
+{
+	int n;
+	for (n=0; n<Np; n++){
+		dist[n] = Den - 0.6666666666666667;
+		dist[Np+n] = 0.055555555555555555;		//double(100*n)+1.f;
+		dist[2*Np+n] = 0.055555555555555555;	//double(100*n)+2.f;
+		dist[3*Np+n] = 0.055555555555555555;	//double(100*n)+3.f;
+		dist[4*Np+n] = 0.055555555555555555;	//double(100*n)+4.f;
+		dist[5*Np+n] = 0.055555555555555555;	//double(100*n)+5.f;
+		dist[6*Np+n] = 0.055555555555555555;	//double(100*n)+6.f;
+		dist[7*Np+n] = 0.0277777777777778;   //double(100*n)+7.f;
+		dist[8*Np+n] = 0.0277777777777778;   //double(100*n)+8.f;
+		dist[9*Np+n] = 0.0277777777777778;   //double(100*n)+9.f;
+		dist[10*Np+n] = 0.0277777777777778;  //double(100*n)+10.f;
+		dist[11*Np+n] = 0.0277777777777778;  //double(100*n)+11.f;
+		dist[12*Np+n] = 0.0277777777777778;  //double(100*n)+12.f;
+		dist[13*Np+n] = 0.0277777777777778;  //double(100*n)+13.f;
+		dist[14*Np+n] = 0.0277777777777778;  //double(100*n)+14.f;
+		dist[15*Np+n] = 0.0277777777777778;  //double(100*n)+15.f;
+		dist[16*Np+n] = 0.0277777777777778;  //double(100*n)+16.f;
+		dist[17*Np+n] = 0.0277777777777778;  //double(100*n)+17.f;
+		dist[18*Np+n] = 0.0277777777777778;  //double(100*n)+18.f;
+	}
+}
+
diff --git a/gpu/D3Q19.cu b/gpu/D3Q19.cu
index ccd125b2..3ae06f3e 100644
--- a/gpu/D3Q19.cu
+++ b/gpu/D3Q19.cu
@@ -267,37 +267,6 @@ __global__ void dvc_ScaLBL_D3Q19_Init(double *dist, int Np)
 	}
 }
 
-__global__ void dvc_ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den)
-{
-	int n;
-	int S = Np/NBLOCKS/NTHREADS + 1;
-	for (int s=0; s<S; s++){
-		//........Get 1-D index for this thread....................
-		n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
-		if (n<Np ){
-			dist[n] = Den - 0.6666666666666667;
-			dist[Np+n] = 0.055555555555555555;		//double(100*n)+1.f;
-			dist[2*Np+n] = 0.055555555555555555;	//double(100*n)+2.f;
-			dist[3*Np+n] = 0.055555555555555555;	//double(100*n)+3.f;
-			dist[4*Np+n] = 0.055555555555555555;	//double(100*n)+4.f;
-			dist[5*Np+n] = 0.055555555555555555;	//double(100*n)+5.f;
-			dist[6*Np+n] = 0.055555555555555555;	//double(100*n)+6.f;
-			dist[7*Np+n] = 0.0277777777777778;   //double(100*n)+7.f;
-			dist[8*Np+n] = 0.0277777777777778;   //double(100*n)+8.f;
-			dist[9*Np+n] = 0.0277777777777778;   //double(100*n)+9.f;
-			dist[10*Np+n] = 0.0277777777777778;  //double(100*n)+10.f;
-			dist[11*Np+n] = 0.0277777777777778;  //double(100*n)+11.f;
-			dist[12*Np+n] = 0.0277777777777778;  //double(100*n)+12.f;
-			dist[13*Np+n] = 0.0277777777777778;  //double(100*n)+13.f;
-			dist[14*Np+n] = 0.0277777777777778;  //double(100*n)+14.f;
-			dist[15*Np+n] = 0.0277777777777778;  //double(100*n)+15.f;
-			dist[16*Np+n] = 0.0277777777777778;  //double(100*n)+16.f;
-			dist[17*Np+n] = 0.0277777777777778;  //double(100*n)+17.f;
-			dist[18*Np+n] = 0.0277777777777778;  //double(100*n)+18.f;
-		}
-	}
-}
-
 //*************************************************************************
 __global__  void dvc_ScaLBL_D3Q19_Swap_Compact(int *neighborList, double *disteven, double *distodd, int Np, int q){
 	int n,nn;
@@ -2358,13 +2327,6 @@ extern "C" void ScaLBL_D3Q19_Init(double *dist, int Np){
 	}
 }
 
-extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den){
-	dvc_ScaLBL_D3Q19_GreyIMRT_Init<<<NBLOCKS,NTHREADS >>>(dist, Np, Den);
-	cudaError_t err = cudaGetLastError();
-	if (cudaSuccess != err){
-		printf("CUDA error in ScaLBL_D3Q19_GreyIMRT_Init: %s \n",cudaGetErrorString(err));
-	}
-}
 
 extern "C" void ScaLBL_D3Q19_Swap(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz){
 	dvc_ScaLBL_D3Q19_Swap<<<NBLOCKS,NTHREADS >>>(ID, disteven, distodd, Nx, Ny, Nz);
diff --git a/gpu/Greyscale.cu b/gpu/Greyscale.cu
index 0a9a63e0..a3fba989 100644
--- a/gpu/Greyscale.cu
+++ b/gpu/Greyscale.cu
@@ -1590,6 +1590,37 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double
 	}
 }
 
+__global__ void dvc_ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den)
+{
+	int n;
+	int S = Np/NBLOCKS/NTHREADS + 1;
+	for (int s=0; s<S; s++){
+		//........Get 1-D index for this thread....................
+		n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
+		if (n<Np ){
+			dist[n] = Den - 0.6666666666666667;
+			dist[Np+n] = 0.055555555555555555;		//double(100*n)+1.f;
+			dist[2*Np+n] = 0.055555555555555555;	//double(100*n)+2.f;
+			dist[3*Np+n] = 0.055555555555555555;	//double(100*n)+3.f;
+			dist[4*Np+n] = 0.055555555555555555;	//double(100*n)+4.f;
+			dist[5*Np+n] = 0.055555555555555555;	//double(100*n)+5.f;
+			dist[6*Np+n] = 0.055555555555555555;	//double(100*n)+6.f;
+			dist[7*Np+n] = 0.0277777777777778;   //double(100*n)+7.f;
+			dist[8*Np+n] = 0.0277777777777778;   //double(100*n)+8.f;
+			dist[9*Np+n] = 0.0277777777777778;   //double(100*n)+9.f;
+			dist[10*Np+n] = 0.0277777777777778;  //double(100*n)+10.f;
+			dist[11*Np+n] = 0.0277777777777778;  //double(100*n)+11.f;
+			dist[12*Np+n] = 0.0277777777777778;  //double(100*n)+12.f;
+			dist[13*Np+n] = 0.0277777777777778;  //double(100*n)+13.f;
+			dist[14*Np+n] = 0.0277777777777778;  //double(100*n)+14.f;
+			dist[15*Np+n] = 0.0277777777777778;  //double(100*n)+15.f;
+			dist[16*Np+n] = 0.0277777777777778;  //double(100*n)+16.f;
+			dist[17*Np+n] = 0.0277777777777778;  //double(100*n)+17.f;
+			dist[18*Np+n] = 0.0277777777777778;  //double(100*n)+18.f;
+		}
+	}
+}
+
 
 extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Fx, double Fy, double Fz,double *Poros,double *Perm, double *Velocity,double *Pressure){
 	
@@ -1631,3 +1662,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dis
 	}
 }
 
+extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den){
+	dvc_ScaLBL_D3Q19_GreyIMRT_Init<<<NBLOCKS,NTHREADS >>>(dist, Np, Den);
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q19_GreyIMRT_Init: %s \n",cudaGetErrorString(err));
+	}
+}
diff --git a/gpu/GreyscaleColor.cu b/gpu/GreyscaleColor.cu
new file mode 100644
index 00000000..c9586d3c
--- /dev/null
+++ b/gpu/GreyscaleColor.cu
@@ -0,0 +1,1590 @@
+#include <stdio.h>
+
+#define NBLOCKS 1024
+#define NTHREADS 256
+
+__global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(double *dist, double *Aq, double *Bq, double *Den,
+                double *DenGradA, double *DenGradB, double *SolidForce, int start, int finish, int Np,
+                double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double Gsc, double Gx, double Gy, double Gz,
+                double *Poros,double *Perm, double *Velocity,double *Pressure){
+
+	int n;
+	double vx,vy,vz,v_mag;
+    double ux,uy,uz,u_mag;
+    double pressure;//defined for this incompressible model
+	// conserved momemnts
+	double jx,jy,jz;
+	// non-conserved moments
+	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
+    double fq;
+    // currently disable 'GeoFun'
+    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 Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
+	double tau,tau_eff,rlx_setA,rlx_setB;
+    double mu_eff;//effective kinematic viscosity for Darcy term
+    double rho0;
+    double phi;
+    double nx,ny,nz,C;
+    double nA,nB;
+	double a1,b1,a2,b2,nAB,delta;
+    double beta=0.95;
+    double nA_gradx,nA_grady,nA_gradz;
+    double nB_gradx,nB_grady,nB_gradz;
+    double Gff_x,Gff_y,Gff_z;
+    double Gfs_x,Gfs_y,Gfs_z;
+
+
+	const double mrt_V1=0.05263157894736842;
+	const double mrt_V2=0.012531328320802;
+	const double mrt_V3=0.04761904761904762;
+	const double mrt_V4=0.004594820384294068;
+	const double mrt_V5=0.01587301587301587;
+	const double mrt_V6=0.0555555555555555555555555;
+	const double mrt_V7=0.02777777777777778;
+	const double mrt_V8=0.08333333333333333;
+	const double mrt_V9=0.003341687552213868;
+	const double mrt_V10=0.003968253968253968;
+	const double mrt_V11=0.01388888888888889;
+	const double mrt_V12=0.04166666666666666;
+
+
+	int S = Np/NBLOCKS/NTHREADS + 1;
+	for (int s=0; s<S; s++){
+	    //........Get 1-D index for this thread....................
+	    n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
+
+		if ( n<finish ){
+
+			// read the component number densities
+			nA = Den[n];
+			nB = Den[Np + n];
+			nA_gradx = DenGradA[n+0*Np];
+			nA_grady = DenGradA[n+1*Np];
+			nA_gradz = DenGradA[n+2*Np];
+			nB_gradx = DenGradB[n+0*Np];
+			nB_grady = DenGradB[n+1*Np];
+			nB_gradz = DenGradB[n+2*Np];
+
+			// compute phase indicator field
+			phi=(nA-nB)/(nA+nB);
+
+			// local density
+			rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
+			// local relaxation time
+			tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
+			rlx_setA = 1.f/tau;
+			rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
+			tau_eff=tauA_eff + 0.5*(1.0-phi)*(tauB_eff-tauA_eff);
+            mu_eff = (tau_eff-0.5)/3.f;//kinematic viscosity
+
+
+            //........................................................................
+            //					READ THE DISTRIBUTIONS
+            //		(read from opposite array due to previous swap operation)
+            //........................................................................
+            // q=0
+            fq = dist[n];
+            m1  = -30.0*fq;
+            m2  = 12.0*fq;
+
+            // q=1
+            fq = dist[2*Np+n];
+            pressure = fq;
+            m1 -= 11.0*fq;
+            m2 -= 4.0*fq;
+            jx = fq;
+            m4 = -4.0*fq;
+            m9 = 2.0*fq;
+            m10 = -4.0*fq;
+
+            // f2 = dist[10*Np+n];
+            fq = dist[1*Np+n];
+            pressure += fq;
+            m1 -= 11.0*(fq);
+            m2 -= 4.0*(fq);
+            jx -= fq;
+            m4 += 4.0*(fq);
+            m9 += 2.0*(fq);
+            m10 -= 4.0*(fq);
+
+            // q=3
+            fq = dist[4*Np+n];
+            pressure += fq;
+            m1 -= 11.0*fq;
+            m2 -= 4.0*fq;
+            jy = fq;
+            m6 = -4.0*fq;
+            m9 -= fq;
+            m10 += 2.0*fq;
+            m11 = fq;
+            m12 = -2.0*fq;
+
+            // q = 4
+            fq = dist[3*Np+n];
+            pressure += fq;
+            m1 -= 11.0*fq;
+            m2 -= 4.0*fq;
+            jy -= fq;
+            m6 += 4.0*fq;
+            m9 -= fq;
+            m10 += 2.0*fq;
+            m11 += fq;
+            m12 -= 2.0*fq;
+
+            // q=5
+            fq = dist[6*Np+n];
+            pressure += fq;
+            m1 -= 11.0*fq;
+            m2 -= 4.0*fq;
+            jz = fq;
+            m8 = -4.0*fq;
+            m9 -= fq;
+            m10 += 2.0*fq;
+            m11 -= fq;
+            m12 += 2.0*fq;
+
+            // q = 6
+            fq = dist[5*Np+n];
+            pressure += fq;
+            m1 -= 11.0*fq;
+            m2 -= 4.0*fq;
+            jz -= fq;
+            m8 += 4.0*fq;
+            m9 -= fq;
+            m10 += 2.0*fq;
+            m11 -= fq;
+            m12 += 2.0*fq;
+
+            // q=7
+            fq = dist[8*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx += fq;
+            m4 += fq;
+            jy += fq;
+            m6 += fq;
+            m9  += fq;
+            m10 += fq;
+            m11 += fq;
+            m12 += fq;
+            m13 = fq;
+            m16 = fq;
+            m17 = -fq;
+
+            // q = 8
+            fq = dist[7*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx -= fq;
+            m4 -= fq;
+            jy -= fq;
+            m6 -= fq;
+            m9 += fq;
+            m10 += fq;
+            m11 += fq;
+            m12 += fq;
+            m13 += fq;
+            m16 -= fq;
+            m17 += fq;
+
+            // q=9
+            fq = dist[10*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx += fq;
+            m4 += fq;
+            jy -= fq;
+            m6 -= fq;
+            m9 += fq;
+            m10 += fq;
+            m11 += fq;
+            m12 += fq;
+            m13 -= fq;
+            m16 += fq;
+            m17 += fq;
+
+            // q = 10
+            fq = dist[9*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx -= fq;
+            m4 -= fq;
+            jy += fq;
+            m6 += fq;
+            m9 += fq;
+            m10 += fq;
+            m11 += fq;
+            m12 += fq;
+            m13 -= fq;
+            m16 -= fq;
+            m17 -= fq;
+
+            // q=11
+            fq = dist[12*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx += fq;
+            m4 += fq;
+            jz += fq;
+            m8 += fq;
+            m9 += fq;
+            m10 += fq;
+            m11 -= fq;
+            m12 -= fq;
+            m15 = fq;
+            m16 -= fq;
+            m18 = fq;
+
+            // q=12
+            fq = dist[11*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx -= fq;
+            m4 -= fq;
+            jz -= fq;
+            m8 -= fq;
+            m9 += fq;
+            m10 += fq;
+            m11 -= fq;
+            m12 -= fq;
+            m15 += fq;
+            m16 += fq;
+            m18 -= fq;
+
+            // q=13
+            fq = dist[14*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx += fq;
+            m4 += fq;
+            jz -= fq;
+            m8 -= fq;
+            m9 += fq;
+            m10 += fq;
+            m11 -= fq;
+            m12 -= fq;
+            m15 -= fq;
+            m16 -= fq;
+            m18 -= fq;
+
+            // q=14
+            fq = dist[13*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx -= fq;
+            m4 -= fq;
+            jz += fq;
+            m8 += fq;
+            m9 += fq;
+            m10 += fq;
+            m11 -= fq;
+            m12 -= fq;
+            m15 -= fq;
+            m16 += fq;
+            m18 += fq;
+
+            // q=15
+            fq = dist[16*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jy += fq;
+            m6 += fq;
+            jz += fq;
+            m8 += fq;
+            m9 -= 2.0*fq;
+            m10 -= 2.0*fq;
+            m14 = fq;
+            m17 += fq;
+            m18 -= fq;
+
+            // q=16
+            fq = dist[15*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jy -= fq;
+            m6 -= fq;
+            jz -= fq;
+            m8 -= fq;
+            m9 -= 2.0*fq;
+            m10 -= 2.0*fq;
+            m14 += fq;
+            m17 -= fq;
+            m18 += fq;
+
+            // q=17
+            fq = dist[18*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jy += fq;
+            m6 += fq;
+            jz -= fq;
+            m8 -= fq;
+            m9 -= 2.0*fq;
+            m10 -= 2.0*fq;
+            m14 -= fq;
+            m17 += fq;
+            m18 += fq;
+
+            // q=18
+            fq = dist[17*Np+n];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jy -= fq;
+            m6 -= fq;
+            jz += fq;
+            m8 += fq;
+            m9 -= 2.0*fq;
+            m10 -= 2.0*fq;
+            m14 -= fq;
+            m17 -= fq;
+            m18 -= fq;
+            //---------------------------------------------------------------------//
+
+            //---------------- Calculate SC fluid-fluid and fluid-solid forces ---------------//
+            // fluid-fluid force
+            Gff_x = -Gsc*nA*nB_gradx*int(phi>0.0)-Gsc*nB*nA_gradx*int(phi<0.0);
+            Gff_y = -Gsc*nA*nB_grady*int(phi>0.0)-Gsc*nB*nA_grady*int(phi<0.0);
+            Gff_z = -Gsc*nA*nB_gradz*int(phi>0.0)-Gsc*nB*nA_gradz*int(phi<0.0);
+            // fluid-solid force
+            Gfs_x = (nA-nB)*SolidForce[n+0*Np];    
+            Gfs_y = (nA-nB)*SolidForce[n+1*Np];    
+            Gfs_z = (nA-nB)*SolidForce[n+2*Np];    
+
+            porosity = Poros[n];
+            // use local saturation as an estimation of effective relperm values
+            perm = Perm[n]*nA/(nA+nB)*int(phi>0.0)+Perm[n]*nB/(nA+nB)*int(phi<0.0);
+
+            c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
+            if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
+            //GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
+            c1 = porosity*0.5*GeoFun/sqrt(perm);
+            if (porosity==1.0) c1 = 0.0;//i.e. apparent pore nodes
+
+            vx = jx/rho0+0.5*(porosity*Gx+Gff_x+Gfs_x);
+            vy = jy/rho0+0.5*(porosity*Gy+Gff_y+Gfs_y);
+            vz = jz/rho0+0.5*(porosity*Gz+Gff_z+Gfs_z);
+            v_mag=sqrt(vx*vx+vy*vy+vz*vz);
+            ux = vx/(c0+sqrt(c0*c0+c1*v_mag));
+            uy = vy/(c0+sqrt(c0*c0+c1*v_mag));
+            uz = vz/(c0+sqrt(c0*c0+c1*v_mag));
+            u_mag=sqrt(ux*ux+uy*uy+uz*uz);
+
+            //Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
+            Fx = rho0*(-porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx + Gff_x + Gfs_x);
+            Fy = rho0*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy + Gff_y + Gfs_y);
+            Fz = rho0*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz + Gff_z + Gfs_z);
+            if (porosity==1.0){
+                Fx=rho0*(Gx + Gff_x + Gfs_x);
+                Fy=rho0*(Gy + Gff_y + Gfs_y);
+                Fz=rho0*(Gz + Gff_z + Gfs_z);
+            }
+
+            //Calculate pressure for Incompressible-MRT model
+            pressure=0.5/porosity*(pressure-0.5*rho0*u_mag*u_mag/porosity);
+
+//            //..............carry out relaxation process...............................................
+//            m1 = m1 + rlx_setA*((-30*rho0+19*(ux*ux+uy*uy+uz*uz)/porosity + 57*pressure*porosity) - m1) 
+//                    + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity;
+//            m2 = m2 + rlx_setA*((12*rho0 - 5.5*(ux*ux+uy*uy+uz*uz)/porosity-27*pressure*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*((rho0*(2*ux*ux-uy*uy-uz*uz)/porosity) - m9)
+//                    + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity;
+//            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*((rho0*(uy*uy-uz*uz)/porosity) - m11)
+//                      + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity;
+//            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*((rho0*ux*uy/porosity) - m13)
+//                      + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity;
+//            m14 = m14 + rlx_setA*((rho0*uy*uz/porosity) - m14)
+//                      + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity;
+//            m15 = m15 + rlx_setA*((rho0*ux*uz/porosity) - 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);
+//            //.......................................................................................................
+
+            //-------------------- IMRT collison where body force has NO higher-order terms -------------//
+            //..............carry out relaxation process...............................................
+            m1 = m1 + rlx_setA*((-30*rho0+19*(ux*ux+uy*uy+uz*uz)/porosity + 57*pressure*porosity) - m1);
+            m2 = m2 + rlx_setA*((12*rho0 - 5.5*(ux*ux+uy*uy+uz*uz)/porosity-27*pressure*porosity) - m2);
+            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*((rho0*(2*ux*ux-uy*uy-uz*uz)/porosity) - m9);
+            m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10);
+            m11 = m11 + rlx_setA*((rho0*(uy*uy-uz*uz)/porosity) - m11);
+            m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12);
+            m13 = m13 + rlx_setA*((rho0*ux*uy/porosity) - m13);
+            m14 = m14 + rlx_setA*((rho0*uy*uz/porosity) - m14);
+            m15 = m15 + rlx_setA*((rho0*ux*uz/porosity) - m15);
+            m16 = m16 + rlx_setB*( - m16);
+            m17 = m17 + rlx_setB*( - m17);
+            m18 = m18 + rlx_setB*( - m18);
+            //.......................................................................................................
+
+            //.................inverse transformation......................................................
+            // q=0
+            fq = mrt_V1*rho0-mrt_V2*m1+mrt_V3*m2;
+            dist[n] = fq;
+
+            // q = 1
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10);
+            dist[1*Np+n] = fq;
+
+            // q=2
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10);
+            dist[2*Np+n] = fq;
+
+            // q = 3
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12);
+            dist[3*Np+n] = fq;
+
+            // q = 4
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12);
+            dist[4*Np+n] = fq;
+
+            // q = 5
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11);
+            dist[5*Np+n] = fq;
+
+            // q = 6
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11);
+            dist[6*Np+n] = fq;
+
+            // q = 7
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17);
+            dist[7*Np+n] = fq;
+
+            // q = 8
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m17-m16);
+            dist[8*Np+n] = fq;
+
+            // q = 9
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17);
+            dist[9*Np+n] = fq;
+
+            // q = 10
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17);
+            dist[10*Np+n] = fq;
+
+            // q = 11
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m18-m16);
+            dist[11*Np+n] = fq;
+
+            // q = 12
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18);
+            dist[12*Np+n] = fq;
+
+            // q = 13
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15-0.125*(m16+m18);
+            dist[13*Np+n] = fq;
+
+            // q= 14
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15+0.125*(m16+m18);
+            dist[14*Np+n] = fq;
+
+            // q = 15
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18);
+            dist[15*Np+n] = fq;
+
+            // q = 16
+            fq =  mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17);
+            dist[16*Np+n] = fq;
+
+            // q = 17
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8)-mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18);
+            dist[17*Np+n] = fq;
+
+            // q = 18
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6)-mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18);
+            dist[18*Np+n] = fq;
+            //........................................................................
+
+            //Update velocity on device
+            Velocity[0*Np+n] = ux;
+            Velocity[1*Np+n] = uy;
+            Velocity[2*Np+n] = uz;
+            //Update pressure on device
+            Pressure[n] = pressure;
+
+            //-----------------------Mass transport------------------------//
+            // Calculate the color gradient
+            nx = (2*nB*nA_gradx-2*nA*nB_gradx)/(nA+nB)/(nA+nB); 
+            ny = (2*nB*nA_grady-2*nA*nB_grady)/(nA+nB)/(nA+nB); 
+            nz = (2*nB*nA_gradz-2*nA*nB_gradz)/(nA+nB)/(nA+nB); 
+			//...........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;		
+			if (C == 0.0)	nx = ny = nz = 0.0;
+
+			// Instantiate mass transport distributions
+			// Stationary value - distribution 0
+			nAB = 1.0/(nA+nB);
+			Aq[n] = 0.3333333333333333*nA;
+			Bq[n] = 0.3333333333333333*nB;
+
+			//...............................................
+			// q = 0,2,4
+			// 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;
+			a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
+			b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
+			a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta;
+			b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta;
+
+			Aq[1*Np+n] = a1;
+			Bq[1*Np+n] = b1;
+			Aq[2*Np+n] = a2;
+			Bq[2*Np+n] = b2;
+
+			//...............................................
+			// q = 2
+			// Cq = {0,1,0}
+			delta = beta*nA*nB*nAB*0.1111111111111111*ny;
+			if (!(nA*nB*nAB>0)) delta=0;
+			a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
+			b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
+			a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta;
+			b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta;
+
+			Aq[3*Np+n] = a1;
+			Bq[3*Np+n] = b1;
+			Aq[4*Np+n] = a2;
+			Bq[4*Np+n] = b2;
+			//...............................................
+			// q = 4
+			// Cq = {0,0,1}
+			delta = beta*nA*nB*nAB*0.1111111111111111*nz;
+			if (!(nA*nB*nAB>0)) delta=0;
+			a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
+			b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
+			a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta;
+			b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta;
+
+			Aq[5*Np+n] = a1;
+			Bq[5*Np+n] = b1;
+			Aq[6*Np+n] = a2;
+			Bq[6*Np+n] = b2;
+			//...............................................
+
+		}
+	}
+}
+
+__global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, double *dist, double *Aq, double *Bq, double *Den,
+                double *DenGradA, double *DenGradB, double *SolidForce, int start, int finish, int Np,
+                double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double Gsc, double Gx, double Gy, double Gz,
+                double *Poros,double *Perm, double *Velocity,double *Pressure){
+
+	int n, nread, nr1,nr2,nr3,nr4,nr5,nr6;
+	double vx,vy,vz,v_mag;
+    double ux,uy,uz,u_mag;
+    double pressure;//defined for this incompressible model
+	// conserved momemnts
+	double jx,jy,jz;
+	// non-conserved moments
+	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
+    double fq;
+    // currently disable 'GeoFun'
+    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 Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
+	double tau,tau_eff,rlx_setA,rlx_setB;
+    double mu_eff;//effective kinematic viscosity for Darcy term
+    double rho0;
+    double phi;
+    double nx,ny,nz,C;
+    double nA,nB;
+	double a1,b1,a2,b2,nAB,delta;
+    double beta=0.95;
+    double nA_gradx,nA_grady,nA_gradz;
+    double nB_gradx,nB_grady,nB_gradz;
+    double Gff_x,Gff_y,Gff_z;
+    double Gfs_x,Gfs_y,Gfs_z;
+
+	const double mrt_V1=0.05263157894736842;
+	const double mrt_V2=0.012531328320802;
+	const double mrt_V3=0.04761904761904762;
+	const double mrt_V4=0.004594820384294068;
+	const double mrt_V5=0.01587301587301587;
+	const double mrt_V6=0.0555555555555555555555555;
+	const double mrt_V7=0.02777777777777778;
+	const double mrt_V8=0.08333333333333333;
+	const double mrt_V9=0.003341687552213868;
+	const double mrt_V10=0.003968253968253968;
+	const double mrt_V11=0.01388888888888889;
+	const double mrt_V12=0.04166666666666666;
+
+	int S = Np/NBLOCKS/NTHREADS + 1;
+	for (int s=0; s<S; s++){
+	    //........Get 1-D index for this thread....................
+	    n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
+
+		if ( n<finish ){		
+
+			// read the component number densities
+			nA = Den[n];
+			nB = Den[Np + n];
+			nA_gradx = DenGradA[n+0*Np];
+			nA_grady = DenGradA[n+1*Np];
+			nA_gradz = DenGradA[n+2*Np];
+			nB_gradx = DenGradB[n+0*Np];
+			nB_grady = DenGradB[n+1*Np];
+			nB_gradz = DenGradB[n+2*Np];
+
+			// compute phase indicator field
+			phi=(nA-nB)/(nA+nB);
+
+			// local density
+			rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
+			// local relaxation time
+			tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
+			rlx_setA = 1.f/tau;
+			rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
+			tau_eff=tauA_eff + 0.5*(1.0-phi)*(tauB_eff-tauA_eff);
+            mu_eff = (tau_eff-0.5)/3.f;//kinematic viscosity
+
+            //........................................................................
+            //					READ THE DISTRIBUTIONS
+            //		(read from opposite array due to previous swap operation)
+            //........................................................................
+            // q=0
+            fq = dist[n];
+            m1  = -30.0*fq;
+            m2  = 12.0*fq;
+
+            // q=1
+            nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
+            fq = dist[nr1]; // reading the f1 data into register fq
+            pressure = fq;
+            m1 -= 11.0*fq;
+            m2 -= 4.0*fq;
+            jx = fq;
+            m4 = -4.0*fq;
+            m9 = 2.0*fq;
+            m10 = -4.0*fq;
+
+            // q=2
+            nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
+            fq = dist[nr2];  // reading the f2 data into register fq
+            pressure += fq;
+            m1 -= 11.0*(fq);
+            m2 -= 4.0*(fq);
+            jx -= fq;
+            m4 += 4.0*(fq);
+            m9 += 2.0*(fq);
+            m10 -= 4.0*(fq);
+
+            // q=3
+            nr3 = neighborList[n+2*Np]; // neighbor 4
+            fq = dist[nr3];
+            pressure += fq;
+            m1 -= 11.0*fq;
+            m2 -= 4.0*fq;
+            jy = fq;
+            m6 = -4.0*fq;
+            m9 -= fq;
+            m10 += 2.0*fq;
+            m11 = fq;
+            m12 = -2.0*fq;
+
+            // q = 4
+            nr4 = neighborList[n+3*Np]; // neighbor 3
+            fq = dist[nr4];
+            pressure += fq;
+            m1 -= 11.0*fq;
+            m2 -= 4.0*fq;
+            jy -= fq;
+            m6 += 4.0*fq;
+            m9 -= fq;
+            m10 += 2.0*fq;
+            m11 += fq;
+            m12 -= 2.0*fq;
+
+            // q=5
+            nr5 = neighborList[n+4*Np];
+            fq = dist[nr5];
+            pressure += fq;
+            m1 -= 11.0*fq;
+            m2 -= 4.0*fq;
+            jz = fq;
+            m8 = -4.0*fq;
+            m9 -= fq;
+            m10 += 2.0*fq;
+            m11 -= fq;
+            m12 += 2.0*fq;
+
+            // q = 6
+            nr6 = neighborList[n+5*Np];
+            fq = dist[nr6];
+            pressure += fq;
+            m1 -= 11.0*fq;
+            m2 -= 4.0*fq;
+            jz -= fq;
+            m8 += 4.0*fq;
+            m9 -= fq;
+            m10 += 2.0*fq;
+            m11 -= fq;
+            m12 += 2.0*fq;
+
+            // q=7
+            nread = neighborList[n+6*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx += fq;
+            m4 += fq;
+            jy += fq;
+            m6 += fq;
+            m9  += fq;
+            m10 += fq;
+            m11 += fq;
+            m12 += fq;
+            m13 = fq;
+            m16 = fq;
+            m17 = -fq;
+
+            // q = 8
+            nread = neighborList[n+7*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx -= fq;
+            m4 -= fq;
+            jy -= fq;
+            m6 -= fq;
+            m9 += fq;
+            m10 += fq;
+            m11 += fq;
+            m12 += fq;
+            m13 += fq;
+            m16 -= fq;
+            m17 += fq;
+
+            // q=9
+            nread = neighborList[n+8*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx += fq;
+            m4 += fq;
+            jy -= fq;
+            m6 -= fq;
+            m9 += fq;
+            m10 += fq;
+            m11 += fq;
+            m12 += fq;
+            m13 -= fq;
+            m16 += fq;
+            m17 += fq;
+
+            // q = 10
+            nread = neighborList[n+9*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx -= fq;
+            m4 -= fq;
+            jy += fq;
+            m6 += fq;
+            m9 += fq;
+            m10 += fq;
+            m11 += fq;
+            m12 += fq;
+            m13 -= fq;
+            m16 -= fq;
+            m17 -= fq;
+
+            // q=11
+            nread = neighborList[n+10*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx += fq;
+            m4 += fq;
+            jz += fq;
+            m8 += fq;
+            m9 += fq;
+            m10 += fq;
+            m11 -= fq;
+            m12 -= fq;
+            m15 = fq;
+            m16 -= fq;
+            m18 = fq;
+
+            // q=12
+            nread = neighborList[n+11*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx -= fq;
+            m4 -= fq;
+            jz -= fq;
+            m8 -= fq;
+            m9 += fq;
+            m10 += fq;
+            m11 -= fq;
+            m12 -= fq;
+            m15 += fq;
+            m16 += fq;
+            m18 -= fq;
+
+            // q=13
+            nread = neighborList[n+12*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx += fq;
+            m4 += fq;
+            jz -= fq;
+            m8 -= fq;
+            m9 += fq;
+            m10 += fq;
+            m11 -= fq;
+            m12 -= fq;
+            m15 -= fq;
+            m16 -= fq;
+            m18 -= fq;
+
+            // q=14
+            nread = neighborList[n+13*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jx -= fq;
+            m4 -= fq;
+            jz += fq;
+            m8 += fq;
+            m9 += fq;
+            m10 += fq;
+            m11 -= fq;
+            m12 -= fq;
+            m15 -= fq;
+            m16 += fq;
+            m18 += fq;
+
+            // q=15
+            nread = neighborList[n+14*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jy += fq;
+            m6 += fq;
+            jz += fq;
+            m8 += fq;
+            m9 -= 2.0*fq;
+            m10 -= 2.0*fq;
+            m14 = fq;
+            m17 += fq;
+            m18 -= fq;
+
+            // q=16
+            nread = neighborList[n+15*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jy -= fq;
+            m6 -= fq;
+            jz -= fq;
+            m8 -= fq;
+            m9 -= 2.0*fq;
+            m10 -= 2.0*fq;
+            m14 += fq;
+            m17 -= fq;
+            m18 += fq;
+
+            // q=17
+            nread = neighborList[n+16*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jy += fq;
+            m6 += fq;
+            jz -= fq;
+            m8 -= fq;
+            m9 -= 2.0*fq;
+            m10 -= 2.0*fq;
+            m14 -= fq;
+            m17 += fq;
+            m18 += fq;
+
+            // q=18
+            nread = neighborList[n+17*Np];
+            fq = dist[nread];
+            pressure += fq;
+            m1 += 8.0*fq;
+            m2 += fq;
+            jy -= fq;
+            m6 -= fq;
+            jz += fq;
+            m8 += fq;
+            m9 -= 2.0*fq;
+            m10 -= 2.0*fq;
+            m14 -= fq;
+            m17 -= fq;
+            m18 -= fq;
+            //---------------------------------------------------------------------//
+
+            //---------------- Calculate SC fluid-fluid and fluid-solid forces ---------------//
+            // fluid-fluid force
+            Gff_x = -Gsc*nA*nB_gradx*int(phi>0.0)-Gsc*nB*nA_gradx*int(phi<0.0);
+            Gff_y = -Gsc*nA*nB_grady*int(phi>0.0)-Gsc*nB*nA_grady*int(phi<0.0);
+            Gff_z = -Gsc*nA*nB_gradz*int(phi>0.0)-Gsc*nB*nA_gradz*int(phi<0.0);
+            // fluid-solid force
+            Gfs_x = (nA-nB)*SolidForce[n+0*Np];    
+            Gfs_y = (nA-nB)*SolidForce[n+1*Np];    
+            Gfs_z = (nA-nB)*SolidForce[n+2*Np];    
+
+            porosity = Poros[n];
+            // use local saturation as an estimation of effective relperm values
+            perm = Perm[n]*nA/(nA+nB)*int(phi>0.0)+Perm[n]*nB/(nA+nB)*int(phi<0.0);
+
+            c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
+            if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
+            //GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
+            c1 = porosity*0.5*GeoFun/sqrt(perm);
+            if (porosity==1.0) c1 = 0.0;//i.e. apparent pore nodes
+
+            vx = jx/rho0+0.5*(porosity*Gx+Gff_x+Gfs_x);
+            vy = jy/rho0+0.5*(porosity*Gy+Gff_y+Gfs_y);
+            vz = jz/rho0+0.5*(porosity*Gz+Gff_z+Gfs_z);
+            v_mag=sqrt(vx*vx+vy*vy+vz*vz);
+            ux = vx/(c0+sqrt(c0*c0+c1*v_mag));
+            uy = vy/(c0+sqrt(c0*c0+c1*v_mag));
+            uz = vz/(c0+sqrt(c0*c0+c1*v_mag));
+            u_mag=sqrt(ux*ux+uy*uy+uz*uz);
+
+            //Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
+            Fx = rho0*(-porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx + Gff_x + Gfs_x);
+            Fy = rho0*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy + Gff_y + Gfs_y);
+            Fz = rho0*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz + Gff_z + Gfs_z);
+            if (porosity==1.0){
+                Fx=rho0*(Gx + Gff_x + Gfs_x);
+                Fy=rho0*(Gy + Gff_y + Gfs_y);
+                Fz=rho0*(Gz + Gff_z + Gfs_z);
+            }
+
+            //Calculate pressure for Incompressible-MRT model
+            pressure=0.5/porosity*(pressure-0.5*rho0*u_mag*u_mag/porosity);
+
+//            //..............carry out relaxation process...............................................
+//            m1 = m1 + rlx_setA*((-30*rho0+19*(ux*ux+uy*uy+uz*uz)/porosity + 57*pressure*porosity) - m1) 
+//                    + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity;
+//            m2 = m2 + rlx_setA*((12*rho0 - 5.5*(ux*ux+uy*uy+uz*uz)/porosity-27*pressure*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*((rho0*(2*ux*ux-uy*uy-uz*uz)/porosity) - m9)
+//                    + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity;
+//            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*((rho0*(uy*uy-uz*uz)/porosity) - m11)
+//                      + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity;
+//            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*((rho0*ux*uy/porosity) - m13)
+//                      + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity;
+//            m14 = m14 + rlx_setA*((rho0*uy*uz/porosity) - m14)
+//                      + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity;
+//            m15 = m15 + rlx_setA*((rho0*ux*uz/porosity) - 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);
+//            //.......................................................................................................
+           
+            //-------------------- IMRT collison where body force has NO higher-order terms -------------//
+            //..............carry out relaxation process...............................................
+            m1 = m1 + rlx_setA*((-30*rho0+19*(ux*ux+uy*uy+uz*uz)/porosity + 57*pressure*porosity) - m1);
+            m2 = m2 + rlx_setA*((12*rho0 - 5.5*(ux*ux+uy*uy+uz*uz)/porosity-27*pressure*porosity) - m2);
+            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*((rho0*(2*ux*ux-uy*uy-uz*uz)/porosity) - m9);
+            m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10);
+            m11 = m11 + rlx_setA*((rho0*(uy*uy-uz*uz)/porosity) - m11);
+            m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12);
+            m13 = m13 + rlx_setA*((rho0*ux*uy/porosity) - m13);
+            m14 = m14 + rlx_setA*((rho0*uy*uz/porosity) - m14);
+            m15 = m15 + rlx_setA*((rho0*ux*uz/porosity) - m15);
+            m16 = m16 + rlx_setB*( - m16);
+            m17 = m17 + rlx_setB*( - m17);
+            m18 = m18 + rlx_setB*( - m18);
+            //.......................................................................................................
+
+
+            //.................inverse transformation......................................................
+            // q=0
+            fq = mrt_V1*rho0-mrt_V2*m1+mrt_V3*m2;
+            dist[n] = fq;
+
+            // q = 1
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10);
+            //nread = neighborList[n+Np];
+            dist[nr2] = fq;
+
+            // q=2
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10);
+            //nread = neighborList[n];
+            dist[nr1] = fq;
+
+            // q = 3
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12);
+            //nread = neighborList[n+3*Np];
+            dist[nr4] = fq;
+
+            // q = 4
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12);
+            //nread = neighborList[n+2*Np];
+            dist[nr3] = fq;
+
+            // q = 5
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11);
+            //nread = neighborList[n+5*Np];
+            dist[nr6] = fq;
+
+            // q = 6
+            fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11);
+            //nread = neighborList[n+4*Np];
+            dist[nr5] = fq;
+
+            // q = 7
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17);
+            nread = neighborList[n+7*Np];
+            dist[nread] = fq;
+
+            // q = 8
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m17-m16);
+            nread = neighborList[n+6*Np];
+            dist[nread] = fq;
+
+            // q = 9
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17);
+            nread = neighborList[n+9*Np];
+            dist[nread] = fq;
+
+            // q = 10
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17);
+            nread = neighborList[n+8*Np];
+            dist[nread] = fq;
+
+            // q = 11
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m18-m16);
+            nread = neighborList[n+11*Np];
+            dist[nread] = fq;
+
+            // q = 12
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18);
+            nread = neighborList[n+10*Np];
+            dist[nread]= fq;
+
+            // q = 13
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15-0.125*(m16+m18);
+            nread = neighborList[n+13*Np];
+            dist[nread] = fq;
+
+            // q= 14
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15+0.125*(m16+m18);
+            nread = neighborList[n+12*Np];
+            dist[nread] = fq;
+
+            // q = 15
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18);
+            nread = neighborList[n+15*Np];
+            dist[nread] = fq;
+
+            // q = 16
+            fq =  mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17);
+            nread = neighborList[n+14*Np];
+            dist[nread] = fq;
+
+            // q = 17
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8)-mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18);
+            nread = neighborList[n+17*Np];
+            dist[nread] = fq;
+
+            // q = 18
+            fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6)-mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18);
+            nread = neighborList[n+16*Np];
+            dist[nread] = fq;
+            //........................................................................
+
+            //Update velocity on device
+            Velocity[0*Np+n] = ux;
+            Velocity[1*Np+n] = uy;
+            Velocity[2*Np+n] = uz;
+            //Update pressure on device
+            Pressure[n] = pressure;
+
+            //-----------------------Mass transport------------------------//
+            // Calculate the color gradient
+            nx = (2*nB*nA_gradx-2*nA*nB_gradx)/(nA+nB)/(nA+nB); 
+            ny = (2*nB*nA_grady-2*nA*nB_grady)/(nA+nB)/(nA+nB); 
+            nz = (2*nB*nA_gradz-2*nA*nB_gradz)/(nA+nB)/(nA+nB); 
+			//...........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;		
+			if (C == 0.0)	nx = ny = nz = 0.0;
+
+			// Instantiate mass transport distributions
+			// Stationary value - distribution 0
+			nAB = 1.0/(nA+nB);
+			Aq[n] = 0.3333333333333333*nA;
+			Bq[n] = 0.3333333333333333*nB;
+
+			//...............................................
+			// q = 0,2,4
+			// 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;
+			a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
+			b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
+			a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta;
+			b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta;
+
+			// q = 1
+			//nread = neighborList[n+Np];
+			Aq[nr2] = a1;
+			Bq[nr2] = b1;
+			// q=2
+			//nread = neighborList[n];
+			Aq[nr1] = a2;
+			Bq[nr1] = b2;
+
+			//...............................................
+			// Cq = {0,1,0}
+			delta = beta*nA*nB*nAB*0.1111111111111111*ny;
+			if (!(nA*nB*nAB>0)) delta=0;
+			a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
+			b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
+			a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta;
+			b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta;
+
+			// q = 3
+			//nread = neighborList[n+3*Np];
+			Aq[nr4] = a1;
+			Bq[nr4] = b1;
+			// q = 4
+			//nread = neighborList[n+2*Np];
+			Aq[nr3] = a2;
+			Bq[nr3] = b2;
+
+			//...............................................
+			// q = 4
+			// Cq = {0,0,1}
+			delta = beta*nA*nB*nAB*0.1111111111111111*nz;
+			if (!(nA*nB*nAB>0)) delta=0;
+			a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
+			b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
+			a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta;
+			b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta;
+
+			// q = 5
+			//nread = neighborList[n+5*Np];
+			Aq[nr6] = a1;
+			Bq[nr6] = b1;
+			// q = 6
+			//nread = neighborList[n+4*Np];
+			Aq[nr5] = a2;
+			Bq[nr5] = b2;
+			//...............................................
+		}
+	}
+}
+
+__global__ void dvc_ScaLBL_D3Q19_GreyColorIMRT_Init(double *dist, double *Den, double rhoA, double rhoB, int Np){
+	int n;
+	int S = Np/NBLOCKS/NTHREADS + 1;
+    double phi;
+    double nA,nB;
+    double Den0;
+	for (int s=0; s<S; s++){
+		//........Get 1-D index for this thread....................
+		n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
+		if (n<Np ){
+            nA = Den[n];
+            nB = Den[n+Np];
+            phi = (nA-nB)/(nA+nB);
+            Den0 = 0.5*(1.f+phi)*rhoA + 0.5*(1.f-phi)*rhoB;
+
+			dist[n] = Den0 - 0.6666666666666667;
+			dist[Np+n] = 0.055555555555555555;		//double(100*n)+1.f;
+			dist[2*Np+n] = 0.055555555555555555;	//double(100*n)+2.f;
+			dist[3*Np+n] = 0.055555555555555555;	//double(100*n)+3.f;
+			dist[4*Np+n] = 0.055555555555555555;	//double(100*n)+4.f;
+			dist[5*Np+n] = 0.055555555555555555;	//double(100*n)+5.f;
+			dist[6*Np+n] = 0.055555555555555555;	//double(100*n)+6.f;
+			dist[7*Np+n] = 0.0277777777777778;   //double(100*n)+7.f;
+			dist[8*Np+n] = 0.0277777777777778;   //double(100*n)+8.f;
+			dist[9*Np+n] = 0.0277777777777778;   //double(100*n)+9.f;
+			dist[10*Np+n] = 0.0277777777777778;  //double(100*n)+10.f;
+			dist[11*Np+n] = 0.0277777777777778;  //double(100*n)+11.f;
+			dist[12*Np+n] = 0.0277777777777778;  //double(100*n)+12.f;
+			dist[13*Np+n] = 0.0277777777777778;  //double(100*n)+13.f;
+			dist[14*Np+n] = 0.0277777777777778;  //double(100*n)+14.f;
+			dist[15*Np+n] = 0.0277777777777778;  //double(100*n)+15.f;
+			dist[16*Np+n] = 0.0277777777777778;  //double(100*n)+16.f;
+			dist[17*Np+n] = 0.0277777777777778;  //double(100*n)+17.f;
+			dist[18*Np+n] = 0.0277777777777778;  //double(100*n)+18.f;
+		}
+	}
+}
+
+__global__ void dvc_ScaLBL_D3Q7_GreyColorIMRT_Init(double *Den, double *Aq, double *Bq, int start, int finish, int Np){
+	int idx;
+    double nA,nB;
+
+	int S = Np/NBLOCKS/NTHREADS + 1;
+	for (int s=0; s<S; s++){
+		//........Get 1-D index for this thread....................
+		idx =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
+		if (idx<finish) {
+            nA = Den[idx];
+            nB = Den[idx+Np];
+
+			Aq[idx]=0.3333333333333333*nA;
+			Aq[Np+idx]=0.1111111111111111*nA;
+			Aq[2*Np+idx]=0.1111111111111111*nA;
+			Aq[3*Np+idx]=0.1111111111111111*nA;
+			Aq[4*Np+idx]=0.1111111111111111*nA;
+			Aq[5*Np+idx]=0.1111111111111111*nA;
+			Aq[6*Np+idx]=0.1111111111111111*nA;
+
+			Bq[idx]=0.3333333333333333*nB;
+			Bq[Np+idx]=0.1111111111111111*nB;
+			Bq[2*Np+idx]=0.1111111111111111*nB;
+			Bq[3*Np+idx]=0.1111111111111111*nB;
+			Bq[4*Np+idx]=0.1111111111111111*nB;
+			Bq[5*Np+idx]=0.1111111111111111*nB;
+			Bq[6*Np+idx]=0.1111111111111111*nB;
+		}
+	}
+}
+
+__global__  void dvc_ScaLBL_D3Q7_AAodd_GreyscaleColorDensity(int *neighborList, double *Aq, double *Bq, double *Den, int start, int finish, int Np){
+	int n,nread;
+	double fq,nA,nB;
+
+	int S = Np/NBLOCKS/NTHREADS + 1;
+	for (int s=0; s<S; s++){
+		//........Get 1-D index for this thread....................
+		n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
+		if (n<finish) {
+			//..........Compute the number density for each component ............
+			// q=0
+			fq = Aq[n];
+			nA = fq;
+			fq = Bq[n];
+			nB = fq;
+			
+			// q=1
+			nread = neighborList[n]; 
+			fq = Aq[nread];
+			nA += fq;
+			fq = Bq[nread]; 
+			nB += fq;
+			
+			// q=2
+			nread = neighborList[n+Np]; 
+			fq = Aq[nread];  
+			nA += fq;
+			fq = Bq[nread]; 
+			nB += fq;
+			
+			// q=3
+			nread = neighborList[n+2*Np]; 
+			fq = Aq[nread];
+			nA += fq;
+			fq = Bq[nread];
+			nB += fq;
+			
+			// q = 4
+			nread = neighborList[n+3*Np]; 
+			fq = Aq[nread];
+			nA += fq;
+			fq = Bq[nread];
+			nB += fq;
+
+			// q=5
+			nread = neighborList[n+4*Np];
+			fq = Aq[nread];
+			nA += fq;
+			fq = Bq[nread];
+			nB += fq;
+			
+			// q = 6
+			nread = neighborList[n+5*Np];
+			fq = Aq[nread];
+			nA += fq;
+			fq = Bq[nread];
+			nB += fq;
+
+			// save the number densities
+			Den[n] = nA;
+			Den[Np+n] = nB;
+		}
+	}
+}
+
+__global__  void dvc_ScaLBL_D3Q7_AAeven_GreyscaleColorDensity(double *Aq, double *Bq, double *Den, int start, int finish, int Np){
+	int n;
+	double fq,nA,nB;
+
+	int S = Np/NBLOCKS/NTHREADS + 1;
+	for (int s=0; s<S; s++){
+		//........Get 1-D index for this thread....................
+		n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
+		if (n<finish) {
+			// compute number density for each component
+			// q=0
+			fq = Aq[n];
+			nA = fq;
+			fq = Bq[n];
+			nB = fq;
+			
+			// q=1
+			fq = Aq[2*Np+n];
+			nA += fq;
+			fq = Bq[2*Np+n];
+			nB += fq;
+
+			// q=2
+			fq = Aq[1*Np+n];
+			nA += fq;
+			fq = Bq[1*Np+n];
+			nB += fq;
+
+			// q=3
+			fq = Aq[4*Np+n];
+			nA += fq;
+			fq = Bq[4*Np+n];
+			nB += fq;
+
+			// q = 4
+			fq = Aq[3*Np+n];
+			nA += fq;
+			fq = Bq[3*Np+n];
+			nB += fq;
+			
+			// q=5
+			fq = Aq[6*Np+n];
+			nA += fq;
+			fq = Bq[6*Np+n];
+			nB += fq;
+			
+			// q = 6
+			fq = Aq[5*Np+n];
+			nA += fq;
+			fq = Bq[5*Np+n];
+			nB += fq;
+
+			// save the number densities
+			Den[n] = nA;
+			Den[Np+n] = nB;
+		}
+	}
+}
+
+__global__ void dvc_ScaLBL_D3Q19_GreyscaleColor_Gradient(int *neighborList, double *Den, double *DenGrad, int start, int finish, int Np){
+
+	int n,nn;
+	// distributions
+	double m1,m2,m3,m4,m5,m6,m7,m8,m9;
+	double m10,m11,m12,m13,m14,m15,m16,m17,m18;
+	double nx,ny,nz;
+
+	int S = Np/NBLOCKS/NTHREADS + 1;
+	for (int s=0; s<S; s++){
+		//........Get 1-D index for this thread....................
+		n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
+		if (n<finish) {
+			nn = neighborList[n+Np]%Np;
+			m1 = Den[nn]*int(n!=nn);
+			nn = neighborList[n]%Np;
+			m2 = Den[nn]*int(n!=nn);
+			nn = neighborList[n+3*Np]%Np;
+			m3 = Den[nn]*int(n!=nn);
+			nn = neighborList[n+2*Np]%Np;
+			m4 = Den[nn]*int(n!=nn);		
+			nn = neighborList[n+5*Np]%Np;
+			m5 = Den[nn]*int(n!=nn);
+			nn = neighborList[n+4*Np]%Np;
+			m6 = Den[nn]*int(n!=nn);		
+			nn = neighborList[n+7*Np]%Np;
+			m7 = Den[nn]*int(n!=nn);
+			nn = neighborList[n+6*Np]%Np;
+			m8 = Den[nn]*int(n!=nn);		
+			nn = neighborList[n+9*Np]%Np;
+			m9 = Den[nn]*int(n!=nn);
+			nn = neighborList[n+8*Np]%Np;
+			m10 = Den[nn]*int(n!=nn);		
+			nn = neighborList[n+11*Np]%Np;
+			m11 = Den[nn]*int(n!=nn);
+			nn = neighborList[n+10*Np]%Np;
+			m12 = Den[nn]*int(n!=nn);		
+			nn = neighborList[n+13*Np]%Np;
+			m13 = Den[nn]*int(n!=nn);
+			nn = neighborList[n+12*Np]%Np;
+			m14 = Den[nn]*int(n!=nn);		
+			nn = neighborList[n+15*Np]%Np;
+			m15 = Den[nn]*int(n!=nn);
+			nn = neighborList[n+14*Np]%Np;
+			m16 = Den[nn]*int(n!=nn);		
+			nn = neighborList[n+17*Np]%Np;
+			m17 = Den[nn]*int(n!=nn);
+			nn = neighborList[n+16*Np]%Np;
+			m18 = Den[nn]*int(n!=nn);					
+			
+			//............Compute the Color Gradient...................................
+			nx = 1.f/18.f*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+			ny = 1.f/18.f*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+			nz = 1.f/18.f*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+			
+			DenGrad[n] = nx;
+			DenGrad[Np+n] = ny;
+			DenGrad[2*Np+n] = nz;
+		}
+	}
+}
+
+
+extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(double *dist, double *Aq, double *Bq, double *Den,
+                double *DenGradA, double *DenGradB, double *SolidForce, int start, int finish, int Np,
+                double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double Gsc, double Gx, double Gy, double Gz,
+                double *Poros,double *Perm, double *Velocity,double *Pressure){
+
+    dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor<<<NBLOCKS,NTHREADS >>>(dist, Aq, Bq, Den, DenGradA, DenGradB, SolidForce, start, finish, Np,
+                                                                 tauA, tauB, tauA_eff, tauB_eff, rhoA, rhoB, Gsc, Gx, Gy, Gz, Poros, Perm, Velocity, Pressure);
+
+    cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q19_AAeven_GreyscaleColor: %s \n",cudaGetErrorString(err));
+	}
+}
+
+extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, double *dist, double *Aq, double *Bq, double *Den,
+                double *DenGradA, double *DenGradB, double *SolidForce, int start, int finish, int Np,
+                double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double Gsc, double Gx, double Gy, double Gz,
+                double *Poros,double *Perm, double *Velocity,double *Pressure){
+
+    dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor<<<NBLOCKS,NTHREADS >>>(neighborList, dist, Aq, Bq, Den, DenGradA, DenGradB, SolidForce, start, finish, Np,
+                                                                 tauA, tauB, tauA_eff, tauB_eff, rhoA, rhoB, Gsc, Gx, Gy, Gz, Poros, Perm, Velocity, Pressure);
+
+    cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q19_AAodd_GreyscaleColor: %s \n",cudaGetErrorString(err));
+	}
+}
+
+extern "C" void ScaLBL_D3Q7_GreyColorIMRT_Init(double *Den, double *Aq, double *Bq, int start, int finish, int Np){
+	dvc_ScaLBL_D3Q7_GreyColorIMRT_Init<<<NBLOCKS,NTHREADS >>>(Den, Aq, Bq, start, finish, Np);
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q7_GreyColorIMRT_Init: %s \n",cudaGetErrorString(err));
+	}
+}
+
+extern "C" void ScaLBL_D3Q19_GreyColorIMRT_Init(double *dist, double *Den, double rhoA, double rhoB, int Np){
+	dvc_ScaLBL_D3Q19_GreyColorIMRT_Init<<<NBLOCKS,NTHREADS >>>(dist,Den,rhoA,rhoB,Np);
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q19_GreyColorIMRT_Init: %s \n",cudaGetErrorString(err));
+	}
+}
+
+extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleColorDensity(int *NeighborList, double *Aq, double *Bq, double *Den, int start, int finish, int Np){
+
+	dvc_ScaLBL_D3Q7_AAodd_GreyscaleColorDensity<<<NBLOCKS,NTHREADS >>>(NeighborList, Aq, Bq, Den, start, finish, Np);
+
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q7_AAodd_GreyscaleColorDensity: %s \n",cudaGetErrorString(err));
+	}
+}
+
+extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleColorDensity(double *Aq, double *Bq, double *Den, int start, int finish, int Np){
+
+	dvc_ScaLBL_D3Q7_AAeven_GreyscaleColorDensity<<<NBLOCKS,NTHREADS >>>(Aq, Bq, Den, start, finish, Np);
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q7_AAeven_GreyscaleColorDensity: %s \n",cudaGetErrorString(err));
+	}
+}
+
+extern "C" void ScaLBL_D3Q19_GreyscaleColor_Gradient(int *neighborList, double *Den, double *DenGrad, int start, int finish, int Np){
+
+	dvc_ScaLBL_D3Q19_GreyscaleColor_Gradient<<<NBLOCKS,NTHREADS >>>(neighborList, Den, DenGrad, start, finish, Np);
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q19_GreyscaleColor_Gradient: %s \n",cudaGetErrorString(err));
+	}
+}
+
diff --git a/models/GreyscaleColorModel.cpp b/models/GreyscaleColorModel.cpp
index 11d92c80..421805c5 100644
--- a/models/GreyscaleColorModel.cpp
+++ b/models/GreyscaleColorModel.cpp
@@ -1,7 +1,7 @@
 /*
 Greyscale lattice boltzmann model
  */
-#include "models/GreyscaleModel.h"
+#include "models/GreyscaleColorModel.h"
 #include "analysis/distance.h"
 #include "analysis/morphology.h"
 #include <stdlib.h>
@@ -13,77 +13,90 @@ void DeleteArray( const TYPE *p )
     delete [] p;
 }
 
-ScaLBL_GreyscaleModel::ScaLBL_GreyscaleModel(int RANK, int NP, MPI_Comm COMM):
-rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tau(0),tau_eff(0),Den(0),Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),GreyPorosity(0),
+ScaLBL_GreyscaleColorModel::ScaLBL_GreyscaleColorModel(int RANK, int NP, MPI_Comm COMM):
+rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),tauA_eff(0),tauB_eff(0),rhoA(0),rhoB(0),Gsc(0),
+Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),GreyPorosity(0),
 Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
 {
 	SignDist.resize(Nx,Ny,Nz);           
     SignDist.fill(0);
 
 }
-ScaLBL_GreyscaleModel::~ScaLBL_GreyscaleModel(){
+ScaLBL_GreyscaleColorModel::~ScaLBL_GreyscaleColorModel(){
 
 }
 
-void ScaLBL_GreyscaleModel::ReadParams(string filename){
+void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
 	// read the input database 
 	db = std::make_shared<Database>( filename );
 	domain_db = db->getDatabase( "Domain" );
-	greyscale_db =  db->getDatabase( "Greyscale" );
+	greyscaleColor_db =  db->getDatabase( "GreyscaleColor" );
 	analysis_db = db->getDatabase( "Analysis" );
 	vis_db = db->getDatabase( "Visualization" );
 
 	// set defaults
 	timestepMax = 100000;
-	tau = 1.0;
-    tau_eff = tau;
-    Den = 1.0;//constant density
+	tauA = 1.0;
+	tauB = 1.0;
+    tauA_eff = tauA;//the effective viscosity of the Darcy term
+    tauB_eff = tauB;
+    rhoA = 1.0;//constant molecular mass (after LB scaling)
+    rhoB = 1.0;
 	tolerance = 0.01;
 	Fx = Fy = Fz = 0.0;
 	Restart=false;
 	din=dout=1.0;
 	flux=0.0;
     dp = 10.0; //unit of 'dp': voxel
-    CollisionType = 1; //1: IMRT; 2: BGK
 	
 	// ---------------------- Greyscale Model parameters -----------------------//
-	if (greyscale_db->keyExists( "timestepMax" )){
-		timestepMax = greyscale_db->getScalar<int>( "timestepMax" );
+	if (greyscaleColor_db->keyExists( "timestepMax" )){
+		timestepMax = greyscaleColor_db->getScalar<int>( "timestepMax" );
 	}
-	if (greyscale_db->keyExists( "tau" )){
-		tau = greyscale_db->getScalar<double>( "tau" );
+	if (greyscaleColor_db->keyExists( "tauA" )){
+		tauA = greyscaleColor_db->getScalar<double>( "tauA" );
 	}
-	tau_eff = greyscale_db->getWithDefault<double>( "tau_eff", tau );
-	if (greyscale_db->keyExists( "Den" )){
-		Den = greyscale_db->getScalar<double>( "Den" );
+	if (greyscaleColor_db->keyExists( "tauB" )){
+		tauB = greyscaleColor_db->getScalar<double>( "tauB" );
 	}
-	if (greyscale_db->keyExists( "dp" )){
-		dp = greyscale_db->getScalar<double>( "dp" );
+	tauA_eff = greyscaleColor_db->getWithDefault<double>( "tauA_eff", tauA);
+	tauB_eff = greyscaleColor_db->getWithDefault<double>( "tauB_eff", tauB);
+	if (greyscaleColor_db->keyExists( "rhoA" )){
+		rhoA = greyscaleColor_db->getScalar<double>( "rhoA" );
 	}
-	if (greyscale_db->keyExists( "F" )){
-		Fx = greyscale_db->getVector<double>( "F" )[0];
-		Fy = greyscale_db->getVector<double>( "F" )[1];
-		Fz = greyscale_db->getVector<double>( "F" )[2];
+	if (greyscaleColor_db->keyExists( "rhoB" )){
+		rhoB = greyscaleColor_db->getScalar<double>( "rhoB" );
 	}
-	if (greyscale_db->keyExists( "Restart" )){
-		Restart = greyscale_db->getScalar<bool>( "Restart" );
+	if (greyscaleColor_db->keyExists( "Gsc" )){
+		Gsc = greyscaleColor_db->getScalar<double>( "Gsc" );
 	}
-	if (greyscale_db->keyExists( "din" )){
-		din = greyscale_db->getScalar<double>( "din" );
+	if (greyscaleColor_db->keyExists( "dp" )){
+		dp = greyscaleColor_db->getScalar<double>( "dp" );
 	}
-	if (greyscale_db->keyExists( "dout" )){
-		dout = greyscale_db->getScalar<double>( "dout" );
+	if (greyscaleColor_db->keyExists( "F" )){
+		Fx = greyscaleColor_db->getVector<double>( "F" )[0];
+		Fy = greyscaleColor_db->getVector<double>( "F" )[1];
+		Fz = greyscaleColor_db->getVector<double>( "F" )[2];
 	}
-	if (greyscale_db->keyExists( "flux" )){
-		flux = greyscale_db->getScalar<double>( "flux" );
+	if (greyscaleColor_db->keyExists( "Restart" )){
+		Restart = greyscaleColor_db->getScalar<bool>( "Restart" );
 	}
-	if (greyscale_db->keyExists( "tolerance" )){
-		tolerance = greyscale_db->getScalar<double>( "tolerance" );
+	if (greyscaleColor_db->keyExists( "din" )){
+		din = greyscaleColor_db->getScalar<double>( "din" );
 	}
-	auto collision = greyscale_db->getWithDefault<std::string>( "collision", "IMRT" );
-	if (collision == "BGK"){
-        CollisionType=2;
+	if (greyscaleColor_db->keyExists( "dout" )){
+		dout = greyscaleColor_db->getScalar<double>( "dout" );
 	}
+	if (greyscaleColor_db->keyExists( "flux" )){
+		flux = greyscaleColor_db->getScalar<double>( "flux" );
+	}
+	if (greyscaleColor_db->keyExists( "tolerance" )){
+		tolerance = greyscaleColor_db->getScalar<double>( "tolerance" );
+	}
+	//auto collision = greyscaleColor_db->getWithDefault<std::string>( "collision", "IMRT" );
+	//if (collision == "BGK"){
+    //    CollisionType=2;
+	//}
 	// ------------------------------------------------------------------------//
     
     //------------------------ Other Domain parameters ------------------------//
@@ -94,7 +107,7 @@ void ScaLBL_GreyscaleModel::ReadParams(string filename){
 	// ------------------------------------------------------------------------//
 }
 
-void ScaLBL_GreyscaleModel::SetDomain(){
+void ScaLBL_GreyscaleColorModel::SetDomain(){
 	Dm  = std::shared_ptr<Domain>(new Domain(domain_db,comm));      // full domain for analysis
 	Mask  = std::shared_ptr<Domain>(new Domain(domain_db,comm));    // mask domain removes immobile phases
 	// domain parameters
@@ -125,7 +138,7 @@ void ScaLBL_GreyscaleModel::SetDomain(){
 	nprocz = Dm->nprocz();
 }
 
-void ScaLBL_GreyscaleModel::ReadInput(){
+void ScaLBL_GreyscaleColorModel::ReadInput(){
 	
 	sprintf(LocalRankString,"%05d",rank);
 	sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
@@ -174,44 +187,151 @@ void ScaLBL_GreyscaleModel::ReadInput(){
 	if (rank == 0) cout << "Domain set." << endl;
 }
 
-/********************************************************
- * AssignComponentLabels                                 *
- ********************************************************/
-void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Permeability)
+void ScaLBL_GreyscaleColorModel::AssignSolidForce(double *SolidPotential, double *SolidForce){
+
+	double *Dst;
+	Dst = new double [3*3*3];
+	for (int kk=0; kk<3; kk++){
+		for (int jj=0; jj<3; jj++){
+			for (int ii=0; ii<3; ii++){
+				int index = kk*9+jj*3+ii;
+				Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
+			}
+		}
+	}
+	double w_face = 1.f/18.f;
+	double w_edge = 1.f/36.f;
+	double w_corner = 0.f;
+	//local 
+	Dst[13] = 0.f;
+	//faces
+	Dst[4] = w_face;
+	Dst[10] = w_face;
+	Dst[12] = w_face;
+	Dst[14] = w_face;
+	Dst[16] = w_face;
+	Dst[22] = w_face;
+	// corners
+	Dst[0] = w_corner;
+	Dst[2] = w_corner;
+	Dst[6] = w_corner;
+	Dst[8] = w_corner;
+	Dst[18] = w_corner;
+	Dst[20] = w_corner;
+	Dst[24] = w_corner;
+	Dst[26] = w_corner;
+	// edges
+	Dst[1] = w_edge;
+	Dst[3] = w_edge;
+	Dst[5] = w_edge;
+	Dst[7] = w_edge;
+	Dst[9] = w_edge;
+	Dst[11] = w_edge;
+	Dst[15] = w_edge;
+	Dst[17] = w_edge;
+	Dst[19] = w_edge;
+	Dst[21] = w_edge;
+	Dst[23] = w_edge;
+	Dst[25] = w_edge;
+
+	for (int k=1; k<Nz-1; k++){
+		for (int j=1; j<Ny-1; j++){
+			for (int i=1; i<Nx-1; i++){
+				int idx=Map(i,j,k);
+				if (!(idx < 0)){
+					double phi_x = 0.f;
+					double phi_y = 0.f;
+					double phi_z = 0.f;
+					for (int kk=0; kk<3; kk++){
+						for (int jj=0; jj<3; jj++){
+							for (int ii=0; ii<3; ii++){
+
+								int index = kk*9+jj*3+ii;
+								double weight= Dst[index];
+
+								int idi=i+ii-1;
+								int idj=j+jj-1;
+								int idk=k+kk-1;
+
+								if (idi < 0) idi=0;
+								if (idj < 0) idj=0;
+								if (idk < 0) idk=0;
+								if (!(idi < Nx)) idi=Nx-1;
+								if (!(idj < Ny)) idj=Ny-1;
+								if (!(idk < Nz)) idk=Nz-1;
+
+								int nn = idk*Nx*Ny + idj*Nx + idi;
+								if ((Mask->id[nn] <= 0)||(Mask->id[nn]>=3)){
+									double vec_x = double(ii-1);
+									double vec_y = double(jj-1);
+									double vec_z = double(kk-1);
+									double GWNS=SolidPotential[nn];
+									phi_x += GWNS*weight*vec_x;
+									phi_y += GWNS*weight*vec_y;
+									phi_z += GWNS*weight*vec_z;
+								}
+							}
+						}
+					}
+					SolidForce[idx] = phi_x;
+					SolidForce[idx+Np] = phi_y;
+					SolidForce[idx+2*Np] = phi_z;
+				}
+			}
+		}
+	}
+	delete [] Dst;
+}
+
+
+void ScaLBL_GreyscaleColorModel::AssignComponentLabels(double *Porosity, double *Permeability, double *SolidPotential)
 {
 	size_t NLABELS=0;
 	signed char VALUE=0;
 	double POROSITY=0.f;
 	double PERMEABILITY=0.f;
+	double AFFINITY=0.f;
 
-	auto LabelList = greyscale_db->getVector<int>( "ComponentLabels" );
-	auto PorosityList = greyscale_db->getVector<double>( "PorosityList" );
-	auto PermeabilityList = greyscale_db->getVector<double>( "PermeabilityList" );
+	auto LabelList = greyscaleColor_db->getVector<int>( "ComponentLabels" );
+	auto AffinityList = greyscaleColor_db->getVector<double>( "ComponentAffinity" );
+	auto PorosityList = greyscaleColor_db->getVector<double>( "PorosityList" );
+	auto PermeabilityList = greyscaleColor_db->getVector<double>( "PermeabilityList" );
+
+    //1. Requirement for "ComponentLabels":
+    //   *labels can be a nagative integer, 0, 1, 2, or a positive integer >= 3
+    //   *label = 1 and 2 are reserved for NW and W phase respectively.
+    //2. Requirement for "ComponentAffinity":
+    //   *should be in the same length as "ComponentLabels"
+    //   *could leave Affinity=0.0 for label=1 and 2
+    //3. Requirement for "PorosityList":
+    //   *for ComponentLables <=0, put porosity value = 0.0;
+    //   *for ComponentLabels >=3, put the corresponding sub-resolution porosity
+    //   *for ComponentLabels =1, 2, put porosity=1 (or if users accidentally put other values it should still be fine)
+    //4. Requirement for "PermeabilityList":
+    //   *for ComponentLabels <=2, does not matter, can leave it as 1.0 
 
 	NLABELS=LabelList.size();
-	if (NLABELS != PorosityList.size()){
-		ERROR("Error: ComponentLabels and PorosityList must be the same length! \n");
+	if (NLABELS != PorosityList.size() || NLABELS != AffinityList.size() || NLABELS != PermeabilityList.size()){
+		ERROR("Error: ComponentLabels, ComponentAffinity, PorosityList and PermeabilityList must all be the same length! \n");
 	}
 
 	double label_count[NLABELS];
 	double label_count_global[NLABELS];
-	// Assign the labels
 
 	for (int idx=0; idx<NLABELS; idx++) label_count[idx]=0;
 
+    //Populate the poroisty map, NOTE only for node_ID > 0, i.e. open or grey nodes
+    //For node_ID <= 0: these are solid nodes of various wettability
 	for (int k=1;k<Nz-1;k++){
 		for (int j=1;j<Ny-1;j++){
 			for (int i=1;i<Nx-1;i++){
 				int n = k*Nx*Ny+j*Nx+i;
 				VALUE=id[n];
-				// Assign the affinity from the paired list
 				for (unsigned int idx=0; idx < NLABELS; idx++){
-				      //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
-					if (VALUE == LabelList[idx]){
+					if ((VALUE>0) && (VALUE == LabelList[idx])){
 						POROSITY=PorosityList[idx];
 						label_count[idx] += 1.0;
 						idx = NLABELS;
-						//Mask->id[n] = 0; // set mask to zero since this is an immobile component
 					}
 				}
 				int idx = Map(i,j,k);
@@ -227,9 +347,8 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
 		}
 	}
 
-	if (NLABELS != PermeabilityList.size()){
-		ERROR("Error: ComponentLabels and PermeabilityList must be the same length! \n");
-	}
+    //Populate the permeability map, NOTE only for node_ID > 0, i.e. open or grey nodes
+    //For node_ID <= 0: these are solid nodes of various wettability
 	for (int k=1;k<Nz-1;k++){
 		for (int j=1;j<Ny-1;j++){
 			for (int i=1;i<Nx-1;i++){
@@ -238,7 +357,7 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
 				// Assign the affinity from the paired list
 				for (unsigned int idx=0; idx < NLABELS; idx++){
 					//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
-					if (VALUE == LabelList[idx]){
+					if ( (VALUE>0) && (VALUE == LabelList[idx])){
 						PERMEABILITY=PermeabilityList[idx];
 						idx = NLABELS;
 						//Mask->id[n] = 0; // set mask to zero since this is an immobile component
@@ -257,6 +376,34 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
 		}
 	}
 
+    //Populate the solid potential map, for ALL range of node_ID except node = 1,2, i.e. NW and W phase
+	for (int k=0;k<Nz;k++){
+		for (int j=0;j<Ny;j++){
+			for (int i=0;i<Nx;i++){
+				int n = k*Nx*Ny+j*Nx+i;
+				VALUE=id[n];
+				// Assign the affinity from the paired list
+				for (unsigned int idx=0; idx < NLABELS; idx++){
+					if (VALUE == LabelList[idx]){
+                        if (VALUE<=0){
+                            AFFINITY=AffinityList[idx];
+                        }
+                        else if (VALUE>=3){
+						    AFFINITY=AffinityList[idx]*(1.0-PorosityList[idx]);//BE CAREFUL! Requires for node_ID<=0, user puts porosity=0.0
+                        }
+                        else{//i.e. label = 1 or 2
+                            AFFINITY=0.0;
+                        } 
+						idx = NLABELS;
+					}
+				}
+				//NOTE: node_ID = 1 and 2 are reserved
+				if ((VALUE == 1)||(VALUE == 2)) AFFINITY=0.0;//NOTE: still need this as users may forget to put label=1,2 in ComponentLabelLists
+				SolidPotential[n] = AFFINITY;
+			}
+		}
+	}
+
 
 	// Set Dm to match Mask
 	for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i]; 
@@ -286,8 +433,61 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
 	}
 }
 
+void ScaLBL_GreyscaleColorModel::DensityField_Init(){
 
-void ScaLBL_GreyscaleModel::Create(){
+	size_t NLABELS=0;
+	signed char VALUE=0;
+    //Read all greyscale node labels and their conrresponding initial saturaiton
+    //TODO make something default, LabelList=1,2, SwList=0.0, 1.0
+	auto LabelList = greyscaleColor_db->getVector<int>( "GreyNodeLabels" );
+	auto SwList = greyscaleColor_db->getVector<double>( "GreyNodeSw" );
+	NLABELS=LabelList.size();
+	if (NLABELS != SwList.size()){
+		ERROR("Error: GreyNodeLabels and GreyNodeSw must be the same length! \n");
+	}
+	
+    double *Den_temp;
+	Den_temp=new double [2*Np];
+    double nA=0.5;//to prevent use may forget to specify all greynodes, then must initialize something to start with, givning just zeros is too risky.
+    double nB=0.5;
+
+	for (int k=1; k<Nz-1; k++){
+		for (int j=1; j<Ny-1; j++){
+			for (int i=1; i<Nx-1; i++){
+				int n = k*Nx*Ny+j*Nx+i;
+				VALUE=Mask->id[n];
+                if (VALUE>0){
+                    for (unsigned int idx=0; idx < NLABELS; idx++){
+                        if (VALUE == LabelList[idx]){
+                            double Sw = SwList[idx];
+                            if ((Sw<0.0) || (Sw>1.0)) ERROR("Error: Initial saturation for grey nodes must be between [0.0, 1.0]! \n");
+                            nB=Sw;
+                            nA=1.0-Sw;
+                            idx = NLABELS;
+                        }
+                    }
+                    if (VALUE==1){//label=1 reserved for NW phase
+                        nA=1.0;
+                        nB=0.0; 
+                    }
+                    else if(VALUE==2){//label=2 reserved for W phase
+                        nA=0.0;
+                        nB=1.0; 
+                    }
+                    int idx = Map(i,j,k);
+                    Den_temp[idx+0*Np] = nA;
+                    Den_temp[idx+1*Np] = nB;
+                }
+			}
+		}
+	}
+    //copy to device
+	ScaLBL_CopyToDevice(Den, Den_temp, 2*Np*sizeof(double));
+	ScaLBL_DeviceBarrier();
+	delete [] Den_temp;
+}
+
+void ScaLBL_GreyscaleColorModel::Create(){
 	/*
 	 *  This function creates the variables needed to run a LBM 
 	 */
@@ -324,99 +524,81 @@ void ScaLBL_GreyscaleModel::Create(){
 	neighborSize=18*(Np*sizeof(int));
 	//...........................................................................
 	ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
-	ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np);
 	ScaLBL_AllocateDeviceMemory((void **) &fq, 19*dist_mem_size);
 	ScaLBL_AllocateDeviceMemory((void **) &Permeability, sizeof(double)*Np);		
 	ScaLBL_AllocateDeviceMemory((void **) &Porosity, sizeof(double)*Np);		
 	ScaLBL_AllocateDeviceMemory((void **) &Pressure_dvc, sizeof(double)*Np);
 	ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &Aq, 7*sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &Bq, 7*sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &Den, 2*sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &SolidForce, 3*sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &DenGradA, 3*sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &DenGradB, 3*sizeof(double)*Np);
 	//...........................................................................
 	// Update GPU data structures
-	if (rank==0)	printf ("Setting up device map and neighbor list \n");
+	if (rank==0)	printf ("Setting up device neighbor list \n");
 	fflush(stdout);
-	int *TmpMap;
-	TmpMap=new int[Np];
-	for (int k=1; k<Nz-1; k++){
-		for (int j=1; j<Ny-1; j++){
-			for (int i=1; i<Nx-1; i++){
-				int idx=Map(i,j,k);
-				if (!(idx < 0))
-					TmpMap[idx] = k*Nx*Ny+j*Nx+i;
-			}
-		}
-	}
-	// check that TmpMap is valid
-	for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
-		int n = TmpMap[idx];
-		if (n > Nx*Ny*Nz){
-			printf("Bad value! idx=%i \n");
-			TmpMap[idx] = Nx*Ny*Nz-1;
-		}
-	}
-	for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
-		int n = TmpMap[idx];
-		if (n > Nx*Ny*Nz){
-			printf("Bad value! idx=%i \n");
-			TmpMap[idx] = Nx*Ny*Nz-1;
-		}
-	}
-	ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
-	ScaLBL_DeviceBarrier();
-	delete [] TmpMap;
-	
 	// copy the neighbor list 
 	ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
 	// initialize phi based on PhaseLabel (include solid component labels)
 	double *Poros, *Perm;
 	Poros = new double[Np];
 	Perm = new double[Np];
-	AssignComponentLabels(Poros,Perm);
+	double *SolidForce_host = new double[3*Np];
+    double *SolidPotential_host = new double [Nx*Ny*Nz];
+	AssignComponentLabels(Poros,Perm,SolidPotential_host);
+    AssignSolidForce(SolidPotential_host,SolidForce_host);
 	ScaLBL_CopyToDevice(Porosity, Poros, Np*sizeof(double));
 	ScaLBL_CopyToDevice(Permeability, Perm, Np*sizeof(double));
+	ScaLBL_CopyToDevice(SolidForce, SolidForce_host, 3*Np*sizeof(double));
+
+	ScaLBL_DeviceBarrier();
+    //TODO make the following smart pointers
+	delete [] SolidForce_host;
+	delete [] SolidPotential_host;
+    delete [] Poros;
+    delete [] Perm;
 }        
 
 
-void ScaLBL_GreyscaleModel::Initialize(){
-	if (rank==0)	printf ("Initializing distributions \n");
-    //TODO: for BGK, you need to consider voxel porosity
-    //      for IMRT, the whole set of feq is different
-    //      if in the future you have different collison mode, need to write two set of initialization functions
-    if (CollisionType==1){
-	    ScaLBL_D3Q19_GreyIMRT_Init(fq, Np, Den);
-        if (rank==0) printf("Collision model: Incompressible MRT.\n");
-    }
-    else if (CollisionType==2){
-	    ScaLBL_D3Q19_Init(fq, Np);
-        if (rank==0) printf("Collision model: BGK.\n");
-    }
-    else{
-        if (rank==0) printf("Unknown collison type! IMRT collision is used.\n"); 
-	    ScaLBL_D3Q19_GreyIMRT_Init(fq, Np, Den);
-        CollisionType=1;
-		greyscale_db->putScalar<std::string>( "collision", "IMRT" );
-    }
-
+void ScaLBL_GreyscaleColorModel::Initialize(){
 	if (Restart == true){
 		if (rank==0){
-			printf("Initializing distributions from Restart! \n");
+			printf("Initializing density field and distributions from Restart! \n");
 		}
 		// Read in the restart file to CPU buffers
         std::shared_ptr<double> cfq;
         cfq = std::shared_ptr<double>(new double[19*Np],DeleteArray<double>);
+        std::shared_ptr<double> cDen;
+        cDen = std::shared_ptr<double>(new double[2*Np],DeleteArray<double>);
         FILE *File;
         File=fopen(LocalRestartFile,"rb");
         fread(cfq.get(),sizeof(double),19*Np,File);
+        fread(cDen.get(),sizeof(double),2*Np,File);
         fclose(File);
 
 		// Copy the restart data to the GPU
 		ScaLBL_CopyToDevice(fq,cfq.get(),19*Np*sizeof(double));
+		ScaLBL_CopyToDevice(Den,cDen.get(),2*Np*sizeof(double));
 		ScaLBL_DeviceBarrier();
-
 		MPI_Barrier(comm);
+
+        ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);//initialize D3Q7 density components
+        ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 	}
+    else{
+        if (rank==0)	printf ("Initializing density field \n");
+        DensityField_Init();//initialize density field
+        ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);//initialize D3Q7 density components
+        ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+        
+        if (rank==0)	printf ("Initializing distributions \n");
+        ScaLBL_D3Q19_GreyColorIMRT_Init(fq, Den, rhoA, rhoB, Np);
+    }
 }
 
-void ScaLBL_GreyscaleModel::Run(){
+void ScaLBL_GreyscaleColorModel::Run(){
 	int nprocs=nprocx*nprocy*nprocz;
 	const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
 	
@@ -432,8 +614,8 @@ void ScaLBL_GreyscaleModel::Run(){
 	if (analysis_db->keyExists( "restart_interval" )){
 		restart_interval = analysis_db->getScalar<int>( "restart_interval" );
 	}
-	if (greyscale_db->keyExists( "timestep" )){
-		timestep = greyscale_db->getScalar<int>( "timestep" );
+	if (greyscaleColor_db->keyExists( "timestep" )){
+		timestep = greyscaleColor_db->getScalar<int>( "timestep" );
 	}
 
 	if (rank==0){
@@ -454,213 +636,217 @@ void ScaLBL_GreyscaleModel::Run(){
 	//************ MAIN ITERATION LOOP ***************************************/
 	PROFILE_START("Loop");
 	auto current_db = db->cloneDatabase();
-	double rlx = 1.0/tau;
-    double rlx_eff = 1.0/tau_eff;
 	double error = 1.0;
 	double flow_rate_previous = 0.0;
 	while (timestep < timestepMax && error > tolerance) {
 		//************************************************************************/
 		// *************ODD TIMESTEP*************//
 		timestep++;
+		// Compute the density field
+		// Read for Aq, Bq happens in this routine (requires communication)
+		ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
+		ScaLBL_D3Q7_AAodd_GreyscaleColorDensity(NeighborList, Aq, Bq, Den, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
+		ScaLBL_DeviceBarrier();
+		ScaLBL_D3Q7_AAodd_GreyscaleColorDensity(NeighborList, Aq, Bq, Den, 0, ScaLBL_Comm->LastExterior(), Np);
+
+        //compute Den gradients - component A
+		ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[0*Np], DenGradA, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_Comm->SendHalo(&Den[0*Np]);
+		ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[0*Np], DenGradA, 0, ScaLBL_Comm->LastExterior(), Np);
+		ScaLBL_Comm->RecvGrad(&Den[0*Np],DenGradA);
+        //compute Den gradients - component B
+		ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[1*Np], DenGradB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_Comm->SendHalo(&Den[1*Np]);
+		ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[1*Np], DenGradB, 0, ScaLBL_Comm->LastExterior(), Np);
+		ScaLBL_Comm->RecvGrad(&Den[1*Np],DenGradB);
+
 		ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
-        switch (CollisionType){
-            case 1: 
-                    ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq,  ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
-                    break;
-            case 2: 
-                    ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq,  ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
-                    break;
-            default: 
-                    ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq,  ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
-                    break;
-        }
+        ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, fq, Aq, Bq, Den, DenGradA, DenGradB, SolidForce, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, 
+                tauA,tauB,tauA_eff,tauB_eff,rhoA,rhoB,Gsc,Fx,Fy,Fz,Porosity,Permeability,Velocity,Pressure_dvc);
 		ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
 		ScaLBL_DeviceBarrier();
-		// Set BCs
-		if (BoundaryCondition == 3){
-			ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
-			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
-		}
-        switch (CollisionType){
-            case 1: 
-		            ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
-                    break;
-            case 2: 
-		            ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
-                    break;
-            default: 
-		            ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
-                    break;
-        }
-		ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+
+//		// Set BCs
+//		if (BoundaryCondition == 3){
+//			ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
+//			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+//		}
+        ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, fq, Aq, Bq, Den, DenGradA, DenGradB, SolidForce, 0, ScaLBL_Comm->LastExterior(), Np, 
+                tauA,tauB,tauA_eff,tauB_eff,rhoA,rhoB,Gsc,Fx,Fy,Fz,Porosity,Permeability,Velocity,Pressure_dvc);
+		ScaLBL_DeviceBarrier(); 
+        MPI_Barrier(comm);
 
 		// *************EVEN TIMESTEP*************//
 		timestep++;
-		ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
-        switch (CollisionType){
-            case 1: 
-		            ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
-                    break;
-            case 2: 
-		            ScaLBL_D3Q19_AAeven_Greyscale(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
-                    break;
-            default: 
-		            ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
-                    break;
-        }
-        ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
+		// Compute the density field
+		// Read for Aq, Bq happens in this routine (requires communication)
+		ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
+		ScaLBL_D3Q7_AAeven_GreyscaleColorDensity(Aq, Bq, Den, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
 		ScaLBL_DeviceBarrier();
-		// Set BCs
-		if (BoundaryCondition == 3){
-			ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
-			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
-		}
-        switch (CollisionType){
-            case 1: 
-		            ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
-                    break;
-            case 2: 
-		            ScaLBL_D3Q19_AAeven_Greyscale(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
-                    break;
-            default: 
-		            ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
-                    break;
-        }
-        ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+		ScaLBL_D3Q7_AAeven_GreyscaleColorDensity(Aq, Bq, Den, 0, ScaLBL_Comm->LastExterior(), Np);
+
+        //compute Den gradients - component A
+		ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[0*Np], DenGradA, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_Comm->SendHalo(&Den[0*Np]);
+		ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[0*Np], DenGradA, 0, ScaLBL_Comm->LastExterior(), Np);
+		ScaLBL_Comm->RecvGrad(&Den[0*Np],DenGradA);
+        //compute Den gradients - component B
+		ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[1*Np], DenGradB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_Comm->SendHalo(&Den[1*Np]);
+		ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[1*Np], DenGradB, 0, ScaLBL_Comm->LastExterior(), Np);
+		ScaLBL_Comm->RecvGrad(&Den[1*Np],DenGradB);
+
+		ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
+        ScaLBL_D3Q19_AAeven_GreyscaleColor(fq, Aq, Bq, Den, DenGradA, DenGradB, SolidForce, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, 
+                tauA,tauB,tauA_eff,tauB_eff,rhoA,rhoB,Gsc,Fx,Fy,Fz,Porosity,Permeability,Velocity,Pressure_dvc);
+		ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
+		ScaLBL_DeviceBarrier();
+
+//		// Set BCs
+//		if (BoundaryCondition == 3){
+//			ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
+//			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+//		}
+        ScaLBL_D3Q19_AAeven_GreyscaleColor(fq, Aq, Bq, Den, DenGradA, DenGradB, SolidForce, 0, ScaLBL_Comm->LastExterior(), Np, 
+                tauA,tauB,tauA_eff,tauB_eff,rhoA,rhoB,Gsc,Fx,Fy,Fz,Porosity,Permeability,Velocity,Pressure_dvc);
+        ScaLBL_DeviceBarrier(); 
+        MPI_Barrier(comm);
 		//************************************************************************/
 		
-		if (timestep%analysis_interval==0){
-			ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Velocity_x);
-			ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Velocity_y);
-			ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Velocity_z);
-			//ScaLBL_Comm->RegularLayout(Map,Porosity,PorosityMap);
-			//ScaLBL_Comm->RegularLayout(Map,Pressure_dvc,Pressure);
-			
-			double count_loc=0;
-			double count;
-			double vax,vay,vaz;
-			double vax_loc,vay_loc,vaz_loc;
-            //double px_loc,py_loc,pz_loc;
-            //double px,py,pz;
-            //double mass_loc,mass_glb;
-            
-            //parameters for domain average
-	        int64_t i,j,k,n,imin,jmin,kmin,kmax;
-            // If external boundary conditions are set, do not average over the inlet and outlet
-            kmin=1; kmax=Nz-1;
-            //In case user forgets to specify the inlet/outlet buffer layers for BC>0
-            if (BoundaryCondition > 0 && Dm->kproc() == 0) kmin=4;
-            if (BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
-
-            imin=jmin=1;
-            // If inlet/outlet layers exist use these as default
-            //if (Dm->inlet_layers_x > 0) imin = Dm->inlet_layers_x;
-            //if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y;
-            if (BoundaryCondition > 0 && Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin = 1 + Dm->inlet_layers_z;//"1" indicates the halo layer
-            if (BoundaryCondition > 0 && Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax = Nz-1 - Dm->outlet_layers_z; 
-
-//			px_loc = py_loc = pz_loc = 0.f;
-//            mass_loc = 0.f;
+//		if (timestep%analysis_interval==0){
+//			ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Velocity_x);
+//			ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Velocity_y);
+//			ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Velocity_z);
+//			//ScaLBL_Comm->RegularLayout(Map,Porosity,PorosityMap);
+//			//ScaLBL_Comm->RegularLayout(Map,Pressure_dvc,Pressure);
+//			
+//			double count_loc=0;
+//			double count;
+//			double vax,vay,vaz;
+//			double vax_loc,vay_loc,vaz_loc;
+//            //double px_loc,py_loc,pz_loc;
+//            //double px,py,pz;
+//            //double mass_loc,mass_glb;
+//            
+//            //parameters for domain average
+//	        int64_t i,j,k,n,imin,jmin,kmin,kmax;
+//            // If external boundary conditions are set, do not average over the inlet and outlet
+//            kmin=1; kmax=Nz-1;
+//            //In case user forgets to specify the inlet/outlet buffer layers for BC>0
+//            if (BoundaryCondition > 0 && Dm->kproc() == 0) kmin=4;
+//            if (BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
+//
+//            imin=jmin=1;
+//            // If inlet/outlet layers exist use these as default
+//            //if (Dm->inlet_layers_x > 0) imin = Dm->inlet_layers_x;
+//            //if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y;
+//            if (BoundaryCondition > 0 && Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin = 1 + Dm->inlet_layers_z;//"1" indicates the halo layer
+//            if (BoundaryCondition > 0 && Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax = Nz-1 - Dm->outlet_layers_z; 
+//
+////			px_loc = py_loc = pz_loc = 0.f;
+////            mass_loc = 0.f;
+////			for (int k=kmin; k<kmax; k++){
+////				for (int j=jmin; j<Ny-1; j++){
+////					for (int i=imin; i<Nx-1; i++){
+////						if (SignDist(i,j,k) > 0){
+////							px_loc   += Velocity_x(i,j,k)*Den*PorosityMap(i,j,k);
+////							py_loc   += Velocity_y(i,j,k)*Den*PorosityMap(i,j,k);
+////							pz_loc   += Velocity_z(i,j,k)*Den*PorosityMap(i,j,k);
+////							mass_loc += Den*PorosityMap(i,j,k);
+////						}
+////					}
+////				}
+////			}
+////			MPI_Allreduce(&px_loc,  &px,      1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+////			MPI_Allreduce(&py_loc,  &py,      1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+////			MPI_Allreduce(&pz_loc,  &pz,      1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+////			MPI_Allreduce(&mass_loc,&mass_glb,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+////			
+////			vax = px/mass_glb;
+////			vay = py/mass_glb;
+////			vaz = pz/mass_glb;
+//            
+//			vax_loc = vay_loc = vaz_loc = 0.f;
 //			for (int k=kmin; k<kmax; k++){
 //				for (int j=jmin; j<Ny-1; j++){
 //					for (int i=imin; i<Nx-1; i++){
 //						if (SignDist(i,j,k) > 0){
-//							px_loc   += Velocity_x(i,j,k)*Den*PorosityMap(i,j,k);
-//							py_loc   += Velocity_y(i,j,k)*Den*PorosityMap(i,j,k);
-//							pz_loc   += Velocity_z(i,j,k)*Den*PorosityMap(i,j,k);
-//							mass_loc += Den*PorosityMap(i,j,k);
+//							vax_loc += Velocity_x(i,j,k);
+//							vay_loc += Velocity_y(i,j,k);
+//							vaz_loc += Velocity_z(i,j,k);
+//							count_loc+=1.0;
 //						}
 //					}
 //				}
 //			}
-//			MPI_Allreduce(&px_loc,  &px,      1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
-//			MPI_Allreduce(&py_loc,  &py,      1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
-//			MPI_Allreduce(&pz_loc,  &pz,      1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
-//			MPI_Allreduce(&mass_loc,&mass_glb,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+//			//MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+//			//MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+//			//MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+//			//MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
 //			
-//			vax = px/mass_glb;
-//			vay = py/mass_glb;
-//			vaz = pz/mass_glb;
-            
-			vax_loc = vay_loc = vaz_loc = 0.f;
-			for (int k=kmin; k<kmax; k++){
-				for (int j=jmin; j<Ny-1; j++){
-					for (int i=imin; i<Nx-1; i++){
-						if (SignDist(i,j,k) > 0){
-							vax_loc += Velocity_x(i,j,k);
-							vay_loc += Velocity_y(i,j,k);
-							vaz_loc += Velocity_z(i,j,k);
-							count_loc+=1.0;
-						}
-					}
-				}
-			}
-			//MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
-			//MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
-			//MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
-			//MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
-			
-            vax = Mask->Comm.sumReduce( vax_loc );
-            vay = Mask->Comm.sumReduce( vay_loc );
-            vaz = Mask->Comm.sumReduce( vaz_loc );
-            count = Mask->Comm.sumReduce( count_loc );
-
-			vax /= count;
-			vay /= count;
-			vaz /= count;
-
-			double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
-			double dir_x = Fx/force_mag;
-			double dir_y = Fy/force_mag;
-			double dir_z = Fz/force_mag;
-			if (force_mag == 0.0){
-				// default to z direction
-				dir_x = 0.0;
-				dir_y = 0.0;
-				dir_z = 1.0;
-				force_mag = 1.0;
-			}
-			//double flow_rate = (px*dir_x + py*dir_y + pz*dir_z)/mass_glb;
-			double flow_rate = (vax*dir_x + vay*dir_y + vaz*dir_z);
-			
-			error = fabs(flow_rate - flow_rate_previous) / fabs(flow_rate);
-			flow_rate_previous = flow_rate;
-			
-			//if (rank==0) printf("Computing Minkowski functionals \n");
-			Morphology.ComputeScalar(SignDist,0.f);
-			//Morphology.PrintAll();
-			double mu = (tau-0.5)/3.f;
-			double Vs = Morphology.V();
-			double As = Morphology.A();
-			double Hs = Morphology.H();
-			double Xs = Morphology.X();
-			Vs = Dm->Comm.sumReduce( Vs);
-			As = Dm->Comm.sumReduce( As);
-			Hs = Dm->Comm.sumReduce( Hs);
-			Xs = Dm->Comm.sumReduce( Xs);
-
-			double h = Dm->voxel_length;
-			//double absperm = h*h*mu*Mask->Porosity()*flow_rate / force_mag;
-			double absperm = h*h*mu*GreyPorosity*flow_rate / force_mag;
-
-            if (rank==0){
-				printf("     AbsPerm = %.5g [micron^2]\n",absperm);
-                bool WriteHeader=false;
-                FILE * log_file = fopen("Permeability.csv","r");
-                if (log_file != NULL)
-                    fclose(log_file);
-                else
-                    WriteHeader=true;
-                log_file = fopen("Permeability.csv","a");
-                if (WriteHeader)
-                    fprintf(log_file,"timestep Fx Fy Fz mu Vs As Hs Xs vax vay vaz AbsPerm \n",
-                            timestep,Fx,Fy,Fz,mu,h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz,absperm);
-
-                fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, 
-                        h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm);
-                fclose(log_file);
-            }
-		}
+//            vax = Mask->Comm.sumReduce( vax_loc );
+//            vay = Mask->Comm.sumReduce( vay_loc );
+//            vaz = Mask->Comm.sumReduce( vaz_loc );
+//            count = Mask->Comm.sumReduce( count_loc );
+//
+//			vax /= count;
+//			vay /= count;
+//			vaz /= count;
+//
+//			double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
+//			double dir_x = Fx/force_mag;
+//			double dir_y = Fy/force_mag;
+//			double dir_z = Fz/force_mag;
+//			if (force_mag == 0.0){
+//				// default to z direction
+//				dir_x = 0.0;
+//				dir_y = 0.0;
+//				dir_z = 1.0;
+//				force_mag = 1.0;
+//			}
+//			//double flow_rate = (px*dir_x + py*dir_y + pz*dir_z)/mass_glb;
+//			double flow_rate = (vax*dir_x + vay*dir_y + vaz*dir_z);
+//			
+//			error = fabs(flow_rate - flow_rate_previous) / fabs(flow_rate);
+//			flow_rate_previous = flow_rate;
+//			
+//			//if (rank==0) printf("Computing Minkowski functionals \n");
+//			Morphology.ComputeScalar(SignDist,0.f);
+//			//Morphology.PrintAll();
+//			double mu = (tau-0.5)/3.f;
+//			double Vs = Morphology.V();
+//			double As = Morphology.A();
+//			double Hs = Morphology.H();
+//			double Xs = Morphology.X();
+//			Vs = Dm->Comm.sumReduce( Vs);
+//			As = Dm->Comm.sumReduce( As);
+//			Hs = Dm->Comm.sumReduce( Hs);
+//			Xs = Dm->Comm.sumReduce( Xs);
+//
+//			double h = Dm->voxel_length;
+//			//double absperm = h*h*mu*Mask->Porosity()*flow_rate / force_mag;
+//			double absperm = h*h*mu*GreyPorosity*flow_rate / force_mag;
+//
+//            if (rank==0){
+//				printf("     AbsPerm = %.5g [micron^2]\n",absperm);
+//                bool WriteHeader=false;
+//                FILE * log_file = fopen("Permeability.csv","r");
+//                if (log_file != NULL)
+//                    fclose(log_file);
+//                else
+//                    WriteHeader=true;
+//                log_file = fopen("Permeability.csv","a");
+//                if (WriteHeader)
+//                    fprintf(log_file,"timestep Fx Fy Fz mu Vs As Hs Xs vax vay vaz AbsPerm \n",
+//                            timestep,Fx,Fy,Fz,mu,h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz,absperm);
+//
+//                fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, 
+//                        h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm);
+//                fclose(log_file);
+//            }
+//		}
 
 		if (timestep%visualization_interval==0){
             VelocityField();
@@ -669,9 +855,9 @@ void ScaLBL_GreyscaleModel::Run(){
 		if (timestep%restart_interval==0){
             //Use rank=0 write out Restart.db
             if (rank==0) {
-                greyscale_db->putScalar<int>("timestep",timestep);    		
-                greyscale_db->putScalar<bool>( "Restart", true );
-                current_db->putDatabase("Greyscale", greyscale_db);
+                greyscaleColor_db->putScalar<int>("timestep",timestep);    		
+                greyscaleColor_db->putScalar<bool>( "Restart", true );
+                current_db->putDatabase("GreyscaleColor", greyscaleColor_db);
                 std::ofstream OutStream("Restart.db");
                 current_db->print(OutStream, "");
                 OutStream.close();
@@ -681,17 +867,21 @@ void ScaLBL_GreyscaleModel::Run(){
             std::shared_ptr<double> cfq;
             cfq = std::shared_ptr<double>(new double[19*Np],DeleteArray<double>);
             ScaLBL_CopyToHost(cfq.get(),fq,19*Np*sizeof(double));// Copy restart data to the CPU
+            std::shared_ptr<double> cDen;
+            cDen = std::shared_ptr<double>(new double[2*Np],DeleteArray<double>);
+            ScaLBL_CopyToHost(cDen.get(),Den,2*Np*sizeof(double));// Copy restart data to the CPU
 
             FILE *RESTARTFILE;
             RESTARTFILE=fopen(LocalRestartFile,"wb");
             fwrite(cfq.get(),sizeof(double),19*Np,RESTARTFILE);
+            fwrite(cDen.get(),sizeof(double),2*Np,RESTARTFILE);
             fclose(RESTARTFILE);
 		    MPI_Barrier(comm);
         }
 	}
 
 	PROFILE_STOP("Loop");
-	PROFILE_SAVE("lbpm_greyscale_simulator",1);
+	PROFILE_SAVE("lbpm_greyscaleColor_simulator",1);
 	//************************************************************************
 	ScaLBL_DeviceBarrier();
 	MPI_Barrier(comm);
@@ -712,7 +902,7 @@ void ScaLBL_GreyscaleModel::Run(){
 	// ************************************************************************
 }
 
-void ScaLBL_GreyscaleModel::VelocityField(){
+void ScaLBL_GreyscaleColorModel::VelocityField(){
 
 /*	Minkowski Morphology(Mask);
 	int SIZE=Np*sizeof(double);
@@ -828,7 +1018,7 @@ void ScaLBL_GreyscaleModel::VelocityField(){
 
 }
 
-void ScaLBL_GreyscaleModel::WriteDebug(){
+void ScaLBL_GreyscaleColorModel::WriteDebug(){
 	// Copy back final phase indicator field and convert to regular layout
 	DoubleArray PhaseField(Nx,Ny,Nz);
 
@@ -840,26 +1030,26 @@ void ScaLBL_GreyscaleModel::WriteDebug(){
 //	fwrite(PhaseField.data(),8,N,OUTFILE);
 //	fclose(OUTFILE);
 //
-//    ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField);
-//	FILE *AFILE;
-//	sprintf(LocalRankFilename,"A.%05i.raw",rank);
-//	AFILE = fopen(LocalRankFilename,"wb");
-//	fwrite(PhaseField.data(),8,N,AFILE);
-//	fclose(AFILE);
-//
-//	ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField);
-//	FILE *BFILE;
-//	sprintf(LocalRankFilename,"B.%05i.raw",rank);
-//	BFILE = fopen(LocalRankFilename,"wb");
-//	fwrite(PhaseField.data(),8,N,BFILE);
-//	fclose(BFILE);
-//
-//	ScaLBL_Comm->RegularLayout(Map,Pressure,PhaseField);
-//	FILE *PFILE;
-//	sprintf(LocalRankFilename,"Pressure.%05i.raw",rank);
-//	PFILE = fopen(LocalRankFilename,"wb");
-//	fwrite(PhaseField.data(),8,N,PFILE);
-//	fclose(PFILE);
+    ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField);
+	FILE *AFILE;
+	sprintf(LocalRankFilename,"A.%05i.raw",rank);
+	AFILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,AFILE);
+	fclose(AFILE);
+
+	ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField);
+	FILE *BFILE;
+	sprintf(LocalRankFilename,"B.%05i.raw",rank);
+	BFILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,BFILE);
+	fclose(BFILE);
+
+	ScaLBL_Comm->RegularLayout(Map,Pressure_dvc,PhaseField);
+	FILE *PFILE;
+	sprintf(LocalRankFilename,"Pressure.%05i.raw",rank);
+	PFILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,PFILE);
+	fclose(PFILE);
 
 	ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
 	FILE *VELX_FILE;
@@ -882,17 +1072,17 @@ void ScaLBL_GreyscaleModel::WriteDebug(){
 	fwrite(PhaseField.data(),8,N,VELZ_FILE);
 	fclose(VELZ_FILE);
 
-	ScaLBL_Comm->RegularLayout(Map,&Porosity[0],PhaseField);
-	FILE *POROS_FILE;
-	sprintf(LocalRankFilename,"Porosity.%05i.raw",rank);
-	POROS_FILE = fopen(LocalRankFilename,"wb");
-	fwrite(PhaseField.data(),8,N,POROS_FILE);
-	fclose(POROS_FILE);
-
-	ScaLBL_Comm->RegularLayout(Map,&Permeability[0],PhaseField);
-	FILE *PERM_FILE;
-	sprintf(LocalRankFilename,"Permeability.%05i.raw",rank);
-	PERM_FILE = fopen(LocalRankFilename,"wb");
-	fwrite(PhaseField.data(),8,N,PERM_FILE);
-	fclose(PERM_FILE);
+//	ScaLBL_Comm->RegularLayout(Map,&Porosity[0],PhaseField);
+//	FILE *POROS_FILE;
+//	sprintf(LocalRankFilename,"Porosity.%05i.raw",rank);
+//	POROS_FILE = fopen(LocalRankFilename,"wb");
+//	fwrite(PhaseField.data(),8,N,POROS_FILE);
+//	fclose(POROS_FILE);
+//
+//	ScaLBL_Comm->RegularLayout(Map,&Permeability[0],PhaseField);
+//	FILE *PERM_FILE;
+//	sprintf(LocalRankFilename,"Permeability.%05i.raw",rank);
+//	PERM_FILE = fopen(LocalRankFilename,"wb");
+//	fwrite(PhaseField.data(),8,N,PERM_FILE);
+//	fclose(PERM_FILE);
 }
diff --git a/models/GreyscaleColorModel.h b/models/GreyscaleColorModel.h
index a99925b1..52841c4a 100644
--- a/models/GreyscaleColorModel.h
+++ b/models/GreyscaleColorModel.h
@@ -1,5 +1,5 @@
 /*
-Implementation of color lattice boltzmann model
+Implementation of multicomponent greyscale lattice boltzmann model
  */
 #include <stdio.h>
 #include <stdlib.h>
@@ -16,10 +16,10 @@ Implementation of color lattice boltzmann model
 #include "ProfilerApp.h"
 #include "threadpool/thread_pool.h"
 
-class ScaLBL_GreyscaleModel{
+class ScaLBL_GreyscaleColorModel{
 public:
-	ScaLBL_GreyscaleModel(int RANK, int NP, MPI_Comm COMM);
-	~ScaLBL_GreyscaleModel();	
+	ScaLBL_GreyscaleColorModel(int RANK, int NP, MPI_Comm COMM);
+	~ScaLBL_GreyscaleColorModel();	
 	
 	// functions in they should be run
 	void ReadParams(string filename);
@@ -35,15 +35,15 @@ public:
 	bool Restart,pBC;
 	int timestep,timestepMax;
 	int BoundaryCondition;
-    int CollisionType;
-	double tau;
-    double tau_eff;
-    double Den;//constant density
+	double tauA,tauB;
+    double tauA_eff,tauB_eff;
+    double rhoA,rhoB;
 	double tolerance;
 	double Fx,Fy,Fz,flux;
 	double din,dout;
     double dp;//solid particle diameter, unit in voxel
     double GreyPorosity;
+    double Gsc;
 	
 	int Nx,Ny,Nz,N,Np;
 	int rank,nprocx,nprocy,nprocz,nprocs;
@@ -56,18 +56,21 @@ public:
     // input database
     std::shared_ptr<Database> db;
     std::shared_ptr<Database> domain_db;
-    std::shared_ptr<Database> greyscale_db;
+    std::shared_ptr<Database> greyscaleColor_db;
     std::shared_ptr<Database> analysis_db;
     std::shared_ptr<Database> vis_db;
 
     signed char *id;    
 	int *NeighborList;
-	int *dvcMap;
-	double *fq;
+	double *fq,*Aq,*Bq;
+    double *Den;
 	double *Permeability;//grey voxel permeability
 	double *Porosity;
 	double *Velocity;
 	double *Pressure_dvc;
+    double *SolidForce;
+    double *DenGradA;
+    double *DenGradB;
     IntArray Map;
     DoubleArray SignDist;
     DoubleArray Velocity_x;
@@ -86,7 +89,9 @@ private:
     char LocalRankFilename[40];
     char LocalRestartFile[40];
    
-    void AssignComponentLabels(double *Porosity, double *Permeablity);
-    
+    void AssignComponentLabels(double *Porosity, double *Permeablity, double *SolidPotential);
+    void AssignSolidForce(double *SolidPotential, double *SolidForce);
+    void DensityField_Init();
+
 };
 
diff --git a/models/GreyscaleModel.cpp b/models/GreyscaleModel.cpp
index 11d92c80..1ef65204 100644
--- a/models/GreyscaleModel.cpp
+++ b/models/GreyscaleModel.cpp
@@ -324,7 +324,6 @@ void ScaLBL_GreyscaleModel::Create(){
 	neighborSize=18*(Np*sizeof(int));
 	//...........................................................................
 	ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
-	ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np);
 	ScaLBL_AllocateDeviceMemory((void **) &fq, 19*dist_mem_size);
 	ScaLBL_AllocateDeviceMemory((void **) &Permeability, sizeof(double)*Np);		
 	ScaLBL_AllocateDeviceMemory((void **) &Porosity, sizeof(double)*Np);		
@@ -332,38 +331,8 @@ void ScaLBL_GreyscaleModel::Create(){
 	ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
 	//...........................................................................
 	// Update GPU data structures
-	if (rank==0)	printf ("Setting up device map and neighbor list \n");
+	if (rank==0)	printf ("Setting up device neighbor list \n");
 	fflush(stdout);
-	int *TmpMap;
-	TmpMap=new int[Np];
-	for (int k=1; k<Nz-1; k++){
-		for (int j=1; j<Ny-1; j++){
-			for (int i=1; i<Nx-1; i++){
-				int idx=Map(i,j,k);
-				if (!(idx < 0))
-					TmpMap[idx] = k*Nx*Ny+j*Nx+i;
-			}
-		}
-	}
-	// check that TmpMap is valid
-	for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
-		int n = TmpMap[idx];
-		if (n > Nx*Ny*Nz){
-			printf("Bad value! idx=%i \n");
-			TmpMap[idx] = Nx*Ny*Nz-1;
-		}
-	}
-	for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
-		int n = TmpMap[idx];
-		if (n > Nx*Ny*Nz){
-			printf("Bad value! idx=%i \n");
-			TmpMap[idx] = Nx*Ny*Nz-1;
-		}
-	}
-	ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
-	ScaLBL_DeviceBarrier();
-	delete [] TmpMap;
-	
 	// copy the neighbor list 
 	ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
 	// initialize phi based on PhaseLabel (include solid component labels)
diff --git a/models/GreyscaleModel.h b/models/GreyscaleModel.h
index a99925b1..3e883b16 100644
--- a/models/GreyscaleModel.h
+++ b/models/GreyscaleModel.h
@@ -62,7 +62,6 @@ public:
 
     signed char *id;    
 	int *NeighborList;
-	int *dvcMap;
 	double *fq;
 	double *Permeability;//grey voxel permeability
 	double *Porosity;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 8b14a9dc..56eda802 100755
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -4,6 +4,7 @@
 ADD_LBPM_EXECUTABLE( lbpm_color_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_greyscale_simulator )
+ADD_LBPM_EXECUTABLE( lbpm_greyscaleColor_simulator )
 #ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator )
 #ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_dfh_simulator )
diff --git a/tests/lbpm_greyscaleColor_simulator.cpp b/tests/lbpm_greyscaleColor_simulator.cpp
new file mode 100644
index 00000000..b836e03c
--- /dev/null
+++ b/tests/lbpm_greyscaleColor_simulator.cpp
@@ -0,0 +1,59 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/stat.h>
+#include <iostream>
+#include <exception>
+#include <stdexcept>
+#include <fstream>
+
+#include "common/ScaLBL.h"
+#include "common/Communication.h"
+#include "common/MPI.h"
+#include "models/GreyscaleColorModel.h"
+//#define WRITE_SURFACES
+
+using namespace std;
+
+
+int main(int argc, char **argv)
+{
+	//*****************************************
+	// ***** MPI STUFF ****************
+	//*****************************************
+	// Initialize MPI
+	int rank,nprocs;
+	MPI_Init(&argc,&argv);
+	MPI_Comm comm = MPI_COMM_WORLD;
+	MPI_Comm_rank(comm,&rank);
+	MPI_Comm_size(comm,&nprocs);
+	{
+		// parallel domain size (# of sub-domains)
+		int nprocx,nprocy,nprocz;
+		int iproc,jproc,kproc;
+
+		if (rank == 0){
+			printf("****************************************\n");
+			printf("Running Greyscale Two-Phase Calculation \n");
+			printf("****************************************\n");
+		}
+		// Initialize compute device
+		int device=ScaLBL_SetDevice(rank);
+		ScaLBL_DeviceBarrier();
+		MPI_Barrier(comm);
+		
+		ScaLBL_GreyscaleColorModel GreyscaleColor(rank,nprocs,comm);
+		auto filename = argv[1];
+		GreyscaleColor.ReadParams(filename);
+		GreyscaleColor.SetDomain();    // this reads in the domain 
+		GreyscaleColor.ReadInput();
+		GreyscaleColor.Create();       // creating the model will create data structure to match the pore structure and allocate variables
+		GreyscaleColor.Initialize();   // initializing the model will set initial conditions for variables
+		GreyscaleColor.Run();	 
+		//GreyscaleColor.VelocityField();
+		GreyscaleColor.WriteDebug();
+	}
+	// ****************************************************
+	MPI_Barrier(comm);
+	MPI_Finalize();
+	// ****************************************************
+}