Merge pull request #66 from EXP-code/pyEXPfix

pyEXP bug fix

CMakeLists.txt

@@ -63,7 +63,6 @@ set(CMAKE_FIND_PACKAGE_SORT_DIRECTION DEC)
 # Package support
 
 find_package(MPI REQUIRED COMPONENTS C CXX)
-find_package(Boost COMPONENTS serialization)
 find_package(OpenMP)
 find_package(FFTW)
 find_package(HDF5 COMPONENTS C CXX HL REQUIRED)
@@ -151,9 +150,6 @@ if(PNG_FOUND AND ENABLE_PNG)
   set(HAVE_LIBPNGPP TRUE)
 endif()
 if(ENABLE_DSMC)
-  if(NOT BOOST_FOUND)
-    message(SEND_ERROR "You need Boost to compile DSMC")
-  endif()
   set(DSMC_ENABLED 1)
 endif()
 if(ENABLE_CUDA_SINGLE)
@@ -230,7 +226,7 @@ add_subdirectory(extern/HighFive)
 
 # Configure the remaining native subdirectories
 add_subdirectory(exputil)
-add_subdirectory(coefs)
+add_subdirectory(expui)
 if (ENABLE_NBODY)
   add_subdirectory(src)
 endif()
@@ -241,6 +237,10 @@ endif()
 
 add_subdirectory(extern/user-modules)
 
+if (ENABLE_DSMC)
+  add_subdirectory(extern/DSMC/src)
+endif()
+
 # Build the tests; set ENABLE_TEST=OFF to disable
 if(ENABLE_TESTS)
   include(CTest)

coefs/BasisFactory.H → expui/BasisFactory.H

@@ -200,13 +200,15 @@ namespace BasisClasses
 
     //! Accumulate coefficient contributions from arrays
     virtual void
-    addFromArray(Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin) = 0;
+    addFromArray(Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin,
+		 bool posvelrows) = 0;
 
     //! Generate coeffients from an array and optional center location
     //! for the expansion
     CoefClasses::CoefStrPtr createFromArray
     (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0,
-     std::vector<double> center={0.0, 0.0, 0.0}, bool roundrobin=true);
+     std::vector<double> center={0.0, 0.0, 0.0},
+     bool roundrobin=true, bool posvelrows=false);
 
     //! Create and the coefficients from the array accumulation with the
     //! provided time value

coefs/BasisFactory.cc → expui/BasisFactory.cc

@@ -267,10 +267,10 @@ namespace BasisClasses
   //
   CoefClasses::CoefStrPtr Basis::createFromArray
   (Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector<double> ctr,
-   bool roundrobin)
+   bool roundrobin, bool posvelrows)
   {
     initFromArray(ctr);
-    addFromArray(m, p, roundrobin);
+    addFromArray(m, p, roundrobin, posvelrows);
     return makeFromArray(time);
   }
 

coefs/BiorthBasis.H → expui/BiorthBasis.H

@@ -96,7 +96,8 @@ namespace BasisClasses
     //! for the expansion
     CoefClasses::CoefStrPtr createFromArray
     (Eigen::VectorXd& m, RowMatrixXd& p, double time=0.0,
-     std::vector<double> center={0.0, 0.0, 0.0}, bool roundrobin=true);
+     std::vector<double> center={0.0, 0.0, 0.0},
+     bool roundrobin=true, bool posvelrows=false);
 
     //! Generate coeffients from an array and optional center location
     //! for the expansion using multiple particle partitions
@@ -110,7 +111,8 @@ namespace BasisClasses
     //! Initialize accumulating coefficients from arrays.  This is
     //! called once to initialize the accumulation.
     void addFromArray
-    (Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin=true);
+    (Eigen::VectorXd& m, RowMatrixXd& p,
+     bool roundrobin=true, bool posvelrows=false);
 
     //! Create and the coefficients from the array accumulation with the
     //! provided time value

coefs/BiorthBasis.cc → expui/BiorthBasis.cc

@@ -90,13 +90,13 @@ namespace BasisClasses
   }
 
   SphericalSL::SphericalSL(const YAML::Node& CONF) :
