Merge pull request #145 from EXP-code/cylheight

Improve scale height consistency

Author: michael-petersen

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features
 
 project(
   EXP
-  VERSION "7.8.6"
+  VERSION "7.9.1"
   HOMEPAGE_URL https://github.com/EXP-code/EXP
   LANGUAGES C CXX Fortran)
 
@@ -268,7 +268,13 @@ endif()
 # try to find pybind11 and build wrapper python module
 find_package(Python3 COMPONENTS Interpreter Development)
 message(STATUS "python3 include dirs: ${Python3_INCLUDE_DIRS}")
-
+if(Python3_FOUND)
+  set(HAVE_PYTHON3 TRUE)
+else()
+  if(ENABLE_PYEXP)
+    message(FATAL_ERROR "You asked for pyEXP but I cannot find a Python3 environment. Please make Python3 available or disable pyEXP. CMake will exit." )
+  endif()
+endif()
 
 # Force installation of the yaml-cpp libraries
 install(TARGETS yaml-cpp DESTINATION lib)

config_cmake.h_in

@@ -19,6 +19,9 @@
 /* Defined if you have HDF5 support */
 #cmakedefine HAVE_HDF5 @HAVE_HDF5@
 
+/* Defined if Python3 runtime exists */
+#cmakedefine HAVE_PYTHON3 1
+
 /* Define to 1 if you have the `cuda' library. */
 #cmakedefine HAVE_LIBCUDA 1
 

doc/exp.cfg

@@ -48,7 +48,7 @@ PROJECT_NAME           = EXP
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = 7.8.0
+PROJECT_NUMBER         = 7.9.1
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a

expui/BasisFactory.H

@@ -137,6 +137,7 @@ namespace BasisClasses
     Eigen::Vector3d currentAccel(double time);
 
   public:
+
     //! The current pseudo acceleration
     Eigen::Vector3d pseudo {0, 0, 0};
 
@@ -249,8 +250,8 @@ namespace BasisClasses
 
     //@{
     //! Initialize non-inertial forces
-    void setNonInertial(int N, Eigen::VectorXd& x, Eigen::MatrixXd& pos);
-    void setNonInertial(int N, const std::string& orient);
+    void setNonInertial(int Naccel, const Eigen::VectorXd& x, const Eigen::MatrixXd& pos);
+    void setNonInertial(int Naccel, const std::string& orient);
     //@}
 
     //! Set the current pseudo acceleration
@@ -259,6 +260,16 @@ namespace BasisClasses
       if (Naccel > 0) pseudo = currentAccel(time);
     }
 
+    //! Returns true if non-inertial forces are active
+    bool usingNonInertial() { return Naccel > 0; }
+
+    //! Reset to inertial coordinates
+    void setInertial()
+    {
+      Naccel = 0;
+      pseudo = {0, 0, 0};
+    }
+
     //! Get the field label vector
     std::vector<std::string> getFieldLabels(void)
     { return getFieldLabels(coordinates); }

expui/BasisFactory.cc

@@ -283,8 +283,16 @@ namespace BasisClasses
     return makeFromArray(time);
   }
 
