diff --git a/models/multiphase/d2q9_pf_velocity/Dynamics.R b/models/multiphase/d2q9_pf_velocity/Dynamics.R
index d2319bb53..602ac8191 100644
--- a/models/multiphase/d2q9_pf_velocity/Dynamics.R
+++ b/models/multiphase/d2q9_pf_velocity/Dynamics.R
@@ -86,20 +86,20 @@ if (Options$RT) {
 	AddStage("BaseInit"  , "Init_distributions" , save=Fields$group %in% c("g","h","Vel"))
 
 	# iteration
-	AddStage("BaseIter"  , "calcHydroIter"      , save=Fields$group %in% c("g","h","Vel","nw") , load=DensityAll$group %in% c("g","h","Vel","nw"))  # TODO: is nw needed here?
+	AddStage("BaseIter"  , "calcHydroIter"      , save=Fields$group %in% c("g","h","Vel","nw") , load=DensityAll$group %in% c("PF","g","h","Vel","nw"))  # TODO: is nw needed here?
 	AddStage("PhaseIter" , "calcPhaseFIter"		, save=Fields$group %in% c("PF")			   , load=DensityAll$group %in% c("g","h","Vel","nw"))
-	AddStage("WallIter"  , "calcWallPhaseIter"	, save=Fields$group %in% c("PF")			   , load=DensityAll$group %in% c("nw"))	
+	AddStage("WallIter"  , "calcWallPhaseIter"	, save=Fields$group %in% c("PF")			   , load=DensityAll$group %in% c("nw","PF"))	# Purposefully read/write of PF for boundary. complex geom may force RACE condition.
 }
 
 AddAction("Iteration", c("BaseIter", "PhaseIter","WallIter"))
 AddAction("Init"     , c("PhaseInit","WallInit", "WallIter","BaseInit"))
 
 # 	Outputs:
-AddQuantity(name="Rho",	  unit="kg/m3")
-AddQuantity(name="PhaseField",unit="1")
-AddQuantity(name="U",	  unit="m/s",vector=T)
-AddQuantity(name="NormalizedPressure",	  unit="Pa")
-AddQuantity(name="Pressure",	  unit="Pa")
+AddQuantity(name="Rho",	unit="kg/m3")
+AddQuantity(name="PhaseField", unit="1")
+AddQuantity(name="U", unit="m/s",vector=T)
+AddQuantity(name="NormalizedPressure", unit="1")
+AddQuantity(name="Pressure", unit="Pa")
 AddQuantity(name="Normal", unit="1", vector=T)
 
 #	Initialisation States
@@ -209,5 +209,4 @@ if (Options$Outflow) {
 }
 AddNodeType(name="Solid", group="BOUNDARY")
 AddNodeType(name="Wall", group="BOUNDARY")
-AddNodeType(name="BGK", group="COLLISION")
 AddNodeType(name="MRT", group="COLLISION")
diff --git a/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt b/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt
index 96e0cd550..d53213371 100755
--- a/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt
+++ b/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt
@@ -40,6 +40,7 @@
 
 #include <math.h>
 #define PI 3.1415926535897
+#define cs2 0.33333333
 
 <?R
 g = PV(Density$name[Density$group == "g"])
@@ -104,8 +105,25 @@ CudaDeviceFunction vector_t getNormal(){
 // 	HELPER FUNCTIONS:
 CudaDeviceFunction vector_t calcGradPhi(){
 	vector_t gradPhi;
-	// eq 34 from "Improved locality of the phase-field lattice Boltzmann 
-	//	model for immiscible fluids at high density ratios"
+	// eq 34 from "Improved locality of the phase-field lattice Boltzmann model for immiscible fluids at high density ratios"
+
+	gradPhi.x = 0.0;
+	gradPhi.y = 0.0;
+
+	real_t phi_a = PhaseF(0,0);
+	real_t phi_b, denom = 0.0;
+
+	for (int i = 1; i<9; i++) {
+		phi_b = PhaseF_dyn(d2q9_ex[i],d2q9_ey[i]);
+		if (!isnan(phi_b)){
+			gradPhi.x += wf[i]*d2q9_ex[i]*(phi_b-phi_a);
+			gradPhi.y += wf[i]*d2q9_ey[i]*(phi_b-phi_a);
+			denom += wf[i];
+		}
+	}
+	gradPhi.x = gradPhi.x/denom/cs2;
+	gradPhi.y = gradPhi.y/denom/cs2;
+
 	#ifdef OPTIONS_Outflow
 		if ((NodeType & NODE_BOUNDARY) == NODE_Neumann_E || (NodeType & NODE_BOUNDARY) == NODE_Convective_E) {
 			gradPhi.x = 0.0;
@@ -113,30 +131,11 @@ CudaDeviceFunction vector_t calcGradPhi(){
 		} else if ((NodeType & NODE_BOUNDARY) == NODE_Convective_N) {
 			gradPhi.x = (PhaseF(1,0) - PhaseF(-1,0))/3.0 + (PhaseF(1,-1) - PhaseF(-1,-1) )/6.0;
 			gradPhi.y = 0.0;
-		} else if ((NodeType & NODE_BOUNDARY) == NODE_Wall){
-			gradPhi.x = (PhaseF(1,0) - PhaseF(-1,0))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(1,-1) - PhaseF(-1,1))/12.0;
-			gradPhi.y = (PhaseF(0,1) - PhaseF(0,-1))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(-1,1) - PhaseF(1,-1))/12.0;
-		} else if (NodeType & NODE_BOUNDARY){
-			gradPhi.x = 0.0;
-			gradPhi.y = 0.0;
-		} else { 
-			gradPhi.x = (PhaseF(1,0) - PhaseF(-1,0))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(1,-1) - PhaseF(-1,1))/12.0;
-			gradPhi.y = (PhaseF(0,1) - PhaseF(0,-1))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(-1,1) - PhaseF(1,-1))/12.0;
-		}
-	#else
-		if ((NodeType & NODE_BOUNDARY) == NODE_Wall){
-			gradPhi.x = (PhaseF(1,0) - PhaseF(-1,0))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(1,-1) - PhaseF(-1,1))/12.0;
-			gradPhi.y = (PhaseF(0,1) - PhaseF(0,-1))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(-1,1) - PhaseF(1,-1))/12.0;
-		} else if (NodeType & NODE_BOUNDARY){
-			gradPhi.x = 0.0;
-			gradPhi.y = 0.0;
-		} else { 
-			gradPhi.x = (PhaseF(1,0) - PhaseF(-1,0))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(1,-1) - PhaseF(-1,1))/12.0;
-			gradPhi.y = (PhaseF(0,1) - PhaseF(0,-1))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(-1,1) - PhaseF(1,-1))/12.0;
-		}
+		} 
 	#endif
