From 81d726030b35cff1b726f20a166ffff45c7979f2 Mon Sep 17 00:00:00 2001
From: James McClure <jemcclur@gmail.com>
Date: Fri, 12 Feb 2021 09:36:50 -0500
Subject: [PATCH] add freelee gpu skeleton

---
 cuda/FreeLee.cu | 2013 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 2013 insertions(+)
 create mode 100644 cuda/FreeLee.cu

diff --git a/cuda/FreeLee.cu b/cuda/FreeLee.cu
new file mode 100644
index 00000000..bc641ed9
--- /dev/null
+++ b/cuda/FreeLee.cu
@@ -0,0 +1,2013 @@
+#include <math.h>
+
+#define STOKES
+
+__global__ void dvc_ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np)
+{
+	int n;
+    double p = 1.0;//NOTE: take initial pressure p=1.0
+    double chem;
+    double cg_x,cg_y,cg_z;
+
+	//for (n=0; n<Np; 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 ){
+			chem = mu_phi[n];//chemical potential
+			cg_x = ColorGrad[0*Np+n];
+			cg_y = ColorGrad[1*Np+n];
+			cg_z = ColorGrad[2*Np+n];
+
+			gqbar[0*Np+n]  = 0.3333333333333333*p;
+			gqbar[1*Np+n]  = 0.055555555555555555*(p - 0.5*(chem*cg_x+Fx));		//double(100*n)+1.f;
+			gqbar[2*Np+n]  = 0.055555555555555555*(p - 0.5*(-chem*cg_x-Fx));	//double(100*n)+2.f;
+			gqbar[3*Np+n]  = 0.055555555555555555*(p - 0.5*(chem*cg_y+Fy));	//double(100*n)+3.f;
+			gqbar[4*Np+n]  = 0.055555555555555555*(p - 0.5*(-chem*cg_y-Fy));	//double(100*n)+4.f;
+			gqbar[5*Np+n]  = 0.055555555555555555*(p - 0.5*(chem*cg_z+Fz));	//double(100*n)+5.f;
+			gqbar[6*Np+n]  = 0.055555555555555555*(p - 0.5*(-chem*cg_z-Fz));	//double(100*n)+6.f;
+
+			gqbar[7*Np+n]  = 0.0277777777777778*(p-0.5*(chem*(cg_x+cg_y)+Fx+Fy));   //double(100*n)+7.f;
+			gqbar[8*Np+n]  = 0.0277777777777778*(p-0.5*(chem*(-cg_x-cg_y)-Fx-Fy));   //double(100*n)+8.f;
+			gqbar[9*Np+n]  = 0.0277777777777778*(p-0.5*(chem*(cg_x-cg_y)+Fx-Fy));   //double(100*n)+9.f;
+			gqbar[10*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_x+cg_y)-Fx+Fy));  //double(100*n)+10.f;
+
+			gqbar[11*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_x+cg_z)+Fx+Fz));  //double(100*n)+11.f;
+			gqbar[12*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_x-cg_z)-Fx-Fz));  //double(100*n)+12.f;
+			gqbar[13*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_x-cg_z)+Fx-Fz));  //double(100*n)+13.f;
+			gqbar[14*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_x+cg_z)-Fx+Fz));  //double(100*n)+14.f;
+
+			gqbar[15*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_y+cg_z)+Fy+Fz)); ;  //double(100*n)+15.f;
+			gqbar[16*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_y-cg_z)-Fy-Fz));;  //double(100*n)+16.f;
+			gqbar[17*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_y-cg_z)+Fy-Fz)); ;  //double(100*n)+17.f;
+			gqbar[18*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_y+cg_z)-Fy+Fz));;  //double(100*n)+18.f;
+		}
+	}
+}
+
+__global__ void dvc_ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np)
+{
+	int n;
+    double p = 1.0;//NOTE: take initial pressure p=1.0
+
+    //	for (n=0; n<Np; 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 ){        
+    		gqbar[0*Np+n]  = 0.3333333333333333*p;
+    		gqbar[1*Np+n]  = 0.055555555555555555*(p - 0.5*(Fx));		//double(100*n)+1.f;
+    		gqbar[2*Np+n]  = 0.055555555555555555*(p - 0.5*(-Fx));	//double(100*n)+2.f;
+    		gqbar[3*Np+n]  = 0.055555555555555555*(p - 0.5*(Fy));	//double(100*n)+3.f;
+    		gqbar[4*Np+n]  = 0.055555555555555555*(p - 0.5*(-Fy));	//double(100*n)+4.f;
+    		gqbar[5*Np+n]  = 0.055555555555555555*(p - 0.5*(Fz));	//double(100*n)+5.f;
+    		gqbar[6*Np+n]  = 0.055555555555555555*(p - 0.5*(-Fz));	//double(100*n)+6.f;
+
+    		gqbar[7*Np+n]  = 0.0277777777777778*(p-0.5*(Fx+Fy));   //double(100*n)+7.f;
+    		gqbar[8*Np+n]  = 0.0277777777777778*(p-0.5*(-Fx-Fy));   //double(100*n)+8.f;
+    		gqbar[9*Np+n]  = 0.0277777777777778*(p-0.5*(Fx-Fy));   //double(100*n)+9.f;
+    		gqbar[10*Np+n] = 0.0277777777777778*(p-0.5*(-Fx+Fy));  //double(100*n)+10.f;
+
+    		gqbar[11*Np+n] = 0.0277777777777778*(p-0.5*(Fx+Fz));  //double(100*n)+11.f;
+    		gqbar[12*Np+n] = 0.0277777777777778*(p-0.5*(-Fx-Fz));  //double(100*n)+12.f;
+    		gqbar[13*Np+n] = 0.0277777777777778*(p-0.5*(Fx-Fz));  //double(100*n)+13.f;
+    		gqbar[14*Np+n] = 0.0277777777777778*(p-0.5*(-Fx+Fz));  //double(100*n)+14.f;
+
+    		gqbar[15*Np+n] = 0.0277777777777778*(p-0.5*(Fy+Fz)); ;  //double(100*n)+15.f;
+    		gqbar[16*Np+n] = 0.0277777777777778*(p-0.5*(-Fy-Fz));;  //double(100*n)+16.f;
+    		gqbar[17*Np+n] = 0.0277777777777778*(p-0.5*(Fy-Fz)); ;  //double(100*n)+17.f;
+    		gqbar[18*Np+n] = 0.0277777777777778*(p-0.5*(-Fy+Fz));;  //double(100*n)+18.f;
+    	}
+    }
+}
+
+__global__ void dvc_ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad, 
+                                                    double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
+	int idx,n;
+	double phi;
+    double nx,ny,nz,cg_mag;
+    double theta;//anti-diffusion term
+    double cs2_inv = 4.5;//inverse of speed of sound for D3Q7
+    double M = 1.0/cs2_inv*(tauM-0.5);//diffusivity (or mobility)
+
+	for (idx=start; idx<finish; idx++){
+
+		n = Map[idx];
+		phi = Phi[n];
+		if (phi > 1.f)   phi = 1.0;
+		if (phi < -1.f)  phi = -1.0;
+		Den[idx] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
+
+        //compute unit normal of color gradient
+        nx = ColorGrad[idx+0*Np];
+        ny = ColorGrad[idx+1*Np];
+        nz = ColorGrad[idx+2*Np];
+        cg_mag = sqrt(nx*nx+ny*ny+nz*nz);
+		double ColorMag_temp = cg_mag;
+		if (cg_mag==0.0) ColorMag_temp=1.0;
+		nx = nx/ColorMag_temp;
+		ny = ny/ColorMag_temp;
+		nz = nz/ColorMag_temp;		
+
+        theta = M*cs2_inv*(1-4.0*phi*phi)/W;
+
+		hq[0*Np+idx]=0.3333333333333333*(phi);
+		hq[1*Np+idx]=0.1111111111111111*(phi+theta*nx);
+		hq[2*Np+idx]=0.1111111111111111*(phi-theta*nx);
+		hq[3*Np+idx]=0.1111111111111111*(phi+theta*ny);
+		hq[4*Np+idx]=0.1111111111111111*(phi-theta*ny);
+		hq[5*Np+idx]=0.1111111111111111*(phi+theta*nz);
+		hq[6*Np+idx]=0.1111111111111111*(phi-theta*nz);
+		
+	}
+}
+
+__global__ void dvc_ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, 
+                                                          double rhoA, double rhoB, int start, int finish, int Np){
+
+	int idx,n,nread;
+	double fq,phi;
+	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 ){
+
+			// q=0
+			fq = hq[n];
+			phi = fq;
+
+			// q=1
+			nread = neighborList[n]; 
+			fq = hq[nread];
+			phi += fq;
+
+			// q=2
+			nread = neighborList[n+Np]; 
+			fq = hq[nread];  
+			phi += fq;
+
+			// q=3
+			nread = neighborList[n+2*Np]; 
+			fq = hq[nread];
+			phi += fq;
+
+			// q = 4
+			nread = neighborList[n+3*Np]; 
+			fq = hq[nread];
+			phi += fq;
+
+			// q=5
+			nread = neighborList[n+4*Np];
+			fq = hq[nread];
+			phi += fq;
+
+			// q = 6
+			nread = neighborList[n+5*Np];
+			fq = hq[nread];
+			phi += fq;
+
+			// save the number densities
+			Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
+
+			// save the phase indicator field
+			idx = Map[n];
+			Phi[idx] = phi; 
+		}
+	}
+}
+
+__global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, 
+			                                               double rhoA, double rhoB, int start, int finish, int Np){
+	int idx,n;
+	double fq,phi;
+	//	for (int n=start; n<finish; 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 + start;
+
+		if ( n<finish ){		
+			// q=0
+			fq = hq[n];
+			phi = fq;
+
+			// q=1
+			fq = hq[2*Np+n];
+			phi += fq;
+
+			// f2 = hq[10*Np+n];
+			fq = hq[1*Np+n];
+			phi += fq;
+
+			// q=3
+			fq = hq[4*Np+n];
+			phi += fq;
+
+			// q = 4
+			fq = hq[3*Np+n];
+			phi += fq;
+
+			// q=5
+			fq = hq[6*Np+n];
+			phi += fq;
+
+			// q = 6
+			fq = hq[5*Np+n];
+			phi += fq;
+
+			// save the number densities
+			Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
+
+			// save the phase indicator field
+			idx = Map[n];
+			Phi[idx] = phi; 	
+		}
+	}	
+}
+
+__global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, double *dist, double *hq, double *Den,	double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, 
+                                                double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 
+                                                int strideY, int strideZ, int start, int finish, int Np){
+	
+	int n,nn,nn2x,ijk;
+	int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
+    double ux,uy,uz;//fluid velocity 
+    double p;//pressure
+    double chem;//chemical potential
+    double phi; //phase field
+    double rho0;//fluid density
+	// distribution functions
+	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
+	double m0,m3,m5,m7;
+	double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
+	double mm3,mm5,mm7;
+    double feq0,feq1,feq2,feq3,feq4,feq5,feq6,feq7,feq8,feq9,feq10,feq11,feq12,feq13,feq14,feq15,feq16,feq17,feq18;
+    double nx,ny,nz;//normal color gradient
+    double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
+
+	//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
+    double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
+	double tau;//position dependent LB relaxation time for fluid
+    double C,theta;
+    double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
+
+    //	for (int n=start; n<finish; 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 + start;
+
+    	if ( n<finish ){
+    		rho0 = Den[n];//load density
+    		phi = Phi[n];// load phase field
+    		// local relaxation time
+    		tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
+
+    		// Get the 1D index based on regular data layout
+    		ijk = Map[n];
+    		//					COMPUTE THE COLOR GRADIENT
+    		//........................................................................
+    		//.................Read Phase Indicator Values............................
+    		//........................................................................
+    		nn = ijk-1;							// neighbor index (get convention)
+    		m1 = Phi[nn];						// get neighbor for phi - 1
+    		//........................................................................
+    		nn = ijk+1;							// neighbor index (get convention)
+    		m2 = Phi[nn];						// get neighbor for phi - 2
+    		//........................................................................
+    		nn = ijk-strideY;							// neighbor index (get convention)
+    		m3 = Phi[nn];					// get neighbor for phi - 3
+    		//........................................................................
+    		nn = ijk+strideY;							// neighbor index (get convention)
+    		m4 = Phi[nn];					// get neighbor for phi - 4
+    		//........................................................................
+    		nn = ijk-strideZ;						// neighbor index (get convention)
+    		m5 = Phi[nn];					// get neighbor for phi - 5
+    		//........................................................................
+    		nn = ijk+strideZ;						// neighbor index (get convention)
+    		m6 = Phi[nn];					// get neighbor for phi - 6
+    		//........................................................................
+    		nn = ijk-strideY-1;						// neighbor index (get convention)
+    		m7 = Phi[nn];					// get neighbor for phi - 7
+    		//........................................................................
+    		nn = ijk+strideY+1;						// neighbor index (get convention)
+    		m8 = Phi[nn];					// get neighbor for phi - 8
+    		//........................................................................
+    		nn = ijk+strideY-1;						// neighbor index (get convention)
+    		m9 = Phi[nn];					// get neighbor for phi - 9
+    		//........................................................................
+    		nn = ijk-strideY+1;						// neighbor index (get convention)
+    		m10 = Phi[nn];					// get neighbor for phi - 10
+    		//........................................................................
+    		nn = ijk-strideZ-1;						// neighbor index (get convention)
+    		m11 = Phi[nn];					// get neighbor for phi - 11
+    		//........................................................................
+    		nn = ijk+strideZ+1;						// neighbor index (get convention)
+    		m12 = Phi[nn];					// get neighbor for phi - 12
+    		//........................................................................
+    		nn = ijk+strideZ-1;						// neighbor index (get convention)
+    		m13 = Phi[nn];					// get neighbor for phi - 13
+    		//........................................................................
+    		nn = ijk-strideZ+1;						// neighbor index (get convention)
+    		m14 = Phi[nn];					// get neighbor for phi - 14
+    		//........................................................................
+    		nn = ijk-strideZ-strideY;					// neighbor index (get convention)
+    		m15 = Phi[nn];					// get neighbor for phi - 15
+    		//........................................................................
+    		nn = ijk+strideZ+strideY;					// neighbor index (get convention)
+    		m16 = Phi[nn];					// get neighbor for phi - 16
+    		//........................................................................
+    		nn = ijk+strideZ-strideY;					// neighbor index (get convention)
+    		m17 = Phi[nn];					// get neighbor for phi - 17
+    		//........................................................................
+    		nn = ijk-strideZ+strideY;					// neighbor index (get convention)
+    		m18 = Phi[nn];					// get neighbor for phi - 18
+
+    		// compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
+    		//........................................................................
+    		nn2x = ijk-2;							// neighbor index (get convention)
+    		mm1 = Phi[nn2x];						// get neighbor for phi - 1
+    		mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
+    		//........................................................................
+    		nn2x = ijk+2;							// neighbor index (get convention)
+    		mm2 = Phi[nn2x];						// get neighbor for phi - 2
+    		mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
+    		//........................................................................
+    		nn2x = ijk-strideY*2;							// neighbor index (get convention)
+    		mm3 = Phi[nn2x];					// get neighbor for phi - 3
+    		mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
+    		//........................................................................
+    		nn2x = ijk+strideY*2;							// neighbor index (get convention)
+    		mm4 = Phi[nn2x];					// get neighbor for phi - 4
+    		mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
+    		//........................................................................
+    		nn2x = ijk-strideZ*2;						// neighbor index (get convention)
+    		mm5 = Phi[nn2x];					// get neighbor for phi - 5
+    		mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
+    		//........................................................................
+    		nn2x = ijk+strideZ*2;						// neighbor index (get convention)
+    		mm6 = Phi[nn2x];					// get neighbor for phi - 6
+    		mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
+    		//........................................................................
+    		nn2x = ijk-strideY*2-2;						// neighbor index (get convention)
+    		mm7 = Phi[nn2x];					// get neighbor for phi - 7
+    		mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
+    		//........................................................................
+    		nn2x = ijk+strideY*2+2;						// neighbor index (get convention)
+    		mm8 = Phi[nn2x];					// get neighbor for phi - 8
+    		mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
+    		//........................................................................
+    		nn2x = ijk+strideY*2-2;						// neighbor index (get convention)
+    		mm9 = Phi[nn2x];					// get neighbor for phi - 9
+    		mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
+    		//........................................................................
+    		nn2x = ijk-strideY*2+2;						// neighbor index (get convention)
+    		mm10 = Phi[nn2x];					// get neighbor for phi - 10
+    		mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
+    		//........................................................................
+    		nn2x = ijk-strideZ*2-2;						// neighbor index (get convention)
+    		mm11 = Phi[nn2x];					// get neighbor for phi - 11
+    		mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
+    		//........................................................................
+    		nn2x = ijk+strideZ*2+2;						// neighbor index (get convention)
+    		mm12 = Phi[nn2x];					// get neighbor for phi - 12
+    		mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
+    		//........................................................................
+    		nn2x = ijk+strideZ*2-2;						// neighbor index (get convention)
+    		mm13 = Phi[nn2x];					// get neighbor for phi - 13
+    		mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
+    		//........................................................................
+    		nn2x = ijk-strideZ*2+2;						// neighbor index (get convention)
+    		mm14 = Phi[nn2x];					// get neighbor for phi - 14
+    		mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
+    		//........................................................................
+    		nn2x = ijk-strideZ*2-strideY*2;					// neighbor index (get convention)
+    		mm15 = Phi[nn2x];					// get neighbor for phi - 15
+    		mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
+    		//........................................................................
+    		nn2x = ijk+strideZ*2+strideY*2;					// neighbor index (get convention)
+    		mm16 = Phi[nn2x];					// get neighbor for phi - 16
+    		mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
+    		//........................................................................
+    		nn2x = ijk+strideZ*2-strideY*2;					// neighbor index (get convention)
+    		mm17 = Phi[nn2x];					// get neighbor for phi - 17
+    		mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
+    		//........................................................................
+    		nn2x = ijk-strideZ*2+strideY*2;					// neighbor index (get convention)
+    		mm18 = Phi[nn2x];					// get neighbor for phi - 18
+    		mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);
+
+
+    		//............Compute the Color Gradient...................................
+    		nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+    		ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+    		nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+    		//............Compute the Chemical Potential...............................
+    		chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
+    		chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
+    		//............Compute the Mixed Gradient...................................
+    		mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14))*0.25;//the factor of 0.25 comes from the denominator of Eq.30
+    		mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18))*0.25;
+    		mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18))*0.25;
+
+    		// q=0
+    		m0 = dist[n];
+    		// q=1
+    		nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
+    		m1 = dist[nr1]; // reading the f1 data into register fq
+
+    		nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
+    		m2 = dist[nr2];  // reading the f2 data into register fq
+
+    		// q=3
+    		nr3 = neighborList[n+2*Np]; // neighbor 4
+    		m3 = dist[nr3];
+
+    		// q = 4
+    		nr4 = neighborList[n+3*Np]; // neighbor 3
+    		m4 = dist[nr4];
+
+    		// q=5
+    		nr5 = neighborList[n+4*Np];
+    		m5 = dist[nr5];
+
+    		// q = 6
+    		nr6 = neighborList[n+5*Np];
+    		m6 = dist[nr6];
+
+    		// q=7
+    		nr7 = neighborList[n+6*Np];
+    		m7 = dist[nr7];
+
+    		// q = 8
+    		nr8 = neighborList[n+7*Np];
+    		m8 = dist[nr8];
+
+    		// q=9
+    		nr9 = neighborList[n+8*Np];
+    		m9 = dist[nr9];
+
+    		// q = 10
+    		nr10 = neighborList[n+9*Np];
+    		m10 = dist[nr10];
+
+    		// q=11
+    		nr11 = neighborList[n+10*Np];
+    		m11 = dist[nr11];
+
+    		// q=12
+    		nr12 = neighborList[n+11*Np];
+    		m12 = dist[nr12];
+
+    		// q=13
+    		nr13 = neighborList[n+12*Np];
+    		m13 = dist[nr13];
+
+    		// q=14
+    		nr14 = neighborList[n+13*Np];
+    		m14 = dist[nr14];
+
+    		// q=15
+    		nr15 = neighborList[n+14*Np];
+    		m15 = dist[nr15];
+
+    		// q=16
+    		nr16 = neighborList[n+15*Np];
+    		m16 = dist[nr16];
+
+    		// q=17
+    		nr17 = neighborList[n+16*Np];
+    		m17 = dist[nr17];
+
+    		// q=18
+    		nr18 = neighborList[n+17*Np];
+    		m18 = dist[nr18];
+
+    		//compute fluid velocity
+    		ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0);
+    		uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0);
+    		uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
+    		//compute pressure
+    		p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17)
+                		  +0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz);
+
+    		//compute equilibrium distributions
+    		feq0 = 0.3333333333333333*p - 0.25*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) - 
+    				0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz) - 0.5*(-(nx*ux) - ny*uy - nz*uz)*
+    				(-0.08333333333333333*(rhoA - rhoB)*(ux*ux + uy*uy + uz*uz) + chem*(0.3333333333333333 - 0.5*(ux*ux + uy*uy + uz*uz)));
+    		feq1 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) - 
+    				0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
+    						0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx - nx*ux - ny*uy - nz*uz)*
+    						(2*chem*ux*ux - 0.3333333333333333*((-rhoA + rhoB)*ux*ux + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz)));
+    		feq2 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) - 
+    				0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
+    						0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx + nx*ux + ny*uy + nz*uz)*
+    						(-2.*chem*ux*ux + 0.1111111111111111*(-4.*chem + rhoB*(-2.*ux - 1.*ux*ux - 1.*uy*uy - 1.*uz*uz) + 
+    								rhoA*(2.*ux + ux*ux + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*ux*ux + 
+    										chem*(4.*ux + 2.*ux*ux + 2.*uy*uy + 2.*uz*uz)));
+    		feq3 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) - 
+    				0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
+    						0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny - nx*ux - ny*uy - nz*uz)*
+    						(2*chem*uy*uy - 0.3333333333333333*((-rhoA + rhoB)*uy*uy + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz)));
+    		feq4 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) - 
+    				0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
+    						0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny + nx*ux + ny*uy + nz*uz)*
+    						(-2.*chem*uy*uy + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 2.*uy - 1.*uy*uy - 1.*uz*uz) + 
+    								rhoA*(ux*ux + 2.*uy + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uy*uy + 
+    										chem*(2.*ux*ux + 4.*uy + 2.*uy*uy + 2.*uz*uz))); 
+    		feq5 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) - 
+    				0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)) - 0.0625*(nx*ux + ny*uy + nz*(-1. + uz))*
+    						(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (2. - 1.*uz)*uz) + 
+    								rhoA*(ux*ux + uy*uy + (-2. + uz)*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
+    										chem*(2.*ux*ux + 2.*uy*uy + uz*(-4. + 2.*uz)))); 
+    		feq6 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) - 
+    				0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))) - 0.0625*(nz + nx*ux + ny*uy + nz*uz)*
+    						(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (-2. - 1.*uz)*uz) + 
+    								rhoA*(ux*ux + uy*uy + uz*(2. + uz))) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
+    										chem*(2.*ux*ux + 2.*uy*uy + uz*(4. + 2.*uz)))); 
+    		feq7 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux + uy)*(ux + uy) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
+    				0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
+    						0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx + ny - nx*ux - ny*uy - nz*uz)*
+    						(2*chem*(ux + uy)*(ux + uy) + 0.3333333333333333*((rhoA - rhoB)*(ux + uy)*(ux + uy) - 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
+    		feq8 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
+    				0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
+    						0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(-(nx*(1 + ux)) - ny*(1 + uy) - nz*uz)*
+    						(2*chem*(ux + uy)*(ux + uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uy)*(ux + uy)) + 
+    								2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
+    		feq9 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
+    				0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
+    						0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx - nx*ux - ny*(1 + uy) - nz*uz)*
+    						(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
+    								2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
+    		feq10 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
+    				0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
+    						0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(ny - nx*(1 + ux) - ny*uy - nz*uz)*
+    						(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
+    								2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
+    		feq11 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux + uz)*(ux + uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
+    				0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
+    						0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nx + nz - nx*ux - ny*uy - nz*uz)*
+    						(2*chem*(ux + uz)*(ux + uz) + 0.3333333333333333*((rhoA - rhoB)*(ux + uz)*(ux + uz) - 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    								0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
+    		feq12 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
+    				0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
+    						0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*(1 + ux)) - ny*uy - nz*(1 + uz))*
+    						(2*chem*(ux + uz)*(ux + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uz)*(ux + uz)) + 
+    								2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
+    		feq13 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
+    				0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
+    						0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(nx - nx*ux - ny*uy - nz*(1 + uz))*
+    						(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
+    								2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
+    		feq14 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
+    				0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
+    						0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*(1 + ux) - ny*uy - nz*uz)*
+    						(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
+    								2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
+    		feq15 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) - 
+    				0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(ny + nz - nx*ux - ny*uy - nz*uz)*
+    						(2*chem*(uy + uz)*(uy + uz) + 0.3333333333333333*((rhoA - rhoB)*(uy + uz)*(uy + uz) - 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
+    								0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))); 
+    		feq16 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) - 
+    				0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*ux) - ny*(1 + uy) - nz*(1 + uz))*
+    						(2*chem*(uy + uz)*(uy + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy + uz)*(uy + uz)) + 
+    								2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))); 
+    		feq17 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) - 
+    				0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(ny - nx*ux - ny*uy - nz*(1 + uz))*
+    						(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
+    								2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))); 
+    		feq18 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) - 
+    				0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*ux - ny*(1 + uy) - nz*uz)*
+    						(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
+    								2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))); 
+
+    		//------------------------------------------------- BCK collison ------------------------------------------------------------//
+    		// q=0
+    		dist[n] = m0 - (m0-feq0)/tau + 0.25*(2*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 
+    				(mgx*ux + mgy*uy + mgz*uz)*(2*chem*(ux*ux + uy*uy + uz*uz) + 
+    						0.3333333333333333*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*uz)))); 
+
+    		// q = 1
+    		dist[nr2] = m1 - (m1-feq1)/tau + 0.125*(2*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
+    				0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
+    				(mgx*(-1 + ux) + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz))));
+
+    		// q=2
+    		dist[nr1] = m2 - (m2-feq2)/tau + 0.125*(2*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
+    				0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) + 
+    				(mgx + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(2*ux + ux*ux + uy*uy + uz*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*uz))));
+
+    		// q = 3
+    		dist[nr4] = m3 - (m3-feq3)/tau + 0.125*(2*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy + 
+    				0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    				(mgx*ux + mgy*(-1 + uy) + mgz*uz)*(-2*chem*(uy*uy) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz))));
+
+    		// q = 4
+    		dist[nr3] = m4 - (m4-feq4)/tau + 0.125*(2*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy + 
+    				0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    				(mgy + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uy*uy) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*uz))));
+
+    		// q = 5
+    		dist[nr6] = m5 - (m5-feq5)/tau + 0.125*(2*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz + 
+    				0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    				(mgx*ux + mgy*uy + mgz*(-1 + uz))*(-2*chem*(uz*uz) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + (-2 + uz)*uz))));
+
+    		// q = 6
+    		dist[nr5] = m6 - (m6-feq6)/tau + 0.125*(2*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz + 
+    				0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) + 
+    				(mgz + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uz*uz) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + uz*(2 + uz))) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*(2 + uz)))));
+
+    		// q = 7
+    		dist[nr8] = m7 - (m7-feq7)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
+    				(0.2222222222222222 + (ux + uy)*(ux + uy) - 
+    						0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    						(mgx*(-1 + ux) + mgy*(-1 + uy) + mgz*uz)*
+    						(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
+
+    		// q = 8
+    		dist[nr7] = m8 - (m8-feq8)/tau + 0.0625*(2*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
+    				(-0.2222222222222222 - (ux + uy)*(ux + uy) + 
+    						0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    						(mgx + mgy + mgx*ux + mgy*uy + mgz*uz)*
+    						(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
+
+    		// q = 9
+    		dist[nr10] = m9 - (m9-feq9)/tau + 0.0625*(2*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
+    				(-0.2222222222222222 - (ux - uy)*(ux - uy) + 
+    						0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    						(mgy + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
+    						(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
+
+    		// q = 10
+    		dist[nr9] = m10 - (m10-feq10)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
+    				(-0.2222222222222222 - (ux - uy)*(ux - uy) + 
+    						0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    						(mgx*(1 + ux) + mgy*(-1 + uy) + mgz*uz)*
+    						(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
+
+    		// q = 11
+    		dist[nr12] = m11 - (m11-feq11)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
+    				(0.2222222222222222 + (ux + uz)*(ux + uz) - 
+    						0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    						(mgx*(-1 + ux) + mgy*uy + mgz*(-1 + uz))*
+    						(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
+
+    		// q = 12
+    		dist[nr11] = m12 - (m12-feq12)/tau + 0.0625*(2*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
+    				(-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))
+    				+ (mgx + mgz + mgx*ux + mgy*uy + mgz*uz)*
+    				(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
+    						(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
+
+    		// q = 13
+    		dist[nr14] = m13 - (m13-feq13)/tau + 0.0625*(2*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
+    				(-0.2222222222222222 - (ux - uz)*(ux - uz) + 
+    						0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
+    						(mgz + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
+    						(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
+
+    		// q= 14
+    		dist[nr13] = m14 - (m14-feq14)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
+    				(-0.2222222222222222 - (ux - uz)*(ux - uz) + 
+    						0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    						(mgx*(1 + ux) + mgy*uy + mgz*(-1 + uz))*
+    						(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
+
+    		// q = 15
+    		dist[nr16] = m15 - (m15-feq15)/tau + 0.0625*(-2*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
+    				(0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))
+    				+ (mgx*ux + mgy*(-1 + uy) + mgz*(-1 + uz))*
+    				(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
+    						(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))));
+
+    		// q = 16
+    		dist[nr15] = m16 - (m16-feq16)/tau + 0.0625*(2*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
+    				(-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))
+    				+ (mgy + mgz + mgx*ux + mgy*uy + mgz*uz)*
+    				(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
+    						(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))));
+
+    		// q = 17
+    		dist[nr18] = m17 - (m17-feq17)/tau + 0.0625*(2*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
+    				(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))
+    				+ (mgz + mgx*ux + mgy*(-1 + uy) + mgz*uz)*
+    				(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
+    						(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))));
+
+    		// q = 18
+    		dist[nr17] = m18 - (m18-feq18)/tau + 0.0625*(2*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
+    				(-0.2222222222222222 - (uy - uz)*(uy - uz) + 
+    						0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
+    						(mgx*ux + mgy*(1 + uy) + mgz*(-1 + uz))*
+    						(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
+    		//----------------------------------------------------------------------------------------------------------------------------------------//
+
+
+    		// ----------------------------- compute phase field evolution ----------------------------------------
+    		//Normalize the Color Gradient
+    		C = sqrt(nx*nx+ny*ny+nz*nz);
+    		double ColorMag = C;
+    		if (C==0.0) ColorMag=1.0;
+    		nx = nx/ColorMag;
+    		ny = ny/ColorMag;
+    		nz = nz/ColorMag;		
+    		//compute surface tension-related parameter
+    		theta = M*4.5*(1-4.0*phi*phi)/W;
+
+    		//load distributions of phase field
+    		//q=0
+    		h0 = hq[n];
+    		//q=1
+    		h1 = hq[nr1]; 
+
+    		//q=2
+    		h2 = hq[nr2];  
+
+    		//q=3
+    		h3 = hq[nr3];
+
+    		//q=4
+    		h4 = hq[nr4];
+
+    		//q=5
+    		h5 = hq[nr5];
+
+    		//q=6
+    		h6 = hq[nr6];
+
+    		//-------------------------------- BGK collison for phase field ---------------------------------//
+    		// q = 0
+    		hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
+
+    		// q = 1
+    		hq[nr2] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
+
+    		// q = 2
+    		hq[nr1] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
+
+    		// q = 3
+    		hq[nr4] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
+
+    		// q = 4
+    		hq[nr3] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
+
+    		// q = 5
+    		hq[nr6] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
+
+    		// q = 6
+    		hq[nr5] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
+    		//........................................................................
+
+    		//Update velocity on device
+    		Vel[0*Np+n] = ux;
+    		Vel[1*Np+n] = uy;
+    		Vel[2*Np+n] = uz;
+    		//Update pressure on device
+    		Pressure[n] = p;
+    		//Update chemical potential on device
+    		mu_phi[n] = chem;
+    		//Update color gradient on device
+    		ColorGrad[0*Np+n] = nx;
+    		ColorGrad[1*Np+n] = ny;
+    		ColorGrad[2*Np+n] = nz;
+    	}
+    }
+}
+
+__global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *hq, double *Den,	double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
+                                                double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 
+                                                int strideY, int strideZ, int start, int finish, int Np){
+	
+	int n,nn,nn2x,ijk;
+	//int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
+    double ux,uy,uz;//fluid velocity 
+    double p;//pressure
+    double chem;//chemical potential
+    double phi; //phase field
+    double rho0;//fluid density
+	// distribution functions
+	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
+	double m0,m3,m5,m7;
+	double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
+	double mm3,mm5,mm7;
+    double feq0,feq1,feq2,feq3,feq4,feq5,feq6,feq7,feq8,feq9,feq10,feq11,feq12,feq13,feq14,feq15,feq16,feq17,feq18;
+    double nx,ny,nz;//normal color gradient
+    double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
+
+	//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
+    double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
+	double tau;//position dependent LB relaxation time for fluid
+    double C,theta;
+    double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
+
+    //	for (int n=start; n<finish; 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 + start;
+
+    	if ( n<finish ){
+    		rho0 = Den[n];//load density
+    		phi = Phi[n];// load phase field
+    		// local relaxation time
+    		tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
+
+    		// Get the 1D index based on regular data layout
+    		ijk = Map[n];
+    		//					COMPUTE THE COLOR GRADIENT
+    		//........................................................................
+    		//.................Read Phase Indicator Values............................
+    		//........................................................................
+    		nn = ijk-1;							// neighbor index (get convention)
+    		m1 = Phi[nn];						// get neighbor for phi - 1
+    		//........................................................................
+    		nn = ijk+1;							// neighbor index (get convention)
+    		m2 = Phi[nn];						// get neighbor for phi - 2
+    		//........................................................................
+    		nn = ijk-strideY;							// neighbor index (get convention)
+    		m3 = Phi[nn];					// get neighbor for phi - 3
+    		//........................................................................
+    		nn = ijk+strideY;							// neighbor index (get convention)
+    		m4 = Phi[nn];					// get neighbor for phi - 4
+    		//........................................................................
+    		nn = ijk-strideZ;						// neighbor index (get convention)
+    		m5 = Phi[nn];					// get neighbor for phi - 5
+    		//........................................................................
+    		nn = ijk+strideZ;						// neighbor index (get convention)
+    		m6 = Phi[nn];					// get neighbor for phi - 6
+    		//........................................................................
+    		nn = ijk-strideY-1;						// neighbor index (get convention)
+    		m7 = Phi[nn];					// get neighbor for phi - 7
+    		//........................................................................
+    		nn = ijk+strideY+1;						// neighbor index (get convention)
+    		m8 = Phi[nn];					// get neighbor for phi - 8
+    		//........................................................................
+    		nn = ijk+strideY-1;						// neighbor index (get convention)
+    		m9 = Phi[nn];					// get neighbor for phi - 9
+    		//........................................................................
+    		nn = ijk-strideY+1;						// neighbor index (get convention)
+    		m10 = Phi[nn];					// get neighbor for phi - 10
+    		//........................................................................
+    		nn = ijk-strideZ-1;						// neighbor index (get convention)
+    		m11 = Phi[nn];					// get neighbor for phi - 11
+    		//........................................................................
+    		nn = ijk+strideZ+1;						// neighbor index (get convention)
+    		m12 = Phi[nn];					// get neighbor for phi - 12
+    		//........................................................................
+    		nn = ijk+strideZ-1;						// neighbor index (get convention)
+    		m13 = Phi[nn];					// get neighbor for phi - 13
+    		//........................................................................
+    		nn = ijk-strideZ+1;						// neighbor index (get convention)
+    		m14 = Phi[nn];					// get neighbor for phi - 14
+    		//........................................................................
+    		nn = ijk-strideZ-strideY;					// neighbor index (get convention)
+    		m15 = Phi[nn];					// get neighbor for phi - 15
+    		//........................................................................
+    		nn = ijk+strideZ+strideY;					// neighbor index (get convention)
+    		m16 = Phi[nn];					// get neighbor for phi - 16
+    		//........................................................................
+    		nn = ijk+strideZ-strideY;					// neighbor index (get convention)
+    		m17 = Phi[nn];					// get neighbor for phi - 17
+    		//........................................................................
+    		nn = ijk-strideZ+strideY;					// neighbor index (get convention)
+    		m18 = Phi[nn];					// get neighbor for phi - 18
+
+    		// compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
+    		//........................................................................
+    		nn2x = ijk-2;							// neighbor index (get convention)
+    		mm1 = Phi[nn2x];						// get neighbor for phi - 1
+    		mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
+    		//........................................................................
+    		nn2x = ijk+2;							// neighbor index (get convention)
+    		mm2 = Phi[nn2x];						// get neighbor for phi - 2
+    		mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
+    		//........................................................................
+    		nn2x = ijk-strideY*2;							// neighbor index (get convention)
+    		mm3 = Phi[nn2x];					// get neighbor for phi - 3
+    		mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
+    		//........................................................................
+    		nn2x = ijk+strideY*2;							// neighbor index (get convention)
+    		mm4 = Phi[nn2x];					// get neighbor for phi - 4
+    		mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
+    		//........................................................................
+    		nn2x = ijk-strideZ*2;						// neighbor index (get convention)
+    		mm5 = Phi[nn2x];					// get neighbor for phi - 5
+    		mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
+    		//........................................................................
+    		nn2x = ijk+strideZ*2;						// neighbor index (get convention)
+    		mm6 = Phi[nn2x];					// get neighbor for phi - 6
+    		mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
+    		//........................................................................
+    		nn2x = ijk-strideY*2-2;						// neighbor index (get convention)
+    		mm7 = Phi[nn2x];					// get neighbor for phi - 7
+    		mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
+    		//........................................................................
+    		nn2x = ijk+strideY*2+2;						// neighbor index (get convention)
+    		mm8 = Phi[nn2x];					// get neighbor for phi - 8
+    		mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
+    		//........................................................................
+    		nn2x = ijk+strideY*2-2;						// neighbor index (get convention)
+    		mm9 = Phi[nn2x];					// get neighbor for phi - 9
+    		mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
+    		//........................................................................
+    		nn2x = ijk-strideY*2+2;						// neighbor index (get convention)
+    		mm10 = Phi[nn2x];					// get neighbor for phi - 10
+    		mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
+    		//........................................................................
+    		nn2x = ijk-strideZ*2-2;						// neighbor index (get convention)
+    		mm11 = Phi[nn2x];					// get neighbor for phi - 11
+    		mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
+    		//........................................................................
+    		nn2x = ijk+strideZ*2+2;						// neighbor index (get convention)
+    		mm12 = Phi[nn2x];					// get neighbor for phi - 12
+    		mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
+    		//........................................................................
+    		nn2x = ijk+strideZ*2-2;						// neighbor index (get convention)
+    		mm13 = Phi[nn2x];					// get neighbor for phi - 13
+    		mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
+    		//........................................................................
+    		nn2x = ijk-strideZ*2+2;						// neighbor index (get convention)
+    		mm14 = Phi[nn2x];					// get neighbor for phi - 14
+    		mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
+    		//........................................................................
+    		nn2x = ijk-strideZ*2-strideY*2;					// neighbor index (get convention)
+    		mm15 = Phi[nn2x];					// get neighbor for phi - 15
+    		mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
+    		//........................................................................
+    		nn2x = ijk+strideZ*2+strideY*2;					// neighbor index (get convention)
+    		mm16 = Phi[nn2x];					// get neighbor for phi - 16
+    		mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
+    		//........................................................................
+    		nn2x = ijk+strideZ*2-strideY*2;					// neighbor index (get convention)
+    		mm17 = Phi[nn2x];					// get neighbor for phi - 17
+    		mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
+    		//........................................................................
+    		nn2x = ijk-strideZ*2+strideY*2;					// neighbor index (get convention)
+    		mm18 = Phi[nn2x];					// get neighbor for phi - 18
+    		mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);
+
+
+    		//............Compute the Color Gradient...................................
+    		nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+    		ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+    		nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+    		//............Compute the Chemical Potential...............................
+    		chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
+    		chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
+    		//............Compute the Mixed Gradient...................................
+    		mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14))*0.25;//the factor of 0.25 comes from the denominator of Eq.30
+    		mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18))*0.25;
+    		mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18))*0.25;
+
+    		// q=0
+    		m0 = dist[n];
+    		// q=1
+    		m1 = dist[2*Np+n]; 
+
+    		// q=2
+    		m2 = dist[1*Np+n];  
+
+    		// q=3
+    		m3 = dist[4*Np+n];
+
+    		// q = 4
+    		m4 = dist[3*Np+n];
+
+    		// q=5
+    		m5 = dist[6*Np+n];
+
+    		// q = 6
+    		m6 = dist[5*Np+n];
+
+    		// q=7
+    		m7 = dist[8*Np+n];
+
+    		// q = 8
+    		m8 = dist[7*Np+n];
+
+    		// q=9
+    		m9 = dist[10*Np+n];
+
+    		// q = 10
+    		m10 = dist[9*Np+n];
+
+    		// q=11
+    		m11 = dist[12*Np+n];
+
+    		// q=12
+    		m12 = dist[11*Np+n];
+
+    		// q=13
+    		m13 = dist[14*Np+n];
+
+    		// q=14
+    		m14 = dist[13*Np+n];
+
+    		// q=15
+    		m15 = dist[16*Np+n];
+
+    		// q=16
+    		m16 = dist[15*Np+n];
+
+    		// q=17
+    		m17 = dist[18*Np+n];
+
+    		// q=18
+    		m18 = dist[17*Np+n];
+
+    		//compute fluid velocity
+    		ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0);
+    		uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0);
+    		uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
+    		//compute pressure
+    		p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17)
+                		  +0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz);
+
+    		//compute equilibrium distributions
+    		feq0 = 0.3333333333333333*p - 0.25*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) - 
+    				0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz) - 0.5*(-(nx*ux) - ny*uy - nz*uz)*
+    				(-0.08333333333333333*(rhoA - rhoB)*(ux*ux + uy*uy + uz*uz) + chem*(0.3333333333333333 - 0.5*(ux*ux + uy*uy + uz*uz)));
+    		feq1 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) - 
+    				0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
+    						0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx - nx*ux - ny*uy - nz*uz)*
+    						(2*chem*ux*ux - 0.3333333333333333*((-rhoA + rhoB)*ux*ux + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz)));
+    		feq2 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) - 
+    				0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
+    						0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx + nx*ux + ny*uy + nz*uz)*
+    						(-2.*chem*ux*ux + 0.1111111111111111*(-4.*chem + rhoB*(-2.*ux - 1.*ux*ux - 1.*uy*uy - 1.*uz*uz) + 
+    								rhoA*(2.*ux + ux*ux + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*ux*ux + 
+    										chem*(4.*ux + 2.*ux*ux + 2.*uy*uy + 2.*uz*uz)));
+    		feq3 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) - 
+    				0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
+    						0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny - nx*ux - ny*uy - nz*uz)*
+    						(2*chem*uy*uy - 0.3333333333333333*((-rhoA + rhoB)*uy*uy + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz)));
+    		feq4 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) - 
+    				0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
+    						0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny + nx*ux + ny*uy + nz*uz)*
+    						(-2.*chem*uy*uy + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 2.*uy - 1.*uy*uy - 1.*uz*uz) + 
+    								rhoA*(ux*ux + 2.*uy + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uy*uy + 
+    										chem*(2.*ux*ux + 4.*uy + 2.*uy*uy + 2.*uz*uz))); 
+    		feq5 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) - 
+    				0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)) - 0.0625*(nx*ux + ny*uy + nz*(-1. + uz))*
+    						(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (2. - 1.*uz)*uz) + 
+    								rhoA*(ux*ux + uy*uy + (-2. + uz)*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
+    										chem*(2.*ux*ux + 2.*uy*uy + uz*(-4. + 2.*uz)))); 
+    		feq6 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) - 
+    				0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))) - 0.0625*(nz + nx*ux + ny*uy + nz*uz)*
+    						(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (-2. - 1.*uz)*uz) + 
+    								rhoA*(ux*ux + uy*uy + uz*(2. + uz))) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
+    										chem*(2.*ux*ux + 2.*uy*uy + uz*(4. + 2.*uz)))); 
+    		feq7 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux + uy)*(ux + uy) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
+    				0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
+    						0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx + ny - nx*ux - ny*uy - nz*uz)*
+    						(2*chem*(ux + uy)*(ux + uy) + 0.3333333333333333*((rhoA - rhoB)*(ux + uy)*(ux + uy) - 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
+    		feq8 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
+    				0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
+    						0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(-(nx*(1 + ux)) - ny*(1 + uy) - nz*uz)*
+    						(2*chem*(ux + uy)*(ux + uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uy)*(ux + uy)) + 
+    								2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
+    		feq9 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
+    				0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
+    						0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx - nx*ux - ny*(1 + uy) - nz*uz)*
+    						(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
+    								2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
+    		feq10 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
+    				0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
+    						0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(ny - nx*(1 + ux) - ny*uy - nz*uz)*
+    						(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
+    								2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
+    		feq11 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux + uz)*(ux + uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
+    				0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
+    						0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nx + nz - nx*ux - ny*uy - nz*uz)*
+    						(2*chem*(ux + uz)*(ux + uz) + 0.3333333333333333*((rhoA - rhoB)*(ux + uz)*(ux + uz) - 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    								0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
+    		feq12 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
+    				0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
+    						0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*(1 + ux)) - ny*uy - nz*(1 + uz))*
+    						(2*chem*(ux + uz)*(ux + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uz)*(ux + uz)) + 
+    								2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
+    		feq13 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
+    				0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
+    						0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(nx - nx*ux - ny*uy - nz*(1 + uz))*
+    						(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
+    								2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
+    		feq14 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
+    				0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
+    						0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*(1 + ux) - ny*uy - nz*uz)*
+    						(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
+    								2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
+    		feq15 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) - 
+    				0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(ny + nz - nx*ux - ny*uy - nz*uz)*
+    						(2*chem*(uy + uz)*(uy + uz) + 0.3333333333333333*((rhoA - rhoB)*(uy + uz)*(uy + uz) - 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
+    								0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))); 
+    		feq16 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) - 
+    				0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*ux) - ny*(1 + uy) - nz*(1 + uz))*
+    						(2*chem*(uy + uz)*(uy + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy + uz)*(uy + uz)) + 
+    								2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))); 
+    		feq17 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) - 
+    				0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(ny - nx*ux - ny*uy - nz*(1 + uz))*
+    						(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
+    								2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))); 
+    		feq18 = 0.027777777777777776*p - 0.041666666666666664*rho0*
+    				(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) - 
+    				0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
+    						0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*ux - ny*(1 + uy) - nz*uz)*
+    						(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
+    								2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
+    								(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))); 
+
+    		//------------------------------------------------- BCK collison ------------------------------------------------------------//
+    		// q=0
+    		dist[n] = m0 - (m0-feq0)/tau + 0.25*(2*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 
+    				(mgx*ux + mgy*uy + mgz*uz)*(2*chem*(ux*ux + uy*uy + uz*uz) + 
+    						0.3333333333333333*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*uz))));
+
+    		// q = 1
+    		dist[1*Np+n] = m1 - (m1-feq1)/tau + 0.125*(2*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
+    				0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
+    				(mgx*(-1 + ux) + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz))));
+
+    		// q=2
+    		dist[2*Np+n] = m2 - (m2-feq2)/tau + 0.125*(2*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
+    				0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) + 
+    				(mgx + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(2*ux + ux*ux + uy*uy + uz*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*uz))));
+
+    		// q = 3
+    		dist[3*Np+n] = m3 - (m3-feq3)/tau + 0.125*(2*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy + 
+    				0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    				(mgx*ux + mgy*(-1 + uy) + mgz*uz)*(-2*chem*(uy*uy) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz))));
+
+    		// q = 4
+    		dist[4*Np+n] = m4 - (m4-feq4)/tau + 0.125*(2*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy + 
+    				0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    				(mgy + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uy*uy) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*uz))));
+
+    		// q = 5
+    		dist[5*Np+n] = m5 - (m5-feq5)/tau + 0.125*(2*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz + 
+    				0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    				(mgx*ux + mgy*uy + mgz*(-1 + uz))*(-2*chem*(uz*uz) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + (-2 + uz)*uz))));
+
+    		// q = 6
+    		dist[6*Np+n] = m6 - (m6-feq6)/tau + 0.125*(2*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz + 
+    				0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) + 
+    				(mgz + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uz*uz) + 
+    						0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + uz*(2 + uz))) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*(2 + uz)))));
+
+    		// q = 7
+    		dist[7*Np+n] = m7 - (m7-feq7)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
+    				(0.2222222222222222 + (ux + uy)*(ux + uy) - 
+    						0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    						(mgx*(-1 + ux) + mgy*(-1 + uy) + mgz*uz)*
+    						(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
+
+    		// q = 8
+    		dist[8*Np+n] = m8 - (m8-feq8)/tau + 0.0625*(2*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
+    				(-0.2222222222222222 - (ux + uy)*(ux + uy) + 
+    						0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    						(mgx + mgy + mgx*ux + mgy*uy + mgz*uz)*
+    						(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
+
+    		// q = 9
+    		dist[9*Np+n] = m9 - (m9-feq9)/tau + 0.0625*(2*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
+    				(-0.2222222222222222 - (ux - uy)*(ux - uy) + 
+    						0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    						(mgy + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
+    						(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
+
+    		// q = 10
+    		dist[10*Np+n] = m10 - (m10-feq10)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
+    				(-0.2222222222222222 - (ux - uy)*(ux - uy) + 
+    						0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    						(mgx*(1 + ux) + mgy*(-1 + uy) + mgz*uz)*
+    						(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
+
+    		// q = 11
+    		dist[11*Np+n] = m11 - (m11-feq11)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
+    				(0.2222222222222222 + (ux + uz)*(ux + uz) - 
+    						0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    						(mgx*(-1 + ux) + mgy*uy + mgz*(-1 + uz))*
+    						(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
+
+    		// q = 12
+    		dist[12*Np+n] = m12 - (m12-feq12)/tau + 0.0625*(2*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
+    				(-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))
+    				+ (mgx + mgz + mgx*ux + mgy*uy + mgz*uz)*
+    				(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
+    						(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
+
+    		// q = 13
+    		dist[13*Np+n] = m13 - (m13-feq13)/tau + 0.0625*(2*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
+    				(-0.2222222222222222 - (ux - uz)*(ux - uz) + 
+    						0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
+    						(mgz + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
+    						(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
+
+    		// q= 14
+    		dist[14*Np+n] = m14 - (m14-feq14)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
+    				(-0.2222222222222222 - (ux - uz)*(ux - uz) + 
+    						0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    						(mgx*(1 + ux) + mgy*uy + mgz*(-1 + uz))*
+    						(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
+
+    		// q = 15
+    		dist[15*Np+n] = m15 - (m15-feq15)/tau + 0.0625*(-2*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
+    				(0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))
+    				+ (mgx*ux + mgy*(-1 + uy) + mgz*(-1 + uz))*
+    				(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
+    						(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))));
+
+    		// q = 16
+    		dist[16*Np+n] = m16 - (m16-feq16)/tau + 0.0625*(2*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
+    				(-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))
+    				+ (mgy + mgz + mgx*ux + mgy*uy + mgz*uz)*
+    				(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
+    						(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))));
+
+    		// q = 17
+    		dist[17*Np+n] = m17 - (m17-feq17)/tau + 0.0625*(2*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
+    				(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))
+    				+ (mgz + mgx*ux + mgy*(-1 + uy) + mgz*uz)*
+    				(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
+    						(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 
+    						0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))));
+
+    		// q = 18
+    		dist[18*Np+n] = m18 - (m18-feq18)/tau + 0.0625*(2*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
+    				(-0.2222222222222222 - (uy - uz)*(uy - uz) + 
+    						0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
+    						(mgx*ux + mgy*(1 + uy) + mgz*(-1 + uz))*
+    						(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
+    								(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
+    								0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
+    		//----------------------------------------------------------------------------------------------------------------------------------------//
+
+
+    		// ----------------------------- compute phase field evolution ----------------------------------------
+    		//Normalize the Color Gradient
+    		C = sqrt(nx*nx+ny*ny+nz*nz);
+    		double ColorMag = C;
+    		if (C==0.0) ColorMag=1.0;
+    		nx = nx/ColorMag;
+    		ny = ny/ColorMag;
+    		nz = nz/ColorMag;		
+    		//compute surface tension-related parameter
+    		theta = M*4.5*(1-4.0*phi*phi)/W;
+
+    		//load distributions of phase field
+    		//q=0
+    		h0 = hq[n];
+    		//q=1
+    		h1 = hq[2*Np+n]; 
+
+    		//q=2
+    		h2 = hq[1*Np+n];  
+
+    		//q=3
+    		h3 = hq[4*Np+n];
+
+    		//q=4
+    		h4 = hq[3*Np+n];
+
+    		//q=5
+    		h5 = hq[6*Np+n];
+
+    		//q=6
+    		h6 = hq[5*Np+n];
+
+    		//-------------------------------- BGK collison for phase field ---------------------------------//
+    		// q = 0
+    		hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
+
+    		// q = 1
+    		hq[1*Np+n] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
+
+    		// q = 2
+    		hq[2*Np+n] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
+
+    		// q = 3
+    		hq[3*Np+n] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
+
+    		// q = 4
+    		hq[4*Np+n] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
+
+    		// q = 5
+    		hq[5*Np+n] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
+
+    		// q = 6
+    		hq[6*Np+n] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
+    		//........................................................................
+
+    		//Update velocity on device
+    		Vel[0*Np+n] = ux;
+    		Vel[1*Np+n] = uy;
+    		Vel[2*Np+n] = uz;
+    		//Update pressure on device
+    		Pressure[n] = p;
+    		//Update chemical potential on device
+    		mu_phi[n] = chem;
+    		//Update color gradient on device
+    		ColorGrad[0*Np+n] = nx;
+    		ColorGrad[1*Np+n] = ny;
+    		ColorGrad[2*Np+n] = nz;
+    	}
+
+	}
+}
+
+__global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure,  
+		double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){
+
+	int n;
+	int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
+	double ux,uy,uz;//fluid velocity 
+	double p;//pressure
+	// distribution functions
+	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
+	double m0,m3,m5,m7;
+
+	//	for (int n=start; n<finish; 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 + start;
+
+		if ( n<finish ){
+			// q=0
+			m0 = dist[n];
+			// q=1
+			nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
+			m1 = dist[nr1]; // reading the f1 data into register fq
+
+			nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
+			m2 = dist[nr2];  // reading the f2 data into register fq
+
+			// q=3
+			nr3 = neighborList[n+2*Np]; // neighbor 4
+			m3 = dist[nr3];
+
+			// q = 4
+			nr4 = neighborList[n+3*Np]; // neighbor 3
+			m4 = dist[nr4];
+
+			// q=5
+			nr5 = neighborList[n+4*Np];
+			m5 = dist[nr5];
+
+			// q = 6
+			nr6 = neighborList[n+5*Np];
+			m6 = dist[nr6];
+
+			// q=7
+			nr7 = neighborList[n+6*Np];
+			m7 = dist[nr7];
+
+			// q = 8
+			nr8 = neighborList[n+7*Np];
+			m8 = dist[nr8];
+
+			// q=9
+			nr9 = neighborList[n+8*Np];
+			m9 = dist[nr9];
+
+			// q = 10
+			nr10 = neighborList[n+9*Np];
+			m10 = dist[nr10];
+
+			// q=11
+			nr11 = neighborList[n+10*Np];
+			m11 = dist[nr11];
+
+			// q=12
+			nr12 = neighborList[n+11*Np];
+			m12 = dist[nr12];
+
+			// q=13
+			nr13 = neighborList[n+12*Np];
+			m13 = dist[nr13];
+
+			// q=14
+			nr14 = neighborList[n+13*Np];
+			m14 = dist[nr14];
+
+			// q=15
+			nr15 = neighborList[n+14*Np];
+			m15 = dist[nr15];
+
+			// q=16
+			nr16 = neighborList[n+15*Np];
+			m16 = dist[nr16];
+
+			// q=17
+			nr17 = neighborList[n+16*Np];
+			m17 = dist[nr17];
+
+			// q=18
+			nr18 = neighborList[n+17*Np];
+			m18 = dist[nr18];
+
+			//compute fluid velocity
+			ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(Fx)/3.0);
+			uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(Fy)/3.0);
+			uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(Fz)/3.0);
+			//compute pressure
+			p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17);
+
+			//------------------------------------------------- BCK collison ------------------------------------------------------------//
+			// q=0
+			dist[n] = m0 + 0.5*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) - 
+					(m0 - 0.3333333333333333*p + 0.25*(Fx*ux + Fy*uy + Fz*uz)*
+							(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz))/
+							tau;
+
+			// q = 1
+			dist[nr2] = m1 + 0.25*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
+					0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) - 
+					(m1 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(ux*ux) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
+							0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) + 0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)))/tau;
+
+			// q=2
+			dist[nr1] = m2 + 0.25*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
+					0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) - 
+					(m2 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(ux*ux) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) + 
+							0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(ux*ux) + 
+									0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)))/tau;
+
+			// q = 3
+			dist[nr4] = m3 + 0.25*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy + 
+					0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) - 
+					(m3 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(uy*uy) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
+							0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*
+							(-0.2222222222222222 - 1.*(uy*uy) + 0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 4
+			dist[nr3] = m4 + 0.25*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy + 
+					0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) - 
+					(m4 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(uy*uy) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
+							0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uy*uy) + 
+									0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 5
+			dist[nr6] = m5 + 0.25*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz + 
+					0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) - 
+					(m5 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
+							0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*
+							(-0.2222222222222222 - 1.*(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
+
+			// q = 6
+			dist[nr5] = m6 + 0.25*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz + 
+					0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) - 
+					(m6 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) + 
+							0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uz*uz) + 
+									0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))))/tau;
+
+			// q = 7
+			dist[nr8] = m7 - 0.125*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
+					(0.2222222222222222 + (ux + uy)*(ux + uy) - 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
+					- (m7 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux + uy)*(ux + uy)) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+							0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) + 
+									0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 8
+			dist[nr7] = m8 + 0.125*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
+					(-0.2222222222222222 - (ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))\
+					- (m8 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux + uy)*(ux + uy)) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+							0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) + 
+									0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 9
+			dist[nr10] = m9 + 0.125*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
+					(-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))
+					- (m9 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux - uy)*(ux - uy)) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+							0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) + 
+									0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 10
+			dist[nr9] = m10 + 0.125*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
+					(-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
+					- (m10 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux - uy)*(ux - uy)) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+							0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) + 
+									0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 11
+			dist[nr12] = m11 - 0.125*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
+					(0.2222222222222222 + (ux + uz)*(ux + uz) - 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
+					- (m11 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux + uz)*(ux + uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+							0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*
+							(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
+
+			// q = 12
+			dist[nr11] = m12 + 0.125*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
+					(-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
+					- (m12 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux + uz)*(ux + uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
+							0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;
+
+			// q = 13
+			dist[nr14] = m13 + 0.125*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
+					(-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
+					- (m13 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux - uz)*(ux - uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
+							0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;
+
+			// q= 14
+			dist[nr13] = m14 + 0.125*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
+					(-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
+					- (m14 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux - uz)*(ux - uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+							0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*
+							(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
+
+			// q = 15
+			dist[nr16] = m15 - 0.125*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
+					(0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))\
+					- (m15 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
+							0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*
+							(-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;
+
+			// q = 16
+			dist[nr15] = m16 + 0.125*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
+					(-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))\
+					- (m16 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 
+							0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))))/tau;
+
+			// q = 17
+			dist[nr18] = m17 + 0.125*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
+					(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))\
+					- (m17 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 
+							0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*
+							(-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))))/tau;
+
+			// q = 18
+			dist[nr17] = m18 + 0.125*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
+					(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))\
+					- (m18 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
+							0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*
+							(-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;
+			//----------------------------------------------------------------------------------------------------------------------------------------//
+
+
+			//Update velocity on device
+			Vel[0*Np+n] = ux;
+			Vel[1*Np+n] = uy;
+			Vel[2*Np+n] = uz;
+			//Update pressure on device
+			Pressure[n] = p;
+		}
+	}
+}
+
+__global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, double *Vel, double *Pressure, 
+		double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){
+
+	int n;
+	double ux,uy,uz;//fluid velocity 
+	double p;//pressure
+	// distribution functions
+	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
+	double m0,m3,m5,m7;
+
+	//	for (int n=start; n<finish; 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 + start;
+
+		if ( n<finish ){	
+			// q=0
+			m0 = dist[n];
+			// q=1
+			m1 = dist[2*Np+n]; 
+
+			// q=2
+			m2 = dist[1*Np+n];  
+
+			// q=3
+			m3 = dist[4*Np+n];
+
+			// q = 4
+			m4 = dist[3*Np+n];
+
+			// q=5
+			m5 = dist[6*Np+n];
+
+			// q = 6
+			m6 = dist[5*Np+n];
+
+			// q=7
+			m7 = dist[8*Np+n];
+
+			// q = 8
+			m8 = dist[7*Np+n];
+
+			// q=9
+			m9 = dist[10*Np+n];
+
+			// q = 10
+			m10 = dist[9*Np+n];
+
+			// q=11
+			m11 = dist[12*Np+n];
+
+			// q=12
+			m12 = dist[11*Np+n];
+
+			// q=13
+			m13 = dist[14*Np+n];
+
+			// q=14
+			m14 = dist[13*Np+n];
+
+			// q=15
+			m15 = dist[16*Np+n];
+
+			// q=16
+			m16 = dist[15*Np+n];
+
+			// q=17
+			m17 = dist[18*Np+n];
+
+			// q=18
+			m18 = dist[17*Np+n];
+
+			//compute fluid velocity
+			ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(Fx)/3.0);
+			uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(Fy)/3.0);
+			uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(Fz)/3.0);
+			//compute pressure
+			p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17);
+
+			//------------------------------------------------- BCK collison ------------------------------------------------------------//
+			// q=0
+			dist[n] = m0 + 0.5*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) - 
+					(m0 - 0.3333333333333333*p + 0.25*(Fx*ux + Fy*uy + Fz*uz)*
+							(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz))/
+							tau; 
+
+			// q = 1
+			dist[1*Np+n] = m1 + 0.25*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
+					0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) - 
+					(m1 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(ux*ux) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
+							0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) + 0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)))/tau;
+
+			// q=2
+			dist[2*Np+n] = m2 + 0.25*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
+					0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) - 
+					(m2 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(ux*ux) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) + 
+							0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(ux*ux) + 
+									0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)))/tau;
+
+			// q = 3
+			dist[3*Np+n] = m3 + 0.25*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy + 
+					0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) - 
+					(m3 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(uy*uy) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
+							0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*
+							(-0.2222222222222222 - 1.*(uy*uy) + 0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 4
+			dist[4*Np+n] = m4 + 0.25*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy + 
+					0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) - 
+					(m4 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(uy*uy) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
+							0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uy*uy) + 
+									0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 5
+			dist[5*Np+n] = m5 + 0.25*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz + 
+					0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) - 
+					(m5 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
+							0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*
+							(-0.2222222222222222 - 1.*(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
+
+			// q = 6
+			dist[6*Np+n] = m6 + 0.25*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz + 
+					0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) - 
+					(m6 - 0.05555555555555555*p + 0.08333333333333333*rho0*
+							(-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) + 
+							0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uz*uz) + 
+									0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))))/tau;
+
+			// q = 7
+			dist[7*Np+n] = m7 - 0.125*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
+					(0.2222222222222222 + (ux + uy)*(ux + uy) - 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
+					- (m7 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux + uy)*(ux + uy)) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+							0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) + 
+									0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 8
+			dist[8*Np+n] = m8 + 0.125*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
+					(-0.2222222222222222 - (ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))\
+					- (m8 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux + uy)*(ux + uy)) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+							0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) + 
+									0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 9
+			dist[9*Np+n] = m9 + 0.125*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
+					(-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))
+					- (m9 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux - uy)*(ux - uy)) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
+							0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) + 
+									0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 10
+			dist[10*Np+n] = m10 + 0.125*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
+					(-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
+					- (m10 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux - uy)*(ux - uy)) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
+							0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) + 
+									0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
+
+			// q = 11
+			dist[11*Np+n] = m11 - 0.125*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
+					(0.2222222222222222 + (ux + uz)*(ux + uz) - 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
+					- (m11 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux + uz)*(ux + uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+							0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*
+							(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
+
+			// q = 12
+			dist[12*Np+n] = m12 + 0.125*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
+					(-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
+					- (m12 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux + uz)*(ux + uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
+							0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;
+
+			// q = 13
+			dist[13*Np+n] = m13 + 0.125*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
+					(-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
+					- (m13 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux - uz)*(ux - uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
+							0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;
+
+			// q= 14
+			dist[14*Np+n] = m14 + 0.125*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
+					(-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
+					- (m14 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((ux - uz)*(ux - uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
+							0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*
+							(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
+
+			// q = 15
+			dist[15*Np+n] = m15 - 0.125*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
+					(0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))\
+					- (m15 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
+							0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*
+							(-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;
+
+			// q = 16
+			dist[16*Np+n] = m16 + 0.125*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
+					(-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))\
+					- (m16 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 
+							0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
+							(-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))))/tau;
+
+			// q = 17
+			dist[17*Np+n] = m17 + 0.125*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
+					(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))\
+					- (m17 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 
+							0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*
+							(-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))))/tau;
+
+			// q = 18
+			dist[18*Np+n] = m18 + 0.125*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
+					(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))\
+					- (m18 - 0.027777777777777776*p + 0.041666666666666664*rho0*
+							(-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
+							0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*
+							(-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) + 
+									0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;
+			//----------------------------------------------------------------------------------------------------------------------------------------//
+
+			//Update velocity on device
+			Vel[0*Np+n] = ux;
+			Vel[1*Np+n] = uy;
+			Vel[2*Np+n] = uz;
+			//Update pressure on device
+			Pressure[n] = p;
+		}
+	}
+}
+
+extern "C" void ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np){
+	
+}
+
+
+extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np){
+	
+}
+
+extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad, 
+                                                    double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
+	
+}
+extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, 
+                                                          double rhoA, double rhoB, int start, int finish, int Np){
+	
+}
+
+extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, 
+			                                               double rhoA, double rhoB, int start, int finish, int Np){
+	
+}
+
+extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, double *dist, double *hq, double *Den,	double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, 
+                                                double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 
+                                                int strideY, int strideZ, int start, int finish, int Np){
+	
+}
+	
+
+extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *hq, double *Den,	double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
+                                                double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 
+                                                int strideY, int strideZ, int start, int finish, int Np){
+	
+}
+
+extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure,  
+                                                                double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){
+	
+
+}
+
+extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, double *Vel, double *Pressure, 
+                                                                 double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){
+	
+}
+