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 e39bfe825..7083fda0c 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; @@ -3169,37 +3169,59 @@ namespace BasisClasses // Sanity check // + if (tfinal == tinit){ + std::ostringstream sout; + sout << "BasicFactor::IntegrateOrbits: tinit cannot be equal to tfinal"; + 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 (h > 0. || tfinal <= tinit){ + std::ostringstream sout; + sout << "BasicFactor::IntegrateOrbits: tfinal should be larger than " + << "tinit when step size is positive"; + throw std::runtime_error(sout.str()); + } + if (nout < 2){ + std::ostringstream sout; + 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::cout << "BasicFactor::IntegrateOrbits: step size is too small or " - << "time interval is too large.\n"; - // Return empty data - // - return {Eigen::VectorXd(), Eigen::Tensor()}; + 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 ); + int numT = floor( (tfinal - tinit)/h ) + 1; // Compute output step + // Use numT and H for calculation // - nout = std::min(numT, nout); - double H = (tfinal - tinit)/nout; + numT = std::min(numT, nout); + double H = (tfinal - tinit)/numT; // 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 +3231,7 @@ namespace BasisClasses // Time array // - Eigen::VectorXd times(nout); + Eigen::VectorXd times(numT); // Do the work // @@ -3218,19 +3240,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 ${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/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index e0b72ab25..89951fd5d 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2004,7 +2004,7 @@ void BasisFactoryClasses(py::module &m) m.def("IntegrateOrbits", [](double tinit, double tfinal, double h, Eigen::MatrixXd ps, std::vector bfe, - BasisClasses::AccelFunc& func, int stride) + 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, h, ps, bfe, F, stride); + BasisClasses::IntegrateOrbits(tinit, tfinal, h, ps, bfe, F, nout); py::array_t ret = make_ndarray3(O); return std::tuple>(T, ret); 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",