From fc365e9aba6837f138981f0b16e33bb710f65271 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 1 May 2018 11:13:26 -0400
Subject: [PATCH] DFH gradient works (no comm yet)

---
 cpu/dfh.cpp                | 11 +++++------
 tests/TestColorGradDFH.cpp | 39 ++++++++++++++++++++++++++++++++------
 2 files changed, 38 insertions(+), 12 deletions(-)

diff --git a/cpu/dfh.cpp b/cpu/dfh.cpp
index 960a51ab..9030aa5f 100644
--- a/cpu/dfh.cpp
+++ b/cpu/dfh.cpp
@@ -1429,7 +1429,6 @@ extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *neighborList, double *Phi, double
 		m1 = Phi[nn];
 		nn = neighborList[n]%Np;
 		m2 = Phi[nn];
-		printf("m1=%f, m2=%f \n",m1,m2);
 		nn = neighborList[n+3*Np]%Np;
 		m3 = Phi[nn];
 		nn = neighborList[n+2*Np]%Np;
@@ -1438,7 +1437,7 @@ extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *neighborList, double *Phi, double
 		m5 = Phi[nn];
 		nn = neighborList[n+4*Np]%Np;
 		m6 = Phi[nn];		
-/*		nn = neighborList[n+7*Np]%Np;
+		nn = neighborList[n+7*Np]%Np;
 		m7 = Phi[nn];
 		nn = neighborList[n+6*Np]%Np;
 		m8 = Phi[nn];		
@@ -1462,12 +1461,12 @@ extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *neighborList, double *Phi, double
 		m17 = Phi[nn];
 		nn = neighborList[n+16*Np]%Np;
 		m18 = Phi[nn];					
-		*/
+		
 		//............Compute the Color Gradient...................................
 		//............Compute the wn fluid Gradient...................................
-		nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
-		ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
-		nz = -(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+		nx = (m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+		ny = (m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+		nz = (m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
 		
 		//...........Normalize the Color Gradient.................................
 		//	C = sqrt(nx*nx+ny*ny+nz*nz);
diff --git a/tests/TestColorGradDFH.cpp b/tests/TestColorGradDFH.cpp
index e8ac71b8..73408382 100644
--- a/tests/TestColorGradDFH.cpp
+++ b/tests/TestColorGradDFH.cpp
@@ -158,7 +158,6 @@ int main(int argc, char **argv)
 		Nz += 2;
 		int N = Nx*Ny*Nz;
 
-		int Np=0;  // number of local pore nodes
 		double *PhaseLabel;
 		PhaseLabel = new double[N];
 		//.......................................................................
@@ -167,7 +166,6 @@ int main(int argc, char **argv)
 				for (i=0;i<Nx;i++){
 					n = k*Nx*Ny+j*Nx+i;
 					Dm.id[n]=1;
-					Np++;
 					// Initialize gradient ColorGrad = (1,2,3)
 					double value=double(3*k+2*j+i);
 					PhaseLabel[n]= value;
@@ -185,13 +183,23 @@ int main(int argc, char **argv)
 
 		// LBM variables
 		if (rank==0)	printf ("Set up the neighborlist \n");
+		
+		int Np=0;  // number of local pore nodes
+		for (k=1;k<Nz-1;k++){
+			for (j=1;j<Ny-1;j++){
+				for (i=1;i<Nx-1;i++){
+					Np++;
+				}
+			}
+		}
+		Np+=32;
 
 		int neighborSize=18*Np*sizeof(int);
 		int *neighborList;
 		IntArray Map(Nx,Ny,Nz);
 		neighborList= new int[18*Np];
 
-		ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np);
+		Np = ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np);
 		MPI_Barrier(comm);
 
 		//......................device distributions.................................
@@ -252,7 +260,20 @@ int main(int argc, char **argv)
     	int SIZE=3*Np*sizeof(double);
     	ScaLBL_CopyToHost(&COLORGRAD[0],&ColorGrad[0],SIZE);
 
-
+    	for (k=1;k<Nz-1;k++){
+    		for (j=1;j<Ny-1;j++){
+    			for (i=1;i<Nx-1;i++){
+    				n = k*Nx*Ny+j*Nx+i;
+    				if (Dm.id[n] > 0){
+    					int idx = Map(i,j,k);
+    					printf("%i ",idx);
+    				}
+    			}
+    			printf("\n");
+    		}
+    		printf("-------\n");
+    	}
+    	
     	double CX,CY,CZ;
     	for (k=2;k<Nz-2;k++){
     		for (j=2;j<Ny-2;j++){
@@ -264,8 +285,14 @@ int main(int argc, char **argv)
     					CY=COLORGRAD[Np+idx];
     					CZ=COLORGRAD[2*Np+idx];
     					double error=sqrt((CX-1.0)*(CX-1.0)+(CY-2.0)*(CY-2.0)+ (CZ-3.0)*(CZ-3.0));
-    					if (error > 1e-8)
-    						printf("i,j,k=%i,%i,%i: Color gradient=%f,%f,%f \n",i,j,k,CX,CY,CZ);
+    					if (error > 1e-8){
+    						printf("i,j,k=%i,%i,%i; idx=%i: Color gradient=%f,%f,%f \n",i,j,k,idx,CX,CY,CZ);
+    						for (int q=0; q<18; q++){
+    							int nn = neighborList[q*Np+idx]%Np;
+    							double value= PHASE[nn];
+    							printf("     q=%i, nn=%i, value=%f \n",q,nn,value);
+    						}
+    					}
     				}
     			}
     		}