From 2336b9af5777c754e12596b0f43827c1ca2a3d31 Mon Sep 17 00:00:00 2001 From: M1ssing-N0 <100515873+M1ssing-N0@users.noreply.github.com> Date: Sun, 15 Dec 2024 21:25:38 +0800 Subject: [PATCH 1/4] Editing IntegrateOrbits --- expui/BiorthBasis.H | 5 ++- expui/BiorthBasis.cc | 86 ++++++++++++++++++++++++++++-------------- pyEXP/BasisWrappers.cc | 12 +++--- 3 files changed, 66 insertions(+), 37 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index a4af27c1a..9b8c4e744 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -3,6 +3,7 @@ #include #include +#include #include #include // For 3d rectangular grids @@ -1033,9 +1034,9 @@ namespace BasisClasses std::tuple> - IntegrateOrbits (double tinit, double tfinal, double h, + IntegrateOrbits (double tinit, double tfinal, Eigen::MatrixXd ps, std::vector bfe, - AccelFunctor F, int nout=std::numeric_limits::max()); + AccelFunctor F, std::optional h=std::nullopt, std::optional nout=std::nullopt); using BiorthBasisPtr = std::shared_ptr; } diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index f69bb0706..442a71e2b 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3147,9 +3147,9 @@ namespace BasisClasses std::tuple> IntegrateOrbits - (double tinit, double tfinal, double h, + (double tinit, double tfinal, Eigen::MatrixXd ps, std::vector bfe, AccelFunctor F, - int nout) + std::optional h, std::optional nout) { int rows = ps.rows(); int cols = ps.cols(); @@ -3169,37 +3169,70 @@ namespace BasisClasses // Sanity check // - if ( (tfinal - tinit)/h > + if (h.has_value()){ + if (*h <= 0.){ + std::ostringstream sout; + sout << "BasicFactor::IntegrateOrbits: unreasonal input"; + throw std::runtime_error(sout.str()); + } + } + if (nout.has_value()){ + if (*nout < 2){ + std::ostringstream sout; + sout << "BasicFactor::IntegrateOrbits: unreasonal input"; + throw std::runtime_error(sout.str()); + } + } + if ( (tfinal - tinit)/ *h > static_cast(std::numeric_limits::max()) ) { - std::cout << "BasicFactor::IntegrateOrbits: step size is too small or " + std::ostringstream sout; + sout << "BasicFactor::IntegrateOrbits: step size is too small or " << "time interval is too large.\n"; - // Return empty data - // - return {Eigen::VectorXd(), Eigen::Tensor()}; + throw std::runtime_error(sout.str()); } + if (tinit < 0 || tfinal < tinit){ + std::ostringstream sout; + sout << "BasicFactor::IntegrateOrbits: unreasonal input"; + throw std::runtime_error(sout.str()); + } - // Number of steps - // - int numT = floor( (tfinal - tinit)/h ); - - // Compute output step - // - nout = std::min(numT, nout); - double H = (tfinal - tinit)/nout; + int numT = 0; + double H = 0.0; + if (h.has_value() && nout.has_value()){ + numT = floor( (tfinal - tinit)/ *h ) + 1; + if (std::abs(numT - *nout) > 0){ + std::ostringstream sout; + sout << "BasicFactor::IntegrateOrbits: inconsistent step number/size input"; + throw std::runtime_error(sout.str()); + } else { + H = *h; + } + } else if (h.has_value()) { + numT = floor( (tfinal - tinit)/ *h ) + 1; + H = *h; + } else if (nout.has_value()) { + numT = *nout; + H = (tfinal - tinit)/ (*nout - 1); + } else { + std::ostringstream sout; + sout << "BasicFactor::IntegrateOrbits: no step number/size input"; + throw std::runtime_error(sout.str()); + } + // Use numT and H for calculation // Return data // Eigen::Tensor ret; try { - ret.resize(rows, 6, nout); + ret.resize(rows, 6, numT); } catch (const std::bad_alloc& e) { std::cout << "BasicFactor::IntegrateOrbits: memory allocation failed: " << e.what() << std::endl << "Your requested number of orbits and time steps requires " - << floor(8.0*rows*6*nout/1e9)+1 << " GB free memory" + << floor(8.0*rows*6*numT/1e9)+1 << " GB free memory" << std::endl; // Return empty data @@ -3209,7 +3242,7 @@ namespace BasisClasses // Time array // - Eigen::VectorXd times(nout); + Eigen::VectorXd times(numT); // Do the work // @@ -3218,19 +3251,14 @@ namespace BasisClasses for (int k=0; k<6; k++) ret(n, k, 0) = ps(n, k); double tnow = tinit; - for (int s=1, cnt=1; s= H*cnt-h*1.0e-8) { - times(cnt) = tnow; - for (int n=0; n= H*cnt-H*1.0e-8) { + times(cnt) = tnow; + for (int n=0; n bfe, - BasisClasses::AccelFunc& func, int stride) + BasisClasses::AccelFunc& func, std::optional h, std::optional nout) { Eigen::VectorXd T; Eigen::Tensor O; @@ -2012,9 +2012,9 @@ void BasisFactoryClasses(py::module &m) AccelFunctor F = [&func](double t, Eigen::MatrixXd& ps, Eigen::MatrixXd& accel, BasisCoef mod)->Eigen::MatrixXd& { return func.F(t, ps, accel, mod);}; std::tie(T, O) = - BasisClasses::IntegrateOrbits(tinit, tfinal, h, ps, bfe, F, stride); + BasisClasses::IntegrateOrbits(tinit, tfinal, ps, bfe, F, h, nout); - py::array_t ret = make_ndarray3(O); + py::array_t ret = make_ndarray(O); return std::tuple>(T, ret); }, R"( @@ -2049,7 +2049,7 @@ void BasisFactoryClasses(py::module &m) tuple(numpy.array, numpy.ndarray) time and phase-space arrays )", - py::arg("tinit"), py::arg("tfinal"), py::arg("h"), + py::arg("tinit"), py::arg("tfinal"), py::arg("h")=std::nullopt, py::arg("ps"), py::arg("basiscoef"), py::arg("func"), - py::arg("nout")=std::numeric_limits::max()); + py::arg("nout")=std::nullopt); } From a47595d0bf14440ae90be8083178a0f80486af9a Mon Sep 17 00:00:00 2001 From: M1ssing-N0 <100515873+M1ssing-N0@users.noreply.github.com> Date: Tue, 17 Dec 2024 00:57:56 +0800 Subject: [PATCH 2/4] Resolve version difference Commit 12176dc: Change names of Eigen::Tensor-->ndarray generators for clarity --- pyEXP/BasisWrappers.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 7c3431315..e1a5f0a2c 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2014,7 +2014,7 @@ void BasisFactoryClasses(py::module &m) std::tie(T, O) = BasisClasses::IntegrateOrbits(tinit, tfinal, ps, bfe, F, h, nout); - py::array_t ret = make_ndarray(O); + py::array_t ret = make_ndarray3(O); return std::tuple>(T, ret); }, R"( From 9a9b873453829054304bdc369a6e2707548596e2 Mon Sep 17 00:00:00 2001 From: M1ssing-N0 <100515873+M1ssing-N0@users.noreply.github.com> Date: Wed, 18 Dec 2024 16:22:13 +0800 Subject: [PATCH 3/4] Revert "Merge branch 'EXP-code:main' into main" This reverts commit c0bfbcd430886dff94e60b50c7d1ed5de11a79cc, reversing changes made to a47595d0bf14440ae90be8083178a0f80486af9a. --- CMakeLists.txt | 3 --- expui/BiorthBasis.cc | 2 +- exputil/CMakeLists.txt | 6 +++--- src/expand.H | 2 +- utils/ICs/cylcache.cc | 4 ++-- utils/ICs/initial.cc | 8 ++++---- utils/SL/diskpot.cc | 2 +- utils/SL/oftest.cc | 2 +- utils/SL/qtest.cc | 2 +- utils/SL/slabchk.cc | 2 +- utils/SL/slshift.cc | 2 +- 11 files changed, 16 insertions(+), 19 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 70a4dc531..1d22f3e97 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,9 +25,6 @@ endif() # Required compiler features add_compile_options(-D_REENTRANT) -# Bake in library paths (esp. useful for HPC sites with modules) -set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) - # Check and enforce that we are a git repository. Necessary for # submodules to work correctly. if(EXISTS "${PROJECT_SOURCE_DIR}/.git") diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 7de4fb202..442a71e2b 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1031,7 +1031,7 @@ namespace BasisClasses nmax = 18; mmax = 6; mlim = std::numeric_limits::max(); - lmaxfid = 72; + lmaxfid = 128; nmaxfid = 64; ncylnx = 256; ncylny = 128; diff --git a/exputil/CMakeLists.txt b/exputil/CMakeLists.txt index c785e5817..9ccb58eca 100644 --- a/exputil/CMakeLists.txt +++ b/exputil/CMakeLists.txt @@ -51,9 +51,9 @@ set(common_INCLUDE_DIRS $ ${DEP_INC} ${EIGEN3_INCLUDE_DIR} ${HDF5_INCLUDE_DIRS} ${FFTW_INCLUDE_DIRS}) -set(common_LINKLIB ${DEP_LIB} OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp - ${VTK_LIBRARIES} ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES} - ${HDF5_HL_LIBRARIES} ${FFTW_DOUBLE_LIB}) +set(common_LINKLIB ${DEP_LIB} OpenMP::OpenMP_CXX MPI::MPI_CXX + yaml-cpp ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES} + ${FFTW_DOUBLE_LIB}) if(ENABLE_CUDA) list(APPEND common_LINKLIB CUDA::toolkit CUDA::cudart) diff --git a/src/expand.H b/src/expand.H index 220b9187c..b8c22d892 100644 --- a/src/expand.H +++ b/src/expand.H @@ -26,7 +26,7 @@ KEYWORD EXPLANATION DEFAULT ---------- ---------------------------- -------- lmax = maximum harmonic order 4 - nmax = maximum radial order 18 + nmax = maxinum radial order 10 nlog = interval between log output 10 nlist = interval between p-s output 50 nsteps = total number of steps to run 500 diff --git a/utils/ICs/cylcache.cc b/utils/ICs/cylcache.cc index fd78cdbf8..9a2422c6f 100644 --- a/utils/ICs/cylcache.cc +++ b/utils/ICs/cylcache.cc @@ -327,7 +327,7 @@ main(int ac, char **av) ("ctype", "DiskHalo radial coordinate scaling type (one of: Linear, Log, Rat)", cxxopts::value(ctype)->default_value("Log")) ("LMAXFID", "Maximum angular order for spherical basis in adaptive construction of the cylindrical basis", - cxxopts::value(LMAXFID)->default_value("72")) + cxxopts::value(LMAXFID)->default_value("32")) ("NMAXFID", "Maximum radial order for the spherical basis in adapative construction of the cylindrical basis", cxxopts::value(NMAXFID)->default_value("32")) ("MMAX", "Maximum azimuthal order for the cylindrical basis", @@ -341,7 +341,7 @@ main(int ac, char **av) ("NCYLODD", "Number of vertically odd basis functions per harmonic order", cxxopts::value(NCYLODD)->default_value("6")) ("NMAX", "Total number of basis functions per harmonic order", - cxxopts::value(NMAX)->default_value("18")) + cxxopts::value(NMAX)->default_value("20")) ("VFLAG", "Diagnostic flag for EmpCylSL", cxxopts::value(VFLAG)->default_value("31")) ("expcond", "Use analytic target density rather than particle distribution", diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index 5a8825aa6..354733a87 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -398,9 +398,9 @@ main(int ac, char **av) ("LMAX", "Harmonic order for halo expansion", cxxopts::value(LMAX)->default_value("18")) ("LMAXFID", "Harmonic order for EOF spherical model", - cxxopts::value(LMAXFID)->default_value("72")) + cxxopts::value(LMAXFID)->default_value("48")) ("NMAXFID", "Radial order for EOF spherical model", - cxxopts::value(NMAXFID)->default_value("64")) + cxxopts::value(NMAXFID)->default_value("48")) ("MMAX", "Aximuthal order for Cylindrical expansion", cxxopts::value(MMAX)->default_value("6")) ("NUMX", "Number of knots in radial dimension of meridional grid", @@ -434,7 +434,7 @@ main(int ac, char **av) ("NOUT", "Number of radial terms in diagnostic basis file for cylinder", cxxopts::value(NOUT)->default_value("18")) ("NODD", "Number of vertically odd terms in cylindrical expansion", - cxxopts::value(NODD)->default_value("9")) + cxxopts::value(NODD)->default_value("6")) ("NMAXH", "Number of radial terms for spherical expansion", cxxopts::value(NMAXH)->default_value("18")) ("NMAXD", "Number of radial terms for cylindrical expansion", @@ -518,7 +518,7 @@ main(int ac, char **av) ("SCSPH", "Scale for Spherical SL coordinate mapping", cxxopts::value(SCSPH)->default_value("1.0")) ("RSPHSL", "Maximum halo expansion radius", - cxxopts::value(RSPHSL)->default_value("1.95")) + cxxopts::value(RSPHSL)->default_value("47.5")) ("ASCALE", "Radial scale length for disk basis construction", cxxopts::value(ASCALE)->default_value("1.0")) ("ASHIFT", "Fraction of scale length for shift in conditioning function", diff --git a/utils/SL/diskpot.cc b/utils/SL/diskpot.cc index c8cfbb3af..d2c930b14 100644 --- a/utils/SL/diskpot.cc +++ b/utils/SL/diskpot.cc @@ -57,7 +57,7 @@ main(int argc, char** argv) int cmap = 0; double scale = 1.0; - int Lmax=16, Nmax=18; + int Lmax=16, Nmax=10; int numr=10000; double rmin=0.0001, rmax=1.0; double delr=0.01, xmax=1.0, zmax=0.1; diff --git a/utils/SL/oftest.cc b/utils/SL/oftest.cc index 1426d3ae3..1c16b36f5 100644 --- a/utils/SL/oftest.cc +++ b/utils/SL/oftest.cc @@ -73,7 +73,7 @@ int main(int argc, char** argv) ("N,mc", "number of particles for Monte Carlo realization", cxxopts::value(N)->default_value("10000")) ("nmax", "maximum number of radial harmonics in the expansion", - cxxopts::value(nmax)->default_value("18")) + cxxopts::value(nmax)->default_value("10")) ("r,rmin", "minimum radius for the SL grid", cxxopts::value(rmin)->default_value("0.0001")) ("R,rmax", "maximum radius for the SL grid", diff --git a/utils/SL/qtest.cc b/utils/SL/qtest.cc index 352777d6e..1a3e1cf40 100644 --- a/utils/SL/qtest.cc +++ b/utils/SL/qtest.cc @@ -35,7 +35,7 @@ int main(int argc, char** argv) ("Lmax", "maximum number of angular harmonics in the expansion", cxxopts::value(Lmax)->default_value("2")) ("nmax", "maximum number of radial harmonics in the expansion", - cxxopts::value(nmax)->default_value("18")) + cxxopts::value(nmax)->default_value("10")) ("numr", "radial knots for the SL grid", cxxopts::value(numr)->default_value("1000")) ("rmin", "minimum radius for the SL grid", diff --git a/utils/SL/slabchk.cc b/utils/SL/slabchk.cc index c5d36eaa9..de3c045ff 100644 --- a/utils/SL/slabchk.cc +++ b/utils/SL/slabchk.cc @@ -29,7 +29,7 @@ main(int argc, char** argv) ("K,kmax", "maximum order of in-plane harmonics", cxxopts::value(kmax)->default_value("4")) ("N,nmax", "maximum number of vertical harmonics", - cxxopts::value(nmax)->default_value("18")) + cxxopts::value(nmax)->default_value("10")) ("n,numz", "size of vertical grid", cxxopts::value(numz)->default_value("1000")) ("Z,zmax", "maximum extent of vertical grid", diff --git a/utils/SL/slshift.cc b/utils/SL/slshift.cc index fabe9d7e4..9356152eb 100644 --- a/utils/SL/slshift.cc +++ b/utils/SL/slshift.cc @@ -353,7 +353,7 @@ main(int argc, char** argv) ("Lmax", "maximum number of angular harmonics in the expansion", cxxopts::value(Lmax)->default_value("2")) ("nmax", "maximum number of radial harmonics in the expansion", - cxxopts::value(nmax)->default_value("18")) + cxxopts::value(nmax)->default_value("10")) ("numr", "radial knots in the shift operator", cxxopts::value(numr)->default_value("1000")) ("rmin", "minimum radius for the shift operator", From 1205a0aef7f0c4e617331c592baae0cfcafcb6d3 Mon Sep 17 00:00:00 2001 From: M1ssing-N0 <100515873+M1ssing-N0@users.noreply.github.com> Date: Wed, 18 Dec 2024 17:03:35 +0800 Subject: [PATCH 4/4] Code clean up This partially reverts commit 2336b9af5777c754e12596b0f43827c1ca2a3d31. Updated sanity check to allow negative tstep --- expui/BiorthBasis.H | 5 ++- expui/BiorthBasis.cc | 77 ++++++++++++++++++------------------------ pyEXP/BasisWrappers.cc | 10 +++--- 3 files changed, 40 insertions(+), 52 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 9b8c4e744..a4af27c1a 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -3,7 +3,6 @@ #include #include -#include #include #include // For 3d rectangular grids @@ -1034,9 +1033,9 @@ namespace BasisClasses std::tuple> - IntegrateOrbits (double tinit, double tfinal, + IntegrateOrbits (double tinit, double tfinal, double h, Eigen::MatrixXd ps, std::vector bfe, - AccelFunctor F, std::optional h=std::nullopt, std::optional nout=std::nullopt); + AccelFunctor F, int nout=std::numeric_limits::max()); using BiorthBasisPtr = std::shared_ptr; } diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 442a71e2b..7083fda0c 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3147,9 +3147,9 @@ namespace BasisClasses std::tuple> IntegrateOrbits - (double tinit, double tfinal, + (double tinit, double tfinal, double h, Eigen::MatrixXd ps, std::vector bfe, AccelFunctor F, - std::optional h, std::optional nout) + int nout) { int rows = ps.rows(); int cols = ps.cols(); @@ -3169,58 +3169,47 @@ namespace BasisClasses // Sanity check // - if (h.has_value()){ - if (*h <= 0.){ - std::ostringstream sout; - sout << "BasicFactor::IntegrateOrbits: unreasonal input"; - throw std::runtime_error(sout.str()); - } + if (tfinal == tinit){ + std::ostringstream sout; + sout << "BasicFactor::IntegrateOrbits: tinit cannot be equal to tfinal"; + throw std::runtime_error(sout.str()); } - if (nout.has_value()){ - if (*nout < 2){ - std::ostringstream sout; - sout << "BasicFactor::IntegrateOrbits: unreasonal input"; - throw std::runtime_error(sout.str()); - } + if (h < 0. || tfinal >= tinit){ + std::ostringstream sout; + sout << "BasicFactor::IntegrateOrbits: tfinal should be smaller than " + << "tinit when step size is negative"; + throw std::runtime_error(sout.str()); } - if ( (tfinal - tinit)/ *h > - static_cast(std::numeric_limits::max()) ) - { - std::ostringstream sout; - sout << "BasicFactor::IntegrateOrbits: step size is too small or " - << "time interval is too large.\n"; - throw std::runtime_error(sout.str()); - } - if (tinit < 0 || tfinal < tinit){ + if (h > 0. || tfinal <= tinit){ std::ostringstream sout; - sout << "BasicFactor::IntegrateOrbits: unreasonal input"; + sout << "BasicFactor::IntegrateOrbits: tfinal should be larger than " + << "tinit when step size is positive"; throw std::runtime_error(sout.str()); } - - int numT = 0; - double H = 0.0; - if (h.has_value() && nout.has_value()){ - numT = floor( (tfinal - tinit)/ *h ) + 1; - if (std::abs(numT - *nout) > 0){ - std::ostringstream sout; - sout << "BasicFactor::IntegrateOrbits: inconsistent step number/size input"; - throw std::runtime_error(sout.str()); - } else { - H = *h; - } - } else if (h.has_value()) { - numT = floor( (tfinal - tinit)/ *h ) + 1; - H = *h; - } else if (nout.has_value()) { - numT = *nout; - H = (tfinal - tinit)/ (*nout - 1); - } else { + if (nout < 2){ std::ostringstream sout; - sout << "BasicFactor::IntegrateOrbits: no step number/size input"; + sout << "BasicFactor::IntegrateOrbits: nout must be larger than 2"; throw std::runtime_error(sout.str()); } + if ( (tfinal - tinit)/h > + static_cast(std::numeric_limits::max()) ) + { + std::ostringstream sout; + sout << "BasicFactor::IntegrateOrbits: step size is too small or " + << "time interval is too large.\n"; + throw std::runtime_error(sout.str()); + } + + // Number of steps + // + int numT = floor( (tfinal - tinit)/h ) + 1; + // Compute output step // Use numT and H for calculation + // + numT = std::min(numT, nout); + double H = (tfinal - tinit)/numT; + // Return data // Eigen::Tensor ret; diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index e1a5f0a2c..89951fd5d 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2002,9 +2002,9 @@ void BasisFactoryClasses(py::module &m) )", py::arg("time"), py::arg("mod")); m.def("IntegrateOrbits", - [](double tinit, double tfinal, Eigen::MatrixXd ps, + [](double tinit, double tfinal, double h, Eigen::MatrixXd ps, std::vector bfe, - BasisClasses::AccelFunc& func, std::optional h, std::optional nout) + BasisClasses::AccelFunc& func, int nout) { Eigen::VectorXd T; Eigen::Tensor O; @@ -2012,7 +2012,7 @@ void BasisFactoryClasses(py::module &m) AccelFunctor F = [&func](double t, Eigen::MatrixXd& ps, Eigen::MatrixXd& accel, BasisCoef mod)->Eigen::MatrixXd& { return func.F(t, ps, accel, mod);}; std::tie(T, O) = - BasisClasses::IntegrateOrbits(tinit, tfinal, ps, bfe, F, h, nout); + BasisClasses::IntegrateOrbits(tinit, tfinal, h, ps, bfe, F, nout); py::array_t ret = make_ndarray3(O); return std::tuple>(T, ret); @@ -2049,7 +2049,7 @@ void BasisFactoryClasses(py::module &m) tuple(numpy.array, numpy.ndarray) time and phase-space arrays )", - py::arg("tinit"), py::arg("tfinal"), py::arg("h")=std::nullopt, + py::arg("tinit"), py::arg("tfinal"), py::arg("h"), py::arg("ps"), py::arg("basiscoef"), py::arg("func"), - py::arg("nout")=std::nullopt); + py::arg("nout")=std::numeric_limits::max()); }