including tfinal when integrating an orbit

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")

expui/BiorthBasis.cc

@@ -1031,7 +1031,7 @@ namespace BasisClasses
     nmax        = 18;
     mmax        = 6;
     mlim        = std::numeric_limits<int>::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<double>(std::numeric_limits<int>::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<float, 3>()};
+	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<int>(numT, nout);
-    double H = (tfinal - tinit)/nout;
+    numT = std::min<int>(numT, nout);
+    double H = (tfinal - tinit)/numT;
 
     // Return data
     //
     Eigen::Tensor<float, 3> 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<numT; s++) {
-      std::tie(tnow, ps) = OneStep(tnow, h, ps, accel, bfe, F);
-      if (tnow >= H*cnt-h*1.0e-8) {
-	times(cnt) = tnow;
-	for (int n=0; n<rows; n++)
-	  for (int k=0; k<6; k++) ret(n, k, cnt) = ps(n, k);
-	cnt += 1;
+    for (int cnt=1; cnt<numT; cnt++) {
+      std::tie(tnow, ps) = OneStep(tnow, H, ps, accel, bfe, F);
+      if (tnow >= H*cnt-H*1.0e-8) {
+        times(cnt) = tnow;
+        for (int n=0; n<rows; n++)
+          for (int k=0; k<6; k++) ret(n, k, cnt) = ps(n, k);
       }
     }
-
-    times(nout-1) = tnow;
-    for (int n=0; n<rows; n++)
-      for (int k=0; k<6; k++) ret(n, k, nout-1) = ps(n, k);
 
     return {times, ret};
   }

exputil/CMakeLists.txt

@@ -51,9 +51,9 @@ set(common_INCLUDE_DIRS $<INSTALL_INTERFACE:include>
   ${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)

pyEXP/BasisWrappers.cc

@@ -2004,15 +2004,15 @@ void BasisFactoryClasses(py::module &m)
   m.def("IntegrateOrbits", 
 	[](double tinit, double tfinal, double h, Eigen::MatrixXd ps,
 	   std::vector<BasisClasses::BasisCoef> bfe,
-	   BasisClasses::AccelFunc& func, int stride)
+	   BasisClasses::AccelFunc& func, int nout)
 	{
 	  Eigen::VectorXd T;
 	  Eigen::Tensor<float, 3> O;
 
 	  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<float> ret = make_ndarray3<float>(O);
 	  return std::tuple<Eigen::VectorXd, py::array_t<float>>(T, ret);

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

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<string>(ctype)->default_value("Log"))
     ("LMAXFID", "Maximum angular order for spherical basis in adaptive construction of the cylindrical basis",
-     cxxopts::value<int>(LMAXFID)->default_value("72"))
+     cxxopts::value<int>(LMAXFID)->default_value("32"))
     ("NMAXFID", "Maximum radial order for the spherical basis in adapative construction of the cylindrical basis",
      cxxopts::value<int>(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<int>(NCYLODD)->default_value("6"))
     ("NMAX", "Total number of basis functions per harmonic order",
-     cxxopts::value<int>(NMAX)->default_value("18"))
+     cxxopts::value<int>(NMAX)->default_value("20"))
     ("VFLAG", "Diagnostic flag for EmpCylSL",
      cxxopts::value<int>(VFLAG)->default_value("31"))
     ("expcond", "Use analytic target density rather than particle distribution",

utils/ICs/initial.cc

@@ -398,9 +398,9 @@ main(int ac, char **av)
     ("LMAX", "Harmonic order for halo expansion",
      cxxopts::value<int>(LMAX)->default_value("18"))
     ("LMAXFID", "Harmonic order for EOF spherical model",
-     cxxopts::value<int>(LMAXFID)->default_value("72"))
+     cxxopts::value<int>(LMAXFID)->default_value("48"))
     ("NMAXFID", "Radial order for EOF spherical model",
-     cxxopts::value<int>(NMAXFID)->default_value("64"))
+     cxxopts::value<int>(NMAXFID)->default_value("48"))
     ("MMAX", "Aximuthal order for Cylindrical expansion",
      cxxopts::value<int>(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<int>(NOUT)->default_value("18"))
     ("NODD", "Number of vertically odd terms in cylindrical expansion",
-     cxxopts::value<int>(NODD)->default_value("9"))
+     cxxopts::value<int>(NODD)->default_value("6"))
     ("NMAXH", "Number of radial terms for spherical expansion",
      cxxopts::value<int>(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<double>(SCSPH)->default_value("1.0"))
     ("RSPHSL", "Maximum halo expansion radius",
-     cxxopts::value<double>(RSPHSL)->default_value("1.95"))
+     cxxopts::value<double>(RSPHSL)->default_value("47.5"))
     ("ASCALE", "Radial scale length for disk basis construction",
      cxxopts::value<double>(ASCALE)->default_value("1.0"))
     ("ASHIFT", "Fraction of scale length for shift in conditioning function",

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;

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<int>(N)->default_value("10000"))
     ("nmax",       "maximum number of radial harmonics in the expansion",
-     cxxopts::value<int>(nmax)->default_value("18"))
+     cxxopts::value<int>(nmax)->default_value("10"))
     ("r,rmin",     "minimum radius for the SL grid",
      cxxopts::value<double>(rmin)->default_value("0.0001"))
     ("R,rmax",     "maximum radius for the SL grid",

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<int>(Lmax)->default_value("2"))
    ("nmax", "maximum number of radial harmonics in the expansion",
-     cxxopts::value<int>(nmax)->default_value("18"))
+     cxxopts::value<int>(nmax)->default_value("10"))
    ("numr", "radial knots for the SL grid",
      cxxopts::value<int>(numr)->default_value("1000"))
    ("rmin", "minimum radius for the SL grid",

utils/SL/slabchk.cc

@@ -29,7 +29,7 @@ main(int argc, char** argv)
     ("K,kmax", "maximum order of in-plane harmonics",
      cxxopts::value<int>(kmax)->default_value("4"))
     ("N,nmax", "maximum number of vertical harmonics",
-     cxxopts::value<int>(nmax)->default_value("18"))
+     cxxopts::value<int>(nmax)->default_value("10"))
     ("n,numz", "size of vertical grid",
      cxxopts::value<int>(numz)->default_value("1000"))
     ("Z,zmax", "maximum extent of vertical grid",

utils/SL/slshift.cc

@@ -353,7 +353,7 @@ main(int argc, char** argv)
    ("Lmax", "maximum number of angular harmonics in the expansion",
      cxxopts::value<int>(Lmax)->default_value("2"))
    ("nmax", "maximum number of radial harmonics in the expansion",
-     cxxopts::value<int>(nmax)->default_value("18"))
+     cxxopts::value<int>(nmax)->default_value("10"))
    ("numr", "radial knots in the shift operator",
      cxxopts::value<int>(numr)->default_value("1000"))
    ("rmin", "minimum radius for the shift operator",