Add opt-in symplectic near-axis pseudo-Cartesian step#401
Conversation
krystophny
changed the title
axis_pcart namelist flag (default off) and the orbit_symplectic_axis_pcart module documenting the planned fix for #398: keep the canonical Boozer residual (symplecticity preserved) and solve the implicit Newton in (X,Y) with field derivatives evaluated natively there, where they stay finite across the axis. Cites Pfefferle et al. arXiv:1412.5464. Full (X,Y) field-derivative evaluator and in-(X,Y) Newton to follow, gated on the idx-12 energy test.
When axis_pcart is on, advance the substep in (X,Y)=(sqrt(s)cos th, sqrt(s)sin th) using the canonical guiding-centre velocity, which is regular across the axis, so no negative-s floor and no (s,theta)->(|s|,theta+pi) reflection are needed. Gate (idx-12 axis-grazing orbit, internal Boozer): max|dH/H| drops from 1.6e-2 (flux-chart Euler1) to 3.9e-4, with the per-near-axis-pass energy random-walk eliminated and the orbit bounded. The current substep is an explicit RK4 in (X,Y) (proof of concept, not yet strictly symplectic and global/slow); the refinement is an implicit symmetric step used only near the axis. Cites Pfefferle arXiv:1412.5464.
Replace the explicit (X,Y) RK4 proof-of-concept with the production form: keep the flux-chart Euler1 residual (same symplectic map, no chart seam) and solve it in the signed radius q (s=q^2) near the axis, with the normalized residual (finite) and a finite-difference Jacobian (no divergent second derivatives). q<0 is the axis crossing (theta+pi); no floor, no reflection. Set pthold at step start. Gate (idx-12, smax=0.1): max|dH/H| 1.6e-2 -> 8.3e-4, orbit bounded. Pfefferle arXiv:1412.5464.
Replace the finite-difference Jacobian in the signed-radius (q,pphi) Newton with the exact 2x2 Jacobian built from the field's second derivatives (get_derivatives2: d2H, d2pth) by quotient and product rule. dr/dq = dr/ds * 2q is finite at the axis because dr/ds ~ s^(-1/2) cancels against 2q ~ s^(1/2). Same Euler1 map and root as the flux chart, so no chart seam. Idx-12 grazing orbit, smax=0.1, 0.011 s: max|dH/H| 1.0e-3, bounded (flat across quarters), vs 1.6e-2 random-walk in the flux chart that ejects the particle. The (X,Y) implicit-midpoint substep seams at the chart boundary (2.8e-3) because it is a different map; the signed-q Euler1 step does not.
The signed-radius Newton in axis_pcart_step exited at maxit silently and proceeded with the last iterate, recording nothing; the flux-chart newton1 increments EVT_NEWTON1_MAXIT on the same condition. Flag convergence and count non-convergence the same way, so runs report near-axis Newton failures instead of hiding them.
krystophny
force-pushed
the
feat/sympl-axis-pseudocartesian
branch
from
8753d68 to
ecffc00
Compare
The implicit Euler1 near-axis Newton stalls at the axis. Symplectic Euler freezes theta during the implicit momentum solve, so the radial solve runs along a fixed ray through the sqrt(s) coordinate singularity; the negative-r floor then oscillates to maxit on crossing substeps. Midpoint and Gauss collocations reduce the floor oscillation but still commit the chart switch near the axis. Advance the optional near-axis shell inside s < axis_pcart_smax with an explicit RK4 of the guiding-centre EOM in (X,Y) = (sqrt(s) cos theta, sqrt(s) sin theta). Then s = X^2 + Y^2 >= 0 and the step needs no negative-s reflection. This path is opt-in and is not a final symplectic fix. It stays stacked after the RK-only correction so it can be reviewed and benchmarked separately. Refs #398.
krystophny
force-pushed
the
feat/sympl-axis-pseudocartesian
branch
from
ecffc00 to
bcbb020
Compare
krystophny
changed the title
## Summary Fix the RK near-axis pseudo-Cartesian transform used by `orbit_timestep_axis`. The old path entered the near-axis chart as `(s cos theta, s sin theta)` but treated it as `(sqrt(s) cos theta, sqrt(s) sin theta)` when evaluating velocity and returning to flux coordinates. That made the axis an artificial sink for some mirror-trapped particles. This PR makes the chart consistent: - enter with `(sqrt(s) cos theta, sqrt(s) sin theta)` - recover `s = x^2 + y^2` - use the matching `0.5` chain-rule factor in `velo_axis` - keep the floor only for field evaluation, not as the physical radial mapping This is the RK-only extraction from PR #401. The symplectic near-axis experiment stays stacked there. ## Verification `fo check` ```text {"ok":true,"stage":"done","target":"","summary":"build and tests passed","hint":"","rerun":"","log_path":"","elapsed_s":2.401,"modules":469,"cached":0,"changed":469} ```
…ocartesian # Conflicts: # src/CMakeLists.txt # src/orbit_symplectic.f90 # src/params.f90
test_axis_pcart: on the QA Boozer field the default flux chart hits 23 near-axis newton1_maxit and 23 r_negative faults at integmode=1; axis_pcart drives both to zero and loses no particle the default kept. Remove the dead si%z(1)>1 recheck (axis_pcart_step returns ierr=1 before updating z on overshoot). State the integmode=1 scope in the module doc.
![slopqueue[bot]](https://avatars.githubusercontent.com/u/48486978?s=60&v=4){kind=small}
There was a problem hiding this comment.
ORCHESTRATE_REVIEW_STATUS: approve
Here's my final review:
PR #401 Review: Opt-in near-axis pseudo-Cartesian step for symplectic integrators
Verdict: Approve. The physics is correct, the test is well-designed, and the scope is minimal and opt-in.
Summary
PR #401 adds an opt-in near-axis pseudo-Cartesian (X, Y) RK4 step for the symplectic Euler1 integrator (integmode=1, the default guiding-centre path). In flux coordinates s = ρ², the m=1 harmonic diverges as s^(-3/2) at the axis, causing the Newton solve to fail and the (s,θ)→(|s|,θ+π) chart switch to inject non-conservative energy. The fix: inside a thin shell s < smax, advance the substep with explicit RK4 in pseudo-Cartesian coordinates where the field is smooth and single-valued across the axis. RK4 is not symplectic, but the shell is thin and the energy error stays below the bulk symplectic level.
Physics verification
pcart_velocity correctly derives the guiding-centre equations in (X, Y, φ, pφ) from the non-canonical Hamiltonian structure. The key equation is sdot (ds/dτ), which uses the chain rule on the non-canonical Hamilton equation dpth/dτ = -∂H/∂θ + dθ/dτ·∂pth/∂θ. The pthd term is not redundant — it's the non-canonical correction to dpth/dτ. The derivation is correct and matches the equations in the existing symplectic integrator (orbit_symplectic_euler1.f90 residual).
The rho = sqrt(max(s, 1.0e-30_dp)) floor prevents division by zero at the axis. The snew > 1.0 check after the RK4 step is consistent with the x(1) > 1.0 check in the Euler1 path. Field state is correctly updated at the end of axis_pcart_step via eval_field/get_derivatives, matching the Euler1 post-Newton update.
Test verification
The test is well-designed:
- Runs the default flux chart first (axis_pcart=False) to generate
start.datwith deterministic particles - Reuses
start.datwith axis_pcart=True (startmode=2) for identical starting conditions - Verifies the default run actually hits near-axis faults (
EVT_R_NEGATIVE,EVT_NEWTON1_MAXIT) — proves the test exercises the code path - Verifies axis_pcart drives both faults to zero
- Verifies no particles are lost that the default kept (on_only = 0)
Findings
1. (medium) OpenACC device copy not updated by set_axis_pcart
!$acc declare copyin at module level copies initial values to the device. set_axis_pcart updates axis_pcart_enabled/axis_pcart_smax on host only — no !$acc update device. Currently harmless: the CPU path has no OpenACC regions, and the GPU path (simple_gpu.f90) doesn't use this module. But if OpenACC offload is added to the CPU path later, the device copy will be stale. Consider adding !$acc update device in set_axis_pcart, or a comment noting the assumption.
2. (low) OpenACC directives in new module are dead code
!$acc routine seq on pcart_velocity/axis_pcart_step and !$acc declare copyin on module variables are not exercised by any current code path (CPU path has no OpenACC regions; GPU path doesn't use this module). Not harmful, but misleading — they suggest the code is device-ready when it isn't fully wired.
3. (info) sdot derivation is correct but non-obvious
The pthd term computes the non-canonical Hamilton equation dpth/dτ = -∂H/∂θ + dθ/dτ·∂pth/∂θ, and the chain-rule cancellation of the ∂pth/∂θ·dθ/dτ terms is the correct derivation of ds/dτ. This is not redundant — it's the non-canonical structure. A brief comment explaining this would help future readers. (Note: the existing RK45 path in orbit_symplectic_quasi.f90:530 uses an approximation zdot(1) = -(∂H/∂θ - (hθ/hφ)·∂H/∂φ) / ∂pth/∂s that drops a dpth(3)·vpar/hφ term; sdot here is more accurate.)
4. (info) Test checks behavioral outcomes, not energy conservation
The test verifies faults → 0 and no lost particles, but doesn't check energy conservation. This is a documented trade-off: RK4 is not symplectic, but the shell is thin and the energy error is expected to be small. Acceptable for an opt-in experimental feature. If energy conservation becomes a concern, a dedicated test can be added later.
Scope and risk
- Opt-in: default
axis_pcart = .False.; only wired into integmode=1 (default GC path) - Other symplectic modes: unaffected (Lobatto, Gauss, etc. still use flux chart throughout)
- GPU path: has its own separate implementation (
simple_gpu.f90), does NOT use this module — correct, not a gap - Stacked on PR #402 (merged): dependency satisfied
- CI: all green (Build and Test, Build Documentation, CGAL Wall Unit Tests, Deploy to GitHub Pages)
- Diff size: small, focused on the new module + wiring + test
1 participant
Summary
Add an opt-in near-axis shell for the symplectic integrator.
This PR is stacked on PR #402. PR #402 fixes the RK pseudo-Cartesian mapping in
orbit_timestep_axis; this PR keeps the separate symplectic experiment isolated for review.The new path is disabled by default:
axis_pcart = .false.axis_pcart_smax = 0.01d0When enabled, symplectic steps with
s < axis_pcart_smaxadvance through a regular pseudo-Cartesian near-axis coordinate,(X,Y) = (sqrt(s) cos(theta), sqrt(s) sin(theta)), then map back withs = X^2 + Y^2.This avoids the flux-chart floor/reflection behavior near the magnetic axis. It is still an opt-in experiment, not the default production path.
Verification
fo check