-  void Basis::setNonInertial(int N, Eigen::VectorXd& t, Eigen::MatrixXd& pos)
+  void Basis::setNonInertial(int N, const Eigen::VectorXd& t, const Eigen::MatrixXd& pos)
   {
+    // Sanity checks
+    if (t.size() < 1)
+      throw std::runtime_error("Basis: setNonInertial: no times in time array");
+
+    if (t.size() != pos.rows())
+      throw std::runtime_error("Basis::setNonInertial: size mismatch in time and position arrays");
+
+    // Set the data
     Naccel  = N;
     t_accel = t;
     p_accel = pos;
@@ -351,16 +359,31 @@ namespace BasisClasses
   {
     Eigen::Vector3d ret;
 
-    if (time < t_accel(0) || time > t_accel(t_accel.size()-1)) {
-      std::cout << "ERROR: " << time << " is outside of range of the non-inertial DB"
-		<< std::endl;
-      ret.setZero();
-      return ret;
+    auto n = t_accel.size();
+
+    // Allow a little bit of buffer in the allowable on-grid range but
+    // otherwise force termination
+    //
+    if ( time < t_accel(0  ) - 0.5*(t_accel(1  ) - t_accel(0  ))  ||
+	 time > t_accel(n-1) + 0.5*(t_accel(n-1) - t_accel(n-2)) ) {
+      
+      std::ostringstream msg;
+      msg << "Basis::currentAccel: " << time
+	  << " is outside the range of the non-inertial DB ["
+	  << t_accel(0) << ", " << t_accel(n-1) << "]";
+
+      throw std::runtime_error(msg.str());
     }
+    // Do the quadratic interpolation
+    //
     else {
       int imin = 0;
-      int imax = std::lower_bound(t_accel.data(), t_accel.data()+t_accel.size(), time) - t_accel.data();
-      if (imax > Naccel) imin = imax - Naccel;
+      int imax = std::lower_bound(t_accel.data(), t_accel.data()+n, time) - t_accel.data();
+
+      // Get a range around the current time of approx size Naccel
+      //
+      imax = std::min<int>(n-1, imax + Naccel/2);
+      imin = std::max<int>(imax - Naccel, 0);
 
       int num = imax - imin + 1;
       Eigen::VectorXd t(num);

expui/BiorthBasis.H

@@ -163,6 +163,10 @@ namespace BasisClasses
 
     //! Clear the particle selector callback
     void clrSelector() { ftor = nullptr; }
+
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    virtual std::vector<double> getAccel(double x, double y, double z) = 0;
   };
 
   /**
@@ -303,10 +307,14 @@ namespace BasisClasses
     {
       auto [ret, worst, lworst] = orthoCompute(orthoCheck(knots));
       // For the CTest log
-      std::cout << "---- Spherical::orthoTest: worst=" << worst << std::endl;
+      if (myid==0)
+	std::cout << "---- Spherical::orthoTest: worst=" << worst << std::endl;
       return ret;
     }
 
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    virtual std::vector<double> getAccel(double x, double y, double z);
   };
 
   /**
@@ -538,10 +546,15 @@ namespace BasisClasses
     {
       auto [ret, worst, lworst] = orthoCompute(orthoCheck());
       // For the CTest log
-      std::cout << "---- FlatDisk::orthoTest: worst=" << worst << std::endl;
+      if (myid==0)
+	std::cout << "---- FlatDisk::orthoTest: worst=" << worst << std::endl;
       return ret;
     }
 
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    std::vector<double> getAccel(double x, double y, double z);
+
   };
 
   /**
@@ -689,10 +702,15 @@ namespace BasisClasses
     {
       auto [ret, worst, lworst] = orthoCompute(orthoCheck());
       // For the CTest log
-      std::cout << "CBDisk::orthoTest: worst=" << worst << std::endl;
+      if (myid==0)
+	std::cout << "CBDisk::orthoTest: worst=" << worst << std::endl;
       return ret;
     }
 
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    std::vector<double> getAccel(double x, double y, double z);
+
   };
 
   /**
@@ -718,7 +736,7 @@ namespace BasisClasses
     int rnum, pnum, tnum;
     double rmin, rmax, rcylmin, rcylmax;
     double acyl, hcyl;
-    bool expcond, logarithmic, density, EVEN_M;
+    bool expcond, logarithmic, density, EVEN_M, sech2 = false;
 
     std::vector<Eigen::MatrixXd> potd, dpot, dpt2, dend;
     std::vector<Eigen::MatrixXd> legs, dlegs, d2legs;
@@ -842,9 +860,15 @@ namespace BasisClasses
     {
       auto [ret, worst, lworst] = orthoCompute(sl->orthoCheck());
       // For the CTest log
-      std::cout << "---- Cylindrical::orthoTest: worst=" << worst << std::endl;
+      if (myid==0)
+	std::cout << "---- Cylindrical::orthoTest: worst=" << worst << std::endl;
       return ret;
     }
+
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    std::vector<double> getAccel(double x, double y, double z);
+
   };
 
   /**
@@ -985,12 +1009,18 @@ namespace BasisClasses
 	  }
 	}
       }
-	
+
+      if (myid==0)
       std::cout << "---- Slab::orthoTest: worst=" << worst << std::endl;
+
       if (worst > __EXP__::orthoTol) return false;
       return true;
     }
 
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    std::vector<double> getAccel(double x, double y, double z);
+
   };
 
   /**
@@ -1108,11 +1138,17 @@ namespace BasisClasses
 	}
       }
 
-      std::cout << "---- Cube::orthoTest: worst=" << worst << std::endl;
+      if (myid==0)
+	std::cout << "---- Cube::orthoTest: worst=" << worst << std::endl;
+
       if (worst > __EXP__::orthoTol) return false;
       return true;
     }
 
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    std::vector<double> getAccel(double x, double y, double z);
+
   };