-		gradPhi.z = sqrt(gradPhi.x*gradPhi.x + gradPhi.y*gradPhi.y + 1e-12); // length
-		return gradPhi;
+
+	gradPhi.z = sqrt(gradPhi.x*gradPhi.x + gradPhi.y*gradPhi.y + 1e-12); // length
+	return gradPhi;
 }
 
 CudaDeviceFunction real_t get_lpPhi(real_t myPhase){
@@ -183,8 +182,7 @@ CudaDeviceFunction real_t calcMu(real_t myPhase){
 		}
 	#endif
 
-	// eq 5 from "Improved locality of the phase-field lattice Boltzmann 
-	// model for immiscible fluids at high density ratios"
+	// eq 5 from "Improved locality of the phase-field lattice Boltzmann model for immiscible fluids at high density ratios"
 	real_t pfavg = 0.5*(PhaseField_l+PhaseField_h);
 	real_t mu = 4.0*(12.0*sigma/W)*(myPhase-PhaseField_l)*(myPhase-PhaseField_h)*(myPhase-pfavg)
 		        - (1.5 *sigma*W) * lpPhi; 
@@ -192,8 +190,7 @@ CudaDeviceFunction real_t calcMu(real_t myPhase){
 }
 
 CudaDeviceFunction real_t calcGamma(int i, real_t u, real_t v, real_t u2mag){
-	// eq 10 from "Improved locality of the phase-field lattice Boltzmann 
-	// model for immiscible fluids at high density ratios"
+	// eq 10 from "Improved locality of the phase-field lattice Boltzmann model for immiscible fluids at high density ratios"
 	real_t gamma, tmp;
 	tmp = (d2q9_ex[i]*u+d2q9_ey[i]*v);
 
@@ -202,8 +199,7 @@ CudaDeviceFunction real_t calcGamma(int i, real_t u, real_t v, real_t u2mag){
 }
 
 CudaDeviceFunction real_t calcF_phi(int i, real_t numerator, real_t nx, real_t ny){
-	// eq 7 from "Improved locality of the phase-field lattice Boltzmann 
-	// model for immiscible fluids at high density ratios"
+	// eq 7 from "Improved locality of the phase-field lattice Boltzmann model for immiscible fluids at high density ratios"
 
 	// nx: normalized PhaseField gradient in x direction
 	// this method calculates F_phi in [0-9] lattice coord
@@ -214,8 +210,7 @@ CudaDeviceFunction real_t calcF_phi(int i, real_t numerator, real_t nx, real_t n
 }
 
 CudaDeviceFunction vector_t calcF_phi_xy(vector_t gradPhi, real_t myPhaseF, real_t pfavg){	
-	// eq 7 from "Improved locality of the phase-field lattice Boltzmann 
-	// model for immiscible fluids at high density ratios"
+	// eq 7 from "Improved locality of the phase-field lattice Boltzmann model for immiscible fluids at high density ratios"
 
 	// this method calculates F_phi in x,y coordinates
 	real_t nx = gradPhi.x/gradPhi.z;  // GradPhi normalized in x, y direction
@@ -285,6 +280,7 @@ CudaDeviceFunction real_t calcRho(real_t myPhaseF)
 CudaDeviceFunction void Init_phase() {
 	PhaseF = PhaseField_init;// read initial PhaseField distribution from config.xml file
 
+	// Special initialisation cases
 	if ( Period > 0 ) {
 		if ( Wave > 0) {
 	    		real_t InterfacePoint = MidPoint + Perturbation*cos(2.0*PI*X/Period + PI);
@@ -304,21 +300,23 @@ CudaDeviceFunction void Init_phase() {
 	    PhaseF = 0.5 * (PhaseField_h + PhaseField_l)
 		   - 0.5 * (PhaseField_h - PhaseField_l) * BubbleType * tanh(2.0*(Ri - Radius)/W);
 	}
-	if ((NodeType & NODE_BOUNDARY) == NODE_Wall) PhaseF = -999;
+
+	// Wall initialisation to calculate normal in second stage of init action
+	if ((NodeType & NODE_BOUNDARY) == NODE_Wall) PhaseF = NAN;
 }
 
 CudaDeviceFunction void Init_wallNorm(){
 	if ((NodeType & NODE_BOUNDARY) == NODE_Wall) {
 		int solidFlags[9], i;
-		solidFlags[0] = PhaseF(0,0)  /-998;
-		solidFlags[1] = PhaseF(1,0)  /-998;
-		solidFlags[2] = PhaseF(0,1)  /-998;
-		solidFlags[3] = PhaseF(-1,0) /-998;
-		solidFlags[4] = PhaseF(0,-1) /-998;
-		solidFlags[5] = PhaseF(1,1)  /-998;
-		solidFlags[6] = PhaseF(-1,1 )/-998;
-		solidFlags[7] = PhaseF(-1,-1)/-998;
-		solidFlags[8] = PhaseF(1,-1) /-998;
+		solidFlags[0] = 1;
+		solidFlags[1] = isnan(PhaseF(1,0)) ? 1 : 0;
+		solidFlags[2] = isnan(PhaseF(0,1)) ? 1 : 0;
+		solidFlags[3] = isnan(PhaseF(-1,0)) ? 1 : 0;
+		solidFlags[4] = isnan(PhaseF(0,-1)) ? 1 : 0;
+		solidFlags[5] = isnan(PhaseF(1,1)) ? 1 : 0;
+		solidFlags[6] = isnan(PhaseF(-1,1)) ? 1 : 0;
+		solidFlags[7] = isnan(PhaseF(-1,-1)) ? 1 : 0;
+		solidFlags[8] = isnan(PhaseF(1,-1)) ? 1 : 0;
 
 		real_t myNorm[2]={0.0,0.0};
 		for (i = 0 ; i < 9; i++){
@@ -350,29 +348,25 @@ CudaDeviceFunction void Init_wallNorm(){
 
 CudaDeviceFunction void Init_distributions(){
 
-	// Find Gradients and normals:
 	real_t myPhaseF = PhaseF(0,0);
-
 	vector_t gradPhi = calcGradPhi();
 	real_t nx, ny;
-	// gradPhi.z stores the magnitude of the vector
-	nx = gradPhi.x/gradPhi.z;
+	nx = gradPhi.x/gradPhi.z; //gradPhi.z stores the magnitude of the vector
 	ny = gradPhi.y/gradPhi.z;
 
-	// Define Equilibrium, then initialise all da things
 	U = VelocityX;
 	V = VelocityY;
 	real_t mag = U*U + V*V;
 
-	real_t Gamma[9], F_phi[9];
+	real_t Gamma, F_phi;
 	real_t pfavg = 0.5*(PhaseField_l+PhaseField_h);
 	real_t tmp1 = (1.0 - 4.0*(myPhaseF - pfavg)*(myPhaseF - pfavg))/W;
 
 	for (int i=0; i< 9; i++){
-		Gamma[i] = calcGamma(i, U, V, mag);
-		F_phi[i] = calcF_phi(i, tmp1, nx, ny);
-		h[i] = myPhaseF * Gamma[i] - 0.5*F_phi[i];   // heq
-		g[i] = Gamma[i] - wf[i]; 				     //geq
+		Gamma = calcGamma(i, U, V, mag);
+		F_phi = calcF_phi(i, tmp1, nx, ny);
+		h[i] = myPhaseF * Gamma - 0.5*F_phi;     // heq
+		g[i] = Gamma - wf[i]; 				     //geq
 	}
 
 	#ifdef OPTIONS_RT
@@ -386,7 +380,6 @@ CudaDeviceFunction void Init_distributions(){
 				C(h_old, h)	}	?>
 		}
 	#endif
-	PhaseF = <?R C(sum(h)) ?>;
 }
 
 //	ITERATION:
@@ -425,26 +418,25 @@ CudaDeviceFunction void BC_Switcher ()
 }
 
 CudaDeviceFunction void calcHydroIter() {
-	real_t tmpPhase = 0.0;
 	if ((NodeType & NODE_ADDITIONALS) == NODE_Smoothing){
 		Init_distributions();
-	}
-	else{
+	} else {
 		PhaseF = PhaseF(0,0);
-		tmpPhase = PhaseF;
+
 		real_t myPhaseF = PhaseF;
-		real_t rho = calcRho(myPhaseF);
 		log_data_from_special_nodes();
+
 		#ifdef OPTIONS_debug	
+			real_t rho = calcRho(myPhaseF);
 			AddToMomentumX(U*rho);
 			AddToMomentumY(V*rho);
 		#endif
 
-		BC_Switcher();		
+		BC_Switcher();	
 
-                if (PhaseF <= 0.5) {
-//			AddToBubbleVelocityX( (1 - PhaseF)*U );
-//			AddToBubbleVelocityY( (1 - PhaseF)*V );
+        if (PhaseF <= 0.5) {
+			AddToBubbleVelocityX( (1 - PhaseF)*U );
+			AddToBubbleVelocityY( (1 - PhaseF)*V );
 			AddToBubbleLocationY( (1 - PhaseF)*Y );
 			AddToSumPhiGas( (1 - PhaseF) );
 		}
@@ -455,15 +447,10 @@ CudaDeviceFunction void calcHydroIter() {
 		#else
 			if (NodeType & NODE_MRT) { CollisionMRT();}
 		#endif
-                if (tmpPhase <= 0.5) {
-			AddToBubbleVelocityX( (1 - tmpPhase)*U );
-			AddToBubbleVelocityY( (1 - tmpPhase)*V );
-		}
 
-		myPhaseF = PhaseF(0,0);
-		rho = calcRho(myPhaseF);
-	
 		#ifdef OPTIONS_debug
+			myPhaseF = PhaseF(0,0);
+			rho = calcRho(myPhaseF);			
 			AddToMomentumX_afterCol(U*rho);
 			AddToMomentumY_afterCol(V*rho);
 		#endif
@@ -505,465 +492,381 @@ CudaDeviceFunction void calcPhaseFIter(){
 }
 
 CudaDeviceFunction void calcWallPhaseIter(){
-// Contact angle is defined with respect to the high density fluid
+	// Contact angle is defined with respect to the high density fluid
 	PhaseF = PhaseF(0,0);
 
 	if ((NodeType & NODE_BOUNDARY) == NODE_Wall) {
 		if (nw_x == 0 && nw_y == 0) { 
-			PhaseF = 0.0; // TODO: why not simply PhaseF = PhaseF(0,0);
+			PhaseF = NAN;
+			//printf("EDGE CASE: Solid found with no normal at (%.1f %.1f)\n", X, Y);
 			return;
 		} 
-		real_t myA, myH, myPhase;
-
+		real_t myA, myPhase; 
 		myPhase = PhaseF_dyn(nw_x, nw_y);
-		myH = sqrt(nw_x*nw_x + nw_y*nw_y );
-
-		if (fabs(radAngle - PI/2) < 1e-4) { PhaseF = myPhase; } 
-		else {   
-			myA = 1.0 - 0.5*myH * (4.0/W)  * cos( radAngle );
+		if (isnan(myPhase)) {
+			//printf("EDGE CASE: Solid normal pointing into solid node at (%.1f %.1f)\n", X, Y);
+			PhaseF=0.0;
+			return;
+		}
+		if (fabs(radAngle - PI/2) < 1e-4) { 
+			PhaseF = myPhase; 
+		} else {   
+			myA = 1.0 - 0.5 * sqrt(nw_x*nw_x + nw_y*nw_y ) * (4.0/W)  * cos( radAngle );
 			PhaseF = (myA - sqrt( myA*myA - 4.0 * (myA-1.0)*myPhase))/(myA-1.0) - myPhase;
 		}
-	} 	 
+	} 	
 }
 
 #ifdef OPTIONS_CM
-	CudaDeviceFunction void relax_CM_hydro(real_t f_in[9], real_t tau, vector_t u)
-	{
-		real_t uxuy = u.x*u.y;
-		real_t ux2 = u.x*u.x;
-		real_t uy2 = u.y*u.y;
-		
-		real_t s_v = 1./tau;
-		real_t s_b = omega_bulk;
-		
-		real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8];
-		
-		real_t temp[9]; real_t cm_eq[9];
-		for (int i = 0; i < 9; i++) {
-			temp[i] = f_in[i];}
-		
-		//raw moments from density-probability functions
-		//[m00, m10, m01, m20, m02, m11, m21, m12, m22]
-		f_in[0] = m00;
-		f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
-		f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
-		f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
-		f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
-		f_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
-		f_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
-		f_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
-		f_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
-		
-		//central moments from raw moments
-		temp[0] = f_in[0];
-		temp[1] = -f_in[0]*u.x + f_in[1];
-		temp[2] = -f_in[0]*u.y + f_in[2];
-		temp[3] = f_in[0]*ux2 - 2*f_in[1]*u.x + f_in[3];
-		temp[4] = f_in[0]*uy2 - 2*f_in[2]*u.y + f_in[4];
-		temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5];
-		temp[6] = -f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 - f_in[3]*u.y - 2*f_in[5]*u.x + f_in[6];
-		temp[7] = -f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy - f_in[4]*u.x - 2*f_in[5]*u.y + f_in[7];
-		temp[8] = f_in[0]*ux2*uy2 - 2*f_in[1]*u.x*uy2 - 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy - 2*f_in[6]*u.y - 2*f_in[7]*u.x + f_in[8];
-		
-		//collision in central moments space
-		//calculate equilibrium distributions in cm space
-		cm_eq[0] = m00;
-		cm_eq[1] = -u.x*(m00 - 1.0);
-		cm_eq[2] = -u.y*(m00 - 1.0);
-		cm_eq[3] = m00*ux2 + 1./3.*m00 - ux2;
-		cm_eq[4] = m00*uy2 + 1./3.*m00 - uy2;
-		cm_eq[5] = uxuy*(m00 - 1.0);
-		cm_eq[6] = -u.y*(m00*ux2 + 1./3.*m00 - ux2 - 1./3.);
-		cm_eq[7] = -u.x*(m00*uy2 + 1./3.*m00 - uy2 - 1./3.);
-		cm_eq[8] = m00*ux2*uy2 + 1./3.*m00*ux2 + 1./3.*m00*uy2 + 1./9.*m00 - ux2*uy2 - 1./3.*ux2 - 1./3.*uy2;
-		//collide eq: -S*(cm - cm_eq)
-		f_in[0] = cm_eq[0] - temp[0];
-		f_in[1] = cm_eq[1] - temp[1];
-		f_in[2] = cm_eq[2] - temp[2];
-		f_in[3] = 0.5*(cm_eq[3] - temp[3])*(s_b + s_v) + 0.5*(cm_eq[4] - temp[4])*(s_b - s_v);
-		f_in[4] = 0.5*(cm_eq[3] - temp[3])*(s_b - s_v) + 0.5*(cm_eq[4] - temp[4])*(s_b + s_v);
-		f_in[5] = s_v*(cm_eq[5] - temp[5]);
-		f_in[6] = cm_eq[6] - temp[6];
-		f_in[7] = cm_eq[7] - temp[7];
-		f_in[8] = cm_eq[8] - temp[8];
-		
-		//back to raw moments
-		temp[0] = f_in[0];
-		temp[1] = f_in[0]*u.x + f_in[1];
-		temp[2] = f_in[0]*u.y + f_in[2];
-		temp[3] = f_in[0]*ux2 + 2*f_in[1]*u.x + f_in[3];
-		temp[4] = f_in[0]*uy2 + 2*f_in[2]*u.y + f_in[4];
-		temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5];
-		temp[6] = f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2*f_in[5]*u.x + f_in[6];
-		temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy + f_in[4]*u.x + 2*f_in[5]*u.y + f_in[7];
-		temp[8] = f_in[0]*ux2*uy2 + 2*f_in[1]*u.x*uy2 + 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy + 2*f_in[6]*u.y + 2*f_in[7]*u.x + f_in[8];
-		
-		//back to density-probability functions
-		f_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
-		f_in[1] = temp[1]/2 + temp[3]/2 - temp[7]/2 - temp[8]/2;
-		f_in[2] = temp[2]/2 + temp[4]/2 - temp[6]/2 - temp[8]/2;
-		f_in[3] = -temp[1]/2 + temp[3]/2 + temp[7]/2 - temp[8]/2;
-		f_in[4] = -temp[2]/2 + temp[4]/2 + temp[6]/2 - temp[8]/2;
-		f_in[5] = temp[5]/4 + temp[6]/4 + temp[7]/4 + temp[8]/4;
-		f_in[6] = -temp[5]/4 + temp[6]/4 - temp[7]/4 + temp[8]/4;
-		f_in[7] = temp[5]/4 - temp[6]/4 - temp[7]/4 + temp[8]/4;
-		f_in[8] = -temp[5]/4 - temp[6]/4 + temp[7]/4 + temp[8]/4;
-	}
+CudaDeviceFunction void relax_CM_hydro(real_t f_in[9], real_t tau, vector_t u)
+{
+	real_t uxuy = u.x*u.y;
+	real_t ux2 = u.x*u.x;
+	real_t uy2 = u.y*u.y;
+	
+	real_t s_v = 1./tau;
+	real_t s_b = omega_bulk;
+	
+	real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8];
+	
+	real_t temp[9]; real_t cm_eq[9];
+	for (int i = 0; i < 9; i++) {
+		temp[i] = f_in[i];}
+	
+	//raw moments from density-probability functions
+	//[m00, m10, m01, m20, m02, m11, m21, m12, m22]
+	f_in[0] = m00;
+	f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
+	f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
+	f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
+	f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
+	f_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
+	f_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
+	f_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
+	f_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
+	
+	//central moments from raw moments
+	temp[0] = f_in[0];
+	temp[1] = -f_in[0]*u.x + f_in[1];
+	temp[2] = -f_in[0]*u.y + f_in[2];
+	temp[3] = f_in[0]*ux2 - 2*f_in[1]*u.x + f_in[3];
+	temp[4] = f_in[0]*uy2 - 2*f_in[2]*u.y + f_in[4];
+	temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5];
+	temp[6] = -f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 - f_in[3]*u.y - 2*f_in[5]*u.x + f_in[6];
+	temp[7] = -f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy - f_in[4]*u.x - 2*f_in[5]*u.y + f_in[7];
+	temp[8] = f_in[0]*ux2*uy2 - 2*f_in[1]*u.x*uy2 - 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy - 2*f_in[6]*u.y - 2*f_in[7]*u.x + f_in[8];
+	
+	//collision in central moments space
+	//calculate equilibrium distributions in cm space
+	cm_eq[0] = m00;
+	cm_eq[1] = -u.x*(m00 - 1.0);
+	cm_eq[2] = -u.y*(m00 - 1.0);
+	cm_eq[3] = m00*ux2 + 1./3.*m00 - ux2;
+	cm_eq[4] = m00*uy2 + 1./3.*m00 - uy2;
+	cm_eq[5] = uxuy*(m00 - 1.0);
+	cm_eq[6] = -u.y*(m00*ux2 + 1./3.*m00 - ux2 - 1./3.);
+	cm_eq[7] = -u.x*(m00*uy2 + 1./3.*m00 - uy2 - 1./3.);
+	cm_eq[8] = m00*ux2*uy2 + 1./3.*m00*ux2 + 1./3.*m00*uy2 + 1./9.*m00 - ux2*uy2 - 1./3.*ux2 - 1./3.*uy2;
+	//collide eq: -S*(cm - cm_eq)
+	f_in[0] = cm_eq[0] - temp[0];
+	f_in[1] = cm_eq[1] - temp[1];
+	f_in[2] = cm_eq[2] - temp[2];
+	f_in[3] = 0.5*(cm_eq[3] - temp[3])*(s_b + s_v) + 0.5*(cm_eq[4] - temp[4])*(s_b - s_v);
+	f_in[4] = 0.5*(cm_eq[3] - temp[3])*(s_b - s_v) + 0.5*(cm_eq[4] - temp[4])*(s_b + s_v);
+	f_in[5] = s_v*(cm_eq[5] - temp[5]);
+	f_in[6] = cm_eq[6] - temp[6];
+	f_in[7] = cm_eq[7] - temp[7];
+	f_in[8] = cm_eq[8] - temp[8];
+	
+	//back to raw moments
+	temp[0] = f_in[0];
+	temp[1] = f_in[0]*u.x + f_in[1];
+	temp[2] = f_in[0]*u.y + f_in[2];
+	temp[3] = f_in[0]*ux2 + 2*f_in[1]*u.x + f_in[3];
+	temp[4] = f_in[0]*uy2 + 2*f_in[2]*u.y + f_in[4];
+	temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5];
+	temp[6] = f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2*f_in[5]*u.x + f_in[6];
+	temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy + f_in[4]*u.x + 2*f_in[5]*u.y + f_in[7];
+	temp[8] = f_in[0]*ux2*uy2 + 2*f_in[1]*u.x*uy2 + 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy + 2*f_in[6]*u.y + 2*f_in[7]*u.x + f_in[8];
+	
+	//back to density-probability functions
+	f_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
+	f_in[1] = temp[1]/2 + temp[3]/2 - temp[7]/2 - temp[8]/2;
+	f_in[2] = temp[2]/2 + temp[4]/2 - temp[6]/2 - temp[8]/2;
+	f_in[3] = -temp[1]/2 + temp[3]/2 + temp[7]/2 - temp[8]/2;
+	f_in[4] = -temp[2]/2 + temp[4]/2 + temp[6]/2 - temp[8]/2;
+	f_in[5] = temp[5]/4 + temp[6]/4 + temp[7]/4 + temp[8]/4;
+	f_in[6] = -temp[5]/4 + temp[6]/4 - temp[7]/4 + temp[8]/4;
+	f_in[7] = temp[5]/4 - temp[6]/4 - temp[7]/4 + temp[8]/4;
+	f_in[8] = -temp[5]/4 - temp[6]/4 + temp[7]/4 + temp[8]/4;
+}
 
-	CudaDeviceFunction void relax_and_collide_CM_hydro(real_t f_in[9], real_t tau, vector_t u, vector_t Fhydro, real_t rho) 
-	{
-		real_t uxuy = u.x*u.y;
-		real_t ux2 = u.x*u.x;
-		real_t uy2 = u.y*u.y;
-
-		real_t s_v = 1./tau;
-		real_t s_b = omega_bulk;
-		real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8];
-		real_t temp[9]; real_t cm_eq[9];
-		for (int i = 0; i < 9; i++) {
-			temp[i] = f_in[i];}
-		//raw moments from density-probability functions
-		//[m00, m10, m01, m20, m02, m11, m21, m12, m22]
-		f_in[0] = m00;
-		f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
-		f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
-		f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
-		f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
-		f_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
-		//central moments from raw moments
-		temp[3] = f_in[0]*ux2 - 2.0*f_in[1]*u.x + f_in[3];
-		temp[4] = f_in[0]*uy2 - 2.0*f_in[2]*u.y + f_in[4];
-		temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5];	
-		//collision in central moments space
-		//calculate equilibrium distributions in cm space
-		cm_eq[0] = m00;
-		cm_eq[1] = -u.x*(m00 - 1.0);
-		cm_eq[2] = -u.y*(m00 - 1.0);
-		cm_eq[3] = m00*ux2 + 1./3.*m00 - ux2;
-		cm_eq[4] = m00*uy2 + 1./3.*m00 - uy2;
-		cm_eq[5] = uxuy*(m00 - 1.0);
-		cm_eq[6] = -u.y*(m00*ux2 + 1./3.*m00 - ux2 - 1./3.);
-		cm_eq[7] = -u.x*(m00*uy2 + 1./3.*m00 - uy2 - 1./3.);
-		cm_eq[8] = m00*ux2*uy2 + 1./3.*m00*ux2 + 1./3.*m00*uy2 + 1./9.*m00 - ux2*uy2 - 1./3.*ux2 - 1./3.*uy2;
-
-		//calculate forces in cm space
-		//He et al. scheme: get_continuous_Maxwellian_DF(dzeta=dzeta, psi=1, u=(ux, uy)) - velocity term
-		//collide eq: (eye(9)-S)*cm + S*cm_eq + (eye(9)-S/2.)*force_in_cm_space // SOI
-		f_in[0] = cm_eq[0];
-		f_in[1] = 0.5*Fhydro.x/rho + cm_eq[1];
-		f_in[2] = 0.5*Fhydro.y/rho + cm_eq[2];
-		f_in[3] = 0.5*cm_eq[3]*(s_b + s_v) + 0.5*cm_eq[4]*(s_b - s_v) - temp[3]*(0.5*s_b + 0.5*s_v - 1.0) - 0.5*temp[4]*(s_b - s_v);
-		f_in[4] = 0.5*cm_eq[3]*(s_b - s_v) + 0.5*cm_eq[4]*(s_b + s_v) - 0.5*temp[3]*(s_b - s_v) - temp[4]*(0.5*s_b + 0.5*s_v - 1.0);
-		f_in[5] = cm_eq[5]*s_v - temp[5]*(s_v - 1.0);
-		f_in[6] = 1./6.*Fhydro.y/rho + cm_eq[6];
-		f_in[7] = 1./6.*Fhydro.x/rho + cm_eq[7];
-		f_in[8] = cm_eq[8];
-
-		//back to raw moments
-		temp[0] = f_in[0];
-		temp[1] = f_in[0]*u.x + f_in[1];
-		temp[2] = f_in[0]*u.y + f_in[2];
-		temp[3] = f_in[0]*ux2 + 2.0*f_in[1]*u.x + f_in[3];
-		temp[4] = f_in[0]*uy2 + 2.0*f_in[2]*u.y + f_in[4];
-		temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5];
-		temp[6] = f_in[0]*ux2*u.y + 2.0*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2.0*f_in[5]*u.x + f_in[6];
-		temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2.0*f_in[2]*uxuy + f_in[4]*u.x + 2.0*f_in[5]*u.y + f_in[7];
-		temp[8] = f_in[0]*ux2*uy2 + 2.0*f_in[1]*u.x*uy2 + 2.0*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4.0*f_in[5]*uxuy + 2.0*f_in[6]*u.y + 2.0*f_in[7]*u.x + f_in[8];
-		//back to density-probability functions
-		f_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
-		f_in[1] = 0.5*temp[1] + 0.5*temp[3] - 0.5*temp[7] - 0.5*temp[8];
-		f_in[2] = 0.5*temp[2] + 0.5*temp[4] - 0.5*temp[6] - 0.5*temp[8];
-		f_in[3] = -0.5*temp[1] + 0.5*temp[3] + 0.5*temp[7] - 0.5*temp[8];
-		f_in[4] = -0.5*temp[2] + 0.5*temp[4] + 0.5*temp[6] - 0.5*temp[8];
-		f_in[5] = 0.25*temp[5] + 0.25*temp[6] + 0.25*temp[7] + 0.25*temp[8];
-		f_in[6] = -0.25*temp[5] + 0.25*temp[6] - 0.25*temp[7] + 0.25*temp[8];
-		f_in[7] = 0.25*temp[5] - 0.25*temp[6] - 0.25*temp[7] + 0.25*temp[8];
-		f_in[8] = -0.25*temp[5] - 0.25*temp[6] + 0.25*temp[7] + 0.25*temp[8];
-	}
+CudaDeviceFunction void relax_and_collide_CM_hydro(real_t f_in[9], real_t tau, vector_t u, vector_t Fhydro, real_t rho) 
+{
+	real_t uxuy = u.x*u.y;
+	real_t ux2 = u.x*u.x;
+	real_t uy2 = u.y*u.y;
+
+	real_t s_v = 1./tau;
+	real_t s_b = omega_bulk;
+	real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8];
+	real_t temp[9]; real_t cm_eq[9];
+	for (int i = 0; i < 9; i++) {
+		temp[i] = f_in[i];}
+	//raw moments from density-probability functions
+	//[m00, m10, m01, m20, m02, m11, m21, m12, m22]
+	f_in[0] = m00;
+	f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
+	f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
+	f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
+	f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
+	f_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
+	//central moments from raw moments
+	temp[3] = f_in[0]*ux2 - 2.0*f_in[1]*u.x + f_in[3];
+	temp[4] = f_in[0]*uy2 - 2.0*f_in[2]*u.y + f_in[4];
+	temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5];	
+	//collision in central moments space
+	//calculate equilibrium distributions in cm space
+	cm_eq[0] = m00;
+	cm_eq[1] = -u.x*(m00 - 1.0);
+	cm_eq[2] = -u.y*(m00 - 1.0);
+	cm_eq[3] = m00*ux2 + 1./3.*m00 - ux2;
+	cm_eq[4] = m00*uy2 + 1./3.*m00 - uy2;
+	cm_eq[5] = uxuy*(m00 - 1.0);
+	cm_eq[6] = -u.y*(m00*ux2 + 1./3.*m00 - ux2 - 1./3.);
+	cm_eq[7] = -u.x*(m00*uy2 + 1./3.*m00 - uy2 - 1./3.);
+	cm_eq[8] = m00*ux2*uy2 + 1./3.*m00*ux2 + 1./3.*m00*uy2 + 1./9.*m00 - ux2*uy2 - 1./3.*ux2 - 1./3.*uy2;
+
+	//calculate forces in cm space
+	//He et al. scheme: get_continuous_Maxwellian_DF(dzeta=dzeta, psi=1, u=(ux, uy)) - velocity term
+	//collide eq: (eye(9)-S)*cm + S*cm_eq + (eye(9)-S/2.)*force_in_cm_space // SOI
+	f_in[0] = cm_eq[0];
+	f_in[1] = 0.5*Fhydro.x/rho + cm_eq[1];
+	f_in[2] = 0.5*Fhydro.y/rho + cm_eq[2];
+	f_in[3] = 0.5*cm_eq[3]*(s_b + s_v) + 0.5*cm_eq[4]*(s_b - s_v) - temp[3]*(0.5*s_b + 0.5*s_v - 1.0) - 0.5*temp[4]*(s_b - s_v);
+	f_in[4] = 0.5*cm_eq[3]*(s_b - s_v) + 0.5*cm_eq[4]*(s_b + s_v) - 0.5*temp[3]*(s_b - s_v) - temp[4]*(0.5*s_b + 0.5*s_v - 1.0);
+	f_in[5] = cm_eq[5]*s_v - temp[5]*(s_v - 1.0);
+	f_in[6] = 1./6.*Fhydro.y/rho + cm_eq[6];
+	f_in[7] = 1./6.*Fhydro.x/rho + cm_eq[7];
+	f_in[8] = cm_eq[8];
+
+	//back to raw moments
+	temp[0] = f_in[0];
+	temp[1] = f_in[0]*u.x + f_in[1];
+	temp[2] = f_in[0]*u.y + f_in[2];
+	temp[3] = f_in[0]*ux2 + 2.0*f_in[1]*u.x + f_in[3];
+	temp[4] = f_in[0]*uy2 + 2.0*f_in[2]*u.y + f_in[4];
+	temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5];
+	temp[6] = f_in[0]*ux2*u.y + 2.0*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2.0*f_in[5]*u.x + f_in[6];
+	temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2.0*f_in[2]*uxuy + f_in[4]*u.x + 2.0*f_in[5]*u.y + f_in[7];
+	temp[8] = f_in[0]*ux2*uy2 + 2.0*f_in[1]*u.x*uy2 + 2.0*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4.0*f_in[5]*uxuy + 2.0*f_in[6]*u.y + 2.0*f_in[7]*u.x + f_in[8];
+	//back to density-probability functions
+	f_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
+	f_in[1] = 0.5*temp[1] + 0.5*temp[3] - 0.5*temp[7] - 0.5*temp[8];
+	f_in[2] = 0.5*temp[2] + 0.5*temp[4] - 0.5*temp[6] - 0.5*temp[8];
+	f_in[3] = -0.5*temp[1] + 0.5*temp[3] + 0.5*temp[7] - 0.5*temp[8];
+	f_in[4] = -0.5*temp[2] + 0.5*temp[4] + 0.5*temp[6] - 0.5*temp[8];
+	f_in[5] = 0.25*temp[5] + 0.25*temp[6] + 0.25*temp[7] + 0.25*temp[8];
+	f_in[6] = -0.25*temp[5] + 0.25*temp[6] - 0.25*temp[7] + 0.25*temp[8];
+	f_in[7] = 0.25*temp[5] - 0.25*temp[6] - 0.25*temp[7] + 0.25*temp[8];
+	f_in[8] = -0.25*temp[5] - 0.25*temp[6] + 0.25*temp[7] + 0.25*temp[8];
+}
 
