From 99f1d9b72741aac7776349ecb9515b9ae305c73c Mon Sep 17 00:00:00 2001
From: Rex Zhe Li <zhe.rex.li@gmail.com>
Date: Thu, 25 Feb 2021 20:30:48 -0500
Subject: [PATCH] add a test routine for mixed gradient

---
 common/ScaLBL.h         |   1 +
 cpu/FreeLee.cpp         | 164 ++++++++++++++++++++++++++++++++++++++++
 models/FreeLeeModel.cpp | 148 ++++++++++++++++++++++++++++++++++++
 models/FreeLeeModel.h   |   3 +
 tests/CMakeLists.txt    |   1 +
 tests/TestMixedGrad.cpp |  85 +++++++++++++++++++++
 6 files changed, 402 insertions(+)
 create mode 100644 tests/TestMixedGrad.cpp

diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index 42c51525..a6a4aeae 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -207,6 +207,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborLis
 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);
 
+extern "C" void ScaLBL_D3Q9_MGTest(int *Map, double *Phi,double *ColorGrad, int start, int finish, Np);
 
 // BOUNDARY CONDITION ROUTINES
 
diff --git a/cpu/FreeLee.cpp b/cpu/FreeLee.cpp
index c77436a2..8d5eeb0d 100644
--- a/cpu/FreeLee.cpp
+++ b/cpu/FreeLee.cpp
@@ -1916,3 +1916,167 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, d
         Pressure[n] = p;
 	}
 }
+
+extern "C" void ScaLBL_D3Q9_MGTest(int *Map, double *Phi,double *ColorGrad, int start, int finish, Np){
+
+	int n,nn,nn2x,ijk;
+	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 nx,ny,nz;//normal color gradient
+    double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
+    double phi;
+
+	for (int n=start; n<finish; n++){
+
+		// Get the 1D index based on regular data layout
+		ijk = Map[n];
+        phi = Phi[ijk];// load phase field
+		//					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 Mixed Gradient...................................
+		mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
+		mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
+		mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
+
+		ColorGrad[0*Np+n] = mgx;
+		ColorGrad[1*Np+n] = mgy;
+		ColorGrad[2*Np+n] = mgz;
+    }
+}
diff --git a/models/FreeLeeModel.cpp b/models/FreeLeeModel.cpp
index 0b7f5e90..588dfb14 100644
--- a/models/FreeLeeModel.cpp
+++ b/models/FreeLeeModel.cpp
@@ -1059,3 +1059,151 @@ void ScaLBL_FreeLeeModel::WriteDebug_SingleFluid(){
 	fwrite(PhaseField.data(),8,N,VELZ_FILE);
 	fclose(VELZ_FILE);
 }
