Commit 76df83cc authored by fabianw's avatar fabianw
Browse files

added MPI'ed version for iso extraction (almost complete)

parent 52ff7678
......@@ -4,13 +4,13 @@ CC = mpic++
# CC = g++
HDR = $(wildcard src/*.h)
SRC = $(wildcard src/*.cpp)
OBJ = ${SRC:.cpp=.o}
.PHONY: clean third
easyIso: $(OBJ)
$(CC) $(CPPFLAGS) $(INC) -o easyIso++ $(OBJ) $(LIB)
easyIso: src/easy_iso.o
$(CC) $(CPPFLAGS) $(INC) -o easyIso++ src/easy_iso.o $(LIB)
easyIsoMPI: src/easy_iso_MPI.o
$(CC) $(CPPFLAGS) $(INC) -o easyIsoMPI++ src/easy_iso_MPI.o $(LIB)
third:
cd third_party; \
......@@ -22,5 +22,5 @@ third:
clean:
find . -iname "*~" -exec rm -f {} \;
rm -f $(OBJ)
rm -f easyIso++
rm -f src/*.o
rm -f easyIso++ easyIsoMPI++
// File : InterpolatorMPI.h
// Date : Mon Nov 28 16:15:41 2016
// Author : Fabian Wermelinger
// Description: MPI'ed Interpolator
// Copyright 2016 ETH Zurich. All Rights Reserved.
#ifndef INTERPOLATORMPI_H_FKS9KH2J
#define INTERPOLATORMPI_H_FKS9KH2J
#include <cassert>
#include <vector>
#include <mpi.h>
#include "Interpolator.h"
#ifdef _FLOAT_PRECISION_
#define MPIREAL MPI_FLOAT
#else
#define MPIREAL MPI_DOUBLE
#endif
class MPInbr
{
public:
MPInbr(const MPI_Comm comm)
{
int myidx[3];
int dims[3];
int periods[3];
MPI_Cart_get(comm, 3, dims, periods, myidx);
for (int i=0; i<27; ++i)
{
m_nbr[i] = -1;
int off[3] = {(i%3)-1, ((i/3)%3)-1, (i/9)-1};
int coords[3];
for (int j=0; j<3; ++j)
coords[j] = myidx[j] + off[j];
if (coords[0] < 0 || coords[0] >= dims[0]) continue;
if (coords[1] < 0 || coords[1] >= dims[1]) continue;
if (coords[2] < 0 || coords[2] >= dims[2]) continue;
int& nbr = m_nbr[(off[0]+1) + (off[1]+1)*3 + (off[2]+1)*9];
MPI_Cart_rank(comm, coords, &nbr);
}
}
inline int operator()(const int ix, const int iy, const int iz) const
{
assert(ix >= -1 && ix <= 1);
assert(iy >= -1 && iy <= 1);
assert(iz >= -1 && iz <= 1);
return m_nbr[(ix+1) + (iy+1)*3 + (iz+1)*9];
}
private:
int m_nbr[27];
};
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);
const int periodic[3] = {false};
int world_size;
MPI_Comm_size(comm, &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);
m_isroot = (0 == myrank);
// cell centered or nodal values?
if (p("data").asString("cell") == "cell") this->_getH = [](const Real e, const int N) { return e/static_cast<Real>(N); };
else if (p("data").asString("cell") == "node") this->_getH = [](const Real e, const int N) { return e/static_cast<Real>(N-1); };
else
{
if (m_isroot)
std::cerr << "ERROR: InterpolatorMPI: Unknown data attribute \"" << p("data").asString() << "\"" << std::endl;
abort();
}
// load the data
if (p("load").asString("h5") == "h5") this->_load_h5_MPI(p);
else if (p("load").asString("h5") == "wavelet") this->_load_wavelet_MPI(p);
else
{
if (m_isroot)
std::cerr << "ERROR: InterpolatorMPI: Undefined loader \"" << p("load").asString() << "\"" << std::endl;
abort();
}
const int Nx = m_PESize[0]*this->m_data.Nx();
const int Ny = m_PESize[1]*this->m_data.Ny();
const int Nz = m_PESize[2]*this->m_data.Nz();
const int Nmax = std::max(Nx, std::max(Ny, Nz));
assert(Nmax > 1);
this->m_h = this->_getH(this->m_extent, Nmax);
this->m_invh = 1.0/this->m_h;
// init coordinate transform
this->m_origin_off = 0.0;
this->m_isNodal = 1;
if (p("data").asString("cell") == "cell")
{
this->m_isNodal = 0;
this->m_origin_off += 0.5*this->m_h;
}
// fetch halos cells
_fetch_halo(cartcomm);
}
virtual ~InterpolatorMPI() {}
inline int getNx() const { return m_PESize[0]*this->m_data.Nx(); }
inline int getNy() const { return m_PESize[1]*this->m_data.Ny(); }
inline int getNz() const { return m_PESize[2]*this->m_data.Nz(); }
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 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;
const int m_PESize[3];
int m_myPEIndex[3];
bool m_isroot;
// loader
void _load_h5_MPI(ArgumentParser& p);
void _load_wavelet_MPI(ArgumentParser& p);
// halo handling
void _fetch_halo(const MPI_Comm comm);
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 _pack_faces(Faces& f, const MPInbr& nbr) const;
void _unpack_faces(const Faces& f, const MPInbr& nbr);
inline void _pack_faceX(Face_t& f, const int _s) const;
inline void _pack_faceY(Face_t& f, const int _s) const;
inline void _pack_faceZ(Face_t& f, const int _s) const;
inline void _upack_faceX(const Face_t& f, const int _d);
inline void _upack_faceY(const Face_t& f, const int _d);
inline void _upack_faceZ(const Face_t& f, const int _d);
inline void _set_halo_faceX(const int _s, const int _d);
inline void _set_halo_faceY(const int _s, const int _d);
inline void _set_halo_faceZ(const int _s, const int _d);
void _send_edges();
void _send_corners();
};
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_load_h5_MPI(ArgumentParser& parser)
{
#ifdef _USE_HDF_
const std::string group = parser("hdf5_group").asString("/data");
const std::string filename = parser("file").asString("");
if (filename == "")
{
std::cerr << "ERROR: Interpolator: -file is not specified. No input file given" << std::endl;
abort();
}
herr_t status;
hid_t file_id, dataset_id, dataspace_id, file_dataspace_id, fspace_id, fapl_id, mspace_id;
int ndims, NCH;
int maxDim[4];
if (m_isroot)
{
hsize_t dims[4];
file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
dataset_id = H5Dopen(file_id, group.c_str(), H5P_DEFAULT);
file_dataspace_id = H5Dget_space(dataset_id);
ndims = H5Sget_simple_extent_dims(file_dataspace_id, dims, NULL);
status = H5Dclose(dataset_id);
status = H5Sclose(file_dataspace_id);
status = H5Fclose(file_id);
maxDim[2] = dims[0];
maxDim[1] = dims[1];
maxDim[0] = dims[2];
maxDim[3] = dims[3];
}
MPI_Bcast(maxDim, 4, MPI_INT, 0, m_comm);
NCH = maxDim[3];
assert(maxDim[0] % m_PESize[0] == 0);
assert(maxDim[1] % m_PESize[1] == 0);
assert(maxDim[2] % m_PESize[2] == 0);
const size_t NX = maxDim[0]/m_PESize[0];
const size_t NY = maxDim[1]/m_PESize[1];
const size_t NZ = maxDim[2]/m_PESize[2];
const size_t num_elem = NX*NY*NZ*NCH;
Real* const data = new Real[num_elem];
hsize_t count[4] = {NZ, NY, NX, NCH};
hsize_t dims[4] = {maxDim[2], maxDim[1], maxDim[0], NCH};
hsize_t offset[4] = {m_myPEIndex[2]*NZ, m_myPEIndex[1]*NY, m_myPEIndex[0]*NX, 0};
H5open();
fapl_id = H5Pcreate(H5P_FILE_ACCESS);
status = H5Pset_fapl_mpio(fapl_id, m_comm, MPI_INFO_NULL);
file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, fapl_id);
status = H5Pclose(fapl_id);
dataset_id = H5Dopen2(file_id, group.c_str(), H5P_DEFAULT);
fapl_id = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(fapl_id, H5FD_MPIO_COLLECTIVE);
fspace_id = H5Dget_space(dataset_id);
H5Sselect_hyperslab(fspace_id, H5S_SELECT_SET, offset, NULL, count, NULL);
mspace_id = H5Screate_simple(4, count, NULL);
status = H5Dread(dataset_id, HDF_PRECISION, mspace_id, fspace_id, fapl_id, data);
if (NCH > 1)
{
const int channel = parser("channel").asInt(0); // data channel
const bool magnitude = parser("magnitude").asBool(false); // vector magnitude (only if NCH == 3)
if (magnitude && NCH == 3)
{
int k = 0;
for (int i = 0; i < num_elem; i += NCH)
{
const Real a = data[i+0];
const Real b = data[i+1];
const Real c = data[i+2];
data[k++] = std::sqrt(a*a + b*b + c*c);
}
}
else
{
assert(channel < NCH);
int k = 0;
for (int i = 0; i < num_elem; i += NCH)
data[k++] = data[i+channel];
}
}
this->m_data = std::move(Matrix_t(NX, NY, NZ, data));
delete [] data;
#endif /* _USE_HDF_ */
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_load_wavelet_MPI(ArgumentParser& p)
{
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_fetch_halo(const MPI_Comm comm)
{
int myrank;
MPI_Comm_rank(comm, &myrank);
MPInbr mynbrs(comm);
// int myrank;
// int size;
// MPI_Comm_rank(comm, &myrank);
// MPI_Comm_size(comm, &size);
// for (int i=0; i<size; ++i)
// {
// if (i==myrank)
// {
// std::cout << "Rank " << myrank << " nbrs" << std::endl;
// for (int i=0; i<27; ++i)
// {
// int off[3] = {(i%3)-1, ((i/3)%3)-1, (i/9)-1};
// std::cout << mynbrs(off[0], off[1], off[2]) << " ";
// }
// std::cout << std::endl << std::endl;
// }
// MPI_Barrier(comm);
// }
// 1.)
const int Nx = this->m_data.Nx();
const int Ny = this->m_data.Ny();
const int Nz = this->m_data.Nz();
// recv faces
std::vector<MPI_Request> fpend_recv;
Faces frecv(Ny,Nz, Nx,Nz, Nx,Ny);
_recv_faces(comm, mynbrs, frecv, fpend_recv);
// send faces
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);
std::vector<MPI_Status> fstat_recv(fpend_recv.size());
MPI_Waitall(fpend_recv.size(), fpend_recv.data(), fstat_recv.data());
_unpack_faces(frecv, mynbrs);
// // send edges
// Edges esend;
// // send corners
// Corners csend;
}
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
{
const int me = nbrs(0,0,0);
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)};
const Real* const data[6] = {f.x0.data, f.x1.data, f.y0.data, f.y1.data, f.z0.data, f.z1.data};
const int Nelem[6] = {f.x0._Nelem, f.x1._Nelem, f.y0._Nelem, f.y1._Nelem, f.z0._Nelem, f.z1._Nelem};
for (int i = 0; i < 6; ++i)
{
if (dst[i] >= 0)
{
MPI_Request req;
MPI_Isend(data[i], Nelem[i], MPIREAL, dst[i], me, comm, &req);
f_req.push_back(req);
}
}
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_recv_faces(const MPI_Comm comm, 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};
const int Nelem[6] = {f.x0._Nelem, f.x1._Nelem, f.y0._Nelem, f.y1._Nelem, f.z0._Nelem, f.z1._Nelem};
for (int i = 0; i < 6; ++i)
{
if (src[i] >= 0)
{
MPI_Request req;
MPI_Irecv(data[i], Nelem[i], MPIREAL, src[i], src[i], comm, &req);
f_req.push_back(req);
}
}
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_pack_faces(Faces& f, const MPInbr& nbr) const
{
const int Nx = this->m_data.Nx();
const int Ny = this->m_data.Ny();
const int Nz = this->m_data.Nz();
if (nbr(-1,0,0) >= 0) _pack_faceX(f.x0, 0);
if (nbr(1,0,0) >= 0) _pack_faceX(f.x1, Nx-Face_t::_P);
if (nbr(0,-1,0) >= 0) _pack_faceY(f.y0, 0);
if (nbr(0,1,0) >= 0) _pack_faceY(f.y1, Ny-Face_t::_P);
if (nbr(0,0,-1) >= 0) _pack_faceZ(f.z0, 0);
if (nbr(0,0,1) >= 0) _pack_faceZ(f.z1, Nz-Face_t::_P);
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_unpack_faces(const Faces& f, const MPInbr& nbr)
{
const int Nx = this->m_data.Nx();
const int Ny = this->m_data.Ny();
const int Nz = this->m_data.Nz();
if (nbr(-1,0,0) >= 0) _upack_faceX(f.x0, -Face_t::_P);
else _set_halo_faceX(0, -Face_t::_P);
if (nbr(1,0,0) >= 0) _upack_faceX(f.x1, Nx);
else _set_halo_faceX(Nx-1, Nx);
if (nbr(0,-1,0) >= 0) _upack_faceY(f.y0, -Face_t::_P);
else _set_halo_faceY(0, -Face_t::_P);
if (nbr(0,1,0) >= 0) _upack_faceY(f.y1, Ny);
else _set_halo_faceY(Ny-1, Ny);
if (nbr(0,0,-1) >= 0) _upack_faceZ(f.z0, -Face_t::_P);
else _set_halo_faceZ(0, -Face_t::_P);
if (nbr(0,0,1) >= 0) _upack_faceZ(f.z1, Nz);
else _set_halo_faceZ(Nz-1, Nz);
}
// PUP kernels
#include "InterpolatorPUPMPI.h"
#endif /* INTERPOLATORMPI_H_FKS9KH2J */
// File : InterpolatorPUPMPI.h
// Date : Tue Nov 29 21:00:01 2016
// Author : Fabian Wermelinger
// Description: pack/unpack kernels
// Copyright 2016 ETH Zurich. All Rights Reserved.
#ifndef INTERPOLATORPUPMPI_H_EPMRY0VT
#define INTERPOLATORPUPMPI_H_EPMRY0VT
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_pack_faceX(Face_t& f, const int _s) const
{
for (int iz=0; iz<this->m_data.Nz(); ++iz)
for (int iy=0; iy<this->m_data.Ny(); ++iy)
for (int ix=0; ix<Face_t::_P; ++ix)
f(iy,iz,ix) = this->m_data(ix+_s,iy,iz);
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_pack_faceY(Face_t& f, const int _s) const
{
for (int iz=0; iz<this->m_data.Nz(); ++iz)
for (int iy=0; iy<Face_t::_P; ++iy)
for (int ix=0; ix<this->m_data.Nx(); ++ix)
f(ix,iz,iy) = this->m_data(ix,iy+_s,iz);
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_pack_faceZ(Face_t& f, const int _s) const
{
for (int iz=0; iz<Face_t::_P; ++iz)
for (int iy=0; iy<this->m_data.Ny(); ++iy)
for (int ix=0; ix<this->m_data.Nx(); ++ix)
f(ix,iy,iz) = this->m_data(ix,iy,iz+_s);
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_upack_faceX(const Face_t& f, const int _d)
{
for (int iz=0; iz<this->m_data.Nz(); ++iz)
for (int iy=0; iy<this->m_data.Ny(); ++iy)
for (int ix=0; ix<Face_t::_P; ++ix)
this->m_data(ix+_d,iy,iz) = f(iy,iz,ix);
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_upack_faceY(const Face_t& f, const int _d)
{
for (int iz=0; iz<this->m_data.Nz(); ++iz)
for (int iy=0; iy<Face_t::_P; ++iy)
for (int ix=0; ix<this->m_data.Nx(); ++ix)
this->m_data(ix,iy+_d,iz) = f(ix,iz,iy);
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_upack_faceZ(const Face_t& f, const int _d)
{
for (int iz=0; iz<Face_t::_P; ++iz)
for (int iy=0; iy<this->m_data.Ny(); ++iy)
for (int ix=0; ix<this->m_data.Nx(); ++ix)
this->m_data(ix,iy,iz+_d) = f(ix,iy,iz);
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_set_halo_faceX(const int _s, const int _d)
{
for (int iz=0; iz<this->m_data.Nz(); ++iz)
for (int iy=0; iy<this->m_data.Ny(); ++iy)
for (int ix=0; ix<Face_t::_P; ++ix)
this->m_data(ix+_d,iy,iz) = this->m_data(_s,iy,iz);
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_set_halo_faceY(const int _s, const int _d)
{
for (int iz=0; iz<this->m_data.Nz(); ++iz)
for (int iy=0; iy<Face_t::_P; ++iy)
for (int ix=0; ix<this->m_data.Nx(); ++ix)
this->m_data(ix,iy+_d,iz) = this->m_data(ix,_s,iz);
}
template <typename Tinterp>
void InterpolatorMPI<Tinterp>::_set_halo_faceZ(const int _s, const int _d)
{
for (int iz=0; iz<Face_t::_P; ++iz)
for (int iy=0; iy<this->m_data.Ny(); ++iy)
for (int ix=0; ix<this->m_data.Nx(); ++ix)
this->m_data(ix,iy,iz+_d) = this->m_data(ix,iy,_s);
}
#endif /* INTERPOLATORPUPMPI_H_EPMRY0VT */
// File : IsoExtractorMPI.h
// Date : Tue Nov 29 10:18:26 2016
// Author : Fabian Wermelinger
// Description: MPI'ed isosurface extractor
// Copyright 2016 ETH Zurich. All Rights Reserved.
#ifndef ISOEXTRACTORMPI_H_XABTB0WN
#define ISOEXTRACTORMPI_H_XABTB0WN
#include <mpi.h>
#include "IsoExtractor.h"
#include "InterpolatorMPI.h"
template <typename Tkernel, template<typename> class Tinterpolator=InterpolatorMPI>
class IsoExtractorMPI : public IsoExtractor<Tkernel,Tinterpolator>
{
public:
IsoExtractorMPI(ArgumentParser& parser) : IsoExtractor<Tkernel,Tinterpolator>(parser) {}
~IsoExtractorMPI() {}
void extract(const Real isoval, const std::string fout)
{
const Real cscale = this->m_parser("cube_scale").asDouble(1.0);
assert(cscale > 0.0);
const int Nx = static_cast<Real>(this->m_interp.getLocalNx()) / cscale;
const int Ny = static_cast<Real>(this->m_interp.getLocalNy()) / cscale;
const int Nz = static_cast<Real>(this->m_interp.getLocalNz()) / cscale;
assert(Nx > 1);
assert(Ny > 1);
assert(Nz > 1);
int pesize[3];
this->m_interp.getPESize(pesize);
// TODO: (fabianw@mavt.ethz.ch; Tue Nov 29 23:51:01 2016) need better
// grid def here
const Real ch = this->m_interp.getExtent() / (std::max(pesize[0]*Nx, std::max(pesize[1]*Ny, pesize[2]*Nz)) - this->m_interp.isNodal());
int peidx[3];
this->m_interp.getPEIndex(peidx);
const Real Ox = ch*Nx*peidx[0];
const Real Oy = ch*Ny*peidx[1];
const Real Oz = ch*Nz*peidx[2];
this->_extractIso(isoval, Nx, Ny, Nz, ch, Ox, Oy, Oz, this->m_interp.isroot());
_writePLY_MPI(fout, this->m_triList);
}
private:
void _writePLY_MPI(const std::string basename, const std::vector<Triangle>& tList)
{
const MPI_Comm comm = this->m_interp.getComm();
const unsigned int mytriangles = tList.size();
unsigned long long alltriangles;
MPI_Reduce(&mytriangles, &alltriangles, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, comm);
const bool isroot = this->m_interp.isroot();
time_t rawtime;
std::time(&rawtime);
struct tm* timeinfo = std::localtime(&rawtime);
char buf[256];
std::strftime(buf, 256, "%A, %h %d %Y, %r %Z %z", timeinfo);
const unsigned int one = 1;