diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp
index 37409c57..f2809d07 100644
--- a/common/ScaLBL.cpp
+++ b/common/ScaLBL.cpp
@@ -4159,7 +4159,7 @@ void ScaLBL_Communicator::Color_BC_z(int *Map, double *Phi, double *Den, double
 	if (kproc == 0) {
 		// Set the phase indicator field and density on the z inlet
 		ScaLBL_Color_BC_z(dvcSendList_z, Map, Phi, Den, vA, vB, sendCount_z, N);
-		ScaLBL_SetSlice_z(Phi,Value,Nx,Ny,Nz,0);
+		//ScaLBL_SetSlice_z(Phi,Value,Nx,Ny,Nz,0);
 	}
 }
 
@@ -4168,7 +4168,7 @@ void ScaLBL_Communicator::Color_BC_Z(int *Map, double *Phi, double *Den, double
 	if (kproc == nprocz-1){
 		// Set the phase indicator field and density on the Z outlet
 		ScaLBL_Color_BC_Z(dvcSendList_Z, Map, Phi, Den, vA, vB, sendCount_Z, N);
-		ScaLBL_SetSlice_z(Phi,Value,Nx,Ny,Nz,Nz-1);
+		//ScaLBL_SetSlice_z(Phi,Value,Nx,Ny,Nz,Nz-1);
 	}
 }
 
diff --git a/gpu/Color.cu b/gpu/Color.cu
index 47251a43..6bd76534 100644
--- a/gpu/Color.cu
+++ b/gpu/Color.cu
@@ -136,8 +136,10 @@ __global__  void dvc_ScaLBL_Color_BC(int *list, int *Map, double *Phi, double *D
 		Den[n] = vA;
 		Den[Np+n] = vB;
 		
-		nm = Map[n];
-		Phi[nm] = (vA-vB)/(vA+vB);
+		//nm = Map[n];
+		//Phi[nm] = (vA-vB)/(vA+vB);
+		Phi[n] = (vA-vB)/(vA+vB);
+
 	}
 }
 __global__  void dvc_ScaLBL_Color_BC_z(int *list, int *Map, double *Phi, double *Den, double vA, double vB, int count, int Np)
@@ -151,8 +153,9 @@ __global__  void dvc_ScaLBL_Color_BC_z(int *list, int *Map, double *Phi, double
 		Den[Np+n] = vB;
 		//double valB = Den[Np+n]; // mass that reaches inlet is conserved
 
-		nm = Map[n];
-		Phi[nm] = (vA-vB)/(vA+vB);
+		//nm = Map[n];
+		//Phi[nm] = (vA-vB)/(vA+vB);
+		Phi[n] = (vA-vB)/(vA+vB);
 	}
 }
 
@@ -167,8 +170,10 @@ __global__  void dvc_ScaLBL_Color_BC_Z(int *list, int *Map, double *Phi, double
 		Den[n] = vA;
 		Den[Np+n] = vB;
 		
-		nm = Map[n];
-		Phi[nm] = (vA-vB)/(vA+vB);
+		//nm = Map[n];
+		//Phi[nm] = (vA-vB)/(vA+vB);
+		Phi[n] = (vA-vB)/(vA+vB);
+
 	}
 }
 //*************************************************************************
diff --git a/tests/lbpm_dfh_simulator.cpp b/tests/lbpm_dfh_simulator.cpp
index d5568b13..a62da5fc 100644
--- a/tests/lbpm_dfh_simulator.cpp
+++ b/tests/lbpm_dfh_simulator.cpp
@@ -492,7 +492,7 @@ int main(int argc, char **argv)
 		ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
 		ScaLBL_DeviceBarrier();
 		delete [] TmpMap;
-		
+
 		// Compute the solid interaction potential and copy result to device
 		if (rank==0) printf("Computing solid interaction potential \n");
 		double *Tmp;
@@ -523,14 +523,14 @@ int main(int argc, char **argv)
 						for (int kk=0; kk<5; kk++){
 							for (int jj=0; jj<5; jj++){
 								for (int ii=0; ii<5; ii++){
-									
+
 									int index = kk*25+jj*5+ii;
 									double distval= Dst[index];
 
 									int idi=i+ii-2;
 									int idj=j+jj-2;
 									int idk=k+kk-2;
-									
+
 									if (idi < 0) idi=0;
 									if (idj < 0) idj=0;
 									if (idk < 0) idk=0;
@@ -540,16 +540,16 @@ int main(int argc, char **argv)
 
 									int nn = idk*Nx*Ny + idj*Nx + idi;
 									if (!(Mask.id[nn] > 0)){
-									double vec_x = double(ii-2);
-									double vec_y = double(jj-2);
-									double vec_z = double(kk-2);
-									
-									double ALPHA=PhaseLabel[nn];
-									double GAMMA=-2.f;
-									if (distval > 2.f) ALPHA=0.f; // symmetric cutoff distance									
-									phi_x += ALPHA*exp(GAMMA*distval)*vec_x/distval;
-									phi_y += ALPHA*exp(GAMMA*distval)*vec_y/distval;
-									phi_z += ALPHA*exp(GAMMA*distval)*vec_z/distval;
+										double vec_x = double(ii-2);
+										double vec_y = double(jj-2);
+										double vec_z = double(kk-2);
+
+										double ALPHA=PhaseLabel[nn];
+										double GAMMA=-2.f;
+										if (distval > 2.f) ALPHA=0.f; // symmetric cutoff distance									
+										phi_x += ALPHA*exp(GAMMA*distval)*vec_x/distval;
+										phi_y += ALPHA*exp(GAMMA*distval)*vec_y/distval;
+										phi_z += ALPHA*exp(GAMMA*distval)*vec_z/distval;
 									}
 								}
 							}
@@ -557,17 +557,17 @@ int main(int argc, char **argv)
 						Tmp[idx] = phi_x;
 						Tmp[idx+Np] = phi_y;
 						Tmp[idx+2*Np] = phi_z;
-						
+
 						/*						double d = Averages->SDs(n);
 												double dx = Averages->SDs_x(n);
 												double dy = Averages->SDs_y(n);
 												double dz = Averages->SDs_z(n);
 												double value=cns*exp(-bns*fabs(d))-cws*exp(-bns*fabs(d));
-												
+
 						Tmp[idx] = value*dx;
 						Tmp[idx+Np] = value*dy;
 						Tmp[idx+2*Np] = value*dz;
-						*/
+						 */
 					}
 				}
 			}
@@ -642,7 +642,7 @@ int main(int argc, char **argv)
 		// Copy the phase from the GPU -> CPU
 		//...........................................................................
 		ScaLBL_DeviceBarrier();
-		ScaLBL_CopyToHost(Averages->Phase.data(),Phi,Np*sizeof(double));
+		ScaLBL_Comm.RegularLayout(Map,Phi,Averages->Phase);
 		ScaLBL_Comm.RegularLayout(Map,Pressure,Averages->Press);
 		ScaLBL_Comm.RegularLayout(Map,&Velocity[0],Averages->Vel_x);
 		ScaLBL_Comm.RegularLayout(Map,&Velocity[Np],Averages->Vel_y);