+
+void ScaLBL_FreeLeeModel::Create_DummyPhase_MGTest(){
+	// Initialize communication structures in averaging domain
+	for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
+	Mask->CommInit();
+	Np=Mask->PoreCount();
+	//...........................................................................
+	if (rank==0)    printf ("Create ScaLBL_Communicator \n");
+	// Create a communicator for the device (will use optimized layout)
+	// ScaLBL_Communicator ScaLBL_Comm(Mask); // original
+	//ScaLBL_Comm  = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
+	//ScaLBL_Comm_Regular  = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
+	ScaLBL_Comm_WideHalo  = std::shared_ptr<ScaLBLWideHalo_Communicator>(new ScaLBLWideHalo_Communicator(Mask,2));
+
+	// create the layout for the LBM
+	int Npad=(Np/16 + 2)*16;
+	if (rank==0)    printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N);
+	Map.resize(Nx,Ny,Nz);       Map.fill(-2);
+	auto neighborList= new int[18*Npad];
+	Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,2);
+	comm.barrier();
+
+	//...........................................................................
+	//                MAIN  VARIABLES ALLOCATED HERE
+	//...........................................................................
+	// LBM variables
+	if (rank==0)    printf ("Allocating distributions \n");
+	//......................device distributions.................................
+	dist_mem_size = Np*sizeof(double);
+	neighborSize=18*(Np*sizeof(int));
+	//...........................................................................
+	//ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
+	ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np);
+	//ScaLBL_AllocateDeviceMemory((void **) &gqbar, 19*dist_mem_size);
+	//ScaLBL_AllocateDeviceMemory((void **) &hq, 7*dist_mem_size);
+	//ScaLBL_AllocateDeviceMemory((void **) &mu_phi, dist_mem_size);
+	//ScaLBL_AllocateDeviceMemory((void **) &Den, dist_mem_size);
+	ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Nh);		
+	//ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np);
+	//ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np);
+	//...........................................................................
+	// Update GPU data structures
+	if (rank==0)	printf ("Setting up device map and neighbor list \n");
+	fflush(stdout);
+	int *TmpMap;
+	TmpMap=new int[Np];
+	for (int k=1; k<Nz-1; k++){
+		for (int j=1; j<Ny-1; j++){
+			for (int i=1; i<Nx-1; i++){
+				int idx=Map(i,j,k);
+				if (!(idx < 0))
+					TmpMap[idx] = ScaLBL_Comm_WideHalo->Map(i,j,k);
+			}
+		}
+	}
+	// check that TmpMap is valid
+	for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
+		auto n = TmpMap[idx];
+		if (n > Nxh*Nyh*Nzh){
+			printf("Bad value! idx=%i \n", n);
+			TmpMap[idx] = Nxh*Nyh*Nzh-1;
+		}
+	}
+	for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
+		auto n = TmpMap[idx];
+		if ( n > Nxh*Nyh*Nzh ){
+			printf("Bad value! idx=%i \n",n);
+			TmpMap[idx] = Nxh*Nyh*Nzh-1;
+		}
+	}
+    // copy the device map
+	ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
+	// copy the neighbor list 
+	//ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
+	comm.barrier();
+
+	double *phase;
+	phase = new double[Nh];
+
+	for (int k=0;k<Nzh;k++){
+		for (int j=0;j<Nyh;j++){
+			for (int i=0;i<Nxh;i++){
+
+                //idx for double-halo array 'phase'
+				int nh = k*Nxh*Nyh+j*Nxh+i;
+
+                //idx for single-halo array Mask->id[n]
+                int x=i-1;
+                int y=j-1;
+                int z=k-1;
+				if (x<0)   x=0;
+				if (y<0)   y=0;
+				if (z<0)   z=0;
+				if (x>=Nx) x=Nx-1;
+				if (y>=Ny) y=Ny-1;
+				if (z>=Nz) z=Nz-1;
+				int n = z*Nx*Ny+y*Nx+x;
+				phase[nh]=id[n];
+			}
+		}
+	}
+	ScaLBL_CopyToDevice(Phi, phase, Nh*sizeof(double));
+	ScaLBL_Comm->Barrier();
+	comm.barrier();
+	delete [] TmpMap;
+	delete [] neighborList;
+    delete [] phase;
+}
+
+void ScaLBL_FreeLeeModel::MGTest(){
+
+	comm.barrier();
+
+	ScaLBL_Comm_WideHalo->Send(Phi);
+    ScaLBL_D3Q9_MGTest(dvcMap,Phi,ColorGrad, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+	ScaLBL_Comm_WideHalo->Send(Phi);
+    ScaLBL_D3Q9_MGTest(dvcMap,Phi,ColorGrad, 0, ScaLBL_Comm->LastExterior(), Np);
+
+    //check the sum of ColorGrad
+    double cgx_loc = 0.0;
+    double cgy_loc = 0.0;
+    double cgz_loc = 0.0;
+    double cgx,cgy,cgz;
+    double *ColorGrad_host;
+	ColorGrad_host = new double [3*Np];
+	ScaLBL_CopyToHost(&ColorGrad_host[0],&ColorGrad[0], 3*Np*sizeof(double));
+    for (int i = ScaLBL_Comm->FirstInterior(), i<ScaLBL_Comm->LastInterior(),i++){
+        cgx_loc+=ColorGrad_host[0*Np+i];
+        cgy_loc+=ColorGrad_host[1*Np+i];
+        cgz_loc+=ColorGrad_host[2*Np+i];
+    }
+    for (int i = 0, i<ScaLBL_Comm->LastExterior(),i++){
+        cgx_loc+=ColorGrad_host[0*Np+i];
+        cgy_loc+=ColorGrad_host[1*Np+i];
+        cgz_loc+=ColorGrad_host[2*Np+i];
+    }
+    cgx=Dm->Comm.sumReduce( cgx_loc);
+    cgy=Dm->Comm.sumReduce( cgy_loc);
+    cgz=Dm->Comm.sumReduce( cgz_loc);
+    if (rank==0){
+        printf("Sum of all x-component of the mixed gradient = %.2g",cgx);
+        printf("Sum of all y-component of the mixed gradient = %.2g",cgy);
+        printf("Sum of all z-component of the mixed gradient = %.2g",cgz);
+    }
+
+    delete [] ColorGrad_host;
+}
diff --git a/models/FreeLeeModel.h b/models/FreeLeeModel.h
index 1e372f50..23afe39b 100644
--- a/models/FreeLeeModel.h
+++ b/models/FreeLeeModel.h
@@ -34,6 +34,9 @@ public:
 	void Initialize_SingleFluid();
 	void Run_SingleFluid();
 	void WriteDebug_SingleFluid();
+    // test utilities
+    void Create_DummyPhase_MGTest();
+    void MGTest();
 	
 	bool Restart,pBC;
 	int timestep,timestepMax;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 8df4e6bd..b32c0b57 100755
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -43,6 +43,7 @@ ADD_LBPM_EXECUTABLE( TestPoissonSolver )
 ADD_LBPM_EXECUTABLE( TestIonModel )
 ADD_LBPM_EXECUTABLE( TestNernstPlanck )
 ADD_LBPM_EXECUTABLE( TestPNP_Stokes )
+ADD_LBPM_EXECUTABLE( TestMixedGrad )
 
 
 
diff --git a/tests/TestMixedGrad.cpp b/tests/TestMixedGrad.cpp
new file mode 100644
index 00000000..821d6465
--- /dev/null
+++ b/tests/TestMixedGrad.cpp
@@ -0,0 +1,85 @@
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <stdexcept>
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/stat.h>
+
+#include "common/Utilities.h"
+#include "models/FreeLeeModel.h"
+
+inline void Initialize_DummyPhaseField(ScaLBL_FreeLeeModel &LeeModel){
+	// initialize a bubble
+	int i,j,k,n;
+	int rank = LeeModel.Mask->rank();
+	int Nx = LeeModel.Mask->Nx;
+	int Ny = LeeModel.Mask->Ny;
+	int Nz = LeeModel.Mask->Nz;
+	if (rank == 0) cout << "Setting up dummy phase field..." << endl;
+	for (k=0;k<Nz;k++){
+		for (j=0;j<Ny;j++){
+			for (i=0;i<Nx;i++){
+				n = k*Nx*Ny + j*Nz + i;
+                LeeModel.Mask->id[n]=1;
+                LeeMOdel.id[n] = LeeModel.Mask->id[n];
+			}
+		}
+	}
+}
+
+
+int main( int argc, char **argv )
+{
+
+    // Initialize
+    Utilities::startup( argc, argv );
+
+    // Load the input database
+    auto db = std::make_shared<Database>( argv[1] );
+
+    { // Limit scope so variables that contain communicators will free before MPI_Finialize
+
+        Utilities::MPI comm( MPI_COMM_WORLD );
+        int rank   = comm.getRank();
+        int nprocs = comm.getSize();
+
+        if ( rank == 0 ) {
+            printf( "********************************************************\n" );
+            printf( "Running Mixed Gradient Test	\n" );
+            printf( "********************************************************\n" );
+        }
+        // Initialize compute device
+        int device = ScaLBL_SetDevice( rank );
+        NULL_USE( device );
+        ScaLBL_DeviceBarrier();
+        comm.barrier();
+
+        PROFILE_ENABLE( 1 );
+        // PROFILE_ENABLE_TRACE();
+        // PROFILE_ENABLE_MEMORY();
+        PROFILE_SYNCHRONIZE();
+        PROFILE_START( "Main" );
+        Utilities::setErrorHandlers();
+
+        auto filename = argv[1];
+        ScaLBL_FreeLeeModel LeeModel( rank, nprocs, comm );
+        LeeModel.ReadParams( filename );
+        LeeModel.SetDomain();
+        Initialize_DummyPhaseField(LeeModel);
+        LeeModel.Create_DummyPhase_MGTest();     
+        LeeModel.MGTest()
+
+        PROFILE_STOP( "Main" );
+        auto file  = db->getWithDefault<std::string>( "TimerFile", "TestMixedGrad" );
+        auto level = db->getWithDefault<int>( "TimerLevel", 1 );
+        PROFILE_SAVE( file, level );
+        // ****************************************************
+
+
+    } // Limit scope so variables that contain communicators will free before MPI_Finialize
+
+    Utilities::shutdown();
+    return 0;
+
+}