Slab force update

expui/BiorthBasis.cc

@@ -3420,7 +3420,9 @@ namespace BasisClasses
     "knots",
     "verbose",
     "check",
-    "method"
+    "method",
+    "self_consistent",
+    "cachename"
   };
 
   Slab::Slab(const YAML::Node& CONF) : BiorthBasis(CONF, "slab")
@@ -3481,6 +3483,8 @@ namespace BasisClasses
       if (conf["knots"])      knots = conf["knots"].as<int>();
 
       if (conf["check"])      check = conf["check"].as<bool>();
+
+      if (conf["cachename"])  cachename = conf["cachename"].as<std::string>();
     } 
     catch (YAML::Exception & error) {
       if (myid==0) std::cout << "Error parsing parameter stanza for <"
@@ -3493,6 +3497,19 @@ namespace BasisClasses
       throw std::runtime_error("Slab: error parsing YAML");
     }
 
+    // Check for non-null cache file name.  This must be specified
+    // to prevent recomputation and unexpected behavior.
+    //
+    if (cachename.size() == 0) {
+      throw std::runtime_error
+	("SlabSL requires a specified cachename in your YAML config\n"
+	 "for consistency with previous invocations and existing coefficient\n"
+	 "sets.  Please add explicitly add 'cachename: name' to your config\n"
+	 "with new 'name' for creating a basis or an existing 'name' for\n"
+	 "reading a previously generated basis cache\n");
+    }
+
+
     // Finally, make the basis
     //
     SLGridSlab::mpi  = 0;
@@ -3502,7 +3519,7 @@ namespace BasisClasses
 
     int nnmax = (nmaxx > nmaxy) ? nmaxx : nmaxy;
 
-    ortho = std::make_shared<SLGridSlab>(nnmax, nmaxz, ngrid, zmax, type);
+    ortho = std::make_shared<SLGridSlab>(nnmax, nmaxz, ngrid, zmax, cachename, type);
 
     // Orthogonality sanity check
     //

exputil/SLGridMP2.cc

@@ -21,6 +21,7 @@
 #include "SLGridMP2.H"
 #include "massmodel.H"
 #include "EXPmath.H"
+#include "libvars.H"		// For parsed version info
 
 #ifdef USE_DMALLOC
 #include <dmalloc.h>
@@ -1844,24 +1845,40 @@ static double poffset=0.0;
 class IsothermalSlab : public SlabModel
 {
 
+private:
+
+  std::string psa =
+    "----     IsothermalSlab NOW uses the traditional density profile proportional to sech^2(z/2H).\n"
+    "----     If you are using the old profile proportional to sech^2(z/H), please update your model\n"
+    "----     to use the new profile and set the scale height H to be half of the old value.  This\n"
+    "----     will ensure that your model has the same density profile and potential as before, but\n"
+    "----     with a more standard functional form.  If you have any questions or concerns about\n"
+    "----     this change, please contact the developers on GitHub.";
+
 public:
 
-  IsothermalSlab() { id = "iso"; }
+  IsothermalSlab() {
+    id = "iso";
+    if (myid==0 and (exp_build.major < 7 or
+                     (exp_build.major == 7 and exp_build.minor < 11)))
+      std::cout << "---- SLGridSlab: IMPORTANT UPDATE for EXP "
+		<< VERSION << '\n' << psa << std::endl;
+  }
 
   double pot(double z)
   {
-    return 2.0*M_PI*SLGridSlab::H*log(cosh(z/SLGridSlab::H)) - poffset;
+    return 4.0*M_PI*SLGridSlab::H*log(cosh(0.5*z/SLGridSlab::H)) - poffset;
   }
 
   double dpot(double z)
   {
-    return 2.0*M_PI*tanh(z/SLGridSlab::H);
+    return 2.0*M_PI*tanh(0.5*z/SLGridSlab::H);
   }
 
   double dens(double z)
   {
-    double tmp = 1.0/cosh(z/SLGridSlab::H);
-    return 4.0*M_PI * 0.5/SLGridSlab::H * tmp*tmp;
+    double tmp = 1.0/cosh(0.5*z/SLGridSlab::H);
+    return 4.0*M_PI*0.25/SLGridSlab::H * tmp*tmp;
   }
 };
 
