Slab force update

expui/BiorthBasis.cc

@@ -3420,7 +3420,8 @@ namespace BasisClasses
     "knots",
     "verbose",
     "check",
-    "method"
+    "method",
+    "self_consistent"
   };
 
   Slab::Slab(const YAML::Node& CONF) : BiorthBasis(CONF, "slab")

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,39 @@ 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(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.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;
   }
 };
 
@@ -2255,6 +2271,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 +2376,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,49 @@
+#pragma once
+
+#include <iostream>
+
+namespace __EXP__
+{
+  // Runtime 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
+  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 {
+	// Add the current digit to the current version component
+	current = current * 10 + (str[i] - '0');
+      }
+    }
+    *target = current;
+    return v;
+  }
+}
+

include/SLGridMP2.H

@@ -292,12 +292,21 @@ 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 versioning
+  inline static const std::string Version = "1.0";
 
   int mpi_myid, mpi_numprocs;
   int mpi_bufsz;

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,38 @@
 #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
+    condtions 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.
+
+*/
 class SlabSL : public PotAccel
 {
 
@@ -47,6 +76,9 @@ 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;
@@ -120,6 +152,12 @@ private:
 
 #endif
 
+  //! Flag self_consitency
+  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;
 

src/SlabSL.cc

@@ -17,7 +17,8 @@ SlabSL::valid_keys = {
   "hslab",
   "zmax",
   "ngrid",
-  "type"
+  "type",
+  "self_consistent"
 };
 
 //@{
@@ -48,6 +49,8 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf)
 
   int nnmax = (nmaxx > nmaxy) ? nmaxx : nmaxy;
 
+  // Make the Sturm-Liouville grid and basis functions
+  //
   grid = std::make_shared<SLGridSlab>(nnmax, nmaxz, ngrid, zmax, type);
 
   // Test for basis consistency (will generate an exception if maximum
@@ -164,6 +167,11 @@ void SlabSL::initialize()
     if (conf["hslab"])          hslab       = conf["hslab"].as<double>();
     if (conf["zmax" ])          zmax        = conf["zmax" ].as<double>();
     if (conf["type" ])          type        = conf["type" ].as<std::string>();
+
+    if (conf["self_consistent"]) {
+      self_consistent = conf["self_consistent"].as<bool>();
+    } else
+      self_consistent = true;
   }
   catch (YAML::Exception & error) {
     if (myid==0) std::cout << "Error parsing parameters in SlabSL: "
@@ -253,6 +261,11 @@ void SlabSL::determine_coefficients(void)
      compute_multistep_coefficients();
   }
 
+  // Only used if self_consistent is false to ensure that coefficients
+  // are only computed once (at the first time step)
+  //
+  firstime_coef = false;
+  if (not self_consistent) expcofF = expccof[0];
 }
 
 void * SlabSL::determine_coefficients_thread(void * arg)
@@ -343,7 +356,7 @@ void SlabSL::get_acceleration_and_potential(Component* C)
 
   if (use_external == false) {
 
-    if (multistep && initializing) {
+    if (multistep && (self_consistent || initializing)) {
       compute_multistep_coefficients();
     }
 
@@ -454,11 +467,17 @@ void * SlabSL::determine_acceleration_and_potential_thread(void * arg)
 
 	  for (int iz=0; iz<imz; iz++) {
 
-	    fac  = facx*facy*zpot[id][iz]*expccof[0](ix, iy, iz);
-	    facf = facx*facy*zfrc[id][iz]*expccof[0](ix, iy, iz);
-
-				// Limit to minimum wave number
+	    // Live for frozen potential
+	    if (self_consistent) {
+	      fac  = facx*facy*zpot[id][iz]*expccof[0](ix, iy, iz);
+	      facf = facx*facy*zfrc[id][iz]*expccof[0](ix, iy, iz);
+	    } else {
+	      fac  = facx*facy*zpot[id][iz]*expcofF(ix, iy, iz);
+	      facf = facx*facy*zfrc[id][iz]*expcofF(ix, iy, iz);
+	    }
 
+	    // Limit to minimum wave number
+	    //
 	    if (abs(ii)<nminx || abs(jj)<nminy) continue;
 
 	    potl += fac;

src/cudaSlabSL.cu

@@ -1047,8 +1047,10 @@ void SlabSL::HtoD_coefs()
   host_coefs.resize(jmax);
 
   // Copy from Slab
-  for (int i=0; i<host_coefs.size(); i++)
-    host_coefs[i] = expccof[0].data()[i];
+  for (int i=0; i<host_coefs.size(); i++) {
+    if (self_consistent) host_coefs[i] = expccof[0].data()[i];
+    else                 host_coefs[i] = expcofF.data()[i];
+  }
 
   // Copy to device
   dev_coefs = host_coefs;

utils/ICs/bonnerebert.cc

@@ -28,8 +28,6 @@
 
 */
 
-#define VERSION 0.1
-
 #include <iostream>
 #include <iomanip>
 #include <sstream>
@@ -366,8 +364,7 @@ decode_switches (int argc, char **argv)
 			   "R:"	/* runit   */
 			   "N:"	/* number  */
 			   "S:"	/* seed    */
-			   "h"	/* help    */
-			   "V",	/* version */
+			   "h", /* help    */
 			   long_options, (int *) 0)) != EOF)
     {
       switch (c)
@@ -399,10 +396,6 @@ decode_switches (int argc, char **argv)
 	case 'S':		/* --seed */
 	  S = atoi(optarg);
 	  break;
-	case 'V':
-	  cout << program_name << " " << VERSION << endl;
-	  exit (0);
-
 	case 'h':
 	  usage (0);
 

utils/Test/CMakeLists.txt

@@ -1,5 +1,5 @@
 
-set(bin_PROGRAMS testBarrier expyaml orthoTest testED)
+set(bin_PROGRAMS testBarrier expyaml orthoTest testED vtest)
 
 set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil
   yaml-cpp ${VTK_LIBRARIES})
@@ -33,6 +33,7 @@ add_executable(testBarrier    test_barrier.cc)
 add_executable(expyaml        test_config.cc)
 add_executable(orthoTest      orthoTest.cc Biorth2Ortho.cc)
 add_executable(testED         testED.cc)
+add_executable(vtest          version_test.cc)
 
 foreach(program ${bin_PROGRAMS})
   target_link_libraries(${program} ${common_LINKLIB})