Minor bug fixes and updates

CMakeLists.txt

@@ -13,6 +13,12 @@ set(CMAKE_CXX_EXTENSIONS OFF)
 # Compiler flags.  Not all tested thoroughly...
 if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
   # using Clang
+  message(STATUS "Clang compiler version is: ${CMAKE_CXX_COMPILER_VERSION}")
+  if(CMAKE_SYSTEM_NAME MATCHES "Linux")
+    if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS 18.0.0)
+      set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -stdlib=libstdc++")
+    endif()
+  endif()
 elseif(CMAKE_CXX_COMPILER_ID MATCHES "GNU")
   # using GCC
 elseif(CMAKE_CXX_COMPILER_ID MATCHES "Intel")
@@ -25,6 +31,9 @@ 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")
@@ -67,6 +76,7 @@ find_package(HDF5 COMPONENTS C CXX HL REQUIRED)
 find_package(TIRPC)	       # Check for alternative Sun rpc support
 find_package(Eigen3 REQUIRED)
 find_package(PNG)
+find_package(Git)
 
 # Check for FE
 include(FEENABLE)
@@ -183,15 +193,15 @@ add_compile_definitions(DATTRIB_CUDA=${CUDA_EXP_DATTRIB})
 
 # Get the current working branch
 execute_process(
-  COMMAND git rev-parse --abbrev-ref HEAD
+  COMMAND ${GIT_EXECUTABLE} rev-parse --abbrev-ref HEAD
   WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
   OUTPUT_VARIABLE GIT_BRANCH
   OUTPUT_STRIP_TRAILING_WHITESPACE
 )
 
 # Git submodule updates
 execute_process(
-  COMMAND git submodule update --init --recursive
+  COMMAND ${GIT_EXECUTABLE} submodule update --init --recursive
   WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
   RESULT_VARIABLE GIT_SUBMOD_RESULT
 )
@@ -204,7 +214,7 @@ endif()
 
 # Get the latest abbreviated commit hash of the working branch
 execute_process(
-  COMMAND git rev-parse HEAD
+  COMMAND ${GIT_EXECUTABLE} rev-parse HEAD
   WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
   OUTPUT_VARIABLE GIT_COMMIT
   OUTPUT_STRIP_TRAILING_WHITESPACE

INSTALL.md

@@ -1,16 +1,25 @@
 # Configuring and building EXP
 
-We are now using git submodules to provide `yaml-cpp`, which is not
-common in the HPC environments.  So, from the top-level directory, do
-the following:
+We are now using git submodules to provide a number of packages that
+are not common in the HPC environments.  These include
+
+- `HighFive`, a C++ wrapper
+class to HDF5
+- `yaml-cpp`, a C++ class for reading and emitting YAML code
+- `pybind11` which is used for Python bindings to the C++ classes
+- `rapidxml`, used to write VTK files for rendering outside of EXP
+- `png++`, C++ wrappers to the png library [Note: png support is optional]
+
+CMake will automatically download and configure these packages on the
+first call. However, if you would to do this manually, from the
+top-level directory, execute the following command:
 
 ```
    git submodule update --init --recursive
 ```
 
-This will install `yaml-cpp` in the `extern` directory.  The png++ C++
-wrappers to the png library are also installed in `extern`.  Note: png
-support is optional.
+This will install the source packages in the `extern` directory.
+
 
 ## EXP uses CMake for building a configuration
 

expui/BasisFactory.H

@@ -9,6 +9,7 @@
 #include <ParticleReader.H>
 #include <OrthoFunction.H>
 #include <Coefficients.H>
+#include <PseudoAccel.H>
 #include <YamlCheck.H>
 #include <localmpi.H>
 #include <exputils.H>
@@ -123,8 +124,22 @@ namespace BasisClasses
     //! Midplane escursion parameter
     double colh = 4.0;
 
-  public:
+    //@
+    //! Pseudo-acceleration db
+    Eigen::VectorXd t_accel;
+    Eigen::MatrixXd p_accel;
+    //@}
 
+    //! Number of center points in acceleration estimator
+    int Naccel = 0;
+
+    //! Get the current pseudo acceleration value
+    Eigen::Vector3d currentAccel(double time);
+
+  public:
+    //! The current pseudo acceleration
+    Eigen::Vector3d pseudo {0, 0, 0};
+
     //! Constructor from YAML node
     Basis(const YAML::Node& conf,     const std::string& name="Basis");
 
@@ -231,6 +246,19 @@ namespace BasisClasses
     //! Height above/below the plane for midplane search in disk scale
     //! lengths
     void setColumnHeight(double value) { colh = value; }
+
+    //@{
+    //! Initialize non-inertial forces
+    void setNonInertial(int N, Eigen::VectorXd& x, Eigen::MatrixXd& pos);
+    void setNonInertial(int N, const std::string& orient);
+    //@}
+
+    //! Set the current pseudo acceleration
+    void setNonInertialAccel(double time)
+    {
+      if (Naccel > 0) pseudo = currentAccel(time);
+    }
+
   };
 
   using BasisPtr = std::shared_ptr<Basis>;

