diff --git a/expui/Centering.cc b/expui/Centering.cc index 8b00b1f07..0afd3f5fa 100644 --- a/expui/Centering.cc +++ b/expui/Centering.cc @@ -143,7 +143,7 @@ namespace Utility int count = 0; for (auto it=combo.rbegin(); it!=combo.rend(); it++) { - if (count > Nsort) break; + if (++count > Nsort) break; for (int k=0; k<3; k++) ctr[k] += it->first * points[it->second].get(k); dentot += it->first; diff --git a/exputil/ParticleReader.cc b/exputil/ParticleReader.cc index ea63c0400..76c68066c 100644 --- a/exputil/ParticleReader.cc +++ b/exputil/ParticleReader.cc @@ -35,7 +35,7 @@ namespace PR { { {"Gas", 0}, {"Halo", 1}, {"Disk", 2}, {"Bulge", 3}, {"Stars", 4}, {"Bndry", 5}}; - GadgetNative::GadgetNative(const std::vector& files, bool verbose) + GadgetNative::GadgetNative(const std::vector& files, bool verbose) : ParticleReader() { _files = files; // Copy file list (bunch) _verbose = verbose; @@ -317,7 +317,7 @@ namespace PR { { {"Gas", 0}, {"Halo", 1}, {"Disk", 2}, {"Bulge", 3}, {"Stars", 4}, {"Bndry", 5}}; - GadgetHDF5::GadgetHDF5(const std::vector& files, bool verbose) + GadgetHDF5::GadgetHDF5(const std::vector& files, bool verbose) : ParticleReader() { _files = files; _verbose = verbose; @@ -705,7 +705,7 @@ namespace PR { } - PSPhdf5::PSPhdf5(const std::vector& files, bool verbose) + PSPhdf5::PSPhdf5(const std::vector& files, bool verbose) : ParticleReader() { _files = files; _verbose = verbose; @@ -1896,7 +1896,7 @@ namespace PR { std::vector ParticleReader::readerTypes {"PSPout", "PSPspl", "GadgetNative", "GadgetHDF5", "PSPhdf5", - "TipsyNative", "TipsyXDR", "Bonsai"}; + "TipsyNative", "TipsyXDR", "Bonsai1", "Bonsai"}; std::vector> @@ -2031,7 +2031,8 @@ namespace PR { exit(1); } #endif - + else if (reader.find("Bonsai1") == 0) + ret = std::make_shared(file, Tipsy::TipsyType::bonsai1, verbose); else if (reader.find("Bonsai") == 0) ret = std::make_shared(file, Tipsy::TipsyType::bonsai, verbose); else { @@ -2071,6 +2072,8 @@ namespace PR { // Native tipsy with ID conversion else if (ttype == TipsyType::bonsai) ps = std::make_shared(file); + else if (ttype == TipsyType::bonsai1) + ps = std::make_shared(file); // Make a tipsy xdr reader else { #ifdef HAVE_XDR @@ -2110,6 +2113,8 @@ namespace PR { ps = std::make_shared(*curfile); else if (ttype == TipsyType::bonsai) ps = std::make_shared(*curfile); + else if (ttype == TipsyType::bonsai1) + ps = std::make_shared(*curfile); else { #ifdef HAVE_XDR ps = std::make_shared(*curfile); @@ -2125,7 +2130,7 @@ namespace PR { } Tipsy::Tipsy(const std::string& file, TipsyType Type, - bool verbose) + bool verbose) : ParticleReader() { ttype = Type; files.push_back(file); @@ -2137,7 +2142,7 @@ namespace PR { } Tipsy::Tipsy(const std::vector& filelist, TipsyType Type, - bool verbose) + bool verbose) : ParticleReader() { files = filelist; @@ -2184,7 +2189,8 @@ namespace PR { P.pos[k] = ps->gas_particles[pcount].pos[k]; P.vel[k] = ps->gas_particles[pcount].vel[k]; } - if (ttype == TipsyType::bonsai) P.indx = ps->gas_particles[pcount].ID(); + P.indx = ps->getIndexOffset(TipsyReader::Ptype::gas) + pcount + 1; + pcount++; return; } @@ -2196,7 +2202,13 @@ namespace PR { P.pos[k] = ps->dark_particles[pcount].pos[k]; P.vel[k] = ps->dark_particles[pcount].vel[k]; } - if (ttype == TipsyType::bonsai) P.indx = ps->dark_particles[pcount].ID(); + if (ttype == TipsyType::bonsai) + P.indx = ps->dark_particles[pcount].ID2(); + else if (ttype == TipsyType::bonsai1) + P.indx = ps->dark_particles[pcount].ID(); + else + P.indx = ps->getIndexOffset(TipsyReader::Ptype::dark) + pcount + 1; + pcount++; return; } @@ -2208,7 +2220,13 @@ namespace PR { P.pos[k] = ps->star_particles[pcount].pos[k]; P.vel[k] = ps->star_particles[pcount].vel[k]; } - if (ttype == TipsyType::bonsai) P.indx = ps->star_particles[pcount].ID(); + if (ttype == TipsyType::bonsai) + P.indx = ps->star_particles[pcount].ID2(); + else if (ttype == TipsyType::bonsai1) + P.indx = ps->star_particles[pcount].ID(); + else + P.indx = ps->getIndexOffset(TipsyReader::Ptype::star) + pcount + 1; + pcount++; return; } diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 7dc2c0e73..5bf18751e 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -2910,6 +2910,7 @@ void SLGridSlab::compute_table(struct TableSlab* table, int KX, int KY) if (VERBOSE && iflag[i] != 0) { + if (iflag[i] > -10) { cout << std::setw(14) << "x" << std::setw(25) << "u(x)" diff --git a/include/ParticleReader.H b/include/ParticleReader.H index 0d150dc45..8dc9a851e 100644 --- a/include/ParticleReader.H +++ b/include/ParticleReader.H @@ -456,7 +456,7 @@ namespace PR public: //! Default constructor - PSP(bool verbose) : VERBOSE(verbose) { init(); } + PSP(bool verbose) : VERBOSE(verbose), ParticleReader() { init(); } //! Destructor virtual ~PSP() { if (in.is_open()) in.close(); } @@ -573,7 +573,7 @@ namespace PR public: - enum class TipsyType { native, xdr, bonsai }; + enum class TipsyType { native, xdr, bonsai1, bonsai }; private: //! List of files diff --git a/include/tipsy.H b/include/tipsy.H index 1f38c03c8..fefe2bcbc 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -1,6 +1,13 @@ #ifndef _tipsy_H #define _tipsy_H +#include + +// For MPI awareness +// +#include + + // Forward declare PR::Tipsy so it can friend TipsyFile // namespace PR { @@ -33,12 +40,6 @@ namespace TipsyReader Real metals ; Real phi ; //@} - - //! Convert phi to index - int ID() const { - union id {Real v; int i;} u; - u.v = phi; return u.i; - } } ; struct dark_particle @@ -52,11 +53,18 @@ namespace TipsyReader Real phi ; //@} - //! Convert phi to index + //! Convert phi to 32-bit index for Bonsai v1 int ID() const { union id {Real v; int i;} u; u.v = phi; return u.i; } + + //! Convert phi to 64-bit index for Bonsai v2 + unsigned long ID2() const { + union id {Real v[2]; uint64_t i;} u; + u.v[0] = eps; u.v[1] = phi; return u.i; + } + } ; struct star_particle @@ -72,11 +80,17 @@ namespace TipsyReader Real phi ; //@} - //! Convert phi to index + //! Convert phi to 32-bit index for Bonsai v1 int ID() const { union id {Real v; int i;} u; u.v = phi; return u.i; } + + //! Convert exp, phi to 64bit index for Bonsai v2 + unsigned long ID2() const { + union id {Real v[2]; uint64_t i;} u; + u.v[0] = eps; u.v[1] = phi; return u.i; + } } ; struct Header @@ -92,6 +106,10 @@ namespace TipsyReader #endif } ; + //! Used to identify particle type for indexing + enum class Ptype { gas, dark, star }; + + //! Base class for Tipsy file reading class TipsyFile { protected: @@ -100,6 +118,9 @@ namespace TipsyReader virtual void read_dark() = 0; virtual void read_star() = 0; + //! Offset for indexing particles in the file + std::map index_offset_map + { {Ptype::gas, 0}, {Ptype::dark, 0}, {Ptype::star, 0} }; public: @@ -111,6 +132,12 @@ namespace TipsyReader virtual int readParticles() = 0; virtual ~TipsyFile() {} + + //! Get the index offset for a given particle type + unsigned long getIndexOffset(const Ptype ptype) + { + return index_offset_map[ptype]; + } }; #ifdef HAVE_XDR @@ -118,13 +145,20 @@ namespace TipsyReader #include #include + //! Tipsy XDR file reader class TipsyXDR : public TipsyFile { private: + //! The input file stream FILE* input; + + //! The XDR stream XDR xdrs; + //! Cache the XDR stream position at beginning of particle group + unsigned long curpos = 0; + int xdr_header() { if (xdr_double(&xdrs, &header.time) != TRUE) return 0; @@ -157,26 +191,95 @@ namespace TipsyReader return header.nbodies; } + //! Generate stream position for parallel read + //! \param nsize The number of particles to read + //! \param ssize The size of each particle structure in bytes + //! \note Using unsigned long to avoid wrapping issues + //! \return The number of particles to read for this process + unsigned long xdr_psize(unsigned long nsize, unsigned long ssize, Ptype ptype) + { + // Number of particles per process + unsigned long psize = nsize/numprocs; + + // Get current stream position + curpos = xdr_getpos(&xdrs); + auto pos = curpos; + + // Set the index offset for indexing particles + index_offset_map[ptype] = psize * myid; + + // Move the stream position to the start of this process' data + pos += index_offset_map[ptype] * ssize; + if (xdr_setpos(&xdrs, pos) == 0) { + throw std::runtime_error("TipsyFile: xdr_setpos failed"); + } + + // If this is the last process, read the remainder + if (myid == numprocs-1) { + index_offset_map[ptype] = (numprocs-1) * psize; + psize = nsize - index_offset_map[ptype]; + }} + + return psize; + } + + //! Set stream position to next group of particles + //! \param nsize The number of particles to read + //! \param ssize The size of each particle structure in bytes + //! \note Using unsigned long to avoid wrapping issues + void xdr_pnext(unsigned long nsize, unsigned long ssize) + { + // Move the stream position to the start of the next particle + // group + if (xdr_setpos(&xdrs, curpos + nsize * ssize) == 0) { + throw std::runtime_error("TipsyFile: xdr_setpos failed"); + } + } + int xdr_read() { int N=0; if (header.nsph != 0) { - gas_particles.resize(header.nsph); + // Get the number of gas particles for this process and set + // stream position + auto psize = xdr_psize(header.nsph, sizeof(gas_particle), Ptype::gas); + + // Do the read + gas_particles.resize(psize); read_gas(); N++; + + // Set stream position to next group of particles + xdr_pnext(header.nsph, sizeof(gas_particle)); } if (header.ndark != 0) { - dark_particles.resize(header.ndark); + // Get the number of dark particles for this process and set + // stream position + auto psize = xdr_psize(header.ndark, sizeof(dark_particle), Ptype::dark); + + // Do the read + dark_particles.resize(psize); read_dark(); N++; + + // Set stream position to next group of particles + xdr_pnext(header.ndark, sizeof(dark_particle)); } if (header.nstar != 0) { - star_particles.resize(header.nstar); + // Get the number of star particles for this process and set + // stream position + auto psize = xdr_psize(header.nstar, sizeof(star_particle), Ptype::star); + + // Do the read + star_particles.resize(psize); read_star(); N++; + + // Set stream position to next group of particles + xdr_pnext(header.nstar, sizeof(star_particle)); } return N; @@ -186,7 +289,7 @@ namespace TipsyReader { if (sizeof(Real) == sizeof(float)) { xdr_vector(&xdrs, (char *) &gas_particles[0], - header.nsph*(sizeof(gas_particle)/sizeof(Real)), + gas_particles.size()*(sizeof(gas_particle)/sizeof(Real)), sizeof(Real), (xdrproc_t) xdr_float); } } @@ -195,7 +298,7 @@ namespace TipsyReader { if (sizeof(Real) == sizeof(float)) { xdr_vector(&xdrs, (char *) &dark_particles[0], - header.ndark*(sizeof(dark_particle)/sizeof(Real)), + dark_particles.size()*(sizeof(dark_particle)/sizeof(Real)), sizeof(Real), (xdrproc_t) xdr_float); } } @@ -204,7 +307,7 @@ namespace TipsyReader { if (sizeof(Real) == sizeof(float)) { xdr_vector(&xdrs, (char *) &star_particles[0], - header.nstar*(sizeof(star_particle)/sizeof(Real)), + star_particles.size()*(sizeof(star_particle)/sizeof(Real)), sizeof(Real), (xdrproc_t) xdr_float); } } @@ -213,6 +316,7 @@ namespace TipsyReader public: + //! Constructor that opens the file and reads the header TipsyXDR(const std::string filename) { // Attempt to read tipsy file @@ -222,8 +326,10 @@ namespace TipsyReader } } + //! Read the particles from the file int readParticles() { return xdr_read(); } + //! Destructor that removes the stream and closes the file ~TipsyXDR() { xdr_destroy(&xdrs); @@ -234,6 +340,7 @@ namespace TipsyReader #endif + //! Tipsy native file reader class TipsyNative : public TipsyFile { private: @@ -241,27 +348,34 @@ namespace TipsyReader //! The input file stream std::ifstream input; + //! Cache stream position at beginning of partcile group + std::istream::pos_type curpos = 0; + //! Read the header int native_header() { try { input.read((char *)&header, sizeof(header)); } catch (std::exception& e) { - std::cerr << "TipsyFile native_header error: " << e.what() << std::endl; - throw; + std::ostringstream s; + s << "TipsyFile native error reading header: " << e.what(); + throw std::runtime_error(s.str()); } return 1; } + //! Open the file and read the header + //! \param filename The name of the file to read int native_init(const std::string& filename) { try { input.open(filename); input.exceptions(ifstream::eofbit | ifstream::failbit | ifstream::badbit); } catch (std::exception& e) { - std::cerr << "TipsyFile native error opening file <" << filename << ">, " - << e.what() << std::endl; - throw; + std::ostringstream s; + s << "TipsyFile native error opening file <" << filename << ">: " + << e.what(); + throw std::runtime_error(s.str()); } if (native_header() != 1) { @@ -271,67 +385,147 @@ namespace TipsyReader return header.nbodies; } + //! Generate stream position for parallel read + //! \param nsize The number of particles to read + //! \param ssize The size of each particle structure in bytes + //! \note Using unsigned long to avoid wrapping issues + //! \return The number of particles to read for this process + unsigned long ios_psize(unsigned long nsize, unsigned long ssize, Ptype ptype) + { + // Number of particles per process + unsigned long psize = nsize/numprocs; + + // Get current stream position + curpos = input.tellg(); + std::streampos pos(curpos); + + // Set the index offset for indexing particles + index_offset_map[ptype] = psize * myid; + + // Move the stream position to the start of this process' data + pos += index_offset_map[ptype] * ssize; + + if (input.seekg(pos).fail()) { + throw std::runtime_error("TipsyFile: seekg failed"); + } + + std::streampos test = input.tellg(); + + // If this is the last process, read the remainder + if (myid == numprocs-1) { + index_offset_map[ptype] = (numprocs-1) * psize; + psize = nsize - index_offset_map[ptype]; + } + + return psize; + } + + //! Move to next particle group + //! \param nsize The number of particles to read + //! \param ssize The size of each particle structure in bytes + //! \note Using unsigned long to avoid wrapping issues + void ios_pnext(unsigned long nsize, unsigned long ssize) + { + // Move the stream position to the start of the next particle + // group + std::streampos pos(curpos); // Make a copy of the current position + pos += nsize * ssize; // Augment the position + + if (input.seekg(pos).fail()) { + throw std::runtime_error("TipsyFile: seekg failed"); + } + } + + //! Read the particles from the file + //! \return The number of particle groups read int native_read() { int N=0; if (header.nsph != 0) { - gas_particles.resize(header.nsph); + // Get the number of gas particles for this process and set + // stream position + auto psize = ios_psize(header.nsph, sizeof(gas_particle), Ptype::gas); + + // Do the read + gas_particles.resize(psize); read_gas(); N++; + + // Set stream position to next group of particles + ios_pnext(header.nsph, sizeof(gas_particle)); } if (header.ndark != 0) { - dark_particles.resize(header.ndark); + // Get the number of dark particles for this process and set + // stream position + auto psize = ios_psize(header.ndark, sizeof(dark_particle), Ptype::dark); + + // Do the read + dark_particles.resize(psize); read_dark(); N++; + + // Set stream position to next group of particles + ios_pnext(header.ndark, sizeof(dark_particle)); } if (header.nstar != 0) { - star_particles.resize(header.nstar); + // Get the number of star particles for this process and set + // stream position + auto psize = ios_psize(header.nstar, sizeof(star_particle), Ptype::star); + + // Do the read + star_particles.resize(psize); read_star(); N++; + + // Set stream position to next group of particles + ios_pnext(header.nstar, sizeof(star_particle)); } return N; } + //! Read the gas particles from the file void read_gas() { try { input.read((char *) &gas_particles[0], - header.nsph*sizeof(gas_particle)); + gas_particles.size()*sizeof(gas_particle)); } catch (std::exception& e) { - std::cerr << "TipsyFile native error reading sph particles: " - << e.what() << std::endl; - throw; + std::ostringstream s; + s << "TipsyFile native error reading gas particles: " << e.what(); + throw std::runtime_error(s.str()); } } + //! Read the dark particles from the file void read_dark() { try { input.read((char *) &dark_particles[0], - header.ndark*sizeof(dark_particle)); + dark_particles.size()*sizeof(dark_particle)); } catch (std::exception& e) { - std::cerr << "TipsyFile native error reading dark particles: " - << e.what() << std::endl; - throw; + std::ostringstream s; + s << "TipsyFile native error reading dark particles: " << e.what(); + throw std::runtime_error(s.str()); } } + //! Read the star particles from the file void read_star() { try { input.read((char *) &star_particles[0], - header.nstar*sizeof(star_particle)); + star_particles.size()*sizeof(star_particle)); } catch (std::exception& e) { - std::cerr << "TipsyFile native error reading star particles: " - << e.what() << std::endl; - throw; + std::ostringstream s; + s << "TipsyFile native error reading star particles: " << e.what(); + throw std::runtime_error(s.str()); } } @@ -339,6 +533,7 @@ namespace TipsyReader public: + //! Constructor that opens the file and reads the header TipsyNative(const std::string filename) { // Attempt to read tipsy file @@ -348,8 +543,10 @@ namespace TipsyReader } } + //! Read the particles from the file int readParticles() { return native_read(); } + //! Destructor that closes the file by scoping ~TipsyNative() { // Nothing diff --git a/pyEXP/UtilWrappers.cc b/pyEXP/UtilWrappers.cc index 04a80ae01..541bf93fb 100644 --- a/pyEXP/UtilWrappers.cc +++ b/pyEXP/UtilWrappers.cc @@ -173,4 +173,74 @@ void UtilityClasses(py::module &m) { Return the version. This is the same version information reported by the EXP N-body code.)"); + + m.def("setMPI", + [](bool verbose) { + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myid); + MPI_Get_processor_name(processor_name, &proc_namelen); + + workers = numprocs - 1; + + if (workers) { + MPI_Comm_group(MPI_COMM_WORLD, &world_group); + std::vector nworkers (workers); + + for (int n=1; n