@@ -1956,8 +1973,12 @@ void SLGridSlab::bomb(string oops)
 				// Constructors
 
 SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX,
-		       const std::string TYPE, bool VERBOSE)
+		       const std::string cachename, const std::string TYPE,
+		       bool VERBOSE)
 {
+  if (cachename.size()) slab_cache_name  = cachename;
+  else throw std::runtime_error("SLGridSlab: you must specify a cachename");
+
   int kx, ky;
 
   numk = NUMK;
@@ -2193,9 +2214,6 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX,
 }
 
 
-const string slab_cache_name = ".slgrid_slab_cache";
-
-
 bool SLGridSlab::ReadH5Cache(void)
 {
   if (!cache) return false;
@@ -2255,6 +2273,18 @@ bool SLGridSlab::ReadH5Cache(void)
     if (not checkStr(geometry, "geometry"))  return false;
     if (not checkStr(forceID,  "forceID"))   return false;
 
+    // Version check
+    //
+    if (h5file.hasAttribute("Version")) {
+      if (not checkStr(Version, "Version"))  return false;
+    } else {
+      if (myid==0)
+	std::cout << "---- SLGridSlab::ReadH5Cache: "
+		  << "recomputing cache for HighFive API change"
+		  << std::endl;
+      return false;
+    }
+
     // Parameter check
     //
     if (not checkStr(type,     "type"))      return false;
@@ -2348,6 +2378,7 @@ void SLGridSlab::WriteH5Cache(void)
 
     file.createAttribute<std::string>("geometry",  HighFive::DataSpace::From(geometry)).write(geometry);
     file.createAttribute<std::string>("forceID",   HighFive::DataSpace::From(forceID)).write(forceID);
+    file.createAttribute<std::string>("Version",   HighFive::DataSpace::From(Version)).write(Version);
 
     // Write parameters
     //

exputil/libvars.cc

@@ -3,6 +3,7 @@
    standalone utilities
  */
 
+#include "config_exp.h"
 #include "libvars.H"
 
 namespace __EXP__
@@ -36,6 +37,7 @@ namespace __EXP__
 
   //! Sanity tolerance for orthogonality
   double           orthoTol         = 1.0e-2;
+
 };
 
 

include/EXPversion.H

@@ -0,0 +1,48 @@
+#pragma once
+
+namespace __EXP__
+{
+  // Compile-time parser for "X.Y.Z" format
+
+  /* Example usage:
+
+     #define VERSION_STR "2.4.12"
+     static constexpr EXPversion current_v = EXP_parse_version(VERSION_STR);
+
+     int main() {
+       static_assert(current_v.major == 2, "Major version mismatch");
+       std::cout << "Parsed: " << current_v.major << "." << current_v.minor << "\n";
+     return 0;}
+  */
+
+
+  struct EXPversion
+  {
+    int major, minor, patch;
+  };
+
+  // Compile-time parser for "X.Y.Z" format using c++-17 features
+  constexpr EXPversion EXP_parse_version(const char* str)
+  {
+    EXPversion v = {0, 0, 0};
+    int* target = &v.major;
+    int current = 0;
+
+    // Scan the version string until we hit a dot, then move to the next target
+    for (int i = 0; str[i] != '\0'; ++i) {
+      if (str[i] == '.') {
+	*target = current;
+	if      (target == &v.major) target = &v.minor;
+	else if (target == &v.minor) target = &v.patch;
+	current = 0;
+      } else {
+	// Append the current decimal digit to the current version
+	current = current * 10 + (str[i] - '0');
+      }
+    }
+    // Store the int of this version string
+    *target = current;		
+    return v;
+  }
+}
+

include/SLGridMP2.H

