Slab force update

[The9Cat]

Updates for consistency with the current conventions

1. Density for isothermal slab is now proportional to sech(z/2h)^2 rather than sech(z/h)^2
2. SlablSL now implements the self_consistent flag for enabling a frozen potential
3. Added a version # to the slab cache along with a user warning about this change

One general development support feature

Added compile -time parsing of the EXP version string into major, minor, patch for internal version control logic. E.g.

#include <iostream>
#include "libvars.H"
using namespace __EXP__;

int main()
{
  std::cout << "Major=" << exp_build.major << '\n';
  std::cout << "Minor=" << exp_build.minor << '\n';
  std::cout << "Patch=" << exp_build.patch << '\n';
  return 0;
}

This will make future version-dependent toggles easier to implement. We have a usage example in both SLGridSlab and in utils/Test.

[Copilot]

Pull request overview

This PR updates the slab force/basis implementation to match current conventions (isothermal slab profile normalization, cache versioning, and a self_consistent “frozen potential” mode) and adds compile-time parsing of the EXP version string for internal version-dependent logic.

Changes:

• Switch isothermal slab model normalization to use a sech^2(z/2H) form and warn users when running older EXP versions.
• Add a self_consistent configuration flag to SlabSL to support freezing the potential coefficients.
• Add compile-time EXP version parsing (EXPversion.H, exp_build) and a small test utility; add slab cache version attribute.

Reviewed changes

Copilot reviewed 12 out of 12 changed files in this pull request and generated 7 comments.

Show a summary per file

File	Description
utils/Test/version_test.cc	Adds a small example program that prints VERSION and parsed major/minor/patch.
utils/Test/CMakeLists.txt	Builds the new vtest utility.
utils/ICs/bonnerebert.cc	Removes the local VERSION macro and the -V version flag handling.
src/SlabSL.H	Documents self_consistent and adds storage for frozen coefficients.
src/SlabSL.cc	Adds YAML parsing for self_consistent and uses frozen coefficients during force evaluation.
src/cudaSlabSL.cu	Uses frozen coefficients when copying coefficients to the GPU if self_consistent=false.
include/SLGridMP2.H	Adds cache versioning for slab cache and improves internal grouping/comments.
include/libvars.H	Exposes exp_build as a compile-time parsed {major,minor,patch} version triple.
include/EXPversion.H	Introduces a constexpr version-string parser.
exputil/SLGridMP2.cc	Updates isothermal slab profile functions and adds slab cache version checks/write.
exputil/libvars.cc	Ensures config header is included in libvars translation unit.
expui/BiorthBasis.cc	Allows self_consistent as a recognized YAML key for the slab basis config.

Comments suppressed due to low confidence (1)

utils/ICs/bonnerebert.cc:404

