Merge branch 'devel' into cylheight for upstream consistency.

Author: Martin D. Weinberg

.gitattributes

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features
 
 project(
   EXP
-  VERSION "7.8.5"
+  VERSION "7.9.1"
   HOMEPAGE_URL https://github.com/EXP-code/EXP
   LANGUAGES C CXX Fortran)
 

CMakeLists.txt

@@ -1,3 +1,7 @@
+See the installation guide in the documentation for the most up-to-date instructions on
+how to install EXP:
+https://exp-docs.readthedocs.io/en/latest/intro/install.html
+
 # Configuring and building EXP
 
 EXP uses CMake for building a configuration.

INSTALL

@@ -4,24 +4,28 @@
 [![Documentation](https://readthedocs.org/projects/exp-docs/badge/?version=latest)](https://exp-docs.readthedocs.io/)
 [![DOI](https://joss.theoj.org/papers/10.21105/joss.07302/status.svg)](https://doi.org/10.21105/joss.07302)
 
+EXP is a Basis Function Expansion C++ library and Python package for running N-body galactic simulations and dynamical discovery.
 
 ## Repo organization
 
-| Files   | Description |
+| Directory   | Description |
 | ---     | ---         |
-| README  | This file |
-| TODO	  | Wish list of features to add, things to fix, big and small |
-| include | Include file for all common classes |
-| src	  | Source for nbody code |
+| Paper   | The JOSS paper for EXP |
+| cmake   | Helpers for `cmake` installation process |
+| doc     | The doxygen documentation for EXP |
+| expui   | Classes used by Python wrappers and by the nbody code and other standalone utilities for coefficient IO support. |
 | exputil | The main EXP function and class library used by both pyEXP and EXP |
-| coefs   | Source for standalone classes used by Python wrappers and other standalone utilities |
+| extern  | The directory for external code needed for installation |
+| include | Include file for all common classes |
 | pyEXP   | Source for Python wrappers |
+| src	  | Source for nbody code |
+| tests   | An assortment of automatic tests for EXP |
 | utils   | Older but still useful standalone C++ utilities |
 
 ## Version reporting
 
 EXP automatically stashes its compile time, git branch, and git commit
-hash when `make` is invoked in the src directory.  You can see this
+hash when `make` is invoked in the src directory (e.g. as part of the `cmake` installation).  You can see this info using `pyEXP` using the `pyEXP.util.getVersionInfo()` command. You can also see this
 info using the -v flag, i.e. `mpirun -np 1 exp -v` or `exp -v`.  Note:
 some MPI implementations require the MPI-aware executable to be run
 using 'mpirun'.  Some recent HPC systems using `slurm` require the use
@@ -30,7 +34,7 @@ resource request.
 
 ## Compile hints
 
-See INSTALL.md for a brief synposis.
+See INSTALL.md for a brief synposis and pointers to more detailed information.
 
 A few quick additional notes. By default, both the n-body code and the
 Python bindings, the pyEXP interface, will be compiled by default.
@@ -54,7 +58,7 @@ We are developing two repositories of examples and tutorials:
 | EXP-examples | Each subdirectory contains a full set of body files and configurations to run EXP with with model galaxy |
 | pyEXP-examples | Tutorials and example workflows for a variety of envisioned use cases |
 
-Both of these are available from the origin as EXP.
+Both of these are available from the [origin as EXP](https://github.com/EXP-code).
 
 ## pyEXP
 
@@ -94,11 +98,9 @@ To provide some context, suppose you want to read some snapshots, make some coef
 5. 'newcoefs' and 'disk' can then be passed to the FieldGenerator to provide density, potential, force fields, etc. for the each principal signal
 
 This is only one example of many possible uses.  There are many
-variants to this work flow, of course, and I expect that you will
+variants to this work flow, of course, and we expect that you will
 invent some interesting ones.
 
 ## Citation
 
-If you use EXP in your research, please cite the following paper:
-
-[![DOI](https://joss.theoj.org/papers/10.21105/joss.07302/status.svg)](https://doi.org/10.21105/joss.07302)
+If you use EXP in your research, please cite the [JOSS paper](https://doi.org/10.21105/joss.07302), and consider including other papers from the list in CITATIONS.bib.

INSTALL.md

@@ -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

README.md

@@ -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); }

doc/exp.cfg

@@ -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/BasisFactory.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;
   };
 
   /**
@@ -307,6 +311,9 @@ namespace BasisClasses
       return ret;
     }
 
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    virtual std::vector<double> getAccel(double x, double y, double z);
   };
 
   /**
@@ -542,6 +549,10 @@ namespace BasisClasses
       return ret;
     }
 
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    std::vector<double> getAccel(double x, double y, double z);
+
   };
 
   /**
@@ -693,6 +704,10 @@ namespace BasisClasses
       return ret;
     }
 
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    std::vector<double> getAccel(double x, double y, double z);
+
   };
 
   /**
@@ -841,10 +856,18 @@ namespace BasisClasses
     bool orthoTest()
     {
       auto [ret, worst, lworst] = orthoCompute(sl->orthoCheck());
-      // For the CTest log
-      std::cout << "---- Cylindrical::orthoTest: worst=" << worst << std::endl;
+      if (myid==0) {
+	// For the CTest log
+	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);
+
   };
 
   /**
@@ -991,6 +1014,10 @@ namespace BasisClasses
       return true;
     }
 
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    std::vector<double> getAccel(double x, double y, double z);
+
   };
 
   /**
@@ -1113,6 +1140,10 @@ namespace BasisClasses
       return true;
     }
 
+    //! Evaluate acceleration in Cartesian coordinates in centered
+    //! coordinate system
+    std::vector<double> getAccel(double x, double y, double z);
+
   };