Add time-dependent support to SCFPotential#1058
Merged
Merged
Conversation
Allow each expansion coefficient Acos/Asin to be a 1D function of time, the SCF analogue of the time-dependent MultipoleExpansionPotential (#847). - __init__ gains a tgrid keyword; Acos/Asin may be callables f(t) returning the (N,L,M) coefficient array or precomputed (Nt,N,L,M) arrays sampled on tgrid. Coefficients are interpolated in time with a cubic spline (CubicSpline, whose PPoly coefficients are passed to C for exact Python/C parity). - from_density gains a tgrid keyword and detects time-dependent densities (a t argument, or a Potential instance), computing coefficients at each time. Static behavior (no tgrid) is unchanged. - Fast C access: the per-coefficient cubic-in-time PPoly blocks are serialized to C, which reconstructs the interpolated coefficients (Horner in t) before the usual summation, for the potential, forces, second derivatives, full 3D Hessian (variational integration), and density. The C cache now keys on t. C and Python agree to machine/integrator precision; new tests in test_scf.py cover callable/array init, reduction to static, isNonAxi detection, from_density time dependence, Python<->C parity (orbits, planar/3D dxdv, dynamical-friction density, direct C evaluators, beyond-tgrid extrapolation), amplitude scaling, array-t broadcasting, and error/warning paths. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01DCoE6U3bdzXnxiNy2QtBps
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #1058 +/- ##
========================================
Coverage 99.93% 99.93%
========================================
Files 227 227
Lines 38188 38509 +321
Branches 840 849 +9
========================================
+ Hits 38164 38485 +321
Misses 24 24 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Follow the same convention as the time-dependent MultipoleExpansionPotential: add a weakly time-dependent, non-axisymmetric SCF test potential to the energy/Jacobi orbit-conservation test (test_orbit via conftest), the actionAngleStaeckel C test (potential evaluation in C), and the dynamical friction C test (density evaluation in C). Generalize the "WeaklyTDMultipole" special-case string to "WeaklyTD" so both the Multipole and SCF mocks are handled. Remove the test_scf.py tests that relied on interpRZPotential's eval_potential_c/eval_force_c/eval_2ndderiv_c helpers. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01DCoE6U3bdzXnxiNy2QtBps
The O(N*L*M) summation loops are inlined into the hot evaluation functions; threading time-dependence through those functions perturbed the compiler's code generation for the inlined non-axisymmetric loop, giving a ~10% slowdown of static non-axisymmetric SCF orbit integration. Force the summation helpers out-of-line (SCF_NOINLINE) so their code generation is isolated from the time-dependence handling, and take the static fast path (use the parsed coefficients directly, skipping the interpolation machinery) explicitly in each evaluation function. Static orbit-integration speed now matches the pre-change code (verified for axisymmetric and non-axisymmetric SCF); time-dependent evaluation is unaffected and C/Python parity is preserved to machine precision. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01DCoE6U3bdzXnxiNy2QtBps
Add a "Time-dependent SCFPotential" section to the potentials/scf_and_multipole notebook, mirroring the time-dependent MultipoleExpansionPotential section (with the same caveats: astropy-Quantity densities unsupported on the time-dependent path; plus a note that the SCF coefficient computation is more expensive per time step than the multipole's). It builds a time-dependent SCFPotential and MultipoleExpansionPotential from the same rotating m=2 bar density and compares the potential and an orbit in both. Executed locally with papermill in a CI-matching environment. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01DCoE6U3bdzXnxiNy2QtBps
Address review feedback on the time-dependent SCFPotential section:
- Overplot a high-order reference ("truth") and add a residual panel, so the
small SCF-vs-multipole difference is shown for what it is: shared angular
truncation error (the theta-independent cos(2phi) bar projects onto
l=2,4,6,..., truncated at finite L), with both expansions converging to the
reference. Bump the compared degree to L=6.
- Add a Tip clarifying that, unlike the multipole, SCFPotential.from_density
treats the density as time-dependent only when a tgrid is passed (so passing
e.g. pot.dens without a tgrid builds the static potential -- no footgun).
- Add tight_layout() to the figures.
Re-executed locally with papermill.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01DCoE6U3bdzXnxiNy2QtBps
Follow up on review: make the time-dependent SCFPotential example show off the time dependence and use the exact wrapped potential as the reference, while keeping the comparison to MultipoleExpansionPotential. - Plot the potential at several times so the m=2 bar is seen to rotate (the previous t=0 snapshot was indistinguishable from the static case). - Keep the SCF-vs-multipole comparison, and separately validate each time-dependent expansion against the exact rigid rotation of its own static expansion (SolidBodyRotationWrapperPotential): each reproduces it to ~1e-5 (the cubic-spline-in-time interpolation error), so the ~1e-3 SCF-vs-multipole offset is shown to be ordinary expansion truncation, not a time-dependence artifact. - The orbit comparison now integrates within the tgrid range (a tgrid spanning only one period previously forced the time-dependent SCF to extrapolate and diverge over the multi-period integration). - tight_layout() on both figures. Re-executed locally with papermill. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01DCoE6U3bdzXnxiNy2QtBps
Vectorize the time-dependent ``SCFPotential.from_density`` coefficient
computation over the ``tgrid`` axis, mirroring the time-dependent
``MultipoleExpansionPotential``: the (time-independent) basis functions
are evaluated once and the density is sampled at all times at once, in a
single Gauss-Legendre quadrature that carries a leading time axis. This
replaces the previous per-timestep loop and gives a large speed-up
(~35x for N=10, L=4 over 41 steps) at machine-precision-identical
coefficients. A density that cannot be evaluated on an array of times
(and one without a ``t`` argument) transparently falls back to the
per-timestep loop / constant-in-time broadcast.
New helpers ``_scf_compute_coeffs_{spherical,axi}_timedep`` and
``_scf_compute_coeffs_timedep`` produce the same ``(Nt,N,L,M)`` coefficient
arrays as before, so the C serialization and evaluation are unchanged.
Tests: add coverage for the vectorized spherical/axi/general paths, the
numOfParam-detection fallbacks, the explicit ``*_order`` overrides, the
non-vectorizable fallback, and the constant-in-time non-axisymmetric
broadcast.
Notebook (scf_and_multipole): reword the build caveat to reflect the now
time-vectorized setup (the remaining consideration is memory, not
per-timestep cost); drop the inaccurate "chaotic bar orbit" wording; and
regenerate outputs.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01DCoE6U3bdzXnxiNy2QtBps
The time-vectorized `from_density` quadrature accumulates a coefficient array with a leading time axis, so its working set grows linearly with the number of time steps. Process the `tgrid` in memory-bounded batches internally (`_batched_timedep` / `_timedep_batch_size`, gated by the module-level `_TIMEDEP_BATCH_BYTES` budget) so peak memory stays bounded for an arbitrarily large `tgrid`. Batching over the independent time axis is exact — results are bitwise-identical to the single-shot build (verified 0 diff) — and there is no public API change. Measured at Nt=12000 (N=10, L=4): peak build memory drops from ~314 MB to ~151 MB (and, from better cache locality, the build is also faster: ~20 s vs ~59 s). The non-vectorizable-density fallback is unaffected (it raises on the first batch and propagates as before). Tests: add `test_tdep_from_density_time_batching`, which forces a tiny budget so the build uses several batches and checks the batched result equals the single-shot build exactly (covering both the Asin-present and Asin=None concatenation branches). Notebook (scf_and_multipole): correct the build caveat — the earlier wording implied a quadrature-grid-scale memory blow-up, which does not happen; memory is now internally batched and bounded, and the cost that grows with the number of time steps is build time, not memory. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01DCoE6U3bdzXnxiNy2QtBps
This was referenced Jul 4, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds time dependence to
SCFPotential, allowing each expansion coefficientAcos/Asinto be a 1D function of time — the SCF analogue of thetime-dependent
MultipoleExpansionPotentialadded in #847.What
SCFPotential.__init__gains atgridkeyword. When given,Acos/Asinare interpreted as time-dependent coefficients, supplied either as
f(t)returning the(N, L, M)coefficient array, or(Nt, N, L, M)arrays sampled ontgrid.The coefficients are sampled on
tgridand interpolated in time with a cubicspline (
scipy.interpolate.CubicSpline, aPPolysubclass whose coefficientsare handed to C directly for exact Python/C parity).
SCFPotential.from_densitygains atgridkeyword and detects atime-dependent density (a
targument, or a galpyPotentialinstance),computing the coefficients at each time in
tgrid. As with the time-dependentMultipole, astropy-
Quantitydensities are not supported on this path; passthe density in internal units. Static behavior (no
tgrid) is unchanged, soexisting
from_density(pot.dens, ...)usage is untouched.Fast C access. The per-coefficient cubic-in-time PPoly blocks are
serialized to C, and the C code reconstructs the interpolated coefficients at
the current time (Horner in
t) before the usual Hernquist–Ostrikersummation. All C entry points are time-aware: potential, forces, second
derivatives, the full 3D Hessian (for
integrate_dxdvvariationalintegration), and density (used by, e.g.,
ChandrasekharDynamicalFrictionForce).The force/derivative cache now keys on
tas well.Parity
C and Python agree to machine / integrator precision (verified for the
potential, forces, densities, planar and 3D variational integration, and time
values extrapolated beyond
tgrid).Tests
New tests in
tests/test_scf.pycover callable and array init, reduction to thestatic potential, grid-node exactness,
isNonAxidetection,from_density(callable-
t,Potential-instance, axisymmetric, and constant-in-time),Python↔C parity for orbits / planar+3D
integrate_dxdv/ dynamical-frictiondensity / direct C evaluators / beyond-
tgridextrapolation, amplitude scaling(deepcopy), array-
tbroadcasting, and the error/warning paths. All new Pythonand C lines are exercised.
🤖 Generated with Claude Code
https://claude.ai/code/session_01DCoE6U3bdzXnxiNy2QtBps