• The -V/--version option was removed from getopt_long parsing, but long_options still defines a version entry and usage() still advertises -V. This makes the CLI help incorrect and --version will likely behave unexpectedly (depending on getopt_long behavior). Either restore handling for -V/--version or remove it consistently from long_options and the usage text.

			   "R:"	/* runit   */
			   "N:"	/* number  */
			   "S:"	/* seed    */
			   "h", /* help    */
			   long_options, (int *) 0)) != EOF)
    {
      switch (c)
	{
	case 'o':		/* --output */
	  oname = string(optarg);
	  break;
	case 'd':		/* --xstep */
	  h = atof(optarg);
	  break;
	case 'r':		/* --ratio */
	  Rratio = atof(optarg);
	  break;
	case 'R':		/* --runit */
	  R = atof(optarg);
	  break;
	case 'M':		/* --munit */
	  Munit = atof(optarg);
	  break;
	case 'm':		/* --mass */
	  M = atof(optarg);
	  break;
	case 'T':		/* --temp */
	  T = atof(optarg);
	  break;
	case 'N':		/* --number */
	  N = atoi(optarg);
	  break;
	case 'S':		/* --seed */
	  S = atoi(optarg);
	  break;
	case 'h':
	  usage (0);

	default:
	  usage (-1);
	}

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[The9Cat]

@michael-petersen I am probably the only one that uses the slab stuff and these changes implement some of the standard features from the other force classes. So no need for you to spend time on most of this. It's all tested and seems to work.

However, note that this PR adds a compile-time function that parses the CMake-generated VERSION string into a structure of three ints (major, minor, patch) that may be used to gate various features rather than reverting to some sort of on-the-fly parse of the VERSION for each case or using CPP macros (yuck). The structure, named EXPversion, is a global library variable compiled into exputil. E.g. this PR uses the EXPversion instance to alert the user to the scale height convention change for the the isothermal disk. Your comments on this bit will be helpful I think.

[Copilot]

@michael-petersen I am probably the only one that uses the slab stuff and these changes implement some of the standard features from the other force classes. So no need for you to spend time on most ...

The version_test.cc file now consistently uses "integer triplet" in both the file comment and the runtime output (commit b7ed7eb). The other resolved review comments were already addressed in earlier commits on this branch.

[The9Cat]

@copilot The slab force is implemented in SlabSL.cc and the Cuda implementation is in cudaSlabSL.cu. There are two parts to this computation: the coefficient accumulation from particles and the force/acceleration computation from the coefficients. Are you able to study these implementations and flag any possible differences that would lead to significantly results?

[Copilot]

@copilot The slab force is implemented in SlabSL.cc and the Cuda implementation is in cudaSlabSL.cu. There are two parts to this computation: the coefficient accumulation from particles and the f...

Here are the three differences I found that could lead to significantly different results:

---

1. Texture lookup: missing (kx, ky) swap when abs(jj) > abs(ii) — both coefficient and force kernels

The textures are initialized for the upper-triangular set (kx, ky) with kx ≥ ky (cudaSLGridMP2.cu:318-350), so the mapping formula 1 + (kx*(kx+1)/2 + ky)*slabNumZ is only correct when kx = abs(ii) ≥ ky = abs(jj). The CPU path explicitly swaps:

if (iix >= iiy) grid->get_pot(zpot, zz, iix, iiy);
else            grid->get_pot(zpot, zz, iiy, iix);

The CUDA kernels omit this swap. For every (ii, jj) pair where abs(jj) > abs(ii), both coefKernelSlab (line 386) and forceKernelSlab (line 504) fetch the wrong texture — one that was stored for a different transverse wavenumber k². For example, (ii=0, jj=2) should use the k²=4 texture but instead fetches the k²=2 texture for (kx=1, ky=1). The fix is:

int kx = max(abs(ii), abs(jj));
int ky = min(abs(ii), abs(jj));
int offset = 1 + (kx*(kx+1)/2 + ky)*slabNumZ;

---

2. Opposite sign convention for the potential term at z < 0 — forceKernelSlab only

coefKernelSlab (line 394) correctly flips sign for odd modes (antisymmetric potential, φ(−z) = −φ(z)):

if (pos[2]<0 && 2*(n/2)!=n) sign = -1;   // odd n

forceKernelSlab (line 511) flips sign for even modes:

if (pos[2]<0 && 2*(n/2)==n) sign = -1;   // even n

That sign is correct for the z-force f (whose even-mode part is antisymmetric), but it is applied to the potential v as well via the same variable. For v, even modes are symmetric (no flip needed), and odd modes are antisymmetric (flip needed). So the potential term in forceKernelSlab has the sign backwards for all modes at z < 0. Since acc[0] and acc[1] are derived from fac = X*Y*v*coef, the transverse accelerations are also wrong for all particles at z < 0 (~half the particles in a slab simulation). The fix requires two separate sign variables:

int signV = 1, signF = 1;
if (pos[2]<0 && 2*(n/2)!=n) signV = -1;  // potential: flip odd modes
if (pos[2]<0 && 2*(n/2)==n) signF = -1;  // z-force:   flip even modes

and apply signV to v, signF to f.

---

3. Missing Adiabatic() factor in coefKernelSlab

CPU (determine_coefficients_thread, line 486):

double adb = cC->Adiabatic();
double mm  = -4.0*M_PI * cC->Mass(i) * adb;

CUDA (coefKernelSlab, line 323):

cuFP_t mm = - p.mass * 2.0*slabDfac;   // = -4π * mass, no Adiabatic()

For standard runs Adiabatic() == 1 so this has no effect, but for any adiabatic ramp the CUDA coefficients will differ from the CPU coefficients by the adiabatic factor — a global multiplicative error in every coefficient.

---

Minor additional discrepancies (bookkeeping/subsampling only):

• sampleMasses in the CUDA path accumulates raw p.mass (line 847 via u_d) while the CPU path accumulates mm (which includes the sign and the adiabatic factor). Only relevant for the covariance/subsampling output.
• cuda_initialize() unconditionally sets fullCovar = false, so the full covariance matrix is never computed on the GPU path.
• The nminx/nminy low-wavenumber cut-off applied in the CPU force loop (line 681) is absent from forceKernelSlab. This is only relevant when nminx or nminy are set to non-zero values.