commit 3d8aa1af3097df8c83054ee82f03f11a1a9845b6
parent cf9b8eab84e6e06c21e638fc69a13e977e2a4663
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date: Thu, 1 Dec 2016 17:35:05 +0100
Merge branch 'master' of gitlab.ethz.ch:fabianw/easy_iso
Diffstat:
5 files changed, 244 insertions(+), 88 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/Interpolator.h b/src/Interpolator.h
@@ -13,6 +13,7 @@
#include <cstring>
#include <string>
#include <functional>
+#include <chrono>
#ifdef _USE_HDF_
#include <hdf5.h>
#endif /* _USE_HDF_ */
@@ -33,6 +34,8 @@ public:
m_extent(p("extent").asDouble(1.0))
{
// load the data
+ std::chrono::time_point<std::chrono::system_clock> start, end;
+ start = std::chrono::system_clock::now();
if (p("load").asString("h5") == "h5") _load_h5(p);
else if (p("load").asString("h5") == "wavelet") _load_wavelet(p);
else
@@ -40,6 +43,9 @@ public:
std::cerr << "ERROR: Interpolator: Undefined loader \"" << p("load").asString() << "\"" << std::endl;
abort();
}
+ end = std::chrono::system_clock::now();
+ std::chrono::duration<double> elapsed_seconds = end-start;
+ std::cout << "Loading data done. Elapsed time = " << elapsed_seconds.count() << "s" << std::endl;
// cell centered or nodal values?
if (p("data").asString("cell") == "cell") _nCells = [](const int N) { return N; };
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,31 +59,49 @@ 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
+ std::chrono::time_point<std::chrono::system_clock> start, end;
+ start = std::chrono::system_clock::now();
if (p("load").asString("h5") == "h5") this->_load_h5_MPI(p);
else if (p("load").asString("h5") == "wavelet") this->_load_wavelet_MPI(p);
else
@@ -91,6 +110,10 @@ public:
std::cerr << "ERROR: InterpolatorMPI: Undefined loader \"" << p("load").asString() << "\"" << std::endl;
abort();
}
+ end = std::chrono::system_clock::now();
+ std::chrono::duration<double> elapsed_seconds = end-start;
+ if (m_isroot)
+ std::cout << "Loading data done. Elapsed time = " << elapsed_seconds.count() << "s" << std::endl;
// cell centered or nodal values?
if (p("data").asString("cell") == "cell") this->_nCells = [](const int N) { return N; };
@@ -134,8 +157,16 @@ public:
for (int i = 0; i < 3; ++i)
this->m_origin_off[i] = origin[i];
- // fetch halos cells
- _fetch_halo(cartcomm);
+ // 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();
+ end = std::chrono::system_clock::now();
+ elapsed_seconds = end-start;
+ if (m_isroot)
+ std::cout << "Exchanging messages data done. Elapsed time = " << elapsed_seconds.count() << "s" << std::endl;
}
virtual ~InterpolatorMPI() {}
@@ -145,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;
@@ -163,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);
@@ -182,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);
@@ -199,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);
@@ -243,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)
@@ -277,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);
@@ -335,8 +368,7 @@ void InterpolatorMPI<Tinterp>::_load_wavelet_MPI(ArgumentParser& p)
abort();
}
- // Reader_WaveletCompressionMPI myreader(m_comm, filename, byte_swap, wavelet_type);
- Reader_WaveletCompressionMPI myreader(MPI::COMM_WORLD, filename, byte_swap, wavelet_type);
+ Reader_WaveletCompressionMPI myreader(m_cartcomm, filename, byte_swap, wavelet_type);
myreader.load_file();
const int NBX = myreader.xblocks();
const int NBY = myreader.yblocks();
@@ -386,7 +418,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
@@ -394,7 +426,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();
@@ -403,51 +435,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);
+ _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};
@@ -460,7 +486,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);
}
}
@@ -468,7 +494,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};
@@ -481,7 +507,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);
}
}
@@ -489,7 +515,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)};
@@ -507,19 +533,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);
}
}
@@ -527,7 +553,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)};
@@ -545,19 +571,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);
}
}
@@ -565,7 +591,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};
@@ -578,7 +604,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);
}
}
@@ -586,7 +612,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};
@@ -599,7 +625,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/IsoExtractor.h b/src/IsoExtractor.h
@@ -11,14 +11,15 @@
#include <string>
#include <ctime>
#include <cstdlib>
+#include <chrono>
#include "common.h"
#include "ArgumentParser.h"
#include "Interpolator.h"
struct Triangle
{
- Real vertices[3][3];
-};
+ float vertices[3][3];
+} __attribute__((packed));
static const Real gVertexOffset[8][3] =
{
@@ -528,6 +529,11 @@ void IsoExtractor<Tkernel,Tinterpolator>::_extractIso(const Real isoval,
abort();
}
+ std::chrono::time_point<std::chrono::system_clock> start, end;
+ start = std::chrono::system_clock::now();
+ if (verbose)
+ std::cout << "Extracting isovalue " << isoval << std::endl;
+
const size_t nCubes = static_cast<size_t>(Nx) * Ny * Nz;
size_t totalSize = 0;
#pragma omp parallel
@@ -578,6 +584,10 @@ void IsoExtractor<Tkernel,Tinterpolator>::_extractIso(const Real isoval,
m_triList.insert(m_triList.end(), myTriList.begin(), myTriList.end());
}
}
+ end = std::chrono::system_clock::now();
+ std::chrono::duration<double> elapsed_seconds = end-start;
+ if (verbose)
+ std::cout << "Extraction done. Elapsed time = " << elapsed_seconds.count() << "s" << std::endl;
}
diff --git a/src/IsoExtractorMPI.h b/src/IsoExtractorMPI.h
@@ -7,9 +7,16 @@
#define ISOEXTRACTORMPI_H_XABTB0WN
#include <mpi.h>
+#include <cstring>
#include "IsoExtractor.h"
#include "InterpolatorMPI.h"
+struct TriangleConnect
+{
+ unsigned char n;
+ int v[3];
+} __attribute__((packed));
+
template <typename Tkernel, template<typename> class Tinterpolator=InterpolatorMPI>
class IsoExtractorMPI : public IsoExtractor<Tkernel,Tinterpolator>
{
@@ -55,6 +62,7 @@ public:
this->m_interp.setOrigin(origin);
this->_extractIso(isoval, Nx, Ny, Nz, chx, chy, chz, Ox, Oy, Oz, this->m_interp.isroot());
+ MPI_Barrier(this->m_interp.getComm());
_writePLY_MPI(fout, this->m_triList);
}
@@ -108,7 +116,8 @@ private:
// write content
// _writePLY_binary_MPI(comm, tList);
- _writePLY_binary_MPI_lazy(tList, basename+".ply");
+ // _writePLY_binary_MPI_lazy(tList, basename+".ply");
+ _writePLY_binary_MPI(comm, tList, basename+".ply");
}
void _writePLY_binary_MPI_lazy(const std::vector<Triangle>& tList, const std::string fname)
@@ -163,28 +172,90 @@ private:
}
}
- // void _writePLY_binary_MPI(const MPI_Comm comm, const std::vector<Triangle>& tList)
- // {
- // for (auto t : tList)
- // {
- // float tc[3][3];
- // for (int j = 0; j < 3; ++j)
- // {
- // tc[j][0] = t.vertices[j][0];
- // tc[j][1] = t.vertices[j][1];
- // tc[j][2] = t.vertices[j][2];
- // }
- // ply.write((char*)&tc[0][0], 9*sizeof(float));
- // }
-
- // const char nVert = 3;
- // for (int i = 0; i < tList.size()*3; i += 3)
- // {
- // int tv[3] = {i, i+1, i+2};
- // ply.write(&nVert, sizeof(char));
- // ply.write((char*)&tv[0], 3*sizeof(int));
- // }
- // }
+ void _writePLY_binary_MPI(const MPI_Comm comm, const std::vector<Triangle>& tList, const std::string fname)
+ {
+ int rank, size;
+ MPI_Status status;
+ MPI_File file_id;
+ MPI_Comm_rank(comm, &rank);
+ MPI_Comm_size(comm, &size);
+
+ 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)
+ std::cerr << "ERROR: IsoExtractorMPI: Can not open file." << std::endl;
+ abort();
+ }
+
+ std::vector<unsigned long long> allsize(size);
+ unsigned long long mytri = tList.size();
+ MPI_Allgather(&mytri, 1, MPI_UNSIGNED_LONG_LONG, allsize.data(), 1, MPI_UNSIGNED_LONG_LONG, comm);
+
+ long long offset = 0;
+ MPI_File_get_position(file_id, &offset);
+
+ assert(sizeof(Triangle) == 9*sizeof(float));
+ for (int i = 0; i < rank; ++i)
+ offset = offset + allsize[i]*sizeof(Triangle);
+
+ Triangle* tarr = nullptr;
+ if (allsize[rank] > 0)
+ {
+ tarr = new Triangle[allsize[rank]];
+ void* const base = &tarr[0].vertices[0][0];
+ for (size_t i=0; i<allsize[rank]; ++i)
+ {
+ void* const dst = base + i*sizeof(Triangle);
+ const void* const src = &tList[i].vertices[0][0];
+ memcpy(dst, src, sizeof(Triangle));
+ }
+ // for (size_t i=0; i<allsize[rank]; ++i)
+ // {
+ // for (int j = 0; j < 3; ++j)
+ // {
+ // tarr[i].vertices[j][0] = tList[i].vertices[j][0];
+ // tarr[i].vertices[j][1] = tList[i].vertices[j][1];
+ // tarr[i].vertices[j][2] = tList[i].vertices[j][2];
+ // }
+ // }
+
+ MPI_File_write_at(file_id, offset, (char *)&tarr[0].vertices[0][0], allsize[rank]*sizeof(Triangle), MPI_CHAR, &status);
+ }
+ offset = offset + allsize[rank]*sizeof(Triangle);
+ MPI_Barrier(comm);
+ MPI_Bcast(&offset, 1, MPI_LONG_LONG, size-1, comm);
+
+ if (tarr != nullptr)
+ delete [] tarr;
+
+ for (int i = 0; i < rank; ++i)
+ offset = offset + allsize[i]*sizeof(TriangleConnect);
+
+ TriangleConnect* tcarr = nullptr;
+ if (allsize[rank] > 0)
+ {
+ tcarr = new TriangleConnect[allsize[rank]];
+ int base = 0;
+ for (int i = 0; i < rank; ++i)
+ base += 3*allsize[i];
+
+ for (size_t i=0; i<allsize[rank]; ++i)
+ {
+ tcarr[i].n = 3;
+ tcarr[i].v[0] = base + 3*i + 0;
+ tcarr[i].v[1] = base + 3*i + 1;
+ tcarr[i].v[2] = base + 3*i + 2;
+ }
+
+ MPI_File_write_at(file_id, offset, (char *)&tcarr[0].n, allsize[rank]*sizeof(TriangleConnect), MPI_CHAR, &status);
+ }
+
+ if (tcarr != nullptr)
+ delete [] tcarr;
+
+ MPI_File_close(&file_id);
+ }
};
#endif /* ISOEXTRACTORMPI_H_XABTB0WN */
