Integrators

This page states the discrete equations each integmode solves, the step size, how a loss crossing is located inside a step, and what SIMPLE does when an implicit step fails. Usage: Running SIMPLE.

Step size

With \(\Delta\tau_\text{macro} = v_0\,\texttt{trace\_time}/(\texttt{ntimstep}-1)\) the microstep is

\[\Delta\tau = \frac{\Delta\tau_\text{macro}}{N_\tau}, \qquad N_\tau = \Bigl\lceil \frac{\Delta\tau_\text{macro}\,\texttt{npoiper2}}{2\pi R_0} \Bigr\rceil ,\]

so \(\Delta\tau \le 2\pi R_0/\texttt{npoiper2}\) and every macrostep holds an integer number of microsteps (params_init in src/params.f90; \(R_0\) from the VMEC or chartmap geometry). The symplectic schemes take \(h = \Delta\tau/\sqrt2\) in the time \(\bar t\) of Guiding-centre model. Losses are resolved within a microstep; confined fractions are recorded per macrostep.

Symplectic schemes

All schemes advance \(z = (r,\vartheta,\varphi,\bar p_\varphi)\) with the canonical equations (4), written with \(H' \equiv \partial_r\bar H/\partial_r\bar p_\vartheta\) and the normalised quantities of Guiding-centre model (bars omitted below). The implicit equations are multiplied by \(\partial_r p_\vartheta\) to avoid the division and are solved by Newton’s method with analytic Jacobians from the second field derivatives.

Scheme (integmode)	Order	Newton unknowns	Field evaluated at
Euler 1 (1)	1	\(r^{n+1}, p_\varphi^{n+1}\)	\((r^{n+1},\vartheta^n,\varphi^n,p_\varphi^{n+1})\)
Euler 2 (2)	1	\(r^{n+1},\vartheta^{n+1},\varphi^{n+1}\)	\((r^{n+1},\vartheta^{n+1},\varphi^{n+1},p_\varphi^n)\)
Midpoint (3)	2	\(z^{n+1}\) and \(r^{n+1/2}\) (5)	averaged angles and \(p_\varphi\), radius \(r^{n+1/2}\)
Gauss–Legendre (4–7)	\(2s\)	all \(s\) stage states (\(4s\))	the \(s = 1..4\) stages
Lobatto IIIA–IIIB (15)	4	\(4s-2 = 10\)	the \(s = 3\) stages

Euler 1 (default) is explicit in the angles and implicit in the momenta. sympl_euler1_residual in src/orbit_symplectic_euler1.f90 solves

\[\partial_r p_\vartheta\,(p_\vartheta - p_\vartheta^n) + h\,(\partial_\vartheta H\,\partial_r p_\vartheta - \partial_r H\,\partial_\vartheta p_\vartheta) = 0, \qquad \partial_r p_\vartheta\,(p_\varphi^{n+1} - p_\varphi^n) + h\,(\partial_\varphi H\,\partial_r p_\vartheta - \partial_r H\,\partial_\varphi p_\vartheta) = 0,\]

with all field quantities at \((r^{n+1},\vartheta^n,\varphi^n,p_\varphi^{n+1})\), then updates the angles explicitly, \(\vartheta^{n+1} = \vartheta^n + hH'\), \(\varphi^{n+1} = \varphi^n + h(v_\parallel - H' h_\vartheta)/h_\varphi\). After convergence the field values are extrapolated to first order in the last Newton correction instead of being re-evaluated. Euler 2 is the adjoint: \(p_\vartheta(r^{n+1},\vartheta^{n+1},\varphi^{n+1},p_\varphi^n) = p_\vartheta^n\) fixes \(r^{n+1}\), the angle equations are implicit, and \(p_\varphi\) is updated explicitly. Midpoint evaluates the field once at \((r^{n+1/2},\bar\vartheta,\bar\varphi,\bar p_\varphi)\) with averaged angles and momentum; the extra unknown \(r^{n+1/2}\) satisfies \(p_\vartheta(z^{n+1/2}) = p_\vartheta^n + \tfrac h2\dot p_\vartheta\), which together with the full-step equation makes \(p_\vartheta^{n+1/2}\) the average of both ends (f_midpoint_part1/part2 in src/orbit_symplectic.f90). Gauss–Legendre solves the stage equations \(z_k = z^n + h\sum_l a_{kl}\dot z_l\) for \((p_\vartheta,\vartheta,\varphi,p_\varphi)\) and updates them with the weights \(b_l\); the stored radius after the step is the last stage radius, and \(r^{n+1}\) is fixed implicitly by \(p_\vartheta^{n+1}\) at the next step (gauss_canfield_residual in src/orbit_rk_core.f90).

