From 3f2432512f17b41acf9a7c753b65b36194004afa Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Apr 2026 10:53:46 -0400 Subject: [PATCH 01/51] Updated isothermal slab scale height with the factor of 2 for consistency with the recent Cylinder update --- exputil/SLGridMP2.cc | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index e13c8cd46..620637cd3 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1844,24 +1844,28 @@ static double poffset=0.0; class IsothermalSlab : public SlabModel { +private: + + std::string psa = "IsothermalSlab NOW uses the traditional density profile proportional to sech^2(z/2H). If you are using the old profile proportional to sech(z/H), please update your model file to use the new profile and set the scale height H to be half of the old value. This will ensure that your model has the same density profile and potential as before, but with a more standard functional form. If you have any questions or concerns about this change, please contact the developers."; + public: - IsothermalSlab() { id = "iso"; } + IsothermalSlab() { id = "iso"; if (myid==0) std::cout << "SLGridSlab: " << psa << std::endl; } double pot(double z) { - return 2.0*M_PI*SLGridSlab::H*log(cosh(z/SLGridSlab::H)) - poffset; + return 4.0*M_PI*SLGridSlab::H*log(cosh(0.5*z/SLGridSlab::H)) - poffset; } double dpot(double z) { - return 2.0*M_PI*tanh(z/SLGridSlab::H); + return 2.0*M_PI*tanh(0.5*z/SLGridSlab::H); } double dens(double z) { - double tmp = 1.0/cosh(z/SLGridSlab::H); - return 4.0*M_PI * 0.5/SLGridSlab::H * tmp*tmp; + double tmp = 1.0/cosh(0.5*z/SLGridSlab::H); + return 0.25/SLGridSlab::H * tmp*tmp; } }; From 85864fae4be67c7c81e5d73d434f4eb9928846ac Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Apr 2026 10:54:36 -0400 Subject: [PATCH 02/51] Updated SlabSL method to have a 'self_consistent' flag --- expui/BiorthBasis.cc | 3 ++- src/SlabSL.H | 44 +++++++++++++++++++++++++++++++++++++++++--- src/SlabSL.cc | 31 +++++++++++++++++++++++++------ src/cudaSlabSL.cu | 6 ++++-- 4 files changed, 72 insertions(+), 12 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index b7723992f..3c79f8209 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3420,7 +3420,8 @@ namespace BasisClasses "knots", "verbose", "check", - "method" + "method", + "self_consistent" }; Slab::Slab(const YAML::Node& CONF) : BiorthBasis(CONF, "slab") diff --git a/src/SlabSL.H b/src/SlabSL.H index 693ded977..d9f5fdd51 100644 --- a/src/SlabSL.H +++ b/src/SlabSL.H @@ -17,9 +17,38 @@ #include #endif -/*! This routine computes the potential, acceleration and density - using expansion periodic in X & Y and outgoing vacuum boundary - condtions in Z */ +/** @class SlabSL + @brief This routine computes the potential, acceleration and density + using expansion periodic in X & Y and outgoing vacuum boundary + condtions in Z + + @details **YAML configuration** + + @param nmaxx is the maximum order of the expansion in x (default 6) + + @param nmaxy is the maximum order of the expansion in y (default 6) + + @param nmaxz is the maximum order of the expansion in z (default 6) + + @param nminx is the minimum order of the expansion in x (default 0) + + @param nminy is the minimum order of the expansion in y (default 0) + + @param hslab is the scale height of the slab (default 0.2) + + @param zmax is the maximum z for the slab (default 10.0) + + @param ngrid is the number of grid points in z for the + Sturm-Liouville solver (default 1000) + + @param type is the type of slab to solve for (default + "isothermal", must be "isothermal", "parabolic", or "constant") + + @param self_consistent set to true allows the particles to evolve + under the time-dependent basis expansion. For a basis fixed in time to + the initial time: set to false. + +*/ class SlabSL : public PotAccel { @@ -47,6 +76,9 @@ private: //! Current coefficient tensor std::vector expccof, expccofP; + //! Coefficient tensor for frozen potential (if self_consistent=false) + coefType expcofF; + int nminx, nminy; int nmaxx, nmaxy, nmaxz; double zmax, hslab; @@ -120,6 +152,12 @@ private: #endif + //! Flag self_consitency + bool self_consistent = true; + + //! Flag whether coefficients have been initialized for the first time + bool firstime_coef = true; + //! Default number of grid points for SLGridSlab int ngrid = 1000; diff --git a/src/SlabSL.cc b/src/SlabSL.cc index 1de35675a..a6c6fd7d8 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -17,7 +17,8 @@ SlabSL::valid_keys = { "hslab", "zmax", "ngrid", - "type" + "type", + "self_consistent" }; //@{ @@ -48,6 +49,8 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) int nnmax = (nmaxx > nmaxy) ? nmaxx : nmaxy; + // Make the Sturm-Liouville grid and basis functions + // grid = std::make_shared(nnmax, nmaxz, ngrid, zmax, type); // Test for basis consistency (will generate an exception if maximum @@ -164,6 +167,11 @@ void SlabSL::initialize() if (conf["hslab"]) hslab = conf["hslab"].as(); if (conf["zmax" ]) zmax = conf["zmax" ].as(); if (conf["type" ]) type = conf["type" ].as(); + + if (conf["self_consistent"]) { + self_consistent = conf["self_consistent"].as(); + } else + self_consistent = true; } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in SlabSL: " @@ -253,6 +261,11 @@ void SlabSL::determine_coefficients(void) compute_multistep_coefficients(); } + // Only used if self_consistent is false to ensure that coefficients + // are only computed once (at the first time step) + // + firstime_coef = false; + if (not self_consistent) expcofF = expccof[0]; } void * SlabSL::determine_coefficients_thread(void * arg) @@ -343,7 +356,7 @@ void SlabSL::get_acceleration_and_potential(Component* C) if (use_external == false) { - if (multistep && initializing) { + if (multistep && (self_consistent || initializing)) { compute_multistep_coefficients(); } @@ -454,11 +467,17 @@ void * SlabSL::determine_acceleration_and_potential_thread(void * arg) for (int iz=0; iz Date: Sun, 5 Apr 2026 12:28:40 -0400 Subject: [PATCH 03/51] Comment updates only --- include/SLGridMP2.H | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index da209403c..caa8976b2 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -292,12 +292,21 @@ private: void compute_table_worker(void); - // Local MPI stuff + //@{ + //! Local MPI stuff void mpi_setup(void); void mpi_unpack_table(void); int mpi_pack_table(TableSlab* table, int kx, int ky); + //@} + + //@{ + //! Cache reading and writing bool ReadH5Cache(void); void WriteH5Cache(void); + //@} + + //! Cache versioning + inline static const std::string Version = "1.0"; int mpi_myid, mpi_numprocs; int mpi_bufsz; From 3d55d7574f263acd07182b18696aec68527025c1 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Apr 2026 12:29:09 -0400 Subject: [PATCH 04/51] Added cache versioning --- exputil/SLGridMP2.cc | 26 +++++++++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 620637cd3..0568c9392 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1846,11 +1846,18 @@ class IsothermalSlab : public SlabModel private: - std::string psa = "IsothermalSlab NOW uses the traditional density profile proportional to sech^2(z/2H). If you are using the old profile proportional to sech(z/H), please update your model file to use the new profile and set the scale height H to be half of the old value. This will ensure that your model has the same density profile and potential as before, but with a more standard functional form. If you have any questions or concerns about this change, please contact the developers."; + std::string psa = + "---- IsothermalSlab NOW uses the traditional density profile proportional to sech^2(z/2H).\n" + "---- If you are using the old profile proportional to sech(z/H), please update your model\n" + "---- to use the new profile and set the scale height H to be half of the old value. This\n" + "---- will ensure that your model has the same density profile and potential as before, but\n" + "---- with a more standard functional form. If you have any questions or concerns about this\n" + "---- change, please contact the developers on GitHUB."; public: - IsothermalSlab() { id = "iso"; if (myid==0) std::cout << "SLGridSlab: " << psa << std::endl; } + IsothermalSlab() { id = "iso"; if (myid==0) std::cout << "---- SLGridSlab: IMPORTANT UPDATE\n" << psa + << std::endl; } double pot(double z) { @@ -1865,7 +1872,7 @@ class IsothermalSlab : public SlabModel double dens(double z) { double tmp = 1.0/cosh(0.5*z/SLGridSlab::H); - return 0.25/SLGridSlab::H * tmp*tmp; + return 4.0*M_PI*0.25/SLGridSlab::H * tmp*tmp; } }; @@ -2259,6 +2266,18 @@ bool SLGridSlab::ReadH5Cache(void) if (not checkStr(geometry, "geometry")) return false; if (not checkStr(forceID, "forceID")) return false; + // Version check + // + if (h5file.hasAttribute("Version")) { + if (not checkStr(Version, "Version")) return false; + } else { + if (myid==0) + std::cout << "---- SLGridSlab::ReadH5Cache: " + << "recomputing cache for HighFive API change" + << std::endl; + return false; + } + // Parameter check // if (not checkStr(type, "type")) return false; @@ -2352,6 +2371,7 @@ void SLGridSlab::WriteH5Cache(void) file.createAttribute("geometry", HighFive::DataSpace::From(geometry)).write(geometry); file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); + file.createAttribute("Version", HighFive::DataSpace::From(Version)).write(Version); // Write parameters // From aeb39a2b7aadb2a436bb70355f5b44b63edd606b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Apr 2026 13:20:05 -0400 Subject: [PATCH 05/51] Added compile-time version string parsing for future development work --- exputil/libvars.cc | 2 ++ include/EXPversion.H | 49 ++++++++++++++++++++++++++++++++++++++++++++ include/libvars.H | 6 ++++++ 3 files changed, 57 insertions(+) create mode 100644 include/EXPversion.H diff --git a/exputil/libvars.cc b/exputil/libvars.cc index 14ff71690..11ee23a56 100644 --- a/exputil/libvars.cc +++ b/exputil/libvars.cc @@ -3,6 +3,7 @@ standalone utilities */ +#include "config_exp.h" #include "libvars.H" namespace __EXP__ @@ -36,6 +37,7 @@ namespace __EXP__ //! Sanity tolerance for orthogonality double orthoTol = 1.0e-2; + }; diff --git a/include/EXPversion.H b/include/EXPversion.H new file mode 100644 index 000000000..41f7625b6 --- /dev/null +++ b/include/EXPversion.H @@ -0,0 +1,49 @@ +#pragma once + +#include + +namespace __EXP__ +{ + // Runtime parser for "X.Y.Z" format + + /* Example usage: + + #define VERSION_STR "2.4.12" + static constexpr EXPversion current_v = EXP_parse_version(VERSION_STR); + + int main() { + static_assert(current_v.major == 2, "Major version mismatch"); + std::cout << "Parsed: " << current_v.major << "." << current_v.minor << "\n"; + return 0;} + */ + + + struct EXPversion + { + int major, minor, patch; + }; + + // Compile-time parser for "X.Y.Z" format + constexpr EXPversion EXP_parse_version(const char* str) + { + EXPversion v = {0, 0, 0}; + int* target = &v.major; + int current = 0; + + // Scan the version string until we hit a dot, then move to the next target + for (int i = 0; str[i] != '\0'; ++i) { + if (str[i] == '.') { + *target = current; + if (target == &v.major) target = &v.minor; + else if (target == &v.minor) target = &v.patch; + current = 0; + } else { + // Add the current digit to the current version component + current = current * 10 + (str[i] - '0'); + } + } + *target = current; + return v; + } +} + diff --git a/include/libvars.H b/include/libvars.H index 12dfd45d3..5dde89c0c 100644 --- a/include/libvars.H +++ b/include/libvars.H @@ -7,6 +7,9 @@ #include +#include "config_exp.h" // For version string +#include "EXPversion.H" // For compile-time version parsing + namespace __EXP__ { //! @name Theading variables @@ -42,6 +45,9 @@ namespace __EXP__ //! Sanity tolerance for orthogonality extern double orthoTol; + //! Compile-time version parsing + static constexpr EXPversion exp_version = EXP_parse_version(VERSION); + }; #endif // END _LIBVARS_H From c8af8b30727ad69754e44c29e3ed4f1992765ab0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Apr 2026 13:20:36 -0400 Subject: [PATCH 06/51] Use compile-time version info to toggle warning message --- exputil/SLGridMP2.cc | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 0568c9392..c8df89b01 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -21,6 +21,7 @@ #include "SLGridMP2.H" #include "massmodel.H" #include "EXPmath.H" +#include "libvars.H" // For parsed version info #ifdef USE_DMALLOC #include @@ -1856,8 +1857,12 @@ class IsothermalSlab : public SlabModel public: - IsothermalSlab() { id = "iso"; if (myid==0) std::cout << "---- SLGridSlab: IMPORTANT UPDATE\n" << psa - << std::endl; } + IsothermalSlab() { + id = "iso"; + if (myid==0 and __EXP__::exp_version.minor<11) + std::cout << "---- SLGridSlab: IMPORTANT UPDATE\n" << psa + << std::endl; + } double pot(double z) { From a5149ed8d941734f31c0338ffc3e537e88c3397c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Apr 2026 13:44:22 -0400 Subject: [PATCH 07/51] Remove unnecessary preprocessor variable --- utils/ICs/bonnerebert.cc | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/utils/ICs/bonnerebert.cc b/utils/ICs/bonnerebert.cc index a7589b6d0..c20c38ff4 100644 --- a/utils/ICs/bonnerebert.cc +++ b/utils/ICs/bonnerebert.cc @@ -28,8 +28,6 @@ */ -#define VERSION 0.1 - #include #include #include @@ -366,8 +364,7 @@ decode_switches (int argc, char **argv) "R:" /* runit */ "N:" /* number */ "S:" /* seed */ - "h" /* help */ - "V", /* version */ + "h", /* help */ long_options, (int *) 0)) != EOF) { switch (c) @@ -399,10 +396,6 @@ decode_switches (int argc, char **argv) case 'S': /* --seed */ S = atoi(optarg); break; - case 'V': - cout << program_name << " " << VERSION << endl; - exit (0); - case 'h': usage (0); From 5cffd92ff3ec15eb6884b82ca4688aca2169fee4 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Apr 2026 13:44:47 -0400 Subject: [PATCH 08/51] Added test example for internal version parser --- utils/Test/CMakeLists.txt | 3 ++- utils/Test/version_test.cc | 15 +++++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) create mode 100644 utils/Test/version_test.cc diff --git a/utils/Test/CMakeLists.txt b/utils/Test/CMakeLists.txt index 5a8a6fc4b..dd3a4c56c 100644 --- a/utils/Test/CMakeLists.txt +++ b/utils/Test/CMakeLists.txt @@ -1,5 +1,5 @@ -set(bin_PROGRAMS testBarrier expyaml orthoTest testED) +set(bin_PROGRAMS testBarrier expyaml orthoTest testED vtest) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil yaml-cpp ${VTK_LIBRARIES}) @@ -33,6 +33,7 @@ add_executable(testBarrier test_barrier.cc) add_executable(expyaml test_config.cc) add_executable(orthoTest orthoTest.cc Biorth2Ortho.cc) add_executable(testED testED.cc) +add_executable(vtest version_test.cc) foreach(program ${bin_PROGRAMS}) target_link_libraries(${program} ${common_LINKLIB}) diff --git a/utils/Test/version_test.cc b/utils/Test/version_test.cc new file mode 100644 index 000000000..5f8ce527e --- /dev/null +++ b/utils/Test/version_test.cc @@ -0,0 +1,15 @@ +// Example of using the version string and parsed octet of integers from libvars.H + +#include +#include "libvars.H" +using namespace __EXP__; + +int main() +{ + std::cout << "Version string: " << VERSION << '\n'; + std::cout << "Parsed octet of integers:\n"; + std::cout << "-- Major=" << exp_build.major << '\n'; + std::cout << "-- Minor=" << exp_build.minor << '\n'; + std::cout << "-- Patch=" << exp_build.patch << '\n'; + return 0; +} From ca731d4a4cc944a5ebce7a44e5680fa8b3a06331 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Apr 2026 13:45:13 -0400 Subject: [PATCH 09/51] Updated name of version struction to prevent variable name conflicts --- exputil/SLGridMP2.cc | 2 +- include/libvars.H | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index c8df89b01..7cc30bd4c 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1859,7 +1859,7 @@ class IsothermalSlab : public SlabModel IsothermalSlab() { id = "iso"; - if (myid==0 and __EXP__::exp_version.minor<11) + if (myid==0 and exp_build.minor<11) std::cout << "---- SLGridSlab: IMPORTANT UPDATE\n" << psa << std::endl; } diff --git a/include/libvars.H b/include/libvars.H index 5dde89c0c..bbbcd6e11 100644 --- a/include/libvars.H +++ b/include/libvars.H @@ -46,7 +46,7 @@ namespace __EXP__ extern double orthoTol; //! Compile-time version parsing - static constexpr EXPversion exp_version = EXP_parse_version(VERSION); + static constexpr EXPversion exp_build = EXP_parse_version(VERSION); }; From 08fc515a5f43d80adb03cb5c7538348d11673736 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Apr 2026 15:36:23 -0400 Subject: [PATCH 10/51] stdout changes only --- exputil/SLGridMP2.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 7cc30bd4c..4d27115a5 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1852,16 +1852,16 @@ class IsothermalSlab : public SlabModel "---- If you are using the old profile proportional to sech(z/H), please update your model\n" "---- to use the new profile and set the scale height H to be half of the old value. This\n" "---- will ensure that your model has the same density profile and potential as before, but\n" - "---- with a more standard functional form. If you have any questions or concerns about this\n" - "---- change, please contact the developers on GitHUB."; + "---- with a more standard functional form. If you have any questions or concerns about\n" + "---- this change, please contact the developers on GitHUB."; public: IsothermalSlab() { id = "iso"; if (myid==0 and exp_build.minor<11) - std::cout << "---- SLGridSlab: IMPORTANT UPDATE\n" << psa - << std::endl; + std::cout << "---- SLGridSlab: IMPORTANT UPDATE for EXP " + << VERSION << '\n' << psa << std::endl; } double pot(double z) From bd75f99f88247569b002645a7b152148f55d619a Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Apr 2026 16:33:36 -0400 Subject: [PATCH 11/51] Update src/SlabSL.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/SlabSL.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SlabSL.cc b/src/SlabSL.cc index a6c6fd7d8..5e2454af6 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -264,8 +264,8 @@ void SlabSL::determine_coefficients(void) // Only used if self_consistent is false to ensure that coefficients // are only computed once (at the first time step) // + if (not self_consistent and firstime_coef) expcofF = expccof[0]; firstime_coef = false; - if (not self_consistent) expcofF = expccof[0]; } void * SlabSL::determine_coefficients_thread(void * arg) From 69a25c4b6bbb636c5213d52deb1e3ba19e262694 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Apr 2026 16:33:55 -0400 Subject: [PATCH 12/51] Update src/SlabSL.H Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/SlabSL.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SlabSL.H b/src/SlabSL.H index d9f5fdd51..8fce83898 100644 --- a/src/SlabSL.H +++ b/src/SlabSL.H @@ -20,7 +20,7 @@ /** @class SlabSL @brief This routine computes the potential, acceleration and density using expansion periodic in X & Y and outgoing vacuum boundary - condtions in Z + conditions in Z @details **YAML configuration** From f23f427ee36f604a887401dbd8d40b40477b1b9b Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Apr 2026 16:34:09 -0400 Subject: [PATCH 13/51] Update src/SlabSL.H Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/SlabSL.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SlabSL.H b/src/SlabSL.H index 8fce83898..c394a9b84 100644 --- a/src/SlabSL.H +++ b/src/SlabSL.H @@ -152,7 +152,7 @@ private: #endif - //! Flag self_consitency + //! Flag self_consistency bool self_consistent = true; //! Flag whether coefficients have been initialized for the first time From 4d7f43b03af41bd0a8e7972909db2d869b06ae53 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Apr 2026 16:35:06 -0400 Subject: [PATCH 14/51] Update exputil/SLGridMP2.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exputil/SLGridMP2.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 4d27115a5..11082a853 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1859,7 +1859,8 @@ class IsothermalSlab : public SlabModel IsothermalSlab() { id = "iso"; - if (myid==0 and exp_build.minor<11) + if (myid==0 and (exp_build.major < 7 or + (exp_build.major == 7 and exp_build.minor < 11))) std::cout << "---- SLGridSlab: IMPORTANT UPDATE for EXP " << VERSION << '\n' << psa << std::endl; } From c84a843616a4fac7b04bb685648d273a4d835807 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Apr 2026 16:35:39 -0400 Subject: [PATCH 15/51] Update exputil/SLGridMP2.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exputil/SLGridMP2.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 11082a853..cd7ef60e5 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1849,11 +1849,11 @@ class IsothermalSlab : public SlabModel std::string psa = "---- IsothermalSlab NOW uses the traditional density profile proportional to sech^2(z/2H).\n" - "---- If you are using the old profile proportional to sech(z/H), please update your model\n" + "---- If you are using the old profile proportional to sech^2(z/H), please update your model\n" "---- to use the new profile and set the scale height H to be half of the old value. This\n" "---- will ensure that your model has the same density profile and potential as before, but\n" "---- with a more standard functional form. If you have any questions or concerns about\n" - "---- this change, please contact the developers on GitHUB."; + "---- this change, please contact the developers on GitHub."; public: From 7e4597528a2b8f032e6e062915d27128669c8d1a Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Apr 2026 16:36:02 -0400 Subject: [PATCH 16/51] Update include/EXPversion.H Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- include/EXPversion.H | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/include/EXPversion.H b/include/EXPversion.H index 41f7625b6..94edc0b18 100644 --- a/include/EXPversion.H +++ b/include/EXPversion.H @@ -1,10 +1,8 @@ #pragma once -#include - namespace __EXP__ { - // Runtime parser for "X.Y.Z" format + // Compile-time parser for "X.Y.Z" format /* Example usage: From a3aaefe2a8cf9afc8ef4c04ebd329082d6dc4777 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Apr 2026 16:57:32 -0400 Subject: [PATCH 17/51] Update output message for parsed version integers --- utils/Test/version_test.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/Test/version_test.cc b/utils/Test/version_test.cc index 5f8ce527e..0d51ae26c 100644 --- a/utils/Test/version_test.cc +++ b/utils/Test/version_test.cc @@ -7,7 +7,7 @@ using namespace __EXP__; int main() { std::cout << "Version string: " << VERSION << '\n'; - std::cout << "Parsed octet of integers:\n"; + std::cout << "Parsed into integer triplet:\n"; std::cout << "-- Major=" << exp_build.major << '\n'; std::cout << "-- Minor=" << exp_build.minor << '\n'; std::cout << "-- Patch=" << exp_build.patch << '\n'; From b7ed7eb5989e7fe7a2e4b33225a730dea865148b Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 5 Apr 2026 21:13:08 +0000 Subject: [PATCH 18/51] Fix remaining 'octet of integers' wording in version_test.cc comment Agent-Logs-Url: https://github.com/EXP-code/EXP/sessions/a9b4ee92-2327-4580-968b-23dded26f6e2 Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- utils/Test/version_test.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/Test/version_test.cc b/utils/Test/version_test.cc index 0d51ae26c..16f2b585b 100644 --- a/utils/Test/version_test.cc +++ b/utils/Test/version_test.cc @@ -1,4 +1,4 @@ -// Example of using the version string and parsed octet of integers from libvars.H +// Example of using the version string and parsed integer triplet from libvars.H #include #include "libvars.H" From 875571d3c9328fdc0703927e1f0dd377b73c294b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Apr 2026 19:49:35 -0400 Subject: [PATCH 19/51] Comments only [no ci] --- include/EXPversion.H | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/include/EXPversion.H b/include/EXPversion.H index 41f7625b6..c0619ab2f 100644 --- a/include/EXPversion.H +++ b/include/EXPversion.H @@ -23,7 +23,7 @@ namespace __EXP__ int major, minor, patch; }; - // Compile-time parser for "X.Y.Z" format + // Compile-time parser for "X.Y.Z" format using c++-17 features constexpr EXPversion EXP_parse_version(const char* str) { EXPversion v = {0, 0, 0}; @@ -34,15 +34,16 @@ namespace __EXP__ for (int i = 0; str[i] != '\0'; ++i) { if (str[i] == '.') { *target = current; - if (target == &v.major) target = &v.minor; + if (target == &v.major) target = &v.minor; else if (target == &v.minor) target = &v.patch; current = 0; } else { - // Add the current digit to the current version component + // Append the current decimal digit to the current version current = current * 10 + (str[i] - '0'); } } - *target = current; + // Store the int of this version string + *target = current; return v; } } From 1c8a6d416aed2116ea10bdbc20252676f452a80c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 10 Apr 2026 10:16:44 -0400 Subject: [PATCH 20/51] Require SLGridSlab to have a non-default cachename for consistency with SLGridSph and provenance policy --- expui/BiorthBasis.cc | 20 ++++- exputil/SLGridMP2.cc | 9 ++- include/SLGridMP2.H | 6 +- src/SlabSL.H | 3 + src/SlabSL.cc | 10 ++- utils/SL/orthochk.cc | 172 ++++++++----------------------------------- 6 files changed, 69 insertions(+), 151 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 3c79f8209..c1d3b433d 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3421,7 +3421,8 @@ namespace BasisClasses "verbose", "check", "method", - "self_consistent" + "self_consistent", + "cachename" }; Slab::Slab(const YAML::Node& CONF) : BiorthBasis(CONF, "slab") @@ -3482,6 +3483,8 @@ namespace BasisClasses if (conf["knots"]) knots = conf["knots"].as(); if (conf["check"]) check = conf["check"].as(); + + if (conf["cachename"]) cachename = conf["cachename"].as(); } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameter stanza for <" @@ -3494,6 +3497,19 @@ namespace BasisClasses throw std::runtime_error("Slab: error parsing YAML"); } + // Check for non-null cache file name. This must be specified + // to prevent recomputation and unexpected behavior. + // + if (cachename.size() == 0) { + throw std::runtime_error + ("SlabSL requires a specified cachename in your YAML config\n" + "for consistency with previous invocations and existing coefficient\n" + "sets. Please add explicitly add 'cachename: name' to your config\n" + "with new 'name' for creating a basis or an existing 'name' for\n" + "reading a previously generated basis cache\n"); + } + + // Finally, make the basis // SLGridSlab::mpi = 0; @@ -3503,7 +3519,7 @@ namespace BasisClasses int nnmax = (nmaxx > nmaxy) ? nmaxx : nmaxy; - ortho = std::make_shared(nnmax, nmaxz, ngrid, zmax, type); + ortho = std::make_shared(nnmax, nmaxz, ngrid, zmax, cachename, type); // Orthogonality sanity check // diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index cd7ef60e5..a6c20e260 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1973,8 +1973,12 @@ void SLGridSlab::bomb(string oops) // Constructors SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, - const std::string TYPE, bool VERBOSE) + const std::string cachename, const std::string TYPE, + bool VERBOSE) { + if (cachename.size()) slab_cache_name = cachename; + else throw std::runtime_error("SLGridSlab: you must specify a cachename"); + int kx, ky; numk = NUMK; @@ -2210,9 +2214,6 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, } -const string slab_cache_name = ".slgrid_slab_cache"; - - bool SLGridSlab::ReadH5Cache(void) { if (!cache) return false; diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index caa8976b2..5b6bc4525 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -305,6 +305,9 @@ private: void WriteH5Cache(void); //@} + //! Cache file name + std::string slab_cache_name; + //! Cache versioning inline static const std::string Version = "1.0"; @@ -412,7 +415,8 @@ public: //! Constructor SLGridSlab(int kmax, int nmax, int numz, double zmax, - const std::string type="isothermal", bool Verbose=false); + const std::string cachename, const std::string type="isothermal", + bool Verbose=false); //! Destructor ~SLGridSlab(); diff --git a/src/SlabSL.H b/src/SlabSL.H index c394a9b84..c2e940465 100644 --- a/src/SlabSL.H +++ b/src/SlabSL.H @@ -48,6 +48,8 @@ under the time-dependent basis expansion. For a basis fixed in time to the initial time: set to false. + @param cachename is the name of the basis cache file. This is a + required parameter. */ class SlabSL : public PotAccel { @@ -86,6 +88,7 @@ private: int imx, imy, imz, jmax, nnmax; double dfac; std::complex kfac; + std::string cachename; std::vector zfrc, zpot; diff --git a/src/SlabSL.cc b/src/SlabSL.cc index 5e2454af6..4a7620f7d 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -18,7 +18,8 @@ SlabSL::valid_keys = { "zmax", "ngrid", "type", - "self_consistent" + "self_consistent", + "cachename" }; //@{ @@ -35,6 +36,7 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) zmax = 10.0; hslab = 0.2; coef_dump = true; + cachename = ""; #if HAVE_LIBCUDA==1 cuda_aware = true; @@ -42,6 +44,9 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) initialize(); + if (cachename.size()==0) + throw std::runtime_error("SlabSL: you must specify a cachename"); + SLGridSlab::mpi = 1; SLGridSlab::ZBEG = 0.0; SLGridSlab::ZEND = 0.1; @@ -51,7 +56,7 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) // Make the Sturm-Liouville grid and basis functions // - grid = std::make_shared(nnmax, nmaxz, ngrid, zmax, type); + grid = std::make_shared(nnmax, nmaxz, ngrid, zmax, cachename, type); // Test for basis consistency (will generate an exception if maximum // error is out of tolerance) @@ -167,6 +172,7 @@ void SlabSL::initialize() if (conf["hslab"]) hslab = conf["hslab"].as(); if (conf["zmax" ]) zmax = conf["zmax" ].as(); if (conf["type" ]) type = conf["type" ].as(); + if (conf["cachename"]) cachename = conf["cachename"].as(); if (conf["self_consistent"]) { self_consistent = conf["self_consistent"].as(); diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index 3e6a86b0c..ce2ad40c8 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -4,69 +4,18 @@ #include #include -#include - #include "biorth1d.H" #include "SLGridMP2.H" #include "gaussQ.H" #include "localmpi.H" +#include "cxxopts.H" -//=========================================================================== - -void usage(char *prog) -{ - cout << "Usage:\n\n" - << prog << " [options]\n\n" - << setw(15) << "Option" << setw(10) << "Argument" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Description" << endl << endl - << resetiosflags(ios::left) - << setw(15) << "-m or --mpi" << setw(10) << "No" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Turn on MPI for SL computation" << endl - << resetiosflags(ios::left) - << setw(15) << "-t or --Trig" << setw(10) << "No" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Use trigonometric basis" << endl - << resetiosflags(ios::left) - << setw(15) << "-s or --SL" << setw(10) << "No" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Use Sturm-Liouville basis" << endl - << setw(15) << "-T or --type" << setw(10) << "string" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Density target (isothermal, constant, parabolic)" << endl - << resetiosflags(ios::left) - << setw(15) << "-n " << setw(10) << "int" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Number of basis functions" << endl - << resetiosflags(ios::left) - << setw(15) << "-H " << setw(10) << "double" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Slab scale height" << endl - << resetiosflags(ios::left) - << setw(15) << "-k " << setw(10) << "double" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Wave number for Trig basis" << endl - << resetiosflags(ios::left) - << setw(15) << "-x " << setw(10) << "double" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Wave number in X for SL basis" << endl - << resetiosflags(ios::left) - << setw(15) << "-y " << setw(10) << "double" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Wave number in Y for SL basis" << endl - << resetiosflags(ios::left) - << "" << endl; - - exit(0); -} enum BioType1d {Trig, SL}; int main(int argc, char** argv) { - bool use_mpi = false; double KX = 0.5; double H = 0.1; double ZMAX = 1.0; @@ -74,101 +23,40 @@ main(int argc, char** argv) int IKX = 1; int IKY = 3; BioType1d Type = Trig; + std::string cachename = ".slab_sl_cache"; std::string slabID = "iso"; + bool use_mpi = false; - int c; - while (1) { - int this_option_optind = optind ? optind : 1; - int option_index = 0; - static struct option long_options[] = { - {"mpi", 0, 0, 0}, - {"Trig", 0, 0, 0}, - {"SL", 0, 0, 0}, - {0, 0, 0, 0} - }; - - c = getopt_long (argc, argv, "msT:tx:y:k:n:z:H:h", - long_options, &option_index); - - if (c == -1) break; - - switch (c) { - case 0: - { - string optname(long_options[option_index].name); - - if (!optname.compare("mpi")) { - use_mpi = true; - } else if (!optname.compare("Trig")) { - Type = Trig; - } else if (!optname.compare("SL")) { - Type = SL; - } else if (!optname.compare("type")) { - slabID = optarg; - } else { - cout << "Option " << long_options[option_index].name; - if (optarg) cout << " with arg " << optarg; - cout << " is not defined " << endl; - exit(0); - } - } - break; - - case 'm': - use_mpi = true; - break; - - case 's': - Type = SL; - break; - - case 'T': - slabID = optarg; - break; - - case 't': - Type = Trig; - break; - - case 'x': - IKX = atoi(optarg); - break; - - case 'y': - IKY = atoi(optarg); - break; - - case 'k': - KX = atof(optarg); - break; - - case 'z': - ZMAX = atof(optarg); - break; - - case 'H': - H = atof(optarg); - break; - - case 'n': - NMAX = atoi(optarg); - break; - - case 'h': - default: - usage(argv[0]); - } - + cxxopts::Options options("orthochk", "Check orthogonality of 1D basis functions"); + + options.add_options() + ("m,mpi", "Use MPI") + ("s,SL", "Use Sturm-Liouville slab basis") + ("T,type", "Slab type (iso, parabolic, or constant)", cxxopts::value()) + ("t,Trig", "Use trigonometric basis") + ("x,ikx", "IKX for SLGridSlab (default: 1)", cxxopts::value()) + ("y,iky", "IKY for SLGridSlab (default: 3)", cxxopts::value()) + ("k,kx", "KX for OneDTrig (default: 0.5)", cxxopts::value()) + ("z,zmax", "ZMAX for OneDTrig and SLGridSlab (default: 1.0)", cxxopts::value()) + ("H,h", "Scale height H for SLGridSlab (default: 0.1)", cxxopts::value()) + ("n,nmax", "NMAX for SLGridSlab (default: 10)", cxxopts::value()) + ("c,cachename", "Cache file name for SLGridSlab (default: .slab_sl_cache)", cxxopts::value()) + ("h,help", "Print usage"); + + auto result = options.parse(argc, argv); + + if (result.count("mpi")) { + local_init_mpi(argc, argv); + use_mpi = true; } - //=================== - // MPI preliminaries - //=================== - - if (use_mpi) { - local_init_mpi(argc, argv); + if (result.count("help")) { + if (myid==0) + std::cout << options.help() << std::endl; + return 0; } + //=================== // Construct ortho //=================== @@ -190,7 +78,7 @@ main(int argc, char** argv) SLGridSlab::H = H; if (use_mpi) SLGridSlab::mpi = 1; - orthoSL = std::make_shared(KMAX, NMAX, NUMZ, ZMAX, slabID, true); + orthoSL = std::make_shared(KMAX, NMAX, NUMZ, ZMAX, cachename, slabID, true); } break; From 0efa0b326012cd116ceac05867296a0ce062b973 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Apr 2026 20:40:09 -0400 Subject: [PATCH 21/51] Preliminary implemenation of covariance for SlabSL; compiles but untested --- expui/BiorthBasis.H | 62 ++++++++++++- expui/BiorthBasis.cc | 61 ++++++++++++ pyEXP/BasisWrappers.cc | 4 + src/SlabSL.H | 39 +++++++- src/SlabSL.cc | 205 ++++++++++++++++++++++++++++++++++++++++- src/cudaSlabSL.cu | 114 ++++++++++++++++++++++- 6 files changed, 478 insertions(+), 7 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index f772ad9ba..060f1bd99 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1256,7 +1256,7 @@ namespace BasisClasses //! SLGridSlab mesh size int ngrid = 1000; - //! Target model type for SLGridSlab + //! Targ*comet model type for SLGridSlab std::string type = "isothermal"; //! Scale height for slab @@ -1317,6 +1317,20 @@ namespace BasisClasses void computeAccel(double x, double y, double z, Eigen::Ref vstore); + //@{ + //! Covariance structures. First index is T, second is the + //! flattened 3-d k vector + std::vector meanV; + std::vector covrV; + std::vector dvarV; + + void init_covariance(); + void zero_covariance(); + void writeCovarH5Params(HighFive::File& file); + int sampT = 100; + bool diagcov = true; + //@} + public: //! Constructor from YAML node @@ -1373,6 +1387,52 @@ namespace BasisClasses return true; } + //! Write coefficient covariance data to an HDF5 file + virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, double time=0.0) + { + if (covarStore) { + SubsampleCovariance::CovarData elem; + std::get<0>(elem) = sampleCounts; + std::get<1>(elem) = sampleMasses; + + int sampT = meanV.size(); + int Nsize = meanV[0].size(); + + std::get<2>(elem) = SubsampleCovariance::CoefType(sampT, 1, Nsize); + std::get<3>(elem) = SubsampleCovariance::CovrType(sampT, 1, Nsize, Nsize); + + for (int T=0; T(elem)(T, 0, n1) = meanV[T](n1); + for (int n2=0; n2(elem)(T, 0, n1, n2) = covrV[T](n1, n2); + } + } + } + covarStore->writeCoefCovariance(compname, runtag, elem, time); + } else { + throw std::runtime_error("Slab::writeCoefCovariance: covariance storage not initialized"); + } + } + + + //! Enable coefficient covariance computation + virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100, + bool ftype=false, + bool covr_tot=false, + bool covar_in=false) + { + pcavar = pcavar_in; + sampT = sampT_in; + scovr = covr_tot; + covar = covar_in; + if (pcavar) { + init_covariance(); + covarStore = std::make_shared + ([this](HighFive::File& file){this->writeCovarH5Params(file);}, BasisID, ftype, scovr, covar); + } + } + }; /** diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index c1d3b433d..49a4595a0 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3991,6 +3991,57 @@ namespace BasisClasses "fullCovar" }; + void Slab::init_covariance() + { + if (pcavar) { + + int jmax = imx * imy * imz; + + meanV.resize(sampT); + for (auto& v : meanV) { + v.resize(jmax); + } + + if (covar) { + if (diagcov) { + dvarV.resize(sampT); + for (auto& v : dvarV) { + v.resize(jmax); + } + } else { + covrV.resize(sampT); + for (auto& v : covrV) { + v.resize(jmax, jmax); + } + } + } else { + covrV.clear(); + dvarV.clear(); + } + + sampleCounts.resize(sampT); + sampleMasses.resize(sampT); + + zero_covariance(); + } + } + + + void Slab::zero_covariance() + { + for (int T=0; T("rcylmax", HighFive::DataSpace::From(rcylmax)).write(rcylmax); } + void Slab::writeCovarH5Params(HighFive::File& file) + { + file.createAttribute("nminx", HighFive::DataSpace::From(nminx)).write(nminx); + file.createAttribute("nminy", HighFive::DataSpace::From(nminy)).write(nminy); + file.createAttribute("nminz", HighFive::DataSpace::From(nminz)).write(nminz); + file.createAttribute("nmaxx", HighFive::DataSpace::From(nmaxx)).write(nmaxx); + file.createAttribute("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy); + file.createAttribute("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz); + } + void Cube::writeCovarH5Params(HighFive::File& file) { file.createAttribute("nminx", HighFive::DataSpace::From(nminx)).write(nminx); diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 1a997a005..22decd37f 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -781,6 +781,10 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Slab, make_coefs,); } + void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) override { + PYBIND11_OVERRIDE(void, Slab, enableCoefCovariance, pcavar, nsamples, ftype, covar); + } + }; diff --git a/src/SlabSL.H b/src/SlabSL.H index c2e940465..634a73e0c 100644 --- a/src/SlabSL.H +++ b/src/SlabSL.H @@ -81,12 +81,12 @@ private: //! Coefficient tensor for frozen potential (if self_consistent=false) coefType expcofF; - int nminx, nminy; + int nminx, nminy, sampT=0, nint=0; int nmaxx, nmaxy, nmaxz; double zmax, hslab; int imx, imy, imz, jmax, nnmax; - double dfac; + double dfac, last=-std::numeric_limits::max(); std::complex kfac; std::string cachename; @@ -131,7 +131,9 @@ private: thrust::device_vector> dw_coef; thrust::device_vector> df_coef; - void resize_coefs(int N, int osize, int gridSize, int stride); + std::vector>> T_samp; + + void resize_coefs(int N, int osize, int gridSize, int stride, int sampT); }; //! A storage instance @@ -167,6 +169,25 @@ private: //! Default slab type (must be "isothermal", "parabolic", or "constant") std::string type = "isothermal"; + //@{ + //! Covariance structures. First index is T, second is the + //! flattened 3-d k vector + std::vector meanV; + std::vector covrV; + + //@{ + //! Per thread covariance structures + std::vector> meanV1; + std::vector> covrV1; + std::vector workV1; + std::vector countV1; + std::vector massV1; + //@} + + void init_covariance(); + void zero_covariance(); + //@} + //@{ //! Usual evaluation interface void determine_coefficients(void); @@ -229,6 +250,9 @@ private: // Biorth ID static const int ID=1; + //! Write basis-specific parameters to HDF5 covariance file + void writeCovarH5Params(HighFive::File& file); + protected: //! Parse parameters and initialize on first call @@ -250,6 +274,15 @@ public: //! Coefficient output void dump_coefs_h5(const std::string& file); + + /** Return a subsample of the basis coefficients and covariance + elements for matrices for analysis. The default implementation + returns empty vectors. This parallels the implementation in + expui. */ + CovarData getSubsample(); + + //! Indicate whether or not subsampling is supported + virtual bool hasSubsample() { return true; } }; diff --git a/src/SlabSL.cc b/src/SlabSL.cc index 4a7620f7d..45f6b4f09 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -18,6 +18,9 @@ SlabSL::valid_keys = { "zmax", "ngrid", "type", + "nint", + "samplesz", + "subsampleFloat", "self_consistent", "cachename" }; @@ -178,6 +181,15 @@ void SlabSL::initialize() self_consistent = conf["self_consistent"].as(); } else self_consistent = true; + + if (conf["nint"]) { + nint = conf["nint"].as(); + if (nint>0) computeSubsample = true; + } + + if (conf["samplesz"]) { + sampT = conf["samplesz"].as(); + } } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in SlabSL: " @@ -190,11 +202,89 @@ void SlabSL::initialize() throw std::runtime_error("SlabSL::initialze: error parsing YAML"); } + // Initialize covariance + // + init_covariance(); + #if HAVE_LIBCUDA==1 + // Cuda initialization (if needed) + // cuda_initialize(); #endif } +void SlabSL::init_covariance() +{ + if (computeSubsample) { + + meanV.resize(sampT); + for (auto& v : meanV) { + v.resize(jmax); + } + + workV1.resize(nthrds); + for (auto& v : workV1) v.resize(jmax); + + if (fullCovar) { + covrV.resize(sampT); + for (auto& v : covrV) { + v.resize(jmax, jmax); + } + } else { + covrV.clear(); + } + + sampleCounts.resize(sampT); + sampleMasses.resize(sampT); + + meanV1.resize(nthrds); + covrV1.resize(nthrds); + countV1.resize(nthrds); + massV1.resize(nthrds); + + for (int n=0; nsetZero(); } + // Determine whether or not to compute a subsample + if (mstep==0 or mstep==std::numeric_limits::max()) { + if (nint>0 && this_step % nint == 0) { + if (tnow > last) { + requestSubsample = true; + last = tnow; + zero_covariance(); + } + } + } else { + subsampleComputed = false; + } + #if HAVE_LIBCUDA==1 (*barrier)("SlabSL::entering cuda coefficients", __FILE__, __LINE__); if (component->cudaDevice>=0 and use_cuda) { @@ -261,6 +364,48 @@ void SlabSL::determine_coefficients(void) MPI_CXX_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD); } + // Accumulate mean and covariance subsample contributions + // + if (requestSubsample) { + + // Only finalize at the last multistep level + // + if ( (multistep and mlevel==multistep) or multistep==0 ) { + + // Sum over threads + // + for (int n=1; nget_pot(zpot[id], zz, iiy, iix); - for (int iz=0; iz("nminx", HighFive::DataSpace::From(nminx)).write(nminx); + file.createAttribute("nminy", HighFive::DataSpace::From(nminy)).write(nminy); + file.createAttribute("nmaxx", HighFive::DataSpace::From(nmaxx)).write(nmaxx); + file.createAttribute("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy); + file.createAttribute("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz); +} + +PotAccel::CovarData SlabSL::getSubsample() +{ + CovarData elem; + + std::get<0>(elem) = sampleCounts; + std::get<1>(elem) = sampleMasses; + std::get<2>(elem) = Eigen::Tensor, 3>(sampT, 1, jmax); + std::get<3>(elem) = Eigen::Tensor, 4>(sampT, 1, jmax, jmax); + + // Fill the covariance structure with subsamples + for (int T=0; T(elem).chip(T, 0).chip(0, 0) = + Eigen::TensorMap, 1>> + (meanV[T].data(), jmax); + + std::get<3>(elem).chip(T, 0).chip(0, 0) = + Eigen::TensorMap, 2>> + (covrV[T].data(), jmax, jmax); + + } + + return elem; +} + diff --git a/src/cudaSlabSL.cu b/src/cudaSlabSL.cu index 877b5d816..003f120ce 100644 --- a/src/cudaSlabSL.cu +++ b/src/cudaSlabSL.cu @@ -171,7 +171,9 @@ cuFP_t cu_d_xi_to_z(cuFP_t xi) // void SlabSL::cuda_initialize() { - // Nothing + // Turn off full covariance + // + fullCovar = false; } // Copy constants to device @@ -516,7 +518,7 @@ public: } }; -void SlabSL::cudaStorage::resize_coefs(int N, int osize, int gridSize, int stride) +void SlabSL::cudaStorage::resize_coefs(int N, int osize, int gridSize, int stride, int sampT) { // Reserve space for coefficient reduction // @@ -531,6 +533,12 @@ void SlabSL::cudaStorage::resize_coefs(int N, int osize, int gridSize, int strid dN_coef.resize(osize*N); dc_coef.resize(osize*gridSize); dw_coef.resize(osize); // This will stay fixed + + // Resize covariance storage, if sampT>0 + if (sampT) { + T_samp.resize(sampT); + for (int T=0; Tstream), cuS.df_coef.begin(), cuS.df_coef.end(), 0.0); + + if (requestSubsample) { + + cuS.T_samp.resize(sampT); + + for (int T=0; Tstream), + cuS.T_samp[T].begin(), cuS.T_samp[T].end(), 0.0); + } + } } void SlabSL::determine_coefficients_cuda() @@ -684,6 +703,13 @@ void SlabSL::determine_coefficients_cuda() // int sMemSize = BLOCK_SIZE * sizeof(CmplxT); + std::vector::iterator> bm; + if (requestSubsample) { + for (int T=0; T()); thrust::advance(beg, jmax); + + if (requestSubsample) { + + int sN = N/sampT; + int nT = sampT; + + if (sN==0) { // Fail-safe underrun + sN = 1; + nT = N; + } + + for (int T=0; T(mend - mbeg); + massV1 [0][T] += static_cast(thrust::reduce(mbeg, mend)); + } + + // Begin the reduction per grid block + // + /* A reminder to consider implementing strides in reduceSum */ + /* + unsigned int stride1 = s/BLOCK_SIZE/deviceProp.maxGridSize[0] + 1; + unsigned int gridSize1 = s/BLOCK_SIZE/stride1; + + if (s > gridSize1*BLOCK_SIZE*stride1) gridSize1++; + */ + + unsigned int gridSize1 = s/BLOCK_SIZE; + if (s > gridSize1*BLOCK_SIZE) gridSize1++; + + // Sum reduce into gridsize1 blocks taking advantage of + // GPU warp structure + // + + // Mean computation only + // + reduceSumS + <<stream>>> + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), jmax, N, k, k+s); + + // Finish the reduction for this order in parallel + // + thrust::counting_iterator index_begin(0); + thrust::counting_iterator index_end(gridSize1*jmax); + + thrust::reduce_by_key + ( + thrust::cuda::par.on(cs->stream), + thrust::make_transform_iterator(index_begin, key_functor(gridSize1)), + thrust::make_transform_iterator(index_end, key_functor(gridSize1)), + cuS.dc_coef.begin(), thrust::make_discard_iterator(), cuS.dw_coef.begin() + ); + + + thrust::transform(thrust::cuda::par.on(cs->stream), + cuS.dw_coef.begin(), cuS.dw_coef.end(), + bm[T], bm[T], thrust::plus()); + + thrust::advance(bm[T], imz); + } + // END: subsample loop } // use1 += N; // Increment particle count @@ -735,6 +833,18 @@ void SlabSL::determine_coefficients_cuda() // host_coefs = cuS.df_coef; + if (requestSubsample) { + // T loop + // + for (int T=0; T retV = cuS.T_samp[T]; + + for (int n=0; n(retV[n]); + } + } + } + // Create a wavenumber tuple from a flattened index // auto indices = [&](int indx) From 626973ef7c53e64eeb44eb41b4e56bb38f11f80d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 22 Apr 2026 14:37:38 -0400 Subject: [PATCH 22/51] Update for main-bug merge --- pyEXP/BasisWrappers.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 772f2605d..2802e93e0 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -881,7 +881,7 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Slab, make_coefs,); } - void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) override { + void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) { PYBIND11_OVERRIDE(void, Slab, enableCoefCovariance, pcavar, nsamples, ftype, covar); } From 37f60718f29d84e9024f0c3ee6c6d53a96f97bf3 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 30 Apr 2026 16:22:20 -0400 Subject: [PATCH 23/51] Clean up covariance tensor storage allocation in Cube and Slab --- src/Cube.cc | 25 ++++++++++--------------- src/SlabSL.cc | 26 +++++++++++++++++--------- 2 files changed, 27 insertions(+), 24 deletions(-) diff --git a/src/Cube.cc b/src/Cube.cc index f4e184dca..5b7b9e56f 100644 --- a/src/Cube.cc +++ b/src/Cube.cc @@ -1038,27 +1038,22 @@ PotAccel::CovarData Cube::getSubsample() std::get<2>(elem) = Eigen::Tensor, 3>(sampT, 1, osize); std::get<3>(elem) = Eigen::Tensor, 4>(sampT, 1, osize, osize); + if (not fullCovar) std::get<3>(elem).setZero(); + // Fill the covariance structure with subsamples + // for (int T=0; T(elem)(T, 0, n1) = meanV[T](n1); - for (int n2=0; n2(elem)(T, 0, n1, n2) = covrV[T](n1, n2); - else - std::get<3>(elem)(T, 0, n1, n2) = 0.0; - } - */ - std::get<2>(elem).chip(T, 0).chip(0, 0) = Eigen::TensorMap, 1>> (meanV[T].data(), osize); - std::get<3>(elem).chip(T, 0).chip(0, 0) = - Eigen::TensorMap, 2>> - (covrV[T].data(), osize, osize); - + // Only fill the covariance if requested + // + if (fullCovar) { + std::get<3>(elem).chip(T, 0).chip(0, 0) = + Eigen::TensorMap, 2>> + (covrV[T].data(), osize, osize); + } } return elem; diff --git a/src/SlabSL.cc b/src/SlabSL.cc index 45f6b4f09..67c808a53 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -104,11 +104,6 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) << worst << std::endl; } - imx = 1 + 2*nmaxx; // Number of x wavenumber - imy = 1 + 2*nmaxy; // Number of y wavenumbers - imz = nmaxz; // Number of vertical functions - jmax = imx * imy * imz; // Total storage in tensor - expccof.resize(nthrds); for (auto & v : expccof) v.resize(imx, imy, imz); @@ -202,6 +197,13 @@ void SlabSL::initialize() throw std::runtime_error("SlabSL::initialze: error parsing YAML"); } + // Set dimensions of coefficient tensor and total storage in tensor + // + imx = 1 + 2*nmaxx; // Number of x wavenumber + imy = 1 + 2*nmaxy; // Number of y wavenumbers + imz = nmaxz; // Number of vertical functions + jmax = imx * imy * imz; // Total storage in tensor + // Initialize covariance // init_covariance(); @@ -887,16 +889,22 @@ PotAccel::CovarData SlabSL::getSubsample() std::get<2>(elem) = Eigen::Tensor, 3>(sampT, 1, jmax); std::get<3>(elem) = Eigen::Tensor, 4>(sampT, 1, jmax, jmax); + if (not fullCovar) std::get<3>(elem).setZero(); + // Fill the covariance structure with subsamples + // for (int T=0; T(elem).chip(T, 0).chip(0, 0) = Eigen::TensorMap, 1>> (meanV[T].data(), jmax); - std::get<3>(elem).chip(T, 0).chip(0, 0) = - Eigen::TensorMap, 2>> - (covrV[T].data(), jmax, jmax); - + // Use the map to copy the covariance if requested + // + if (fullCovar) { + std::get<3>(elem).chip(T, 0).chip(0, 0) = + Eigen::TensorMap, 2>> + (covrV[T].data(), jmax, jmax); + } } return elem; From c883db15d7dfe5eea759e8a90903e05ed6d21eda Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 1 May 2026 14:57:41 -0400 Subject: [PATCH 24/51] Comments only --- src/SlabSL.cc | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/SlabSL.cc b/src/SlabSL.cc index 67c808a53..be9d4629f 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -34,6 +34,9 @@ static bool cudaAccelOverride = false; SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) { id = "Slab (Sturm-Liouville)"; + + // Default values + // nminx = nminy = 0; nmaxx = nmaxy = nmaxz = 6; zmax = 10.0; @@ -45,6 +48,8 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) cuda_aware = true; #endif + // Parse the YAML and initialize CUDA if needed + // initialize(); if (cachename.size()==0) @@ -68,6 +73,8 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) int kxw = 0, kyw = 0; auto test = grid->orthoCheck(10000); for (int kx=0, indx=0; kx<=nnmax; kx++) { + // We only need the lower triangle of the kx-ky plane since the + // basis is symmetric in kx and ky for (int ky=0; ky<=kx; ky++, indx++) { for (int n1=0; n1 __EXP__::orthoTol) { if (myid==0) std::cout << "SlabSL: orthogonality failure, worst=" << worst @@ -104,12 +113,17 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) << worst << std::endl; } + // Allocate temporary storage for coefficients + // expccof.resize(nthrds); for (auto & v : expccof) v.resize(imx, imy, imz); + // Lx and Ly are both 1.0, so the wavenumbers are 2*pi*n with n an + // integer dfac = 2.0*M_PI; kfac = std::complex(0.0, dfac); + // Arrays for multithreading force evaluation zpot.resize(nthrds); zfrc.resize(nthrds); @@ -315,6 +329,7 @@ void SlabSL::determine_coefficients(void) } // Determine whether or not to compute a subsample + // if (mstep==0 or mstep==std::numeric_limits::max()) { if (nint>0 && this_step % nint == 0) { if (tnow > last) { @@ -481,6 +496,8 @@ void * SlabSL::determine_coefficients_thread(void * arg) std::cerr << "Out of bounds: iiy=" << jj << std::endl; } + // We only populate the lower triangle of the kx-ky plane + // since the basis is symmetric in kx and ky if (iix>=iiy) grid->get_pot(zpot[id], zz, iix, iiy); else From 7c487bc1fb2c0b2451a2d77febd6eee134ebd7e3 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 1 May 2026 14:58:05 -0400 Subject: [PATCH 25/51] Fix option tags --- utils/SL/orthochk.cc | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index ce2ad40c8..19457064c 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -30,18 +30,18 @@ main(int argc, char** argv) cxxopts::Options options("orthochk", "Check orthogonality of 1D basis functions"); options.add_options() - ("m,mpi", "Use MPI") - ("s,SL", "Use Sturm-Liouville slab basis") - ("T,type", "Slab type (iso, parabolic, or constant)", cxxopts::value()) - ("t,Trig", "Use trigonometric basis") - ("x,ikx", "IKX for SLGridSlab (default: 1)", cxxopts::value()) - ("y,iky", "IKY for SLGridSlab (default: 3)", cxxopts::value()) - ("k,kx", "KX for OneDTrig (default: 0.5)", cxxopts::value()) - ("z,zmax", "ZMAX for OneDTrig and SLGridSlab (default: 1.0)", cxxopts::value()) - ("H,h", "Scale height H for SLGridSlab (default: 0.1)", cxxopts::value()) - ("n,nmax", "NMAX for SLGridSlab (default: 10)", cxxopts::value()) + ("m,mpi", "Use MPI") + ("s,SL", "Use Sturm-Liouville slab basis") + ("T,type", "Slab type (iso, parabolic, or constant)", cxxopts::value()) + ("t,Trig", "Use trigonometric basis") + ("x,ikx", "IKX for SLGridSlab (default: 1)", cxxopts::value()) + ("y,iky", "IKY for SLGridSlab (default: 3)", cxxopts::value()) + ("k,kx", "KX for OneDTrig (default: 0.5)", cxxopts::value()) + ("z,zmax", "ZMAX for OneDTrig and SLGridSlab (default: 1.0)", cxxopts::value()) + ("H,height", "Scale height H for SLGridSlab (default: 0.1)", cxxopts::value()) + ("n,nmax", "NMAX for SLGridSlab (default: 10)", cxxopts::value()) ("c,cachename", "Cache file name for SLGridSlab (default: .slab_sl_cache)", cxxopts::value()) - ("h,help", "Print usage"); + ("h,help", "Print usage"); auto result = options.parse(argc, argv); From dc2ae45f22001207dfabec79bca6d57dbd6806e4 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 1 May 2026 16:45:17 -0400 Subject: [PATCH 26/51] Use linear map on table by default --- exputil/SLGridMP2.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index a6c20e260..30a886e86 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1996,7 +1996,7 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, // This could be controlled by a parameter...but at this point, this // is a fixed tuning. - mM = CoordMap::factory(CoordMapTypes::Sech, H); + mM = CoordMap::factory(CoordMapTypes::Linear, H); init_table(); From f7e7be434f6c0b71f9420db3de14170aebfd433f Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 1 May 2026 16:45:30 -0400 Subject: [PATCH 27/51] Update orthogonality and Poisson check --- utils/SL/orthochk.cc | 222 ++++++++++++++++++++++++++----------------- 1 file changed, 133 insertions(+), 89 deletions(-) diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index 19457064c..fa417c54d 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -13,6 +13,16 @@ enum BioType1d {Trig, SL}; +std::map bioTypeMap = { + {"trig", Trig}, + {"sl", SL} +}; + +std::map bioStrMap = { + {Trig, "trig"}, + {SL, "sl"} +}; + int main(int argc, char** argv) { @@ -22,7 +32,7 @@ main(int argc, char** argv) int NMAX = 10; int IKX = 1; int IKY = 3; - BioType1d Type = Trig; + std::string bioTypeStr = "trig"; std::string cachename = ".slab_sl_cache"; std::string slabID = "iso"; bool use_mpi = false; @@ -31,15 +41,15 @@ main(int argc, char** argv) options.add_options() ("m,mpi", "Use MPI") - ("s,SL", "Use Sturm-Liouville slab basis") - ("T,type", "Slab type (iso, parabolic, or constant)", cxxopts::value()) + ("T,type", "Slab type (iso, parabolic, or constant)", cxxopts::value(bioTypeStr)->default_value("trig")) ("t,Trig", "Use trigonometric basis") - ("x,ikx", "IKX for SLGridSlab (default: 1)", cxxopts::value()) - ("y,iky", "IKY for SLGridSlab (default: 3)", cxxopts::value()) - ("k,kx", "KX for OneDTrig (default: 0.5)", cxxopts::value()) - ("z,zmax", "ZMAX for OneDTrig and SLGridSlab (default: 1.0)", cxxopts::value()) - ("H,height", "Scale height H for SLGridSlab (default: 0.1)", cxxopts::value()) - ("n,nmax", "NMAX for SLGridSlab (default: 10)", cxxopts::value()) + ("M,matrix", "Print orthogonality matrix for a particular kx, ky choice") + ("x,ikx", "IKX for SLGridSlab (default: 1)", cxxopts::value(IKX)) + ("y,iky", "IKY for SLGridSlab (default: 3)", cxxopts::value(IKY)) + ("k,kx", "KX for OneDTrig (default: 0.5)", cxxopts::value(KX)) + ("z,zmax", "ZMAX for OneDTrig and SLGridSlab (default: 1.0)", cxxopts::value(ZMAX)) + ("H,height", "Scale height H for SLGridSlab (default: 0.1)", cxxopts::value(H)) + ("n,nmax", "NMAX for SLGridSlab (default: 10)", cxxopts::value(NMAX)) ("c,cachename", "Cache file name for SLGridSlab (default: .slab_sl_cache)", cxxopts::value()) ("h,help", "Print usage"); @@ -57,6 +67,13 @@ main(int argc, char** argv) } + // Determine biorthogonal function type + std::transform(bioTypeStr.begin(), bioTypeStr.end(), bioTypeStr.begin(), + [](unsigned char c){ return std::tolower(c); }); + + BioType1d Type = std::find_if(bioTypeMap.begin(), bioTypeMap.end(), + [bioTypeStr](const std::pair& pair) { return pair.first == bioTypeStr; }) != bioTypeMap.end() ? bioTypeMap[bioTypeStr] : Trig; + //=================== // Construct ortho //=================== @@ -71,7 +88,7 @@ main(int argc, char** argv) case SL: { - const int NUMZ=800; + const int NUMZ=2000; int KMAX = max(IKX+1, IKY+1); SLGridSlab::ZBEG = 0.0; SLGridSlab::ZEND = 0.1; @@ -169,64 +186,66 @@ main(int argc, char** argv) int num; cin >> num; - cout << "N? "; - int N; - cin >> N; - - cout << "dz? "; - double dz; - cin >> dz; + cout << "eps? "; + double eps; + cin >> eps; + double dz = 2.0*ZMAX/num; + double h = dz*eps; double x, z, d1, d2, d3; for (int i=0; ir_to_rb(z); - d1 = ( - ortho->potl(N, i, z+dz) - - ortho->potl(N, i, z)*2.0 + - ortho->potl(N, i, z-dz) - ) / (dz*dz); + for (int n=0; nforce(N, i, z-0.5*dz) - - ortho->force(N, i, z+0.5*dz) - ) / dz; + d1 = ( + ortho->potl(n, i, z+h) - + ortho->potl(n, i, z)*2.0 + + ortho->potl(n, i, z-h) + ) / (dz*dz); + + d2 = ( + ortho->force(n, i, z-0.5*h) - + ortho->force(n, i, z+0.5*h) + ) / dz; - d3 = -KX*KX*ortho->potl(N, i, z); + d3 = -KX*KX*ortho->potl(n, i, z); - out << setw(15) << z - << setw(15) << d1+d3 - << setw(15) << d2+d3 - << setw(15) << -ortho->dens(N, i, z) - << endl; + out << setw(15) << d1+d3 + << setw(15) << d2+d3 + << setw(15) << -ortho->dens(n, i, z); + } + } else { - x = orthoSL->z_to_xi(z); - d1 = ( - orthoSL->get_pot(x+dz, IKX, IKY, N) - - orthoSL->get_pot(x, IKX, IKY, N)*2.0 + - orthoSL->get_pot(x-dz, IKX, IKY, N) - ) / (dz*dz); + Eigen::VectorXd pot0(NMAX), potN(NMAX), potP(NMAX), den0(NMAX); + Eigen::VectorXd frcP(NMAX), frcN(NMAX); - d2 = ( - orthoSL->get_force(x+0.5*dz, IKX, IKY, N) - - orthoSL->get_force(x-0.5*dz, IKX, IKY, N) - ) / dz; + orthoSL->get_pot (pot0, z , IKX, IKY); + orthoSL->get_pot (potN, z-h, IKX, IKY); + orthoSL->get_pot (potP, z+h, IKX, IKY); + orthoSL->get_force(frcN, z-0.5*h, IKX, IKY); + orthoSL->get_force(frcP, z+0.5*h, IKX, IKY); + orthoSL->get_dens (den0, z, IKX, IKY); - d3 = -4.0*M_PI*M_PI*(IKX*IKX+IKY*IKY)* - orthoSL->get_pot(x, IKX, IKY, N); + for (int n=0; nget_dens(x, IKX, IKY, N) - << endl; + d1 = (potP(n) - 2.0*pot0(n) + potN(n)) / (h*h); + d2 = (frcP(n) - frcN(n)) / h; + d3 = -4.0*M_PI*M_PI*(IKX*IKX+IKY*IKY) * pot0(n); + out << setw(15) << d1+d3 + << setw(15) << d2+d3 + << setw(15) << den0(n); + } } + out << endl; } } @@ -240,54 +259,79 @@ main(int argc, char** argv) LegeQuad lw(num); - std::cout << "N1, N2? "; - int N1, N2; - std::cin >> N1; - std::cin >> N2; - - double ximin, ximax; - switch (Type) { - case Trig: - ximin = ortho->r_to_rb(-ZMAX); - ximax = ortho->r_to_rb( ZMAX); - break; - case SL: - ximin = orthoSL->z_to_xi(-ZMAX); - ximax = orthoSL->z_to_xi( ZMAX); - break; + if (result.count("matrix")) { + + Eigen::MatrixXd ans(NMAX, NMAX); + ans.setZero(); + + for (int i=0; ipotl(n1, i, z) * ortho->dens(n2, i, z) * W; + break; + case SL: + ans(n1, n2) += + orthoSL->get_pot(z, IKX, IKY, n1) * orthoSL->get_dens(z, IKX, IKY, n2) * W; + break; + } + } + } + } + + std::cout << std::endl << ans << std::endl; } + else { + std::cout << "N1, N2? "; + int N1, N2; + std::cin >> N1; + std::cin >> N2; - double x, r, ans=0.0; - for (int i=0; ipotl(N1, i, x); - double tmp2 = ortho->dens(N1, i, x); - double tmp3 = ortho->d_r_to_rb(x); - - ans += ortho->potl(N1, i, x)* - ortho->dens(N2, i, x) * - ortho->d_r_to_rb(x) * (ximax - ximin)*lw.weight(i); - } - + ximin = ortho->r_to_rb(-ZMAX); + ximax = ortho->r_to_rb( ZMAX); break; - case SL: + ximin = orthoSL->z_to_xi(-ZMAX); + ximax = orthoSL->z_to_xi( ZMAX); + break; + } + + double z, r, ans=0.0; + for (int i=0; ipotl(N1, i, z); + double tmp2 = ortho->dens(N1, i, z); + + ans += ortho->potl(N1, i, z)* + ortho->dens(N2, i, z) * 2.0*ZMAX*lw.weight(i); + } - ans += orthoSL->get_pot(x, IKX, IKY, N1)* - orthoSL->get_dens(x, IKX, IKY, N2) / - orthoSL->d_xi_to_z(x) * (ximax - ximin)*lw.weight(i); + break; - break; + case SL: + + ans += orthoSL->get_pot(z, IKX, IKY, N1)* + orthoSL->get_dens(z, IKX, IKY, N2) * 2.0*ZMAX*lw.weight(i); + + break; + } } + std::cout << "<" << N1 << "|" << N2 << "> = " << ans << std::endl; } - - cout << "<" << N1 << "|" << N2 << "> = " << ans << endl; } break; From a5a780baa0b98a0f04e0178a3d93ea75039b51ab Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 1 May 2026 17:20:47 -0400 Subject: [PATCH 28/51] Added cmap variable to Slab basis in both EXP and pyEXP --- expui/BiorthBasis.H | 5 ++++- expui/BiorthBasis.cc | 5 ++++- exputil/SLGridMP2.cc | 20 ++++++++++++++++---- include/SLGridMP2.H | 17 ++++++++++++++++- src/SlabSL.H | 5 +++++ src/SlabSL.cc | 4 +++- src/cudaSlabSL.cu | 2 +- utils/SL/orthochk.cc | 4 +++- 8 files changed, 52 insertions(+), 10 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index f8c4ba574..314df3971 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1194,7 +1194,10 @@ namespace BasisClasses //! SLGridSlab mesh size int ngrid = 1000; - //! Targ*comet model type for SLGridSlab + //! Coordinate map type for SLGridSlab + std::string cmap = "linear"; + + //! Target model type for SLGridSlab std::string type = "isothermal"; //! Scale height for slab diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 49a4595a0..654d335f9 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3416,6 +3416,7 @@ namespace BasisClasses "hslab", "zmax", "ngrid", + "cmap", "type", "knots", "verbose", @@ -3478,6 +3479,7 @@ namespace BasisClasses if (conf["hslab"]) hslab = conf["hslab"].as(); if (conf["zmax" ]) zmax = conf["zmax" ].as(); if (conf["ngrid"]) ngrid = conf["ngrid"].as(); + if (conf["cmap" ]) cmap = conf["cmap" ].as(); if (conf["type" ]) type = conf["type" ].as(); if (conf["knots"]) knots = conf["knots"].as(); @@ -3519,7 +3521,8 @@ namespace BasisClasses int nnmax = (nmaxx > nmaxy) ? nmaxx : nmaxy; - ortho = std::make_shared(nnmax, nmaxz, ngrid, zmax, cachename, type); + ortho = std::make_shared(nnmax, nmaxz, ngrid, zmax, cachename, + cmap, type); // Orthogonality sanity check // diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 30a886e86..6b988c0cb 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1973,7 +1973,9 @@ void SLGridSlab::bomb(string oops) // Constructors SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, - const std::string cachename, const std::string TYPE, + const std::string cachename, + const std::string CMAP, + const std::string TYPE, bool VERBOSE) { if (cachename.size()) slab_cache_name = cachename; @@ -1984,6 +1986,7 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, numk = NUMK; nmax = NMAX; numz = NUMZ; + cmap = CMAP; type = TYPE; zmax = ZMAX; @@ -1994,9 +1997,16 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, tbdbg = VERBOSE; - // This could be controlled by a parameter...but at this point, this - // is a fixed tuning. - mM = CoordMap::factory(CoordMapTypes::Linear, H); + + std::transform(cmap.begin(), cmap.end(), cmap.begin(), + [](unsigned char c){ return std::tolower(c); }); + + cmtype = + std::find_if(CoordMapNames.begin(), CoordMapNames.end(), + [&](const std::pair& pair) + { return pair.second == cmap; })->first; + + mM = CoordMap::factory(cmtype, H); init_table(); @@ -2287,6 +2297,7 @@ bool SLGridSlab::ReadH5Cache(void) // Parameter check // + if (not checkStr(cmap, "cmap")) return false; if (not checkStr(type, "type")) return false; if (not checkInt(numk, "numk")) return false; if (not checkInt(nmax, "nmax")) return false; @@ -2382,6 +2393,7 @@ void SLGridSlab::WriteH5Cache(void) // Write parameters // + file.createAttribute ("cmap", HighFive::DataSpace::From(cmap)).write(cmap); file.createAttribute ("type", HighFive::DataSpace::From(type)).write(type); file.createAttribute ("numk", HighFive::DataSpace::From(numk)).write(numk); file.createAttribute ("nmax", HighFive::DataSpace::From(nmax)).write(nmax); diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index 5b6bc4525..543a6fd7e 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -327,6 +327,13 @@ private: //! to finite interval enum CoordMapTypes {Tanh, Sech, Linear}; + //! Coordinate map chooser + const std::map CoordMapNames = { + {Tanh, "tanh" }, + {Sech, "sech" }, + {Linear, "linear"} + }; + //! The base class for all coordinate maps class CoordMap { @@ -386,6 +393,12 @@ private: virtual double d_xi_to_z(double z); }; + //! Coordinate map type (enum) + CoordMapTypes cmtype; + + //! Coordinate map type (string) + std::string cmap = "linear"; + //! Model type std::string type; @@ -415,7 +428,9 @@ public: //! Constructor SLGridSlab(int kmax, int nmax, int numz, double zmax, - const std::string cachename, const std::string type="isothermal", + const std::string cachename, + const std::string coordmap="linear", + const std::string type="isothermal", bool Verbose=false); //! Destructor diff --git a/src/SlabSL.H b/src/SlabSL.H index 634a73e0c..0198aa527 100644 --- a/src/SlabSL.H +++ b/src/SlabSL.H @@ -41,6 +41,8 @@ @param ngrid is the number of grid points in z for the Sturm-Liouville solver (default 1000) + @param cmap is the coordinate mapping for the SL solver (must be one of "tanh", "sech", or "linear", default "linear") + @param type is the type of slab to solve for (default "isothermal", must be "isothermal", "parabolic", or "constant") @@ -166,6 +168,9 @@ private: //! Default number of grid points for SLGridSlab int ngrid = 1000; + //! Default coordinate mapping for SLGridSlab (must be "tanh", "sech", or "linear") + std::string cmap = "linear"; + //! Default slab type (must be "isothermal", "parabolic", or "constant") std::string type = "isothermal"; diff --git a/src/SlabSL.cc b/src/SlabSL.cc index be9d4629f..20923e2a4 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -17,6 +17,7 @@ SlabSL::valid_keys = { "hslab", "zmax", "ngrid", + "cmap", "type", "nint", "samplesz", @@ -64,7 +65,8 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) // Make the Sturm-Liouville grid and basis functions // - grid = std::make_shared(nnmax, nmaxz, ngrid, zmax, cachename, type); + grid = std::make_shared(nnmax, nmaxz, ngrid, zmax, cachename, + cmap, type); // Test for basis consistency (will generate an exception if maximum // error is out of tolerance) diff --git a/src/cudaSlabSL.cu b/src/cudaSlabSL.cu index 003f120ce..1f6363b86 100644 --- a/src/cudaSlabSL.cu +++ b/src/cudaSlabSL.cu @@ -216,7 +216,7 @@ void SlabSL::initialize_constants() __FILE__, __LINE__, "Error copying slabNum"); // Sech map - int Cmap = 1; + int Cmap = static_cast(cmtype); cuda_safe_call(cudaMemcpyToSymbol(slabCmap, &Cmap, sizeof(int), size_t(0), cudaMemcpyHostToDevice), diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index fa417c54d..08c5f67eb 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -34,6 +34,7 @@ main(int argc, char** argv) int IKY = 3; std::string bioTypeStr = "trig"; std::string cachename = ".slab_sl_cache"; + std::string cmap = "linear"; std::string slabID = "iso"; bool use_mpi = false; @@ -95,7 +96,8 @@ main(int argc, char** argv) SLGridSlab::H = H; if (use_mpi) SLGridSlab::mpi = 1; - orthoSL = std::make_shared(KMAX, NMAX, NUMZ, ZMAX, cachename, slabID, true); + orthoSL = std::make_shared(KMAX, NMAX, NUMZ, ZMAX, cachename, + cmap, slabID, true); } break; From bab8d36ec8e5a51aad35bd329cd68bba18d1240a Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 1 May 2026 19:08:43 -0400 Subject: [PATCH 29/51] Some updates for covariance cuda code; untested but compiling --- src/SlabSL.H | 4 +- src/cudaSlabSL.cu | 122 +++++++++++++++++++++++++++------------------- 2 files changed, 74 insertions(+), 52 deletions(-) diff --git a/src/SlabSL.H b/src/SlabSL.H index 0198aa527..fb0a6200c 100644 --- a/src/SlabSL.H +++ b/src/SlabSL.H @@ -133,9 +133,11 @@ private: thrust::device_vector> dw_coef; thrust::device_vector> df_coef; + thrust::device_vector u_d; + std::vector>> T_samp; - void resize_coefs(int N, int osize, int gridSize, int stride, int sampT); + void resize_coefs(int N, int osize, int gridSize, int sampT); }; //! A storage instance diff --git a/src/cudaSlabSL.cu b/src/cudaSlabSL.cu index 1f6363b86..3213b38e4 100644 --- a/src/cudaSlabSL.cu +++ b/src/cudaSlabSL.cu @@ -1,5 +1,6 @@ // -*- C++ -*- +#include #include #include @@ -215,8 +216,15 @@ void SlabSL::initialize_constants() size_t(0), cudaMemcpyHostToDevice), __FILE__, __LINE__, "Error copying slabNum"); - // Sech map - int Cmap = static_cast(cmtype); + // Coordinate map + std::map CoordMap = + { + {"tanh" , 0}, + {"sech" , 1}, + {"linear", 2} + }; + + int Cmap = CoordMap[cmap]; cuda_safe_call(cudaMemcpyToSymbol(slabCmap, &Cmap, sizeof(int), size_t(0), cudaMemcpyHostToDevice), @@ -246,7 +254,7 @@ void SlabSL::initialize_constants() __global__ void coefKernelSlab (dArray P, dArray I, dArray coef, - dArray tex, int stride, PII lohi) + dArray tex, dArray used, int stride, PII lohi) { // Thread ID // @@ -271,6 +279,8 @@ __global__ void coefKernelSlab cuFP_t pos[3] = {p.pos[0], p.pos[1], p.pos[2]}; cuFP_t mm = - p.mass * 2.0*slabDfac; + used._v[i] = mm; // Mark particle as used + // Restore particles to the unit slab // for (int k=0; k<2; k++) { @@ -518,7 +528,7 @@ public: } }; -void SlabSL::cudaStorage::resize_coefs(int N, int osize, int gridSize, int stride, int sampT) +void SlabSL::cudaStorage::resize_coefs(int N, int osize, int gridSize, int sampT) { // Reserve space for coefficient reduction // @@ -534,6 +544,9 @@ void SlabSL::cudaStorage::resize_coefs(int N, int osize, int gridSize, int strid dc_coef.resize(osize*gridSize); dw_coef.resize(osize); // This will stay fixed + u_d.resize(N); + + // Resize covariance storage, if sampT>0 // Resize covariance storage, if sampT>0 if (sampT) { T_samp.resize(sampT); @@ -640,6 +653,11 @@ void SlabSL::determine_coefficients_cuda() // cuda_zero_coefs(); + // Zero mass accumulation array + // + thrust::fill(cuS.u_d.begin(), cuS.u_d.end(), 0.0); + + // Get sorted particle range for mlevel // Get sorted particle range for mlevel // PII lohi = component->CudaGetLevelRange(mlevel, mlevel), cur; @@ -712,7 +730,7 @@ void SlabSL::determine_coefficients_cuda() // Adjust cached storage, if necessary // - cuS.resize_coefs(N, jmax, gridSize, stride); + cuS.resize_coefs(N, jmax, gridSize, sampT); // Compute the coefficient contribution for each order // @@ -720,7 +738,8 @@ void SlabSL::determine_coefficients_cuda() coefKernelSlab<<stream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), toKernel(t_d), stride, cur); + toKernel(cuS.dN_coef), toKernel(t_d), toKernel(cuS.u_d), + stride, cur); // Begin the reduction by blocks [perhaps this should use a // stride?] @@ -778,52 +797,52 @@ void SlabSL::determine_coefficients_cuda() countV1[0][T] += static_cast(mend - mbeg); massV1 [0][T] += static_cast(thrust::reduce(mbeg, mend)); - } - // Begin the reduction per grid block - // - /* A reminder to consider implementing strides in reduceSum */ - /* - unsigned int stride1 = s/BLOCK_SIZE/deviceProp.maxGridSize[0] + 1; - unsigned int gridSize1 = s/BLOCK_SIZE/stride1; - - if (s > gridSize1*BLOCK_SIZE*stride1) gridSize1++; - */ + // Begin the reduction per grid block + // + /* A reminder to consider implementing strides in reduceSum */ + /* + unsigned int stride1 = s/BLOCK_SIZE/deviceProp.maxGridSize[0] + 1; + unsigned int gridSize1 = s/BLOCK_SIZE/stride1; + + if (s > gridSize1*BLOCK_SIZE*stride1) gridSize1++; + */ - unsigned int gridSize1 = s/BLOCK_SIZE; - if (s > gridSize1*BLOCK_SIZE) gridSize1++; + unsigned int gridSize1 = s/BLOCK_SIZE; + if (s > gridSize1*BLOCK_SIZE) gridSize1++; - // Sum reduce into gridsize1 blocks taking advantage of - // GPU warp structure - // - - // Mean computation only - // - reduceSumS - <<stream>>> - (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), jmax, N, k, k+s); - - // Finish the reduction for this order in parallel - // - thrust::counting_iterator index_begin(0); - thrust::counting_iterator index_end(gridSize1*jmax); - - thrust::reduce_by_key - ( - thrust::cuda::par.on(cs->stream), - thrust::make_transform_iterator(index_begin, key_functor(gridSize1)), - thrust::make_transform_iterator(index_end, key_functor(gridSize1)), - cuS.dc_coef.begin(), thrust::make_discard_iterator(), cuS.dw_coef.begin() - ); - - - thrust::transform(thrust::cuda::par.on(cs->stream), - cuS.dw_coef.begin(), cuS.dw_coef.end(), - bm[T], bm[T], thrust::plus()); - - thrust::advance(bm[T], imz); + // Sum reduce into gridsize1 blocks taking advantage of + // GPU warp structure + // + + // Mean computation only + // + reduceSumS + <<stream>>> + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), jmax, N, k, k+s); + + // Finish the reduction for this order in parallel + // + thrust::counting_iterator index_begin(0); + thrust::counting_iterator index_end(gridSize1*jmax); + + thrust::reduce_by_key + ( + thrust::cuda::par.on(cs->stream), + thrust::make_transform_iterator(index_begin, key_functor(gridSize1)), + thrust::make_transform_iterator(index_end, key_functor(gridSize1)), + cuS.dc_coef.begin(), thrust::make_discard_iterator(), cuS.dw_coef.begin() + ); + + + thrust::transform(thrust::cuda::par.on(cs->stream), + cuS.dw_coef.begin(), cuS.dw_coef.end(), + bm[T], bm[T], thrust::plus()); + + thrust::advance(bm[T], imz); + } + // END: subsample loop } - // END: subsample loop } // use1 += N; // Increment particle count @@ -839,7 +858,7 @@ void SlabSL::determine_coefficients_cuda() for (int T=0; T retV = cuS.T_samp[T]; - for (int n=0; n(retV[n]); } } @@ -1242,7 +1261,7 @@ void SlabSL::multistep_update_cuda() // Adjust cached storage, if necessary // - cuS.resize_coefs(N, jmax, gridSize, stride); + cuS.resize_coefs(N, jmax, gridSize, sampT); // Compute the coefficient contribution for each order // @@ -1252,7 +1271,8 @@ void SlabSL::multistep_update_cuda() // coefKernelSlab<<stream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), toKernel(t_d), stride, cur); + toKernel(cuS.dN_coef), toKernel(t_d), toKernel(cuS.u_d), + stride, cur); unsigned int gridSize1 = N/BLOCK_SIZE; if (N > gridSize1*BLOCK_SIZE) gridSize1++; From 86baa995339a34fc36b7cd61d627863bc344fc71 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 1 May 2026 19:23:28 -0400 Subject: [PATCH 30/51] Use mass, not mass scaled for biorthogonal function accumulation, for covariance analysis --- src/Component.cc | 2 +- src/cudaCube.cu | 3 +-- src/cudaSlabSL.cu | 6 ++---- 3 files changed, 4 insertions(+), 7 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index c50884ce6..af9109217 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1393,7 +1393,7 @@ void Component::initialize(void) cout << "---- Component <" << name << ">: "; - if (nlevel<0) + if (nlevel<=0) std::cout << "no multistep level reporting"; else std::cout << "multistep level reporting every " << nlevel << " steps"; diff --git a/src/cudaCube.cu b/src/cudaCube.cu index 9f3098066..0d13a190e 100644 --- a/src/cudaCube.cu +++ b/src/cudaCube.cu @@ -1,5 +1,3 @@ -// -*- C++ -*- - #include #include @@ -1462,4 +1460,5 @@ void Cube::destroy_cuda() // Nothing } +// -*- C++ -*- diff --git a/src/cudaSlabSL.cu b/src/cudaSlabSL.cu index 3213b38e4..6e7eda3e8 100644 --- a/src/cudaSlabSL.cu +++ b/src/cudaSlabSL.cu @@ -1,5 +1,3 @@ -// -*- C++ -*- - #include #include @@ -279,7 +277,7 @@ __global__ void coefKernelSlab cuFP_t pos[3] = {p.pos[0], p.pos[1], p.pos[2]}; cuFP_t mm = - p.mass * 2.0*slabDfac; - used._v[i] = mm; // Mark particle as used + used._v[i] = p.mass; // Mark particle as used // Restore particles to the unit slab // @@ -1328,4 +1326,4 @@ void SlabSL::destroy_cuda() // Nothing } - +// -*- C++ -*- From 3766f97aef55e2eb7a590aecf1f7cf1d6fa13b98 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 2 May 2026 10:12:58 -0400 Subject: [PATCH 31/51] Allow check to use all the coordinate mappings for testing --- utils/SL/orthochk.cc | 46 +++++++++++++++++++++++++++++++++++++------- 1 file changed, 39 insertions(+), 7 deletions(-) diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index 08c5f67eb..a2f36af78 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -42,8 +42,7 @@ main(int argc, char** argv) options.add_options() ("m,mpi", "Use MPI") - ("T,type", "Slab type (iso, parabolic, or constant)", cxxopts::value(bioTypeStr)->default_value("trig")) - ("t,Trig", "Use trigonometric basis") + ("T,type", "Slab type (trig or SL)", cxxopts::value(bioTypeStr)->default_value("trig")) ("M,matrix", "Print orthogonality matrix for a particular kx, ky choice") ("x,ikx", "IKX for SLGridSlab (default: 1)", cxxopts::value(IKX)) ("y,iky", "IKY for SLGridSlab (default: 3)", cxxopts::value(IKY)) @@ -51,6 +50,8 @@ main(int argc, char** argv) ("z,zmax", "ZMAX for OneDTrig and SLGridSlab (default: 1.0)", cxxopts::value(ZMAX)) ("H,height", "Scale height H for SLGridSlab (default: 0.1)", cxxopts::value(H)) ("n,nmax", "NMAX for SLGridSlab (default: 10)", cxxopts::value(NMAX)) + ("C,cmap", "Coordinate map for SLGridSlab (linear, tanh, or sech; default: linear)", cxxopts::value(cmap)->default_value("linear")) + ("s,slabid", "Slab model ID for SLGridSlab (iso, parabolic, or constant; default: iso)", cxxopts::value(slabID)->default_value("iso")) ("c,cachename", "Cache file name for SLGridSlab (default: .slab_sl_cache)", cxxopts::value()) ("h,help", "Print usage"); @@ -75,6 +76,18 @@ main(int argc, char** argv) BioType1d Type = std::find_if(bioTypeMap.begin(), bioTypeMap.end(), [bioTypeStr](const std::pair& pair) { return pair.first == bioTypeStr; }) != bioTypeMap.end() ? bioTypeMap[bioTypeStr] : Trig; + // Check for valid coordinate map choice + std::transform(cmap.begin(), cmap.end(), cmap.begin(), + [](unsigned char c){ return std::tolower(c); }); + + std::vector valid_cmaps = { "tanh", "sech", "linear" }; + if (valid_cmaps.end() == + std::find_if(valid_cmaps.begin(), valid_cmaps.end(), + [cmap](const std::string& vmap) { return vmap == cmap; })) { + std::cerr << "Invalid coordinate map choice: " << cmap << std::endl; + return 1; + } + //=================== // Construct ortho //=================== @@ -123,7 +136,8 @@ main(int argc, char** argv) cout << "1: Print out density, potential pairs" << endl; cout << "2: Check density" << endl; cout << "3: Check orthogonality" << endl; - cout << "4: Quit" << endl; + cout << "4: Internal ortho check" << endl; + cout << "5: Quit" << endl; cout << "?? "; cin >> iwhich; @@ -335,9 +349,27 @@ main(int argc, char** argv) std::cout << "<" << N1 << "|" << N2 << "> = " << ans << std::endl; } } - - break; - + case 4: + { + std::vector orthoMat; + switch (Type) { + case Trig: + std::cout << "No internal orthogonality check for OneDTrig" + << std::endl; + break; + case SL: + { + int knots; + std::cout << "Number of knots? "; + std::cin >> knots; + auto orthoMat = orthoSL->orthoCheck(knots); + std::cout << "Orthogonality check for SLGridSlab" << std::endl; + for (size_t i=0; i Date: Sat, 2 May 2026 12:01:51 -0400 Subject: [PATCH 32/51] Add numz as a variable --- utils/SL/orthochk.cc | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index a2f36af78..9f1bdebd6 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -32,6 +32,7 @@ main(int argc, char** argv) int NMAX = 10; int IKX = 1; int IKY = 3; + int numz = 2000; std::string bioTypeStr = "trig"; std::string cachename = ".slab_sl_cache"; std::string cmap = "linear"; @@ -50,6 +51,7 @@ main(int argc, char** argv) ("z,zmax", "ZMAX for OneDTrig and SLGridSlab (default: 1.0)", cxxopts::value(ZMAX)) ("H,height", "Scale height H for SLGridSlab (default: 0.1)", cxxopts::value(H)) ("n,nmax", "NMAX for SLGridSlab (default: 10)", cxxopts::value(NMAX)) + ("N,numz", "Number of z points for SLGridSlab (default: 2000)", cxxopts::value(numz)->default_value("2000")) ("C,cmap", "Coordinate map for SLGridSlab (linear, tanh, or sech; default: linear)", cxxopts::value(cmap)->default_value("linear")) ("s,slabid", "Slab model ID for SLGridSlab (iso, parabolic, or constant; default: iso)", cxxopts::value(slabID)->default_value("iso")) ("c,cachename", "Cache file name for SLGridSlab (default: .slab_sl_cache)", cxxopts::value()) @@ -102,14 +104,13 @@ main(int argc, char** argv) case SL: { - const int NUMZ=2000; int KMAX = max(IKX+1, IKY+1); SLGridSlab::ZBEG = 0.0; SLGridSlab::ZEND = 0.1; SLGridSlab::H = H; if (use_mpi) SLGridSlab::mpi = 1; - orthoSL = std::make_shared(KMAX, NMAX, NUMZ, ZMAX, cachename, + orthoSL = std::make_shared(KMAX, NMAX, numz, ZMAX, cachename, cmap, slabID, true); } break; @@ -141,7 +142,7 @@ main(int argc, char** argv) cout << "?? "; cin >> iwhich; - if (iwhich < 1 || iwhich > 4) iwhich = 4; + if (iwhich < 1 || iwhich > 5) iwhich = 5; switch(iwhich) { case 1: From 8e25b530a50312216b27ee5369cb10d919a580df Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 2 May 2026 12:35:00 -0400 Subject: [PATCH 33/51] Missing parameter parsing for cmap --- src/SlabSL.cc | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/SlabSL.cc b/src/SlabSL.cc index 20923e2a4..aad87f62e 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -185,6 +185,7 @@ void SlabSL::initialize() if (conf["ngrid"]) ngrid = conf["ngrid"].as(); if (conf["hslab"]) hslab = conf["hslab"].as(); if (conf["zmax" ]) zmax = conf["zmax" ].as(); + if (conf["cmap" ]) cmap = conf["cmap" ].as(); if (conf["type" ]) type = conf["type" ].as(); if (conf["cachename"]) cachename = conf["cachename"].as(); @@ -446,9 +447,9 @@ void * SlabSL::determine_coefficients_thread(void * arg) std::complex stepx, stepy; unsigned nbodies = component->levlist[mlevel].size(); - int id = *((int*)arg); - int nbeg = nbodies*id/nthrds; - int nend = nbodies*(id+1)/nthrds; + int id = *((int*)arg); + int nbeg = nbodies*id/nthrds; + int nend = nbodies*(id+1)/nthrds; double adb = cC->Adiabatic(); for (int q=nbeg; q Date: Sat, 2 May 2026 12:35:23 -0400 Subject: [PATCH 34/51] Add periodic box wrapping to cuda implementation for testing --- src/cudaSlabSL.cu | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/src/cudaSlabSL.cu b/src/cudaSlabSL.cu index 6e7eda3e8..3de4045a4 100644 --- a/src/cudaSlabSL.cu +++ b/src/cudaSlabSL.cu @@ -250,6 +250,44 @@ void SlabSL::initialize_constants() } +__global__ void wrapKernelSlab +(dArray P, dArray I, int stride, PII lohi) +{ + // Thread ID + // + const int tid = blockDim.x * blockIdx.x + threadIdx.x; + const int N = lohi.second - lohi.first; + + for (int n=0; n=P._s) printf("out of bounds: %s:%d\n", __FILE__, __LINE__); +#endif + cudaParticle & p = P._v[I._v[npart]]; + + // Truncate to box with sides in [0,1] + // + for (int k=0; k<2; k++) { + if (p.pos[k]<0.0) + p.pos[k] += floor(-p.pos[k]) + 1.0; + else + p.pos[k] += -floor(p.pos[k]); + } + } + // END: particle index limit + } + // END: stride loop +} + + __global__ void coefKernelSlab (dArray P, dArray I, dArray coef, dArray tex, dArray used, int stride, PII lohi) @@ -730,6 +768,11 @@ void SlabSL::determine_coefficients_cuda() // cuS.resize_coefs(N, jmax, gridSize, sampT); + // Wrap the coordinates to the unit slab + // + wrapKernelSlab<<stream>>> + (toKernel(cs->cuda_particles), toKernel(cs->indx1), stride, cur); + // Compute the coefficient contribution for each order // auto beg = cuS.df_coef.begin(); From e39cfc8da654d6e584ce1043435d6f10f9760465 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 3 May 2026 10:51:43 -0400 Subject: [PATCH 35/51] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 2802e93e0..b81ffc4e1 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -881,8 +881,9 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Slab, make_coefs,); } - void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) { - PYBIND11_OVERRIDE(void, Slab, enableCoefCovariance, pcavar, nsamples, ftype, covar); + void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covr_tot, bool covar) override + { + PYBIND11_OVERRIDE(void, Slab, enableCoefCovariance, pcavar, nsamples, ftype, covr_tot, covar); } }; From 6590d598bf61f17b6ef94925a283bf167c1fc89c Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 3 May 2026 10:52:28 -0400 Subject: [PATCH 36/51] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/SL/orthochk.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index 9f1bdebd6..30a54a6e0 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -54,7 +54,7 @@ main(int argc, char** argv) ("N,numz", "Number of z points for SLGridSlab (default: 2000)", cxxopts::value(numz)->default_value("2000")) ("C,cmap", "Coordinate map for SLGridSlab (linear, tanh, or sech; default: linear)", cxxopts::value(cmap)->default_value("linear")) ("s,slabid", "Slab model ID for SLGridSlab (iso, parabolic, or constant; default: iso)", cxxopts::value(slabID)->default_value("iso")) - ("c,cachename", "Cache file name for SLGridSlab (default: .slab_sl_cache)", cxxopts::value()) + ("c,cachename", "Cache file name for SLGridSlab (default: .slab_sl_cache)", cxxopts::value(cachename)->default_value(".slab_sl_cache")) ("h,help", "Print usage"); auto result = options.parse(argc, argv); From a111a904857c1082dd0465a61b12a80f75d65e16 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 3 May 2026 10:52:59 -0400 Subject: [PATCH 37/51] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/SL/orthochk.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index 30a54a6e0..b24b8bc9c 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -224,12 +224,12 @@ main(int argc, char** argv) ortho->potl(n, i, z+h) - ortho->potl(n, i, z)*2.0 + ortho->potl(n, i, z-h) - ) / (dz*dz); + ) / (h*h); d2 = ( ortho->force(n, i, z-0.5*h) - ortho->force(n, i, z+0.5*h) - ) / dz; + ) / h; d3 = -KX*KX*ortho->potl(n, i, z); From 7b4f5d715ed8f45cae98a7bf1e538734a11f96b5 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 3 May 2026 14:53:24 +0000 Subject: [PATCH 38/51] Document nint and samplesz subsampling params in SlabSL YAML docstring Agent-Logs-Url: https://github.com/EXP-code/EXP/sessions/61daa681-8beb-4cbd-877f-f58ff1d18b1a Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- src/SlabSL.H | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/SlabSL.H b/src/SlabSL.H index fb0a6200c..fee9220ec 100644 --- a/src/SlabSL.H +++ b/src/SlabSL.H @@ -50,6 +50,12 @@ under the time-dependent basis expansion. For a basis fixed in time to the initial time: set to false. + @param nint integer number of steps between subsample computations + (default: 0, no subsampling) + + @param samplesz integer number of subsample bins (default: 0, + no subsampling). Must be greater than 0 when `nint>0`. + @param cachename is the name of the basis cache file. This is a required parameter. */ From d87a355d6f35c78509dedb72093aeda79b45144d Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 3 May 2026 10:53:37 -0400 Subject: [PATCH 39/51] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.H | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 314df3971..96dbfbbbe 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1346,7 +1346,13 @@ namespace BasisClasses for (int n1=0; n1(elem)(T, 0, n1) = meanV[T](n1); for (int n2=0; n2(elem)(T, 0, n1, n2) = covrV[T](n1, n2); + double cov = 0.0; + if (diagcov) { + if (n1 == n2) cov = dvarV[T](n1); + } else if (covar && T < covrV.size()) { + cov = covrV[T](n1, n2); + } + std::get<3>(elem)(T, 0, n1, n2) = cov; } } } From 5e9e60bc207c05bc761eae7497e25fb09e1e094d Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 3 May 2026 10:54:36 -0400 Subject: [PATCH 40/51] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- include/SLGridMP2.H | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index 543a6fd7e..c4de61250 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -433,6 +433,12 @@ public: const std::string type="isothermal", bool Verbose=false); + //! Backward-compatible constructor for legacy call sites + SLGridSlab(int kmax, int nmax, int numz, double zmax, + const std::string type, + bool Verbose=false) + : SLGridSlab(kmax, nmax, numz, zmax, default_cache, "linear", type, Verbose) {} + //! Destructor ~SLGridSlab(); From 5dbe30af4efcee1a28cecca8b79ea2fa597ee34a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 3 May 2026 15:05:54 +0000 Subject: [PATCH 41/51] Remove unused outer orthoMat declaration in orthochk.cc case 4 Agent-Logs-Url: https://github.com/EXP-code/EXP/sessions/61aca825-1c82-4e34-af64-6de5d2314370 Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- utils/SL/orthochk.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index b24b8bc9c..86a3ccdbc 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -352,7 +352,6 @@ main(int argc, char** argv) } case 4: { - std::vector orthoMat; switch (Type) { case Trig: std::cout << "No internal orthogonality check for OneDTrig" From 0cc1bc5c5428443bc82b80395de8bbe4dabe64b6 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 3 May 2026 15:09:14 +0000 Subject: [PATCH 42/51] Fix SlabSL.ortho diagnostic gate, CoordMap safety, and bonnerebert -V option Agent-Logs-Url: https://github.com/EXP-code/EXP/sessions/437a8908-e2a0-48c1-9224-8263a765ef4e Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- src/SlabSL.cc | 7 +++++-- src/cudaSlabSL.cu | 9 ++++++++- utils/ICs/bonnerebert.cc | 5 +++++ 3 files changed, 18 insertions(+), 3 deletions(-) diff --git a/src/SlabSL.cc b/src/SlabSL.cc index aad87f62e..192ba348c 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -92,8 +92,11 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) } } - // Diagnostic output for debugging - if (true) { + // Diagnostic output for debugging. + // +--- Set to 'true' only for deep debugging; 'false' for production. + // | + // v + if (false) { std::ofstream tmp("SlabSL.ortho"); for (int kx=0, indx=0; kx<=nnmax; kx++) { for (int ky=0; ky<=kx; ky++, indx++) { diff --git a/src/cudaSlabSL.cu b/src/cudaSlabSL.cu index 3de4045a4..806dfd914 100644 --- a/src/cudaSlabSL.cu +++ b/src/cudaSlabSL.cu @@ -222,7 +222,14 @@ void SlabSL::initialize_constants() {"linear", 2} }; - int Cmap = CoordMap[cmap]; + auto it = CoordMap.find(cmap); + if (it == CoordMap.end()) { + std::ostringstream sout; + sout << "cudaSlabSL: unknown coordinate map type '" << cmap + << "'. Valid values are: tanh, sech, linear"; + throw std::runtime_error(sout.str()); + } + int Cmap = it->second; cuda_safe_call(cudaMemcpyToSymbol(slabCmap, &Cmap, sizeof(int), size_t(0), cudaMemcpyHostToDevice), diff --git a/utils/ICs/bonnerebert.cc b/utils/ICs/bonnerebert.cc index c20c38ff4..8d7cceb00 100644 --- a/utils/ICs/bonnerebert.cc +++ b/utils/ICs/bonnerebert.cc @@ -364,6 +364,7 @@ decode_switches (int argc, char **argv) "R:" /* runit */ "N:" /* number */ "S:" /* seed */ + "V" /* version */ "h", /* help */ long_options, (int *) 0)) != EOF) { @@ -399,6 +400,10 @@ decode_switches (int argc, char **argv) case 'h': usage (0); + case 'V': + cout << program_name << " version " << VERSION << endl; + exit (0); + default: usage (-1); } From 227e8594c36d29facb45115030a1c55bb174a1d7 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 3 May 2026 12:14:27 -0400 Subject: [PATCH 43/51] Consistency fixes between CPU and GPU --- src/SlabSL.cc | 18 +++++++++--------- src/cudaSlabSL.cu | 26 +++++++++++++++++--------- 2 files changed, 26 insertions(+), 18 deletions(-) diff --git a/src/SlabSL.cc b/src/SlabSL.cc index 192ba348c..9c9e88116 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -313,8 +313,8 @@ void SlabSL::determine_coefficients(void) // n=-nmax,-nmax+1,...,0,...,nmax-1,nmax in a single array for each // dimension with z dimension changing most rapidly - // Clean - + // Clean acculumulation arrays + // for (int i=0; i(nmaxx)*kfac*cC->Pos(i, 0)); starty = exp(static_cast(nmaxy)*kfac*cC->Pos(i, 1)); - double zz = cC->Pos(i, 2), mm = -4.0*M_PI * cC->Mass(i) * adb; + double zz = cC->Pos(i, 2), ms = cC->Mass(i); + double mm = -4.0*M_PI * ms * adb; for (facx=startx, ix=0; ix P, dArray I, // Flip sign for antisymmetric basis functions // - int sign = 1; - if (pos[2]<0 && 2*(n/2)==n) sign = -1; + int signV = 1, signF = 1; + + // potential: flip odd modes + if (pos[2]<0 && 2*(n/2)!=n) signV = -1; + + // z-force: flip even modes + if (pos[2]<0 && 2*(n/2)==n) signF = -1; int k = offset + n; @@ -523,7 +531,7 @@ forceKernelSlab(dArray P, dArray I, a*int2_as_double(tex1D(tex._v[k], in0 )) + b*int2_as_double(tex1D(tex._v[k], in0+1)) #endif - ) * p0 * sign; + ) * p0 * signV; cuFP_t f = ( #if cuREAL == 4 @@ -535,7 +543,7 @@ forceKernelSlab(dArray P, dArray I, -2.0*s*int2_as_double(tex1D(tex._v[k], jn0)) * int2_as_double(tex1D(tex._v[0], jn0)) + (s + 0.5)*int2_as_double(tex1D(tex._v[k], jn0+1)) * int2_as_double(tex1D(tex._v[0], jn0+1)) #endif - ) * sign * dx; + ) * signF * dx; fac = X * Y * v * coef._v[slabIndex(ii, jj, n)]; facf = X * Y * f * coef._v[slabIndex(ii, jj, n)]; From 8ff45d18f2e0a75ca1f6e2e6de88b6847d9891f8 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 3 May 2026 12:28:22 -0400 Subject: [PATCH 44/51] Fix up some missing variables and type mismatches introduced by signatures changes proposed by Copilot --- expui/BiorthBasis.H | 2 +- include/SLGridMP2.H | 6 +++++- utils/SL/slabchk.cc | 4 +++- 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 96dbfbbbe..1be1e269f 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1346,7 +1346,7 @@ namespace BasisClasses for (int n1=0; n1(elem)(T, 0, n1) = meanV[T](n1); for (int n2=0; n2 cov = 0.0; if (diagcov) { if (n1 == n2) cov = dvarV[T](n1); } else if (covar && T < covrV.size()) { diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index c4de61250..fa4091efc 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -405,6 +405,9 @@ private: //! The current map created by the SLGridSlab constructor std::unique_ptr mM; + //! Default cache file name + const std::string default_cache = ".slgrid_slab_cache"; + public: //! Global MPI flag, default: 0=off @@ -437,7 +440,8 @@ public: SLGridSlab(int kmax, int nmax, int numz, double zmax, const std::string type, bool Verbose=false) - : SLGridSlab(kmax, nmax, numz, zmax, default_cache, "linear", type, Verbose) {} + : SLGridSlab(kmax, nmax, numz, zmax, + default_cache, "linear", type, Verbose) {} //! Destructor ~SLGridSlab(); diff --git a/utils/SL/slabchk.cc b/utils/SL/slabchk.cc index 3baf1589f..ee0db910b 100644 --- a/utils/SL/slabchk.cc +++ b/utils/SL/slabchk.cc @@ -89,7 +89,9 @@ main(int argc, char** argv) // Generate Sturm-Liouville grid // - auto ortho = std::make_shared(kmax, nmax, numz, zmax, "isothermal"); + auto ortho = std::make_shared(kmax, nmax, numz, zmax, + ".slgrid_slab_cache", + "linear", "isothermal"); LegeQuad lw(knots); From 5e3b42f27675495e3cf914d441f5d53c0d676af3 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 3 May 2026 16:53:06 -0400 Subject: [PATCH 45/51] Miswritten 'assign' members for Slab and Cube --- pyEXP/CoefWrappers.cc | 65 ++++++++++++++++++++++++++----------------- 1 file changed, 40 insertions(+), 25 deletions(-) diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index 1164ebf6b..8ed67fe1f 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -1000,8 +1000,13 @@ void CoefficientClasses(py::module &m) { py::class_, CoefStruct> (m, "SlabStruct") .def(py::init<>(), "Slab coefficient data structure object") - .def("assign", &SlabStruct::assign, - R"( + .def("assign", + [](CoefClasses::SlabStruct& A, py::array_t> mat) + { + auto M = make_tensor3>(mat); + A.assign(M); + }, + R"( Assign a coefficient matrix to CoefStruct. Parameters @@ -1022,8 +1027,13 @@ void CoefficientClasses(py::module &m) { py::class_, CoefStruct> (m, "CubeStruct") .def(py::init<>(), "Cube coefficient data structure object") - .def("assign", &CubeStruct::assign, - R"( + .def("assign", + [](CoefClasses::SlabStruct& A, py::array_t> mat) + { + auto M = make_tensor3>(mat); + A.assign(M); + }, + R"( Assign a coefficient matrix to CoefStruct. Parameters @@ -1045,13 +1055,13 @@ void CoefficientClasses(py::module &m) { (m, "TblStruct") .def(py::init<>(), "Multicolumn table data structure object") .def("assign", &TblStruct::assign, - R"( + R"( Assign a coefficient matrix to CoefStruct. Parameters ---------- mat : numpy.ndarray, complex - complex-valued NumPy array of table values + complex-valued NumPy array of table values Returns ------- @@ -1585,9 +1595,9 @@ void CoefficientClasses(py::module &m) { Parameters ---------- time : float - snapshot time corresponding to the the coefficient matrix - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient matrix + mat : numpy.ndarray + the new coefficient array. Returns ------- @@ -1660,9 +1670,9 @@ void CoefficientClasses(py::module &m) { Parameters ---------- time : float - snapshot time corresponding to the the coefficient matrix - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient matrix + mat : numpy.ndarray + the new coefficient array. Returns ------- @@ -1772,9 +1782,9 @@ void CoefficientClasses(py::module &m) { Parameters ---------- time : float - snapshot time corresponding to the the coefficient matrix - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient matrix + mat : numpy.ndarray + the new coefficient array. Returns ------- @@ -1856,9 +1866,9 @@ void CoefficientClasses(py::module &m) { Parameters ---------- time : float - snapshot time corresponding to the the coefficient matrix - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient matrix + mat : numpy.ndarray + the new coefficient array. Returns ------- @@ -1919,16 +1929,21 @@ void CoefficientClasses(py::module &m) { )", py::arg("time")) .def("setTensor", - &CoefClasses::SlabCoefs::setTensor, + [](CoefClasses::SlabCoefs& A, double time, + py::array_t> mat) + { + auto M = make_tensor3>(mat); + A.setTensor(time, M); + }, R"( Enter and/or rewrite the coefficient tensor at the provided time Parameters ---------- time : float - snapshot time corresponding to the the coefficient tensor - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient matrix + mat : numpy.ndarray of complex values + the new coefficient array. Returns ------- @@ -2020,9 +2035,9 @@ void CoefficientClasses(py::module &m) { Parameters ---------- time : float - snapshot time corresponding to the the coefficient tensor - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient tensor + mat : numpy.ndarray + the new coefficient array. Returns ------- From c0a233a1f84db51cac6cf639b3c2c982e216137d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 4 May 2026 10:22:59 -0400 Subject: [PATCH 46/51] Sign error in force table folding --- exputil/SLGridMP2.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 6b988c0cb..edcd34800 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -2496,7 +2496,7 @@ double SLGridSlab::get_dens(double x, int kx, int ky, int n, int which) int hold; int sign=1; - if (x<0 && 2*(n/2)!=n) sign=-1; + if (x<0 && 2*(n/2)==n) sign=-1; x = fabs(x); if (which) // Convert from z to x @@ -2557,10 +2557,10 @@ double SLGridSlab::get_force(double x, int kx, int ky, int n, int which) // Point 1: indx+1 return mM->d_xi_to_z(x)/dxi * ( - (p - 0.5)*table[kx][ky].ef(n, indx-1)*p0[indx-1] - -2.0*p*table[kx][ky].ef(n, indx)*p0[indx] - + (p + 0.5)*table[kx][ky].ef(n, indx+1)*p0[indx+1] - ) / sqrt(table[kx][ky].ev[n]) * sign; + (p - 0.5)*table[kx][ky].ef(n, indx-1)*p0[indx-1] + -2.0*p*table[kx][ky].ef(n, indx)*p0[indx] + + (p + 0.5)*table[kx][ky].ef(n, indx+1)*p0[indx+1] + ) / sqrt(table[kx][ky].ev[n]) * sign; } From ec6e8537195cf75a9e01d127cf354852f6bce3a0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 4 May 2026 10:23:55 -0400 Subject: [PATCH 47/51] Correct the mass prefactor for the biorthogonal projection --- src/cudaSlabSL.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cudaSlabSL.cu b/src/cudaSlabSL.cu index 838e51b25..1278e9289 100644 --- a/src/cudaSlabSL.cu +++ b/src/cudaSlabSL.cu @@ -320,7 +320,7 @@ __global__ void coefKernelSlab cudaParticle & p = P._v[I._v[npart]]; cuFP_t pos[3] = {p.pos[0], p.pos[1], p.pos[2]}; - cuFP_t mm = - p.mass * 2.0*slabDfac; + cuFP_t mm = - p.mass; // Biorthogonality gives a minus here used._v[i] = p.mass; // Mark particle as used From f871e2692d9fa511aec636c6b2b6507652f5bbd8 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 4 May 2026 10:24:59 -0400 Subject: [PATCH 48/51] Correct the biorthogonal prefactor and update ui for Cube and Slab coefficient creation in pyEXP --- expui/BiorthBasis.cc | 12 ++++++++---- src/SlabSL.cc | 4 +++- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 654d335f9..e8f983001 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3669,10 +3669,8 @@ namespace BasisClasses ortho->get_pot(zpot, z, iiy, iix); for (int iz=0; iz(); else if (name.compare("flatdisk") == 0) coef = std::make_shared(); + else if (name.compare("slab") == 0) + coef = std::make_shared(); else if (name.compare("cube") == 0) coef = std::make_shared(); else { @@ -4661,6 +4661,10 @@ namespace BasisClasses coefret = std::make_shared(); else if (name.compare("flatdisk") == 0) coefret = std::make_shared(); + else if (name.compare("slab") == 0) + coefret = std::make_shared(); + else if (name.compare("cube") == 0) + coefret = std::make_shared(); else { std::ostringstream sout; sout << "Basis::createCoefficients: basis <" << name << "> not recognized" diff --git a/src/SlabSL.cc b/src/SlabSL.cc index 9c9e88116..be54b6921 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -483,7 +483,9 @@ void * SlabSL::determine_coefficients_thread(void * arg) starty = exp(static_cast(nmaxy)*kfac*cC->Pos(i, 1)); double zz = cC->Pos(i, 2), ms = cC->Mass(i); - double mm = -4.0*M_PI * ms * adb; + // Biorthogonaity gives a minus sign + // here + double mm = - ms * adb; for (facx=startx, ix=0; ix Date: Mon, 4 May 2026 15:45:57 -0400 Subject: [PATCH 49/51] Fix 4*pi factors in slab normalization --- exputil/SLGridMP2.cc | 6 +++--- utils/SL/orthochk.cc | 41 ++++++++++++++++++++++------------------- 2 files changed, 25 insertions(+), 22 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index edcd34800..cfdc7986b 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1878,7 +1878,7 @@ class IsothermalSlab : public SlabModel double dens(double z) { double tmp = 1.0/cosh(0.5*z/SLGridSlab::H); - return 4.0*M_PI*0.25/SLGridSlab::H * tmp*tmp; + return 0.25/SLGridSlab::H * tmp*tmp; } }; @@ -1903,7 +1903,7 @@ class ConstantSlab : public SlabModel double dens(double z) { - return 4.0*M_PI / (2.0 * SLGridSlab::H); + return 1.0 / (2.0 * SLGridSlab::H); } }; @@ -1935,7 +1935,7 @@ class ParabolicSlab : public SlabModel { double h = SLGridSlab::H; double h2 = h*h; - return 4.0*M_PI * 3.0*(1.0 - z*z/h2)/(4.0*h); + return 3.0*(1.0 - z*z/h2)/(4.0*h); } }; diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index 86a3ccdbc..39775ff6d 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -122,6 +122,8 @@ main(int argc, char** argv) } + int ikx = max(IKX, IKY), iky = min(IKX, IKY); + //=================== // Get info //=================== @@ -178,9 +180,9 @@ main(int argc, char** argv) } else { x = orthoSL->z_to_xi(z); out << setw(15) << z - << setw(15) << orthoSL->get_pot(x, IKX, IKY, N) - << setw(15) << orthoSL->get_force(x, IKX, IKY, N) - << setw(15) << orthoSL->get_dens(x, IKX, IKY, N) + << setw(15) << orthoSL->get_pot (x, ikx, iky, N) + << setw(15) << orthoSL->get_force(x, ikx, iky, N) + << setw(15) << orthoSL->get_dens (x, ikx, iky, N) << endl; } } @@ -244,12 +246,12 @@ main(int argc, char** argv) Eigen::VectorXd pot0(NMAX), potN(NMAX), potP(NMAX), den0(NMAX); Eigen::VectorXd frcP(NMAX), frcN(NMAX); - orthoSL->get_pot (pot0, z , IKX, IKY); - orthoSL->get_pot (potN, z-h, IKX, IKY); - orthoSL->get_pot (potP, z+h, IKX, IKY); - orthoSL->get_force(frcN, z-0.5*h, IKX, IKY); - orthoSL->get_force(frcP, z+0.5*h, IKX, IKY); - orthoSL->get_dens (den0, z, IKX, IKY); + orthoSL->get_pot (pot0, z , ikx, iky); + orthoSL->get_pot (potN, z-h, ikx, iky); + orthoSL->get_pot (potP, z+h, ikx, iky); + orthoSL->get_force(frcN, z-0.5*h, ikx, iky); + orthoSL->get_force(frcP, z+0.5*h, ikx, iky); + orthoSL->get_dens (den0, z, ikx, iky); for (int n=0; npotl(i, i, z, pot); + ortho->dens(i, i, z, den); + } else { + orthoSL->get_pot (pot, z, ikx, iky); + orthoSL->get_dens(den, z, ikx, iky); + } + for (int n1=0; n1potl(n1, i, z) * ortho->dens(n2, i, z) * W; - break; - case SL: - ans(n1, n2) += - orthoSL->get_pot(z, IKX, IKY, n1) * orthoSL->get_dens(z, IKX, IKY, n2) * W; - break; - } + ans(n1, n2) += pot(n1) * den(n2) * W; } } } @@ -349,6 +351,7 @@ main(int argc, char** argv) } std::cout << "<" << N1 << "|" << N2 << "> = " << ans << std::endl; } + break; } case 4: { From ac21815216b09f5d97bccc96ac4746aea8d7f322 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 5 May 2026 10:21:14 -0400 Subject: [PATCH 50/51] Restored previous 4*pi factors which were correct --- exputil/SLGridMP2.cc | 15 ++++++++++++--- utils/SL/orthochk.cc | 4 ++-- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index cfdc7986b..57b6d62b4 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1878,7 +1878,10 @@ class IsothermalSlab : public SlabModel double dens(double z) { double tmp = 1.0/cosh(0.5*z/SLGridSlab::H); - return 0.25/SLGridSlab::H * tmp*tmp; + return 4.0*M_PI * 0.25/SLGridSlab::H * tmp*tmp; + // ^ + // | + // +--- The 4*pi factor simplifies the SL solution } }; @@ -1903,7 +1906,10 @@ class ConstantSlab : public SlabModel double dens(double z) { - return 1.0 / (2.0 * SLGridSlab::H); + return 4.0*M_PI / (2.0 * SLGridSlab::H); + // ^ + // | + // +--- The 4*pi factor simplifies the SL solution } }; @@ -1935,7 +1941,10 @@ class ParabolicSlab : public SlabModel { double h = SLGridSlab::H; double h2 = h*h; - return 3.0*(1.0 - z*z/h2)/(4.0*h); + return 4.0*M_PI * 3.0*(1.0 - z*z/h2)/(4.0*h); + // ^ + // | + // +--- The 4*pi factor simplifies the SL solution } }; diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index 39775ff6d..72295fffe 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -173,9 +173,9 @@ main(int argc, char** argv) if (Type == Trig) { x = ortho->r_to_rb(z); out << setw(15) << z - << setw(15) << ortho->potl(N, i, z) + << setw(15) << ortho->potl (N, i, z) << setw(15) << ortho->force(N, i, z) - << setw(15) << ortho->dens(N, i, z) + << setw(15) << ortho->dens (N, i, z) << endl; } else { x = orthoSL->z_to_xi(z); From 11beb18c4b20b20527e9e5144ce16bd83165034b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 5 May 2026 12:59:31 -0400 Subject: [PATCH 51/51] Fix 4*pi factors in slab normalization, 2nd try --- src/SlabSL.cc | 2 +- src/cudaSlabSL.cu | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/SlabSL.cc b/src/SlabSL.cc index be54b6921..48b5a219c 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -485,7 +485,7 @@ void * SlabSL::determine_coefficients_thread(void * arg) double zz = cC->Pos(i, 2), ms = cC->Mass(i); // Biorthogonaity gives a minus sign // here - double mm = - ms * adb; + double mm = - 4.0*M_PI * ms * adb; for (facx=startx, ix=0; ix