@@ -292,12 +292,24 @@ private:
   void compute_table_worker(void);
 
 
-				// Local MPI stuff
+  //@{
+  //! Local MPI stuff
   void mpi_setup(void);
   void mpi_unpack_table(void);
   int  mpi_pack_table(TableSlab* table, int kx, int ky);
+  //@}
+
+  //@{
+  //! Cache reading and writing
   bool ReadH5Cache(void);
   void WriteH5Cache(void);
+  //@}
+
+  //! Cache file name
+  std::string slab_cache_name;
+
+  //! Cache versioning
+  inline static const std::string Version = "1.0";
 
   int mpi_myid, mpi_numprocs;
   int mpi_bufsz;
@@ -403,7 +415,8 @@ public:
 
   //! Constructor
   SLGridSlab(int kmax, int nmax, int numz, double zmax,
-	     const std::string type="isothermal", bool Verbose=false);
+	     const std::string cachename, const std::string type="isothermal",
+	     bool Verbose=false);
 
   //! Destructor
   ~SLGridSlab();

include/libvars.H

@@ -7,6 +7,9 @@
 
 #include <mpi.h>
 
+#include "config_exp.h"		// For version string
+#include "EXPversion.H"		// For compile-time version parsing
+
 namespace __EXP__
 {
   //! @name Theading variables
@@ -42,6 +45,9 @@ namespace __EXP__
   //! Sanity tolerance for orthogonality
   extern double orthoTol;
 
+  //! Compile-time version parsing
+  static constexpr EXPversion exp_build = EXP_parse_version(VERSION);
+
 };
 
 #endif	// END _LIBVARS_H

src/SlabSL.H

@@ -17,9 +17,40 @@
 #include <cudaMappingConstants.cuH>
 #endif
 
-/*! This routine computes the potential, acceleration and density
-  using expansion periodic in X & Y and outgoing vacuum boundary
-  condtions in Z */
+/** @class SlabSL
+    @brief This routine computes the potential, acceleration and density
+    using expansion periodic in X & Y and outgoing vacuum boundary
+    conditions in Z
+
+    @details **YAML configuration**
+
+    @param nmaxx is the maximum order of the expansion in x (default 6)
+
+    @param nmaxy is the maximum order of the expansion in y (default 6)
+
+    @param nmaxz is the maximum order of the expansion in z (default 6)
+
+    @param nminx is the minimum order of the expansion in x (default 0)
+
+    @param nminy is the minimum order of the expansion in y (default 0)
+
+    @param hslab is the scale height of the slab (default 0.2)
+
+    @param zmax is the maximum z for the slab (default 10.0)
+
+    @param ngrid is the number of grid points in z for the
+    Sturm-Liouville solver (default 1000)
+
+    @param type is the type of slab to solve for (default
+    "isothermal", must be "isothermal", "parabolic", or "constant")
+
+    @param self_consistent set to true allows the particles to evolve
+    under the time-dependent basis expansion. For a basis fixed in time to
+    the initial time: set to false.
+
+    @param cachename is the name of the basis cache file.  This is a
+    required parameter.
+*/
 class SlabSL : public PotAccel
 {
 
@@ -47,13 +78,17 @@ private:
   //! Current coefficient tensor
   std::vector<coefType> expccof, expccofP;
 
+  //! Coefficient tensor for frozen potential (if self_consistent=false)
+  coefType expcofF;
+
   int nminx, nminy;
   int nmaxx, nmaxy, nmaxz;
   double zmax, hslab;
 
   int imx, imy, imz, jmax, nnmax;
   double  dfac;
   std::complex<double> kfac;
+  std::string cachename;
 
   std::vector<Eigen::VectorXd> zfrc, zpot;
 
@@ -120,6 +155,12 @@ private:
 
 #endif
 
+  //! Flag self_consistency
+  bool self_consistent = true;
+
+  //! Flag whether coefficients have been initialized for the first time
+  bool firstime_coef = true;
+
   //! Default number of grid points for SLGridSlab
   int ngrid = 1000;