easy-iso

Iso-surface extraction from volume data
git clone https://git.0xfab.ch/easy-iso.git
Log | Files | Refs | Submodules | README | LICENSE

commit a620016e1d48a7e2446fdfea028cee5c0368c69e
parent bd3523e7ea4204278a69d2bea376074868e6f2da
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date:   Thu,  1 Dec 2016 17:14:14 +0100

only use cartesian communicator

Diffstat:
Minclude/common.h | 53++++++++++++++++++++++++++++++++++++++++++++++++-----
Msrc/InterpolatorMPI.h | 130++++++++++++++++++++++++++++++++++++++++++++-----------------------------------
Msrc/IsoExtractorMPI.h | 2+-
3 files changed, 121 insertions(+), 64 deletions(-)

diff --git a/include/common.h b/include/common.h @@ -33,16 +33,28 @@ template <typename T, size_t P> struct Face { T* data; - const int _N; - const int _M; - const int _Nelem; + int _N; + int _M; + int _Nelem; static constexpr int _P = P; + Face() : data(nullptr) {} Face(const size_t N, const size_t M) : data(new T[N*M*P]), _N(N), _M(M), _Nelem(N*M*P) {} ~Face() { delete [] data; } + void alloc(const size_t N, const size_t M) + { + if (data == nullptr) + { + _N = N; + _M = M; + _Nelem = N*M*P; + data = new T[_Nelem]; + } + } + inline T& operator()(const size_t ix, const size_t iy, const size_t iz) { assert(ix < _N); @@ -65,16 +77,27 @@ template <typename T, size_t P, size_t Q> struct Edge { T* data; - const int _N; - const int _Nelem; + int _N; + int _Nelem; static constexpr int _P = P; static constexpr int _Q = Q; + Edge() : data(nullptr) {} Edge(const size_t N) : data(new T[N*P*Q]), _N(N), _Nelem(N*P*Q) {} ~Edge() { delete [] data; } + void alloc(const size_t N) + { + if (data == nullptr) + { + _N = N; + _Nelem = N*P*Q; + data = new T[_Nelem]; + } + } + inline T& operator()(const size_t ix, const size_t iy, const size_t iz) { assert(ix < _N); @@ -133,11 +156,23 @@ struct Faces // indices 0 and 1 correspond to the near and far faces relative to the // origin. Face_t x0, x1, y0, y1, z0, z1; + Faces() {} Faces(const size_t XN, const size_t XM, const size_t YN, const size_t YM, const size_t ZN, const size_t ZM) : x0(XN,XM), x1(XN,XM), y0(YN,YM), y1(YN,YM), z0(ZN,ZM), z1(ZN,ZM) {} + void alloc(const size_t XN, const size_t XM, + const size_t YN, const size_t YM, + const size_t ZN, const size_t ZM) + { + x0.alloc(XN,XM); + x1.alloc(XN,XM); + y0.alloc(YN,YM); + y1.alloc(YN,YM); + z0.alloc(ZN,ZM); + z1.alloc(ZN,ZM); + } }; struct Edges @@ -146,11 +181,19 @@ struct Edges // using right hand rule for the indicated direction. The 0-edge is the // one that goes through the origin. Edge_t x0, x1, x2, x3, y0, y1, y2, y3, z0, z1, z2, z3; + Edges() {} Edges(const size_t XN, const size_t YN, const size_t ZN) : x0(XN), x1(XN), x2(XN), x3(XN), y0(YN), y1(YN), y2(YN), y3(YN), z0(ZN), z1(ZN), z2(ZN), z3(ZN) {} + void alloc(const size_t XN, const size_t YN, const size_t ZN) + { + x0.alloc(XN); x1.alloc(XN); x2.alloc(XN); x3.alloc(XN); + y0.alloc(YN); y1.alloc(YN); y2.alloc(YN); y3.alloc(YN); + z0.alloc(ZN); z1.alloc(ZN); z2.alloc(ZN); z3.alloc(ZN); + } + }; struct Corners diff --git a/src/InterpolatorMPI.h b/src/InterpolatorMPI.h @@ -24,6 +24,7 @@ class MPInbr public: MPInbr(const MPI_Comm comm) { + // assumes comm is a communicator created with MPI_Cart_create int myidx[3]; int dims[3]; int periods[3]; @@ -58,28 +59,44 @@ private: int m_nbr[27]; }; +struct CubeHalo +{ + Faces faces; + Edges edges; + Corners corners; + + CubeHalo() {} + CubeHalo(const int Nx, const int Ny, const int Nz) : + faces(Ny,Nz, Nx,Nz, Nx,Ny), + edges(Nx, Ny, Nz) + {} + void alloc(const int Nx, const int Ny, const int Nz) + { + faces.alloc(Ny,Nz, Nx,Nz, Nx,Ny); + edges.alloc(Nx, Ny, Nz); + } +}; + template <typename Tinterp> class InterpolatorMPI : public Interpolator<Tinterp> { public: InterpolatorMPI(ArgumentParser& p, const MPI_Comm comm=MPI_COMM_WORLD) : Interpolator<Tinterp>(), - m_comm(comm), m_PESize{p("xpesize").asInt(1),p("ypesize").asInt(1),p("zpesize").asInt(1)} { this->m_extent = p("extent").asDouble(1.0); int periodic[3] = {false}; + MPI_Cart_create(comm, 3, m_PESize, periodic, true, &m_cartcomm); + int world_size; - MPI_Comm_size(comm, &world_size); + MPI_Comm_size(m_cartcomm, &world_size); assert(m_PESize[0]*m_PESize[1]*m_PESize[2] == world_size); - MPI_Comm cartcomm; - MPI_Cart_create(comm, 3, m_PESize, periodic, true, &cartcomm); - int myrank; - MPI_Comm_rank(cartcomm, &myrank); - MPI_Cart_coords(cartcomm, myrank, 3, m_myPEIndex); + MPI_Comm_rank(m_cartcomm, &myrank); + MPI_Cart_coords(m_cartcomm, myrank, 3, m_myPEIndex); m_isroot = (0 == myrank); // load the data @@ -141,8 +158,11 @@ public: this->m_origin_off[i] = origin[i]; // fetch halo cells + m_sendbuf.alloc(this->m_data.Nx(), this->m_data.Ny(), this->m_data.Nz()); + m_recvbuf.alloc(this->m_data.Nx(), this->m_data.Ny(), this->m_data.Nz()); + start = std::chrono::system_clock::now(); - _fetch_halo(cartcomm); + _fetch_halo(); end = std::chrono::system_clock::now(); elapsed_seconds = end-start; if (m_isroot) @@ -156,14 +176,14 @@ public: inline int getLocalNx() const { return this->m_data.Nx(); } inline int getLocalNy() const { return this->m_data.Ny(); } inline int getLocalNz() const { return this->m_data.Nz(); } - inline MPI_Comm getComm() const { return m_comm; } + inline MPI_Comm getComm() const { return m_cartcomm; } inline void getPESize(int ps[3]) const { for (int i = 0; i < 3; ++i) ps[i] = m_PESize[i]; } inline void getPEIndex(int pi[3]) const { for (int i = 0; i < 3; ++i) pi[i] = m_myPEIndex[i]; } inline bool isroot() const {return m_isroot; } private: - const MPI_Comm m_comm; + MPI_Comm m_cartcomm; int m_PESize[3]; int m_myPEIndex[3]; bool m_isroot; @@ -174,10 +194,12 @@ private: // halo handling using func_t = std::function<void(const int, const int, const int)>; - void _fetch_halo(const MPI_Comm comm); + CubeHalo m_sendbuf; + CubeHalo m_recvbuf; + void _fetch_halo(); - void _send_faces(const MPI_Comm comm, const MPInbr& nbrs, const Faces& f, std::vector<MPI_Request>& f_req) const; - void _recv_faces(const MPI_Comm comm, const MPInbr& nbrs, Faces& f, std::vector<MPI_Request>& f_req) const; + void _send_faces(const MPInbr& nbrs, const Faces& f, std::vector<MPI_Request>& f_req) const; + void _recv_faces(const MPInbr& nbrs, Faces& f, std::vector<MPI_Request>& f_req) const; void _pack_faces(Faces& f, const MPInbr& nbr); void _unpack_faces(const Faces& f, const MPInbr& nbr); inline void _process_faceX(func_t& eval); @@ -193,8 +215,8 @@ private: inline void _halo_faceY(const int _s0, const int _d0); inline void _halo_faceZ(const int _s0, const int _d0); - void _send_edges(const MPI_Comm comm, const MPInbr& nbrs, const Edges& e, std::vector<MPI_Request>& e_req) const; - void _recv_edges(const MPI_Comm comm, const MPInbr& nbrs, Edges& e, std::vector<MPI_Request>& e_req) const; + void _send_edges(const MPInbr& nbrs, const Edges& e, std::vector<MPI_Request>& e_req) const; + void _recv_edges(const MPInbr& nbrs, Edges& e, std::vector<MPI_Request>& e_req) const; void _pack_edges(Edges& f, const MPInbr& nbr); void _unpack_edges(const Edges& e, const MPInbr& nbr); inline void _process_edgeX(func_t& eval); @@ -210,8 +232,8 @@ private: inline void _halo_edgeY(const int _s0, const int _s1, const int _d0, const int _d1); inline void _halo_edgeZ(const int _s0, const int _s1, const int _d0, const int _d1); - void _send_corners(const MPI_Comm comm, const MPInbr& nbrs, const Corners& c, std::vector<MPI_Request>& c_req) const; - void _recv_corners(const MPI_Comm comm, const MPInbr& nbrs, Corners& c, std::vector<MPI_Request>& c_req) const; + void _send_corners(const MPInbr& nbrs, const Corners& c, std::vector<MPI_Request>& c_req) const; + void _recv_corners(const MPInbr& nbrs, Corners& c, std::vector<MPI_Request>& c_req) const; void _pack_corners(Corners& c, const MPInbr& nbr); void _unpack_corners(const Corners& c, const MPInbr& nbr); inline void _process_corner(func_t& eval); @@ -254,7 +276,7 @@ void InterpolatorMPI<Tinterp>::_load_h5_MPI(ArgumentParser& parser) maxDim[0] = dims[2]; maxDim[3] = dims[3]; } - MPI_Bcast(maxDim, 4, MPI_INT, 0, m_comm); + MPI_Bcast(maxDim, 4, MPI_INT, 0, m_cartcomm); NCH = maxDim[3]; if (maxDim[0] % m_PESize[0] != 0) @@ -288,7 +310,7 @@ void InterpolatorMPI<Tinterp>::_load_h5_MPI(ArgumentParser& parser) H5open(); fapl_id = H5Pcreate(H5P_FILE_ACCESS); - status = H5Pset_fapl_mpio(fapl_id, m_comm, MPI_INFO_NULL); + status = H5Pset_fapl_mpio(fapl_id, m_cartcomm, MPI_INFO_NULL); file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, fapl_id); status = H5Pclose(fapl_id); @@ -346,7 +368,7 @@ void InterpolatorMPI<Tinterp>::_load_wavelet_MPI(ArgumentParser& p) abort(); } - // Reader_WaveletCompressionMPI myreader(m_comm, filename, byte_swap, wavelet_type); + // Reader_WaveletCompressionMPI myreader(m_cartcomm, filename, byte_swap, wavelet_type); Reader_WaveletCompressionMPI myreader(MPI::COMM_WORLD, filename, byte_swap, wavelet_type); myreader.load_file(); const int NBX = myreader.xblocks(); @@ -397,7 +419,7 @@ void InterpolatorMPI<Tinterp>::_load_wavelet_MPI(ArgumentParser& p) template <typename Tinterp> -void InterpolatorMPI<Tinterp>::_fetch_halo(const MPI_Comm comm) +void InterpolatorMPI<Tinterp>::_fetch_halo() { // concept: // 0.) prepare neighbor info @@ -405,7 +427,7 @@ void InterpolatorMPI<Tinterp>::_fetch_halo(const MPI_Comm comm) // 2.) unpack buffers // 0.) - MPInbr mynbrs(comm); + MPInbr mynbrs(m_cartcomm); // 1.) const int Nx = this->m_data.Nx(); @@ -414,53 +436,45 @@ void InterpolatorMPI<Tinterp>::_fetch_halo(const MPI_Comm comm) // faces std::vector<MPI_Request> fpend_recv; - Faces frecv(Ny,Nz, Nx,Nz, Nx,Ny); - _recv_faces(comm, mynbrs, frecv, fpend_recv); + _recv_faces(mynbrs, m_recvbuf.faces, fpend_recv); std::vector<MPI_Request> fpend_send; - Faces fsend(Ny,Nz, Nx,Nz, Nx,Ny); - _pack_faces(fsend, mynbrs); - _send_faces(comm, mynbrs, fsend, fpend_send); + _pack_faces(m_sendbuf.faces, mynbrs); + _send_faces(mynbrs, m_sendbuf.faces, fpend_send); // edges std::vector<MPI_Request> epend_recv; - Edges erecv(Nx, Ny, Nz); - _recv_edges(comm, mynbrs, erecv, epend_recv); + _recv_edges(mynbrs, m_recvbuf.edges, epend_recv); std::vector<MPI_Request> epend_send; - Edges esend(Nx, Ny, Nz); - _pack_edges(esend, mynbrs); - _send_edges(comm, mynbrs, esend, epend_send); + _pack_edges(m_sendbuf.edges, mynbrs); + _send_edges(mynbrs, m_sendbuf.edges, epend_send); // corners std::vector<MPI_Request> cpend_recv; - Corners crecv; - _recv_corners(comm, mynbrs, crecv, cpend_recv); + _recv_corners(mynbrs, m_recvbuf.corners, cpend_recv); std::vector<MPI_Request> cpend_send; - Corners csend; - _pack_corners(csend, mynbrs); - _send_corners(comm, mynbrs, csend, cpend_send); + _pack_corners(m_sendbuf.corners, mynbrs); + _send_corners(mynbrs, m_sendbuf.corners, cpend_send); // 2.) std::vector<MPI_Status> fstat_recv(fpend_recv.size()); MPI_Waitall(fpend_recv.size(), fpend_recv.data(), fstat_recv.data()); - _unpack_faces(frecv, mynbrs); + _unpack_faces(m_recvbuf.faces, mynbrs); std::vector<MPI_Status> estat_recv(epend_recv.size()); MPI_Waitall(epend_recv.size(), epend_recv.data(), estat_recv.data()); - _unpack_edges(erecv, mynbrs); + _unpack_edges(m_recvbuf.edges, mynbrs); std::vector<MPI_Status> cstat_recv(cpend_recv.size()); MPI_Waitall(cpend_recv.size(), cpend_recv.data(), cstat_recv.data()); - _unpack_corners(crecv, mynbrs); - - MPI_Barrier(comm); + _unpack_corners(m_recvbuf.corners, mynbrs); } template <typename Tinterp> -void InterpolatorMPI<Tinterp>::_send_faces(const MPI_Comm comm, const MPInbr& nbrs, const Faces& f, std::vector<MPI_Request>& f_req) const +void InterpolatorMPI<Tinterp>::_send_faces(const MPInbr& nbrs, const Faces& f, std::vector<MPI_Request>& f_req) const { const int dst[6] = {nbrs(-1,0,0), nbrs(1,0,0), nbrs(0,-1,0), nbrs(0,1,0), nbrs(0,0,-1), nbrs(0,0,1)}; Real* const data[6] = {f.x0.data, f.x1.data, f.y0.data, f.y1.data, f.z0.data, f.z1.data}; @@ -473,7 +487,7 @@ void InterpolatorMPI<Tinterp>::_send_faces(const MPI_Comm comm, const MPInbr& nb if (dst[i] >= 0) { MPI_Request req; - MPI_Isend(data[i], Nelem[i], MPIREAL, dst[i], tagbase+dir*2+side, comm, &req); + MPI_Isend(data[i], Nelem[i], MPIREAL, dst[i], tagbase+dir*2+side, m_cartcomm, &req); f_req.push_back(req); } } @@ -481,7 +495,7 @@ void InterpolatorMPI<Tinterp>::_send_faces(const MPI_Comm comm, const MPInbr& nb template <typename Tinterp> -void InterpolatorMPI<Tinterp>::_recv_faces(const MPI_Comm comm, const MPInbr& nbrs, Faces& f, std::vector<MPI_Request>& f_req) const +void InterpolatorMPI<Tinterp>::_recv_faces(const MPInbr& nbrs, Faces& f, std::vector<MPI_Request>& f_req) const { const int src[6] = {nbrs(-1,0,0), nbrs(1,0,0), nbrs(0,-1,0), nbrs(0,1,0), nbrs(0,0,-1), nbrs(0,0,1)}; Real* const data[6] = {f.x0.data, f.x1.data, f.y0.data, f.y1.data, f.z0.data, f.z1.data}; @@ -494,7 +508,7 @@ void InterpolatorMPI<Tinterp>::_recv_faces(const MPI_Comm comm, const MPInbr& nb if (src[i] >= 0) { MPI_Request req; - MPI_Irecv(data[i], Nelem[i], MPIREAL, src[i], tagbase+dir*2+side, comm, &req); + MPI_Irecv(data[i], Nelem[i], MPIREAL, src[i], tagbase+dir*2+side, m_cartcomm, &req); f_req.push_back(req); } } @@ -502,7 +516,7 @@ void InterpolatorMPI<Tinterp>::_recv_faces(const MPI_Comm comm, const MPInbr& nb template <typename Tinterp> -void InterpolatorMPI<Tinterp>::_send_edges(const MPI_Comm comm, const MPInbr& nbrs, const Edges& e, std::vector<MPI_Request>& e_req) const +void InterpolatorMPI<Tinterp>::_send_edges(const MPInbr& nbrs, const Edges& e, std::vector<MPI_Request>& e_req) const { const int dstX[4] = {nbrs(0,-1,-1), nbrs(0,1,-1), nbrs(0,1,1), nbrs(0,-1,1)}; const int dstY[4] = {nbrs(-1,0,-1), nbrs(-1,0,1), nbrs(1,0,1), nbrs(1,0,-1)}; @@ -520,19 +534,19 @@ void InterpolatorMPI<Tinterp>::_send_edges(const MPI_Comm comm, const MPInbr& nb if (dstX[i] >= 0) { MPI_Request req; - MPI_Isend(dataX[i], NelemX[i], MPIREAL, dstX[i], tagbase+0*4+edge, comm, &req); + MPI_Isend(dataX[i], NelemX[i], MPIREAL, dstX[i], tagbase+0*4+edge, m_cartcomm, &req); e_req.push_back(req); } if (dstY[i] >= 0) { MPI_Request req; - MPI_Isend(dataY[i], NelemY[i], MPIREAL, dstY[i], tagbase+1*4+edge, comm, &req); + MPI_Isend(dataY[i], NelemY[i], MPIREAL, dstY[i], tagbase+1*4+edge, m_cartcomm, &req); e_req.push_back(req); } if (dstZ[i] >= 0) { MPI_Request req; - MPI_Isend(dataZ[i], NelemZ[i], MPIREAL, dstZ[i], tagbase+2*4+edge, comm, &req); + MPI_Isend(dataZ[i], NelemZ[i], MPIREAL, dstZ[i], tagbase+2*4+edge, m_cartcomm, &req); e_req.push_back(req); } } @@ -540,7 +554,7 @@ void InterpolatorMPI<Tinterp>::_send_edges(const MPI_Comm comm, const MPInbr& nb template <typename Tinterp> -void InterpolatorMPI<Tinterp>::_recv_edges(const MPI_Comm comm, const MPInbr& nbrs, Edges& e, std::vector<MPI_Request>& e_req) const +void InterpolatorMPI<Tinterp>::_recv_edges(const MPInbr& nbrs, Edges& e, std::vector<MPI_Request>& e_req) const { const int srcX[4] = {nbrs(0,-1,-1), nbrs(0,1,-1), nbrs(0,1,1), nbrs(0,-1,1)}; const int srcY[4] = {nbrs(-1,0,-1), nbrs(-1,0,1), nbrs(1,0,1), nbrs(1,0,-1)}; @@ -558,19 +572,19 @@ void InterpolatorMPI<Tinterp>::_recv_edges(const MPI_Comm comm, const MPInbr& nb if (srcX[i] >= 0) { MPI_Request req; - MPI_Irecv(dataX[i], NelemX[i], MPIREAL, srcX[i], tagbase+0*4+edge, comm, &req); + MPI_Irecv(dataX[i], NelemX[i], MPIREAL, srcX[i], tagbase+0*4+edge, m_cartcomm, &req); e_req.push_back(req); } if (srcY[i] >= 0) { MPI_Request req; - MPI_Irecv(dataY[i], NelemY[i], MPIREAL, srcY[i], tagbase+1*4+edge, comm, &req); + MPI_Irecv(dataY[i], NelemY[i], MPIREAL, srcY[i], tagbase+1*4+edge, m_cartcomm, &req); e_req.push_back(req); } if (srcZ[i] >= 0) { MPI_Request req; - MPI_Irecv(dataZ[i], NelemZ[i], MPIREAL, srcZ[i], tagbase+2*4+edge, comm, &req); + MPI_Irecv(dataZ[i], NelemZ[i], MPIREAL, srcZ[i], tagbase+2*4+edge, m_cartcomm, &req); e_req.push_back(req); } } @@ -578,7 +592,7 @@ void InterpolatorMPI<Tinterp>::_recv_edges(const MPI_Comm comm, const MPInbr& nb template <typename Tinterp> -void InterpolatorMPI<Tinterp>::_send_corners(const MPI_Comm comm, const MPInbr& nbrs, const Corners& c, std::vector<MPI_Request>& c_req) const +void InterpolatorMPI<Tinterp>::_send_corners(const MPInbr& nbrs, const Corners& c, std::vector<MPI_Request>& c_req) const { const int dst[8] = {nbrs(-1,-1,-1), nbrs(-1,1,-1), nbrs(-1,1,1), nbrs(-1,-1,1), nbrs(1,-1,-1), nbrs(1,1,-1), nbrs(1,1,1), nbrs(1,-1,1)}; Real* const data[8] = {c.x00.data, c.x01.data, c.x02.data, c.x03.data, c.x10.data, c.x11.data, c.x12.data, c.x13.data}; @@ -591,7 +605,7 @@ void InterpolatorMPI<Tinterp>::_send_corners(const MPI_Comm comm, const MPInbr& if (dst[i] >= 0) { MPI_Request req; - MPI_Isend(data[i], Nelem[i], MPIREAL, dst[i], tagbase+side*4+cycl, comm, &req); + MPI_Isend(data[i], Nelem[i], MPIREAL, dst[i], tagbase+side*4+cycl, m_cartcomm, &req); c_req.push_back(req); } } @@ -599,7 +613,7 @@ void InterpolatorMPI<Tinterp>::_send_corners(const MPI_Comm comm, const MPInbr& template <typename Tinterp> -void InterpolatorMPI<Tinterp>::_recv_corners(const MPI_Comm comm, const MPInbr& nbrs, Corners& c, std::vector<MPI_Request>& c_req) const +void InterpolatorMPI<Tinterp>::_recv_corners(const MPInbr& nbrs, Corners& c, std::vector<MPI_Request>& c_req) const { const int src[8] = {nbrs(-1,-1,-1), nbrs(-1,1,-1), nbrs(-1,1,1), nbrs(-1,-1,1), nbrs(1,-1,-1), nbrs(1,1,-1), nbrs(1,1,1), nbrs(1,-1,1)}; Real* const data[8] = {c.x00.data, c.x01.data, c.x02.data, c.x03.data, c.x10.data, c.x11.data, c.x12.data, c.x13.data}; @@ -612,7 +626,7 @@ void InterpolatorMPI<Tinterp>::_recv_corners(const MPI_Comm comm, const MPInbr& if (src[i] >= 0) { MPI_Request req; - MPI_Irecv(data[i], Nelem[i], MPIREAL, src[i], tagbase+side*4+cycl, comm, &req); + MPI_Irecv(data[i], Nelem[i], MPIREAL, src[i], tagbase+side*4+cycl, m_cartcomm, &req); c_req.push_back(req); } } diff --git a/src/IsoExtractorMPI.h b/src/IsoExtractorMPI.h @@ -180,7 +180,7 @@ private: MPI_Comm_rank(comm, &rank); MPI_Comm_size(comm, &size); - int rc = MPI_File_open(MPI_COMM_SELF, fname.c_str(), MPI_MODE_APPEND | MPI_MODE_WRONLY, MPI_INFO_NULL, &file_id); + int rc = MPI_File_open(MPI_COMM_SELF, (char*)fname.c_str(), MPI_MODE_APPEND | MPI_MODE_WRONLY, MPI_INFO_NULL, &file_id); if (rc) { if (0 == rank)