From 6d82dc83a2655cbc7b7c53a9315d6702b0e2d115 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 28 Dec 2021 16:12:49 -0500 Subject: [PATCH] add membrane unit test --- analysis/Membrane.cpp | 33 +++++-- analysis/Membrane.h | 6 +- tests/CMakeLists.txt | 1 + tests/TestMembrane.cpp | 198 +++++++++++++++++++++++++++++++++++++++++ 4 files changed, 230 insertions(+), 8 deletions(-) create mode 100644 tests/TestMembrane.cpp diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index b63bab20..6b9a99d8 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -343,12 +343,22 @@ 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 - const int *originalSendList, const int shift_x, const int shift_y, const int shift_z, +int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int rank_q, int rank_Q, + 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; @@ -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; - membraneSendList[sendCount-membraneCount] = n; + RecvList_q[sendCount-membraneCount] = nn; + SendList[sendCount-membraneCount] = n; } else { - membraneRecvList[regularCount] = nn; - membraneSendList[regularCount++] = n; + RecvList_q[regularCount] = nn; + 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; } diff --git a/analysis/Membrane.h b/analysis/Membrane.h index fe207614..92432c86 100644 --- a/analysis/Membrane.h +++ b/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, - DoubleArray &Distance, double *recvDistance, int *d3q19_recvlist); + int D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int rank_q, int rank_Q, + 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 //...................................................................................... diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5e4441e3..c391613b 100755 --- a/tests/CMakeLists.txt +++ b/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 ) diff --git a/tests/TestMembrane.cpp b/tests/TestMembrane.cpp new file mode 100644 index 00000000..bd7e8234 --- /dev/null +++ b/tests/TestMembrane.cpp @@ -0,0 +1,198 @@ + +//************************************************************************* +// Lattice Boltzmann Simulator for Single Phase Flow in Porous Media +// James E. McCLure +//************************************************************************* +#include +#include +#include +#include "common/ScaLBL.h" +#include "common/MPI.h" + +using namespace std; + +std::shared_ptr loadInputs( int nprocs ) +{ + //auto db = std::make_shared( "Domain.in" ); + auto db = std::make_shared(); + db->putScalar( "BC", 0 ); + db->putVector( "nproc", { 1, 1, 1 } ); + db->putVector( "n", { 5, 5, 5 } ); + db->putScalar( "nspheres", 1 ); + db->putVector( "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( "n" )[0]; + int Ny = db->getVector( "n" )[1]; + int Nz = db->getVector( "n" )[2]; + auto Dm = std::make_shared(db,comm); + + Nx += 2; + Ny += 2; + Nz += 2; + int N = Nx*Ny*Nz; + //....................................................................... + int Np = 0; + for (k=1;kid[n] = 1; + Np++; + } + } + } + Dm->CommInit(); + + // Create a communicator for the device (will use optimized layout) + std::shared_ptr ScaLBL_Comm(new ScaLBL_Communicator(Dm)); + //Create a second communicator based on the regular data layout + std::shared_ptr 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; idxLastExterior(); 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(); idxLastInterior(); 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; kfirst_interior; idxlast_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; +} +