Pvar-dissip-struct: velocity-dependent variational equations in C (+ Chandrasekhar Jacobians)

[jobovy]

Summary

Extends the 3D variational equations to velocity-dependent (dissipative) forces, per the design doc: the variational Jacobian becomes the general $J = [[0, I], [\partial F/\partial x,\ \partial F/\partial v]]$ (non-symmetric position block + velocity block), with the conservative case as the special case $\partial F/\partial v = 0$.

C machinery

• struct potentialArg gains RectForceJacobian(t, q, jac_x, jac_v, ...) filling the rectangular 3×3 $\partial F/\partial x$ and $\partial F/\partial v$ blocks (NULL for conservative; NULL-init centralized in init_potentialArgs, verified used by every parser).
• evalRectDeriv_dxdv generalized to $d(\delta v)/dt = (K + \Sigma,\partial F/\partial x)\delta x + (\Sigma,\partial F/\partial v)\delta v$; new calcRforceConservative/calcphitorqueConservative keep friction out of the cylindrical→Cartesian K-assembly terms. Conservative path verified bit-identical (pre/post builds compared on MN/LogHalo/DehnenBar × 3 integrators).
• ChandrasekharDynamicalFrictionForce Jacobians fully analytic (F = A(x,|v|)·v ⇒ closed forms incl. the $d\ln\Lambda/dv$ term and clamped-spline σ′ — nothing dropped), except one documented term: $\partial\rho/\partial x$ uses central FD of calcDensity (no density-gradient C interface; ~1e-8 accurate).

Python gating

• DissipativeForce.hasC_dxdv3d = False default; Chandrasekhar → True. FDMDynamicalFrictionForce explicitly overridden back to False — it subclasses Chandrasekhar and would otherwise silently use wrong Jacobians (hole caught and closed).
• Pure-Python dxdv methods raise NotImplementedError for dissipative pots (the Python _EOM_dxdv is conservative-only) — loud failure instead of silently wrong results.

Validation (all committed as tests)