-    BiorthBasis(CONF, "SphericalSL")
+    BiorthBasis(CONF, "sphereSL")
   {
     initialize();
   }
 
   SphericalSL::SphericalSL(const std::string& confstr) :
-    BiorthBasis(confstr, "SphericalSL")
+    BiorthBasis(confstr, "sphereSL")
   {
     initialize();
   }
@@ -279,7 +279,7 @@ namespace BasisClasses
     int ldim = (lmax+1)*(lmax+2)/2;
 
     // Allocate the coefficient storage
-    cf->store.resize((lmax+1)*(lmax+2)/2, nmax);
+    cf->store.resize((lmax+1)*(lmax+2)/2*nmax);
 
     // Make the coefficient map
     cf->coefs = std::make_shared<CoefClasses::SphStruct::coefType>
@@ -826,13 +826,13 @@ namespace BasisClasses
   };
 
   Cylindrical::Cylindrical(const YAML::Node& CONF) :
-    BiorthBasis(CONF, "Cylindrical")
+    BiorthBasis(CONF, "cylinder")
   {
     initialize();
   }
 
   Cylindrical::Cylindrical(const std::string& confstr) :
-    BiorthBasis(confstr, "Cylindrical")
+    BiorthBasis(confstr, "cylinder")
   {
     initialize();
   }
@@ -1147,8 +1147,7 @@ namespace BasisClasses
 			       << mtype << ">, valid types are: "
 			       << "Exponential, Gaussian, Plummer, Power "
 			       << "(not case sensitive)" << std::endl;
-	if (use_mpi) MPI_Finalize();
-	return;
+	throw std::runtime_error("Cylindrical:initialize: EmpCylSL bad parameter");
       }
 
       // Convert dtype string to lower case
@@ -1172,8 +1171,7 @@ namespace BasisClasses
 	  for (auto v : dtlookup) std::cout << v.first << " ";
 	  std::cout << std::endl;
 	}
-	if (use_mpi) MPI_Finalize();
-	return;
+	throw std::runtime_error("Cylindrical::initialize: invalid DiskType");
       }
 
       // Use these user models to deproject for the EOF spherical basis
@@ -1303,7 +1301,14 @@ namespace BasisClasses
 
     Eigen::VectorXd cos1(nmax), sin1(nmax);
 
-    cf->store((mmax+1)*nmax);
+    // Initialize the values
+    cos1.setZero();
+    sin1.setZero();
+
+    // Allocate the coefficient storage
+    cf->store.resize((mmax+1)*nmax);
+
+    // Create a new instance
     cf->coefs = std::make_shared<CoefClasses::CylStruct::coefType>
       (cf->store.data(), mmax+1, nmax);
 
@@ -1428,12 +1433,14 @@ namespace BasisClasses
     "cachename"
   };
 
-  FlatDisk::FlatDisk(const YAML::Node& CONF) : BiorthBasis(CONF)
+  FlatDisk::FlatDisk(const YAML::Node& CONF) :
+    BiorthBasis(CONF, "flatdisk")
   {
     initialize();
   }
 
-  FlatDisk::FlatDisk(const std::string& confstr) : BiorthBasis(confstr)
+  FlatDisk::FlatDisk(const std::string& confstr) :
+    BiorthBasis(confstr, "flatdisk")
   {
     initialize();
   }
@@ -1895,12 +1902,12 @@ namespace BasisClasses
     "method"
   };
 
-  Cube::Cube(const YAML::Node& CONF) : BiorthBasis(CONF, "Cube")
+  Cube::Cube(const YAML::Node& CONF) : BiorthBasis(CONF, "cube")
   {
     initialize();
   }
 
-  Cube::Cube(const std::string& confstr) : BiorthBasis(confstr, "Cube")
+  Cube::Cube(const std::string& confstr) : BiorthBasis(confstr, "cube")
   {
     initialize();
   }
@@ -2198,6 +2205,8 @@ namespace BasisClasses
       coef = std::make_shared<CoefClasses::CylStruct>();
     else if (name.compare("flatdisk") == 0)
       coef = std::make_shared<CoefClasses::CylStruct>();
+    else if (name.compare("cube") == 0)
+      coef = std::make_shared<CoefClasses::CubeStruct>();
     else {
       std::ostringstream sout;
       sout << "Basis::createCoefficients: basis <" << name << "> not recognized"
@@ -2279,7 +2288,8 @@ namespace BasisClasses
   }
 
   // Accumulate coefficient contributions from arrays