Newton stops when all residuals are below \(10^{-15}\) or all corrections below relerr times a reference scale (Euler 1: 1 for \(r\), \(10\,\psi_\text{tor}^\text{edge}/\bar\rho_0\) for \(p_\varphi\)), after at most 32 iterations. An outward correction past the loss boundary \(r_\text{max}\) is cut to 0.8 of the remaining distance; a step ending at \(r < 0\) is mapped through the axis, \((r,\vartheta)\to(-r,\vartheta+\pi)\); Euler and midpoint then re-evaluate the field at the new endpoint, because the old-chart increments are not a first-order extrapolation across the flip. The schemes conserve the canonical two-form to the Newton tolerance, so backward-error analysis bounds the long-time energy error for small steps in smooth regions (Hairer, Lubich and Wanner 2006); with early classification, loss computations cost about one order of magnitude less than with adaptive Runge–Kutta (Albert, Kasilov and Kernbichler 2020, JPP).

Adaptive Runge–Kutta (integmode = 0)

The Cash–Karp RK5(4) integrator of libneo (odeint_allroutines) solves the non-canonical guiding-centre equations in any coordinates, with relative tolerance relerr and \(w = v/v_0\) (velo_can in src/sub_alpha_lifetime_can.f90):

\[\frac{dx^i}{d\tau} = \frac{w\lambda\,h^{*i} + \dfrac{\rho_0 w^2(1-\lambda^2)}{2B\sqrt g}\,\varepsilon^{ijk}h_j\partial_k\ln B} {h^*_\parallel}, \qquad \frac{d\lambda}{d\tau} = -\frac{(1-\lambda^2)\,w}{2h^*_\parallel}\, \mathbf h^*\!\cdot\nabla\ln B, \qquad \frac{dw}{d\tau} = 0,\]

with \(\mathbf h^* = \mathbf h + w\lambda(\rho_0/B)\nabla\times\mathbf h\) and \(h^*_\parallel = 1 + w\lambda(\rho_0/B)\,\mathbf h\cdot\nabla\times \mathbf h\). It is not symplectic. Below \(s = 0.01\) the step switches to \((\sqrt s\cos\vartheta, \sqrt s\sin\vartheta, \varphi)\), which is regular at the axis (orbit_timestep_axis). The same equations serve as the first recovery method below.

Loss crossing inside a step

If a symplectic trial step ends beyond \(r_\text{max}\) (LCFS, outer chartmap boundary, or \(M_\text{vol}\) for SPECTRE), locate_symplectic_boundary repeats it from the accepted state with a fraction \(f\) of \(h\), chosen by a secant estimate of \(r(f) = r_\text{max}\) damped to 0.8 of the way, else by bisection; trials outside the domain become the upper bracket. It converges when the bracket width is at most boundary_event_fraction_tolerance and \(|r - r_\text{max}|\) at most boundary_event_radial_tolerance (defaults \(\max(10^{-12}, 10\,\texttt{relerr})\), \(\max(10^{-10}, 10\,\texttt{relerr})\)), within 64 trials. The loss time is \((n+f)\Delta\tau\); residual and bracket width go to orbit_exit_code.dat.

Failure recovery

With symplectic_newton_warning_mode = .true. (default) a failed microstep, typically near the axis, goes through these stages (src/orbit_symplectic.f90, src/simple_main.f90):

1. At the iteration limit, accept a final correction of at most 100 (SPECTRE: 10) times the convergence threshold.

2. Retry as two half steps, recursively to depth 8.

3. RK equations above from the last accepted state, at most 10,000 steps; then, at reference radius \(< 0.1\) with \(|B^*_\parallel/B| \ge 0.2\), with axis floor \(10^{-3}\); then a toroidally regularised RK step (at most 100,000 steps); then, below reference radius 0.02 with a chartmap field, a Boris full-orbit bridge. Each must keep \(v/v_0\) within 0.5 %.

4. Resume the symplectic scheme (reseed_sympl) after 256 consecutive RK steps with \(|B^*_\parallel/B| \ge 0.3\) and reference radius \(\ge 0.12\) (recovery started below 0.1) or \(\le 0.98\) (started above 0.98), otherwise after \(2^{7+k}\) such steps for the \(k\)-th recovery (\(k\le9\)).

5. Otherwise roll back and skip the microstep; a second consecutive failure ends the marker with exit code 101–105, or 4 if collisionless at \(s < 0.01\) (Outputs).

With the warning mode off, the first failure ends the marker. Stages 2–5 are not a fixed symplectic map.

References

• C. G. Albert, S. V. Kasilov, W. Kernbichler, Symplectic integration with non-canonical quadrature for guiding-center orbits in magnetic confinement devices, J. Comput. Phys. 403, 109065 (2020), https://doi.org/10.1016/j.jcp.2019.109065

• C. G. Albert, S. V. Kasilov, W. Kernbichler, Accelerated methods for direct computation of fusion alpha particle losses within stellarator optimization, J. Plasma Phys. 86, 815860201 (2020), https://doi.org/10.1017/S0022377820000203

• E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration, 2nd ed., Springer (2006), https://doi.org/10.1007/3-540-30666-8