From caa6a1925cb155914e8b8c9042c2c70dd9410435 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 14 Aug 2026 12:14:20 -0400 Subject: [PATCH 01/28] Added new ParticleHDF class to generalize EXP HDF5 particle writing and applied this to all of the native particle generators. --- tests/CMakeLists.txt | 31 ++++++++++++++++++++++ utils/ICs/cubeICs.cc | 49 +++++++++++++++++++++++++---------- utils/ICs/genslab.cc | 40 +++++++++++++++++++++++++++- utils/ICs/gensph.cc | 37 +++++++++++++++++++++----- utils/ICs/initial.cc | 59 +++++++++++++++++++++++++++++------------- utils/ICs/initial2d.cc | 23 +++++++++++++--- 6 files changed, 196 insertions(+), 43 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 45225dc3e..afd7a2807 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -2,6 +2,16 @@ include(CTest) set(CTEST_OUTPUT_ON_FAILURE ON) +# HDF5 particle IC smoke checks use h5py for concise schema inspection. Keep +# them optional so a Python installation without h5py does not disable CTest. +find_package(Python3 COMPONENTS Interpreter QUIET) +if(Python3_Interpreter_FOUND) + execute_process( + COMMAND ${Python3_EXECUTABLE} -c "import h5py, numpy" + RESULT_VARIABLE EXP_H5PY_STATUS + OUTPUT_QUIET ERROR_QUIET) +endif() + # Is EXP configured for pyEXP? If yes, run pyEXP tests... # if(ENABLE_PYEXP) @@ -150,6 +160,27 @@ if(ENABLE_NBODY) set_tests_properties(removeCubeFiles PROPERTIES DEPENDS expCubeCheckPos REQUIRED_FILES "config.runS.yml;current.processor.rates.runS;cube.bods;OUTLOG.runS;runS.levels;") + if(Python3_Interpreter_FOUND AND EXP_H5PY_STATUS EQUAL 0) + add_test(NAME makeCubeHDF5ICTest + COMMAND ${CMAKE_BINARY_DIR}/utils/ICs/cubeics --hdf5 --number 16 + --seed 17 --zerovel --file hdf5-cube.bods + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube) + + add_test(NAME checkCubeHDF5ICTest + COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py + hdf5-cube.bods.h5 --count 16 + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube) + set_tests_properties(checkCubeHDF5ICTest PROPERTIES DEPENDS makeCubeHDF5ICTest) + + add_test(NAME removeCubeHDF5ICFiles + COMMAND ${CMAKE_COMMAND} -E remove hdf5-cube.bods.h5 + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube) + set_tests_properties(removeCubeHDF5ICFiles PROPERTIES DEPENDS checkCubeHDF5ICTest + REQUIRED_FILES "hdf5-cube.bods.h5") + set_tests_properties(makeCubeHDF5ICTest checkCubeHDF5ICTest removeCubeHDF5ICFiles + PROPERTIES LABELS "quick") + endif() + # Set labels for pyEXP tests set_tests_properties(expExecuteTest PROPERTIES LABELS "quick") set_tests_properties(makeICTest expNbodyTest expNbodyCheck2TW diff --git a/utils/ICs/cubeICs.cc b/utils/ICs/cubeICs.cc index fff4586ed..787fee237 100644 --- a/utils/ICs/cubeICs.cc +++ b/utils/ICs/cubeICs.cc @@ -15,6 +15,7 @@ #include #include "cxxopts.H" +#include "ParticleHDF5.H" int main(int ac, char **av) @@ -28,6 +29,8 @@ main(int ac, char **av) std::string bodyfile; // Output file unsigned seed; // Will be inialized by /dev/random if // not set on the command line + unsigned hdf5_filter = 1; + bool hdf5_output = false, hdf5_double = false; // Default values for the velocity dispersion and bulk velocity // @@ -55,6 +58,10 @@ main(int ac, char **av) cxxopts::value(seed)) ("o,file", "Output body file", cxxopts::value(bodyfile)->default_value("cube.bods")) + ("5,hdf5", "Write HDF5 phase-space output instead of ASCII") + ("8,double", "Use float64 HDF5 output (default is float32)") + ("f,filter", "HDF5 filter ID (default: 1 = GZIP)", + cxxopts::value(hdf5_filter)->default_value("1")) ("w,wave", "Perturbation wave vector", cxxopts::value>(pert)) ("p,pamp", "Perturbation amplitude", @@ -76,6 +83,8 @@ main(int ac, char **av) std::cout << options.help() << std::endl << std::endl; return 1; } + hdf5_output = vm.count("hdf5") > 0; + hdf5_double = vm.count("double") > 0; // Set from /dev/random if not specified if (vm.count("seed")==0) { @@ -122,13 +131,15 @@ main(int ac, char **av) // Open the output file // - std::ofstream out(bodyfile); + std::ofstream out; + if (!hdf5_output) out.open(bodyfile); - if (out) { + if (hdf5_output || out) { // Header // - out << std::setw(10) << N << std::setw(10) << 0 << std::setw(10) << 0 - << std::endl << std::setprecision(10); + if (!hdf5_output) + out << std::setw(10) << N << std::setw(10) << 0 << std::setw(10) << 0 + << std::endl << std::setprecision(10); std::mt19937 gen(seed); std::uniform_real_distribution<> uniform(0.0, 1.0); @@ -184,19 +195,31 @@ main(int ac, char **av) double mass = M/N; - for (int i=0; i> particles; + particles.reserve(N); + for (int i=0; i(bodyfile + ".h5", particles, 0, 0, hdf5_filter); + else + EXP::ParticleHDF5::write(bodyfile + ".h5", particles, 0, 0, hdf5_filter); + } else { + for (int i=0; i +struct ParticleData +{ + std::optional> index; // optional particle index + std::vector m; // mass + std::vector x, y, z; // position + std::vector u, v, w; // velocity + std::vector> aux_ints; // auxiliary integer fields + std::vector> aux_floats; // auxiliary float fields + + size_t num_particles = 0; + size_t num_aux_ints = 0; + size_t num_aux_floats = 0; +}; + +// Write particle data to HDF5 file with specified filter and precision +template +void write_hdf5_data(const std::string& outfile, + const ParticleData& data, + unsigned filter_id, + bool verbose) +{ + std::string hdf5_file = outfile + ".h5"; + EXP::ParticleHDF5::write(hdf5_file, data.m, data.x, data.y, data.z, + data.u, data.v, data.w, data.aux_ints, + data.aux_floats, data.index, filter_id); + + if (verbose) { + std::string precision_str = "float32"; + if (sizeof(T) == 8) precision_str = "float64"; + std::cout << "Successfully wrote " << data.num_particles << " particles to " + << hdf5_file << " (" << precision_str << ")" << std::endl; + } +} int main(int argc, char **argv) @@ -214,4 +253,3 @@ main(int argc, char **argv) << " 2T/VC=" << KE/VC << std::endl; } - diff --git a/utils/ICs/gensph.cc b/utils/ICs/gensph.cc index 4d72a2e70..37f373713 100644 --- a/utils/ICs/gensph.cc +++ b/utils/ICs/gensph.cc @@ -49,6 +49,7 @@ #include "fpetrap.h" #include "cxxopts.H" #include "EXPini.H" +#include "ParticleHDF5.H" // Reference to n-body globals // @@ -76,6 +77,8 @@ main(int argc, char **argv) double Emin0, Emax0, Kmin0, Kmax0, RBAR, MBAR, BRATIO, CRATIO, SMOOTH; bool LOGR, ELIMIT, VERBOSE, GRIDPOT, MODELS, EBAR, zeropos, zerovel; bool VTEST; + bool hdf5_output = false, hdf5_double = false; + unsigned hdf5_filter = 1; std::string INFILE, MMFILE, OUTFILE, OUTPS, config; const double goldenRatio = 1.618033988749895; @@ -108,6 +111,10 @@ main(int argc, char **argv) cxxopts::value(MMFILE)) ("o,PSFILE", "Phase-space output file", cxxopts::value(OUTPS)->default_value("new.bods")) + ("5,hdf5", "Write HDF5 phase-space output instead of ASCII") + ("8,double", "Use float64 HDF5 output (default is float32)") + ("f,filter", "HDF5 filter ID (default: 1 = GZIP)", + cxxopts::value(hdf5_filter)->default_value("1")) ("p,prefix", "Diagnostic output file prefix", cxxopts::value(OUTFILE)->default_value("gensph")) ("zeropos", "Set the origin at the center of mass", @@ -259,6 +266,8 @@ main(int argc, char **argv) return 0; } } + hdf5_output = vm.count("hdf5") > 0; + hdf5_double = vm.count("double") > 0; if (vm.count("verbose")) VERBOSE = true; else VERBOSE = false; @@ -272,9 +281,10 @@ main(int argc, char **argv) // std::ostringstream sout; sout << OUTPS << "." << myid; - std::ofstream out(sout.str()); + std::ofstream out; int bad = 0; - if (!out) { + if (!hdf5_output) out.open(sout.str()); + if (!hdf5_output && !out) { std::cerr << "[" << myid << "] Couldn't open <" << sout.str() << "> for output" << std::endl; bad = 1; @@ -285,8 +295,10 @@ main(int argc, char **argv) exit(-1); } - out.precision(11); - out.setf(ios::scientific, ios::floatfield); + if (!hdf5_output) { + out.precision(11); + out.setf(ios::scientific, ios::floatfield); + } // Begin integration // @@ -701,6 +713,7 @@ main(int argc, char **argv) int ierr; Eigen::VectorXd ps(7), ps0(7); + std::vector> hdf5_particles; ps0[0] = 0.0; ps0[1] = X0; @@ -710,7 +723,7 @@ main(int argc, char **argv) ps0[5] = V0; ps0[6] = W0; - if (myid==0) { + if (!hdf5_output && myid==0) { out << std::setw(12) << N << std::setw( 6) << NI << std::setw( 6) << ND << std::endl; } @@ -852,7 +865,10 @@ main(int argc, char **argv) if (zeropos) for (int i=1; i<3; i++) zz[i] -= ps[0]*ps[i]; if (zerovel) for (int i=4; i<7; i++) zz[i] -= ps[0]*ps[i]; } - else { + else if (hdf5_output) { + hdf5_particles.push_back({mass * ps[0], ps[1] + ps0[1], ps[2] + ps0[2], ps[3] + ps0[3], + ps[4] + ps0[4], ps[5] + ps0[5], ps[6] + ps0[6]}); + } else { out << std::setw(20) << mass * ps[0]; for (int i=1; i<=6; i++) out << std::setw(20) << ps[i]+ps0[i]; @@ -919,6 +935,11 @@ main(int argc, char **argv) } for (auto ps : PS) { + if (hdf5_output) { + hdf5_particles.push_back({ps[0], ps[1] + ps0[1], ps[2] + ps0[2], ps[3] + ps0[3], + ps[4] + ps0[4], ps[5] + ps0[5], ps[6] + ps0[6]}); + continue; + } out << std::setw(20) << ps[0]; for (int i=1; i<=6; i++) out << std::setw(20) << ps[i]+ps0[i]; @@ -966,11 +987,13 @@ main(int argc, char **argv) } out.close(); + if (hdf5_output) + EXP::ParticleHDF5::gather_and_write(OUTPS + ".h5", hdf5_particles, NI, ND, hdf5_filter, hdf5_double); MPI_Barrier(MPI_COMM_WORLD); // Make the final phase-space file and clean up // - if (myid==0) { + if (!hdf5_output && myid==0) { std::ostringstream sout; sout << "cat " << OUTPS << ".* > " << OUTPS; sout << "; rm " << OUTPS << ".*"; diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index c94929b46..8e480078a 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -103,6 +103,7 @@ #include "libvars.H" // Library globals #include "cxxopts.H" // Command-line parsing #include "EXPini.H" // Ini-style config +#include "ParticleHDF5.H" #include "norminv.H" @@ -403,6 +404,8 @@ main(int ac, char **av) int nhalo, ndisk, ngas, ngparam; std::string hbods, dbods, gbods, outtag, runtag, centerfile, halofile1, halofile2; std::string cachefile, config, gentype, dtype, dmodel, mtype, ctype; + unsigned hdf5_filter = 1; + bool hdf5_output = false, hdf5_double = false; const std::string mesg("Generates a Monte Carlo realization of a halo with an\n embedded disk using Jeans' equations\n"); @@ -427,6 +430,10 @@ main(int ac, char **av) cxxopts::value(gbods)->default_value("gas.bods")) ("dbods", "The output bodyfile for the stellar disc", cxxopts::value(dbods)->default_value("disk.bods")) + ("5,hdf5", "Write HDF5 phase-space output instead of ASCII") + ("8,double", "Use float64 HDF5 output (default is float32)") + ("f,filter", "HDF5 filter ID (default: 1 = GZIP)", + cxxopts::value(hdf5_filter)->default_value("1")) ("cachefile", "The cache file for the cylindrical basis", cxxopts::value(cachefile)->default_value(".eof.cache.file")) ("ctype", "DiskHalo radial coordinate scaling type (one of: Linear, Log,Rat)", @@ -670,6 +677,8 @@ main(int ac, char **av) return 0; } } + hdf5_output = vm.count("hdf5") > 0; + hdf5_double = vm.count("double") > 0; if (vm.count("spline")) { SphericalModelTable::linear = 0; @@ -1067,7 +1076,7 @@ main(int ac, char **av) // before realizing a large phase space) std::ofstream out_halo, out_disk; if (myid==0) { - if (not evolved and n_particlesH) { + if (!hdf5_output && not evolved && n_particlesH) { out_halo.open(hbods); if (!out_halo) { cout << "Could not open <" << hbods << "> for output\n"; @@ -1076,7 +1085,7 @@ main(int ac, char **av) } } - if (ndisk) { + if (!hdf5_output && ndisk) { out_disk.open(dbods); if (!out_disk) { std::cout << "Could not open <" << dbods << "> for output" << std::endl; @@ -1477,13 +1486,19 @@ main(int ac, char **av) if (not evolved) { if (myid==0) std::cout << "Writing phase space file for halo . . . " << std::flush; - diskhalo->write_file(out_halo, hparticles); + if (hdf5_output) + EXP::ParticleHDF5::gather_and_write(hbods + ".h5", EXP::ParticleHDF5::records(hparticles), 0, 0, hdf5_filter, hdf5_double); + else + diskhalo->write_file(out_halo, hparticles); if (myid==0) std::cout << "done" << std::endl; out_halo.close(); } if (myid==0) std::cout << "Writing phase space file for disk . . . " << std::flush; - diskhalo->write_file(out_disk, dparticles); + if (hdf5_output) + EXP::ParticleHDF5::gather_and_write(dbods + ".h5", EXP::ParticleHDF5::records(dparticles), 0, 0, hdf5_filter, hdf5_double); + else + diskhalo->write_file(out_disk, dparticles); if (myid==0) std::cout << "done" << std::endl; out_disk.close(); // Diagnostic . . . @@ -1693,9 +1708,11 @@ main(int ac, char **av) // // Prepare output stream // - ofstream outps("gas.bods"); - if (!outps) { - cerr << "Couldn't open <" << "gas.bods" << "> for output\n"; + ofstream outps; + std::vector> gas_particles; + if (!hdf5_output) outps.open(gbods); + if (!hdf5_output && !outps) { + cerr << "Couldn't open <" << gbods << "> for output\n"; exit (-1); } @@ -1726,8 +1743,8 @@ main(int ac, char **av) gmass = gmass0; fr = fz = potr = 0.0; - outps << setw(8) << ngas - << setw(6) << 0 << setw(6) << ngparam << endl; + if (!hdf5_output) + outps << setw(8) << ngas << setw(6) << 0 << setw(6) << ngparam << endl; for (int n=0; naccumulated_eval(R, z, phi, p0, p, fr, fz, fp); @@ -1828,6 +1844,13 @@ main(int ac, char **av) std::cout << endl << "Done!" << std::endl; + if (hdf5_output) { + if (hdf5_double) + EXP::ParticleHDF5::write(gbods + ".h5", gas_particles, 0, ngparam, hdf5_filter); + else + EXP::ParticleHDF5::write(gbods + ".h5", gas_particles, 0, ngparam, hdf5_filter); + } + std::cout << "****************************" << std::endl << " Gas disk" << std::endl << "----------------------------" << std::endl diff --git a/utils/ICs/initial2d.cc b/utils/ICs/initial2d.cc index 5acbda48b..588d2f7ab 100644 --- a/utils/ICs/initial2d.cc +++ b/utils/ICs/initial2d.cc @@ -50,6 +50,7 @@ #include "libvars.H" // Library globals #include "cxxopts.H" // Command-line parsing #include "EXPini.H" // Ini-style config +#include "ParticleHDF5.H" #include "norminv.H" @@ -209,6 +210,8 @@ main(int ac, char **av) std::string hbods, dbods, suffix, centerfile, halofile1, halofile2; std::string cachefile, config, gentype, dtype, dmodel, mtype, ctype; std::string diskconf; + unsigned hdf5_filter = 1; + bool hdf5_output = false, hdf5_double = false; const std::string mesg("Generates a Monte Carlo realization of a halo with an\n embedded disk using Jeans' equations\n"); @@ -223,6 +226,10 @@ main(int ac, char **av) cxxopts::value(hbods)->default_value("halo.bods")) ("dbods", "The output bodyfile for the stellar disc", cxxopts::value(dbods)->default_value("disk.bods")) + ("5,hdf5", "Write HDF5 phase-space output instead of ASCII") + ("8,double", "Use float64 HDF5 output (default is float32)") + ("f,filter", "HDF5 filter ID (default: 1 = GZIP)", + cxxopts::value(hdf5_filter)->default_value("1")) ("cachefile", "The cache file for the cylindrical basis", cxxopts::value(cachefile)->default_value(".eof_2d_cache")) ("ctype", "DiskHalo radial coordinate scaling type (one of: Linear, Log,Rat)", @@ -429,6 +436,8 @@ main(int ac, char **av) return 0; } } + hdf5_output = vm.count("hdf5") > 0; + hdf5_double = vm.count("double") > 0; if (vm.count("spline")) { SphericalModelTable::linear = 0; @@ -688,7 +697,7 @@ main(int ac, char **av) // before realizing a large phase space) std::ofstream out_halo, out_disk; if (myid==0) { - if (not evolved and n_particlesH) { + if (!hdf5_output && not evolved && n_particlesH) { out_halo.open(hbods); if (!out_halo) { cout << "Could not open <" << hbods << "> for output\n"; @@ -697,7 +706,7 @@ main(int ac, char **av) } } - if (ndisk) { + if (!hdf5_output && ndisk) { out_disk.open(dbods); if (!out_disk) { std::cout << "Could not open <" << dbods << "> for output" << std::endl; @@ -927,13 +936,19 @@ main(int ac, char **av) if (not evolved) { if (myid==0) std::cout << "Writing phase space file for halo . . . " << std::flush; - diskhalo->write_file(out_halo, hparticles); + if (hdf5_output) + EXP::ParticleHDF5::gather_and_write(hbods + ".h5", EXP::ParticleHDF5::records(hparticles), 0, 0, hdf5_filter, hdf5_double); + else + diskhalo->write_file(out_halo, hparticles); if (myid==0) std::cout << "done" << std::endl; out_halo.close(); } if (myid==0) std::cout << "Writing phase space file for disk . . . " << std::flush; - diskhalo->write_file(out_disk, dparticles); + if (hdf5_output) + EXP::ParticleHDF5::gather_and_write(dbods + ".h5", EXP::ParticleHDF5::records(dparticles), 0, 0, hdf5_filter, hdf5_double); + else + diskhalo->write_file(out_disk, dparticles); if (myid==0) std::cout << "done" << std::endl; out_disk.close(); // Diagnostic . . . From 0b48448c753696ccaa08b237dc91b487d44f7383 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 14 Aug 2026 12:25:22 -0400 Subject: [PATCH 02/28] Added new ParticleHDF5 class header --- utils/ICs/ParticleHDF5.H | 172 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 172 insertions(+) create mode 100644 utils/ICs/ParticleHDF5.H diff --git a/utils/ICs/ParticleHDF5.H b/utils/ICs/ParticleHDF5.H new file mode 100644 index 000000000..dc7fe26d7 --- /dev/null +++ b/utils/ICs/ParticleHDF5.H @@ -0,0 +1,172 @@ +#ifndef EXP_PARTICLE_HDF5_H +#define EXP_PARTICLE_HDF5_H + +// Shared writer for the simple EXP phase-space HDF5 schema. The +// schema unchanged from the Component update and is the same one tested +// by the hdf5bods and slabics implementations. + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace EXP +{ + + namespace ParticleHDF5 + { + + //! Collects the mass, position, and velocity of a particle container into a vector of 7-element arrays. + template + inline std::vector> records(const Particles& particles) + { + std::vector> result; + result.reserve(particles.size()); + for (const auto& p : particles) + result.push_back({p.mass, p.pos[0], p.pos[1], p.pos[2], + p.vel[0], p.vel[1], p.vel[2]}); + return result; + } + + //! Create a HighFive::DataSetCreateProps object with chunking and optional filter. + inline HighFive::DataSetCreateProps properties(hsize_t size, unsigned filter) + { + if (size == 0) throw std::invalid_argument("cannot write an empty particle file"); + + HighFive::DataSetCreateProps props; + // chunk size is min(size, 256kB) but at least 1kB + props.add(HighFive::Chunking({std::min(size, std::min(262144, std::max(1024, size)))})); + if (filter != 3 && filter != 4 && filter != 32001) props.add(HighFive::Shuffle()); + + // The filter IDs are from HDF5 1.12.0 and later, see https://support.hdfgroup.org/HDF5/doc/Advanced/Filters.html + // Filter types: + // 1: GZIP (DEFLATE) + // 2: SZIP + // 3: SHUFFLE + // 4: FLETCHER32 + // 307: LZ4 + // 32004: ZSTD + // 32001: BLOSC + // + // The filter options are specific to each filter type, see the HDF5 documentation for details. + // + hid_t plist = props.getId(); + switch (filter) { + case 1: + case 307: { const unsigned level = 5; H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 1, &level); break; } + case 2: { const unsigned options[] = {141, 16}; H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 2, options); break; } + case 3: return props; // shuffle only + case 4: H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 0, nullptr); return props; + case 32004: { const unsigned options[] = {0}; H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 1, options); break; } + case 32001: { const unsigned options[] = {0, 0, 0, 0, 5, 1, 1}; H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 7, options); break; } + default: + throw std::invalid_argument("unsupported HDF5 filter ID"); + } + return props; + } + + //! Write a particle phase-space file with optional auxiliary integer and floating-point data. + template + inline void write(const std::string& filename, + const std::vector& m, const std::vector& x, + const std::vector& y, const std::vector& z, + const std::vector& u, const std::vector& v, + const std::vector& w, + const std::vector>& aux_ints, + const std::vector>& aux_floats, + const std::optional>& index, + unsigned filter) + { + const auto count = m.size(); + if (count == 0) return; + if (x.size() != count || y.size() != count || z.size() != count || + u.size() != count || v.size() != count || w.size() != count) + throw std::invalid_argument("inconsistent phase-space vector sizes"); + + HighFive::File file(filename, HighFive::File::ReadWrite | HighFive::File::Create | HighFive::File::Truncate); + + const int n = static_cast(count); + file.createAttribute("num_particles", HighFive::DataSpace::From(n)).write(n); + + const int num_aux_ints = static_cast(aux_ints.size()); + const int num_aux_floats = static_cast(aux_floats.size()); + file.createAttribute("num_aux_ints", HighFive::DataSpace::From(num_aux_ints)).write(num_aux_ints); + file.createAttribute("num_aux_floats", HighFive::DataSpace::From(num_aux_floats)).write(num_aux_floats); + + auto group = file.createGroup("particles"); + auto props = properties(count, filter); + group.createDataSet("m", m, props); group.createDataSet("x", x, props); group.createDataSet("y", y, props); + group.createDataSet("z", z, props); group.createDataSet("u", u, props); group.createDataSet("v", v, props); + group.createDataSet("w", w, props); + if (index) group.createDataSet("index", *index, props); + for (size_t j = 0; j < aux_ints.size(); ++j) + group.createDataSet("aux_int_" + std::to_string(j), aux_ints[j], props); + for (size_t j = 0; j < aux_floats.size(); ++j) + group.createDataSet("aux_float_" + std::to_string(j), aux_floats[j], props); + } + + //! Write a particle phase-space file from a vector of 7-element arrays with optional auxiliary integer and floating-point data. + template + inline void write(const std::string& filename, + const std::vector>& particles, + int num_aux_ints, int num_aux_floats, unsigned filter) + { + std::vector m, x, y, z, u, v, w; + m.reserve(particles.size()); x.reserve(particles.size()); y.reserve(particles.size()); z.reserve(particles.size()); + u.reserve(particles.size()); v.reserve(particles.size()); w.reserve(particles.size()); + for (const auto& p : particles) { + m.push_back(p[0]); x.push_back(p[1]); y.push_back(p[2]); z.push_back(p[3]); + u.push_back(p[4]); v.push_back(p[5]); w.push_back(p[6]); + } + std::vector> auxi(num_aux_ints, std::vector(particles.size(), 0)); + std::vector> auxf(num_aux_floats, std::vector(particles.size(), 0)); + write(filename, m, x, y, z, u, v, w, auxi, auxf, std::nullopt, filter); + } + + //! Collect rank-local phase space records and have rank zero write one file. + inline void gather_and_write(const std::string& filename, + const std::vector>& local, + int num_aux_ints, int num_aux_floats, + unsigned filter, bool double_precision, + MPI_Comm comm = MPI_COMM_WORLD) + { + int rank, ranks, local_count = static_cast(local.size()); + MPI_Comm_rank(comm, &rank); MPI_Comm_size(comm, &ranks); + std::vector counts(rank == 0 ? ranks : 0); + MPI_Gather(&local_count, 1, MPI_INT, rank == 0 ? counts.data() : nullptr, 1, MPI_INT, 0, comm); + + std::vector send(7 * local.size()); + for (size_t i = 0; i < local.size(); ++i) + std::copy(local[i].begin(), local[i].end(), send.begin() + 7*i); + std::vector counts7, offsets; + std::vector received; + if (rank == 0) { + counts7.resize(ranks); offsets.resize(ranks); + int total = 0; + for (int i = 0; i < ranks; ++i) { offsets[i] = total; counts7[i] = 7 * counts[i]; total += counts7[i]; } + received.resize(total); + } + MPI_Gatherv(send.data(), static_cast(send.size()), MPI_DOUBLE, + rank == 0 ? received.data() : nullptr, + rank == 0 ? counts7.data() : nullptr, rank == 0 ? offsets.data() : nullptr, + MPI_DOUBLE, 0, comm); + if (rank != 0) return; + + std::vector> all(received.size()/7); + for (size_t i = 0; i < all.size(); ++i) + std::copy(received.begin() + 7*i, received.begin() + 7*(i+1), all[i].begin()); + if (double_precision) write(filename, all, num_aux_ints, num_aux_floats, filter); + else write(filename, all, num_aux_ints, num_aux_floats, filter); + } + } // namespace ParticleHDF5 +} // namespace EXP + +#endif From 369a77856dd7137235359eca6d37b691a735a0f7 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 14 Aug 2026 13:46:53 -0400 Subject: [PATCH 03/28] Cherry pick the refactoring of 1d mass models from slabUpdate --- exputil/CMakeLists.txt | 2 +- {utils/ICs => exputil}/massmodel1d.cc | 209 ++++- include/massmodel1d.H | 1072 +++++++++++++++++++++++++ utils/ICs/massmodel1d.H | 442 ---------- 4 files changed, 1241 insertions(+), 484 deletions(-) rename {utils/ICs => exputil}/massmodel1d.cc (58%) mode change 100755 => 100644 create mode 100644 include/massmodel1d.H delete mode 100755 utils/ICs/massmodel1d.H diff --git a/exputil/CMakeLists.txt b/exputil/CMakeLists.txt index ba785f48c..191d0f7f8 100644 --- a/exputil/CMakeLists.txt +++ b/exputil/CMakeLists.txt @@ -28,7 +28,7 @@ set(PHASE_SRC phase.cc ensemble.cc io_ensemble.cc move_ensemble.cc set(SPECFUNC_SRC gammln.cc bessel.cc OrthoPoly.cc CauchyPV.cc) # modbessel.cc set(INTERP_SRC Spline.cc SplintE.cc Vodd2.cc Vlocate.cc levsurf.cc Interp1d.cc Cheby1d.cc MonotCubicInterpolator.cc) set(MASSMODEL_SRC massmodel.cc massmodel_dist.cc embedded.cc isothermal.cc realize_model.cc GenPoly.cc mestel.cc - toomre.cc exponential.cc) + massmodel1d.cc toomre.cc exponential.cc) set(ORBIT_SRC orbit.cc orbit_trans.cc FindOrb.cc) set(BIORTH_SRC biorth_wake.cc biorth.cc biorth2d.cc biorth_grid.cc diff --git a/utils/ICs/massmodel1d.cc b/exputil/massmodel1d.cc old mode 100755 new mode 100644 similarity index 58% rename from utils/ICs/massmodel1d.cc rename to exputil/massmodel1d.cc index 81e66db04..a9f22294b --- a/utils/ICs/massmodel1d.cc +++ b/exputil/massmodel1d.cc @@ -1,50 +1,19 @@ -/***************************************************************************** - * Description: - * ----------- - * - * These routines computes massmodels for 1-d slab - * - * - * Call sequence: - * ------------- - * - * Parameters: - * ---------- - * - * x as above - * - * Returns: - * ------- - * - * Value - * - * Notes: - * ----- - * - * - * By: - * -- - * - * MDW 11/13/88 - * - ***************************************************************************/ - - #include -#include +#include #include #include -#include #include "gaussQ.H" #include "interp.H" #include "massmodel1d.H" +#include "biorth1d.H" #include "RK4.H" -double Sech2::HMAX=1.0e6; -double Sech2mu::HMAX=1.0e6; -double Sech2Halo::HMAX=3.0e1; +double Sech2::HMAX=10.0; +double Sech2u::HMAX=10.0; +double Sech2mu::HMAX=10.0; +double Sech2Halo::HMAX=30.0; OneDModelTable::OneDModelTable(std::string filename, int PARM) { @@ -64,7 +33,7 @@ OneDModelTable::OneDModelTable(std::string filename, int PARM) // Read header in.getline((char *)line, 144); - while (string(line).find_first_of("!#") != string::npos) + while (std::string(line).find_first_of("!#") != std::string::npos) in.getline((char *)line, 144); @@ -76,7 +45,7 @@ OneDModelTable::OneDModelTable(std::string filename, int PARM) for (int i=0; i> z[i]; iline >> d[i]; @@ -87,8 +56,8 @@ OneDModelTable::OneDModelTable(std::string filename, int PARM) if (PARM) { in.getline(line, 144); - istringstream ins(line); - string word; + std::istringstream ins(line); + std::string word; while (1) { ins >> word; if (ins.good()) params.push_back(atof(word.c_str())); @@ -365,3 +334,161 @@ void Sech2Halo::reset() model_computed = true; dist_defined = true; } + +void Cosine::compute_model(void) +{ + double Emin = get_pot(0.0); + double Emax = get_pot(get_max_radius()); + + Etab.resize(numdf+1); + Ftab.resize(numdf+1); + + const int numz = 10000; + std::vector pot(numz+1), z(numz+1); + for (int i=0; i<=numz; i++) { + z[i] = H * i/numz; + pot[i] = get_pot(z[i]); + } + +#ifdef DEBUG + { + std::ofstream out("cosine_pot.dat"); + if (out) { + out << "#" + << std::setw(15) << std::right << "Z" + << std::setw(16) << std::right << "Potential" + << std::endl; + for (int i=0; i<=numz; i++) { + out << std::setw(16) << std::right << z[i] + << std::setw(16) << std::right << pot[i] + << std::endl; + } + } else { + std::cerr << "Could not open cosine_pot.dat for writing" << std::endl; + } + } +#endif + + Linear1d getZ(pot, z); + + const int numL = 400; + double fac = rho0/(2.0*sqrt(2.0)*H); + LegeQuad lege(numL); + + for (int i=0; i<=numdf; i++) { + Etab[i] = Emin + (Emax - Emin)*i/numdf; + Ftab[i] = 0.0; + double umax = sqrt(Emax - Etab[i]); + for (int j=0; j Z(numZ+1), Rho(numZ+1, 0.0); + const int nint = 800; + LegeQuad lg(nint); + + for (int i=0; i<=numZ; i++) { + Z[i] = get_max_radius() * i/numZ; + double P = get_pot(Z[i]); + double vmax = sqrt(2.0*(Emax - P)); + for (int j=0; j maxdif) { + maxdif = dif; + maxz = Z[i]; + } + } + std::cout << "---- Max density difference from DF integration: " << maxdif + << " at z=" << maxz << std::endl; + + std::ofstream out("cosine_df_check_density.dat"); + if (out) { + out << "#" + << std::setw(15) << std::right << "Z" + << std::setw(16) << std::right << "Density" + << std::setw(16) << std::right << "DF Density" + << std::setw(16) << std::right << "Difference" + << std::endl; + for (int i=0; i<=numZ; i++) { + out << std::setw(16) << std::right << Z[i] + << std::setw(16) << std::right << get_density(Z[i]) + << std::setw(16) << std::right << Rho[i] + << std::setw(16) << std::right << Rho[i] - get_density(Z[i]) + << std::endl; + } + + } else { + std::cerr << "Could not open cosine_df_check.dat for writing" << std::endl; + } + + out.close(); + out.open("cosine_df_check_harmonic.dat"); + if (out) { + out << "#" + << std::setw(15) << std::right << "Z" + << std::setw(16) << std::right << "E" + << std::setw(16) << std::right << "DF" + << std::setw(16) << std::right << "Harmonic DF" + << std::setw(16) << std::right << "Difference" + << std::endl; + + double Omega = sqrt(4.0*M_PI*rho0); + double Sigma2 = 8.0*rho0*H*H/M_PI; + + for (int i=0; i<=numdf; i++) { + double Z = getZ(Etab[i]); + double E = get_pot(Z) + Phi0; + double DF = df_interp(Etab[i]); + double DFh = rho0*Omega/(2.0*M_PI*Sigma2)*exp(-E/Sigma2); + + out << std::setw(16) << std::right << getZ(Z) + << std::setw(16) << std::right << Etab[i] + << std::setw(16) << std::right << DF + << std::setw(16) << std::right << DFh + << std::setw(16) << std::right << DF - DFh + << std::endl; + } + } else { + std::cerr << "Could not open cosine_df_check_harmonic.dat for writing" << std::endl; + } + } +#endif + // END DEBUG +} diff --git a/include/massmodel1d.H b/include/massmodel1d.H new file mode 100644 index 000000000..1d5f42982 --- /dev/null +++ b/include/massmodel1d.H @@ -0,0 +1,1072 @@ +#ifndef _massmodel1d_H +#define _massmodel1d_H + +#include +#include + +#include "massmodel.H" +#include "biorth1d.H" + +class OneDModel : public MassModel +{ + using Complex = std::complex; + using Vector = Eigen::VectorXd; + using CVector = Eigen::VectorXcd; + using Matrix = Eigen::MatrixXd; + using CMatrix = Eigen::MatrixXcd; + +protected: + + int Nfgrid = 800; // Number of points for tabulating + // distribution function for + // interpolation in getE(M) + + Linear1d dfgrid; // Interpolation object for + // distribution function grid (for + // getE(M)) + + std::vector Egrid, Mgrid; + + void compute_dfgrid(void); + +public: + + bool dist_defined; + + //! Null constructor + // AxiSymModel(void) {}; + + //@{ + //! Required members of mass model + + virtual double get_mass (const double) = 0; + virtual double get_density (const double) = 0; + virtual double get_pot (const double) = 0; + virtual double get_dpot (const double) = 0; + virtual double get_dpot2 (const double) = 0; + virtual std::tuple + get_pot_dpot(const double) = 0; + + double get_mass(const double x, const double y, const double z) + { return get_mass(z); } + + double get_density(const double x, const double y, const double z) + { return get_density(z); } + + double get_pot(const double x, const double y, const double z) + { return get_pot(z); } + //@} + + + //@{ + //! Additional member functions + virtual double get_min_radius(void) = 0; + virtual double get_max_radius(void) = 0; + virtual double get_scale_height(void) = 0; + virtual double distf(const double, const double V=0.0) = 0; + virtual double dfde (const double, const double V=0.0) = 0; + virtual double dfdv (const double, const double V=0.0) = 0; + //@} + + + //! Model types + enum class ModelType {table, lowiso, sech2, sech2u, sech2halo, sech2mu, + uniform, cosine}; + + + //! Model names for sysout display (avoid temporaries and dynamic + //! allocations on every class instantiation) + inline static const std::map model_names = { + {ModelType::table, "file"}, + {ModelType::lowiso, "LowIso"}, + {ModelType::sech2, "Sech2"}, + {ModelType::sech2u, "Sech2u"}, + {ModelType::sech2halo, "Sech2Halo"}, + {ModelType::sech2mu, "Sech2(mu)"}, + {ModelType::uniform, "Uniform"}, + {ModelType::cosine, "Cosine"} + }; + + //! Model types for command line parsing + //! 'inline static const' allows initialization inside the class definition. + //! 'std::string_view' keys completely prevent temporary string creations. + inline static const std::map model_map = { + {"table", ModelType::table}, + {"lowiso", ModelType::lowiso}, + {"sech2", ModelType::sech2}, + {"sech2u", ModelType::sech2u}, + {"sech2halo", ModelType::sech2halo}, + {"sech2mu", ModelType::sech2mu}, + {"uniform", ModelType::uniform}, + {"cosine", ModelType::cosine} + }; + + //! Model type for sysout display + static std::string get_model_name(ModelType type) { + return std::string(model_names.at(type)); + } + + //! Command line parsing for model types + static ModelType parse_model_type(const std::string& str) { + std::string MODEL = str; + std::transform(MODEL.begin(), MODEL.end(), MODEL.begin(), + [](unsigned char c){ return std::tolower(c); }); + + // A std::string implicitly and cheaply converts to a std::string_view, + // so finding it in our map works cleanly and creates zero temporaries. + auto it = model_map.find(MODEL); + if (it != model_map.end()) { + return it->second; + } else { + std::ostringstream oss; + oss << "Invalid model type: " << str << ". Valid types are: "; + for (const auto& pair : model_map) { + oss << pair.first << " "; + } + throw std::invalid_argument(oss.str()); + } + } +}; + +class OneDModelTable : public OneDModel +{ +protected: + + Spline1d mass, dens, pot; + + int even; + int num; + int numdf; + double half_height; + std::vector params; + +public: + + OneDModelTable() {}; + + OneDModelTable(std::string filename, int PARM=0); + + OneDModelTable(int num, double *r, double *d, + double *m, double *p, std::string ID = "" ); + + // Required member functions + + double get_mass (const double); + double get_density(const double); + double get_pot (const double); + double get_dpot (const double); + double get_dpot2 (const double); + std::tuple get_pot_dpot (const double); + + // Additional member functions + + const int get_num_param(void) { return params.size(); } + const double get_param(int i) { return params[i-1]; } + double get_scale_height(void) { return half_height; } + + double get_min_radius(void) { return mass.xlo(); } + double get_max_radius(void) { return mass.xhi(); } + int grid_size(void) { return num; } + +// double distf(double E, double V); +// double dfde(double E, double V); +// double dfdv(double E, double V); +}; + +class LowIso : public OneDModelTable +{ +private: + double w0, Bfac, betak, gammak; + double dispx, normx; + + void setup_model(void); + +public: + LowIso(std::string filename, double DISPX=0.159154943091895335768) : + OneDModelTable(filename, 1) { + dispx = DISPX; + setup_model(); + } + + double get_pot(const double); + double get_dpot(const double); + double get_dpot2(const double); + std::tuple get_pot_dpot(const double); + + double distf(const double E, const double V=0.0); + double dfde(const double E, const double V=0.0); + double dfdv(const double E, const double V=0.0); +}; + + +class Sech2 : public OneDModel +{ +private: + double h; + double dispz, dispx; + double norm; + + static double HMAX; + +public: + + Sech2(void) + { + dispz = 1.0; + dispx = 1.0; + // + // Units: G = rho_o = 1 + // + h = sqrt(dispz/(2.0*M_PI)); + norm = 1.0/( sqrt(2.0*M_PI*dispz) ); + dist_defined = true; + } + + Sech2(const double DISPZ, const double DISPX=1.0) + { + dispz = DISPZ; + dispx = DISPX; + // + // Units: G = rho_o = 1 + // + h = sqrt(dispz/(2.0*M_PI)); + norm = 1.0/( sqrt(2.0*M_PI*dispz) ); + dist_defined = true; + } + + double get_mass(const double z) + { + return 2.0*h/(1.0 + exp(-2.0*z/h)); + } + + double get_density(const double z) + { + double zz = fabs(z); + double ret = 2.0*exp(-zz/h)/(1.0 + exp(-2.0*zz/h)); + return ret*ret; + } + + double get_pot(const double z) + { + double zz = fabs(z); + return 4.0*M_PI*h*(zz + h*log(1.0 + exp(-2.0*zz/h)) - h*M_LN2); + } + + double get_dpot(const double z) + { + double zz = fabs(z); + double ret = (1.0 - exp(-2.0*zz/h))/(1.0 + exp(-2.0*zz/h)); + return 4.0*M_PI*h* ret * z/(zz+1.0e-18); + } + + double get_dpot2(const double z) + { + double zz = fabs(z); + double ret = 2.0*exp(-zz/h)/(1.0 + exp(-2*zz/h)); + return 4.0*M_PI*ret*ret; + } + + std::tuple get_pot_dpot(const double z) + { + double zz = fabs(z); + double p = 4.0*M_PI*h*(zz + h*log(1.0 + exp(-2.0*zz/h)) - h*M_LN2); + double ret = (1.0 - exp(-2.0*zz/h))/(1.0 + exp(-2.0*zz/h)); + double dp = 4.0*M_PI*h* ret * z/(zz+1.0e-18); + return {p, dp}; + } + + double get_mass(const double x, const double y, const double z) + { + return get_mass(z); + } + + double get_density(const double x, const double y, const double z) + { + return get_density(z); + } + + double get_pot(const double x, const double y, const double z) + { + return get_pot(z); + } + + double get_min_radius(void) { return 0.0; } + double get_max_radius(void) { return HMAX*h; } + double get_scale_height(void) { return h; } + + static void set_hmax(double hmax) { HMAX = hmax; } + + double distf(const double E, const double p) + { + return exp(-E/dispz - 0.5*p*p/dispx) * norm; + } + + double dfde(const double E, const double p=0.0) + { + return -exp(-E/dispz - 0.5*p*p/dispx)/dispz * norm; + } + + double dfdv(const double E, const double p=0.0) + { + return -exp(-E/dispz - 0.5*p*p/dispx)*p/dispx * norm; + } + +}; + + +class Sech2u : public OneDModel +{ +private: + double h, rho0, S0; + double dispz, dispx; + double norm; + + static double HMAX; + +public: + + Sech2u(void) + { + S0 = 1.0; + dispz = 1.0; + dispx = 1.0; + + // Units: G = Sigma_o = 1 + // + rho0 = S0/(4.0*h); + h = dispz/(2.0*M_PI*S0); + norm = S0/( sqrt(2.0*M_PI*dispz) ) / (4.0*h); + dist_defined = true; + } + + Sech2u(const double DISPZ, const double DISPX) + { + S0 = 1.0; + dispz = DISPZ; + dispx = DISPX; + + // Units: G = Sigma_o = 1 + // + h = dispz/(2.0*M_PI*S0); + norm = S0/( sqrt(2.0*M_PI*dispz) ) / (4.0*h); + dist_defined = true; + } + + Sech2u(const double S, const double DISPZ, const double DISPX) + { + S0 = S; + dispz = DISPZ; + dispx = DISPX; + + // Units: G = Sigma_o = 1 + // + h = dispz/(2.0*M_PI*S0); + norm = S0/( sqrt(2.0*M_PI*dispz) ) / (4.0*h); + dist_defined = true; + } + + double get_mass(const double z) + { + // To avoid overflow, split into two cases. + if (z<0.0) + return S0*exp(z/h)/(1.0 + exp(z/h)); + else + return S0/(1.0 + exp(-z/h)); + } + + double get_density(const double z) + { + double zz = fabs(z); // Even function + double ret = 1.0 + exp(-zz/h); + return S0*exp(-zz/h)/(ret*ret)/h; + } + + double get_pot(const double z) + { + double zz = fabs(z); // Even function + return 2.0*M_PI*S0*(zz + 2.0*h*log(1.0 + exp(-zz/h)) - 2.0*h*M_LN2); + } + + double get_dpot(const double z) + { + double zz = fabs(z); // Odd function + double ret = (1.0 - exp(-zz/h))/(1.0 + exp(-zz/h)); + return 2.0*M_PI*S0 * ret * std::copysign(1.0, z); + } + + double get_dpot2(const double z) + { + double zz = fabs(z); // Even function + double ret = exp(-zz/(2.0*h))/(1.0 + exp(-zz/h)); + return 4.0*M_PI*S0 * ret*ret/h; + } + + std::tuple get_pot_dpot(const double z) + { + double zz = fabs(z); + double p = 2.0*M_PI*S0 * (zz + 2.0*h*log(1.0 + exp(-zz/h)) - 2.0*h*M_LN2); + double ret = (1.0 - exp(-zz/h))/(1.0 + exp(-zz/h)); + double dp = 2.0*M_PI*S0 * ret * std::copysign(1.0, z); + return {p, dp}; + } + + double get_mass(const double x, const double y, const double z) + { + return get_mass(z); + } + + double get_density(const double x, const double y, const double z) + { + return get_density(z); + } + + double get_pot(const double x, const double y, const double z) + { + return get_pot(z); + } + + double get_min_radius(void) { return 0.0; } + double get_max_radius(void) { return HMAX*h; } + double get_scale_height(void) { return h; } + + void set_vdisp(double DISPZ, double DISPX=1.0) + { + dispz = DISPZ; + dispx = DISPX; + h = dispz/(2.0*M_PI*S0); + norm = S0/( sqrt(2.0*M_PI*dispz) ) / (4.0*h); + dist_defined = true; + } + + static void set_hmax(double hmax) { HMAX = hmax; } + + double distf(const double E, const double p=0.0) + { + return exp(-E/dispz - 0.5*p*p/dispx) * norm; + } + + double dfde(const double E, const double p=0.0) + { + return -exp(-E/dispz - 0.5*p*p/dispx)/dispz * norm; + } + + double dfdv(const double E, const double p=0.0) + { + return -exp(-E/dispz - 0.5*p*p/dispx)*p/dispx * norm; + } + +}; + + + +class Sech2mu : public OneDModel +{ +private: + double mu, h; + double dispz, dispx; + double dnorm, fnorm; + + static double HMAX; + +public: + + Sech2mu(void) + { + // Units: G = 1 + // + mu = 1.0; + dispz = 1.0; + dispx = 1.0; + h = dispz/(2.0*M_PI*mu); + + // Normalizations + // + dnorm = 0.25*mu/h; + fnorm = dnorm/sqrt(2.0*M_PI*dispz); + dist_defined = true; + } + + Sech2mu(const double DISPZ, const double DISPX=1.0) + { + mu = 1.0; + dispz = DISPZ; + dispx = DISPX; + h = dispz/(2.0*M_PI*mu); + // + // Units: G = 1 + // + dnorm = 0.25*mu/h; + fnorm = dnorm/sqrt(2.0*M_PI*dispz); + dist_defined = true; + } + + void setMu(double Sigma0) + { + mu = Sigma0; + h = dispz/(2.0*M_PI*mu); + // + // Units: G = 1, mu = Sigma0 + // + dnorm = 0.25*mu/h; + fnorm = dnorm/sqrt(2.0*M_PI*dispz); + dist_defined = true; + } + + double get_mass(const double z) { + return mu*exp(z/h)/(1.0 + exp(z/h)); + } + + double get_density(const double z) { + double zz = fabs(z); + double fac = exp(zz/h); + return 4.0*dnorm/(fac + 1.0/fac + 2.0); + } + + double get_pot(const double z) { + double zz = fabs(z); + // Define potential to be zero at z=0, so subtract off the constant term. + return 2.0*dispz*(0.5*zz/h + log(1.0 + exp(-zz/h)) - M_LN2); + } + + double get_dpot(const double z) { + double zz = fabs(z); + double ret = (1.0 - exp(-zz/h))/(1.0 + exp(-zz/h)); + return dispz*ret/h * z/(zz+1.0e-18); + } + + double get_dpot2(const double z) { + double zz = fabs(z); + double ret = 2.0*exp(-0.5*zz/h)/(1.0 + exp(-zz/h)); + return 0.5*dispz*ret*ret/h/h; + } + + std::tuple get_pot_dpot(const double z) + { + double zz = fabs(z); + double p = 2.0*dispz*(0.5*zz/h + log(1.0 + exp(-zz/h)) - M_LN2); + double ret = (1.0 - exp(-zz/h))/(1.0 + exp(-zz/h)); + double dp = dispz*ret/h * z/(zz+1.0e-18); + return {p, dp}; + } + + double get_mass(const double x, const double y, const double z) { + return get_mass(z); + } + + double get_density(const double x, const double y, const double z) { + return get_density(z); + } + + double get_pot(const double x, const double y, const double z) { + return get_pot(z); + } + + double get_min_radius(void) { return 0.0; } + double get_max_radius(void) { return HMAX*h; } + double get_scale_height(void) { return h; } + + static void set_hmax(double hmax) { HMAX = hmax; } + + double distf(const double E, const double p) { + return exp(-E/dispz - 0.5*p*p/dispx) * fnorm; + } + + double dfde(const double E, const double p=0.0) { + return -exp(-E/dispz - 0.5*p*p/dispx)/dispz * fnorm; + } + + double dfdv(const double E, const double p=0.0) { + return -exp(-E/dispz - 0.5*p*p/dispx)*p/dispx * fnorm; + } + +}; + + +class Sech2Halo : public OneDModelTable +{ +private: + double h, rho0; + double dispz, dispx; + double dratio, hratio; + double hh, rho0h; + double norm, hmax; + + bool model_computed; + + static double HMAX; + + void reset(); + +public: + + static int NTABLE; + static double OFFSET; + static bool MU; + + Sech2Halo(void) { + dispz = 1.0; + dispx = 1.0; + dratio = 0.0; + hratio = 1.0; + + dist_defined = false; + model_computed = false; + } + + Sech2Halo(const double DISPZ, const double DRATIO, const double HRATIO, + const double DISPX=1.0); + + double get_pot(const double z) { + if (!model_computed) reset(); + return OneDModelTable::get_pot(z) + 4.0*M_PI*hh*hh*rho0h*log(cosh(z/hh)); + } + + double get_dpot(const double z) { + if (!model_computed) reset(); + return OneDModelTable::get_dpot(z) + 4.0*M_PI*hh*rho0h*tanh(z/hh); + } + + double get_dpot2(const double z) { + if (!model_computed) reset(); + double sech = 1.0/cosh(z/hh); + return OneDModelTable::get_dpot2(z) + 4.0*M_PI*rho0h*sech*sech; + } + + std::tuple get_pot_dpot(const double z) { + if (!model_computed) reset(); + auto [p, dp] = OneDModelTable::get_pot_dpot(z); + p += 4.0*M_PI*hh*hh*rho0h*log(cosh(z/hh)); + dp += 4.0*M_PI*hh*rho0h*tanh(z/hh); + return {p, dp}; + } + + + double get_min_radius(void) { return 0.0; } + double get_max_radius(void) { return hmax*h; } + double get_scale_height(void) { return h; } + double get_scale_height_halo(void) { return hh; } + double get_rho0(void) { return rho0; } + double get_rho0_halo(void) { return rho0h; } + + static void set_hmax(double hmax) { HMAX = hmax; } + + double distf(const double E, const double p=0.0) { + if (!model_computed) reset(); + return exp(-E/dispz - 0.5*p*p/dispx) * norm; + } + + double dfde(const double E, const double p=0.0) { + if (!model_computed) reset(); + return -exp(-E/dispz - 0.5*p*p/dispx)/dispz * norm; + } + + double dfdv(const double E, const double p=0.0) { + if (!model_computed) reset(); + return -exp(-E/dispz - 0.5*p*p/dispx)*p/dispx * norm; + } + +}; + +//! A simple uniform slab model with finite thickness and constant density +class Uniform : public OneDModel +{ +private: + //! Half-thickness of the uniform slab + double H; + + //! Velocity widths in x and y directions (for distribution function) + double Vx, Vy; + + //! Vertical frequency for the uniform slab potential + double Omega; + + //! Normalization constant for the distribution function + double norm; + + //! Maximum action + double Jmax; + + //! Regularization parameter for the distribution function (if needed) + double epsilon; + +public: + + //! Null constructor with default parameters + Uniform(void) + { + Vx = 1.0; + Vy = 1.0; + // + // Units: G = Sigma_o = 1 + // + H = 1.0; + Omega = std::sqrt(2.0*M_PI/H); + norm = std::sqrt(2/Omega)/(M_PI*2.0*H); + Jmax = 0.5*Omega*H*H; + ModelID = "Uniform"; + } + + //! Constructor with specified half-thickness and velocity width + //! (same for x and y) + Uniform(double h, double Vx, double epsilon) + { + this->H = h; + this->Vx = Vx; + this->Vy = Vx; + this->epsilon = epsilon; + // + // Units: G = Sigma_o = 1 + // + Omega = std::sqrt(2.0*M_PI/H); + norm = std::sqrt(2/Omega)/(M_PI*2.0*H); + Jmax = 0.5*Omega*H*H; + ModelID = "Uniform"; + dist_defined = true; + } + + //! Constructor with specified half-thickness and velocity widths in + //! x and y + Uniform(double h, double Vx, double Vy, double epsilon) + { + this->H = h; + this->Vx = Vx; + this->Vy = Vy; + this->epsilon = epsilon; + // + // Units: G = Sigma_o = 1 + // + Omega = std::sqrt(2.0*M_PI/H); + norm = std::sqrt(2/Omega)/(M_PI*2.0*H); + Jmax = 0.5*Omega*H*H; + ModelID = "Uniform"; + dist_defined = true; + } + + double get_mass(const double z) + { + if (z<-H) + return 0.0; + else if (z > H) + return 1.0; + else { + return (z + H)/(2.0*H); + } + } + + double get_density(const double z) + { + if (fabs(z) > H) + return 0.0; + else + return 1.0/(2.0*H); + } + + double get_pot(const double z) + { + // Offset defined so that potential is zero at |z|=H, and negative + // for |z| H) + return 2.0*M_PI*zz - 2.0*M_PI*H; + else + return M_PI*zz*zz/H - M_PI*H; + } + + double get_dpot(const double z) + { + double zz = fabs(z); + if (zz > H) + return 2.0*M_PI*z/zz; + else + return 2.0*M_PI*z/H; + } + + double get_dpot2(const double z) + { + double zz = fabs(z); + if (zz > H) + return 0.0; + else + return 2.0*M_PI/H; + } + + std::tuple get_pot_dpot(const double z) + { + double zz = fabs(z), p, dp; + if (zz > H) { + p = 2.0*M_PI*zz - 2.0*M_PI*H; + dp = 2.0*M_PI*z/zz; + } else { + p = M_PI*zz*zz/H - M_PI*H; + dp = 2.0*M_PI*z/H; + } + return {p, dp}; + } + + double get_mass(const double x, const double y, const double z) + { + return get_mass(z); + } + + double get_density(const double x, const double y, const double z) + { + return get_density(z); + } + + double get_pot(const double x, const double y, const double z) + { + return get_pot(z); + } + + double get_min_radius(void) { return 0.0; } + double get_max_radius(void) { return H; } + double get_scale_height(void) { return H; } + + double distf(const double E, const double p=0.0) + { + // Convert energy to potential, so that the distribution function + // is nonzero only. That is, assume E = Omega*Ja, not zero at +/- + // H. + double J = fabs(E + M_PI*H)/Omega; + return norm / std::sqrt(Jmax - J + epsilon); + } + + double dfde(const double E, const double p=0.0) + { + // Convert energy to potential, so that the distribution function + // is nonzero only. That is, assume E = Omega*Ja, not zero at +/- + // H. + double J = fabs(E + M_PI*H)/Omega; + return -0.5*Omega*norm / std::pow(Jmax - J + epsilon, 1.5); + } + + double dfdv(const double E, const double p=0.0) + { + // The distribution function is independent of velocity, so this + // derivative is zero. There are delta functions at v=-Vx and + // v=+Vx, but those are not captured by this simple functional + // form. + return 0.0; + } +}; + +//! A simple consine-bell slab model with finite thickness +class Cosine : public OneDModel +{ +private: + //! Half-thickness of the uniform slab + double H; + + //! Central density of the cosine slab + double rho0; + + //! Potential zero point is defined at |z|=H, so that potential is + //! negative inside the slab and zero at the edges + double Phi0; + + //! Velocity widths in x and y directions (for distribution function) + double Vx, Vy; + + //! Vertical frequency for the uniform slab potential + double Omega; + + //! Normalization constant for the distribution function + double norm; + + //! Regularization parameter for the distribution function (if needed) + double epsilon; + + //! Number of points in the tabulated distribution function + int numdf; + + //! Number of quadrature points for computing DF + int nint; + + //! Tabulated action and distribution function for interpolation + std::vector Etab, Ftab; + + //! DF table interpolation object + Linear1d df_interp; + + //! Precompute the distribution function table for interpolation + void compute_model(); + +public: + + //! Null constructor with default parameters + Cosine(void) + { + Vx = 1.0; + Vy = 1.0; + // + // Units: G = Sigma_o = 1 + // + H = 1.0; + rho0 = 1.0/H; + Phi0 = M_PI*H*(4.0/(M_PI*M_PI) + 1.0); + ModelID = "Cosine"; + numdf = 800; + nint = 800; + dist_defined = true; + compute_model(); + } + + //! Constructor with specified half-thickness and velocity width + //! (same for x and y) + Cosine(double h, double Vx, double epsilon) + { + this->H = h; + this->Vx = Vx; + this->Vy = Vx; + this->epsilon = epsilon; + // + // Units: G = Sigma_o = 1 + // + rho0 = 1.0/H; + Phi0 = M_PI*H*(4.0/(M_PI*M_PI) + 1.0); + ModelID = "Cosine"; + numdf = 800; + nint = 800; + dist_defined = true; + compute_model(); + } + + //! Constructor with specified half-thickness and velocity widths in + //! x and y + Cosine(double h, double Vx, double Vy, double epsilon) + { + this->H = h; + this->Vx = Vx; + this->Vy = Vy; + this->epsilon = epsilon; + // + // Units: G = Sigma_o = 1 + // + rho0 = 1.0/H; + Phi0 = M_PI*H*(4.0/(M_PI*M_PI) + 1.0); + ModelID = "Cosine"; + numdf = 800; + nint = 800; + dist_defined = true; + compute_model(); + } + + void setDF(int numdf, int nint=800) + { + this->numdf = numdf; + this->nint = nint; + compute_model(); + } + + double get_mass(const double z) + { + if (z<-H) + return 0.0; + else if (z > H) + return 1.0; + else { + return 0.5*rho0*(z + H + H/M_PI*sin(M_PI*z/H)); + } + } + + double get_density(const double z) + { + if (fabs(z) > H) + return 0.0; + else { + double cosfac = cos(M_PI*z/(2.0*H)); + return rho0*cosfac*cosfac; + } + } + + double get_pot(const double z) + { + // Offset defined so that potential is zero at |z|=H, and negative + // for |z| H) + return 2.0*M_PI*rho0*H*(zz - H); + else + return 2.0*M_PI*rho0*(0.5*z*z + H*H/(M_PI*M_PI)*(1 - cos(M_PI*z/H))) - Phi0; + } + + double get_dpot(const double z) + { + double zz = fabs(z); + if (zz > H) + return 2.0*M_PI*rho0*H*z/zz; + else + return 2.0*M_PI*rho0*(z + H/M_PI*sin(M_PI*z/H)); + } + + double get_dpot2(const double z) + { + double zz = fabs(z); + if (zz > H) + return 0.0; + else + return 2.0*M_PI*rho0*(1 + cos(M_PI*z/H)); + } + + std::tuple get_pot_dpot(const double z) + { + double zz = fabs(z), p, dp; + if (zz > H) { + p = 2.0*M_PI*rho0*H*(zz - H); + dp = 2.0*M_PI*rho0*H*z/zz; + } else { + p = 2.0*M_PI*rho0*(0.5*z*z + H*H/(M_PI*M_PI)*(1 - cos(M_PI*z/H))) - Phi0; + dp = 2.0*M_PI*rho0*(z + H/M_PI*sin(M_PI*z/H)); + } + return {p, dp}; + } + + double get_mass(const double x, const double y, const double z) + { + return get_mass(z); + } + + double get_density(const double x, const double y, const double z) + { + return get_density(z); + } + + double get_pot(const double x, const double y, const double z) + { + return get_pot(z); + } + + double get_min_radius(void) { return 0.0; } + double get_max_radius(void) { return H; } + double get_scale_height(void) { return H; } + + double distf(const double E, const double p=0.0) + { + if (E>=0.0) + return 0.0; + else + return df_interp.eval(E); + } + + double dfde(const double E, const double p=0.0) + { + if (E>=0.0) + return 0.0; + else + return df_interp.deriv(E); + } + + double dfdv(const double E, const double p=0.0) + { + // The distribution function is independent of velocity, so this + // derivative is zero. There are delta functions at v=-Vx and + // v=+Vx, but those are not captured by this simple functional + // form. + return 0.0; + } + +}; + +#endif + +// -*- C++ -*- diff --git a/utils/ICs/massmodel1d.H b/utils/ICs/massmodel1d.H deleted file mode 100755 index 06033f522..000000000 --- a/utils/ICs/massmodel1d.H +++ /dev/null @@ -1,442 +0,0 @@ -// -*- C++ -*- - -#ifndef _massmodel1d_H -#define _massmodel1d_H - -#include -#include - -#include "massmodel.H" - -class OneDModel : public MassModel -{ - -public: - bool dist_defined; - -// AxiSymModel(void) {}; - - virtual double get_mass(const double) = 0; - virtual double get_density(const double) = 0; - virtual double get_pot(const double) = 0; - virtual double get_dpot(const double) = 0; - virtual double get_dpot2 (const double) = 0; - virtual tuple get_pot_dpot (const double) = 0; - - // Required members of mass model - - double get_mass(const double x, const double y, const double z) { - return get_mass(z);} - - double get_density(const double x, const double y, const double z) { - return get_density(z);} - - double get_pot(const double x, const double y, const double z) { - return get_pot(z);} - - - // Addiional member functions - - virtual double get_min_radius(void) = 0; - virtual double get_max_radius(void) = 0; - virtual double get_scale_height(void) = 0; - virtual double distf(const double, const double V=0.0) = 0; - virtual double dfde (const double, const double V=0.0) = 0; - virtual double dfdv (const double, const double V=0.0) = 0; -}; - -class OneDModelTable : public OneDModel -{ -protected: - - Spline1d mass, dens, pot; - - int even; - int num; - int numdf; - double half_height; - std::vector params; - -public: - - OneDModelTable() {}; - - OneDModelTable(string filename, int PARM=0); - - OneDModelTable(int num, double *r, double *d, - double *m, double *p, string ID = "" ); - - // Required member functions - - double get_mass (const double); - double get_density(const double); - double get_pot (const double); - double get_dpot (const double); - double get_dpot2 (const double); - tuple get_pot_dpot (const double); - - // Additional member functions - - const int get_num_param(void) { return params.size(); } - const double get_param(int i) { return params[i-1]; } - double get_scale_height(void) { return half_height; } - - double get_min_radius(void) { return mass.xlo(); } - double get_max_radius(void) { return mass.xhi(); } - int grid_size(void) { return num; } - -// double distf(double E, double V); -// double dfde(double E, double V); -// double dfdv(double E, double V); -}; - -class LowIso : public OneDModelTable -{ -private: - double w0, Bfac, betak, gammak; - double dispx, normx; - - void setup_model(void); - -public: - LowIso(string filename, double DISPX=0.159154943091895335768) : - OneDModelTable(filename, 1) { - dispx = DISPX; - setup_model(); - } - - double get_pot(const double); - double get_dpot(const double); - double get_dpot2(const double); - tuple get_pot_dpot(const double); - - double distf(const double E, const double V=0.0); - double dfde(const double E, const double V=0.0); - double dfdv(const double E, const double V=0.0); -}; - - -class Sech2 : public OneDModel -{ -private: - double h; - double dispz, dispx; - double norm; - - static double HMAX; - -public: - - Sech2(void) - { - dispz = 1.0; - dispx = 1.0; - // - // Units: G = rho_o = 1 - // - h = sqrt(dispz/(2.0*M_PI)); - norm = 1.0/( sqrt(2.0*M_PI*dispz) ); - dist_defined = true; - } - - Sech2(const double DISPZ, const double DISPX=1.0) - { - dispz = DISPZ; - dispx = DISPX; - // - // Units: G = rho_o = 1 - // - h = sqrt(dispz/(2.0*M_PI)); - norm = 1.0/( sqrt(2.0*M_PI*dispz) ); - dist_defined = true; - } - - double get_mass(const double z) - { - return 2.0*h/(1.0 + exp(-2.0*z/h)); - } - - double get_density(const double z) - { - double zz = fabs(z); - double ret = 2.0*exp(-zz/h)/(1.0 + exp(-2.0*zz/h)); - return ret*ret; - } - - double get_pot(const double z) - { - double zz = fabs(z); - return 4.0*M_PI*h*(zz + h*log(1.0 + exp(-2.0*zz/h)) - h*M_LN2); - } - - double get_dpot(const double z) - { - double zz = fabs(z); - double ret = (1.0 - exp(-2.0*zz/h))/(1.0 + exp(-2.0*zz/h)); - return 4.0*M_PI*h* ret * z/(zz+1.0e-18); - } - - double get_dpot2(const double z) - { - double zz = fabs(z); - double ret = 2.0*exp(-zz/h)/(1.0 + exp(-2*zz/h)); - return 4.0*M_PI*ret*ret; - } - - tuple get_pot_dpot(const double z) - { - double zz = fabs(z); - double p = 4.0*M_PI*h*(zz + h*log(1.0 + exp(-2.0*zz/h)) - h*M_LN2); - double ret = (1.0 - exp(-2.0*zz/h))/(1.0 + exp(-2.0*zz/h)); - double dp = 4.0*M_PI*h* ret * z/(zz+1.0e-18); - return {p, dp}; - } - - double get_mass(const double x, const double y, const double z) - { - return get_mass(z); - } - - double get_density(const double x, const double y, const double z) - { - return get_density(z); - } - - double get_pot(const double x, const double y, const double z) - { - return get_pot(z); - } - - double get_min_radius(void) { return 0.0; } - double get_max_radius(void) { return HMAX*h; } - double get_scale_height(void) { return h; } - - static void set_hmax(double hmax) { HMAX = hmax; } - - double distf(const double E, const double p) - { - return exp(-E/dispz - 0.5*p*p/dispx) * norm; - } - - double dfde(const double E, const double p) - { - return -exp(-E/dispz - 0.5*p*p/dispx)/dispz * norm; - } - - double dfdv(const double E, const double p) - { - return -exp(-E/dispz - 0.5*p*p/dispx)*p/dispx * norm; - } - -}; - - - -class Sech2mu : public OneDModel -{ -private: - double mu, h; - double dispz, dispx; - double dnorm, knorm, fnorm; - - static double HMAX; - -public: - - Sech2mu(void) - { - mu = 1.0; - dispz = 1.0; - dispx = 1.0; - h = 1.0; - // - // Units: G = mu = 1 - // - dnorm = 0.25*mu/h; - knorm = 2.0*M_PI*mu*h/dispz; - fnorm = mu/(4.0*h*knorm*sqrt(2.0*M_PI*dispz)); - dist_defined = true; - } - - Sech2mu(const double DISPZ, const double H, const double DISPX=1.0) - { - mu = 1.0; - dispz = DISPZ; - dispx = DISPX; - h = H; - // - // Units: G = mu = 1 - // - dnorm = 0.25*mu/h; - knorm = 2.0*M_PI*mu*h/dispz; - fnorm = mu/(4.0*h*knorm*sqrt(2.0*M_PI*dispz)); - dist_defined = true; - } - - void setMu(double Sigma0) - { - mu = Sigma0; - // - // Units: G = 1, mu = Sigma0 - // - dnorm = 0.25*mu/h; - knorm = 2.0*M_PI*mu*h/dispz; - fnorm = mu/(4.0*h*knorm*sqrt(2.0*M_PI*dispz)); - dist_defined = true; - } - - double get_mass(const double z) { - return mu*exp(z/h)/(1.0 + exp(z/h)); - } - - double get_density(const double z) { - double zz = fabs(z); - double fac = exp(zz/h); - return 4.0*dnorm/(fac + 1.0/fac + 2.0); - } - - double get_pot(const double z) { - double zz = fabs(z); - return 2.0*dispz*(-M_LN2 + 0.5*zz/h + log(1.0 + exp(-zz/h)) ) - dispz*log(knorm); - } - - double get_dpot(const double z) { - double zz = fabs(z); - double ret = (1.0 - exp(-zz/h))/(1.0 + exp(-zz/h)); - return 0.5*dispz*ret/h * z/(zz+1.0e-18); - } - - double get_dpot2(const double z) { - double zz = fabs(z); - double ret = 2.0*exp(-0.5*zz/h)/(1.0 + exp(-zz/h)); - return 0.5*dispz*ret*ret/h/h; - } - - void get_pot_dpot(const double z, double& p, double& dp) { - double zz = fabs(z); - p = 2.0*dispz*(-M_LN2 + 0.5*zz/h + log(1.0 + exp(-zz/h)) ) - dispz*log(knorm); - double ret = (1.0 - exp(-zz/h))/(1.0 + exp(-zz/h)); - dp = 0.5*dispz*ret/h * z/(zz+1.0e-18); - } - - double get_mass(const double x, const double y, const double z) { - return get_mass(z); - } - - double get_density(const double x, const double y, const double z) { - return get_density(z); - } - - double get_pot(const double x, const double y, const double z) { - return get_pot(z); - } - - double get_min_radius(void) { return 0.0; } - double get_max_radius(void) { return HMAX*h; } - double get_scale_height(void) { return h; } - - static void set_hmax(double hmax) { HMAX = hmax; } - - double distf(const double E, const double p) { - return exp(-E/dispz - 0.5*p*p/dispx) * fnorm; - } - - double dfde(const double E, const double p) { - return -exp(-E/dispz - 0.5*p*p/dispx)/dispz * fnorm; - } - - double dfdv(const double E, const double p) { - return -exp(-E/dispz - 0.5*p*p/dispx)*p/dispx * fnorm; - } - -}; - - -class Sech2Halo : public OneDModelTable -{ -private: - double h, rho0; - double dispz, dispx; - double dratio, hratio; - double hh, rho0h; - double norm, hmax; - - bool model_computed; - - static double HMAX; - - void reset(); - -public: - - static int NTABLE; - static double OFFSET; - static bool MU; - - Sech2Halo(void) { - dispz = 1.0; - dispx = 1.0; - dratio = 0.0; - hratio = 1.0; - - dist_defined = false; - model_computed = false; - } - - Sech2Halo(const double DISPZ, const double DRATIO, const double HRATIO, - const double DISPX=1.0); - - double get_pot(const double z) { - if (!model_computed) reset(); - return OneDModelTable::get_pot(z) + 4.0*M_PI*hh*hh*rho0h*log(cosh(z/hh)); - } - - double get_dpot(const double z) { - if (!model_computed) reset(); - return OneDModelTable::get_dpot(z) + 4.0*M_PI*hh*rho0h*tanh(z/hh); - } - - double get_dpot2(const double z) { - if (!model_computed) reset(); - double sech = 1.0/cosh(z/hh); - return OneDModelTable::get_dpot2(z) + 4.0*M_PI*rho0h*sech*sech; - } - - std::tuple get_pot_dpot(const double z) { - if (!model_computed) reset(); - auto [p, dp] = OneDModelTable::get_pot_dpot(z); - p += 4.0*M_PI*hh*hh*rho0h*log(cosh(z/hh)); - dp += 4.0*M_PI*hh*rho0h*tanh(z/hh); - return {p, dp}; - } - - - double get_min_radius(void) { return 0.0; } - double get_max_radius(void) { return hmax*h; } - double get_scale_height(void) { return h; } - double get_scale_height_halo(void) { return hh; } - double get_rho0(void) { return rho0; } - double get_rho0_halo(void) { return rho0h; } - - static void set_hmax(double hmax) { HMAX = hmax; } - - double distf(const double E, const double p=0.0) { - if (!model_computed) reset(); - return exp(-E/dispz - 0.5*p*p/dispx) * norm; - } - - double dfde(const double E, const double p=0.0) { - if (!model_computed) reset(); - return -exp(-E/dispz - 0.5*p*p/dispx)/dispz * norm; - } - - double dfdv(const double E, const double p=0.0) { - if (!model_computed) reset(); - return -exp(-E/dispz - 0.5*p*p/dispx)*p/dispx * norm; - } - -}; - - -#endif - From 5852d9acfce5d9e50e0a453c1a6542c82471f25d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 14 Aug 2026 13:48:15 -0400 Subject: [PATCH 04/28] Remove a stray character --- utils/ICs/genslab.cc | 247 ++++++++++++++++++++++++++++++------------- 1 file changed, 172 insertions(+), 75 deletions(-) mode change 100755 => 100644 utils/ICs/genslab.cc diff --git a/utils/ICs/genslab.cc b/utils/ICs/genslab.cc old mode 100755 new mode 100644 index e8798197c..fc2282a39 --- a/utils/ICs/genslab.cc +++ b/utils/ICs/genslab.cc @@ -1,37 +1,13 @@ -/***************************************************************************** - * Description: - * ----------- - * - * Generate slab initial conditions in a unit sqaure - * - * - * Call sequence: - * ------------- - * - * Parameters: - * ---------- - * - * - * Returns: - * ------- - * - * - * Notes: - * ----- - * - * - * By: - * -- - * - * MDW 11/20/91 - * - ***************************************************************************/ +// Generate slab initial conditions for varous slab models. The +// output can be in ASCII or HDF5 format, with optional compression +// filters for HDF5. #include #include #include #include #include +#include #include #include #include @@ -40,7 +16,6 @@ #include "massmodel1d.H" #include "interp.H" - #include "cxxopts.H" #include "ParticleHDF5.H" @@ -60,6 +35,7 @@ struct ParticleData size_t num_aux_floats = 0; }; + // Write particle data to HDF5 file with specified filter and precision template void write_hdf5_data(const std::string& outfile, @@ -84,24 +60,36 @@ int main(int argc, char **argv) { unsigned int seed; - int Ntable, Number; - double Dratio, Hratio, R, Hmax, DispX, DispZ, fJ; + int Ntable, Num_particles, Num_aux_ints, Num_aux_floats; + double Dratio, Hratio, R, Hmax, DispX, DispZ, fJ, Lx, Ly; std::string outfile, config, modfile, modelType; - bool Mu; + unsigned filter_id = 1; + bool Mu, HDF5 = true; // Parse command line // - std::string message = "Generate unit-box slab initial conditions\n"; + std::string message = "Generate slab initial conditions for varous slab models. The\n" + "output can be in ASCII or HDF5 format, with optional compression\n" + "filters for HDF5.\n"; cxxopts::Options options(argv[0], message); options.add_options() ("h,help", "Print this help message") + ("5,hdf5", "Write HDF5 output (default)") + ("A,ascii", "Write old-style ASCII output (default is HDF5)") + ("v,verbose", "Verbose output") + ("8,double", "Use double precision for HDF5 output (default is float)") + ("f,filter", "HDF5 filter ID to use (default: 1 = GZIP)", cxxopts::value(filter_id)->default_value("1")) ("N,number", "Number of bodies", - cxxopts::value(Number)->default_value("10000")) + cxxopts::value(Num_particles)->default_value("10000")) + ("nauxint", "Number of auxiliary integer fields", + cxxopts::value(Num_aux_ints)->default_value("0")) + ("nauxfloat", "Number of auxiliary float fields", + cxxopts::value(Num_aux_floats)->default_value("0")) ("n,ntable", "Number of points in model table", cxxopts::value(Ntable)->default_value("400")) - ("t,model", "Model type (LowIso, Sech2, Sech2Halo)", + ("m,model", "Model type (LowIso, Sech2, Sech2mu, Sech2Halo)", cxxopts::value(modelType)->default_value("Sech2")) ("d,dratio", "Ratio of disk to halo density", cxxopts::value(Dratio)->default_value("3.0")) @@ -111,6 +99,10 @@ main(int argc, char **argv) cxxopts::value(R)->default_value("1.0")) ("H,hmax", "Maximum vertical size in scale heights", cxxopts::value(Hmax)->default_value("10.0")) + ("x,Lx", "Slab length in the x-direction", + cxxopts::value(Lx)->default_value("1.0")) + ("y,Ly", "Slab length in the y-direction", + cxxopts::value(Ly)->default_value("1.0")) ("X,DispX", "In-plane velocity variance", cxxopts::value(DispX)->default_value("1.0")) ("F,fJ", "Ratio of Jeans length to box scale", @@ -129,6 +121,8 @@ main(int argc, char **argv) cxxopts::ParseResult vm; + // Parse the command line + // try { vm = options.parse(argc, argv); } catch (cxxopts::OptionException& e) { @@ -143,15 +137,11 @@ main(int argc, char **argv) return 0; } - std::ofstream out(outfile); - if (!out) { - std::cerr << "Can't open <" << outfile << ">" << std::endl; - exit(-1); - } - out.precision(6); - out.setf(ios::scientific); - + // Select output format // + if (vm.count("hdf5")) HDF5 = true; + if (vm.count("ascii")) HDF5 = false; + // Define model // double h = 1.0; @@ -161,9 +151,13 @@ main(int argc, char **argv) if (modelType.compare("LowIso") == 0) { model = std::make_shared(modfile); } + else if (modelType.compare("Sech2mu") == 0) { + auto tmp = std::make_shared(DispZ, DispX); + h = tmp->get_scale_height(); + if (Hmax>0) tmp->set_hmax(Hmax); + model = tmp; + } else if (modelType.compare("Sech2") == 0) { - DispZ = 2.0/(M_PI*fJ*fJ); - DispX = DispZ/(R*R); auto tmp = std::make_shared(DispZ); h = tmp->get_scale_height(); if (Hmax>0) tmp->set_hmax(Hmax); @@ -182,10 +176,8 @@ main(int argc, char **argv) exit(-1); } - // // Jeans' length for selecting scale height // - double maxZ = model->get_max_radius(); double mu = model->get_mass(maxZ); double KJ = 2.0*M_PI*mu/DispX; @@ -194,18 +186,19 @@ main(int argc, char **argv) std::cout.setf(ios::left); char prev = cout.fill('.'); - std::cout << std::setw(40) << "Model type" << modelType << std::endl; - std::cout << std::setw(40) << "Surface mass density" << mu << std::endl; - std::cout << std::setw(40) << "Jeans' length" << LJ << std::endl; - std::cout << std::setw(40) << "Scale height" << h << std::endl; - std::cout << std::setw(40) << "Maximum thickness" << maxZ << std::endl; + if (vm.count("verbose")) { + std::cout << std::endl << "Slab model parameters:" << std::endl; + std::cout << std::setw(40) << "Model type" << modelType << std::endl; + std::cout << std::setw(40) << "Surface mass density" << mu << std::endl; + std::cout << std::setw(40) << "Jeans' length" << LJ << std::endl; + std::cout << std::setw(40) << "Scale height" << h << std::endl; + std::cout << std::setw(40) << "Maximum thickness" << maxZ << std::endl; + } cout.fill(prev); - // // Make mass table // - std::vector Z(Ntable); std::vector M(Ntable); double z, dz = 2.0*maxZ/(Ntable-1.0); @@ -222,34 +215,138 @@ main(int argc, char **argv) std::normal_distribution Vv{0.0, sqrt(DispZ)}; std::normal_distribution Vh{0.0, sqrt(DispX)}; - // Header line - out << std::setw(10) << Number << std::setw(15) << 0 << std::setw(15) << 0 << std::endl; - double KE = 0.0; double VC = 0.0; - double mass = mu/Number; + double mass = mu/Num_particles; - for (int n=0; n; + + ParticleData data; - KE += vel[2]*vel[2]; - VC += pos[2]*model->get_dpot(pos[2]); + data.m.resize(Num_particles); + data.x.resize(Num_particles); + data.y.resize(Num_particles); + data.z.resize(Num_particles); + data.u.resize(Num_particles); + data.v.resize(Num_particles); + data.w.resize(Num_particles); + + data.num_particles = Num_particles; + data.num_aux_ints = Num_aux_ints; + data.num_aux_floats = Num_aux_floats; + + data.aux_ints.resize(Num_aux_ints); + for (auto& vec : data.aux_ints) { + vec.resize(Num_particles); + vec.assign(Num_particles, 0); // Initialize auxiliary integer + // fields to zero + } + + data.aux_floats.resize(Num_aux_floats); + for (auto& vec : data.aux_floats) { + vec.resize(Num_particles); + vec.assign(Num_particles, 0.0f); // Initialize auxiliary float + // fields to zero + } - out << std::setw(15) << mass - << std::setw(15) << pos[0] - << std::setw(15) << pos[1] - << std::setw(15) << pos[2] - << std::setw(15) << vel[0] - << std::setw(15) << vel[1] - << std::setw(15) << vel[2] - << std::endl; + // This is fast enough that multithreading is not strictly + // necessary, but we can use OpenMP to parallelize the particle + // generation loop. The reduction clause is used to accumulate + // KE and VC across threads. +#pragma omp parallel for schedule(dynamic, 256) reduction(+:KE, VC) num_threads(omp_get_max_threads()) + for (int n=0; nget_dpot(data.z[n]); + } + + write_hdf5_data(outfile, data, filter_id, vm.count("verbose") > 0); + }; + + // Write HDF5 with float32 precision. This induces the compiler + // to build the template for float and double precision, allowing + // the user to specifies --double, we instantiate with double type + // instead. + if (vm.count("double")) { + double precision = 1.0; // Placeholder to indicate double precision + create_and_write(precision); + } else { + float precision = 1.0f; // Placeholder to indicate float precision + create_and_write(precision); + } + } + // Old-style ascii output for backwards compatibility. This is not + // recommended at this point (it storage inefficient and slower to + // read), but is retained for users who may have scripts that use + // original ascii format. + else { + + if (vm.count("verbose")) { + std::cout << "Writing ASCII output to " << outfile << std::endl; + } + + std::ofstream out(outfile); + if (!out) { + std::cerr << "Can't open <" << outfile << ">" << std::endl; + exit(-1); + } + out.precision(6); + out.setf(ios::scientific); + + // Header line + out << std::setw(10) << Num_particles << std::setw(15) << 0 << std::setw(15) << 0 << std::endl; + + Eigen::MatrixXd posvel(Num_particles, 6); + + // Generate particles +#pragma omp parallel for schedule(dynamic, 256) reduction(+:KE, VC) num_threads(omp_get_max_threads()) + for (int n=0; nget_dpot(posvel(n, 2)); + } - std::cout << std::endl - << "Virial parameters: KE=" << 0.5*mass*KE - << " VC=" << mass*VC - << " 2T/VC=" << KE/VC << std::endl; + std::cout << "Done generating particles" << std::endl; + + for (int n=0; n Date: Fri, 14 Aug 2026 13:48:55 -0400 Subject: [PATCH 05/28] Add a Python HDF5 IC check script --- tests/check_hdf5_particles.py | 48 +++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 tests/check_hdf5_particles.py diff --git a/tests/check_hdf5_particles.py b/tests/check_hdf5_particles.py new file mode 100644 index 000000000..f1cadac9d --- /dev/null +++ b/tests/check_hdf5_particles.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python3 +"""Check the portable EXP particle HDF5 schema produced by IC generators.""" + +import argparse +import sys + +import h5py +import numpy as np + + +def fail(message): + print(f"HDF5 particle check failed: {message}", file=sys.stderr) + raise SystemExit(1) + + +parser = argparse.ArgumentParser() +parser.add_argument("file") +parser.add_argument("--count", type=int, required=True) +parser.add_argument("--float-bytes", type=int, default=4) +args = parser.parse_args() + +with h5py.File(args.file, "r") as h5: + for name, expected in (("num_particles", args.count), ("num_aux_ints", 0), + ("num_aux_floats", 0)): + if name not in h5.attrs or h5.attrs[name] != expected: + fail(f"attribute {name!r} is not {expected}") + + if "particles" not in h5: + fail("missing /particles group") + particles = h5["particles"] + required = {"m", "x", "y", "z", "u", "v", "w"} + if set(particles) != required: + fail(f"unexpected datasets: {sorted(particles)}") + + for name in required: + data = particles[name] + if data.shape != (args.count,): + fail(f"/particles/{name} has shape {data.shape}, expected ({args.count},)") + if data.dtype.kind != "f" or data.dtype.itemsize != args.float_bytes: + fail(f"/particles/{name} has dtype {data.dtype}, expected float{8 * args.float_bytes}") + if not np.isfinite(data[:]).all(): + fail(f"/particles/{name} contains non-finite values") + + if not np.isclose(particles["m"][:].sum(), 1.0, rtol=2e-6, atol=2e-6): + fail("particle masses do not sum to one") + for name in ("x", "y", "z"): + if (particles[name][:] < 0.0).any() or (particles[name][:] > 1.0).any(): + fail(f"/particles/{name} is outside the unit cube") From ffcf140c5b2e88ef57d8fc47b8f104546b4975d0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 14 Aug 2026 14:10:27 -0400 Subject: [PATCH 06/28] Tweaked the test profiles to allow user-specified total mass and make enclosing cube test optional --- tests/CMakeLists.txt | 50 ++++++++++++++++++++++++++++++++--- tests/check_hdf5_particles.py | 11 +++++--- utils/ICs/CMakeLists.txt | 2 +- 3 files changed, 54 insertions(+), 9 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index afd7a2807..ac9d41406 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -66,7 +66,7 @@ if(ENABLE_NBODY) add_test(NAME expExecuteTest COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/src/exp -v) - if (ENABLE_UTILS) + if(ENABLE_UTILS) # Makes some spherical ICs using utils/ICs/gensph add_test(NAME makeICTest COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph -N 10000 -i SLGridSph.model @@ -160,10 +160,52 @@ if(ENABLE_NBODY) set_tests_properties(removeCubeFiles PROPERTIES DEPENDS expCubeCheckPos REQUIRED_FILES "config.runS.yml;current.processor.rates.runS;cube.bods;OUTLOG.runS;runS.levels;") - if(Python3_Interpreter_FOUND AND EXP_H5PY_STATUS EQUAL 0) + # These are basic sanity checks for writing and reading EXP HDF5 IC files + if(Python3_Interpreter_FOUND AND EXP_H5PY_STATUS EQUAL 0 AND ENABLE_UTILS) + + # Create the banner / message test + add_test(NAME "--- Description ---" + COMMAND ${CMAKE_COMMAND} -E echo "The following tests check for consistent HDF5 IC generation") + + # Check spherical body file generation + add_test(NAME makeSphHDF5ICTest + COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph --hdf5 -N 16 -i SLGridSph.model --SEED 17 --zerovel -o hdf5-sph.bods + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + + add_test(NAME checkSphHDF5ICTest + COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py hdf5-sph.bods.h5 --count 16 + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + set_tests_properties(checkSphHDF5ICTest PROPERTIES DEPENDS makeSphHDF5ICTest) + + add_test(NAME removeSphHDF5ICFiles + COMMAND ${CMAKE_COMMAND} -E remove hdf5-sph.bods.h5 + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + set_tests_properties(removeSphHDF5ICFiles PROPERTIES DEPENDS checkSphHDF5ICTest + REQUIRED_FILES "hdf5-sph.bods.h5") + set_tests_properties(makeSphHDF5ICTest checkSphHDF5ICTest removeSphHDF5ICFiles + PROPERTIES LABELS "quick") + + # Check slab body file generation + add_test(NAME makeSlabHDF5ICTest + COMMAND ${CMAKE_BINARY_DIR}/utils/ICs/slabics --hdf5 --number 16 --seed 17 --cube --outfile hdf5-slab.bods + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Slab) + + add_test(NAME checkSlabHDF5ICTest + COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py hdf5-slab.bods.h5 --count 16 + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Slab) + set_tests_properties(checkSlabHDF5ICTest PROPERTIES DEPENDS makeSlabHDF5ICTest) + + add_test(NAME removeSlabHDF5ICFiles + COMMAND ${CMAKE_COMMAND} -E remove hdf5-slab.bods.h5 + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Slab) + set_tests_properties(removeSlabHDF5ICFiles PROPERTIES DEPENDS checkSlabHDF5ICTest + REQUIRED_FILES "hdf5-slab.bods.h5") + set_tests_properties(makeSlabHDF5ICTest checkSlabHDF5ICTest removeSlabHDF5ICFiles + PROPERTIES LABELS "quick") + + # Check cube body file generation add_test(NAME makeCubeHDF5ICTest - COMMAND ${CMAKE_BINARY_DIR}/utils/ICs/cubeics --hdf5 --number 16 - --seed 17 --zerovel --file hdf5-cube.bods + COMMAND ${CMAKE_BINARY_DIR}/utils/ICs/cubeics --hdf5 --number 16 --seed 17 --cube --zerovel --file hdf5-cube.bods WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube) add_test(NAME checkCubeHDF5ICTest diff --git a/tests/check_hdf5_particles.py b/tests/check_hdf5_particles.py index f1cadac9d..4b432b523 100644 --- a/tests/check_hdf5_particles.py +++ b/tests/check_hdf5_particles.py @@ -17,6 +17,8 @@ def fail(message): parser.add_argument("file") parser.add_argument("--count", type=int, required=True) parser.add_argument("--float-bytes", type=int, default=4) +parser.add_argument("--mass", type=float, default=1.0) +parser.add_argument("--cube", action="store_true", help="check that particles are in the unit cube") args = parser.parse_args() with h5py.File(args.file, "r") as h5: @@ -41,8 +43,9 @@ def fail(message): if not np.isfinite(data[:]).all(): fail(f"/particles/{name} contains non-finite values") - if not np.isclose(particles["m"][:].sum(), 1.0, rtol=2e-6, atol=2e-6): + if not np.isclose(particles["m"][:].sum(), args.mass, rtol=2e-6, atol=2e-6): fail("particle masses do not sum to one") - for name in ("x", "y", "z"): - if (particles[name][:] < 0.0).any() or (particles[name][:] > 1.0).any(): - fail(f"/particles/{name} is outside the unit cube") + if args.cube: + for name in ("x", "y", "z"): + if (particles[name][:] < 0.0).any() or (particles[name][:] > 1.0).any(): + fail(f"/particles/{name} is outside the unit cube") diff --git a/utils/ICs/CMakeLists.txt b/utils/ICs/CMakeLists.txt index 7ac2198cc..c354438b3 100644 --- a/utils/ICs/CMakeLists.txt +++ b/utils/ICs/CMakeLists.txt @@ -87,7 +87,7 @@ add_executable(cubeics cubeICs.cc) add_executable(zangics ZangICs.cc) -add_executable(slabics genslab.cc massmodel1d.cc) +add_executable(slabics genslab.cc) add_executable(addring addring.cc) From 4fe6b8053d62cf6aa138c2dd662258d1e1f6d757 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 14 Aug 2026 14:21:37 -0400 Subject: [PATCH 07/28] Added mass on spherical table as a particle check --- tests/CMakeLists.txt | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index ac9d41406..f4902e188 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -164,16 +164,19 @@ if(ENABLE_NBODY) if(Python3_Interpreter_FOUND AND EXP_H5PY_STATUS EQUAL 0 AND ENABLE_UTILS) # Create the banner / message test - add_test(NAME "--- Description ---" + add_test(NAME HDF5IC COMMAND ${CMAKE_COMMAND} -E echo "The following tests check for consistent HDF5 IC generation") # Check spherical body file generation add_test(NAME makeSphHDF5ICTest - COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph --hdf5 -N 16 -i SLGridSph.model --SEED 17 --zerovel -o hdf5-sph.bods + COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph + --hdf5 -N 16 -i SLGridSph.model --SEED 17 --zerovel -o hdf5-sph.bods WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) add_test(NAME checkSphHDF5ICTest - COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py hdf5-sph.bods.h5 --count 16 + COMMAND ${Python3_EXECUTABLE} + ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py + hdf5-sph.bods.h5 --count 16 --mass 1.008 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) set_tests_properties(checkSphHDF5ICTest PROPERTIES DEPENDS makeSphHDF5ICTest) @@ -187,11 +190,14 @@ if(ENABLE_NBODY) # Check slab body file generation add_test(NAME makeSlabHDF5ICTest - COMMAND ${CMAKE_BINARY_DIR}/utils/ICs/slabics --hdf5 --number 16 --seed 17 --cube --outfile hdf5-slab.bods + COMMAND ${CMAKE_BINARY_DIR}/utils/ICs/slabics --hdf5 --number + 16 --seed 17 --outfile hdf5-slab.bods WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Slab) add_test(NAME checkSlabHDF5ICTest - COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py hdf5-slab.bods.h5 --count 16 + COMMAND ${Python3_EXECUTABLE} + ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py + hdf5-slab.bods.h5 --count 16 --cube WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Slab) set_tests_properties(checkSlabHDF5ICTest PROPERTIES DEPENDS makeSlabHDF5ICTest) @@ -210,7 +216,7 @@ if(ENABLE_NBODY) add_test(NAME checkCubeHDF5ICTest COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py - hdf5-cube.bods.h5 --count 16 + hdf5-cube.bods.h5 --count 16 --cube WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube) set_tests_properties(checkCubeHDF5ICTest PROPERTIES DEPENDS makeCubeHDF5ICTest) From df779e16dea0ac82a258abfc6a08a90feeaad301 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 14 Aug 2026 15:07:18 -0400 Subject: [PATCH 08/28] Add HDF5 file generation to Zang ICs; update the mass tolerance --- tests/CMakeLists.txt | 8 ++--- tests/check_hdf5_particles.py | 2 +- utils/ICs/ZangICs.cc | 61 ++++++++++++++++++++++++++++------- 3 files changed, 54 insertions(+), 17 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f4902e188..a89330a27 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -164,7 +164,7 @@ if(ENABLE_NBODY) if(Python3_Interpreter_FOUND AND EXP_H5PY_STATUS EQUAL 0 AND ENABLE_UTILS) # Create the banner / message test - add_test(NAME HDF5IC + add_test(NAME TestHDF5InitialConditionGeneration COMMAND ${CMAKE_COMMAND} -E echo "The following tests check for consistent HDF5 IC generation") # Check spherical body file generation @@ -191,13 +191,13 @@ if(ENABLE_NBODY) # Check slab body file generation add_test(NAME makeSlabHDF5ICTest COMMAND ${CMAKE_BINARY_DIR}/utils/ICs/slabics --hdf5 --number - 16 --seed 17 --outfile hdf5-slab.bods + 16 --seed 17 --outfile hdf5-slab.bods --model Sech2mu WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Slab) add_test(NAME checkSlabHDF5ICTest COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py - hdf5-slab.bods.h5 --count 16 --cube + hdf5-slab.bods.h5 --count 16 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Slab) set_tests_properties(checkSlabHDF5ICTest PROPERTIES DEPENDS makeSlabHDF5ICTest) @@ -211,7 +211,7 @@ if(ENABLE_NBODY) # Check cube body file generation add_test(NAME makeCubeHDF5ICTest - COMMAND ${CMAKE_BINARY_DIR}/utils/ICs/cubeics --hdf5 --number 16 --seed 17 --cube --zerovel --file hdf5-cube.bods + COMMAND ${CMAKE_BINARY_DIR}/utils/ICs/cubeics --hdf5 --number 16 --seed 17 --zerovel --file hdf5-cube.bods WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube) add_test(NAME checkCubeHDF5ICTest diff --git a/tests/check_hdf5_particles.py b/tests/check_hdf5_particles.py index 4b432b523..6e4a03330 100644 --- a/tests/check_hdf5_particles.py +++ b/tests/check_hdf5_particles.py @@ -43,7 +43,7 @@ def fail(message): if not np.isfinite(data[:]).all(): fail(f"/particles/{name} contains non-finite values") - if not np.isclose(particles["m"][:].sum(), args.mass, rtol=2e-6, atol=2e-6): + if not np.isclose(particles["m"][:].sum(), args.mass, rtol=2e-4, atol=2e-4): fail("particle masses do not sum to one") if args.cube: for name in ("x", "y", "z"): diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index dcc286d0b..6b8d023ff 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -18,6 +18,7 @@ #include "Progress.H" // Progress bar #include "cxxopts.H" // Option parsing +#include "ParticleHDF5.H" // Particle HDF5 output int main(int ac, char **av) @@ -28,12 +29,17 @@ main(int ac, char **av) int N; // Number of particles int Nrepl; // Number of particle replicates per orbit + int NI, ND; // Number of integer and double attributes double mu, nu, Ri, Ro; // Taper paramters double Rmin, Rmax; // Radial range double sigma; // Velocity dispersion std::string bodyfile; // Output file unsigned seed; // Will be inialized by /dev/random if // not set on the command line + unsigned hdf5_filter = 1; // HDF5 filter type (1 = GZIP, 2 = SZIP, 3 = LZF) + bool hdf5_output = false; // Write HDF5 output instead of ASCII + bool hdf5_double = false; // Write HDF5 output in double precision + cxxopts::Options options(av[0], "Ideal tapered Mestel IC generator"); @@ -42,6 +48,10 @@ main(int ac, char **av) ("V,nozerovel", "Do not zero the mean velocity") ("P,nozeropos", "Do not zero the center of mass") ("d,debug", "Print debug grid") + ("5,hdf5", "Write HDF5 phase-space output instead of ASCII") + ("8,double", "Use float64 HDF5 output (default is float32)") + ("Z,filter", "HDF5 filter ID (default: 1 = GZIP)", + cxxopts::value(hdf5_filter)->default_value("1")) ("N,number", "Number of particles to generate", cxxopts::value(N)->default_value("100000")) ("n,nu", "Inner taper exponent (0 for no taper)", @@ -64,6 +74,10 @@ main(int ac, char **av) cxxopts::value(Nrepl)->default_value("1")) ("f,file", "Output body file", cxxopts::value(bodyfile)->default_value("zang.bods")) + ("NI", "Number of integer attributes for PSP file", + cxxopts::value(NI)->default_value("0")) + ("ND", "Number of double attributes for PSP file", + cxxopts::value(ND)->default_value("0")) ; cxxopts::ParseResult vm; @@ -82,6 +96,9 @@ main(int ac, char **av) return 1; } + hdf5_output = vm.count("hdf5") > 0; + hdf5_double = vm.count("double") > 0; + // Set from /dev/random if not specified if (vm.count("seed")==0) { seed = std::random_device{}(); @@ -99,11 +116,14 @@ main(int ac, char **av) // Open the output file // - std::ofstream out(bodyfile); - if (not out) { - std::string msg(av[0]); - msg += ": output file <" + bodyfile + "> can not be opened"; - throw std::runtime_error(msg); + std::ofstream out; + if (not hdf5_output) { + out.open(bodyfile); + if (not out) { + std::string msg(av[0]); + msg += ": output file <" + bodyfile + "> can not be opened"; + throw std::runtime_error(msg); + } } SphericalOrbit::ZFRAC=0.3; // TEST @@ -294,16 +314,27 @@ main(int ac, char **av) std::cout << "** " << over << " particles failed acceptance" << std::endl << "** Particle mass=" << mass << std::endl; - out << std::setw(8) << N << std::setw(8) << 0 << std::setw(8) << 0 - << std::endl; + if (not hdf5_output) { + out << std::setw(8) << N << std::setw(8) << 0 << std::setw(8) << 0 + << std::endl; + } double ektot = 0.0, clausius = 0.0; - for (int n=0; n> hdf5_particles; + for (int n=0; n Date: Fri, 14 Aug 2026 16:09:42 -0400 Subject: [PATCH 09/28] Allow use of non-MPI applications --- utils/ICs/ParticleHDF5.H | 96 +++++++++++++++++++++++++++++----------- 1 file changed, 70 insertions(+), 26 deletions(-) diff --git a/utils/ICs/ParticleHDF5.H b/utils/ICs/ParticleHDF5.H index dc7fe26d7..312dcc674 100644 --- a/utils/ICs/ParticleHDF5.H +++ b/utils/ICs/ParticleHDF5.H @@ -138,35 +138,79 @@ namespace EXP unsigned filter, bool double_precision, MPI_Comm comm = MPI_COMM_WORLD) { - int rank, ranks, local_count = static_cast(local.size()); - MPI_Comm_rank(comm, &rank); MPI_Comm_size(comm, &ranks); - std::vector counts(rank == 0 ? ranks : 0); - MPI_Gather(&local_count, 1, MPI_INT, rank == 0 ? counts.data() : nullptr, 1, MPI_INT, 0, comm); - - std::vector send(7 * local.size()); - for (size_t i = 0; i < local.size(); ++i) - std::copy(local[i].begin(), local[i].end(), send.begin() + 7*i); - std::vector counts7, offsets; - std::vector received; - if (rank == 0) { - counts7.resize(ranks); offsets.resize(ranks); - int total = 0; - for (int i = 0; i < ranks; ++i) { offsets[i] = total; counts7[i] = 7 * counts[i]; total += counts7[i]; } - received.resize(total); + // Check if MPI is initialized + // + int is_init = 0; + MPI_Initialized(&is_init); + + // Gather the local counts to rank 0 only if MPI is initialized, otherwise just use the local data + // + if (is_init) { + int myid=0, ranks=1, local_count = static_cast(local.size()); + // Get the rank and size of the communicator + MPI_Comm_rank(comm, &myid ); + MPI_Comm_size(comm, &ranks); + + // Gather the local counts to rank 0 + std::vector counts(myid == 0 ? ranks : 0); + + MPI_Gather(&local_count, 1, MPI_INT, + myid == 0 ? counts.data() : nullptr, + 1, MPI_INT, 0, comm); + + // Gather the local data to rank 0 + std::vector send(7 * local.size()); + for (size_t i = 0; i < local.size(); ++i) + std::copy(local[i].begin(), local[i].end(), send.begin() + 7*i); + + std::vector counts7, offsets; + std::vector received; + + // Only Rank 0 needs to know the counts and offsets for the gather + // + if (myid == 0) { + counts7.resize(ranks); + offsets.resize(ranks); + int total = 0; + for (int i = 0; i < ranks; ++i) { + offsets[i] = total; + counts7[i] = 7 * counts[i]; + total += counts7[i]; + } + received.resize(total); + } + + // Do the gather of the 7-element arrays to Rank 0 + // + MPI_Gatherv(send.data(), static_cast(send.size()), MPI_DOUBLE, + myid == 0 ? received.data() : nullptr, + myid == 0 ? counts7. data() : nullptr, + myid == 0 ? offsets. data() : nullptr, + MPI_DOUBLE, 0, comm); + + if (myid != 0) return; + + // Rank 0 reconstructs the all vector from the received data + // + std::vector> all(received.size()/7); + for (size_t i = 0; i < all.size(); ++i) + std::copy(received.begin() + 7*i, received.begin() + 7*(i+1), all[i].begin()); + + // Write the gathered data to the HDF5 file + // + if (double_precision) write(filename, all, num_aux_ints, num_aux_floats, filter); + else write(filename, all, num_aux_ints, num_aux_floats, filter); + + } else { + // Write the local/only data to the HDF5 file + // + if (double_precision) write(filename, local, num_aux_ints, num_aux_floats, filter); + else write(filename, local, num_aux_ints, num_aux_floats, filter); } - MPI_Gatherv(send.data(), static_cast(send.size()), MPI_DOUBLE, - rank == 0 ? received.data() : nullptr, - rank == 0 ? counts7.data() : nullptr, rank == 0 ? offsets.data() : nullptr, - MPI_DOUBLE, 0, comm); - if (rank != 0) return; - - std::vector> all(received.size()/7); - for (size_t i = 0; i < all.size(); ++i) - std::copy(received.begin() + 7*i, received.begin() + 7*(i+1), all[i].begin()); - if (double_precision) write(filename, all, num_aux_ints, num_aux_floats, filter); - else write(filename, all, num_aux_ints, num_aux_floats, filter); } + } // namespace ParticleHDF5 + } // namespace EXP #endif From 18a0ce9d1bf1192c759838123a6116d1084302cc Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 14 Aug 2026 16:36:22 -0400 Subject: [PATCH 10/28] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/ICs/ParticleHDF5.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/ICs/ParticleHDF5.H b/utils/ICs/ParticleHDF5.H index 312dcc674..768e0dcbc 100644 --- a/utils/ICs/ParticleHDF5.H +++ b/utils/ICs/ParticleHDF5.H @@ -44,7 +44,7 @@ namespace EXP HighFive::DataSetCreateProps props; // chunk size is min(size, 256kB) but at least 1kB props.add(HighFive::Chunking({std::min(size, std::min(262144, std::max(1024, size)))})); - if (filter != 3 && filter != 4 && filter != 32001) props.add(HighFive::Shuffle()); + if (filter != 4 && filter != 32001) props.add(HighFive::Shuffle()); // The filter IDs are from HDF5 1.12.0 and later, see https://support.hdfgroup.org/HDF5/doc/Advanced/Filters.html // Filter types: From 828b00a46f0fd82033c881ddf6bb08277c8cba81 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 14 Aug 2026 16:36:50 -0400 Subject: [PATCH 11/28] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/ICs/ParticleHDF5.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/ICs/ParticleHDF5.H b/utils/ICs/ParticleHDF5.H index 768e0dcbc..53794ac94 100644 --- a/utils/ICs/ParticleHDF5.H +++ b/utils/ICs/ParticleHDF5.H @@ -52,8 +52,8 @@ namespace EXP // 2: SZIP // 3: SHUFFLE // 4: FLETCHER32 - // 307: LZ4 - // 32004: ZSTD + // 307: BZIP2 + // 32004: LZ4 // 32001: BLOSC // // The filter options are specific to each filter type, see the HDF5 documentation for details. From 4d2c233086da65090f1ebfe6dc4503e150b98e8a Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 14 Aug 2026 16:37:31 -0400 Subject: [PATCH 12/28] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/ICs/ParticleHDF5.H | 2 ++ 1 file changed, 2 insertions(+) diff --git a/utils/ICs/ParticleHDF5.H b/utils/ICs/ParticleHDF5.H index 53794ac94..4fdf2fb00 100644 --- a/utils/ICs/ParticleHDF5.H +++ b/utils/ICs/ParticleHDF5.H @@ -126,6 +126,8 @@ namespace EXP m.push_back(p[0]); x.push_back(p[1]); y.push_back(p[2]); z.push_back(p[3]); u.push_back(p[4]); v.push_back(p[5]); w.push_back(p[6]); } + if (num_aux_ints < 0 || num_aux_floats < 0) + throw std::invalid_argument("auxiliary field counts must be non-negative"); std::vector> auxi(num_aux_ints, std::vector(particles.size(), 0)); std::vector> auxf(num_aux_floats, std::vector(particles.size(), 0)); write(filename, m, x, y, z, u, v, w, auxi, auxf, std::nullopt, filter); From 658a79942bd7c08479984007fe9ea0bb54efb6e1 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 14 Aug 2026 16:38:43 -0400 Subject: [PATCH 13/28] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/ICs/genslab.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/utils/ICs/genslab.cc b/utils/ICs/genslab.cc index fc2282a39..0c7ab896c 100644 --- a/utils/ICs/genslab.cc +++ b/utils/ICs/genslab.cc @@ -217,6 +217,8 @@ main(int argc, char **argv) double KE = 0.0; double VC = 0.0; + if (Num_particles <= 0) + throw std::invalid_argument("number of particles must be positive"); double mass = mu/Num_particles; if (vm.count("verbose")) { From f4bd1b70d9693f1676238478f4a909a81f8372b3 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 14 Aug 2026 16:40:52 -0400 Subject: [PATCH 14/28] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/ICs/genslab.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/utils/ICs/genslab.cc b/utils/ICs/genslab.cc index 0c7ab896c..f02ea0e33 100644 --- a/utils/ICs/genslab.cc +++ b/utils/ICs/genslab.cc @@ -251,7 +251,9 @@ main(int argc, char **argv) data.num_particles = Num_particles; data.num_aux_ints = Num_aux_ints; data.num_aux_floats = Num_aux_floats; - + + if (Num_aux_ints < 0 || Num_aux_floats < 0) + throw std::invalid_argument("auxiliary field counts must be non-negative"); data.aux_ints.resize(Num_aux_ints); for (auto& vec : data.aux_ints) { vec.resize(Num_particles); From bff80dde63a3ab49182852241ec332e6bab44626 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 14 Aug 2026 16:41:22 -0400 Subject: [PATCH 15/28] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- tests/check_hdf5_particles.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/check_hdf5_particles.py b/tests/check_hdf5_particles.py index 6e4a03330..3c3c00aea 100644 --- a/tests/check_hdf5_particles.py +++ b/tests/check_hdf5_particles.py @@ -44,7 +44,7 @@ def fail(message): fail(f"/particles/{name} contains non-finite values") if not np.isclose(particles["m"][:].sum(), args.mass, rtol=2e-4, atol=2e-4): - fail("particle masses do not sum to one") + fail(f"particle masses do not sum to {args.mass}") if args.cube: for name in ("x", "y", "z"): if (particles[name][:] < 0.0).any() or (particles[name][:] > 1.0).any(): From 103f3e09eb89895f53e6be6db7753a96c41a99c7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 14 Aug 2026 21:11:18 +0000 Subject: [PATCH 16/28] Add gendisk HDF5 schema CTest coverage Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- tests/CMakeLists.txt | 35 ++++++++++++++++++ .../check_hdf5_particles.cpython-312.pyc | Bin 0 -> 4339 bytes tests/check_hdf5_particles.py | 16 ++++++-- 3 files changed, 47 insertions(+), 4 deletions(-) create mode 100644 tests/__pycache__/check_hdf5_particles.cpython-312.pyc diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index a89330a27..192db3d0f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -227,6 +227,41 @@ if(ENABLE_NBODY) REQUIRED_FILES "hdf5-cube.bods.h5") set_tests_properties(makeCubeHDF5ICTest checkCubeHDF5ICTest removeCubeHDF5ICFiles PROPERTIES LABELS "quick") + + # Check gendisk body file generation, including gas auxiliary fields + add_test(NAME makeGendiskHDF5ICTest + COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gendisk + --hdf5 --nhalo 8 --ndisk 8 --ngas 8 --ngparam 2 --SEED 17 + --hbods hdf5-gendisk-halo.bods --dbods hdf5-gendisk-disk.bods + --gbods hdf5-gendisk-gas.bods --cachefile hdf5-gendisk.cache + --suffix hdf5-gendisk + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + + add_test(NAME checkGendiskDiskHDF5ICTest + COMMAND ${Python3_EXECUTABLE} + ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py + hdf5-gendisk-disk.bods.h5 --count 8 + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + set_tests_properties(checkGendiskDiskHDF5ICTest PROPERTIES DEPENDS makeGendiskHDF5ICTest) + + add_test(NAME checkGendiskGasHDF5ICTest + COMMAND ${Python3_EXECUTABLE} + ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py + hdf5-gendisk-gas.bods.h5 --count 8 --aux-floats 2 + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + set_tests_properties(checkGendiskGasHDF5ICTest PROPERTIES DEPENDS makeGendiskHDF5ICTest) + + add_test(NAME removeGendiskHDF5ICFiles + COMMAND ${CMAKE_COMMAND} -E remove + hdf5-gendisk-halo.bods.h5 hdf5-gendisk-disk.bods.h5 hdf5-gendisk-gas.bods.h5 + hdf5-gendisk.cache hdf5-gendisk.ortho_check + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + set_tests_properties(removeGendiskHDF5ICFiles PROPERTIES DEPENDS + "checkGendiskDiskHDF5ICTest;checkGendiskGasHDF5ICTest" + REQUIRED_FILES "hdf5-gendisk-disk.bods.h5;hdf5-gendisk-gas.bods.h5") + set_tests_properties(makeGendiskHDF5ICTest checkGendiskDiskHDF5ICTest + checkGendiskGasHDF5ICTest removeGendiskHDF5ICFiles + PROPERTIES LABELS "quick") endif() # Set labels for pyEXP tests diff --git a/tests/__pycache__/check_hdf5_particles.cpython-312.pyc b/tests/__pycache__/check_hdf5_particles.cpython-312.pyc new file mode 100644 index 0000000000000000000000000000000000000000..b4e3e7fef8d1b186131a2552b44c73da6e849950 GIT binary patch literal 4339 zcma(TTWk|o_Re@5JGK)W;yiEyojeFJ!K8%%Aq7H6nuQi6VF_Ju%UXl)O>E;wIx~iN z`hY&9CR!vyf9h3*)={NhM<5ZtV5OC|tKC*qmA2y~kTC5ByL_SjiQuDZzxLd5CX;1D z%OrEpx#w}u^PcfHWn~ToPy0s~VxN~H^pI|h$GMC=iP#YO6k!xZn89qE2{JllgRBlM zL5mK#pcPv_LctOi5hhk{9Am<3%V^L}?6C@yZeW|gY%f@kSl*qIXCv;b=TMbkAEACJeY7LV`Y$2ORUc17#JywZ8C5Z z;)-i*akmQvJ=hCU4JOyXhFM$*oF-#$Lf}*Zr`e@R2!P4hGt4!r^mku);g0%y)@*gZ@p-eIbBaH zLm;Dib~ciYv3*7^O!$nRV02r}x~&(WCJ6Iu>jS5L?K6PWaMiZTD{rjJv}<6v&X5py zf`x+y-Gs42CT?K3QTOk?RsJ=ttqopmzG^jWFer)FmRr+Gv)-GH&(P>KyfGnUqR-qJ znC@%KbNRXjx2%1JVnwtUOUy{%JZtQ*VZ(%EBW5w{Tj|x6Ri62; zT1*STYvPb@=cu7+LQ;*l(#|Sciw&3@12e2ZM1uy+gt!&AUf*_&D^jFN-HX85vQPDQ zS?t&5IR+;7;q5S!x>t97(Kp<77N%6}gdu4{!$Zku+>>Pu*IB6U&uTE2cCX6Yi>vc? zsRyPN)Bv-_{AEyv4YQ0sL(ao$V2$V{n0@DAe+n`1mtn^&r_g!E|B}+z6CvTVyc{9? zXiAiYGjYPda&nL#6+}52hRxyLzFoW&j*x`FkBTXr4in6u$?!*d_z{vMq9CV4Y5OA@ zY=3z9i>gOXGfuypPljOm6v-<*%oJyl4B09^vE}L*dkU|wib{q#S$yzH-+aLMxqhWi zcW=d6a1(+V7NT*2J9$5&al?QyRb4Il>U^o{*zuQ&zxw{{yUX={@9<6@`&~fgs-J)q5@{+FTObfl zrIYfIpVc@yGfFhONZv|EMS_*GKwvnY66C;{3^?(K%KKT7{sPi$m<$W)xU7@}0trEq zMt}&)S?c5!r+ILj;^ZTO#7hxj6#Q{*<5v_~6n2|mW1`Cv zm==Xse5)jcBxFY>BShr&iZR2AXBAs_;*lMER%-chDk%%mB#ky17>+_V34czAr->x` zKw>1&Q}>^L@7~%)%P-2QxW>UXiz}Oq6!5gVMB;Hu&zF==@N$axTSWTvOniYZJUs_q z{ZTg1)TYw16vgDJ3btytPza~Op%COv903y&5%t`1kzLTsR;U;NHaV&{v{v3Nj-(SL zDG%zBT1inua|sxS-~gScS?Eo(>OpG`U4rr?2t;4|;9BafuCb9_npJQ1m9j~z$_5%q zAp!kwR6?`Tj*&$AO%UC%XeH@UC@0aZqL3URaJ`8P*Bgcz!5T*kty%Sg(iVrQCY_JU z5zVSsP~*-1A>ai*>E}p&P;Jpjz8;`@o&{Z})$_U+p}#Shp>Y zZ=|MEYV83?LhV83F|t@YC(6gU@yG+Gds6&#?8;cd`9f~TpE{;G^1V|BW=_sE?!WKy zOm=2*_GIp4K9Ucr>pJeb_Ac2#6*6Y8$Tk#gb@%NqwQS40y=Bp}e)8oj?u$n>_l}E) z7dh9*1Md$^pO|UCwg2Y+xj^?l@1Z<1IhJFx(kH*2^&Xn%dJLiT%+RgS&Cp!iEBCy8 z`T9vESD&kxt=v58?VIO*b>HR~-}~Wv*_v5fJ+V0^L|3+R-r2n9uDGmBD3>oxT*x)fx;NkVcqiY^wdb1i z&Y9j@uikuhW}n)1O05ju^}O-mDcQXDhHu)p;NMs9@B4yVI5=21I5_JcT=R$lb?}%v zbaG(`7vO&os{=%>9G>@#EV_N;;^ns|-p;m9G1)zH?#2c8mV$fBJ-2@)8nvlQU3c)F zt9z-MreMj3Tvg*cXzY{GoaIV#>X=%!WudCIP}MqD70AmsGSivc@|~u;&H+7}0cefr zbTt2lTDEuI-nnRVew_Irv*>MnjM&m%X1sUNvmx7b#XZHSo`&4VI*o!KA zm>)H7-w)L$`6uZh+}FQv?Dq7xSpKoSn*;FW4p;vs$Cn410spGW)!*p&>c!r@u=%F- VP%VIe+REv=$wAi^YyV5O{{p|x3H1N~ literal 0 HcmV?d00001 diff --git a/tests/check_hdf5_particles.py b/tests/check_hdf5_particles.py index 3c3c00aea..9dbb1a18b 100644 --- a/tests/check_hdf5_particles.py +++ b/tests/check_hdf5_particles.py @@ -18,12 +18,14 @@ def fail(message): parser.add_argument("--count", type=int, required=True) parser.add_argument("--float-bytes", type=int, default=4) parser.add_argument("--mass", type=float, default=1.0) +parser.add_argument("--aux-ints", type=int, default=0) +parser.add_argument("--aux-floats", type=int, default=0) parser.add_argument("--cube", action="store_true", help="check that particles are in the unit cube") args = parser.parse_args() with h5py.File(args.file, "r") as h5: - for name, expected in (("num_particles", args.count), ("num_aux_ints", 0), - ("num_aux_floats", 0)): + for name, expected in (("num_particles", args.count), ("num_aux_ints", args.aux_ints), + ("num_aux_floats", args.aux_floats)): if name not in h5.attrs or h5.attrs[name] != expected: fail(f"attribute {name!r} is not {expected}") @@ -31,6 +33,8 @@ def fail(message): fail("missing /particles group") particles = h5["particles"] required = {"m", "x", "y", "z", "u", "v", "w"} + required.update({f"aux_int_{i}" for i in range(args.aux_ints)}) + required.update({f"aux_float_{i}" for i in range(args.aux_floats)}) if set(particles) != required: fail(f"unexpected datasets: {sorted(particles)}") @@ -38,8 +42,12 @@ def fail(message): data = particles[name] if data.shape != (args.count,): fail(f"/particles/{name} has shape {data.shape}, expected ({args.count},)") - if data.dtype.kind != "f" or data.dtype.itemsize != args.float_bytes: - fail(f"/particles/{name} has dtype {data.dtype}, expected float{8 * args.float_bytes}") + if name.startswith("aux_int_"): + if data.dtype.kind != "i": + fail(f"/particles/{name} has dtype {data.dtype}, expected an integer type") + else: + if data.dtype.kind != "f" or data.dtype.itemsize != args.float_bytes: + fail(f"/particles/{name} has dtype {data.dtype}, expected float{8 * args.float_bytes}") if not np.isfinite(data[:]).all(): fail(f"/particles/{name} contains non-finite values") From 86fc6c71f72d2b60be9852fe908649da33707fff Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 14 Aug 2026 22:26:54 +0000 Subject: [PATCH 17/28] Add gendisk2d HDF5 CTest coverage Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- tests/CMakeLists.txt | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 192db3d0f..4185ceb81 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -262,6 +262,40 @@ if(ENABLE_NBODY) set_tests_properties(makeGendiskHDF5ICTest checkGendiskDiskHDF5ICTest checkGendiskGasHDF5ICTest removeGendiskHDF5ICFiles PROPERTIES LABELS "quick") + + # Check gendisk2d body file generation + add_test(NAME makeGendisk2dHDF5ICTest + COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gendisk2d + --hdf5 --nhalo 8 --ndisk 8 --SEED 17 + --hbods hdf5-gendisk2d-halo.bods --dbods hdf5-gendisk2d-disk.bods + --cachefile hdf5-gendisk2d.cache + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + + add_test(NAME checkGendisk2dHaloHDF5ICTest + COMMAND ${Python3_EXECUTABLE} + ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py + hdf5-gendisk2d-halo.bods.h5 --count 8 + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + set_tests_properties(checkGendisk2dHaloHDF5ICTest PROPERTIES DEPENDS makeGendisk2dHDF5ICTest) + + add_test(NAME checkGendisk2dDiskHDF5ICTest + COMMAND ${Python3_EXECUTABLE} + ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py + hdf5-gendisk2d-disk.bods.h5 --count 8 + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + set_tests_properties(checkGendisk2dDiskHDF5ICTest PROPERTIES DEPENDS makeGendisk2dHDF5ICTest) + + add_test(NAME removeGendisk2dHDF5ICFiles + COMMAND ${CMAKE_COMMAND} -E remove + hdf5-gendisk2d-halo.bods.h5 hdf5-gendisk2d-disk.bods.h5 + hdf5-gendisk2d.cache hdf5-gendisk2d.ortho_check + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + set_tests_properties(removeGendisk2dHDF5ICFiles PROPERTIES DEPENDS + "checkGendisk2dHaloHDF5ICTest;checkGendisk2dDiskHDF5ICTest" + REQUIRED_FILES "hdf5-gendisk2d-halo.bods.h5;hdf5-gendisk2d-disk.bods.h5") + set_tests_properties(makeGendisk2dHDF5ICTest checkGendisk2dHaloHDF5ICTest + checkGendisk2dDiskHDF5ICTest removeGendisk2dHDF5ICFiles + PROPERTIES LABELS "quick") endif() # Set labels for pyEXP tests From 1f52901d45c03b09582ca347b34c7adbff222b89 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 14 Aug 2026 22:29:33 +0000 Subject: [PATCH 18/28] test: add gendisk halo hdf5 schema check Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- tests/CMakeLists.txt | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4185ceb81..575ac8d9c 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -244,6 +244,13 @@ if(ENABLE_NBODY) WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) set_tests_properties(checkGendiskDiskHDF5ICTest PROPERTIES DEPENDS makeGendiskHDF5ICTest) + add_test(NAME checkGendiskHaloHDF5ICTest + COMMAND ${Python3_EXECUTABLE} + ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py + hdf5-gendisk-halo.bods.h5 --count 8 + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) + set_tests_properties(checkGendiskHaloHDF5ICTest PROPERTIES DEPENDS makeGendiskHDF5ICTest) + add_test(NAME checkGendiskGasHDF5ICTest COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py @@ -257,9 +264,9 @@ if(ENABLE_NBODY) hdf5-gendisk.cache hdf5-gendisk.ortho_check WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) set_tests_properties(removeGendiskHDF5ICFiles PROPERTIES DEPENDS - "checkGendiskDiskHDF5ICTest;checkGendiskGasHDF5ICTest" - REQUIRED_FILES "hdf5-gendisk-disk.bods.h5;hdf5-gendisk-gas.bods.h5") - set_tests_properties(makeGendiskHDF5ICTest checkGendiskDiskHDF5ICTest + "checkGendiskDiskHDF5ICTest;checkGendiskHaloHDF5ICTest;checkGendiskGasHDF5ICTest" + REQUIRED_FILES "hdf5-gendisk-halo.bods.h5;hdf5-gendisk-disk.bods.h5;hdf5-gendisk-gas.bods.h5") + set_tests_properties(makeGendiskHDF5ICTest checkGendiskDiskHDF5ICTest checkGendiskHaloHDF5ICTest checkGendiskGasHDF5ICTest removeGendiskHDF5ICFiles PROPERTIES LABELS "quick") From 148025382f9dc61881ccb4eba50298670b59bd36 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 14 Aug 2026 22:30:05 +0000 Subject: [PATCH 19/28] Reduce root memory overhead in gathered HDF5 writes Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- Testing/Temporary/LastTest.log | 3 ++ utils/ICs/ParticleHDF5.H | 50 ++++++++++++++++++++++++++++------ 2 files changed, 45 insertions(+), 8 deletions(-) create mode 100644 Testing/Temporary/LastTest.log diff --git a/Testing/Temporary/LastTest.log b/Testing/Temporary/LastTest.log new file mode 100644 index 000000000..bae840da5 --- /dev/null +++ b/Testing/Temporary/LastTest.log @@ -0,0 +1,3 @@ +Start testing: Aug 14 22:26 UTC +---------------------------------------------------------- +End testing: Aug 14 22:26 UTC diff --git a/utils/ICs/ParticleHDF5.H b/utils/ICs/ParticleHDF5.H index 4fdf2fb00..adb841ab8 100644 --- a/utils/ICs/ParticleHDF5.H +++ b/utils/ICs/ParticleHDF5.H @@ -133,6 +133,46 @@ namespace EXP write(filename, m, x, y, z, u, v, w, auxi, auxf, std::nullopt, filter); } + //! Write a particle phase-space file from an interleaved [m,x,y,z,u,v,w] buffer. + template + inline void write_interleaved(const std::string& filename, + const std::vector& particles, + int num_aux_ints, int num_aux_floats, + unsigned filter) + { + if (particles.empty()) return; + if (particles.size() % 7 != 0) + throw std::invalid_argument("interleaved phase-space buffer must have a multiple-of-7 length"); + if (num_aux_ints < 0 || num_aux_floats < 0) + throw std::invalid_argument("auxiliary field counts must be non-negative"); + + const auto count = particles.size() / 7; + const int n = static_cast(count); + + HighFive::File file(filename, HighFive::File::ReadWrite | HighFive::File::Create | HighFive::File::Truncate); + file.createAttribute("num_particles", HighFive::DataSpace::From(n)).write(n); + file.createAttribute("num_aux_ints", HighFive::DataSpace::From(num_aux_ints)).write(num_aux_ints); + file.createAttribute("num_aux_floats", HighFive::DataSpace::From(num_aux_floats)).write(num_aux_floats); + + auto group = file.createGroup("particles"); + auto props = properties(count, filter); + + static constexpr std::array names = {"m", "x", "y", "z", "u", "v", "w"}; + std::vector field(count); + for (size_t col = 0; col < names.size(); ++col) { + for (size_t i = 0; i < count; ++i) field[i] = static_cast(particles[7*i + col]); + group.createDataSet(names[col], field, props); + } + + std::vector auxi(count, 0); + for (int j = 0; j < num_aux_ints; ++j) + group.createDataSet("aux_int_" + std::to_string(j), auxi, props); + + std::vector auxf(count, 0); + for (int j = 0; j < num_aux_floats; ++j) + group.createDataSet("aux_float_" + std::to_string(j), auxf, props); + } + //! Collect rank-local phase space records and have rank zero write one file. inline void gather_and_write(const std::string& filename, const std::vector>& local, @@ -192,16 +232,10 @@ namespace EXP if (myid != 0) return; - // Rank 0 reconstructs the all vector from the received data - // - std::vector> all(received.size()/7); - for (size_t i = 0; i < all.size(); ++i) - std::copy(received.begin() + 7*i, received.begin() + 7*(i+1), all[i].begin()); - // Write the gathered data to the HDF5 file // - if (double_precision) write(filename, all, num_aux_ints, num_aux_floats, filter); - else write(filename, all, num_aux_ints, num_aux_floats, filter); + if (double_precision) write_interleaved(filename, received, num_aux_ints, num_aux_floats, filter); + else write_interleaved(filename, received, num_aux_ints, num_aux_floats, filter); } else { // Write the local/only data to the HDF5 file From fd116a9bb5073e2db49eede9bb0e5f81bd2e4f9c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 15 Aug 2026 10:58:51 -0400 Subject: [PATCH 20/28] Removed gas particle tests, decrease basis resolution to speed up tests, added corrected total mass checks --- tests/CMakeLists.txt | 54 ++++++++++++++++++++++---------------------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 575ac8d9c..8e3ac5e91 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -69,7 +69,7 @@ if(ENABLE_NBODY) if(ENABLE_UTILS) # Makes some spherical ICs using utils/ICs/gensph add_test(NAME makeICTest - COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph -N 10000 -i SLGridSph.model + COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph -N 10000 -i SLGridSph.model --NUMG 100 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) # Runs those ICs using exp @@ -171,6 +171,7 @@ if(ENABLE_NBODY) add_test(NAME makeSphHDF5ICTest COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph --hdf5 -N 16 -i SLGridSph.model --SEED 17 --zerovel -o hdf5-sph.bods + --NUMG 100 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) add_test(NAME checkSphHDF5ICTest @@ -186,7 +187,7 @@ if(ENABLE_NBODY) set_tests_properties(removeSphHDF5ICFiles PROPERTIES DEPENDS checkSphHDF5ICTest REQUIRED_FILES "hdf5-sph.bods.h5") set_tests_properties(makeSphHDF5ICTest checkSphHDF5ICTest removeSphHDF5ICFiles - PROPERTIES LABELS "quick") + PROPERTIES LABELS "long") # Check slab body file generation add_test(NAME makeSlabHDF5ICTest @@ -228,52 +229,51 @@ if(ENABLE_NBODY) set_tests_properties(makeCubeHDF5ICTest checkCubeHDF5ICTest removeCubeHDF5ICFiles PROPERTIES LABELS "quick") - # Check gendisk body file generation, including gas auxiliary fields + # Check gendisk body file generation with a pared-down + # low-resolution basis add_test(NAME makeGendiskHDF5ICTest COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gendisk - --hdf5 --nhalo 8 --ndisk 8 --ngas 8 --ngparam 2 --SEED 17 + --nhalo 10000 --ndisk 10000 + --hdf5 --SEED 17 --ignore --MMAX 2 --NUMX 128 --NUMY 64 + --LMAXFID 32 --NMAXFID 24 --PNUM 0 --TNUM 32 --hbods hdf5-gendisk-halo.bods --dbods hdf5-gendisk-disk.bods - --gbods hdf5-gendisk-gas.bods --cachefile hdf5-gendisk.cache - --suffix hdf5-gendisk + --cachefile hdf5-gendisk.cache --suffix hdf5-gendisk WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) add_test(NAME checkGendiskDiskHDF5ICTest COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py - hdf5-gendisk-disk.bods.h5 --count 8 + hdf5-gendisk-disk.bods.h5 --count 10000 --mass 0.05 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) set_tests_properties(checkGendiskDiskHDF5ICTest PROPERTIES DEPENDS makeGendiskHDF5ICTest) add_test(NAME checkGendiskHaloHDF5ICTest COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py - hdf5-gendisk-halo.bods.h5 --count 8 + hdf5-gendisk-halo.bods.h5 --count 10000 --mass 1.0079 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) set_tests_properties(checkGendiskHaloHDF5ICTest PROPERTIES DEPENDS makeGendiskHDF5ICTest) - add_test(NAME checkGendiskGasHDF5ICTest - COMMAND ${Python3_EXECUTABLE} - ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py - hdf5-gendisk-gas.bods.h5 --count 8 --aux-floats 2 - WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) - set_tests_properties(checkGendiskGasHDF5ICTest PROPERTIES DEPENDS makeGendiskHDF5ICTest) - add_test(NAME removeGendiskHDF5ICFiles COMMAND ${CMAKE_COMMAND} -E remove - hdf5-gendisk-halo.bods.h5 hdf5-gendisk-disk.bods.h5 hdf5-gendisk-gas.bods.h5 + hdf5-gendisk-halo.bods.h5 hdf5-gendisk-disk.bods.h5 hdf5-gendisk.cache hdf5-gendisk.ortho_check WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) set_tests_properties(removeGendiskHDF5ICFiles PROPERTIES DEPENDS - "checkGendiskDiskHDF5ICTest;checkGendiskHaloHDF5ICTest;checkGendiskGasHDF5ICTest" - REQUIRED_FILES "hdf5-gendisk-halo.bods.h5;hdf5-gendisk-disk.bods.h5;hdf5-gendisk-gas.bods.h5") - set_tests_properties(makeGendiskHDF5ICTest checkGendiskDiskHDF5ICTest checkGendiskHaloHDF5ICTest - checkGendiskGasHDF5ICTest removeGendiskHDF5ICFiles - PROPERTIES LABELS "quick") - - # Check gendisk2d body file generation + "checkGendiskDiskHDF5ICTest;checkGendiskHaloHDF5ICTest" + REQUIRED_FILES "hdf5-gendisk-halo.bods.h5;hdf5-gendisk-disk.bods.h5") + set_tests_properties(makeGendiskHDF5ICTest + checkGendiskDiskHDF5ICTest checkGendiskHaloHDF5ICTest + removeGendiskHDF5ICFiles + PROPERTIES LABELS "long") + + # Check gendisk2d body file generation with a pared-down + # low-resolution basis add_test(NAME makeGendisk2dHDF5ICTest COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gendisk2d - --hdf5 --nhalo 8 --ndisk 8 --SEED 17 + --hdf5 --nhalo 10000 --ndisk 10000 --SEED 17 --NUMX 128 --NUMY 64 + --PNUM 0 --TNUM 40 --LMAX 4 --MMAX 2 --NMAXH 4 --NMAXD 4 + # --RMIN 3.333e-05 --hbods hdf5-gendisk2d-halo.bods --dbods hdf5-gendisk2d-disk.bods --cachefile hdf5-gendisk2d.cache WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) @@ -281,14 +281,14 @@ if(ENABLE_NBODY) add_test(NAME checkGendisk2dHaloHDF5ICTest COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py - hdf5-gendisk2d-halo.bods.h5 --count 8 + hdf5-gendisk2d-halo.bods.h5 --count 10000 --mass 1.0079 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) set_tests_properties(checkGendisk2dHaloHDF5ICTest PROPERTIES DEPENDS makeGendisk2dHDF5ICTest) add_test(NAME checkGendisk2dDiskHDF5ICTest COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_LIST_DIR}/check_hdf5_particles.py - hdf5-gendisk2d-disk.bods.h5 --count 8 + hdf5-gendisk2d-disk.bods.h5 --count 10000 --mass 0.05 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) set_tests_properties(checkGendisk2dDiskHDF5ICTest PROPERTIES DEPENDS makeGendisk2dHDF5ICTest) @@ -302,7 +302,7 @@ if(ENABLE_NBODY) REQUIRED_FILES "hdf5-gendisk2d-halo.bods.h5;hdf5-gendisk2d-disk.bods.h5") set_tests_properties(makeGendisk2dHDF5ICTest checkGendisk2dHaloHDF5ICTest checkGendisk2dDiskHDF5ICTest removeGendisk2dHDF5ICFiles - PROPERTIES LABELS "quick") + PROPERTIES LABELS "long") endif() # Set labels for pyEXP tests From 97231fca35236e3334a916011759dfb3e0b4a7c4 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 15 Aug 2026 12:07:05 -0400 Subject: [PATCH 21/28] gendisk2d needs a smaller default RMIN to match the table --- tests/CMakeLists.txt | 2 +- utils/ICs/ParticleHDF5.H | 73 +++++++++++++++++++++++++++++----------- 2 files changed, 54 insertions(+), 21 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8e3ac5e91..a36a1a8d6 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -273,7 +273,7 @@ if(ENABLE_NBODY) COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gendisk2d --hdf5 --nhalo 10000 --ndisk 10000 --SEED 17 --NUMX 128 --NUMY 64 --PNUM 0 --TNUM 40 --LMAX 4 --MMAX 2 --NMAXH 4 --NMAXD 4 - # --RMIN 3.333e-05 + --RMIN 3.333e-05 --hbods hdf5-gendisk2d-halo.bods --dbods hdf5-gendisk2d-disk.bods --cachefile hdf5-gendisk2d.cache WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) diff --git a/utils/ICs/ParticleHDF5.H b/utils/ICs/ParticleHDF5.H index adb841ab8..0af0418f5 100644 --- a/utils/ICs/ParticleHDF5.H +++ b/utils/ICs/ParticleHDF5.H @@ -1,21 +1,23 @@ #ifndef EXP_PARTICLE_HDF5_H #define EXP_PARTICLE_HDF5_H -// Shared writer for the simple EXP phase-space HDF5 schema. The -// schema unchanged from the Component update and is the same one tested -// by the hdf5bods and slabics implementations. +// Help class for writing the EXP phase-space HDF5 schema from a +// phase-space array. The schema unchanged from the Component update +// and is the same one tested by the hdf5bods and slabics +// implementations. #include -#include -#include #include +#include #include #include +#include #include #include #include #include + #include namespace EXP @@ -24,7 +26,8 @@ namespace EXP namespace ParticleHDF5 { - //! Collects the mass, position, and velocity of a particle container into a vector of 7-element arrays. + //! Collects the mass, position, and velocity of from an EXP + //! Particle container into a vector of 7-element arrays. template inline std::vector> records(const Particles& particles) { @@ -36,17 +39,24 @@ namespace EXP return result; } - //! Create a HighFive::DataSetCreateProps object with chunking and optional filter. + //! Create a HighFive::DataSetCreateProps object with chunking and + //! optional filter. inline HighFive::DataSetCreateProps properties(hsize_t size, unsigned filter) { - if (size == 0) throw std::invalid_argument("cannot write an empty particle file"); + if (size == 0) throw std::invalid_argument("ParticleHDF5: cannot write an empty particle file"); HighFive::DataSetCreateProps props; // chunk size is min(size, 256kB) but at least 1kB props.add(HighFive::Chunking({std::min(size, std::min(262144, std::max(1024, size)))})); if (filter != 4 && filter != 32001) props.add(HighFive::Shuffle()); - // The filter IDs are from HDF5 1.12.0 and later, see https://support.hdfgroup.org/HDF5/doc/Advanced/Filters.html + // We could use HighFive::Deflate(level) for GZIP, but we want + // to support other filters as well, so we use the HDF5 C API + // directly. + + // The filter IDs are from HDF5 1.12.0 and later, see + // https://support.hdfgroup.org/HDF5/doc/Advanced/Filters.html + // // Filter types: // 1: GZIP (DEFLATE) // 2: SZIP @@ -61,14 +71,37 @@ namespace EXP hid_t plist = props.getId(); switch (filter) { case 1: - case 307: { const unsigned level = 5; H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 1, &level); break; } - case 2: { const unsigned options[] = {141, 16}; H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 2, options); break; } - case 3: return props; // shuffle only - case 4: H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 0, nullptr); return props; - case 32004: { const unsigned options[] = {0}; H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 1, options); break; } - case 32001: { const unsigned options[] = {0, 0, 0, 0, 5, 1, 1}; H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 7, options); break; } + case 307: + { + const unsigned level = 5; + H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 1, &level); + break; + } + case 2: + { + const unsigned options[] = {141, 16}; + H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 2, options); + break; + } + case 3: + return props; // shuffle only, no compression + case 4: + H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 0, nullptr); + return props; + case 32004: + { + const unsigned options[] = {0}; + H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 1, options); + break; + } + case 32001: + { + const unsigned options[] = {0, 0, 0, 0, 5, 1, 1}; + H5Pset_filter(plist, filter, H5Z_FLAG_OPTIONAL, 7, options); + break; + } default: - throw std::invalid_argument("unsupported HDF5 filter ID"); + throw std::invalid_argument("ParticleHDF5: unsupported HDF5 filter ID"); } return props; } @@ -89,7 +122,7 @@ namespace EXP if (count == 0) return; if (x.size() != count || y.size() != count || z.size() != count || u.size() != count || v.size() != count || w.size() != count) - throw std::invalid_argument("inconsistent phase-space vector sizes"); + throw std::invalid_argument("ParticleHDF5: inconsistent phase-space vector sizes"); HighFive::File file(filename, HighFive::File::ReadWrite | HighFive::File::Create | HighFive::File::Truncate); @@ -127,7 +160,7 @@ namespace EXP u.push_back(p[4]); v.push_back(p[5]); w.push_back(p[6]); } if (num_aux_ints < 0 || num_aux_floats < 0) - throw std::invalid_argument("auxiliary field counts must be non-negative"); + throw std::invalid_argument("ParticleHDF5: auxiliary field counts must be non-negative"); std::vector> auxi(num_aux_ints, std::vector(particles.size(), 0)); std::vector> auxf(num_aux_floats, std::vector(particles.size(), 0)); write(filename, m, x, y, z, u, v, w, auxi, auxf, std::nullopt, filter); @@ -142,9 +175,9 @@ namespace EXP { if (particles.empty()) return; if (particles.size() % 7 != 0) - throw std::invalid_argument("interleaved phase-space buffer must have a multiple-of-7 length"); + throw std::invalid_argument("ParticleHDF5: interleaved phase-space buffer must have a multiple-of-7 length"); if (num_aux_ints < 0 || num_aux_floats < 0) - throw std::invalid_argument("auxiliary field counts must be non-negative"); + throw std::invalid_argument("ParticleHDF5: auxiliary field counts must be non-negative"); const auto count = particles.size() / 7; const int n = static_cast(count); From 3e8d002b58cb0df8030cba49fc37bd5ef0e95818 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 16 Aug 2026 09:27:49 -0400 Subject: [PATCH 22/28] Apply suggestions from code review Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/ICs/genslab.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/ICs/genslab.cc b/utils/ICs/genslab.cc index f02ea0e33..6ae9838f8 100644 --- a/utils/ICs/genslab.cc +++ b/utils/ICs/genslab.cc @@ -272,7 +272,7 @@ main(int argc, char **argv) // necessary, but we can use OpenMP to parallelize the particle // generation loop. The reduction clause is used to accumulate // KE and VC across threads. -#pragma omp parallel for schedule(dynamic, 256) reduction(+:KE, VC) num_threads(omp_get_max_threads()) + // Generate serially so the seeded RNG state is not shared across threads. for (int n=0; n Date: Sun, 16 Aug 2026 11:59:24 -0400 Subject: [PATCH 23/28] Added compression to save some disk space --- exputil/EmpCylSL.cc | 62 +++++++++++++++++++++---- include/EmpCylSL.H | 6 +++ utils/ICs/initial.cc | 107 ++++++++++++++++++++++++++++++++++++------- 3 files changed, 149 insertions(+), 26 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 5030c42df..723f5213f 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -68,6 +68,8 @@ double EmpCylSL::RMAX = 20.0; double EmpCylSL::HFAC = 0.2; double EmpCylSL::PPOW = 4.0; bool EmpCylSL::NewCoefs = true; +int EmpCylSL::H5compress = 5; +bool EmpCylSL::H5shuffle = false; EmpCylSL::EmpModel EmpCylSL::mtype = EmpCylSL::EmpModel::Exponential; @@ -2403,6 +2405,22 @@ void EmpCylSL::generate_eof(int numr, int nump, int numt, Timer timer; if (VFLAG & 16) timer.start(); + // Sanity check on quadrature parameters + // + numr = std::max(1, numr); + nump = std::max(1, nump); + numt = std::max(1, numt); + + if (myid==0 and numr < 10) { + std::cerr << "EmpCylSL: Warning, numr=" << numr + << " is very small for quadrature" << std::endl; + } + + if (myid==0 and numt < 10) { + std::cerr << "EmpCylSL: Warning, numt=" << numt + << " is very small for quadrature" << std::endl; + } + // Create spherical orthogonal basis if necessary // if (not ortho) { @@ -2667,7 +2685,7 @@ void EmpCylSL::generate_eof(int numr, int nump, int numt, } // *** r quadrature loop - if (VFLAG & 8) { + if (VFLAG & 16) { auto t = timer.stop(); if (myid==0) { std::cout << std::endl @@ -2700,7 +2718,7 @@ void EmpCylSL::generate_eof(int numr, int nump, int numt, // make_eof(); - if (VFLAG & 8) { + if (VFLAG & 16) { cout << "Process " << setw(4) << myid << ": completed basis in " << timer.stop() << " seconds" << endl; @@ -7463,6 +7481,30 @@ void EmpCylSL::WriteH5Cache() file.createAttribute ("eof_numt", HighFive::DataSpace::From(eof_vars.numt)).write(eof_vars.numt); + // To write an Eigen::MatrixXd with compression in HighFive, we we + // create the dataset with explicit chunk sizes and a deflation + // (gzip) filter before writing, because shortcut functions do not + // pass compression parameters. + + auto dcpl = HighFive::DataSetCreateProps{}; + + // Define dimensions + std::vector dims + {static_cast(NUMX+1), static_cast(NUMY+1)}; + + // Dataspace size definition + HighFive::DataSpace ds(dims); + + // Define chunk size (e.g., matching matrix size or split into blocks) + std::vector chunk_dims = + {static_cast(NUMX+1), static_cast(NUMY+1)}; + + if (H5compress>0) { + dcpl.add(HighFive::Chunking(chunk_dims)); + if (H5shuffle) dcpl.add(HighFive::Shuffle()); + dcpl.add(HighFive::Deflate(H5compress)); + } + // Cosine functions auto cosine = file.createGroup("Cosine"); @@ -7480,10 +7522,10 @@ void EmpCylSL::WriteH5Cache() sout << n; auto order = harmonic.createGroup(sout.str()); - order.createDataSet("potC", potC [m][n]); - order.createDataSet("rforceC", rforceC[m][n]); - order.createDataSet("zforceC", zforceC[m][n]); - order.createDataSet("densC", densC [m][n]); + order.createDataSet("potC", ds, dcpl).write(potC [m][n]); + order.createDataSet("rforceC", ds, dcpl).write(rforceC[m][n]); + order.createDataSet("zforceC", ds, dcpl).write(zforceC[m][n]); + order.createDataSet("densC", ds, dcpl).write(densC [m][n]); } } @@ -7504,10 +7546,10 @@ void EmpCylSL::WriteH5Cache() sout << n; auto order = harmonic.createGroup(sout.str()); - order.createDataSet("potS", potS [m][n]); - order.createDataSet("rforceS", rforceS[m][n]); - order.createDataSet("zforceS", zforceS[m][n]); - order.createDataSet("densS", densS [m][n]); + order.createDataSet("potS", ds, dcpl).write(potS [m][n]); + order.createDataSet("rforceS", ds, dcpl).write(rforceS[m][n]); + order.createDataSet("zforceS", ds, dcpl).write(zforceS[m][n]); + order.createDataSet("densS", ds, dcpl).write(densS [m][n]); } } diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index 3b92175fc..a66ec3b6b 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -448,6 +448,12 @@ public: //! No extrapolating beyond grid (default: false) static bool enforce_limits; + //! H5compression level (default: 5) + static int H5compress; + + //! H5shuffle (default: true) + static bool H5shuffle; + /** @brief Density model type Available options: diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index 8e480078a..007139cd2 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -80,8 +80,14 @@ #include #include +#include +#include +#include +#include + #include +#include "quickdigest5.hpp" #include "config_exp.h" #ifdef HAVE_OMP_H #include @@ -401,11 +407,11 @@ main(int ac, char **av) double disk_mass, gas_mass, gscal_length, ToomreQ, Temp, Tmin; bool const_height, images, multi, basis, zeropos, zerovel; bool report, ignore, evolved, diskmodel; - int nhalo, ndisk, ngas, ngparam; + int nhalo, ndisk, ngas, ngparam, H5compress; std::string hbods, dbods, gbods, outtag, runtag, centerfile, halofile1, halofile2; std::string cachefile, config, gentype, dtype, dmodel, mtype, ctype; unsigned hdf5_filter = 1; - bool hdf5_output = false, hdf5_double = false; + bool hdf5_output = false, hdf5_double = false, H5shuffle = false; const std::string mesg("Generates a Monte Carlo realization of a halo with an\n embedded disk using Jeans' equations\n"); @@ -620,9 +626,15 @@ main(int ac, char **av) cxxopts::value(runtag)) ("threads", "Number of threads to run", cxxopts::value(nthrds)->default_value("1")) + ("H5compress", "HDF5 compression level (0-9)", + cxxopts::value(H5compress)->default_value("0")) + ("H5shuffle", "HDF5 shuffle filter (true/false)", + cxxopts::value(H5shuffle)->default_value("false")) ("allow", "Allow multimass algorithm to generature negative masses for testing") ("nomono", "Allow non-monotonic mass interpolation") ("diskmodel", "Table describing the model for the disk plane") + ("ortho", "Check the orthogonality of the basis functions") + ("spline", "Use spline interpolation for the spherical model table"); ("pyname", "Name of module with the user-specified target disk density", cxxopts::value(pyname)) ; @@ -917,6 +929,8 @@ main(int ac, char **av) EmpCylSL::CMAPZ = CMAPZ; EmpCylSL::VFLAG = VFLAG; EmpCylSL::logarithmic = LOGR; + EmpCylSL::H5compress = H5compress; + EmpCylSL::H5shuffle = H5shuffle; // Create expansion only if needed . . . std::shared_ptr expandd; @@ -927,19 +941,19 @@ main(int ac, char **av) expandd = std::make_shared(NMAXFID, LMAXFID, MMAX, NMAXD, ASCALE, HSCALE, NODD, cachefile); #ifdef DEBUG - std::cout << "Process " << myid << ": " - << " rmin=" << EmpCylSL::RMIN - << " rmax=" << EmpCylSL::RMAX - << " a=" << ASCALE - << " h=" << HSCALE - << " nmaxfid=" << NMAXFID - << " lmaxfid=" << LMAXFID - << " mmax=" << MMAX - << " nmax=" << NMAXD - << " nodd=" << NODD - << std::endl << std::flush; + std::cout << "Process " << myid << ": " + << " rmin=" << EmpCylSL::RMIN + << " rmax=" << EmpCylSL::RMAX + << " a=" << ASCALE + << " h=" << HSCALE + << " nmaxfid=" << NMAXFID + << " lmaxfid=" << LMAXFID + << " mmax=" << MMAX + << " nmax=" << NMAXD + << " nodd=" << NODD + << std::endl << std::flush; #endif - + // Try to read existing cache to get EOF // if (not ignore) { @@ -956,7 +970,7 @@ main(int ac, char **av) return 0; } } - + // Use these user models to deproject for the EOF spherical basis // if (vm.count("deproject")) { @@ -964,7 +978,7 @@ main(int ac, char **av) // height relative to the length // double H = scale_height/scale_length; - + // The model instance (you can add others in DiskModels.H). // It's MN or Exponential if not MN. // @@ -992,6 +1006,67 @@ main(int ac, char **av) if (expcond and not save_eof) { expandd->generate_eof(RNUM, PNUM, TNUM, dcond); save_eof = true; + + if (myid==0) { + + // Orthogonalith check for the cylindrical basis functions + // + if (vm.count("ortho")) { + auto oc = expandd->orthoCheck(); + for (int M=0; M to cache file <" << cachefile << ">" << std::endl; + + // Write the md5sum for the Python module + if (DTYPE == DiskType::python) { + try { + std::vector pyinfo = + {pyname, QuickDigest5::fileToHash(pyname + ".py")}; + + file.createAttribute("pythonDiskType", pyinfo); + + std::cout << "---- Cylindrical: writing pythonDiskType <" << pyname + ".py" + << "> to cache file <" << cachefile << ">" << std::endl; + + + } catch (const std::runtime_error& e) { + if (myid==0) { + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() + << ", can not write the pyname and md5 hash to HDF5" + << std::endl; + } + } + } + } catch (const HighFive::Exception& err) { + if (myid==0) { + std::cerr << err.what() << std::endl; + std::cerr << "Error writing metadata to cache file <" << cachefile + << std::endl; + } + } + // Errors will prevent metadata from being written to the cache + } + // Only the root process should be updating the cache } // Basis orthgonality check From 43bbc9321f738cb9cd68f8421bb80c279b359d3c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 16 Aug 2026 12:08:43 -0400 Subject: [PATCH 24/28] Change hsize_t to size_t to try to keep Clang happy --- exputil/EmpCylSL.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 723f5213f..553140fbb 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -7489,15 +7489,15 @@ void EmpCylSL::WriteH5Cache() auto dcpl = HighFive::DataSetCreateProps{}; // Define dimensions - std::vector dims - {static_cast(NUMX+1), static_cast(NUMY+1)}; + std::vector dims + {static_cast(NUMX+1), static_cast(NUMY+1)}; // Dataspace size definition HighFive::DataSpace ds(dims); // Define chunk size (e.g., matching matrix size or split into blocks) - std::vector chunk_dims = - {static_cast(NUMX+1), static_cast(NUMY+1)}; + std::vector chunk_dims = + {static_cast(NUMX+1), static_cast(NUMY+1)}; if (H5compress>0) { dcpl.add(HighFive::Chunking(chunk_dims)); From 581199f54c5c6e18bc914b1e6e2d1738caf1ab93 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 16 Aug 2026 12:16:53 -0400 Subject: [PATCH 25/28] Another try at type matching --- exputil/EmpCylSL.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 553140fbb..373b6883c 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -7489,15 +7489,15 @@ void EmpCylSL::WriteH5Cache() auto dcpl = HighFive::DataSetCreateProps{}; // Define dimensions - std::vector dims - {static_cast(NUMX+1), static_cast(NUMY+1)}; + std::vector dims {static_cast(NUMX+1), + static_cast(NUMY+1)}; // Dataspace size definition HighFive::DataSpace ds(dims); // Define chunk size (e.g., matching matrix size or split into blocks) - std::vector chunk_dims = - {static_cast(NUMX+1), static_cast(NUMY+1)}; + std::vector chunk_dims {static_cast(NUMX+1), + static_cast(NUMY+1)}; if (H5compress>0) { dcpl.add(HighFive::Chunking(chunk_dims)); From a51ac72541ad0485571a65193ebc375f59249732 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 16 Aug 2026 12:18:37 -0400 Subject: [PATCH 26/28] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/ICs/initial.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index 007139cd2..2568e0627 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -634,7 +634,7 @@ main(int ac, char **av) ("nomono", "Allow non-monotonic mass interpolation") ("diskmodel", "Table describing the model for the disk plane") ("ortho", "Check the orthogonality of the basis functions") - ("spline", "Use spline interpolation for the spherical model table"); + ("spline", "Use spline interpolation for the spherical model table") ("pyname", "Name of module with the user-specified target disk density", cxxopts::value(pyname)) ; From f7d63fdd5bbd00dc2725e8341d9118bcc0754ab5 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 16 Aug 2026 13:30:00 -0400 Subject: [PATCH 27/28] Missing commit of pyEXP updates for Cylindrical, oops --- expui/BiorthBasis.H | 3 +++ expui/BiorthBasis.cc | 13 ++++++++++--- include/EmpCylSL.H | 6 ++++-- 3 files changed, 17 insertions(+), 5 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 6b8bc09eb..a9b2dfe62 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -985,6 +985,9 @@ namespace BasisClasses double rcylmin, rcylmax, acyl, hcyl, bias; bool expcond, logarithmic, density, EVEN_M, sech2 = true; + int H5compress = 0; + bool H5shuffle = false; + std::vector potd, dpot, dpt2, dend; std::vector legs, dlegs, d2legs; diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 0367512c5..fd9e0cc1b 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1279,7 +1279,9 @@ namespace BasisClasses "pyproj", "nint", "totalCovar", - "fullCovar" + "fullCovar", + "H5compress", + "H5shuffle" }; Cylindrical::Cylindrical(const YAML::Node& CONF) : @@ -1529,10 +1531,13 @@ namespace BasisClasses if (conf["vflag" ]) vflag = conf["vflag" ].as(); if (conf["pyname" ]) pyname = conf["pyname" ].as(); if (conf["pyproj" ]) pyproj = conf["pyproj" ].as(); - if (conf["pcavar"] ) pcavar = conf["pcavar" ].as(); - if (conf["subsamp"] ) sampT = conf["subsamp" ].as(); + if (conf["pcavar" ]) pcavar = conf["pcavar" ].as(); + if (conf["subsamp" ]) sampT = conf["subsamp" ].as(); if (conf["samplesz" ]) sampT = conf["samplesz" ].as(); + if (conf["H5compress"]) H5compress = conf["H5compress"].as(); + if (conf["H5shuffle" ]) H5shuffle = conf["H5shuffle" ].as(); + // Sanity sampT = std::max(1, sampT); @@ -1604,6 +1609,8 @@ namespace BasisClasses EmpCylSL::CMAPZ = cmapZ; EmpCylSL::logarithmic = logarithmic; EmpCylSL::VFLAG = vflag; + EmpCylSL::H5compress = H5compress; + EmpCylSL::H5shuffle = H5shuffle; // Convert dmodel string to lower case (deprojection model for EOF // basis construction) diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index a66ec3b6b..df81afcce 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -448,10 +448,12 @@ public: //! No extrapolating beyond grid (default: false) static bool enforce_limits; - //! H5compression level (default: 5) + //! H5compression level (default: 5; good compromise between speed + //! and compression factor) static int H5compress; - //! H5shuffle (default: true) + //! H5shuffle (default: false; tests suggests that shuffling makes + //! cache compression worse) static bool H5shuffle; /** From 4d1654f3c5a59be72ed2488dfbbb78fe401c48d1 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 16 Aug 2026 17:47:12 +0000 Subject: [PATCH 28/28] Change IC generation defaults to HDF5 and H5compress=5 Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- tests/CMakeLists.txt | 8 ++++---- tests/Cube/config.yml | 2 +- tests/Halo/config.yml | 2 +- utils/ICs/ZangICs.cc | 8 +++++--- utils/ICs/cubeICs.cc | 8 +++++--- utils/ICs/gensph.cc | 8 +++++--- utils/ICs/initial.cc | 10 ++++++---- utils/ICs/initial2d.cc | 8 +++++--- 8 files changed, 32 insertions(+), 22 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index a36a1a8d6..bee275c7b 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -118,14 +118,14 @@ if(ENABLE_NBODY) # perhaps there is a better way? add_test(NAME removeTempFiles COMMAND ${CMAKE_COMMAND} -E remove - config.run0.yml current.processor.rates.run0 new.bods + config.run0.yml current.processor.rates.run0 new.bods.h5 OUTLOG.run0 run0.levels SLGridSph.cache.run0 test.grid outcoef.halo.run0 SLGridSph.cache.run0 coefcovar.halo.test_covar.h5 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) # Remove the temporary files set_tests_properties(removeTempFiles PROPERTIES DEPENDS expNbodyCheck2TW - REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods;run0.levels;SLGridSph.cache.run0;test.grid;" + REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods.h5;run0.levels;SLGridSph.cache.run0;test.grid;" ) # Makes some cube ICs using utils/ICs/cubeics @@ -152,13 +152,13 @@ if(ENABLE_NBODY) # A separate test to remove the generated files if they all exist add_test(NAME removeCubeFiles COMMAND ${CMAKE_COMMAND} -E remove - config.runS.yml current.processor.rates.runS cube.bods + config.runS.yml current.processor.rates.runS cube.bods.h5 OUTLOG.runS runS.levels WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Cube) # Remove the temporary files set_tests_properties(removeCubeFiles PROPERTIES DEPENDS expCubeCheckPos - REQUIRED_FILES "config.runS.yml;current.processor.rates.runS;cube.bods;OUTLOG.runS;runS.levels;") + REQUIRED_FILES "config.runS.yml;current.processor.rates.runS;cube.bods.h5;OUTLOG.runS;runS.levels;") # These are basic sanity checks for writing and reading EXP HDF5 IC files if(Python3_Interpreter_FOUND AND EXP_H5PY_STATUS EQUAL 0 AND ENABLE_UTILS) diff --git a/tests/Cube/config.yml b/tests/Cube/config.yml index 148e3a228..4872da17d 100644 --- a/tests/Cube/config.yml +++ b/tests/Cube/config.yml @@ -30,7 +30,7 @@ Global: Components: - name : cube parameters : {nlevel: 1, indexing: true} - bodyfile : cube.bods + bodyfile : cube.bods.h5 force : id : cube parameters : diff --git a/tests/Halo/config.yml b/tests/Halo/config.yml index 596d51f36..b236c74fc 100644 --- a/tests/Halo/config.yml +++ b/tests/Halo/config.yml @@ -28,7 +28,7 @@ Global: Components: - name : halo parameters : {nlevel: 1, indexing: true} - bodyfile : new.bods + bodyfile : new.bods.h5 force : id : sphereSL parameters : diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index 6b8d023ff..ef59f0fc0 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -37,7 +37,7 @@ main(int ac, char **av) unsigned seed; // Will be inialized by /dev/random if // not set on the command line unsigned hdf5_filter = 1; // HDF5 filter type (1 = GZIP, 2 = SZIP, 3 = LZF) - bool hdf5_output = false; // Write HDF5 output instead of ASCII + bool hdf5_output = true; // Write HDF5 output by default bool hdf5_double = false; // Write HDF5 output in double precision @@ -48,7 +48,8 @@ main(int ac, char **av) ("V,nozerovel", "Do not zero the mean velocity") ("P,nozeropos", "Do not zero the center of mass") ("d,debug", "Print debug grid") - ("5,hdf5", "Write HDF5 phase-space output instead of ASCII") + ("5,hdf5", "Write HDF5 phase-space output (default)") + ("A,ascii", "Write old-style ASCII output (default is HDF5)") ("8,double", "Use float64 HDF5 output (default is float32)") ("Z,filter", "HDF5 filter ID (default: 1 = GZIP)", cxxopts::value(hdf5_filter)->default_value("1")) @@ -96,7 +97,8 @@ main(int ac, char **av) return 1; } - hdf5_output = vm.count("hdf5") > 0; + if (vm.count("hdf5")) hdf5_output = true; + if (vm.count("ascii")) hdf5_output = false; hdf5_double = vm.count("double") > 0; // Set from /dev/random if not specified diff --git a/utils/ICs/cubeICs.cc b/utils/ICs/cubeICs.cc index 787fee237..24122f8ab 100644 --- a/utils/ICs/cubeICs.cc +++ b/utils/ICs/cubeICs.cc @@ -30,7 +30,7 @@ main(int ac, char **av) unsigned seed; // Will be inialized by /dev/random if // not set on the command line unsigned hdf5_filter = 1; - bool hdf5_output = false, hdf5_double = false; + bool hdf5_output = true, hdf5_double = false; // Default values for the velocity dispersion and bulk velocity // @@ -58,7 +58,8 @@ main(int ac, char **av) cxxopts::value(seed)) ("o,file", "Output body file", cxxopts::value(bodyfile)->default_value("cube.bods")) - ("5,hdf5", "Write HDF5 phase-space output instead of ASCII") + ("5,hdf5", "Write HDF5 phase-space output (default)") + ("A,ascii", "Write old-style ASCII output (default is HDF5)") ("8,double", "Use float64 HDF5 output (default is float32)") ("f,filter", "HDF5 filter ID (default: 1 = GZIP)", cxxopts::value(hdf5_filter)->default_value("1")) @@ -83,7 +84,8 @@ main(int ac, char **av) std::cout << options.help() << std::endl << std::endl; return 1; } - hdf5_output = vm.count("hdf5") > 0; + if (vm.count("hdf5")) hdf5_output = true; + if (vm.count("ascii")) hdf5_output = false; hdf5_double = vm.count("double") > 0; // Set from /dev/random if not specified diff --git a/utils/ICs/gensph.cc b/utils/ICs/gensph.cc index 37f373713..5f3c53211 100644 --- a/utils/ICs/gensph.cc +++ b/utils/ICs/gensph.cc @@ -77,7 +77,7 @@ main(int argc, char **argv) double Emin0, Emax0, Kmin0, Kmax0, RBAR, MBAR, BRATIO, CRATIO, SMOOTH; bool LOGR, ELIMIT, VERBOSE, GRIDPOT, MODELS, EBAR, zeropos, zerovel; bool VTEST; - bool hdf5_output = false, hdf5_double = false; + bool hdf5_output = true, hdf5_double = false; unsigned hdf5_filter = 1; std::string INFILE, MMFILE, OUTFILE, OUTPS, config; @@ -111,7 +111,8 @@ main(int argc, char **argv) cxxopts::value(MMFILE)) ("o,PSFILE", "Phase-space output file", cxxopts::value(OUTPS)->default_value("new.bods")) - ("5,hdf5", "Write HDF5 phase-space output instead of ASCII") + ("5,hdf5", "Write HDF5 phase-space output (default)") + ("A,ascii", "Write old-style ASCII output (default is HDF5)") ("8,double", "Use float64 HDF5 output (default is float32)") ("f,filter", "HDF5 filter ID (default: 1 = GZIP)", cxxopts::value(hdf5_filter)->default_value("1")) @@ -266,7 +267,8 @@ main(int argc, char **argv) return 0; } } - hdf5_output = vm.count("hdf5") > 0; + if (vm.count("hdf5")) hdf5_output = true; + if (vm.count("ascii")) hdf5_output = false; hdf5_double = vm.count("double") > 0; if (vm.count("verbose")) VERBOSE = true; diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index 2568e0627..28618d423 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -411,7 +411,7 @@ main(int ac, char **av) std::string hbods, dbods, gbods, outtag, runtag, centerfile, halofile1, halofile2; std::string cachefile, config, gentype, dtype, dmodel, mtype, ctype; unsigned hdf5_filter = 1; - bool hdf5_output = false, hdf5_double = false, H5shuffle = false; + bool hdf5_output = true, hdf5_double = false, H5shuffle = false; const std::string mesg("Generates a Monte Carlo realization of a halo with an\n embedded disk using Jeans' equations\n"); @@ -436,7 +436,8 @@ main(int ac, char **av) cxxopts::value(gbods)->default_value("gas.bods")) ("dbods", "The output bodyfile for the stellar disc", cxxopts::value(dbods)->default_value("disk.bods")) - ("5,hdf5", "Write HDF5 phase-space output instead of ASCII") + ("5,hdf5", "Write HDF5 phase-space output (default)") + ("A,ascii", "Write old-style ASCII output (default is HDF5)") ("8,double", "Use float64 HDF5 output (default is float32)") ("f,filter", "HDF5 filter ID (default: 1 = GZIP)", cxxopts::value(hdf5_filter)->default_value("1")) @@ -627,7 +628,7 @@ main(int ac, char **av) ("threads", "Number of threads to run", cxxopts::value(nthrds)->default_value("1")) ("H5compress", "HDF5 compression level (0-9)", - cxxopts::value(H5compress)->default_value("0")) + cxxopts::value(H5compress)->default_value("5")) ("H5shuffle", "HDF5 shuffle filter (true/false)", cxxopts::value(H5shuffle)->default_value("false")) ("allow", "Allow multimass algorithm to generature negative masses for testing") @@ -689,7 +690,8 @@ main(int ac, char **av) return 0; } } - hdf5_output = vm.count("hdf5") > 0; + if (vm.count("hdf5")) hdf5_output = true; + if (vm.count("ascii")) hdf5_output = false; hdf5_double = vm.count("double") > 0; if (vm.count("spline")) { diff --git a/utils/ICs/initial2d.cc b/utils/ICs/initial2d.cc index 588d2f7ab..c2b688e7b 100644 --- a/utils/ICs/initial2d.cc +++ b/utils/ICs/initial2d.cc @@ -211,7 +211,7 @@ main(int ac, char **av) std::string cachefile, config, gentype, dtype, dmodel, mtype, ctype; std::string diskconf; unsigned hdf5_filter = 1; - bool hdf5_output = false, hdf5_double = false; + bool hdf5_output = true, hdf5_double = false; const std::string mesg("Generates a Monte Carlo realization of a halo with an\n embedded disk using Jeans' equations\n"); @@ -226,7 +226,8 @@ main(int ac, char **av) cxxopts::value(hbods)->default_value("halo.bods")) ("dbods", "The output bodyfile for the stellar disc", cxxopts::value(dbods)->default_value("disk.bods")) - ("5,hdf5", "Write HDF5 phase-space output instead of ASCII") + ("5,hdf5", "Write HDF5 phase-space output (default)") + ("A,ascii", "Write old-style ASCII output (default is HDF5)") ("8,double", "Use float64 HDF5 output (default is float32)") ("f,filter", "HDF5 filter ID (default: 1 = GZIP)", cxxopts::value(hdf5_filter)->default_value("1")) @@ -436,7 +437,8 @@ main(int ac, char **av) return 0; } } - hdf5_output = vm.count("hdf5") > 0; + if (vm.count("hdf5")) hdf5_output = true; + if (vm.count("ascii")) hdf5_output = false; hdf5_double = vm.count("double") > 0; if (vm.count("spline")) {