Merge pull request #113 from EXP-code/WrapperFix

`orthoCheck` Python wrapper is missing its argument. This is fixed in this PR.

Plus: bug fixes for velocity field calculations.

Author: michael-petersen

expui/BasisFactory.H

@@ -259,6 +259,10 @@ namespace BasisClasses
       if (Naccel > 0) pseudo = currentAccel(time);
     }
 
+    //! Get the field label vector
+    std::vector<std::string> getFieldLabels(void)
+    { return getFieldLabels(coordinates); }
+
   };
 
   using BasisPtr = std::shared_ptr<Basis>;

expui/Coefficients.cc

@@ -2792,7 +2792,7 @@ namespace CoefClasses
 
     for (int t=0; t<ntim; t++) {
       auto & cof = *(coefs[roundTime(times[t])]->coefs);
-      for (int i=0; i<4; i++) {
+      for (int i=0; i<Nfld; i++) {
 	for (int l=0; l<(Lmax+2)*(Lmax+1)/2; l++) {
 	  for (int n=0; n<Nmax; n++) {
 	    ret(i, l, n, t) = cof(i, l, n);

expui/FieldBasis.H

@@ -98,9 +98,13 @@ namespace BasisClasses
     virtual std::vector<double>
     crt_eval(double x, double y, double z);
 
+    //! Get the field labels
+    std::vector<std::string> getFieldLabels(const Coord ctype)
+    { return fieldLabels; }
+  
   public:
 
-    //! Constructor from YAML node
+   //! Constructor from YAML node
     FieldBasis(const YAML::Node& conf,
 	       const std::string& name="FieldBasis") : Basis(conf, name)
     { configure(); }
@@ -158,10 +162,6 @@ namespace BasisClasses
     {
     }
 
-    //! Get the field labels
-    std::vector<std::string> getFieldLabels(const Coord ctype)
-    { return fieldLabels; }
-  
     //! Return current maximum harmonic order in expansion
     int getLmax() { return lmax; }
 

expui/FieldBasis.cc

@@ -731,8 +731,8 @@ namespace BasisClasses
     double r   = sqrt(R*R + z*z);
 
     double vr  = (u*x + v*y + w*z)/r;
-    double vt  = (u*z*x + v*z*y - w*R)/R/r;
-    double vp  = (u*y - v*x)/R;
+    double vt  = (u*z*x + v*z*y - w*R*R)/R/r;
+    double vp  = (v*x - u*y)/R;
 
     return {vr, vt, vp, vr*vr, vt*vt, vp*vp};
   }

exputil/EmpCylSL.cc

@@ -16,6 +16,7 @@
 #include <thread>
 #include <exp_thread.h>
 #include <EXPException.H>
+#include <exputils.H>
 
 #include <Eigen/Eigenvalues>
 
@@ -391,6 +392,8 @@ void EmpCylSL::reset(int numr, int lmax, int mmax, int nord,
   ortho = std::make_shared<SLGridSph>(make_sl(), LMAX, NMAX, NUMR,
 				      RMIN, RMAX*0.99, false, 1, 1.0);
 
+  orthoTest(ortho->orthoCheck(std::max<int>(NMAX*50, 200)), "EmpCylSL[SLGridSph]", "l");
+
   // Resize (should not be necessary) but just in case some future
   // feature changes mulitstep on the fly
   //
@@ -868,6 +871,8 @@ int EmpCylSL::read_eof_file(const string& eof_file)
   ortho = std::make_shared<SLGridSph>(make_sl(), LMAX, NMAX, NUMR,
 					RMIN, RMAX*0.99, false, 1, 1.0);
 
+  orthoTest(ortho->orthoCheck(std::max<int>(NMAX*50, 200)), "EmpCylSL[SLGridSph]", "l");
+
   setup_eof();
   setup_accumulation();
 
@@ -1434,10 +1439,13 @@ void EmpCylSL::compute_eof_grid(int request_id, int m)
 {
   // Check for existence of ortho and create if necessary
   //
-  if (not ortho)
+  if (not ortho) {
     ortho = std::make_shared<SLGridSph>(make_sl(), LMAX, NMAX, NUMR,
 					  RMIN, RMAX*0.99, false, 1, 1.0);
 
+    orthoTest(ortho->orthoCheck(std::max<int>(NMAX*50, 200)), "EmpCylSL[SLGridSph]", "l");
+  }
+
 
   //  Read in coefficient matrix or
   //  make grid if needed
@@ -1613,11 +1621,14 @@ void EmpCylSL::compute_even_odd(int request_id, int m)
 {
   // check for ortho
   //
-  if (not ortho)
+  if (not ortho) {
     ortho = std::make_shared<SLGridSph>(make_sl(),
 					  LMAX, NMAX, NUMR, RMIN, RMAX*0.99,
 					  false, 1, 1.0);
 
+    orthoTest(ortho->orthoCheck(std::max<int>(NMAX*50, 200)), "EmpCylSL[SLGridSph]", "l");
+  }
+
   double dens, potl, potr, pott;
 
   int icnt, off;
@@ -2293,9 +2304,14 @@ void EmpCylSL::generate_eof(int numr, int nump, int numt,
 
   // Create spherical orthogonal basis if necessary
   //
-  if (not ortho)
+  if (not ortho) {
     ortho = std::make_shared<SLGridSph>(make_sl(), LMAX, NMAX, NUMR,
 					  RMIN, RMAX*0.99, false, 1, 1.0);
+
+    orthoTest(ortho->orthoCheck(std::max<int>(NMAX*50, 200)), "EmpCylSL[SLGridSph]", "l");
+  }
+
+
   // Initialize fixed variables and storage
   //
   setup_eof();
@@ -2585,9 +2601,13 @@ void EmpCylSL::generate_eof(int numr, int nump, int numt,
 void EmpCylSL::accumulate_eof(double r, double z, double phi, double mass, 
 			      int id, int mlevel)
 {
-  if (not ortho)
+  if (not ortho) {
     ortho = std::make_shared<SLGridSph>
       (make_sl(), LMAX, NMAX, NUMR, RMIN, RMAX*0.99, false, 1, 1.0);
+
+    orthoTest(ortho->orthoCheck(std::max<int>(NMAX*50, 200)), "EmpCylSL[SLGridSph]", "l");
+  }
+
   if (eof_made) {
     if (VFLAG & 2)
       cerr << "accumulate_eof: Process " << setw(4) << myid << ", Thread " 

pyEXP/BasisWrappers.cc

@@ -928,8 +928,25 @@ void BasisFactoryClasses(py::module &m)
          Returns
          -------
          None
-         )");
+         )")
+    .def("getFieldLabels",
+	 [](BasisClasses::Basis& A)
+	 {
+	   return A.getFieldLabels();
+	 },
+	 R"(
+         Provide the field labels for the basis functions
 
+         Parameters
+         ----------
+         None
+
+         Returns
+         -------
+         list: str
+           list of basis function labels
+      )"
+    );
 
   py::class_<BasisClasses::BiorthBasis, std::shared_ptr<BasisClasses::BiorthBasis>, PyBiorthBasis, BasisClasses::Basis>
     (m, "BiorthBasis")
@@ -1279,9 +1296,9 @@ void BasisFactoryClasses(py::module &m)
     // orthoCheck is not in the base class and needs to have different
     // parameters depending on the basis type.  Here, the quadrature
     // is determined by the scale of the meridional grid.
-    .def("orthoCheck", [](BasisClasses::Cylindrical& A)
+    .def("orthoCheck", [](BasisClasses::Cylindrical& A, int knots)
 	 {
-	   return A.orthoCheck();
+	   return A.orthoCheck(knots);
 	 },
 	R"(
         Check orthgonality of basis functions by quadrature
@@ -1298,7 +1315,7 @@ void BasisFactoryClasses(py::module &m)
         -------
         list(numpy.ndarray)
 	    list of numpy.ndarrays from [0, ... , Mmax]
-        )")
+        )", py::arg("knots")=400)
     .def_static("cacheInfo", [](std::string cachefile)
     {
       return BasisClasses::Cylindrical::cacheInfo(cachefile);
@@ -2044,11 +2061,11 @@ void BasisFactoryClasses(py::module &m)
       // orthoCheck is not in the base class and needs to have
       // different parameters depending on the basis type.  Here the
       // user can and will often need to specify a quadrature value.
-      .def("orthoCheck", [](BasisClasses::FieldBasis& A)
+    .def("orthoCheck", [](BasisClasses::FieldBasis& A)
       {
 	return A.orthoCheck();
       },
-	R"(
+      R"(
         Check orthgonality of basis functions by quadrature
 
         Inner-product matrix of orthogonal functions
@@ -2062,10 +2079,10 @@ void BasisFactoryClasses(py::module &m)
         numpy.ndarray
 	    orthogonality matrix
         )"
-	);
+      );
 
   py::class_<BasisClasses::VelocityBasis, std::shared_ptr<BasisClasses::VelocityBasis>, BasisClasses::FieldBasis>(m, "VelocityBasis")
-      .def(py::init<const std::string&>(),
+    .def(py::init<const std::string&>(),
 	 R"(
          Create a orthogonal velocity-field basis
 

pyEXP/TensorToArray.H

@@ -11,7 +11,7 @@ py::array_t<T> make_ndarray3(Eigen::Tensor<T, 3>& mat)
   // Check rank
   if (dims.size() != 3) {
     std::ostringstream sout;
-    sout << "make_ndarray: tensor rank must be 3, found " << dims.size();
+    sout << "make_ndarray3: tensor rank must be 3, found " << dims.size();
     throw std::runtime_error(sout.str());
   }
 
@@ -37,10 +37,17 @@ py::array_t<T> make_ndarray4(Eigen::Tensor<T, 4>& mat)
   // Check rank
   if (dims.size() != 4) {
     std::ostringstream sout;
-    sout << "make_ndarray: tensor rank must be 4, found " << dims.size();
+    sout << "make_ndarray4: tensor rank must be 4, found " << dims.size();
     throw std::runtime_error(sout.str());
   }
 
+  // Sanity check
+  for (int i=0; i<mat.size(); i++) {
+    if (isnan(std::abs(mat.data()[i]))) {
+      throw std::runtime_error("make_ndarray4: NaN encountered");
+    }
+  }
+
   // Make the memory mapping
   return py::array_t<T>
     (
@@ -107,11 +114,11 @@ Eigen::Tensor<T, 4> make_tensor4(py::array_t<T> array)
 
   // Build result tensor with col-major ordering
   Eigen::Tensor<T, 4> tensor(shape[0], shape[1], shape[2], shape[3]);
-  for (int i=0, l=0; i < shape[0]; i++) {
+  for (int i=0, c=0; i < shape[0]; i++) {
     for (int j=0; j < shape[1]; j++) {
       for (int k=0; k < shape[2]; k++) {
-	for (int l=0; l < shape[3]; k++) {
-	  tensor(i, j, k, l) = data[l++];
+	for (int l=0; l < shape[3]; l++, c++) {
+	  tensor(i, j, k, l) = data[c];
 	}
       }
     }