expui/BasisFactory.cc

@@ -184,6 +184,9 @@ namespace BasisClasses
       else if ( !name.compare("flatdisk") ) {
 	basis = std::make_shared<FlatDisk>(conf);
       }
+      else if ( !name.compare("CBDisk") ) {
+	basis = std::make_shared<CBDisk>(conf);
+      }
       else if ( !name.compare("slabSL") ) {
 	basis = std::make_shared<Slab>(conf);
       }
@@ -280,5 +283,98 @@ namespace BasisClasses
     return makeFromArray(time);
   }
 
+  void Basis::setNonInertial(int N, Eigen::VectorXd& t, Eigen::MatrixXd& pos)
+  {
+    Naccel  = N;
+    t_accel = t;
+    p_accel = pos;
+  }
+
+  void Basis::setNonInertial(int N, const std::string& orient)
+  {
+    std::ifstream in(orient);
+
+    if (not in) {
+      throw std::runtime_error("Cannot open Orient file with centering data: " + orient);
+    }
+
+    const int cbufsiz = 16384;
+    std::unique_ptr<char[]> cbuf(new char [cbufsiz]);
+
+    // Look for data and write it while
+    // accumlating data for averaging
+    Eigen::Vector3d testread;
+    double time, dummy;
+
+    std::vector<double> times;
+    std::vector<Eigen::Vector3d> centers;
+
+    while (in) {
+
+      in.getline(cbuf.get(), cbufsiz);
+      if (in.rdstate() & (ios::failbit | ios::eofbit)) break;
+
+      // Skip comment lines
+      //
+      if (cbuf[0] == '#') continue;
+
+      std::istringstream line(cbuf.get());
+
+      // Read until current time is reached
+      line >> time;		// 
+      line >> dummy;
+      line >> dummy;
+
+      bool allRead = true;
+      for (int i=0; i<8; i++) {
+	if (line.eof()) allRead = false;
+	for (int k; k<3; k++) line >> testread(k);
+      }
+      if (allRead) {
+	times.push_back(time);
+	centers.push_back(testread);
+      }
+    }
+
+    // Repack data
+    Naccel = N;
+    t_accel.resize(times.size());
+    p_accel.resize(times.size(), 3);
+    for (int i=0; i<times.size(); i++) {
+      t_accel(i) = times[i];
+      for (int k=0; k<3; k++) p_accel(i, k) = centers[i](k);
+    }
+  }
+
+
+  Eigen::Vector3d Basis::currentAccel(double time)
+  {
+    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;
+    }
+    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 num = imax - imin + 1;
+      Eigen::VectorXd t(num);
+      Eigen::MatrixXd p(num, 3);
+      for (int i=imin; i<=imax; i++) {
+	t(i-imin) = t_accel(i);
+	for (int k=0; k<3; k++)	p(i-imin, k) = p_accel(i, k);
+      }
+
+      for (int k=0; k<3; k++)
+	ret(k) = 2.0*std::get<0>(QuadLS<Eigen::VectorXd>(t, p.col(k)).coefs());
+    }
+    return ret;
+  }
+
 }
 // END namespace BasisClasses