-	CudaDeviceFunction void relax_and_collide_CM_phase_field(real_t f_in[9], real_t tau, vector_t u, vector_t F_phi, real_t rho) 
-	{
-		real_t uxuy = u.x*u.y;
-		real_t ux2 = u.x*u.x;
-		real_t uy2 = u.y*u.y;
-
-		real_t s_v = 1./tau;
-		real_t s_b = omega_bulk;
-
-		real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8];
-
-		real_t temp[9];
-
-		for (int i = 0; i < 9; i++) {
-			temp[i] = f_in[i];}
-
-		//raw moments from density-probability functions
-		//[m00, m10, m01, m20, m02, m11, m21, m12, m22]
-		f_in[0] = m00;
-		f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
-		f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
-
-		// f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
-		// f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
-		// f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
-		// f_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
-		// f_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
-		// f_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
-		// f_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
-
-		// Central moments from raw moments
-		// temp[0] = m00;
-		temp[1] = -f_in[0]*u.x + f_in[1];
-		temp[2] = -f_in[0]*u.y + f_in[2];
-
-		// temp[3] = f_in[0]*ux2 - 2*f_in[1]*u.x + f_in[3];
-		// temp[4] = f_in[0]*uy2 - 2*f_in[2]*u.y + f_in[4];
-		// temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5];
-		// temp[6] = -f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 - f_in[3]*u.y - 2*f_in[5]*u.x + f_in[6];
-		// temp[7] = -f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy - f_in[4]*u.x - 2*f_in[5]*u.y + f_in[7];
-		// temp[8] = f_in[0]*ux2*uy2 - 2*f_in[1]*u.x*uy2 - 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy - 2*f_in[6]*u.y - 2*f_in[7]*u.x + f_in[8];
-
-		// Collision in central moments space
-		// Default is second order integration (trapezoidal discretization)
-		// collide eq: (eye(9)-S)*cm + S*cm_eq + (eye(9)-S/2.)*F_phi_cm
-		// relax 1st moments
-		f_in[0] = m00;
-		f_in[1] = -F_phi.x*(0.5*s_v - 1.0) - temp[1]*(s_v - 1.0);
-		f_in[2] = -F_phi.y*(0.5*s_v - 1.0) - temp[2]*(s_v - 1.0);
-		f_in[3] = 1./3.*m00;
-		f_in[4] = 1./3.*m00;
-		f_in[5] = 0;
-		f_in[6] = 1./6.*F_phi.y;
-		f_in[7] = 1./6.*F_phi.x;
-		f_in[8] = 1./9.*m00;
-
-		//back to raw moments
-		temp[0] = f_in[0];
-		temp[1] = f_in[0]*u.x + f_in[1];
-		temp[2] = f_in[0]*u.y + f_in[2];
-		temp[3] = f_in[0]*ux2 + 2*f_in[1]*u.x + f_in[3];
-		temp[4] = f_in[0]*uy2 + 2*f_in[2]*u.y + f_in[4];
-		temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5];
-		temp[6] = f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2*f_in[5]*u.x + f_in[6];
-		temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy + f_in[4]*u.x + 2*f_in[5]*u.y + f_in[7];
-		temp[8] = f_in[0]*ux2*uy2 + 2*f_in[1]*u.x*uy2 + 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy + 2*f_in[6]*u.y + 2*f_in[7]*u.x + f_in[8];
-
-		//back to density-probability functions
-		f_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
-		f_in[1] = temp[1]/2 + temp[3]/2 - temp[7]/2 - temp[8]/2;
-		f_in[2] = temp[2]/2 + temp[4]/2 - temp[6]/2 - temp[8]/2;
-		f_in[3] = -temp[1]/2 + temp[3]/2 + temp[7]/2 - temp[8]/2;
-		f_in[4] = -temp[2]/2 + temp[4]/2 + temp[6]/2 - temp[8]/2;
-		f_in[5] = temp[5]/4 + temp[6]/4 + temp[7]/4 + temp[8]/4;
-		f_in[6] = -temp[5]/4 + temp[6]/4 - temp[7]/4 + temp[8]/4;
-		f_in[7] = temp[5]/4 - temp[6]/4 - temp[7]/4 + temp[8]/4;
-		f_in[8] = -temp[5]/4 - temp[6]/4 + temp[7]/4 + temp[8]/4;
-	}
+CudaDeviceFunction void relax_and_collide_CM_phase_field(real_t f_in[9], real_t tau, vector_t u, vector_t F_phi, real_t rho) 
+{
+	real_t uxuy = u.x*u.y;
+	real_t ux2 = u.x*u.x;
+	real_t uy2 = u.y*u.y;
+
+	real_t s_v = 1./tau;
+	real_t s_b = omega_bulk;
+
+	real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8];
+
+	real_t temp[9];
+
+	for (int i = 0; i < 9; i++) {
+		temp[i] = f_in[i];}
+
+	//raw moments from density-probability functions
+	//[m00, m10, m01, m20, m02, m11, m21, m12, m22]
+	f_in[0] = m00;
+	f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8];
+	f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
+
+	// f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8];
+	// f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8];
+	// f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8];
+	// f_in[5] = temp[5] - temp[6] + temp[7] - temp[8];
+	// f_in[6] = temp[5] + temp[6] - temp[7] - temp[8];
+	// f_in[7] = temp[5] - temp[6] - temp[7] + temp[8];
+	// f_in[8] = temp[5] + temp[6] + temp[7] + temp[8];
+
+	// Central moments from raw moments
+	// temp[0] = m00;
+	temp[1] = -f_in[0]*u.x + f_in[1];
+	temp[2] = -f_in[0]*u.y + f_in[2];
+
+	// temp[3] = f_in[0]*ux2 - 2*f_in[1]*u.x + f_in[3];
+	// temp[4] = f_in[0]*uy2 - 2*f_in[2]*u.y + f_in[4];
+	// temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5];
+	// temp[6] = -f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 - f_in[3]*u.y - 2*f_in[5]*u.x + f_in[6];
+	// temp[7] = -f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy - f_in[4]*u.x - 2*f_in[5]*u.y + f_in[7];
+	// temp[8] = f_in[0]*ux2*uy2 - 2*f_in[1]*u.x*uy2 - 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy - 2*f_in[6]*u.y - 2*f_in[7]*u.x + f_in[8];
+
+	// Collision in central moments space
+	// Default is second order integration (trapezoidal discretization)
+	// collide eq: (eye(9)-S)*cm + S*cm_eq + (eye(9)-S/2.)*F_phi_cm
+	// relax 1st moments
+	f_in[0] = m00;
+	f_in[1] = -F_phi.x*(0.5*s_v - 1.0) - temp[1]*(s_v - 1.0);
+	f_in[2] = -F_phi.y*(0.5*s_v - 1.0) - temp[2]*(s_v - 1.0);
+	f_in[3] = 1./3.*m00;
+	f_in[4] = 1./3.*m00;
+	f_in[5] = 0;
+	f_in[6] = 1./6.*F_phi.y;
+	f_in[7] = 1./6.*F_phi.x;
+	f_in[8] = 1./9.*m00;
+
+	//back to raw moments
+	temp[0] = f_in[0];
+	temp[1] = f_in[0]*u.x + f_in[1];
+	temp[2] = f_in[0]*u.y + f_in[2];
+	temp[3] = f_in[0]*ux2 + 2*f_in[1]*u.x + f_in[3];
+	temp[4] = f_in[0]*uy2 + 2*f_in[2]*u.y + f_in[4];
+	temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5];
+	temp[6] = f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2*f_in[5]*u.x + f_in[6];
+	temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy + f_in[4]*u.x + 2*f_in[5]*u.y + f_in[7];
+	temp[8] = f_in[0]*ux2*uy2 + 2*f_in[1]*u.x*uy2 + 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy + 2*f_in[6]*u.y + 2*f_in[7]*u.x + f_in[8];
+
+	//back to density-probability functions
+	f_in[0] = temp[0] - temp[3] - temp[4] + temp[8];
+	f_in[1] = temp[1]/2 + temp[3]/2 - temp[7]/2 - temp[8]/2;
+	f_in[2] = temp[2]/2 + temp[4]/2 - temp[6]/2 - temp[8]/2;
+	f_in[3] = -temp[1]/2 + temp[3]/2 + temp[7]/2 - temp[8]/2;
+	f_in[4] = -temp[2]/2 + temp[4]/2 + temp[6]/2 - temp[8]/2;
+	f_in[5] = temp[5]/4 + temp[6]/4 + temp[7]/4 + temp[8]/4;
+	f_in[6] = -temp[5]/4 + temp[6]/4 - temp[7]/4 + temp[8]/4;
+	f_in[7] = temp[5]/4 - temp[6]/4 - temp[7]/4 + temp[8]/4;
+	f_in[8] = -temp[5]/4 - temp[6]/4 + temp[7]/4 + temp[8]/4;
+}
 
