OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

merging multi-halo

Browse Source
This commit is contained in:
JamesEMcclure
2021-01-13 21:44:23 -05:00
parent 29e486a1a0 64b49f720e
commit 0fe3125833
184 changed files with 22054 additions and 7588 deletions
Download Patch File Download Diff File Expand all files Collapse all files

230
common/Communication.h
View File

@@ -1,7 +1,7 @@
#ifndef COMMUNICATION_H_INC
#define COMMUNICATION_H_INC
#include "common/MPI_Helpers.h"
#include "common/MPI.h"
#include "common/Utilities.h"
#include "common/Array.h"
@@ -38,7 +38,7 @@ struct RankInfoStruct {
//! Redistribute domain data (dst may be smaller than the src)
template<class TYPE>
Array<TYPE> redistribute( const RankInfoStruct& src_rank, const Array<TYPE>& src_data,
const RankInfoStruct& dst_rank, std::array<int,3> dst_size, MPI_Comm comm );
const RankInfoStruct& dst_rank, std::array<int,3> dst_size, const Utilities::MPI& comm );
/*!
@@ -59,7 +59,7 @@ public:
* @param[in] fill Fill {faces,edges,corners}
* @param[in] periodic Periodic dimensions
*/
fillHalo( MPI_Comm comm, const RankInfoStruct& info,
fillHalo( const Utilities::MPI& comm, const RankInfoStruct& info,
std::array<int,3> n, std::array<int,3> ng, int tag, int depth,
std::array<bool,3> fill = {true,true,true},
std::array<bool,3> periodic = {true,true,true} );
@@ -83,7 +83,7 @@ public:
private:
MPI_Comm comm;
Utilities::MPI comm;
RankInfoStruct info;
std::array<int,3> n, ng;
int depth;
@@ -93,8 +93,6 @@ private:
TYPE *mem;
TYPE *send[3][3][3], *recv[3][3][3];
MPI_Request send_req[3][3][3], recv_req[3][3][3];
size_t N_type;
MPI_Datatype datatype;
fillHalo(); // Private empty constructor
fillHalo(const fillHalo&); // Private copy constructor
fillHalo& operator=(const fillHalo&); // Private assignment operator
@@ -104,7 +102,7 @@ private:
//***************************************************************************************
inline void PackMeshData(int *list, int count, double *sendbuf, double *data){
inline void PackMeshData(const int *list, int count, double *sendbuf, double *data){
// Fill in the phase ID values from neighboring processors
// This packs up the values that need to be sent from one processor to another
int idx,n;
@@ -113,7 +111,7 @@ inline void PackMeshData(int *list, int count, double *sendbuf, double *data){
sendbuf[idx] = data[n];
}
}
inline void UnpackMeshData(int *list, int count, double *recvbuf, double *data){
inline void UnpackMeshData(const int *list, int count, double *recvbuf, double *data){
// Fill in the phase ID values from neighboring processors
// This unpacks the values once they have been recieved from neighbors
int idx,n;
@@ -136,7 +134,7 @@ void InitializeRanks( const int rank, const int nprocx, const int nprocy, const
//***************************************************************************************
inline void CommunicateSendRecvCounts( MPI_Comm Communicator, int sendtag, int recvtag,
inline void CommunicateSendRecvCounts( const Utilities::MPI& comm, int sendtag, int recvtag,
int rank_x, int rank_y, int rank_z,
int rank_X, int rank_Y, int rank_Z,
int rank_xy, int rank_XY, int rank_xY, int rank_Xy,
@@ -154,54 +152,53 @@ inline void CommunicateSendRecvCounts( MPI_Comm Communicator, int sendtag, int r
int& recvCount_yz, int& recvCount_YZ, int& recvCount_yZ, int& recvCount_Yz )
{
MPI_Request req1[18], req2[18];
MPI_Status stat1[18],stat2[18];
MPI_Isend(&sendCount_x, 1,MPI_INT,rank_x,sendtag+0,Communicator,&req1[0]);
MPI_Irecv(&recvCount_X, 1,MPI_INT,rank_X,recvtag+0,Communicator,&req2[0]);
MPI_Isend(&sendCount_X, 1,MPI_INT,rank_X,sendtag+1,Communicator,&req1[1]);
MPI_Irecv(&recvCount_x, 1,MPI_INT,rank_x,recvtag+1,Communicator,&req2[1]);
MPI_Isend(&sendCount_y, 1,MPI_INT,rank_y,sendtag+2,Communicator,&req1[2]);
MPI_Irecv(&recvCount_Y, 1,MPI_INT,rank_Y,recvtag+2,Communicator,&req2[2]);
MPI_Isend(&sendCount_Y, 1,MPI_INT,rank_Y,sendtag+3,Communicator,&req1[3]);
MPI_Irecv(&recvCount_y, 1,MPI_INT,rank_y,recvtag+3,Communicator,&req2[3]);
MPI_Isend(&sendCount_z, 1,MPI_INT,rank_z,sendtag+4,Communicator,&req1[4]);
MPI_Irecv(&recvCount_Z, 1,MPI_INT,rank_Z,recvtag+4,Communicator,&req2[4]);
MPI_Isend(&sendCount_Z, 1,MPI_INT,rank_Z,sendtag+5,Communicator,&req1[5]);
MPI_Irecv(&recvCount_z, 1,MPI_INT,rank_z,recvtag+5,Communicator,&req2[5]);
req1[0] = comm.Isend(&sendCount_x,1,rank_x,sendtag+0);
req2[0] = comm.Irecv(&recvCount_X,1,rank_X,recvtag+0);
req1[1] = comm.Isend(&sendCount_X,1,rank_X,sendtag+1);
req2[1] = comm.Irecv(&recvCount_x,1,rank_x,recvtag+1);
req1[2] = comm.Isend(&sendCount_y,1,rank_y,sendtag+2);
req2[2] = comm.Irecv(&recvCount_Y,1,rank_Y,recvtag+2);
req1[3] = comm.Isend(&sendCount_Y,1,rank_Y,sendtag+3);
req2[3] = comm.Irecv(&recvCount_y,1,rank_y,recvtag+3);
req1[4] = comm.Isend(&sendCount_z,1,rank_z,sendtag+4);
req2[4] = comm.Irecv(&recvCount_Z,1,rank_Z,recvtag+4);
req1[5] = comm.Isend(&sendCount_Z,1,rank_Z,sendtag+5);
req2[5] = comm.Irecv(&recvCount_z,1,rank_z,recvtag+5);
MPI_Isend(&sendCount_xy, 1,MPI_INT,rank_xy,sendtag+6,Communicator,&req1[6]);
MPI_Irecv(&recvCount_XY, 1,MPI_INT,rank_XY,recvtag+6,Communicator,&req2[6]);
MPI_Isend(&sendCount_XY, 1,MPI_INT,rank_XY,sendtag+7,Communicator,&req1[7]);
MPI_Irecv(&recvCount_xy, 1,MPI_INT,rank_xy,recvtag+7,Communicator,&req2[7]);
MPI_Isend(&sendCount_Xy, 1,MPI_INT,rank_Xy,sendtag+8,Communicator,&req1[8]);
MPI_Irecv(&recvCount_xY, 1,MPI_INT,rank_xY,recvtag+8,Communicator,&req2[8]);
MPI_Isend(&sendCount_xY, 1,MPI_INT,rank_xY,sendtag+9,Communicator,&req1[9]);
MPI_Irecv(&recvCount_Xy, 1,MPI_INT,rank_Xy,recvtag+9,Communicator,&req2[9]);
req1[6] = comm.Isend(&sendCount_xy,1,rank_xy,sendtag+6);
req2[6] = comm.Irecv(&recvCount_XY,1,rank_XY,recvtag+6);
req1[7] = comm.Isend(&sendCount_XY,1,rank_XY,sendtag+7);
req2[7] = comm.Irecv(&recvCount_xy,1,rank_xy,recvtag+7);
req1[8] = comm.Isend(&sendCount_Xy,1,rank_Xy,sendtag+8);
req2[8] = comm.Irecv(&recvCount_xY,1,rank_xY,recvtag+8);
req1[9] = comm.Isend(&sendCount_xY,1,rank_xY,sendtag+9);
req2[9] = comm.Irecv(&recvCount_Xy,1,rank_Xy,recvtag+9);
MPI_Isend(&sendCount_xz, 1,MPI_INT,rank_xz,sendtag+10,Communicator,&req1[10]);
MPI_Irecv(&recvCount_XZ, 1,MPI_INT,rank_XZ,recvtag+10,Communicator,&req2[10]);
MPI_Isend(&sendCount_XZ, 1,MPI_INT,rank_XZ,sendtag+11,Communicator,&req1[11]);
MPI_Irecv(&recvCount_xz, 1,MPI_INT,rank_xz,recvtag+11,Communicator,&req2[11]);
MPI_Isend(&sendCount_Xz, 1,MPI_INT,rank_Xz,sendtag+12,Communicator,&req1[12]);
MPI_Irecv(&recvCount_xZ, 1,MPI_INT,rank_xZ,recvtag+12,Communicator,&req2[12]);
MPI_Isend(&sendCount_xZ, 1,MPI_INT,rank_xZ,sendtag+13,Communicator,&req1[13]);
MPI_Irecv(&recvCount_Xz, 1,MPI_INT,rank_Xz,recvtag+13,Communicator,&req2[13]);
req1[10] = comm.Isend(&sendCount_xz,1,rank_xz,sendtag+10);
req2[10] = comm.Irecv(&recvCount_XZ,1,rank_XZ,recvtag+10);
req1[11] = comm.Isend(&sendCount_XZ,1,rank_XZ,sendtag+11);
req2[11] = comm.Irecv(&recvCount_xz,1,rank_xz,recvtag+11);
req1[12] = comm.Isend(&sendCount_Xz,1,rank_Xz,sendtag+12);
req2[12] = comm.Irecv(&recvCount_xZ,1,rank_xZ,recvtag+12);
req1[13] = comm.Isend(&sendCount_xZ,1,rank_xZ,sendtag+13);
req2[13] = comm.Irecv(&recvCount_Xz,1,rank_Xz,recvtag+13);
MPI_Isend(&sendCount_yz, 1,MPI_INT,rank_yz,sendtag+14,Communicator,&req1[14]);
MPI_Irecv(&recvCount_YZ, 1,MPI_INT,rank_YZ,recvtag+14,Communicator,&req2[14]);
MPI_Isend(&sendCount_YZ, 1,MPI_INT,rank_YZ,sendtag+15,Communicator,&req1[15]);
MPI_Irecv(&recvCount_yz, 1,MPI_INT,rank_yz,recvtag+15,Communicator,&req2[15]);
MPI_Isend(&sendCount_Yz, 1,MPI_INT,rank_Yz,sendtag+16,Communicator,&req1[16]);
MPI_Irecv(&recvCount_yZ, 1,MPI_INT,rank_yZ,recvtag+16,Communicator,&req2[16]);
MPI_Isend(&sendCount_yZ, 1,MPI_INT,rank_yZ,sendtag+17,Communicator,&req1[17]);
MPI_Irecv(&recvCount_Yz, 1,MPI_INT,rank_Yz,recvtag+17,Communicator,&req2[17]);
MPI_Waitall(18,req1,stat1);
MPI_Waitall(18,req2,stat2);
MPI_Barrier(Communicator);
req1[14] = comm.Isend(&sendCount_yz,1,rank_yz,sendtag+14);
req2[14] = comm.Irecv(&recvCount_YZ,1,rank_YZ,recvtag+14);
req1[15] = comm.Isend(&sendCount_YZ,1,rank_YZ,sendtag+15);
req2[15] = comm.Irecv(&recvCount_yz,1,rank_yz,recvtag+15);
req1[16] = comm.Isend(&sendCount_Yz,1,rank_Yz,sendtag+16);
req2[16] = comm.Irecv(&recvCount_yZ,1,rank_yZ,recvtag+16);
req1[17] = comm.Isend(&sendCount_yZ,1,rank_yZ,sendtag+17);
req2[17] = comm.Irecv(&recvCount_Yz,1,rank_Yz,recvtag+17);
comm.waitAll( 18, req1 );
comm.waitAll( 18, req2 );
comm.barrier();
}
//***************************************************************************************
inline void CommunicateRecvLists( MPI_Comm Communicator, int sendtag, int recvtag,
inline void CommunicateRecvLists( const Utilities::MPI& comm, int sendtag, int recvtag,
int *sendList_x, int *sendList_y, int *sendList_z, int *sendList_X, int *sendList_Y, int *sendList_Z,
int *sendList_xy, int *sendList_XY, int *sendList_xY, int *sendList_Xy,
int *sendList_xz, int *sendList_XZ, int *sendList_xZ, int *sendList_Xz,
@@ -222,53 +219,52 @@ inline void CommunicateRecvLists( MPI_Comm Communicator, int sendtag, int recvta
int rank_Xy, int rank_xz, int rank_XZ, int rank_xZ, int rank_Xz, int rank_yz, int rank_YZ, int rank_yZ, int rank_Yz)
{
MPI_Request req1[18], req2[18];
MPI_Status stat1[18],stat2[18];
MPI_Isend(sendList_x, sendCount_x,MPI_INT,rank_x,sendtag,Communicator,&req1[0]);
MPI_Irecv(recvList_X, recvCount_X,MPI_INT,rank_X,recvtag,Communicator,&req2[0]);
MPI_Isend(sendList_X, sendCount_X,MPI_INT,rank_X,sendtag,Communicator,&req1[1]);
MPI_Irecv(recvList_x, recvCount_x,MPI_INT,rank_x,recvtag,Communicator,&req2[1]);
MPI_Isend(sendList_y, sendCount_y,MPI_INT,rank_y,sendtag,Communicator,&req1[2]);
MPI_Irecv(recvList_Y, recvCount_Y,MPI_INT,rank_Y,recvtag,Communicator,&req2[2]);
MPI_Isend(sendList_Y, sendCount_Y,MPI_INT,rank_Y,sendtag,Communicator,&req1[3]);
MPI_Irecv(recvList_y, recvCount_y,MPI_INT,rank_y,recvtag,Communicator,&req2[3]);
MPI_Isend(sendList_z, sendCount_z,MPI_INT,rank_z,sendtag,Communicator,&req1[4]);
MPI_Irecv(recvList_Z, recvCount_Z,MPI_INT,rank_Z,recvtag,Communicator,&req2[4]);
MPI_Isend(sendList_Z, sendCount_Z,MPI_INT,rank_Z,sendtag,Communicator,&req1[5]);
MPI_Irecv(recvList_z, recvCount_z,MPI_INT,rank_z,recvtag,Communicator,&req2[5]);
req1[0] = comm.Isend(sendList_x,sendCount_x,rank_x,sendtag);
req2[0] = comm.Irecv(recvList_X,recvCount_X,rank_X,recvtag);
req1[1] = comm.Isend(sendList_X,sendCount_X,rank_X,sendtag);
req2[1] = comm.Irecv(recvList_x,recvCount_x,rank_x,recvtag);
req1[2] = comm.Isend(sendList_y,sendCount_y,rank_y,sendtag);
req2[2] = comm.Irecv(recvList_Y,recvCount_Y,rank_Y,recvtag);
req1[3] = comm.Isend(sendList_Y,sendCount_Y,rank_Y,sendtag);
req2[3] = comm.Irecv(recvList_y,recvCount_y,rank_y,recvtag);
req1[4] = comm.Isend(sendList_z,sendCount_z,rank_z,sendtag);
req2[4] = comm.Irecv(recvList_Z,recvCount_Z,rank_Z,recvtag);
req1[5] = comm.Isend(sendList_Z,sendCount_Z,rank_Z,sendtag);
req2[5] = comm.Irecv(recvList_z,recvCount_z,rank_z,recvtag);
MPI_Isend(sendList_xy, sendCount_xy,MPI_INT,rank_xy,sendtag,Communicator,&req1[6]);
MPI_Irecv(recvList_XY, recvCount_XY,MPI_INT,rank_XY,recvtag,Communicator,&req2[6]);
MPI_Isend(sendList_XY, sendCount_XY,MPI_INT,rank_XY,sendtag,Communicator,&req1[7]);
MPI_Irecv(recvList_xy, recvCount_xy,MPI_INT,rank_xy,recvtag,Communicator,&req2[7]);
MPI_Isend(sendList_Xy, sendCount_Xy,MPI_INT,rank_Xy,sendtag,Communicator,&req1[8]);
MPI_Irecv(recvList_xY, recvCount_xY,MPI_INT,rank_xY,recvtag,Communicator,&req2[8]);
MPI_Isend(sendList_xY, sendCount_xY,MPI_INT,rank_xY,sendtag,Communicator,&req1[9]);
MPI_Irecv(recvList_Xy, recvCount_Xy,MPI_INT,rank_Xy,recvtag,Communicator,&req2[9]);
req1[6] = comm.Isend(sendList_xy,sendCount_xy,rank_xy,sendtag);
req2[6] = comm.Irecv(recvList_XY,recvCount_XY,rank_XY,recvtag);
req1[7] = comm.Isend(sendList_XY,sendCount_XY,rank_XY,sendtag);
req2[7] = comm.Irecv(recvList_xy,recvCount_xy,rank_xy,recvtag);
req1[8] = comm.Isend(sendList_Xy,sendCount_Xy,rank_Xy,sendtag);
req2[8] = comm.Irecv(recvList_xY,recvCount_xY,rank_xY,recvtag);
req1[9] = comm.Isend(sendList_xY,sendCount_xY,rank_xY,sendtag);
req2[9] = comm.Irecv(recvList_Xy,recvCount_Xy,rank_Xy,recvtag);
MPI_Isend(sendList_xz, sendCount_xz,MPI_INT,rank_xz,sendtag,Communicator,&req1[10]);
MPI_Irecv(recvList_XZ, recvCount_XZ,MPI_INT,rank_XZ,recvtag,Communicator,&req2[10]);
MPI_Isend(sendList_XZ, sendCount_XZ,MPI_INT,rank_XZ,sendtag,Communicator,&req1[11]);
MPI_Irecv(recvList_xz, recvCount_xz,MPI_INT,rank_xz,recvtag,Communicator,&req2[11]);
MPI_Isend(sendList_Xz, sendCount_Xz,MPI_INT,rank_Xz,sendtag,Communicator,&req1[12]);
MPI_Irecv(recvList_xZ, recvCount_xZ,MPI_INT,rank_xZ,recvtag,Communicator,&req2[12]);
MPI_Isend(sendList_xZ, sendCount_xZ,MPI_INT,rank_xZ,sendtag,Communicator,&req1[13]);
MPI_Irecv(recvList_Xz, recvCount_Xz,MPI_INT,rank_Xz,recvtag,Communicator,&req2[13]);
req1[10] = comm.Isend(sendList_xz,sendCount_xz,rank_xz,sendtag);
req2[10] = comm.Irecv(recvList_XZ,recvCount_XZ,rank_XZ,recvtag);
req1[11] = comm.Isend(sendList_XZ,sendCount_XZ,rank_XZ,sendtag);
req2[11] = comm.Irecv(recvList_xz,recvCount_xz,rank_xz,recvtag);
req1[12] = comm.Isend(sendList_Xz,sendCount_Xz,rank_Xz,sendtag);
req2[12] = comm.Irecv(recvList_xZ,recvCount_xZ,rank_xZ,recvtag);
req1[13] = comm.Isend(sendList_xZ,sendCount_xZ,rank_xZ,sendtag);
req2[13] = comm.Irecv(recvList_Xz,recvCount_Xz,rank_Xz,recvtag);
MPI_Isend(sendList_yz, sendCount_yz,MPI_INT,rank_yz,sendtag,Communicator,&req1[14]);
MPI_Irecv(recvList_YZ, recvCount_YZ,MPI_INT,rank_YZ,recvtag,Communicator,&req2[14]);
MPI_Isend(sendList_YZ, sendCount_YZ,MPI_INT,rank_YZ,sendtag,Communicator,&req1[15]);
MPI_Irecv(recvList_yz, recvCount_yz,MPI_INT,rank_yz,recvtag,Communicator,&req2[15]);
MPI_Isend(sendList_Yz, sendCount_Yz,MPI_INT,rank_Yz,sendtag,Communicator,&req1[16]);
MPI_Irecv(recvList_yZ, recvCount_yZ,MPI_INT,rank_yZ,recvtag,Communicator,&req2[16]);
MPI_Isend(sendList_yZ, sendCount_yZ,MPI_INT,rank_yZ,sendtag,Communicator,&req1[17]);
MPI_Irecv(recvList_Yz, recvCount_Yz,MPI_INT,rank_Yz,recvtag,Communicator,&req2[17]);
MPI_Waitall(18,req1,stat1);
MPI_Waitall(18,req2,stat2);
req1[14] = comm.Isend(sendList_yz,sendCount_yz,rank_yz,sendtag);
req2[14] = comm.Irecv(recvList_YZ,recvCount_YZ,rank_YZ,recvtag);
req1[15] = comm.Isend(sendList_YZ,sendCount_YZ,rank_YZ,sendtag);
req2[15] = comm.Irecv(recvList_yz,recvCount_yz,rank_yz,recvtag);
req1[16] = comm.Isend(sendList_Yz,sendCount_Yz,rank_Yz,sendtag);
req2[16] = comm.Irecv(recvList_yZ,recvCount_yZ,rank_yZ,recvtag);
req1[17] = comm.Isend(sendList_yZ,sendCount_yZ,rank_yZ,sendtag);
req2[17] = comm.Irecv(recvList_Yz,recvCount_Yz,rank_Yz,recvtag);
comm.waitAll( 18, req1 );
comm.waitAll( 18, req2 );
}
//***************************************************************************************
inline void CommunicateMeshHalo(DoubleArray &Mesh, MPI_Comm Communicator,
inline void CommunicateMeshHalo(DoubleArray &Mesh, const Utilities::MPI& comm,
double *sendbuf_x,double *sendbuf_y,double *sendbuf_z,double *sendbuf_X,double *sendbuf_Y,double *sendbuf_Z,
double *sendbuf_xy,double *sendbuf_XY,double *sendbuf_xY,double *sendbuf_Xy,
double *sendbuf_xz,double *sendbuf_XZ,double *sendbuf_xZ,double *sendbuf_Xz,
@@ -318,42 +314,24 @@ inline void CommunicateMeshHalo(DoubleArray &Mesh, MPI_Comm Communicator,
PackMeshData(sendList_yZ, sendCount_yZ ,sendbuf_yZ, MeshData);
PackMeshData(sendList_YZ, sendCount_YZ ,sendbuf_YZ, MeshData);
//......................................................................................
MPI_Sendrecv(sendbuf_x,sendCount_x,MPI_DOUBLE,rank_x,sendtag,
recvbuf_X,recvCount_X,MPI_DOUBLE,rank_X,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_X,sendCount_X,MPI_DOUBLE,rank_X,sendtag,
recvbuf_x,recvCount_x,MPI_DOUBLE,rank_x,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_y,sendCount_y,MPI_DOUBLE,rank_y,sendtag,
recvbuf_Y,recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_Y,sendCount_Y,MPI_DOUBLE,rank_Y,sendtag,
recvbuf_y,recvCount_y,MPI_DOUBLE,rank_y,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_z,sendCount_z,MPI_DOUBLE,rank_z,sendtag,
recvbuf_Z,recvCount_Z,MPI_DOUBLE,rank_Z,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_Z,sendCount_Z,MPI_DOUBLE,rank_Z,sendtag,
recvbuf_z,recvCount_z,MPI_DOUBLE,rank_z,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_xy,sendCount_xy,MPI_DOUBLE,rank_xy,sendtag,
recvbuf_XY,recvCount_XY,MPI_DOUBLE,rank_XY,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_XY,sendCount_XY,MPI_DOUBLE,rank_XY,sendtag,
recvbuf_xy,recvCount_xy,MPI_DOUBLE,rank_xy,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_Xy,sendCount_Xy,MPI_DOUBLE,rank_Xy,sendtag,
recvbuf_xY,recvCount_xY,MPI_DOUBLE,rank_xY,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_xY,sendCount_xY,MPI_DOUBLE,rank_xY,sendtag,
recvbuf_Xy,recvCount_Xy,MPI_DOUBLE,rank_Xy,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_xz,sendCount_xz,MPI_DOUBLE,rank_xz,sendtag,
recvbuf_XZ,recvCount_XZ,MPI_DOUBLE,rank_XZ,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_XZ,sendCount_XZ,MPI_DOUBLE,rank_XZ,sendtag,
recvbuf_xz,recvCount_xz,MPI_DOUBLE,rank_xz,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_Xz,sendCount_Xz,MPI_DOUBLE,rank_Xz,sendtag,
recvbuf_xZ,recvCount_xZ,MPI_DOUBLE,rank_xZ,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_xZ,sendCount_xZ,MPI_DOUBLE,rank_xZ,sendtag,
recvbuf_Xz,recvCount_Xz,MPI_DOUBLE,rank_Xz,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_yz,sendCount_yz,MPI_DOUBLE,rank_yz,sendtag,
recvbuf_YZ,recvCount_YZ,MPI_DOUBLE,rank_YZ,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_YZ,sendCount_YZ,MPI_DOUBLE,rank_YZ,sendtag,
recvbuf_yz,recvCount_yz,MPI_DOUBLE,rank_yz,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_Yz,sendCount_Yz,MPI_DOUBLE,rank_Yz,sendtag,
recvbuf_yZ,recvCount_yZ,MPI_DOUBLE,rank_yZ,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_yZ,sendCount_yZ,MPI_DOUBLE,rank_yZ,sendtag,
recvbuf_Yz,recvCount_Yz,MPI_DOUBLE,rank_Yz,recvtag,Communicator,MPI_STATUS_IGNORE);
comm.sendrecv(sendbuf_x,sendCount_x,rank_x,sendtag,recvbuf_X,recvCount_X,rank_X,recvtag);
comm.sendrecv(sendbuf_X,sendCount_X,rank_X,sendtag,recvbuf_x,recvCount_x,rank_x,recvtag);
comm.sendrecv(sendbuf_y,sendCount_y,rank_y,sendtag,recvbuf_Y,recvCount_Y,rank_Y,recvtag);
comm.sendrecv(sendbuf_Y,sendCount_Y,rank_Y,sendtag,recvbuf_y,recvCount_y,rank_y,recvtag);
comm.sendrecv(sendbuf_z,sendCount_z,rank_z,sendtag,recvbuf_Z,recvCount_Z,rank_Z,recvtag);
comm.sendrecv(sendbuf_Z,sendCount_Z,rank_Z,sendtag,recvbuf_z,recvCount_z,rank_z,recvtag);
comm.sendrecv(sendbuf_xy,sendCount_xy,rank_xy,sendtag,recvbuf_XY,recvCount_XY,rank_XY,recvtag);
comm.sendrecv(sendbuf_XY,sendCount_XY,rank_XY,sendtag,recvbuf_xy,recvCount_xy,rank_xy,recvtag);
comm.sendrecv(sendbuf_Xy,sendCount_Xy,rank_Xy,sendtag,recvbuf_xY,recvCount_xY,rank_xY,recvtag);
comm.sendrecv(sendbuf_xY,sendCount_xY,rank_xY,sendtag,recvbuf_Xy,recvCount_Xy,rank_Xy,recvtag);
comm.sendrecv(sendbuf_xz,sendCount_xz,rank_xz,sendtag,recvbuf_XZ,recvCount_XZ,rank_XZ,recvtag);
comm.sendrecv(sendbuf_XZ,sendCount_XZ,rank_XZ,sendtag,recvbuf_xz,recvCount_xz,rank_xz,recvtag);
comm.sendrecv(sendbuf_Xz,sendCount_Xz,rank_Xz,sendtag,recvbuf_xZ,recvCount_xZ,rank_xZ,recvtag);
comm.sendrecv(sendbuf_xZ,sendCount_xZ,rank_xZ,sendtag,recvbuf_Xz,recvCount_Xz,rank_Xz,recvtag);
comm.sendrecv(sendbuf_yz,sendCount_yz,rank_yz,sendtag,recvbuf_YZ,recvCount_YZ,rank_YZ,recvtag);
comm.sendrecv(sendbuf_YZ,sendCount_YZ,rank_YZ,sendtag,recvbuf_yz,recvCount_yz,rank_yz,recvtag);
comm.sendrecv(sendbuf_Yz,sendCount_Yz,rank_Yz,sendtag,recvbuf_yZ,recvCount_yZ,rank_yZ,recvtag);
comm.sendrecv(sendbuf_yZ,sendCount_yZ,rank_yZ,sendtag,recvbuf_Yz,recvCount_Yz,rank_Yz,recvtag);
//........................................................................................
UnpackMeshData(recvList_x, recvCount_x ,recvbuf_x, MeshData);
UnpackMeshData(recvList_X, recvCount_X ,recvbuf_X, MeshData);

53
common/Communication.hpp
View File

@@ -2,9 +2,8 @@
#define COMMUNICATION_HPP_INC
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "common/MPI.h"
#include "common/Utilities.h"
//#include "ProfilerApp.h"
/********************************************************
@@ -12,17 +11,19 @@
********************************************************/
template<class TYPE>
Array<TYPE> redistribute( const RankInfoStruct& src_rank, const Array<TYPE>& src_data,
const RankInfoStruct& dst_rank, std::array<int,3> dst_size, MPI_Comm comm )
const RankInfoStruct& dst_rank, std::array<int,3> dst_size, const Utilities::MPI& comm )
{
#ifdef USE_MPI
if ( comm.getSize() == 1 ) {
return src_data.subset( { 0, (size_t) dst_size[0]-1, 0, (size_t) dst_size[1]-1, 0, (size_t) dst_size[2]-1 } );
}
// Get the src size
std::array<int,3> src_size;
int size0[3] = { (int) src_data.size(0), (int) src_data.size(1), (int) src_data.size(2) };
MPI_Allreduce( size0, src_size.data(), 3, MPI_INT, MPI_MAX, comm );
comm.maxReduce( size0, src_size.data(), 3 );
if ( !src_data.empty() )
ASSERT( src_size[0] == size0[0] && src_size[1] == size0[1] && src_size[2] == size0[2] );
// Check that dst_size matches on all ranks
MPI_Allreduce( dst_size.data(), size0, 3, MPI_INT, MPI_MAX, comm );
comm.maxReduce( dst_size.data(), size0, 3 );
ASSERT( dst_size[0] == size0[0] && dst_size[1] == size0[1] && dst_size[2] == size0[2] );
// Function to get overlap range
auto calcOverlap = []( int i1[3], int i2[3], int j1[3], int j2[3] ) {
@@ -60,7 +61,7 @@ Array<TYPE> redistribute( const RankInfoStruct& src_rank, const Array<TYPE>& src
}
std::vector<MPI_Request> send_request( send_rank.size() );
for (size_t i=0; i<send_rank.size(); i++)
MPI_Isend( send_data[i].data(), sizeof(TYPE)*send_data[i].length(), MPI_BYTE, send_rank[i], 5462, comm, &send_request[i]);
send_request[i] = comm.Isend( send_data[i].data(), send_data[i].length(), send_rank[i], 5462 );
// Unpack data from the appropriate ranks (including myself)
Array<TYPE> dst_data( dst_size[0], dst_size[1], dst_size[2] );
int i1[3] = { dst_size[0] * dst_rank.ix, dst_size[1] * dst_rank.jy, dst_size[2] * dst_rank.kz };
@@ -75,17 +76,14 @@ Array<TYPE> redistribute( const RankInfoStruct& src_rank, const Array<TYPE>& src
continue;
int rank = src_rank.getRankForBlock(i,j,k);
Array<TYPE> data( index[1] - index[0] + 1, index[3] - index[2] + 1, index[5] - index[4] + 1 );
MPI_Recv( data.data(), sizeof(TYPE)*data.length(), MPI_BYTE, rank, 5462, comm, MPI_STATUS_IGNORE );
comm.recv( data.data(), data.length(), rank, 5462 );
dst_data.copySubset( index, data );
}
}
}
// Free data
MPI_Waitall( send_request.size(), send_request.data(), MPI_STATUSES_IGNORE );
comm.waitAll( send_request.size(), send_request.data() );
return dst_data;
#else
return src_data.subset( { 0, dst_size[0]-1, 0, dst_size[1]-1, 0, dst_size[2]-1 );
#endif
}
@@ -94,27 +92,11 @@ Array<TYPE> redistribute( const RankInfoStruct& src_rank, const Array<TYPE>& src
* Structure to fill halo cells *
********************************************************/
template<class TYPE>
fillHalo<TYPE>::fillHalo( MPI_Comm comm_, const RankInfoStruct& info_,
fillHalo<TYPE>::fillHalo( const Utilities::MPI& comm_, const RankInfoStruct& info_,
std::array<int,3> n_, std::array<int,3> ng_, int tag0, int depth_,
std::array<bool,3> fill, std::array<bool,3> periodic ):
comm(comm_), info(info_), n(n_), ng(ng_), depth(depth_)
{
if ( std::is_same<TYPE,double>() ) {
N_type = 1;
datatype = MPI_DOUBLE;
} else if ( std::is_same<TYPE,float>() ) {
N_type = 1;
datatype = MPI_FLOAT;
} else if ( sizeof(TYPE)%sizeof(double)==0 ) {
N_type = sizeof(TYPE) / sizeof(double);
datatype = MPI_DOUBLE;
} else if ( sizeof(TYPE)%sizeof(float)==0 ) {
N_type = sizeof(TYPE) / sizeof(float);
datatype = MPI_FLOAT;
} else {
N_type = sizeof(TYPE);
datatype = MPI_BYTE;
}
// Set the fill pattern
memset(fill_pattern,0,sizeof(fill_pattern));
if ( fill[0] ) {
@@ -251,8 +233,8 @@ void fillHalo<TYPE>::fill( Array<TYPE>& data )
for (int k=0; k<3; k++) {
if ( !fill_pattern[i][j][k] )
continue;
MPI_Irecv( recv[i][j][k], N_type*depth2*N_send_recv[i][j][k], datatype,
info.rank[i][j][k], tag[2-i][2-j][2-k], comm, &recv_req[i][j][k] );
recv_req[i][j][k] = comm.Irecv( recv[i][j][k], depth2*N_send_recv[i][j][k],
info.rank[i][j][k], tag[2-i][2-j][2-k] );
}
}
}
@@ -263,19 +245,18 @@ void fillHalo<TYPE>::fill( Array<TYPE>& data )
if ( !fill_pattern[i][j][k] )
continue;
pack( data, i-1, j-1, k-1, send[i][j][k] );
MPI_Isend( send[i][j][k], N_type*depth2*N_send_recv[i][j][k], datatype,
info.rank[i][j][k], tag[i][j][k], comm, &send_req[i][j][k] );
send_req[i][j][k] = comm.Isend( send[i][j][k], depth2*N_send_recv[i][j][k],
info.rank[i][j][k], tag[i][j][k] );
}
}
}
// Recv the dst data and unpack (we recive in reverse order to match the sends)
MPI_Status status;
for (int i=2; i>=0; i--) {
for (int j=2; j>=0; j--) {
for (int k=2; k>=0; k--) {
if ( !fill_pattern[i][j][k] )
continue;
MPI_Wait(&recv_req[i][j][k],&status);
comm.wait( recv_req[i][j][k] );
unpack( data, i-1, j-1, k-1, recv[i][j][k] );
}
}
@@ -286,7 +267,7 @@ void fillHalo<TYPE>::fill( Array<TYPE>& data )
for (int k=0; k<3; k++) {
if ( !fill_pattern[i][j][k] )
continue;
MPI_Wait(&send_req[i][j][k],&status);
comm.wait( send_req[i][j][k] );
}
}
}

843
common/Domain.cpp
View File

File diff suppressed because it is too large Load Diff

60
common/Domain.h
View File

@@ -5,6 +5,7 @@
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>
#include <time.h>
#include <exception>
@@ -12,7 +13,7 @@
#include "common/Array.h"
#include "common/Utilities.h"
#include "common/MPI_Helpers.h"
#include "common/MPI.h"
#include "common/Communication.h"
#include "common/Database.h"
@@ -63,7 +64,7 @@ private:
class Domain{
public:
//! Default constructor
Domain( std::shared_ptr<Database> db, MPI_Comm Communicator);
Domain( std::shared_ptr<Database> db, const Utilities::MPI& Communicator);
//! Obsolete constructor
Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz,
@@ -116,12 +117,10 @@ public: // Public variables (need to create accessors instead)
double porosity;
RankInfoStruct rank_info;
MPI_Comm Comm; // MPI Communicator for this domain
Utilities::MPI Comm; // MPI Communicator for this domain
int BoundaryCondition;
MPI_Group Group; // Group of processors associated with this domain
//**********************************
// MPI ranks for all 18 neighbors
//**********************************
@@ -157,31 +156,22 @@ public: // Public variables (need to create accessors instead)
// Get the actual D3Q19 communication counts (based on location of solid phase)
// Discrete velocity set symmetry implies the sendcount = recvcount
//......................................................................................
int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z;
int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ;
int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ;
//......................................................................................
int *sendList_x, *sendList_y, *sendList_z, *sendList_X, *sendList_Y, *sendList_Z;
int *sendList_xy, *sendList_yz, *sendList_xz, *sendList_Xy, *sendList_Yz, *sendList_xZ;
int *sendList_xY, *sendList_yZ, *sendList_Xz, *sendList_XY, *sendList_YZ, *sendList_XZ;
//......................................................................................
int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z;
int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ;
int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ;
//......................................................................................
int *recvList_x, *recvList_y, *recvList_z, *recvList_X, *recvList_Y, *recvList_Z;
int *recvList_xy, *recvList_yz, *recvList_xz, *recvList_Xy, *recvList_Yz, *recvList_xZ;
int *recvList_xY, *recvList_yZ, *recvList_Xz, *recvList_XY, *recvList_YZ, *recvList_XZ;
inline int recvCount( const char* dir ) const { return getRecvList( dir ).size(); }
inline int sendCount( const char* dir ) const { return getSendList( dir ).size(); }
inline const int* recvList( const char* dir ) const { return getRecvList( dir ).data(); }
inline const int* sendList( const char* dir ) const { return getSendList( dir ).data(); }
//......................................................................................
// Solid indicator function
signed char *id;
std::vector<signed char> id;
void ReadIDs();
void Decomp( const std::string& filename );
void ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData);
void CommunicateMeshHalo(DoubleArray &Mesh);
void CommInit();
int PoreCount();
void ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData);
void AggregateLabels( const std::string& filename );
void AggregateLabels( const std::string& filename, DoubleArray &UserData );
@@ -193,22 +183,18 @@ private:
//......................................................................................
MPI_Request req1[18], req2[18];
MPI_Status stat1[18],stat2[18];
//......................................................................................
std::vector<int> sendList_x, sendList_y, sendList_z, sendList_X, sendList_Y, sendList_Z;
std::vector<int> sendList_xy, sendList_yz, sendList_xz, sendList_Xy, sendList_Yz, sendList_xZ;
std::vector<int> sendList_xY, sendList_yZ, sendList_Xz, sendList_XY, sendList_YZ, sendList_XZ;
//......................................................................................
std::vector<int> recvList_x, recvList_y, recvList_z, recvList_X, recvList_Y, recvList_Z;
std::vector<int> recvList_xy, recvList_yz, recvList_xz, recvList_Xy, recvList_Yz, recvList_xZ;
std::vector<int> recvList_xY, recvList_yZ, recvList_Xz, recvList_XY, recvList_YZ, recvList_XZ;
//......................................................................................
const std::vector<int>& getRecvList( const char* dir ) const;
const std::vector<int>& getSendList( const char* dir ) const;
int *sendBuf_x, *sendBuf_y, *sendBuf_z, *sendBuf_X, *sendBuf_Y, *sendBuf_Z;
int *sendBuf_xy, *sendBuf_yz, *sendBuf_xz, *sendBuf_Xy, *sendBuf_Yz, *sendBuf_xZ;
int *sendBuf_xY, *sendBuf_yZ, *sendBuf_Xz, *sendBuf_XY, *sendBuf_YZ, *sendBuf_XZ;
//......................................................................................
int *recvBuf_x, *recvBuf_y, *recvBuf_z, *recvBuf_X, *recvBuf_Y, *recvBuf_Z;
int *recvBuf_xy, *recvBuf_yz, *recvBuf_xz, *recvBuf_Xy, *recvBuf_Yz, *recvBuf_xZ;
int *recvBuf_xY, *recvBuf_yZ, *recvBuf_Xz, *recvBuf_XY, *recvBuf_YZ, *recvBuf_XZ;
//......................................................................................
double *sendData_x, *sendData_y, *sendData_z, *sendData_X, *sendData_Y, *sendData_Z;
double *sendData_xy, *sendData_yz, *sendData_xz, *sendData_Xy, *sendData_Yz, *sendData_xZ;
double *sendData_xY, *sendData_yZ, *sendData_Xz, *sendData_XY, *sendData_YZ, *sendData_XZ;
double *recvData_x, *recvData_y, *recvData_z, *recvData_X, *recvData_Y, *recvData_Z;
double *recvData_xy, *recvData_yz, *recvData_xz, *recvData_Xy, *recvData_Yz, *recvData_xZ;
double *recvData_xY, *recvData_yZ, *recvData_Xz, *recvData_XY, *recvData_YZ, *recvData_XZ;
};

1176
common/MPI.I Normal file
View File

File diff suppressed because it is too large Load Diff

3831
common/MPI.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

1166
common/MPI.h Normal file
View File

File diff suppressed because it is too large Load Diff

266
common/MPI_Helpers.cpp
View File

@@ -1,266 +0,0 @@
#include "common/MPI_Helpers.h"
#include "common/Utilities.h"
/********************************************************
* Return the MPI data type *
********************************************************/
template<> MPI_Datatype getMPItype<char>() {
return MPI_CHAR;
}
template<> MPI_Datatype getMPItype<unsigned char>() {
return MPI_UNSIGNED_CHAR;
}
template<> MPI_Datatype getMPItype<int>() {
return MPI_INT;
}
template<> MPI_Datatype getMPItype<long>() {
return MPI_LONG;
}
template<> MPI_Datatype getMPItype<unsigned long>() {
return MPI_UNSIGNED_LONG;
}
template<> MPI_Datatype getMPItype<long long>() {
return MPI_LONG_LONG;
}
template<> MPI_Datatype getMPItype<float>() {
return MPI_FLOAT;
}
template<> MPI_Datatype getMPItype<double>() {
return MPI_DOUBLE;
}
/********************************************************
* Concrete implimentations for packing/unpacking *
********************************************************/
// unsigned char
template<>
size_t packsize<unsigned char>( const unsigned char& )
{
return sizeof(unsigned char);
}
template<>
void pack<unsigned char>( const unsigned char& rhs, char *buffer )
{
memcpy(buffer,&rhs,sizeof(unsigned char));
}
template<>
void unpack<unsigned char>( unsigned char& data, const char *buffer )
{
memcpy(&data,buffer,sizeof(unsigned char));
}
// char
template<>
size_t packsize<char>( const char& )
{
return sizeof(char);
}
template<>
void pack<char>( const char& rhs, char *buffer )
{
memcpy(buffer,&rhs,sizeof(char));
}
template<>
void unpack<char>( char& data, const char *buffer )
{
memcpy(&data,buffer,sizeof(char));
}
// int
template<>
size_t packsize<int>( const int& )
{
return sizeof(int);
}
template<>
void pack<int>( const int& rhs, char *buffer )
{
memcpy(buffer,&rhs,sizeof(int));
}
template<>
void unpack<int>( int& data, const char *buffer )
{
memcpy(&data,buffer,sizeof(int));
}
// unsigned int
template<>
size_t packsize<unsigned int>( const unsigned int& )
{
return sizeof(unsigned int);
}
template<>
void pack<unsigned int>( const unsigned int& rhs, char *buffer )
{
memcpy(buffer,&rhs,sizeof(int));
}
template<>
void unpack<unsigned int>( unsigned int& data, const char *buffer )
{
memcpy(&data,buffer,sizeof(int));
}
// size_t
template<>
size_t packsize<size_t>( const size_t& )
{
return sizeof(size_t);
}
template<>
void pack<size_t>( const size_t& rhs, char *buffer )
{
memcpy(buffer,&rhs,sizeof(size_t));
}
template<>
void unpack<size_t>( size_t& data, const char *buffer )
{
memcpy(&data,buffer,sizeof(size_t));
}
// std::string
template<>
size_t packsize<std::string>( const std::string& rhs )
{
return rhs.size()+1;
}
template<>
void pack<std::string>( const std::string& rhs, char *buffer )
{
memcpy(buffer,rhs.c_str(),rhs.size()+1);
}
template<>
void unpack<std::string>( std::string& data, const char *buffer )
{
data = std::string(buffer);
}
/********************************************************
* Fake MPI routines *
********************************************************/
#ifndef USE_MPI
int MPI_Init(int*,char***)
{
return 0;
}
int MPI_Init_thread(int*,char***, int required, int *provided )
{
*provided = required;
return 0;
}
int MPI_Finalize()
{
return 0;
}
int MPI_Comm_size( MPI_Comm, int *size )
{
*size = 1;
return 0;
}
int MPI_Comm_rank( MPI_Comm, int *rank )
{
*rank = 0;
return 0;
}
int MPI_Barrier( MPI_Comm )
{
return 0;
}
int MPI_Waitall( int, MPI_Request[], MPI_Status[] )
{
return 0;
}
int MPI_Wait( MPI_Request*, MPI_Status* )
{
return 0;
}
int MPI_Bcast( void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm )
{
return 0;
}
int MPI_Send(const void *buf, int count, MPI_Datatype datatype, int dest, int tag,
MPI_Comm comm)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Recv(void *buf, int count, MPI_Datatype datatype, int source, int tag,
MPI_Comm comm, MPI_Status *status)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Isend(const void *buf, int count, MPI_Datatype datatype, int dest, int tag,
MPI_Comm comm, MPI_Request *request)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Irecv(void *buf, int count, MPI_Datatype datatype, int source,
int tag, MPI_Comm comm, MPI_Request *request)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count,
MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
void *recvbuf, int recvcount, MPI_Datatype recvtype,
MPI_Comm comm)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Allgatherv(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
void *recvbuf, const int *recvcounts, const int *displs,
MPI_Datatype recvtype, MPI_Comm comm)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Sendrecv(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
int dest, int sendtag,
void *recvbuf, int recvcount, MPI_Datatype recvtype,
int source, int recvtag,
MPI_Comm comm, MPI_Status *status)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Reduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype,
MPI_Op op, int root, MPI_Comm comm)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Comm_group(MPI_Comm comm, MPI_Group *group)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Comm_create(MPI_Comm comm, MPI_Group group, MPI_Comm *newcomm)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Comm_dup(MPI_Comm comm, MPI_Comm *newcomm)
{
*newcomm = comm;
return 0;
}
double MPI_Wtime( void )
{
return 0.0;
}
int MPI_Comm_free(MPI_Comm *group)
{
return 0;
}
int MPI_Group_free(MPI_Group *group)
{
return 0;
}
#endif

239
common/MPI_Helpers.h
View File

@@ -1,239 +0,0 @@
// This file contains wrappers for MPI routines and functions to pack/unpack data structures
#ifndef MPI_WRAPPERS_INC
#define MPI_WRAPPERS_INC
#include <string.h>
#include <vector>
#include <set>
#include <map>
#ifdef USE_MPI
// Inlcude MPI
#include "mpi.h"
#else
// Create fake MPI types
typedef int MPI_Comm;
typedef int MPI_Request;
typedef int MPI_Status;
#define MPI_COMM_WORLD 0
#define MPI_COMM_SELF 0
#define MPI_COMM_NULL -1
#define MPI_GROUP_NULL -2
#define MPI_STATUS_IGNORE NULL
enum MPI_Datatype { MPI_LOGICAL, MPI_CHAR, MPI_UNSIGNED_CHAR, MPI_INT,
MPI_UNSIGNED, MPI_LONG, MPI_UNSIGNED_LONG, MPI_LONG_LONG, MPI_FLOAT, MPI_DOUBLE };
enum MPI_Op { MPI_MIN, MPI_MAX, MPI_SUM };
typedef int MPI_Group;
#define MPI_THREAD_SINGLE 0
#define MPI_THREAD_FUNNELED 1
#define MPI_THREAD_SERIALIZED 2
#define MPI_THREAD_MULTIPLE 3
// Fake MPI functions
int MPI_Init(int*,char***);
int MPI_Init_thread( int *argc, char ***argv, int required, int *provided );
int MPI_Finalize();
int MPI_Comm_size( MPI_Comm, int *size );
int MPI_Comm_rank( MPI_Comm, int *rank );
int MPI_Barrier(MPI_Comm);
int MPI_Wait(MPI_Request*,MPI_Status*);
int MPI_Waitall(int,MPI_Request[],MPI_Status[]);
int MPI_Bcast(void*,int,MPI_Datatype,int,MPI_Comm);
int MPI_Send(const void *buf, int count, MPI_Datatype datatype, int dest, int tag,
MPI_Comm comm);
int MPI_Recv(void *buf, int count, MPI_Datatype datatype, int source, int tag,
MPI_Comm comm, MPI_Status *status);
int MPI_Isend(const void *buf, int count, MPI_Datatype datatype, int dest, int tag,
MPI_Comm comm, MPI_Request *request);
int MPI_Irecv(void *buf, int count, MPI_Datatype datatype, int source,
int tag, MPI_Comm comm, MPI_Request *request);
int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count,
MPI_Datatype datatype, MPI_Op op, MPI_Comm comm);
int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
void *recvbuf, int recvcount, MPI_Datatype recvtype,
MPI_Comm comm);
int MPI_Allgatherv(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
void *recvbuf, const int *recvcounts, const int *displs,
MPI_Datatype recvtype, MPI_Comm comm);
int MPI_Sendrecv(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
int dest, int sendtag,
void *recvbuf, int recvcount, MPI_Datatype recvtype,
int source, int recvtag,
MPI_Comm comm, MPI_Status *status);
int MPI_Reduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype,
MPI_Op op, int root, MPI_Comm comm);
double MPI_Wtime( void );
int MPI_Comm_group(MPI_Comm comm, MPI_Group *group);
int MPI_Comm_create(MPI_Comm comm, MPI_Group group, MPI_Comm *newcomm);
int MPI_Comm_free(MPI_Comm *group);
int MPI_Group_free(MPI_Group *group);
int MPI_Comm_dup(MPI_Comm comm, MPI_Comm *newcomm);
#endif
//! Get the size of the MPI_Comm
// Note: this is a thread and interrupt safe function
inline int comm_size( MPI_Comm comm ) {
int size = 1;
MPI_Comm_size( comm, &size );
return size;
}
//! Get the rank of the MPI_Comm
// Note: this is a thread and interrupt safe function
inline int comm_rank( MPI_Comm comm ) {
int rank = 1;
MPI_Comm_rank( comm, &rank );
return rank;
}
//! Get the size of MPI_COMM_WORLD
inline int MPI_WORLD_SIZE( ) {
return comm_size( MPI_COMM_WORLD );
}
//! Get the size of MPI_COMM_WORLD
inline int MPI_WORLD_RANK( ) {
return comm_rank( MPI_COMM_WORLD );
}
//! Return the appropriate MPI datatype for a class
template<class TYPE>
MPI_Datatype getMPItype();
//! Template function to return the buffer size required to pack a class
template<class TYPE>
size_t packsize( const TYPE& rhs );
//! Template function to pack a class to a buffer
template<class TYPE>
void pack( const TYPE& rhs, char *buffer );
//! Template function to unpack a class from a buffer
template<class TYPE>
void unpack( TYPE& data, const char *buffer );
//! Template function to return the buffer size required to pack a std::vector
template<class TYPE>
size_t packsize( const std::vector<TYPE>& rhs );
//! Template function to pack a class to a buffer
template<class TYPE>
void pack( const std::vector<TYPE>& rhs, char *buffer );
//! Template function to pack a class to a buffer
template<class TYPE>
void unpack( std::vector<TYPE>& data, const char *buffer );
//! Template function to return the buffer size required to pack a std::pair
template<class TYPE1, class TYPE2>
size_t packsize( const std::pair<TYPE1,TYPE2>& rhs );
//! Template function to pack a class to a buffer
template<class TYPE1, class TYPE2>
void pack( const std::pair<TYPE1,TYPE2>& rhs, char *buffer );
//! Template function to pack a class to a buffer
template<class TYPE1, class TYPE2>
void unpack( std::pair<TYPE1,TYPE2>& data, const char *buffer );
//! Template function to return the buffer size required to pack a std::map
template<class TYPE1, class TYPE2>
size_t packsize( const std::map<TYPE1,TYPE2>& rhs );
//! Template function to pack a class to a buffer
template<class TYPE1, class TYPE2>
void pack( const std::map<TYPE1,TYPE2>& rhs, char *buffer );
//! Template function to pack a class to a buffer
template<class TYPE1, class TYPE2>
void unpack( std::map<TYPE1,TYPE2>& data, const char *buffer );
//! Template function to return the buffer size required to pack a std::set
template<class TYPE>
size_t packsize( const std::set<TYPE>& rhs );
//! Template function to pack a class to a buffer
template<class TYPE>
void pack( const std::set<TYPE>& rhs, char *buffer );
//! Template function to pack a class to a buffer
template<class TYPE>
void unpack( std::set<TYPE>& data, const char *buffer );
// Helper functions
inline double sumReduce( MPI_Comm comm, double x )
{
double y = 0;
MPI_Allreduce(&x,&y,1,MPI_DOUBLE,MPI_SUM,comm);
return y;
}
inline float sumReduce( MPI_Comm comm, float x )
{
float y = 0;
MPI_Allreduce(&x,&y,1,MPI_FLOAT,MPI_SUM,comm);
return y;
}
inline int sumReduce( MPI_Comm comm, int x )
{
int y = 0;
MPI_Allreduce(&x,&y,1,MPI_INT,MPI_SUM,comm);
return y;
}
inline long long sumReduce( MPI_Comm comm, long long x )
{
long long y = 0;
MPI_Allreduce(&x,&y,1,MPI_LONG_LONG,MPI_SUM,comm);
return y;
}
inline bool sumReduce( MPI_Comm comm, bool x )
{
int y = sumReduce( comm, x?1:0 );
return y>0;
}
inline std::vector<float> sumReduce( MPI_Comm comm, const std::vector<float>& x )
{
auto y = x;
MPI_Allreduce(x.data(),y.data(),x.size(),MPI_FLOAT,MPI_SUM,comm);
return y;
}
inline std::vector<int> sumReduce( MPI_Comm comm, const std::vector<int>& x )
{
auto y = x;
MPI_Allreduce(x.data(),y.data(),x.size(),MPI_INT,MPI_SUM,comm);
return y;
}
inline double maxReduce( MPI_Comm comm, double x )
{
double y = 0;
MPI_Allreduce(&x,&y,1,MPI_DOUBLE,MPI_MAX,comm);
return y;
}
inline float maxReduce( MPI_Comm comm, float x )
{
float y = 0;
MPI_Allreduce(&x,&y,1,MPI_FLOAT,MPI_MAX,comm);
return y;
}
inline int maxReduce( MPI_Comm comm, int x )
{
int y = 0;
MPI_Allreduce(&x,&y,1,MPI_INT,MPI_MAX,comm);
return y;
}
#endif
#include "common/MPI_Helpers.hpp"

154
common/MPI_Helpers.hpp
View File

@@ -1,154 +0,0 @@
// This file contains wrappers for MPI routines and functions to pack/unpack data structures
#ifndef MPI_WRAPPERS_HPP
#define MPI_WRAPPERS_HPP
#include "common/MPI_Helpers.h"
#include <string.h>
#include <vector>
#include <set>
#include <map>
/********************************************************
* Default instantiations for std::vector *
********************************************************/
template<class TYPE>
size_t packsize( const std::vector<TYPE>& rhs )
{
size_t bytes = sizeof(size_t);
for (size_t i=0; i<rhs.size(); i++)
bytes += packsize(rhs[i]);
return bytes;
}
template<class TYPE>
void pack( const std::vector<TYPE>& rhs, char *buffer )
{
size_t size = rhs.size();
memcpy(buffer,&size,sizeof(size_t));
size_t pos = sizeof(size_t);
for (size_t i=0; i<rhs.size(); i++) {
pack(rhs[i],&buffer[pos]);
pos += packsize(rhs[i]);
}
}
template<class TYPE>
void unpack( std::vector<TYPE>& data, const char *buffer )
{
size_t size;
memcpy(&size,buffer,sizeof(size_t));
data.clear();
data.resize(size);
size_t pos = sizeof(size_t);
for (size_t i=0; i<data.size(); i++) {
unpack(data[i],&buffer[pos]);
pos += packsize(data[i]);
}
}
/********************************************************
* Default instantiations for std::pair *
********************************************************/
template<class TYPE1, class TYPE2>
size_t packsize( const std::pair<TYPE1,TYPE2>& rhs )
{
return packsize(rhs.first)+packsize(rhs.second);
}
template<class TYPE1, class TYPE2>
void pack( const std::pair<TYPE1,TYPE2>& rhs, char *buffer )
{
pack(rhs.first,buffer);
pack(rhs.second,&buffer[packsize(rhs.first)]);
}
template<class TYPE1, class TYPE2>
void unpack( std::pair<TYPE1,TYPE2>& data, const char *buffer )
{
unpack(data.first,buffer);
unpack(data.second,&buffer[packsize(data.first)]);
}
/********************************************************
* Default instantiations for std::map *
********************************************************/
template<class TYPE1, class TYPE2>
size_t packsize( const std::map<TYPE1,TYPE2>& rhs )
{
size_t bytes = sizeof(size_t);
typename std::map<TYPE1,TYPE2>::const_iterator it;
for (it=rhs.begin(); it!=rhs.end(); ++it) {
bytes += packsize(it->first);
bytes += packsize(it->second);
}
return bytes;
}
template<class TYPE1, class TYPE2>
void pack( const std::map<TYPE1,TYPE2>& rhs, char *buffer )
{
size_t N = rhs.size();
pack(N,buffer);
size_t pos = sizeof(size_t);
typename std::map<TYPE1,TYPE2>::const_iterator it;
for (it=rhs.begin(); it!=rhs.end(); ++it) {
pack(it->first,&buffer[pos]); pos+=packsize(it->first);
pack(it->second,&buffer[pos]); pos+=packsize(it->second);
}
}
template<class TYPE1, class TYPE2>
void unpack( std::map<TYPE1,TYPE2>& data, const char *buffer )
{
size_t N = 0;
unpack(N,buffer);
size_t pos = sizeof(size_t);
data.clear();
for (size_t i=0; i<N; i++) {
std::pair<TYPE1,TYPE2> tmp;
unpack(tmp.first,&buffer[pos]); pos+=packsize(tmp.first);
unpack(tmp.second,&buffer[pos]); pos+=packsize(tmp.second);
data.insert(tmp);
}
}
/********************************************************
* Default instantiations for std::set *
********************************************************/
template<class TYPE>
size_t packsize( const std::set<TYPE>& rhs )
{
size_t bytes = sizeof(size_t);
typename std::set<TYPE>::const_iterator it;
for (it=rhs.begin(); it!=rhs.end(); ++it) {
bytes += packsize(*it);
}
return bytes;
}
template<class TYPE>
void pack( const std::set<TYPE>& rhs, char *buffer )
{
size_t N = rhs.size();
pack(N,buffer);
size_t pos = sizeof(size_t);
typename std::set<TYPE>::const_iterator it;
for (it=rhs.begin(); it!=rhs.end(); ++it) {
pack(*it); pos+=packsize(*it);
}
}
template<class TYPE>
void unpack( std::set<TYPE>& data, const char *buffer )
{
size_t N = 0;
unpack(N,buffer);
size_t pos = sizeof(size_t);
data.clear();
for (size_t i=0; i<N; i++) {
TYPE tmp;
unpack(tmp,&buffer[pos]); pos+=packsize(tmp);
data.insert(tmp);
}
}
#endif

4
common/ReadMicroCT.cpp
View File

@@ -64,11 +64,11 @@ Array<uint8_t> readMicroCT( const std::string& filename )
// Read the compressed micro CT data and distribute
Array<uint8_t> readMicroCT( const Database& domain, MPI_Comm comm )
Array<uint8_t> readMicroCT( const Database& domain, const Utilities::MPI& comm )
{
// Get the local problem info
auto n = domain.getVector<int>( "n" );
int rank = comm_rank(MPI_COMM_WORLD);
int rank = comm.getRank();
auto nproc = domain.getVector<int>( "nproc" );
RankInfoStruct rankInfo( rank, nproc[0], nproc[1], nproc[2] );

3
common/ReadMicroCT.h
View File

@@ -5,11 +5,12 @@
#include "common/Array.h"
#include "common/Communication.h"
#include "common/Database.h"
#include "common/MPI.h"
Array<uint8_t> readMicroCT( const std::string& filename );
Array<uint8_t> readMicroCT( const Database& domain, MPI_Comm comm );
Array<uint8_t> readMicroCT( const Database& domain, const Utilities::MPI& comm );
#endif

916
common/ScaLBL.cpp
View File

File diff suppressed because it is too large Load Diff

107
common/ScaLBL.h
View File

@@ -70,6 +70,49 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Fx, double Fy, double Fz,
double *Poros,double *Perm, double *Velocity,double Den,double *Pressure);
// ION TRANSPORT MODEL
extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *dist, double *Den, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit, int Np);
extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np);
extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np);
// LBM Poisson solver
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,
int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,
int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(int *neighborList,int *Map, double *dist, double *Psi, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *dist, double *Psi, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, int start, int finish, int Np);
//maybe deprecated
//extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
// int strideY, int strideZ,int start, int finish, int Np);
// LBM Stokes Model (adapted from MRT model)
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_MRT(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Fx, double Fy, double Fz,
double *Poros,double *Perm, double *Velocity,double Den,double *Pressure);
@@ -78,12 +121,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist
double *Poros,double *Perm, double *Velocity,double Den,double *Pressure);
extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel,
double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel,double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel,
double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel,double *Pressure,
double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
// ION TRANSPORT MODEL
@@ -130,6 +173,8 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, do
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np);
extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, double *Bq, int start, int finish, int Np);
// MRT MODEL
extern "C" void ScaLBL_D3Q19_AAeven_MRT(double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx,
double Fy, double Fz);
@@ -231,6 +276,14 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion_Concentration_BC_z(int *d_neighborList, in
extern "C" void ScaLBL_D3Q7_AAodd_Ion_Concentration_BC_Z(int *d_neighborList, int *list, double *dist, double Cout, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Ion_Flux_BC_z(int *list, double *dist, double Cin, double tau, double *VelocityZ, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Ion_Flux_BC_Z(int *list, double *dist, double Cout, double tau, double *VelocityZ, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Ion_Flux_BC_z(int *d_neighborList, int *list, double *dist, double Cin, double tau, double *VelocityZ, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Ion_Flux_BC_Z(int *d_neighborList, int *list, double *dist, double Cout, double tau, double *VelocityZ, int count, int Np);
class ScaLBL_Communicator{
public:
//......................................................................................
@@ -239,7 +292,6 @@ public:
//ScaLBL_Communicator(Domain &Dm, IntArray &Map);
~ScaLBL_Communicator();
//......................................................................................
MPI_Comm MPI_COMM_SCALBL; // MPI Communicator
unsigned long int CommunicationCount,SendCount,RecvCount;
int Nx,Ny,Nz,N;
int n_bb_d3q7, n_bb_d3q19;
@@ -266,6 +318,10 @@ public:
int LastInterior();
int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np, int width);
void Barrier(){
ScaLBL_DeviceBarrier();
MPI_COMM_SCALBL.barrier();
};
void SendD3Q19AA(double *dist);
void RecvD3Q19AA(double *dist);
void SendD3Q7AA(double *fq, int Component);
@@ -282,31 +338,31 @@ public:
void SolidDirichletD3Q7(double *fq, double *BoundaryValue);
void SolidNeumannD3Q7(double *fq, double *BoundaryValue);
// Routines to set boundary conditions
void Color_BC_z(int *Map, double *Phi, double *Den, double vA, double vB);
void Color_BC_Z(int *Map, double *Phi, double *Den, double vA, double vB);
void D3Q19_Pressure_BC_z(int *neighborList, double *fq, double din, int time);
void D3Q19_Pressure_BC_Z(int *neighborList, double *fq, double dout, int time);
void D3Q19_Reflection_BC_z(double *fq);
void D3Q19_Reflection_BC_Z(double *fq);
double D3Q19_Flux_BC_z(int *neighborList, double *fq, double flux, int time);
void GreyscaleSC_BC_z(int *Map, double *DenA, double *DenB, double vA, double vB);
void GreyscaleSC_BC_Z(int *Map, double *DenA, double *DenB, double vA, double vB);
// Routines to set boundary conditions
void Color_BC_z(int *Map, double *Phi, double *Den, double vA, double vB);
void Color_BC_Z(int *Map, double *Phi, double *Den, double vA, double vB);
void D3Q19_Pressure_BC_z(int *neighborList, double *fq, double din, int time);
void D3Q19_Pressure_BC_Z(int *neighborList, double *fq, double dout, int time);
void D3Q19_Reflection_BC_z(double *fq);
void D3Q19_Reflection_BC_Z(double *fq);
double D3Q19_Flux_BC_z(int *neighborList, double *fq, double flux, int time);
void D3Q7_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time);
void D3Q7_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time);
void Poisson_D3Q7_BC_z(int *Map, double *Psi, double Vin);
void Poisson_D3Q7_BC_Z(int *Map, double *Psi, double Vout);
void D3Q7_Ion_Concentration_BC_z(int *neighborList, double *fq, double Cin, int time);
void D3Q7_Ion_Concentration_BC_Z(int *neighborList, double *fq, double Cout, int time);
void D3Q7_Ion_Flux_BC_z(int *neighborList, double *fq, double Cin, double tau, double *VelocityZ, int time);
void D3Q7_Ion_Flux_BC_Z(int *neighborList, double *fq, double Cout, double tau, double *VelocityZ, int time);
void GreyscaleSC_BC_z(int *Map, double *DenA, double *DenB, double vA, double vB);
void GreyscaleSC_BC_Z(int *Map, double *DenA, double *DenB, double vA, double vB);
void GreyscaleSC_Pressure_BC_z(int *neighborList, double *fqA, double *fqB, double dinA, double dinB, int time);
void GreyscaleSC_Pressure_BC_Z(int *neighborList, double *fqA, double *fqB, double doutA, double doutB, int time);
void D3Q7_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time);
void D3Q7_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time);
void Poisson_D3Q7_BC_z(int *Map, double *Psi, double Vin);
void Poisson_D3Q7_BC_Z(int *Map, double *Psi, double Vout);
void D3Q7_Ion_Concentration_BC_z(int *neighborList, double *fq, double Cin, int time);
void D3Q7_Ion_Concentration_BC_Z(int *neighborList, double *fq, double Cout, int time);
// Debugging and unit testing functions
void PrintD3Q19();
// Debugging and unit testing functions
void PrintD3Q19();
private:
//void D3Q19_MapRecv_OLD(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count, int *d3q19_recvlist);
void D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, int *list, int start, int count, int *d3q19_recvlist);
void D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, const int *list, int start, int count, int *d3q19_recvlist);
bool Lock; // use Lock to make sure only one call at a time to protect data in transit
// only one set of Send requests can be active at any time (per instance)
@@ -317,9 +373,8 @@ private:
int sendtag,recvtag;
// Give the object it's own MPI communicator
RankInfoStruct rank_info;
MPI_Group Group; // Group of processors associated with this domain
Utilities::MPI MPI_COMM_SCALBL; // MPI Communicator for this domain
MPI_Request req1[18],req2[18];
MPI_Status stat1[18],stat2[18];
//......................................................................................
// MPI ranks for all 18 neighbors
//......................................................................................

