add membrane unit test

analysis/Membrane.cpp

@@ -343,13 +343,23 @@ Membrane::~Membrane() {
 	ScaLBL_FreeDeviceMemory( dvcRecvDist_YZ );
 }
 
-int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int rank_p, int rank_q
+int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int rank_q, int rank_Q,
-		const int *originalSendList, const int shift_x, const int shift_y, const int shift_z, 
+		const int shift_x, const int shift_y, const int shift_z, 
+		const int *originalSendList, const int sendCount, const int recvCount,
 		const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList){
 
+	int sendtag = 2389;
+	int recvtag = 2389;
 	int i,j,k,n,nn,idx;
 	double dist,locdist;
 	
+	int *SendList;
+	int *RecvList;
+	int *RecvList_q;
+	SendList = new int [sendCount]; // for this process
+	RecvList = new int [recvCount];   // filled in by the neighbor process
+	RecvList_q = new int [sendCount]; // goes to the neighbor process
+
 	int regularCount = 0;
 	int membraneCount = 0;
 	for (idx=0; idx<sendCount; idx++){
@@ -377,16 +387,27 @@ int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int
 		if (dist*locdist < 0.0){
 			/* store membrane values at the end */
 			membraneCount++;
-			membraneRecvList[sendCount-membraneCount] = nn;
+			RecvList_q[sendCount-membraneCount] = nn;
-			membraneSendList[sendCount-membraneCount] = n;
+			SendList[sendCount-membraneCount] = n;
 		}
 		else {
-			membraneRecvList[regularCount] = nn;
+			RecvList_q[regularCount] = nn;
-			membraneSendList[regularCount++] = n;
+			SendList[regularCount++] = n;
 		}
 	}
+	// send the recv list to the neighbor process
+    req1[0] = Dm->Comm.Isend(RecvList_q, sendCount, rank_q, sendtag);
+    req2[0] = Dm->Comm.Irecv(RecvList, recvCount, rank_Q, recvtag);
+    Dm->Barrier();
 	
+	// Return updated version to the device
+	ScaLBL_CopyToDevice(&membraneSendList[start], SendList, sendCount*sizeof(int));
+	ScaLBL_CopyToDevice(&membraneRecvList[start], RecvList, recvCount*sizeof(int));
 
+	// clean up the work arrays
+	delete [] SendList;
+	delete [] RecvList;
+	delete [] RecvList_q;
 	return membraneCount;
 }
 

analysis/Membrane.h

@@ -63,8 +63,10 @@ private:
     * @param recvDistance   - distance values from neighboring processor
     * @param d3q19_recvlist - device array with the saved list
     * */
-    int D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, const int *list,  int start, int count,
+    int D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int rank_q, int rank_Q,
-    		DoubleArray &Distance, double *recvDistance, int *d3q19_recvlist);
+    		const int shift_x, const int shift_y, const int shift_z, 
+    		const int *originalSendList, const int sendCount, const int recvCount,
+    		const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList);
 	//......................................................................................
 	// MPI ranks for all 18 neighbors
 	//......................................................................................

tests/CMakeLists.txt

@@ -59,6 +59,7 @@ ADD_LBPM_TEST( TestTopo3D )
 ADD_LBPM_TEST( TestFluxBC )
 ADD_LBPM_TEST( TestFlowAdaptor )
 ADD_LBPM_TEST( TestMap )
+ADD_LBPM_TEST( TestMembrane )
 #ADD_LBPM_TEST( TestMRT )
 #ADD_LBPM_TEST( TestColorGrad )
 ADD_LBPM_TEST( TestWideHalo )

tests/TestMembrane.cpp