Check	Result
Conservative regression (dxdv/liouville battery)	121/121, bit-identical
C Jacobian vs FD-of-C-force (18 entries, 3 lnΛ configs)	≤ 2.1e-8
FD-of-flow STM (MW2014 + friction, decaying orbit)	≤ 2.8e-5
Phase-volume law $\det M(t) = e^{\int \mathrm{tr}(\partial F/\partial v)dt'}$	≤ 3.1e-6 rel; $\det M = 0.518 < 1$ (friction contracts phase volume)
Dissipative excluded from det(M)=1/symplecticity registries	tripwire test

Flagged numpy-path change (standing rule)

CompositePotential now aggregates hasC_dxdv3d (it didn't) — so conservative list inputs like MWPotential2014 previously fell back to odeint silently; they now take the C path (outputs differ 1e-11–4e-8 = odeint-vs-C, both correct). Note: #926 contains the same one-line fix independently — whichever merges second rebases and dedupes (trivial).

Follow-ups (separate PRs)

FDM Jacobians; planar (4D) dissipative dxdv; NonInertialFrameForce.

🤖 Generated with Claude Code

[codecov]

Codecov Report

❌ Patch coverage is 0% with 142 lines in your changes missing coverage. Please review.
✅ Project coverage is 15.85%. Comparing base (6d48e49) to head (2e23404).

Files with missing lines	Patch %	Lines
...ential_c_ext/ChandrasekharDynamicalFrictionForce.c	0.00%	76 Missing ⚠️
galpy/potential/potential_c_ext/galpy_potentials.c	0.00%	31 Missing ⚠️
galpy/orbit/orbit_c_ext/integrateFullOrbit.c	0.00%	30 Missing ⚠️
galpy/orbit/integrateFullOrbit.py	0.00%	2 Missing ⚠️
...y/potential/ChandrasekharDynamicalFrictionForce.py	0.00%	1 Missing ⚠️
galpy/potential/DissipativeForce.py	0.00%	1 Missing ⚠️
galpy/potential/FDMDynamicalFrictionForce.py	0.00%	1 Missing ⚠️

❗ There is a different number of reports uploaded between BASE (6d48e49) and HEAD (2e23404). Click for more details.

HEAD has 12 uploads less than BASE

Flag	BASE (6d48e49)	HEAD (2e23404)
	13	1

Additional details and impacted files

@@                   Coverage Diff                   @@
##           feat/variational3d     #933       +/-   ##
=======================================================
- Coverage               99.92%   15.85%   -84.08%     
=======================================================
  Files                     225      225               
  Lines                   36207    36319      +112     
  Branches                  785      787        +2     
=======================================================
- Hits                    36180     5757    -30423     
- Misses                     27    30562    +30535     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[jobovy]

Review addressed (branch reworked in place, force-pushed):

1. Rename: RectForceJacobian → RectDissipativeForceJacobian everywhere (struct member, calcRectDissipativeForceJacobian aggregator, per-force functions, wiring) — the name now signals it exists only for the velocity-dependent (DissipativeForce-subclass) forces; conservative potentials leave it NULL (documented at the struct member).
2. Direct C accessors removed: evalRectForceJacobian/evalRectForceVelocity and the ctypes per-entry unit test are gone. Validation is now purely through regular galpy orbit use: the FD-of-flow STM test, and the phase-volume law with tr(∂F/∂v) computed by central FD of the Python forces — which strengthens the test (the C STM's det M(t) is checked against a fully independent Python-derived trace; agreement still ≤1.6e-6).
3. evalRectDeriv_dxdv simplified to the galpy aggregator pattern: one straight-line RHS, with the dissipative blocks supplied by the NULL-safe aggregator (zero-filled when absent) — no parallel code paths. Conservative path re-verified bit-identical vs a fresh feat/variational3d build (18/18 byte-identical: 3 potentials × 3 integrators × 2 seeds).
4. Stale force enumerations removed from all warnings/docstrings (grep-gated to zero).
5. Per-force renames cascaded to Pvar-dissip: FDM dynamical-friction rectangular Jacobians in C #934/Pvar-dissip: NonInertialFrameForce rectangular Jacobians in C #935 (rebased + reworked identically).

Batteries: #933 126 passed (the −1 vs before is exactly the deleted ctypes test), #934 128, #935 130 + test_noninertial 53; all physics numbers reproduce (phase-volume 1.59e-6/det M=0.677 for FDM; det M−1 ≤ 9.2e-11 for the rotating frame).

[jobovy]

Round-2 review addressed (amended in place):

• calcRforceConservative/calcphitorqueConservative removed — calcRforce/calczforce/calcphitorque now take a final include_dissipative int (the existing macros supply 1, so every current call site is unchanged; the variational RHS calls the functions directly with 0). Conservative path re-verified bit-identical vs the base branch (hash-equal dxdv outputs, 3 potentials × 3 integrators × 2 seeds).
• Struct-member comment shortened to 5 lines; HISTORY entry shortened to 6.
• The 11 uncovered Chandrasekhar Jacobian lines are now covered by orbit-level tests only (test_chandrasekhar_dxdv_fd_of_flow_branches, parametrized over const-lnΛ / rhm-Coulomb-log / minr-gate / σ-clamp configs, each with a non-vacuity guard asserting the branch is genuinely selected; mutation-validated 1:1 — each locally-perturbed branch term is caught by exactly its test). One physics note: the minr-gate test keeps the orbit entirely inside minr — a crossing orbit's flow has a saltation term at the force discontinuity that the gate's zero Jacobian correctly omits, so FD-of-flow is only well-posed in the non-crossing regime (documented in the test).
• FDM hasC_dxdv3d=False override: confirmed removed in Pvar-dissip: FDM dynamical-friction rectangular Jacobians in C #934 (its Jacobian is wired there).
  Battery: 130 passed (126 + 4 branch tests); cascades re-run green (Pvar-dissip: FDM dynamical-friction rectangular Jacobians in C #934: 8, Pvar-dissip: NonInertialFrameForce rectangular Jacobians in C #935: 5+53).

[jobovy]

Round-3 refinements (amended; #934/#935 deliberately NOT rebased pending agreement here):

• The include_dissipative int now sits before the optional velocities, so every caller uses the plain calcRforce(...) macro — no parenthesized direct calls. The macro defaults it to 1 for the existing arities; the six existing velocity call sites pass 1 explicitly; the variational K-assembly is now simply calcRforce(R,z,phi,t,nargs,potentialArgs,0,0.,0.,0.).
• Conditional flattened to if (!requiresVelocity) add Rforce; else if (include_dissipative) add RforceVelocity; — note the exact form if (include_dissipative && requiresVelocity) with a plain else would wrongly route a dissipative component to its (unset) conservative Rforce pointer when the flag is 0, hence the else-if form.
• This rebase also dropped the duplicated CompositePotential.hasC_dxdv3d fix (now inherited from Pvar-lyapunov: largest Lyapunov exponent from the variational equations (Orbit.lyapunov) #926's merge) and merged the HISTORY entries.

Conservative path bit-identical (hash-equal vs base: 51a75aa8…); battery 130 passed.