2
common/SpherePack.cpp
View File

@@ -9,7 +9,7 @@
#include "common/Array.h"
#include "common/Utilities.h"
#include "common/MPI_Helpers.h"
#include "common/MPI.h"
#include "common/Communication.h"
#include "common/Database.h"
#include "common/SpherePack.h"

2
common/SpherePack.h
View File

@@ -12,7 +12,7 @@
#include "common/Array.h"
#include "common/Utilities.h"
#include "common/MPI_Helpers.h"
#include "common/MPI.h"
#include "common/Communication.h"
#include "common/Database.h"

24
common/Utilities.cpp
View File

@@ -8,7 +8,7 @@
#endif
#ifdef USE_MPI
#include "mpi.h"
#include "common/MPI.h"
#endif
#include <algorithm>
@@ -37,7 +37,7 @@ static std::mutex Utilities_mutex;
/****************************************************************************
* Function to perform the default startup/shutdown sequences *
****************************************************************************/
void Utilities::startup( int argc, char **argv )
void Utilities::startup( int argc, char **argv, bool multiple )
{
NULL_USE( argc );
NULL_USE( argv );
@@ -46,15 +46,19 @@ void Utilities::startup( int argc, char **argv )
Utilities::setenv( "MKL_NUM_THREADS", "1" );
// Start MPI
#ifdef USE_MPI
int provided;
MPI_Init_thread( &argc, &argv, MPI_THREAD_MULTIPLE, &provided );
if ( provided < MPI_THREAD_MULTIPLE ) {
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
if ( rank == 0 )
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
if ( multiple ) {
int provided;
MPI_Init_thread( &argc, &argv, MPI_THREAD_MULTIPLE, &provided );
if ( provided < MPI_THREAD_MULTIPLE ) {
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
if ( rank == 0 )
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
}
StackTrace::globalCallStackInitialize( MPI_COMM_WORLD );
} else {
MPI_Init( &argc, &argv );
}
StackTrace::globalCallStackInitialize( MPI_COMM_WORLD );
#endif
// Set the error handlers
Utilities::setAbortBehavior( true, 3 );

2
common/Utilities.h
View File

@@ -31,7 +31,7 @@ using StackTrace::Utilities::sleep_s;
* \param argc argc from main
* \param argv argv from main
*/
void startup( int argc, char **argv );
void startup( int argc, char **argv, bool multiple=true );
/*!
* \brief Stop MPI, error handlers

234
common/Utilities.hpp Normal file
View File

@@ -0,0 +1,234 @@
#ifndef included_Utilities_hpp
#define included_Utilities_hpp
#include "Utilities.h"
#include <vector>
namespace Utilities {
/************************************************************************
* templated quicksort routines *
************************************************************************/
template<class T>
void quicksort( std::vector<T> &x )
{
if ( x.size() <= 1u )
return;
T *arr = &x[0];
bool test;
long int i, ir, j, jstack, k, l, istack[100];
T a, tmp_a;
jstack = 0;
l = 0;
ir = x.size() - 1;
while ( 1 ) {
if ( ir - l < 7 ) { // Insertion sort when subarray small enough.
for ( j = l + 1; j <= ir; j++ ) {
a = arr[j];
test = true;
for ( i = j - 1; i >= 0; i-- ) {
if ( arr[i] < a ) {
arr[i + 1] = a;
test = false;
break;
}
arr[i + 1] = arr[i];
}
if ( test ) {
i = l - 1;
arr[i + 1] = a;
}
}
if ( jstack == 0 )
return;
ir = istack[jstack]; // Pop stack and begin a new round of partitioning.
l = istack[jstack - 1];
jstack -= 2;
} else {
k = ( l + ir ) / 2; // Choose median of left, center and right elements as partitioning
// element a. Also rearrange so that a(l) < a(l+1) < a(ir).
tmp_a = arr[k];
arr[k] = arr[l + 1];
arr[l + 1] = tmp_a;
if ( arr[l] > arr[ir] ) {
tmp_a = arr[l];
arr[l] = arr[ir];
arr[ir] = tmp_a;
}
if ( arr[l + 1] > arr[ir] ) {
tmp_a = arr[l + 1];
arr[l + 1] = arr[ir];
arr[ir] = tmp_a;
}
if ( arr[l] > arr[l + 1] ) {
tmp_a = arr[l];
arr[l] = arr[l + 1];
arr[l + 1] = tmp_a;
}
// Scan up to find element > a
j = ir;
a = arr[l + 1]; // Partitioning element.
for ( i = l + 2; i <= ir; i++ ) {
if ( arr[i] < a )
continue;
while ( arr[j] > a ) // Scan down to find element < a.
j--;
if ( j < i )
break; // Pointers crossed. Exit with partitioning complete.
tmp_a = arr[i]; // Exchange elements of both arrays.
arr[i] = arr[j];
arr[j] = tmp_a;
}
arr[l + 1] = arr[j]; // Insert partitioning element in both arrays.
arr[j] = a;
jstack += 2;
// Push pointers to larger subarray on stack, process smaller subarray immediately.
if ( ir - i + 1 >= j - l ) {
istack[jstack] = ir;
istack[jstack - 1] = i;
ir = j - 1;
} else {
istack[jstack] = j - 1;
istack[jstack - 1] = l;
l = i;
}
}
}
}
template<class T1, class T2>
void quicksort( std::vector<T1> &x, std::vector<T2> &y )
{
if ( x.size() <= 1u )
return;
T1 *arr = &x[0];
T2 *brr = &y[0];
bool test;
long int i, ir, j, jstack, k, l, istack[100];
T1 a, tmp_a;
T2 b, tmp_b;
jstack = 0;
l = 0;
ir = x.size() - 1;
while ( 1 ) {
if ( ir - l < 7 ) { // Insertion sort when subarray small enough.
for ( j = l + 1; j <= ir; j++ ) {
a = arr[j];
b = brr[j];
test = true;
for ( i = j - 1; i >= 0; i-- ) {
if ( arr[i] < a ) {
arr[i + 1] = a;
brr[i + 1] = b;
test = false;
break;
}
arr[i + 1] = arr[i];
brr[i + 1] = brr[i];
}
if ( test ) {
i = l - 1;
arr[i + 1] = a;
brr[i + 1] = b;
}
}
if ( jstack == 0 )
return;
ir = istack[jstack]; // Pop stack and begin a new round of partitioning.
l = istack[jstack - 1];
jstack -= 2;
} else {
k = ( l + ir ) / 2; // Choose median of left, center and right elements as partitioning
// element a. Also rearrange so that a(l) ? a(l+1) ? a(ir).
tmp_a = arr[k];
arr[k] = arr[l + 1];
arr[l + 1] = tmp_a;
tmp_b = brr[k];
brr[k] = brr[l + 1];
brr[l + 1] = tmp_b;
if ( arr[l] > arr[ir] ) {
tmp_a = arr[l];
arr[l] = arr[ir];
arr[ir] = tmp_a;
tmp_b = brr[l];
brr[l] = brr[ir];
brr[ir] = tmp_b;
}
if ( arr[l + 1] > arr[ir] ) {
tmp_a = arr[l + 1];
arr[l + 1] = arr[ir];
arr[ir] = tmp_a;
tmp_b = brr[l + 1];
brr[l + 1] = brr[ir];
brr[ir] = tmp_b;
}
if ( arr[l] > arr[l + 1] ) {
tmp_a = arr[l];
arr[l] = arr[l + 1];
arr[l + 1] = tmp_a;
tmp_b = brr[l];
brr[l] = brr[l + 1];
brr[l + 1] = tmp_b;
}
// Scan up to find element > a
j = ir;
a = arr[l + 1]; // Partitioning element.
b = brr[l + 1];
for ( i = l + 2; i <= ir; i++ ) {
if ( arr[i] < a )
continue;
while ( arr[j] > a ) // Scan down to find element < a.
j--;
if ( j < i )
break; // Pointers crossed. Exit with partitioning complete.
tmp_a = arr[i]; // Exchange elements of both arrays.
arr[i] = arr[j];
arr[j] = tmp_a;
tmp_b = brr[i];
brr[i] = brr[j];
brr[j] = tmp_b;
}
arr[l + 1] = arr[j]; // Insert partitioning element in both arrays.
arr[j] = a;
brr[l + 1] = brr[j];
brr[j] = b;
jstack += 2;
// Push pointers to larger subarray on stack, process smaller subarray immediately.
if ( ir - i + 1 >= j - l ) {
istack[jstack] = ir;
istack[jstack - 1] = i;
ir = j - 1;
} else {
istack[jstack] = j - 1;
istack[jstack - 1] = l;
l = i;
}
}
}
}
template<class T>
void unique( std::vector<T> &x )
{
if ( x.size() <= 1 )
return;
// First perform a quicksort
quicksort( x );
// Next remove duplicate entries
size_t pos = 1;
for ( size_t i = 1; i < x.size(); i++ ) {
if ( x[i] != x[pos - 1] ) {
x[pos] = x[i];
pos++;
}
}
if ( pos < x.size() )
x.resize( pos );
}
}
#endif