@@ -0,0 +1,198 @@
+
+//*************************************************************************
+// Lattice Boltzmann Simulator for Single Phase Flow in Porous Media
+// James E. McCLure
+//*************************************************************************
+#include <stdio.h>
+#include <iostream>
+#include <fstream>
+#include "common/ScaLBL.h"
+#include "common/MPI.h"
+
+using namespace std;
+
+std::shared_ptr<Database> loadInputs( int nprocs )
+{
+  //auto db = std::make_shared<Database>( "Domain.in" );
+    auto db = std::make_shared<Database>();
+    db->putScalar<int>( "BC", 0 );
+    db->putVector<int>( "nproc", { 1, 1, 1 } );
+    db->putVector<int>( "n", { 5, 5, 5 } );
+    db->putScalar<int>( "nspheres", 1 );
+    db->putVector<double>( "L", { 1, 1, 1 } );
+    return db;
+}
+
+//***************************************************************************************
+int main(int argc, char **argv)
+{
+	// Initialize MPI
+    Utilities::startup( argc, argv );
+	Utilities::MPI comm( MPI_COMM_WORLD );
+	int check=0;
+	{
+
+		int i,j,k,n;
+		
+		static int D3Q19[18][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1},
+				{1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0},
+				{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1},
+				{0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}};
+
+        int rank = comm.getRank();
+		if (rank == 0){
+			printf("********************************************************\n");
+			printf("Running unit test: TestMap	\n");
+			printf("********************************************************\n");
+		}
+		
+	    // Load inputs
+	    auto db = loadInputs( comm.getSize() );
+	    int Nx = db->getVector<int>( "n" )[0];
+	    int Ny = db->getVector<int>( "n" )[1];
+	    int Nz = db->getVector<int>( "n" )[2];
+		auto Dm = std::make_shared<Domain>(db,comm);
+
+		Nx += 2;
+		Ny += 2;
+		Nz += 2;
+		int N = Nx*Ny*Nz;
+		//.......................................................................
+		int Np = 0;
+		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;
+					Dm->id[n] = 1;
+					Np++;
+				}
+			}
+		}
+		Dm->CommInit();
+
+		// Create a communicator for the device (will use optimized layout)
+		std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm(new ScaLBL_Communicator(Dm));
+		//Create a second communicator based on the regular data layout
+		std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular(new ScaLBL_Communicator(Dm));
+
+		if (rank==0){
+			printf("Total domain size = %i \n",N);
+			printf("Reduced domain size = %i \n",Np);
+		}
+
+		// LBM variables
+		if (rank==0)	printf ("Set up the neighborlist \n");
+		int Npad=Np+32;
+		int neighborSize=18*Npad*sizeof(int);
+		int *neighborList;
+		IntArray Map(Nx,Ny,Nz);
+		neighborList= new int[18*Npad];
+
+		Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,1);
+		comm.barrier();
+		
+		// Check the neighborlist
+		printf("Check neighborlist: exterior %i, first interior %i last interior %i \n",ScaLBL_Comm->LastExterior(),ScaLBL_Comm->FirstInterior(),ScaLBL_Comm->LastInterior());
+		for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
+			for (int q=0; q<18; q++){
+				int nn = neighborList[q*Np+idx]%Np;
+				if (nn>Np) printf("neighborlist error (exterior) at q=%i, idx=%i \n",q,idx);
+		      
+			}
+		}
+		for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
+			for (int q=0; q<18; q++){
+				int nn = neighborList[q*Np+idx]%Np;
+				if (nn>Np) printf("neighborlist error (exterior) at q=%i, idx=%i \n",q,idx);
+		      
+			}
+		}
+
+		//......................device distributions.................................
+		int *NeighborList;
+		int *dvcMap;
+		//...........................................................................
+		ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
+		ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Npad);
+		
+		//...........................................................................
+		// Update GPU data structures
+		if (rank==0)	printf ("Setting up device map and neighbor list \n");
+		int *TmpMap;
+		TmpMap=new int[Np*sizeof(int)];
+		for (k=1; k<Nz-1; k++){
+			for (j=1; j<Ny-1; j++){
+				for (i=1; i<Nx-1; i++){
+					int idx=Map(i,j,k);
+					if (!(idx < 0))
+						TmpMap[idx] = k*Nx*Ny+j*Nx+i;
+				}
+			}
+		}
+		ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
+		ScaLBL_DeviceBarrier();
+		
+		// Create a dummy distribution data structure
+		double *fq;
+		fq = new double[19*Np];
+		if (rank==0)	printf ("Setting up distributions \n");
+		for (k=1; k<Nz-1; k++){
+			for (j=1; j<Ny-1; j++){
+				for (i=1; i<Nx-1; i++){
+					int idx=Map(i,j,k);
+					if (!(idx<0)){
+					  for (int q=0; q<19; q++){
+					    fq[q*Np+idx]=k*100.f+j*10.f+i*1.f+0.01*q;
+					  }
+					}
+				}
+			}
+		}
+		/*		for (int idx=0; idx<Np; idx++){
+			n = TmpMap[idx];
+			// back out the 3D indices
+			k = n/(Nx*Ny);
+			j = (n-Nx*Ny*k)/Nx;
+			i = n-Nx*Ny*k-Nx*j;
+			for (int q=0; q<19; q++){
+				fq[q*Np+idx]=k*100.f+j*10.f+i*1.f+0.01*q;
+			}
+		}
+				if (rank==0)	printf ("Setting up distributions \n");
+		
+		*/
+		// Loop over the distributions for interior lattice sites
+		if (rank==0)	printf ("Loop over distributions \n");
+
+		for (int idx=ScaLBL_Comm->first_interior; idx<ScaLBL_Comm->last_interior; idx++){
+			n = TmpMap[idx];
+			k = n/(Nx*Ny);
+			j = (n-Nx*Ny*k)/Nx;
+			i = n-Nx*Ny*k-Nx*j;
+			for (int q=1; q<19; q++){
+				int nn = neighborList[(q-1)*Np+idx];
+				double value=fq[nn];
+				// 3D index of neighbor
+				int iq=i-D3Q19[q-1][0];
+				int jq=j-D3Q19[q-1][1];
+				int kq=k-D3Q19[q-1][2];
+				if (iq==0) iq=1;
+				if (jq==0) jq=1;
+				if (kq==0) kq=1;
+				if (iq==Nx-1) iq=Nx-2;
+				if (jq==Ny-1) jq=Ny-2;
+				if (kq==Nz-1) kq=Nz-2;
+				double check = kq*100.f+jq*10.f+iq*1.f+q*0.01;
+				if (value != check)
+				  printf("Neighbor q=%i, i=%i,j=%i,k=%i: %f \n",q,iq,jq,kq,value);
+			}
+		}
+
+		delete [] TmpMap;
+
+	}
+    Utilities::shutdown();
+
+	return check;
+}
+