Added coefficient-origin support to exp_fields.cc

.github/workflows/benchmarks.yml

@@ -43,7 +43,7 @@ jobs:
       github.event.pull_request.draft == false }}
     env:
       # Run tests against a tagged EXP version, unless this is the devel branch or a PR into devel, in which case test against EXP devel
-      EXP_REF: ${{ (github.ref == 'refs/heads/devel' || (github.event_name == 'pull_request' && github.event.pull_request.base.ref == 'devel')) && 'devel' || 'v7.9.1' }}
+      EXP_REF: ${{ (github.ref == 'refs/heads/devel' || (github.event_name == 'pull_request' && github.event.pull_request.base.ref == 'devel')) && 'devel' || 'v7.10.0' }}
     steps:
       - uses: actions/checkout@v6
         with:

.github/workflows/tests.yml

@@ -40,7 +40,7 @@ jobs:
       !contains(needs.check_skip_flags.outputs.head-commit-message, '[skip tests]') }}
     env:
       # Run tests against a tagged EXP version, unless this is the devel branch or a PR into devel, in which case test against EXP devel
-      EXP_REF: ${{ (github.ref == 'refs/heads/devel' || (github.event_name == 'pull_request' && github.event.pull_request.base.ref == 'devel')) && 'devel' || 'v7.9.1' }}
+      EXP_REF: ${{ (github.ref == 'refs/heads/devel' || (github.event_name == 'pull_request' && github.event.pull_request.base.ref == 'devel')) && 'devel' || 'v7.10.0' }}
       MACOSX_DEPLOYMENT_TARGET: "15.0"
     strategy:
       fail-fast: true

CHANGES.rst

@@ -26,6 +26,10 @@ Bug fixes
   This was caused by floating-point differences between the requested time and the
   limited-precision snapshot times stored in the coefficient file.
 
+- EXP: fixes a bug where the center set by calling ``setCoefCenter()`` on a pyEXP
+  ``Coefs`` object was ignored. Gala now requires EXP >= v7.10.0 when building with EXP
+  support.
+
 API changes
 -----------
 

docs/tutorials/exp.rst

@@ -27,7 +27,7 @@ is the best place to read about how to build EXP. Gala doesn't have any special
 instructions for the EXP build, except that the user must actually "install" EXP,
 rather than just build it. This is demonstrated below.
 
-Gala is compatible with EXP version >= 7.9.1. If you encounter build issues, double
+Gala is compatible with EXP version >= 7.10.0. If you encounter build issues, double
 check the EXP version.
 
 To install EXP's dependencies, here is one recipe that we have found to work on Ubuntu 24.04::

src/gala/potential/potential/builtin/exp_fields.cc

@@ -294,7 +294,7 @@ double exp_value(double t, double *pars, double *q, int n_dim, void* state) {
   // Get the field quantities
   // TODO: ask Martin/Mike for a way to compute only the potential - we're wasting
   // computation time here by computing all fields
-  auto field = exp_state->basis->getFields(q[0], q[1], q[2]);
+  auto field = exp_state->basis->getFieldsOrigin(q[0], q[1], q[2]);
 
   return field[5];
 }
@@ -315,6 +315,7 @@ void exp_gradient(double t, double *__restrict__ pars, double *__restrict__ q_in
   Eigen::Map<Eigen::VectorXd> eigen_y(q.y, N);
   Eigen::Map<Eigen::VectorXd> eigen_z(q.z, N);
 
+  // getAccel works with respect to the user's origin, same as getFieldsOrigin
   auto& allaccel = exp_state->basis->getAccel(eigen_x, eigen_y, eigen_z);
 
   for(size_t i = 0; i < N; i++) {
@@ -336,7 +337,7 @@ double exp_density(double t, double *pars, double *q, int n_dim, void* state) {
 
   // TODO: ask Martin/Mike for a way to compute only the density - we're wasting
   // computation time here by computing all fields
-  auto field = exp_state->basis->getFields(q[0], q[1], q[2]);
+  auto field = exp_state->basis->getFieldsOrigin(q[0], q[1], q[2]);
 
   return field[2];
 }
@@ -351,7 +352,7 @@ double exp_density(double t, double *pars, double *q, int n_dim, void* state) {
 //     );
 //   }
 
-//   auto field = exp_state->basis->getFields(q[0], q[1], q[2]);
+//   auto field = exp_state->basis->getFieldsOrigin(q[0], q[1], q[2]);
 
 //   for(int i=0; i<9; i++) {
 //     hess[i] += NAN;  // TODO: get hessian from EXP

tests/potential/potential/test_exp.py

@@ -343,6 +343,52 @@ def test_pyexp_unit_tests():
     assert u.allclose(pot_multi.tmax_exp, 2.0 * u.Gyr)
 
 
+@pytest.mark.skipif(not HAVE_PYEXP, reason="requires pyEXP")
+def test_pyexp_moving_centers():
+    """Test that setting per-snapshot centers shifts the EXP potential.
+
+    See scripts/moving_EXPPotential.py for the underlying pyEXP recipe.
+    """
+    units = EXP_UNITS
+
+    with open(EXP_CONFIG_FILE) as fp, chdir(EXP_CONFIG_FILE.parent):
+        basis = pyEXP.basis.Basis.factory(fp.read())
+
+    # Two independent Coefs objects -- setCoefCenter mutates in place
+    coefs_ref = pyEXP.coefs.Coefs.factory(str(EXP_MULTI_COEF_FILE))
+    coefs_moved = pyEXP.coefs.Coefs.factory(str(EXP_MULTI_COEF_FILE))
+
+    times = np.array(coefs_moved.Times())
+    t1, t2 = times[0], times[-1]
+
+    c1 = np.array([10.0, 0.0, 0.0])
+    c2 = np.array([0.0, 15.0, 0.0])
+    coefs_moved.getCoefStruct(t1).setCoefCenter(c1)
+    coefs_moved.getCoefStruct(t2).setCoefCenter(c2)
+
+    pot_ref = PyEXPPotential(basis=basis, coefs=coefs_ref, units=units)
+    pot_moved = PyEXPPotential(basis=basis, coefs=coefs_moved, units=units)
+
+    t1_q = t1 * units["time"]
+    t2_q = t2 * units["time"]
+    c1_q = c1 * units["length"]
+    c2_q = c2 * units["length"]
+    origin = [0.0, 0.0, 0.0] * units["length"]
+    fixed_pt = [8.0, 0.0, 0.0] * u.kpc
+
+    # At the moved center, the moved potential matches the ref at origin
+    assert u.allclose(pot_moved.energy(c1_q, t=t1_q), pot_ref.energy(origin, t=t1_q))
+    assert u.allclose(pot_moved.energy(c2_q, t=t2_q), pot_ref.energy(origin, t=t2_q))
+
+    # At a fixed point, the moved potential differs from the unmoved one
+    assert not u.allclose(
+        pot_moved.energy(fixed_pt, t=t1_q), pot_ref.energy(fixed_pt, t=t1_q)
+    )
+    assert not u.allclose(
+        pot_moved.energy(fixed_pt, t=t2_q), pot_ref.energy(fixed_pt, t=t2_q)
+    )
+
+
 def test_multi_different_snapshot_time_unit():
     pot_multi = EXPPotential(
         config_file=EXP_CONFIG_FILE,