-	CudaDeviceFunction vector_t calcTotalHydrodynamicForceCM(vector_t gradPhi, real_t myPhaseF, real_t rho, real_t tau, real_t mu, real_t p)
-	{
-		// arguments:
-		// gradPhi - Phase field gradients
-		// myPhaseF - PhaseField value at cell of interest (0,0 or BC like inlet_PhaseF, moving_wall_PhaseF)
-
-		// Fluid Properties: //
-		// real_t rho - fluid density
-		// real_t tau - relaxation parameter
-		// real_t mu - chemical potential 
-		// real_t p - normalized pressure
-
-		real_t Gamma[9], geq[9], mag;	// equilibrium, pressure equilibrium, velocity magnitude
-		real_t F_pressure[2], F_body[2], F_mu[2], F_surf_tension[2];   // Forces
-		vector_t F_total_hydro;	
-
-		real_t stress[3];				// Stress tensor calculation
-		real_t R[9];				// Populations for MRT relaxation
-
-		// Force Calculations
-		// eq 19
-		real_t density_coeff = (Density_h-Density_l)/(PhaseField_h-PhaseField_l);
-		F_pressure[0] = (-1.0/3.0) * p * density_coeff * gradPhi.x;
-		F_pressure[1] = (-1.0/3.0) * p * density_coeff * gradPhi.y;	
-
-		F_body[0] = -1.0*(rho-Density_l)*BuoyancyX + rho*GravitationX;
-		F_body[1] = -1.0*(rho-Density_l)*BuoyancyY + rho*GravitationY;
-		
-		// eq 4
-		F_surf_tension[0] = mu*gradPhi.x;
-		F_surf_tension[1] = mu*gradPhi.y;
-
-		// Calculate viscous force
-		for (int j=0;j<fixedIterator;j++){
-			// User has the choice in attempt to improve stability to use
-			// a fixed point iteration style update for the viscous force.
-			for (int i=0; i< 9; i++){
-				R[i] = g[i];
-			}
-			
-			vector_t u;
-			u.x = U;
-			u.y = V;
-			relax_CM_hydro(R, tau, u); 
-
-			// Stress/strain Tensor
-			stress[0] = 0.0;stress[1] = 0.0;stress[2] = 0.0;
-			for (int i=0; i< 9; i++){
-				stress[0] -= R[i]*d2q9_ex[i]*d2q9_ex[i];
-				stress[1] -= R[i]*d2q9_ex[i]*d2q9_ey[i];
-				stress[2] -= R[i]*d2q9_ey[i]*d2q9_ey[i];    
-			}
-			
-			// Viscous force
-			F_mu[0] = (0.5-tau) * (stress[0]*gradPhi.x + stress[1]*gradPhi.y) * density_coeff;
-			F_mu[1] = (0.5-tau) * (stress[1]*gradPhi.x + stress[2]*gradPhi.y) * density_coeff;
+CudaDeviceFunction vector_t calcTotalHydrodynamicForceCM(vector_t gradPhi, real_t myPhaseF, real_t rho, real_t tau, real_t mu, real_t p)
+{
+	// arguments:
+	// gradPhi - Phase field gradients
+	// myPhaseF - PhaseField value at cell of interest (0,0 or BC like inlet_PhaseF, moving_wall_PhaseF)
 
-			// Total Hydrodynamic Force
-			F_total_hydro.x = F_surf_tension[0] + F_pressure[0] + F_body[0] + F_mu[0];
-			F_total_hydro.y = F_surf_tension[1] + F_pressure[1] + F_body[1] + F_mu[1];
-			F_total_hydro.z = 0;
+	// Fluid Properties: //
+	// real_t rho - fluid density
+	// real_t tau - relaxation parameter
+	// real_t mu - chemical potential 
+	// real_t p - normalized pressure
 
-			// update velocity and re-iterate
-			U = g[1]-g[3]+g[5]-g[6]-g[7]+g[8];
-			V = g[2]-g[4]+g[5]+g[6]-g[7]-g[8];
+	real_t Gamma[9], geq[9], mag;	// equilibrium, pressure equilibrium, velocity magnitude
+	real_t F_pressure[2], F_body[2], F_mu[2], F_surf_tension[2];   // Forces
+	vector_t F_total_hydro;	
 
-			U += 0.5*F_total_hydro.x/rho;
-			V += 0.5*F_total_hydro.y/rho;
-		}
-		#ifdef OPTIONS_debug
-			AddToF_pressureX(F_pressure[0]);
-			AddToF_pressureY(F_pressure[1]);
-			AddToF_bodyX(F_body[0]);
-			AddToF_bodyY(F_body[1]);
-			AddToF_surf_tensionX(F_surf_tension[0]);
-			AddToF_surf_tensionY(F_surf_tension[1]);
-			AddToF_muX(F_mu[0]);
-			AddToF_muY(F_mu[1]);
-			AddToF_total_hydroX(F_total_hydro.x);
-			AddToF_total_hydroY(F_total_hydro.y);
-		#endif
+	real_t stress[3];				// Stress tensor calculation
+	real_t R[9];				// Populations for MRT relaxation
 
-		return F_total_hydro;
-	}
+	// Force Calculations
+	// eq 19
+	real_t density_coeff = (Density_h-Density_l)/(PhaseField_h-PhaseField_l);
+	F_pressure[0] = (-1.0/3.0) * p * density_coeff * gradPhi.x;
+	F_pressure[1] = (-1.0/3.0) * p * density_coeff * gradPhi.y;	
 
-	CudaDeviceFunction void CollisionCM(){
-		#ifdef OPTIONS_RT
-			real_t uold = U;
-			real_t vold = V;
-		#endif
+	F_body[0] = -1.0*(rho-Density_l)*BuoyancyX + rho*GravitationX;
+	F_body[1] = -1.0*(rho-Density_l)*BuoyancyY + rho*GravitationY;
+	
+	// eq 4
+	F_surf_tension[0] = mu*gradPhi.x;
+	F_surf_tension[1] = mu*gradPhi.y;
 
-		real_t myPhaseF = PhaseF(0,0);
+	// Calculate viscous force
+	for (int j=0;j<fixedIterator;j++){
+		// User has the choice in attempt to improve stability to use
+		// a fixed point iteration style update for the viscous force.
+		for (int i=0; i< 9; i++){
+			R[i] = g[i];
+		}
+		
+		vector_t u;
+		u.x = U;
+		u.y = V;
+		relax_CM_hydro(R, tau, u); 
 
-		// Find Macroscopic Details
-		real_t rho = calcRho(myPhaseF);
-		real_t tau = calcTau(myPhaseF, rho);
-		real_t pfavg = 0.5*(PhaseField_l+PhaseField_h);
+		// Stress/strain Tensor
+		stress[0] = 0.0;stress[1] = 0.0;stress[2] = 0.0;
+		for (int i=0; i< 9; i++){
+			stress[0] -= R[i]*d2q9_ex[i]*d2q9_ex[i];
+			stress[1] -= R[i]*d2q9_ex[i]*d2q9_ey[i];
+			stress[2] -= R[i]*d2q9_ey[i]*d2q9_ey[i];    
+		}
+		
+		// Viscous force
+		F_mu[0] = (0.5-tau) * (stress[0]*gradPhi.x + stress[1]*gradPhi.y) * density_coeff;
+		F_mu[1] = (0.5-tau) * (stress[1]*gradPhi.x + stress[2]*gradPhi.y) * density_coeff;
 
-		// Gradients & Normals
-		vector_t gradPhi = calcGradPhi(); // Phase field gradient
-		real_t p = getNormalizedPressure(); // normalized pressure
+		// Total Hydrodynamic Force
+		F_total_hydro.x = F_surf_tension[0] + F_pressure[0] + F_body[0] + F_mu[0];
+		F_total_hydro.y = F_surf_tension[1] + F_pressure[1] + F_body[1] + F_mu[1];
+		F_total_hydro.z = 0;
 
-		vector_t Fhydro = calcTotalHydrodynamicForceCM(gradPhi, myPhaseF, rho, tau,  calcMu(myPhaseF), p);
+		// update velocity and re-iterate
+		U = g[1]-g[3]+g[5]-g[6]-g[7]+g[8];
+		V = g[2]-g[4]+g[5]+g[6]-g[7]-g[8];
 
-		// PHASE FIELD POPULATION UPDATE:
-		vector_t u; u.x = U; u.y = V;
-		vector_t F_phi_vect = calcF_phi_xy(gradPhi,  myPhaseF,  pfavg);
-		#ifdef OPTIONS_debug
-			AddToF_phiX(F_phi_vect.x);
-			AddToF_phiY(F_phi_vect.y);
-		#endif
-		relax_and_collide_CM_phase_field(h, 1./omega_phi, u, F_phi_vect, rho);
-		
-		// PRESSURE EVOLUTION UPDATE:
-		relax_and_collide_CM_hydro(g, tau, u, Fhydro, rho);
+		U += 0.5*F_total_hydro.x/rho;
+		V += 0.5*F_total_hydro.y/rho;
 	}
-#elif OPTIONS_BGK
-	CudaDeviceFunction void CollisionBGK(){
-		int i, j;
-		real_t C = PhaseF(0,0), mu;
-		real_t tau, DynVisc, rho, p;		// Macroscopic Properties
-		vector_t gradPhi;					// Phase field gradients
-		real_t nx, ny;						// Normals
-		real_t Gamma[9], geq[9], mag;		// equilibrium, pressure equilibrium, velocity magnitude
-		real_t F_pressure[2], F_body[2], F_mu[2], F_total[2]; // Forces
-		real_t tmp1, stress[3];				// Stress tensor calculation
-		real_t F_phi[9], heq[9];			// Phase field collision terms
-		real_t F_i[9];						// Momentum distribution forcing term
-
-		// Find Macroscopic Details
-		mu = calcMu( C );
-		rho = Density_l + (C - PhaseField_l)*(Density_h - Density_l);
-		p = g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6]+g[7]+g[8];
-
-		// Updating of tau:
-		//    Updating via the kinematic viscosity here gives more stable solutions
-		//    Updating via the dynVisc section of the code is more accurate
-		//        but will go unstable for low viscosities
-		if ( C < PhaseField_l){
-			tau = tau_l + 0.5;
-		} else if (C > PhaseField_h) {
-			tau = tau_h + 0.5;
-		} else {
-			// Inverse update:
-			//tau = (C - PhaseField_l)/(PhaseField_h - PhaseField_l) * (1.0/tau_h - 1.0/tau_l) + 1.0/tau_l;
-			//tau = 1.0/tau + 0.5;
-			// Linear update:
-			tau = 0.5 + tau_l + C*(tau_h - tau_l);
-			// Viscosity update:
-			//DynVisc = Density_l*Viscosity_l + C * (Density_h*Viscosity_h - Density_l*Viscosity_l);
-			//tau = 3.0 * DynVisc / rho + 0.5;
-		}
+	#ifdef OPTIONS_debug
+		AddToF_pressureX(F_pressure[0]);
+		AddToF_pressureY(F_pressure[1]);
+		AddToF_bodyX(F_body[0]);
+		AddToF_bodyY(F_body[1]);
+		AddToF_surf_tensionX(F_surf_tension[0]);
+		AddToF_surf_tensionY(F_surf_tension[1]);
+		AddToF_muX(F_mu[0]);
+		AddToF_muY(F_mu[1]);
+		AddToF_total_hydroX(F_total_hydro.x);
+		AddToF_total_hydroY(F_total_hydro.y);
+	#endif
 
-		// GRADIENTS AND NORMALS
-		gradPhi = calcGradPhi();
-		nx = gradPhi.x/gradPhi.z;
-		ny = gradPhi.y/gradPhi.z;
-
-		// CALCULATE FORCES:
-		real_t density_coeff = (Density_h-Density_l)/(PhaseField_h-PhaseField_l);
-		F_pressure[0] = (-1.0/3.0) * p * density_coeff * gradPhi.x;
-		F_pressure[1] = (-1.0/3.0) * p * density_coeff * gradPhi.y;
-		F_body[0] = -1.0*(rho-Density_l)*BuoyancyX + rho*GravitationX;
-		F_body[1] = -1.0*(rho-Density_l)*BuoyancyY + rho*GravitationY;
-
-		// Finding viscous force:
-		for (j=0;j<fixedIterator;j++) {
-			// GAMMA AND EQUILIBRIUM
-			mag = U*U + V*V;
-			for (i=0; i< 9; i++){
-				Gamma[i] = calcGamma(i, U, V, mag);
-				geq[i] = wf[i]*p + Gamma[i] - wf[i];
-			}
+	return F_total_hydro;
+}
 
-				// Stress/strain Tensor
-			stress[0] = 0.0;stress[1] = 0.0;stress[2] = 0.0;
-			for (i=0; i< 9; i++){
-				tmp1 = g[i] - geq[i];
-				stress[0] += tmp1*d2q9_ex[i]*d2q9_ex[i];
-				stress[1] += tmp1*d2q9_ex[i]*d2q9_ey[i];
-				stress[2] += tmp1*d2q9_ey[i]*d2q9_ey[i];    }
-
-			F_mu[0] = (0.5-tau)/tau * (stress[0]*gradPhi.x + stress[1]*gradPhi.y) * density_coeff;
-			F_mu[1] = (0.5-tau)/tau * (stress[1]*gradPhi.x + stress[2]*gradPhi.y) * density_coeff;
-			F_total[0] = mu*gradPhi.x + F_pressure[0] + F_body[0] + F_mu[0];
-			F_total[1] = mu*gradPhi.y + F_pressure[1] + F_body[1] + F_mu[1];
-			U = g[1]-g[3]+g[5]-g[6]-g[7]+g[8] + 0.5*F_total[0]/rho;
-			V = g[2]-g[4]+g[5]+g[6]-g[7]-g[8] + 0.5*F_total[1]/rho;
-		}
-		// PHASE FIELD POPULATION UPDATE:
-
-		real_t pfavg = 0.5*(PhaseField_l+PhaseField_h);
-		tmp1 = (1.0 - 4.0*(C - pfavg)*(C - pfavg))/W;
-		for (i=0; i< 9; i++){
-			F_phi[i] = calcF_phi(i, tmp1, nx, ny); // Forcing Terms
-			heq[i] = C * Gamma[i] - 0.5*F_phi[i];  // heq
-			h[i] = h[i] - omega_phi * (h[i]-heq[i]) + F_phi[i];	// collision
-		}
-		// PRESSURE EVOLUTION UPDATE:
-		for (i=0; i< 9; i++) {
-			F_i[i] = 3.0*wf[i] * (F_total[0]*d2q9_ex[i] + F_total[1]*d2q9_ey[i])/rho;
-			g[i] = g[i]-(g[i]-(geq[i]-0.5*F_i[i]))/tau+F_i[i];
-		}
-	}
+CudaDeviceFunction void CollisionCM(){
+	#ifdef OPTIONS_RT
+		real_t uold = U;
+		real_t vold = V;
+	#endif
+
+	real_t myPhaseF = PhaseF(0,0);
+
+	// Find Macroscopic Details
+	real_t rho = calcRho(myPhaseF);
+	real_t tau = calcTau(myPhaseF, rho);
+	real_t pfavg = 0.5*(PhaseField_l+PhaseField_h);
+
+	// Gradients & Normals
+	vector_t gradPhi = calcGradPhi(); // Phase field gradient
+	real_t p = getNormalizedPressure(); // normalized pressure
+
+	vector_t Fhydro = calcTotalHydrodynamicForceCM(gradPhi, myPhaseF, rho, tau,  calcMu(myPhaseF), p);
+
+	// PHASE FIELD POPULATION UPDATE:
+	vector_t u; u.x = U; u.y = V;
+	vector_t F_phi_vect = calcF_phi_xy(gradPhi,  myPhaseF,  pfavg);
+	#ifdef OPTIONS_debug
+		AddToF_phiX(F_phi_vect.x);
+		AddToF_phiY(F_phi_vect.y);
+	#endif
+	relax_and_collide_CM_phase_field(h, 1./omega_phi, u, F_phi_vect, rho);
+	
+	// PRESSURE EVOLUTION UPDATE:
+	relax_and_collide_CM_hydro(g, tau, u, Fhydro, rho);
+}
 #else
 CudaDeviceFunction void relaxMRT(real_t R[9], real_t tau)
 {
@@ -1099,8 +1002,8 @@ CudaDeviceFunction void CollisionMRT(){
 	real_t pfavg = 0.5*(PhaseField_l+PhaseField_h);
 
 	//  Gradients & Normals
-	vector_t gradPhi = calcGradPhi(); // Phase field gradient
-	real_t nx = gradPhi.x/gradPhi.z;  // GradPhi normalized in x, y direction
+	vector_t gradPhi = calcGradPhi(); // Phase field gradient, .z stores magnitude
+	real_t nx = gradPhi.x/gradPhi.z;
 	real_t ny = gradPhi.y/gradPhi.z;
 
 	real_t p = getNormalizedPressure(); // normalized pressure
@@ -1120,41 +1023,35 @@ CudaDeviceFunction void CollisionMRT(){
 
 	real_t tmp1 = (1.0 - 4.0*(myPhaseF - pfavg)*(myPhaseF - pfavg))/W;
 	for (int i=0; i< 9; i++){
-		#ifdef OPTIONS_RT
-			Req[i] = 3.0 * wf[i] * (d2q9_ex[i] * (myPhaseF*U - PhaseOld*uold)
-				+ d2q9_ey[i] * (myPhaseF*V - PhaseOld*vold));
-			F_phi[i] = calcF_phi(i, tmp1, nx, ny) + Req[i];
-			heq[i] = myPhaseF * wf[i] * (1 + 3.0 * (d2q9_ex[i]*U + d2q9_ey[i]*V) ) - 0.5*F_phi[i];
-		#else
-			F_phi[i] = calcF_phi(i, tmp1, nx, ny); // Forcing Terms:
-			//calcF_phi: F_phi[i] = wf[i] * tmp1 * (d2q9_ex[i]*nx + d2q9_ey[i]*ny);
-			heq[i] = myPhaseF * Gamma[i] - 0.5*F_phi[i];  // heq
-		#endif
-		#ifdef OPTIONS_debug
-			AddToF_phiX(F_phi[i]*d2q9_ex[i]);
-			AddToF_phiY(F_phi[i]*d2q9_ey[i]);
-		#endif
-			h[i] = h[i] - omega_phi * (h[i]-heq[i]) + F_phi[i];	// collision
+			#ifdef OPTIONS_RT
+				Req[i] = 3.0 * wf[i] * (d2q9_ex[i] * (myPhaseF*U - PhaseOld*uold) + d2q9_ey[i] * (myPhaseF*V - PhaseOld*vold));
+				F_phi[i] = calcF_phi(i, tmp1, nx, ny) + Req[i];
+				heq[i] = myPhaseF * wf[i] * (1 + 3.0 * (d2q9_ex[i]*U + d2q9_ey[i]*V) ) - 0.5*F_phi[i];
+			#else
+				F_phi[i] = calcF_phi(i, tmp1, nx, ny); // Forcing Terms: F_phi[i] = wf[i] * tmp1 * (d2q9_ex[i]*nx + d2q9_ey[i]*ny);
+				heq[i] = myPhaseF * Gamma[i] - 0.5*F_phi[i];  // heq
+			#endif
+			#ifdef OPTIONS_debug
+				AddToF_phiX(F_phi[i]*d2q9_ex[i]);
+				AddToF_phiY(F_phi[i]*d2q9_ey[i]);
+			#endif
+			h[i] = h[i] - omega_phi * (h[i]-heq[i]) + F_phi[i];
 		}
 
 	//  PRESSURE EVOLUTION UPDATE:
 	// ---------------------- MRT ---------------------------
 	#ifdef OPTIONS_GF
-		real_t Falpha_x, Falpha_y; // FOR EXTENDED FORCE SCHEME
+		real_t Falpha_x, Falpha_y;
 	#endif
 	for (int i=0; i< 9; i++){
-	#ifdef OPTIONS_GF
-		Falpha_x = d2q9_ex[i] - U + 3.0 * d2q9_ex[i] * (d2q9_ex[i]*U + d2q9_ey[i]*V); //FOR EXTENDED FORCE SCHEME
-		Falpha_y = d2q9_ey[i] - V + 3.0 * d2q9_ey[i] * (d2q9_ex[i]*U + d2q9_ey[i]*V); //FOR EXTENDED FORCE SCHEME
-		F_i[i] = 3.0*wf[i] * (Falpha_x * Fhydro.x + Falpha_y * Fhydro.y )/rho;        //FOR EXTENDED FORCE SCHEME
+	#ifdef OPTIONS_GF 
+		Falpha_x = d2q9_ex[i] - U + 3.0 * d2q9_ex[i] * (d2q9_ex[i]*U + d2q9_ey[i]*V);
+		Falpha_y = d2q9_ey[i] - V + 3.0 * d2q9_ey[i] * (d2q9_ex[i]*U + d2q9_ey[i]*V);
+		F_i[i] = 3.0*wf[i] * (Falpha_x * Fhydro.x + Falpha_y * Fhydro.y )/rho;        
 	#else
-		// eq 15 from "Improved locality of the phase-field lattice Boltzmann 
-		//	model for immiscible fluids at high density ratios"
-		F_i[i] = 3.0*wf[i] * (Fhydro.x*d2q9_ex[i] + Fhydro.y*d2q9_ey[i])/rho;
+		F_i[i] = 3.0*wf[i] * (Fhydro.x*d2q9_ex[i] + Fhydro.y*d2q9_ey[i])/rho; // eq 15 from "Improved locality of the phase-field lattice Boltzmann model for immiscible fluids at high density ratios"
 	#endif
-		// eq 16 from "Improved locality of the phase-field lattice Boltzmann 
-		// model for immiscible fluids at high density ratios"
-		R[i] = g[i] - (geq[i] - 0.5*F_i[i]);
+		R[i] = g[i] - (geq[i] - 0.5*F_i[i]); // eq 16 from "Improved locality of the phase-field lattice Boltzmann model for immiscible fluids at high density ratios"
 	}
 
 	relaxMRT(R, tau);
diff --git a/models/multiphase/d2q9_pf_velocity/conf.mk b/models/multiphase/d2q9_pf_velocity/conf.mk
index 10d4bb3b5..d4cb8d1ca 100755
--- a/models/multiphase/d2q9_pf_velocity/conf.mk
+++ b/models/multiphase/d2q9_pf_velocity/conf.mk
@@ -1,6 +1,6 @@
 ADJOINT=0
 TEST=FALSE
-OPT="(GF+RT+Outflow+GuoCM+debug+BGK+CM)*autosym"
+OPT="(GF+RT+Outflow+debug+CM)*autosym"
 
 # GF: Guo Forcing - feature in MRT model;
 #	This is using a higher order Forcing scheme
@@ -10,7 +10,6 @@ OPT="(GF+RT+Outflow+GuoCM+debug+BGK+CM)*autosym"
 #	phase field equilibrium distribution function by
 #	Ren et al. (2016)
 # debug: Enables tracking of momentum and force globals 
-# BGK: Applies single relaxation time - not recommended
 # CM: Applies the Central moments relaxation
 
 # Boundary Conditions: