From 86fde7a67ca552483eb6709925dba073d8ea0e46 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 7 Aug 2025 11:23:53 -0400 Subject: [PATCH 01/15] Enable parallel reading from a Tipsy file --- include/tipsy.H | 130 ++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 104 insertions(+), 26 deletions(-) diff --git a/include/tipsy.H b/include/tipsy.H index 1f38c03c8..6c46baf5c 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -1,6 +1,11 @@ #ifndef _tipsy_H #define _tipsy_H +// For MPI awareness +// +#include + + // Forward declare PR::Tipsy so it can friend TipsyFile // namespace PR { @@ -157,24 +162,59 @@ 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 + unsigned xdr_psize(int nsize, int ssize) + { + // Number of particles per process + int psize = nsize/numprocs; + + // Get current stream position + auto pos = xdr_getpos(&xdrs); + + // Move the stream position to the start of this process' data + pos += psize * myid * 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) psize = nsize - (numprocs-1)*psize; + + return psize; + } + int xdr_read() { int N=0; if (header.nsph != 0) { - gas_particles.resize(header.nsph); + // Get the number of gas particles for this process + auto psize = xdr_psize(header.nsph, sizeof(gas_particle)); + + // Do the read + gas_particles.resize(psize); read_gas(); N++; } if (header.ndark != 0) { - dark_particles.resize(header.ndark); + // Get the number of dark particles for this process + auto psize = xdr_psize(header.ndark, sizeof(dark_particle)); + + // Do the read + dark_particles.resize(psize); read_dark(); N++; } if (header.nstar != 0) { - star_particles.resize(header.nstar); + // Get the number of star particles for this process + auto psize = xdr_psize(header.nstar, sizeof(star_particle)); + + // Do the read + star_particles.resize(psize); read_star(); N++; } @@ -186,7 +226,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 +235,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 +244,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); } } @@ -247,8 +287,9 @@ namespace TipsyReader 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; } @@ -259,9 +300,10 @@ namespace TipsyReader 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,24 +313,60 @@ 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 + unsigned ios_psize(int nsize, int ssize) + { + // Number of particles per process + int psize = nsize/numprocs; + + // Get current stream position + auto pos = input.tellg(); + + // Move the stream position to the start of this process' data + pos += psize * myid * ssize; + input.seekg(pos); + if (input.fail()) { + throw std::runtime_error("TipsyFile: seekg failed"); + } + + // If this is the last process, read the remainder + if (myid == numprocs-1) psize = nsize - (numprocs-1)*psize; + + return psize; + } + int native_read() { int N=0; if (header.nsph != 0) { - gas_particles.resize(header.nsph); + // Get the number of gas particles for this process + auto psize = ios_psize(header.nsph, sizeof(gas_particle)); + + // Do the read + gas_particles.resize(psize); read_gas(); N++; } if (header.ndark != 0) { - dark_particles.resize(header.ndark); + // Get the number of dark particles for this process + auto psize = ios_psize(header.ndark, sizeof(dark_particle)); + + // Do the read + dark_particles.resize(psize); read_dark(); N++; } if (header.nstar != 0) { - star_particles.resize(header.nstar); + // Get the number of star particles for this process + int psize = ios_psize(header.nstar, sizeof(star_particle)); + + // Do the read + star_particles.resize(psize); read_star(); N++; } @@ -300,12 +378,12 @@ namespace TipsyReader { 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()); } } @@ -313,12 +391,12 @@ namespace TipsyReader { 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()); } } @@ -326,12 +404,12 @@ namespace TipsyReader { 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()); } } From 38bae87482421501bcf00827a17ae2f9460b56fe Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Thu, 7 Aug 2025 11:35:51 -0400 Subject: [PATCH 02/15] Update exputil/SLGridMP2.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exputil/SLGridMP2.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index bb358fd0b..a50bc5a74 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -120,7 +120,7 @@ double sphdens(double r) return 4.0*M_PI * model->get_density(r); } -// Use entire grid for for l=Lswitch: rmin=rmap/rfac, rmax=rmap*rfac with rfac=pow(10, From b4379ba217d0ef76f59d2f54da478c8b3863292f Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Thu, 7 Aug 2025 11:35:59 -0400 Subject: [PATCH 03/15] Update exputil/SLGridMP2.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exputil/SLGridMP2.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index a50bc5a74..983c5507f 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -2907,7 +2907,7 @@ void SLGridSlab::compute_table(struct TableSlab* table, int KX, int KY) << std::setw( 5) << iflag[i] << std::endl; - if (VERBOSE and iflag[i] != 0) { + if (VERBOSE && iflag[i] != 0) { if (iflag[i] > -10) { cout << std::setw(14) << "x" From 56c71caf18b62b76df2ed31cd6c86aed2bb044b5 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 7 Aug 2025 13:44:07 -0400 Subject: [PATCH 04/15] Need to reposition file pointer after reading each particle type --- include/tipsy.H | 90 +++++++++++++++++++++++++++++++++++-------- pyEXP/UtilWrappers.cc | 65 +++++++++++++++++++++++++++++++ 2 files changed, 139 insertions(+), 16 deletions(-) diff --git a/include/tipsy.H b/include/tipsy.H index 6c46baf5c..32a6308c9 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -39,7 +39,7 @@ namespace TipsyReader Real phi ; //@} - //! Convert phi to index + //! Convert phi to index for Bonsai int ID() const { union id {Real v; int i;} u; u.v = phi; return u.i; @@ -57,7 +57,7 @@ namespace TipsyReader Real phi ; //@} - //! Convert phi to index + //! Convert phi to index for Bonsai int ID() const { union id {Real v; int i;} u; u.v = phi; return u.i; @@ -77,7 +77,7 @@ namespace TipsyReader Real phi ; //@} - //! Convert phi to index + //! Convert phi to index for Bonsai int ID() const { union id {Real v; int i;} u; u.v = phi; return u.i; @@ -127,9 +127,15 @@ namespace TipsyReader { private: + //! The input file stream FILE* input; + + //! The XDR stream XDR xdrs; + //! Cache the XDR stream position at beginning of particle group + unsigned curpos = 0; + int xdr_header() { if (xdr_double(&xdrs, &header.time) != TRUE) return 0; @@ -165,13 +171,14 @@ namespace TipsyReader //! Generate stream position for parallel read //! \param nsize The number of particles to read //! \param ssize The size of each particle structure in bytes - unsigned xdr_psize(int nsize, int ssize) + unsigned xdr_psize(unsigned int nsize, unsigned int ssize) { // Number of particles per process - int psize = nsize/numprocs; + unsigned int psize = nsize/numprocs; // Get current stream position - auto pos = xdr_getpos(&xdrs); + curpos = xdr_getpos(&xdrs); + auto pos = curpos // Move the stream position to the start of this process' data pos += psize * myid * ssize; @@ -185,38 +192,60 @@ namespace TipsyReader return psize; } + //! Set stream position to next group of particles + void xdr_pnext(unsigned int nsize, unsigned int 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) { - // Get the number of gas particles for this process + // Get the number of gas particles for this process and set + // stream position auto psize = xdr_psize(header.nsph, sizeof(gas_particle)); // 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) { - // Get the number of dark particles for this process + // Get the number of dark particles for this process and set + // stream position auto psize = xdr_psize(header.ndark, sizeof(dark_particle)); // 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) { - // Get the number of star particles for this process + // Get the number of star particles for this process and set + // stream position auto psize = xdr_psize(header.nstar, sizeof(star_particle)); // 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; @@ -281,6 +310,9 @@ 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() { @@ -316,13 +348,14 @@ namespace TipsyReader //! Generate stream position for parallel read //! \param nsize The number of particles to read //! \param ssize The size of each particle structure in bytes - unsigned ios_psize(int nsize, int ssize) + unsigned ios_psize(unsigned int nsize, unsigned int ssize) { // Number of particles per process - int psize = nsize/numprocs; + unsigned int psize = nsize/numprocs; // Get current stream position - auto pos = input.tellg(); + curpos = input.tellg(); + auto pos = curpos; // Move the stream position to the start of this process' data pos += psize * myid * ssize; @@ -337,38 +370,63 @@ namespace TipsyReader return psize; } + //! Move to next particle group + void ios_pnext(unsigned int nsize, unsigned int ssize) + { + // Move the stream position to the start of the next particle + // group + auto pos = curpos; // Make a copy of the current position + pos += nsize * ssize; // Augment the position + input.seekg(pos); + if (input.fail()) { + throw std::runtime_error("TipsyFile: seekg failed"); + } + } + int native_read() { int N=0; if (header.nsph != 0) { - // Get the number of gas particles for this process + // Get the number of gas particles for this process and set + // stream position auto psize = ios_psize(header.nsph, sizeof(gas_particle)); // 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) { - // Get the number of dark particles for this process + // Get the number of dark particles for this process and set + // stream position auto psize = ios_psize(header.ndark, sizeof(dark_particle)); // 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) { - // Get the number of star particles for this process - int psize = ios_psize(header.nstar, sizeof(star_particle)); + // Get the number of star particles for this process and set + // stream position + auto psize = ios_psize(header.nstar, sizeof(star_particle)); // 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; diff --git a/pyEXP/UtilWrappers.cc b/pyEXP/UtilWrappers.cc index 04a80ae01..f04032bdd 100644 --- a/pyEXP/UtilWrappers.cc +++ b/pyEXP/UtilWrappers.cc @@ -173,4 +173,69 @@ 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 Date: Thu, 7 Aug 2025 16:39:41 -0400 Subject: [PATCH 05/15] Promote stream arithmetic to unsigned long --- include/tipsy.H | 34 +++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/include/tipsy.H b/include/tipsy.H index 32a6308c9..e7ba0f804 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -134,7 +134,7 @@ namespace TipsyReader XDR xdrs; //! Cache the XDR stream position at beginning of particle group - unsigned curpos = 0; + unsigned long curpos = 0; int xdr_header() { @@ -171,14 +171,14 @@ namespace TipsyReader //! Generate stream position for parallel read //! \param nsize The number of particles to read //! \param ssize The size of each particle structure in bytes - unsigned xdr_psize(unsigned int nsize, unsigned int ssize) + unsigned long xdr_psize(unsigned long nsize, unsigned long ssize) { // Number of particles per process - unsigned int psize = nsize/numprocs; + unsigned long psize = nsize/numprocs; // Get current stream position curpos = xdr_getpos(&xdrs); - auto pos = curpos + auto pos = curpos; // Move the stream position to the start of this process' data pos += psize * myid * ssize; @@ -193,7 +193,7 @@ namespace TipsyReader } //! Set stream position to next group of particles - void xdr_pnext(unsigned int nsize, unsigned int ssize) + void xdr_pnext(unsigned long nsize, unsigned long ssize) { // Move the stream position to the start of the next particle // group @@ -348,22 +348,24 @@ namespace TipsyReader //! Generate stream position for parallel read //! \param nsize The number of particles to read //! \param ssize The size of each particle structure in bytes - unsigned ios_psize(unsigned int nsize, unsigned int ssize) + unsigned long ios_psize(unsigned long nsize, unsigned long ssize) { // Number of particles per process - unsigned int psize = nsize/numprocs; + unsigned long psize = nsize/numprocs; // Get current stream position curpos = input.tellg(); - auto pos = curpos; + std::streampos pos(curpos); // Move the stream position to the start of this process' data pos += psize * myid * ssize; - input.seekg(pos); - if (input.fail()) { + + 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) psize = nsize - (numprocs-1)*psize; @@ -371,14 +373,14 @@ namespace TipsyReader } //! Move to next particle group - void ios_pnext(unsigned int nsize, unsigned int ssize) + void ios_pnext(unsigned long nsize, unsigned long ssize) { // Move the stream position to the start of the next particle // group - auto pos = curpos; // Make a copy of the current position - pos += nsize * ssize; // Augment the position - input.seekg(pos); - if (input.fail()) { + 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"); } } @@ -409,6 +411,8 @@ namespace TipsyReader // Do the read dark_particles.resize(psize); read_dark(); + // Stream position debug + std::streampos test = input.tellg(); N++; // Set stream position to next group of particles From 40ee042daa697c71f20619733d12a96896e2b343 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 7 Aug 2025 16:49:28 -0400 Subject: [PATCH 06/15] Add losts more comments [no CI] --- include/tipsy.H | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/include/tipsy.H b/include/tipsy.H index e7ba0f804..4b05c6b16 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -171,6 +171,8 @@ namespace TipsyReader //! 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) { // Number of particles per process @@ -193,6 +195,9 @@ namespace TipsyReader } //! 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 @@ -282,6 +287,7 @@ namespace TipsyReader public: + //! Constructor that opens the file and reads the header TipsyXDR(const std::string filename) { // Attempt to read tipsy file @@ -291,8 +297,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); @@ -326,6 +334,8 @@ namespace TipsyReader 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 { @@ -348,6 +358,8 @@ namespace TipsyReader //! 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) { // Number of particles per process @@ -373,6 +385,9 @@ namespace TipsyReader } //! 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 @@ -385,6 +400,8 @@ namespace TipsyReader } } + //! Read the particles from the file + //! \return The number of particle groups read int native_read() { int N=0; @@ -436,6 +453,7 @@ namespace TipsyReader return N; } + //! Read the gas particles from the file void read_gas() { try { @@ -449,6 +467,7 @@ namespace TipsyReader } } + //! Read the dark particles from the file void read_dark() { try { @@ -462,6 +481,7 @@ namespace TipsyReader } } + //! Read the star particles from the file void read_star() { try { @@ -479,6 +499,7 @@ namespace TipsyReader public: + //! Constructor that opens the file and reads the header TipsyNative(const std::string filename) { // Attempt to read tipsy file @@ -488,8 +509,10 @@ namespace TipsyReader } } + //! Read the particles from the file int readParticles() { return native_read(); } + //! Destructor that closes the file by scoping ~TipsyNative() { // Nothing From 0bd7607321bee45a62e9bc5b851ef96339e53754 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 7 Aug 2025 20:36:15 -0400 Subject: [PATCH 07/15] Added more explanation for setMPI() --- exputil/ParticleReader.cc | 10 +++++----- include/ParticleReader.H | 2 +- pyEXP/UtilWrappers.cc | 13 +++++++++---- 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/exputil/ParticleReader.cc b/exputil/ParticleReader.cc index ea63c0400..fa6b7baa6 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; @@ -2125,7 +2125,7 @@ namespace PR { } Tipsy::Tipsy(const std::string& file, TipsyType Type, - bool verbose) + bool verbose) : ParticleReader() { ttype = Type; files.push_back(file); @@ -2137,7 +2137,7 @@ namespace PR { } Tipsy::Tipsy(const std::vector& filelist, TipsyType Type, - bool verbose) + bool verbose) : ParticleReader() { files = filelist; diff --git a/include/ParticleReader.H b/include/ParticleReader.H index 0d150dc45..8a22cada5 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(); } diff --git a/pyEXP/UtilWrappers.cc b/pyEXP/UtilWrappers.cc index f04032bdd..541bf93fb 100644 --- a/pyEXP/UtilWrappers.cc +++ b/pyEXP/UtilWrappers.cc @@ -224,14 +224,19 @@ void UtilityClasses(py::module &m) { } }, R"( - Initialize EXP for MPI from Python. + Initialize EXP's C++ MPI communicators and variables from pyEXP - Sets internal MPI-specific library values necessary to MPI from Python scripts. - + Sets all of the internal MPI-specific library values necessary + to use the full complement of the EXP C++ classes from Python + scripts. For example, this is needed for ParticleReader, which + does not provide its own MPI initialization. The BasisFactory, + Basis, ParticleReader, and FieldGenerator classes can do this + automatically. + Parameters ---------- verbose : bool - If true, print the node and process assignments + Print the node and process assignments. Default: false. Returns ------- From 5e2d557f6c030cbfe8f19b9f96546e0998a12896 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 10 Aug 2025 08:53:48 -0400 Subject: [PATCH 08/15] Save native Tipsy particle offsets per node for creating unique sequential indices --- include/tipsy.H | 48 ++++++++++++++++++++++++++++++++++++------------ 1 file changed, 36 insertions(+), 12 deletions(-) diff --git a/include/tipsy.H b/include/tipsy.H index 4b05c6b16..fa36aadda 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -1,6 +1,8 @@ #ifndef _tipsy_H #define _tipsy_H +#include + // For MPI awareness // #include @@ -97,6 +99,8 @@ namespace TipsyReader #endif } ; + enum class Ptype { gas, dark, star }; + class TipsyFile { protected: @@ -105,6 +109,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: @@ -116,6 +123,11 @@ namespace TipsyReader virtual int readParticles() = 0; virtual ~TipsyFile() {} + + unsigned long getIndexOffset(const Ptype ptype) + { + return index_offset_map[ptype]; + } }; #ifdef HAVE_XDR @@ -173,7 +185,7 @@ namespace TipsyReader //! \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) + unsigned long xdr_psize(unsigned long nsize, unsigned long ssize, Ptype ptype) { // Number of particles per process unsigned long psize = nsize/numprocs; @@ -182,14 +194,20 @@ namespace TipsyReader 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 += psize * myid * ssize; + 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) psize = nsize - (numprocs-1)*psize; + if (myid == numprocs-1) { + index_offset_map[ptype] = (numprocs-1) * psize; + psize = nsize - index_offset; + }} return psize; } @@ -214,7 +232,7 @@ namespace TipsyReader if (header.nsph != 0) { // Get the number of gas particles for this process and set // stream position - auto psize = xdr_psize(header.nsph, sizeof(gas_particle)); + auto psize = xdr_psize(header.nsph, sizeof(gas_particle), Ptype::gas); // Do the read gas_particles.resize(psize); @@ -228,7 +246,7 @@ namespace TipsyReader if (header.ndark != 0) { // Get the number of dark particles for this process and set // stream position - auto psize = xdr_psize(header.ndark, sizeof(dark_particle)); + auto psize = xdr_psize(header.ndark, sizeof(dark_particle), Ptype::dark); // Do the read dark_particles.resize(psize); @@ -242,7 +260,7 @@ namespace TipsyReader if (header.nstar != 0) { // Get the number of star particles for this process and set // stream position - auto psize = xdr_psize(header.nstar, sizeof(star_particle)); + auto psize = xdr_psize(header.nstar, sizeof(star_particle), Ptype::star); // Do the read star_particles.resize(psize); @@ -360,7 +378,7 @@ namespace TipsyReader //! \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) + unsigned long ios_psize(unsigned long nsize, unsigned long ssize, Ptype ptype) { // Number of particles per process unsigned long psize = nsize/numprocs; @@ -369,8 +387,11 @@ namespace TipsyReader 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 += psize * myid * ssize; + pos += index_offset_map[ptype] * ssize; if (input.seekg(pos).fail()) { throw std::runtime_error("TipsyFile: seekg failed"); @@ -379,7 +400,10 @@ namespace TipsyReader std::streampos test = input.tellg(); // If this is the last process, read the remainder - if (myid == numprocs-1) psize = nsize - (numprocs-1)*psize; + if (myid == numprocs-1) { + index_offset_map[ptype] = (numprocs-1) * psize; + psize = nsize - index_offset_map[ptype]; + } return psize; } @@ -409,7 +433,7 @@ namespace TipsyReader if (header.nsph != 0) { // Get the number of gas particles for this process and set // stream position - auto psize = ios_psize(header.nsph, sizeof(gas_particle)); + auto psize = ios_psize(header.nsph, sizeof(gas_particle), Ptype::gas); // Do the read gas_particles.resize(psize); @@ -423,7 +447,7 @@ namespace TipsyReader if (header.ndark != 0) { // Get the number of dark particles for this process and set // stream position - auto psize = ios_psize(header.ndark, sizeof(dark_particle)); + auto psize = ios_psize(header.ndark, sizeof(dark_particle), Ptype::dark); // Do the read dark_particles.resize(psize); @@ -439,7 +463,7 @@ namespace TipsyReader if (header.nstar != 0) { // Get the number of star particles for this process and set // stream position - auto psize = ios_psize(header.nstar, sizeof(star_particle)); + auto psize = ios_psize(header.nstar, sizeof(star_particle), Ptype::star); // Do the read star_particles.resize(psize); From 85945c85b6fc8319bb71aacd6b2799e1a9e4883f Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 10 Aug 2025 09:03:46 -0400 Subject: [PATCH 09/15] Assign EXP index from native Tipsy beginning with 1 --- exputil/ParticleReader.cc | 6 ++++++ include/tipsy.H | 7 ++++++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/exputil/ParticleReader.cc b/exputil/ParticleReader.cc index fa6b7baa6..71e33d2f0 100644 --- a/exputil/ParticleReader.cc +++ b/exputil/ParticleReader.cc @@ -2185,6 +2185,8 @@ namespace PR { P.vel[k] = ps->gas_particles[pcount].vel[k]; } if (ttype == TipsyType::bonsai) P.indx = ps->gas_particles[pcount].ID(); + else P.indx = ps->getIndexOffset(TipsyReader::Ptype::gas) + pcount + 1; + pcount++; return; } @@ -2197,6 +2199,8 @@ namespace PR { P.vel[k] = ps->dark_particles[pcount].vel[k]; } if (ttype == TipsyType::bonsai) P.indx = ps->dark_particles[pcount].ID(); + else P.indx = ps->getIndexOffset(TipsyReader::Ptype::dark) + pcount + 1; + pcount++; return; } @@ -2209,6 +2213,8 @@ namespace PR { P.vel[k] = ps->star_particles[pcount].vel[k]; } if (ttype == TipsyType::bonsai) P.indx = ps->star_particles[pcount].ID(); + else P.indx = ps->getIndexOffset(TipsyReader::Ptype::star) + pcount + 1; + pcount++; return; } diff --git a/include/tipsy.H b/include/tipsy.H index fa36aadda..71c6917a6 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -99,8 +99,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: @@ -109,7 +111,7 @@ namespace TipsyReader virtual void read_dark() = 0; virtual void read_star() = 0; - // Offset for indexing particles in the file + //! Offset for indexing particles in the file std::map index_offset_map { {Ptype::gas, 0}, {Ptype::dark, 0}, {Ptype::star, 0} }; @@ -124,6 +126,7 @@ namespace TipsyReader virtual ~TipsyFile() {} + //! Get the index offset for a given particle type unsigned long getIndexOffset(const Ptype ptype) { return index_offset_map[ptype]; @@ -135,6 +138,7 @@ namespace TipsyReader #include #include + //! Tipsy XDR file reader class TipsyXDR : public TipsyFile { private: @@ -329,6 +333,7 @@ namespace TipsyReader #endif + //! Tipsy native file reader class TipsyNative : public TipsyFile { private: From 4eaa424eb1a2ed37e66b49a492bf4f8e0a8ead73 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 10 Aug 2025 10:46:17 -0400 Subject: [PATCH 10/15] Typo in Centering routine; missing count variable update --- expui/Centering.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/Centering.cc b/expui/Centering.cc index 8b00b1f07..4cff61094 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; From 9db73e1df6cb5bc5072b016a8145c40d553b6062 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 13 Aug 2025 10:28:54 -0400 Subject: [PATCH 11/15] Updated Bonsai reader variant to v2 --- exputil/ParticleReader.cc | 24 ++++++++++++++++-------- include/ParticleReader.H | 2 +- include/tipsy.H | 23 +++++++++++++++-------- 3 files changed, 32 insertions(+), 17 deletions(-) diff --git a/exputil/ParticleReader.cc b/exputil/ParticleReader.cc index 71e33d2f0..5a85dced2 100644 --- a/exputil/ParticleReader.cc +++ b/exputil/ParticleReader.cc @@ -1896,7 +1896,7 @@ namespace PR { std::vector ParticleReader::readerTypes {"PSPout", "PSPspl", "GadgetNative", "GadgetHDF5", "PSPhdf5", - "TipsyNative", "TipsyXDR", "Bonsai"}; + "TipsyNative", "TipsyXDR", "ChaNGa", "Bonsai"}; std::vector> @@ -2031,7 +2031,8 @@ namespace PR { exit(1); } #endif - + else if (reader.find("ChaNGa") == 0) + ret = std::make_shared(file, Tipsy::TipsyType::changa, verbose); else if (reader.find("Bonsai") == 0) ret = std::make_shared(file, Tipsy::TipsyType::bonsai, verbose); else { @@ -2184,8 +2185,7 @@ 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(); - else P.indx = ps->getIndexOffset(TipsyReader::Ptype::gas) + pcount + 1; + P.indx = ps->getIndexOffset(TipsyReader::Ptype::gas) + pcount + 1; pcount++; return; @@ -2198,8 +2198,12 @@ 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(); - else P.indx = ps->getIndexOffset(TipsyReader::Ptype::dark) + pcount + 1; + if (ttype == TipsyType::bonsai) + P.indx = ps->dark_particles[pcount].ID2(); + else if (ttype == TipsyType::changa) + P.indx = ps->dark_particles[pcount].ID(); + else + P.indx = ps->getIndexOffset(TipsyReader::Ptype::dark) + pcount + 1; pcount++; return; @@ -2212,8 +2216,12 @@ 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(); - else P.indx = ps->getIndexOffset(TipsyReader::Ptype::star) + pcount + 1; + if (ttype == TipsyType::bonsai) + P.indx = ps->star_particles[pcount].ID2(); + else if (ttype == TipsyType::changa) + P.indx = ps->star_particles[pcount].ID(); + else + P.indx = ps->getIndexOffset(TipsyReader::Ptype::star) + pcount + 1; pcount++; return; diff --git a/include/ParticleReader.H b/include/ParticleReader.H index 8a22cada5..8852504de 100644 --- a/include/ParticleReader.H +++ b/include/ParticleReader.H @@ -573,7 +573,7 @@ namespace PR public: - enum class TipsyType { native, xdr, bonsai }; + enum class TipsyType { native, xdr, changa, bonsai }; private: //! List of files diff --git a/include/tipsy.H b/include/tipsy.H index 71c6917a6..ee2420fc6 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -40,12 +40,6 @@ namespace TipsyReader Real metals ; Real phi ; //@} - - //! Convert phi to index for Bonsai - int ID() const { - union id {Real v; int i;} u; - u.v = phi; return u.i; - } } ; struct dark_particle @@ -59,11 +53,18 @@ namespace TipsyReader Real phi ; //@} - //! Convert phi to index for Bonsai + //! Convert phi to 32-bit index for ChaNGa int ID() const { union id {Real v; int i;} u; u.v = phi; return u.i; } + + //! Convert phi to 64-bit index for Bonsai + 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 @@ -79,11 +80,17 @@ namespace TipsyReader Real phi ; //@} - //! Convert phi to index for Bonsai + //! Convert phi to 32-bit index for ChaNGa int ID() const { union id {Real v; int i;} u; u.v = phi; return u.i; } + + //! Convert exp, phi to 64bit index for Bonsai + 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 From 83ef0bde28e7e3ec949943c175aa9d61e0437369 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 13 Aug 2025 13:08:53 -0400 Subject: [PATCH 12/15] Change Tipsy tag from ChaNGa to Bonsai1; since I can't figure out what ChaNGa is doing for particle ID. --- exputil/ParticleReader.cc | 14 +++++++++----- include/ParticleReader.H | 2 +- include/tipsy.H | 8 ++++---- 3 files changed, 14 insertions(+), 10 deletions(-) diff --git a/exputil/ParticleReader.cc b/exputil/ParticleReader.cc index 5a85dced2..76c68066c 100644 --- a/exputil/ParticleReader.cc +++ b/exputil/ParticleReader.cc @@ -1896,7 +1896,7 @@ namespace PR { std::vector ParticleReader::readerTypes {"PSPout", "PSPspl", "GadgetNative", "GadgetHDF5", "PSPhdf5", - "TipsyNative", "TipsyXDR", "ChaNGa", "Bonsai"}; + "TipsyNative", "TipsyXDR", "Bonsai1", "Bonsai"}; std::vector> @@ -2031,8 +2031,8 @@ namespace PR { exit(1); } #endif - else if (reader.find("ChaNGa") == 0) - ret = std::make_shared(file, Tipsy::TipsyType::changa, verbose); + 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 { @@ -2072,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 @@ -2111,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); @@ -2200,7 +2204,7 @@ namespace PR { } if (ttype == TipsyType::bonsai) P.indx = ps->dark_particles[pcount].ID2(); - else if (ttype == TipsyType::changa) + else if (ttype == TipsyType::bonsai1) P.indx = ps->dark_particles[pcount].ID(); else P.indx = ps->getIndexOffset(TipsyReader::Ptype::dark) + pcount + 1; @@ -2218,7 +2222,7 @@ namespace PR { } if (ttype == TipsyType::bonsai) P.indx = ps->star_particles[pcount].ID2(); - else if (ttype == TipsyType::changa) + else if (ttype == TipsyType::bonsai1) P.indx = ps->star_particles[pcount].ID(); else P.indx = ps->getIndexOffset(TipsyReader::Ptype::star) + pcount + 1; diff --git a/include/ParticleReader.H b/include/ParticleReader.H index 8852504de..8dc9a851e 100644 --- a/include/ParticleReader.H +++ b/include/ParticleReader.H @@ -573,7 +573,7 @@ namespace PR public: - enum class TipsyType { native, xdr, changa, bonsai }; + enum class TipsyType { native, xdr, bonsai1, bonsai }; private: //! List of files diff --git a/include/tipsy.H b/include/tipsy.H index ee2420fc6..06d37403a 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -53,13 +53,13 @@ namespace TipsyReader Real phi ; //@} - //! Convert phi to 32-bit index for ChaNGa + //! 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 + //! 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; @@ -80,13 +80,13 @@ namespace TipsyReader Real phi ; //@} - //! Convert phi to 32-bit index for ChaNGa + //! 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 + //! 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; From b609af3055f310aeb5dbb6d4d73b618ec3b2df7e Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 17 Aug 2025 09:45:05 -0400 Subject: [PATCH 13/15] Update include/tipsy.H Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- include/tipsy.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/tipsy.H b/include/tipsy.H index 06d37403a..655325bc1 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -217,7 +217,7 @@ namespace TipsyReader // If this is the last process, read the remainder if (myid == numprocs-1) { index_offset_map[ptype] = (numprocs-1) * psize; - psize = nsize - index_offset; + psize = nsize - index_offset_map[ptype]; }} return psize; From b58315671b8914450053ad8f0563d501a515b5dc Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 17 Aug 2025 09:46:11 -0400 Subject: [PATCH 14/15] Update expui/Centering.cc Sure. I agree that pre-increment better represents the intent. Although post-increment is not wrong either. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/Centering.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/Centering.cc b/expui/Centering.cc index 4cff61094..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; From c9a8370c334c9ed9568b9b9762376682e96f5e2d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 17 Aug 2025 10:01:45 -0400 Subject: [PATCH 15/15] Remove debug cruft --- include/tipsy.H | 2 -- 1 file changed, 2 deletions(-) diff --git a/include/tipsy.H b/include/tipsy.H index 655325bc1..fefe2bcbc 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -464,8 +464,6 @@ namespace TipsyReader // Do the read dark_particles.resize(psize); read_dark(); - // Stream position debug - std::streampos test = input.tellg(); N++; // Set stream position to next group of particles