-  void BiorthBasis::addFromArray(Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin)
+  void BiorthBasis::addFromArray(Eigen::VectorXd& m, RowMatrixXd& p,
+				 bool RoundRobin, bool PosVelRows)
   {
     // Sanity check: is coefficient instance created?  This is not
     // foolproof.  It is really up the user to make sure that a call
@@ -2292,13 +2302,35 @@ namespace BasisClasses
       throw std::runtime_error(msg);
     }
 
-    std::vector<double> p1(3), v1(3, 0);
+    // Assume position arrays in rows by default
+    //
+    int rows = p.rows();
+    int cols = p.cols(); 
+
+    bool ambiguous = false;
 
-    if (p.rows() < 10 and p.cols() > p.rows()) {
-      std::cout << "Basis::addFromArray: interpreting your "
-		<< p.rows() << "X" << p.cols() << " input array as "
-		<< p.cols() << "X" << p.rows() << "." << std::endl;
+    if (cols==3 or cols==6) {
+      if (rows != 3 and rows != 6) PosVelRows = false;
+      else ambiguous = true;
+    }
+
+    if (rows==3 or rows==6) {
+      if (cols != 3 and cols != 6) PosVelRows = true;
+      else ambiguous = true;
+    }
+
+    if (ambiguous and myid==0) {
+      std::cout << "---- BiorthBasis::addFromArray: dimension deduction "
+		<< "is ambiguous.  Assuming that ";
+      if (PosVelRows) std::cout << "positions are in rows" << std::endl;
+      else std::cout << "positions are in columns" << std::endl;
+      std::cout << "---- BiorthBasis::addFromArray: reset 'posvelrows' flag "
+		<< "if this assumption is wrong." << std::endl;
+    }
+
+    std::vector<double> p1(3), v1(3, 0);
 
+    if (PosVelRows) {
       if (p.rows()<3) {
 	std::ostringstream msg;
 	msg << "Basis::addFromArray: you must pass a position array with at "
@@ -2308,7 +2340,7 @@ namespace BasisClasses
 
       for (int n=0; n<p.cols(); n++) {
 
-	if (n % numprocs==myid or not roundrobin) {
+	if (n % numprocs==myid or not RoundRobin) {
 
 	  bool use = true;
 	  if (ftor) {
@@ -2336,7 +2368,7 @@ namespace BasisClasses
 
       for (int n=0; n<p.rows(); n++) {
 
-	if (n % numprocs==myid or not roundrobin) {
+	if (n % numprocs==myid or not RoundRobin) {
 
 	  bool use = true;
 	  if (ftor) {
@@ -2367,10 +2399,10 @@ namespace BasisClasses
   //
   CoefClasses::CoefStrPtr BiorthBasis::createFromArray
   (Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector<double> ctr,
-   bool roundrobin)
+   bool RoundRobin, bool PosVelRows)
   {
     initFromArray(ctr);
-    addFromArray(m, p, roundrobin);
+    addFromArray(m, p, RoundRobin, PosVelRows);
     return makeFromArray(time);
   }
 

coefs/CMakeLists.txt → expui/CMakeLists.txt

@@ -29,18 +29,18 @@ if(ENABLE_XDR AND TIRPC_FOUND)
   list(APPEND common_LINKLIB ${TIRPC_LIBRARIES})
 endif()
 
-# Make the expcoef shared library
+# Make the expui shared library
 #
-set(expcoef_SOURCES BasisFactory.cc BiorthBasis.cc FieldBasis.cc
+set(expui_SOURCES BasisFactory.cc BiorthBasis.cc FieldBasis.cc
   CoefContainer.cc CoefStruct.cc FieldGenerator.cc expMSSA.cc
   Coefficients.cc KMeans.cc Centering.cc ParticleIterator.cc
   Koopman.cc)
-add_library(expcoef ${expcoef_SOURCES})
-set_target_properties(expcoef PROPERTIES OUTPUT_NAME expcoef)
-target_include_directories(expcoef PUBLIC ${common_INCLUDE})
-target_link_libraries(expcoef PUBLIC ${common_LINKLIB})
+add_library(expui ${expui_SOURCES})
+set_target_properties(expui PROPERTIES OUTPUT_NAME expui)
+target_include_directories(expui PUBLIC ${common_INCLUDE})
+target_link_libraries(expui PUBLIC ${common_LINKLIB})
 
-install(TARGETS expcoef DESTINATION lib)
+install(TARGETS expui DESTINATION lib)
 
 # Configure and build the test routines
 #
@@ -52,7 +52,7 @@ add_executable(makecoefs  makecoefs.cc)
 add_executable(testread   testread.cc)
 
 foreach(program ${bin_PROGRAMS})
-  target_link_libraries(${program} expcoef exputil ${common_LINKLIB})
+  target_link_libraries(${program} expui exputil ${common_LINKLIB})
   target_include_directories(${program} PUBLIC ${common_INCLUDE})
   target_compile_options(${program} PUBLIC ${OpenMP_CXX_FLAGS})
   install(TARGETS ${program} DESTINATION bin)

coefs/Coefficients.cc → expui/Coefficients.cc

@@ -743,6 +743,8 @@ namespace CoefClasses
   void SphCoefs::add(CoefStrPtr coef)
   {
     auto p = std::dynamic_pointer_cast<SphStruct>(coef);
+    if (not p) throw std::runtime_error("SphCoefs::add: Null coefficient structure, nothing added!");
+
     Lmax = p->lmax;
     Nmax = p->nmax;
     coefs[roundTime(coef->time)] = p;
@@ -787,6 +789,10 @@ namespace CoefClasses
       Eigen::MatrixXcd in(Mmax+1, Nmax);
       stanza.getDataSet("coefficients").read(in);
 
+      // Work around for previous unitiaized data bug; enforces real data
+      //
+      for (int n=0; n<Nmax; n++) in(0, n) = std::real(in(0, n));
+
       // Pack the data into the coefficient variable
       //
       auto coef = std::make_shared<CylStruct>();
@@ -1949,6 +1955,8 @@ namespace CoefClasses
   void CylCoefs::add(CoefStrPtr coef)
   {
     auto p = std::dynamic_pointer_cast<CylStruct>(coef);
+    if (not p) throw std::runtime_error("CylCoefs::add: Null coefficient structure, nothing added!");
+
     Mmax = p->mmax;
     Nmax = p->nmax;
     coefs[roundTime(coef->time)] = p;
@@ -1957,6 +1965,8 @@ namespace CoefClasses
   void CubeCoefs::add(CoefStrPtr coef)
   {
     auto p = std::dynamic_pointer_cast<CubeStruct>(coef);
+    if (not p) throw std::runtime_error("CubeCoefs::add: Null coefficient structure, nothing added!");
+
     NmaxX = p->nmaxx;
     NmaxY = p->nmaxy;
     NmaxZ = p->nmaxz;
@@ -1965,7 +1975,10 @@ namespace CoefClasses
 
   void TableData::add(CoefStrPtr coef)
   {
-    coefs[roundTime(coef->time)] = std::dynamic_pointer_cast<TblStruct>(coef);
+    auto p = std::dynamic_pointer_cast<TblStruct>(coef);
+    if (not p) throw std::runtime_error("TableData::add: Null coefficient structure, nothing added!");
+
+    coefs[roundTime(coef->time)] = p;
   }
 
 
@@ -1981,6 +1994,8 @@ namespace CoefClasses
   void CylFldCoefs::add(CoefStrPtr coef)
   {
     auto p = std::dynamic_pointer_cast<CylFldStruct>(coef);
+    if (not p) throw std::runtime_error("CylFldCoefs::add: Null coefficient structure, nothing added!");
+
     Nfld = p->nfld;
     Mmax = p->mmax;
     Nmax = p->nmax;

coefs/FieldBasis.H → expui/FieldBasis.H

@@ -146,7 +146,7 @@ namespace BasisClasses
     //! Initialize accumulating coefficients from arrays.  This is
     //! called once to initialize the accumulation.
     void addFromArray
-    (Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin=true);
+    (Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin=true, bool posvelrows=false);
 
     //! Accumulate new coefficients
     virtual void accumulate(double mass,