Add opt-in symplectic near-axis pseudo-Cartesian step#401
Open
krystophny wants to merge 9 commits into
Open
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
krystophny
added a commit
that referenced
this pull request
Jun 19, 2026
## 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.
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.
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