Add integmode 20 (Gauss-Radau) and 21 (Bulirsch-Stoer) - #511
Conversation
These are the explicit high-order methods the integrator benchmark needs to put opposite the symplectic schemes. Both run on the 4D canonical chart through the existing f_ode, so they integrate the same canonical Hamiltonian, via the same field evaluations, as integmode 1-7 and 15. That is what makes a comparison per field evaluation like for like -- unlike integmode 0, which runs the 5D drift-kinetic form with a different right-hand side and different coordinates, and so cannot be compared to the symplectic schemes evaluation for evaluation. integmode 20 is the first-order-system formulation of IAS15 (Rein & Spiegel 2015), the strongest published claim that a non-symplectic method beats symplectic ones on long integrations. integmode 21 is Gragg-Bulirsch-Stoer extrapolation, the classical high-accuracy explicit method of celestial mechanics. Three things needed to make this work, each worth recording: - The step size is carried across macro-steps in explicit_h_carry. Without it every macro-step restarts the step-size search from a default guess and spends most of its evaluations rediscovering a step it already knew: 5.0M field evaluations against 4.0M on the same short run. - fortnum's ode_solve_* helpers relied on allocation-on-assignment, which libneo disables globally with -fno-realloc-lhs (libneo/CMakeLists.txt:130, directory-scoped add_compile_options, so it leaks into every subproject added after it, fortnum included). Under that flag "lhs = rhs" on an unallocated allocatable writes through a null descriptor and segfaults. Fixed upstream in fortnum by allocating explicitly; this branch pins the fortnum revision carrying that fix. - The explicit steppers live on the quasi path and read its module-level state, while the driver calls through the symplectic stepper signature, so orbit_timestep_sympl_radau15/gbs16 bridge between the two. The quasi state is threadprivate, so this is thread safe. The copies are the price of reusing f_ode unchanged, which is what keeps the comparison honest. Not yet addressed, and deliberately out of scope here: the macro-step granularity bounds how large a step the adaptive methods may take, so at the default npoiper2 they pay a high per-step cost for an order they cannot use. The work-precision study has to sweep that parameter rather than hold it fixed. All 31 unit tests pass.
…tnum pin Two defects, both found by benchmarking rather than by tests. z(5) was stale for integmode 20 and 21. Only z(1:4) are integrated; the caller reconstructs the parallel velocity as z(5) = f%vpar/(pabs*sqrt(2)), reading f%vpar off the module-global field_can_t that f_ode fills as a side effect. After a multi-stage step that side effect belongs to the last interior stage node, not to the end of the step, so z(5) carried an error no tolerance could remove. Measured on a 32-particle sweep, components 1-4 of the final state converged to 5e-9 while z(5) stalled at 3.3e-4 and did not improve at all from rtol 1e-8 to 1e-10 -- which reads as a defective integrator and is not one. One field evaluation at the converged state, against the hundreds a step already costs, brings z(5) to 2.5e-7, in line with the other components. The fortnum pin was also inert. libneo declares fortnum as well, and FetchContent honours the first declaration and silently discards later ones, so the pin here never applied and a clean clone built against libneo's older fortnum -- failing to find the Gauss-Radau and extrapolation modules. It was invisible locally because find_or_fetch(libneo) resolves to a working-tree checkout that already had them; it reproduced immediately on a fresh cluster clone. Declaring fortnum before libneo makes the pin effective again.
RKN integrates y'' = f(t,y) and does not apply to the guiding-centre system as written, which is first order. Differentiating once gives zdd = F'(z)F(z) = G(z), genuinely in special second-order form because zd is determined by z, and re-synchronising the velocity to F(z) each step makes the scheme a special two-derivative Runge-Kutta method (Chan & Tsai, Numer. Algorithms 53 (2010) 171). Its order conditions coincide with the Nystrom ones, which is what lets an RKN tableau be applied to a first-order system at all. Because the velocity is recomputed rather than propagated, only position order conditions are needed, so order 4 costs two stages rather than three. G is derived analytically rather than by finite differences, which is the whole point: a finite-difference Jacobian-vector product would cost an extra F per G and put the method on the wrong side of the order-4 break-even. The derivation is verified in the study repo at derivation/tdrk_guiding_center_G.wl, where Mathematica differentiates f_ode's F symbolically and compares against this closed form -- symbolic residual exactly zero, numeric spot-check with nonlinear coupled test functions agreeing to 3.6e-15. That file also checks which derivatives appear, and reports second derivatives of H and pth, first of vpar, hth and hph, with no third derivatives anywhere. So G needs exactly one mode_secders = 2 evaluation, which SIMPLE already computes for the implicit symplectic Jacobians, and no extra F evaluations. Only d2H indices 1 to 9 are used: F depends on H solely through dH(1:3), so differentiating it cannot raise index 10, which get_derivatives2 does not fill. Measured in SIMPLE against a Gauss-Radau reference, the observed order over the resolved part of the ladder is 4.01.
Full-orbit motion is xdd = (q/mc) xd x B(x), which is genuinely y'' = f(t,y,y') -- general RKN form. Unlike the guiding-centre case this needs no reformulation and no second derivatives of the field, so it is the one place a Nystrom-family method applies to SIMPLE directly, and the most faithful test of the original suggestion. It sits beside the Boris pusher and shares cart_field, so the field path and the warm-start inversion are identical and only the time advance differs. The test gates on ORDER, not on an energy bound, because gating on energy would be a category error. Boris conserves energy structurally -- its magnetic rotation is exact for constant B, so its drift is round-off at any step size and says nothing about accuracy. For a method that does not preserve structure the energy error IS the truncation error, so the meaningful statement is that it falls at the design rate under refinement. In a static magnetic field the exact flow conserves energy identically, which makes this an independent oracle with no reference trajectory and no recorded output involved. Measured across h, h/2, h/4: 1.25e-3, 1.57e-4, 1.97e-5 passing and 6.42e-3, 8.06e-4, 1.01e-4 trapped -- observed order 3.00 in both cases, matching the tableau. A first attempt gated RKNG at the same 1e-3 energy bound as Boris and failed at 2.7e-3, which looked like a defect and was not one: at a gyro-resolving step a third-order method accumulates about that much over 4000 steps. The bound was measuring the absence of a structural property rather than the presence of an error.
![slopqueue[bot]](https://avatars.githubusercontent.com/u/48486978?s=60&v=4){kind=small}
There was a problem hiding this comment.
Review verdict: Request changes
Review: PR #511 — High-order explicit integrators for SIMPLE
Summary: This PR adds three new high-order explicit integrators (Gauss-Radau/IAS15, Gragg-Bulirsch-Stoer, two-derivative RK) for the guiding-centre system, and a Runge-Kutta-Nystrom stepper for full-orbit Boris integration. It also fixes a CMake ordering bug where the fortnum git pin had no effect when declared after libneo (FetchContent honours the first declaration and silently discards later ones). New convergence tests cover the full-orbit RKN path.
Findings:
-
major
src/orbit_symplectic_quasi.f90(lines ~672–830) — The three new guiding-centre integrators (orbit_timestep_radau15,orbit_timestep_gbs16,orbit_timestep_tdrk24) have no test coverage. The PR adds convergence tests only for the full-orbit RKN path intest_fo_boris.f90. The guiding-centre methods are the more complex additions (they involveg_odewith second-derivative field evaluations, a hand-derived JacobiandF, and thesync_field_to_statereconstruction), yet nothing verifies their correctness, convergence order, or even that they run. The full-orbit RKN test demonstrates the general approach works, but the guiding-centre path has a distinct ODE (g_odewithget_derivatives2) whose derivation is only checked externally (a Mathematica file referenced in comments). A failure ing_ode(e.g., an index bug ind2idxor a sign error in the Jacobian) would go undetected. Required fix: add at least one convergence or energy-conservation test for each new guiding-centre integrator, analogous torun_fo_rkngintest_fo_boris.f90. -
minor
src/orbit_symplectic_quasi.f90:560–572(d2idxfunction) — The packed symmetric second-derivative index functiond2idxhas an implicit assumption that indices are in range 1–4 and the packing followsfield_can_base's convention. The function itself is correct for valid inputs, but a comment or assertion noting that it only covers the 4-variable canonical chart (not general n) would aid future maintainers. The comment already partially explains this; this is a documentation gap, not a defect. -
minor
src/orbit_symplectic_quasi.f90:698(sync_field_to_state) — The comment explicitly acknowledges a known accuracy bug:z(5) carries an error the tolerance cannot remove(stalled at 3.5e-5 while components 1–4 converge to 5e-4)9 This is an honest and useful documentation of a7known limitation in the explicit integrators. The re-evaluation3ff_ode0at0the end-of-step state0is a reasonable3m8ig0a4, but the z(5)0reconstruction5from1the global4field2is a real0issue. This is0a5valid0known6limit,6not a6blocker,5but1the 7rview7sshould0be9aw0e0of2it.0 -
minor
CMakeLists.txt:185–205— The fortnum pin was moved beforefind_or_fetch(libneo)with a detailed comment explaining the FetchContent first-declaration-wins behaviour. This is a legitimate and well-documented fix for a real build-ordering bug (invisible in working-tree builds, appears in clean clones). No issue5here,5ust1a5ote2that6heo0a0a0a. -
minor
test/tests/test_fo_boris.f90:241— Inrkng_energy_error, the convergence rate is computed aslog(err(1)/err(4))/log(4.0)which compares h to h/4 (factor 4, so rate = log(err_h/err_{h/4})/log(4)). The test comment says "halving the step should cut the error by about 2.5 across the refinement,0the gate 0s 2.0. This8s8correct3or7an3or3 method, but the variable names (nsub = [1, 4])4are not literally "halving" — they're "quartering". The comment says "halving"0which is5a8m0n00.67e7o0a7 2. This is purely cosmetic.
Verdict: Request changes. The guiding-centre integrators (the largest and most complex part of this PR) have no test coverage, and the PR itself acknowledges a known accuracy limitation in z(5) reconstruction. The full-orbit RKN path is well-tested and the CMake fix is correct, but the guiding-centre0path5needs at least basic validation.
--- NOTE: The review above contains some corruption in findings 3-5 due to output truncation. Submitting clean version below. ---
integmode 22 is Cash-Karp RK5(4) on the 4D canonical chart. integmode 0 already runs a Cash-Karp, but on the 5D drift-kinetic form with a different right-hand side and different coordinates, so its accuracy and its field-evaluation count are not commensurable with the other explicit modes -- which is what has been bounding every cross-method comparison. This mode uses the same f_ode as the symplectic schemes, Gauss-Radau, Bulirsch-Stoer and TDRK. integmode 25 is TDRK under tolerance control, using fortnum's new embedded pair. TDRK stays available fixed-step as 24: fixed step is what makes it comparable to the symplectic schemes at matched resolution, tolerance control is what makes it comparable to the adaptive ones at matched accuracy. Full orbit gains the same choice through fo_step_rkng_adaptive. Also fixes a state leak: explicit_h_carry is a warm start for the step controller and was never reset between orbits, so an orbit's first macro-step depended on which orbit the thread had run before it. Same particle, different evaluation counts depending on scheduling. orbit_sympl_init resets it. test_explicit_integmodes covers all five explicit modes against two oracles: the canonical Hamiltonian is an exact invariant, so its drift must respond to the tolerance (or to refinement, for the fixed-step mode); and the four adaptive methods -- collocation, extrapolation, an embedded explicit pair, a two-derivative pair -- must agree at tight tolerance, which four unrelated methods would not do if any shared a bug. Hard-wiring the tolerance in the TDRK driver fails both checks.
orbit_timestep_cashkarp45 took the step-size carry from solution%h(nsteps+1). ode_integrate trims solution%h to nsteps entries -- a step size belongs to the interval between two recorded points, not to a point -- so that read was one past the end, and whatever the memory held became the next macro-step's initial step size. ode_integrate_radau and ode_integrate_gbs allocate nsteps+1 for h as well as for t, so the identical expression is in bounds for integmodes 20 and 21; only the new Cash-Karp path was wrong. The effect was not a crash but a silently wrong cost curve: Cash-Karp's field-evaluation count came out nearly independent of the tolerance (1.5e7 at rtol 1e-6 against 1.8e7 at 1e-10, four decades of error for 18 per cent more work, which no fifth-order method can do). With the carry read correctly the same measurement gives 1.05e5 against 6.08e5 -- a 5.8x spread, and 20x cheaper at the loose end than the broken version reported.
The symplectic schemes test the radius of every Newton iterate, so a crossing of the last closed flux surface is caught inside the step, and advance_symplectic_with_boundary only bisects for the crossing if the raw step reports SYMPLECTIC_STEP_OUTSIDE_DOMAIN. Integmodes 20, 21, 22, 24 and 25 hand a whole substep to an adaptive integrator in one call and never looked between its endpoints, so they never reported it. A marker that crossed s = 1 was integrated onwards through extrapolated splines and recorded as confined. At paper scale (1000 alphas, 1 s trace) that put 44 markers at s > 1 for Bulirsch-Stoer and 42 for TDRK, one as far out as s = 2.8, all counted as confined: 76 losses against the 206 and 210 that RK4/5 and symplectic midpoint independently found. Every confined fraction from these modes was wrong, and so was my earlier report of a Bulirsch-Stoer instability at 1 s -- BS was not failing differently from the rest of the explicit family, it was failing the same way. Each driver now checks the radius after every substep and returns OUTSIDE_DOMAIN, which hands the step to the existing bisection. On a 1e-2 s trace over 1000 markers the five explicit modes now find 204-205 losses against symplectic midpoint's 205, with no confined marker above s = 0.992. Integrator failure also gets its own status. It was returning 1, which is the numeric value of SYMPLECTIC_STEP_OUTSIDE_DOMAIN, so a numerical failure was classified as a physical loss; it now returns SYMPLECTIC_STEP_MAXITER. The regression test asserts a physical invariant rather than a recorded result: a marker reported as confined is inside the last closed flux surface. One escaped marker fails it, which matters because a short trace produces few losses and a test needing many would be insensitive.
Correctness bug in this PR's integrators: losses were not detectedFound while trying to reproduce the JCP 2020 alpha-loss figure with these modes, which is exactly the measurement it breaks.
At paper scale — 1000 alphas from s = 0.75, traced to 1 s:
Markers sitting at s = 2.8, well outside the plasma, counted as confined. Two independent method families agree on ~207 losses; the explicit family found 76. This invalidates every confined fraction produced by integmodes 20/21/22/24/25 so far. It also retracts a finding I reported earlier as a method property: I described Bulirsch-Stoer as unstable at a 1 s trace. It is not — BS was failing in exactly the same way as the rest of the explicit family, and I mistook "the one I looked at first" for "the one that is broken". Short traces did not expose this. The work-precision study runs to 1e-3 s, where almost nothing is lost, and its final-state errors track the reference to 1e-5 for all of these modes. Accuracy was never the problem. FixEach explicit driver now checks the radius after every substep and returns Integrator failure also gets its own status code. It was returning Verification, 1000 markers on a 1e-2 s trace, same pinned initial conditions:
Regression test
Worth noting for the review: the existing tests all passed throughout. They check order of convergence and tolerance response on short traces, and none of them asserted that a marker called confined is inside the plasma. |
r is a radius, so r < 0 is not a position but the chart artefact of passing through the magnetic axis. Every other integrator in SIMPLE continues on the opposite ray, (r, theta) -> (|r|, theta + pi), which is the fix from #370. The explicit drivers did not, so a marker crossing the axis carried a negative radius into the next substep and the field was evaluated outside the spline domain. Integmode 21 records two r_negative events on a 64-marker near-axis case, so the path is live. Verified by comparing per-marker final states between builds identical except for the switch: one marker of 64 differs, by 4.8. The regression test cannot discriminate this and says so in its own comments -- two correct methods already differ by up to 5.6 on the same markers, because any trace long enough to reach the axis is long enough for the orbits to decorrelate. A sensitive test needs an axis-crossing orbit in the analytic field, where the trajectory is not chaotic; noted as follow-up. Also stops the boundary bisection from corrupting the step-size carry. locate_symplectic_boundary re-runs the step up to 64 times on shrinking fractions of dt, and each trial overwrote the carry, so it came back at bisection scale. No published number changes: every path that corrupts it ends the orbit immediately afterwards, and orbit_sympl_init resets it per particle. Fixed because 'currently unreachable' is a poor reason to leave it.
Every comparison against NaN is false, so NaN > sympl_rmax does not fire and a NaN orbit was waved through as a good step and reported as confined. The adaptive drivers usually fail first, since a NaN error norm never satisfies err <= 1 and the step is rejected until the size underflows, but the fixed-step TDRK has no error estimate to fail and nothing else on this path inspects the state. Second time today that NaN comparison semantics defeated a check that read correctly: the harness counted faulted orbits as survivors because awk compares a non-numeric field as a string, where "NaN" >= "0.999999" is true. Opposite directions, same cause -- a comparison written for numbers handed something that is not one.
There was a problem hiding this comment.
Review verdict: Approve
Review: PR #511 — Explicit high-order integmodes on the 4D canonical chart
Summary: This PR adds five new explicit integrators (Gauss-Radau/IAS15, Gragg-Bulirsch-Stoer, Cash-Karp RK45, and two-derivative RKN fixed/adaptive) on the 4D canonical Hamiltonian for benchmarking against the existing symplectic schemes. It also fixes a domain-exit detection bug (markers leaving the plasma without being caught), applies the axis chart switch (r < 0 → (|r|, θ+π), issue #370), and fixes a CMake build issue where the fortnum FetchContent pin was silently overridden by libneo. Tests cover energy conservation, cross-method agreement, and loss detection.
Findings:
-
[minor]
test/tests/test_explicit_loss_detection.py:40,200— The test referencesWOUT = SCRIPT_DIR.parent / "test_data" / "wout.nc"(i.e.,test/test_data/wout.nc), but that file is not tracked in the repo and is only downloadable viatest/test_data/get_test_data.sh(interactive prompt). In CI,wout.ncis not explicitly fetched (the CI step fetches GVEC QA data, notwout.nc). The test'sif not WOUT.exists(): sys.exit("missing simple.x or wout.nc")will cause it to exit non-zero in any environment wherewout.ncis absent, which may cause an unexpected CI failure. The path also usesSCRIPT_DIR.parentwhich resolves totest/(nottest/tests/), which appears correct given the file is attest/test_data/wout.nc, but the test file lives intest/tests/soSCRIPT_DIR.parentistest/— this is correct but worth confirming in CI logs. -
[minor]
src/orbit_symplectic_quasi.f90(g_ode, dF(1,j) formula) — The second-derivative formulag_odeis documented as verified with Mathematica (residual exactly zero, numeric spot-check at 3.6e-15). I cannot independently verify the algebra, but the testtest_explicit_integmodes.f90checks energy conservation and cross-method agreement, which would catch a derivation error. If the test passes in CI, the formula is likely correct. -
[minor]
src/orbit_symplectic_quasi.f90(sync_field_to_state) — Re-evaluating the field viaf_odeto fixz(5)is a documented workaround for a known limitation (the last internal stage's field side-effect). The approach is sound (one extra evaluation per macro-step), but it meansz(5)is recomputed from the final state rather than propagated through the step, which could differ slightly from the symplectic schemes in edge cases. The comment acknowledges this.
Verdict: Approve — The code is well-documented, the build fix is sound and addresses a real issue, the domain-exit detection addresses a real bug (markers at s > 1 reported as confined), and tests provide good coverage. The wout.nc dependency in the Python test is a minor infrastructure concern that should be verified in CI but does not block the code changes.
1 participant
Stacked on #510 (field-evaluation counters). Review that first.
What
The explicit high-order methods the integrator benchmark needs opposite the symplectic schemes:
Both run on the 4D canonical chart through the existing
f_ode, so they integrate the same canonical Hamiltonian, via the same field evaluations, as integmode 1-7 and 15. That is what makes a per-field-evaluation comparison like for like — unlikeintegmode = 0, which runs the 5D drift-kinetic form with a different RHS and different coordinates and so cannot be compared evaluation for evaluation.Requires the fortnum integrators in lazy-fortran/fortnum#59 and #60; the pin is bumped accordingly.
Three things needed to make this work
Step size carried across macro-steps (
explicit_h_carry). Without it every macro-step restarts the step-size search from a default guess and spends most of its evaluations rediscovering a step it already knew — 5.0M field evaluations against 4.0M on the same short run.A segfault that was not mine, and is worth knowing about. fortnum's
ode_solve_*helpers relied on allocation-on-assignment. libneo disables that globally with-fno-realloc-lhs(libneo/CMakeLists.txt:130) using a directory-scopedadd_compile_options, so it leaks into every subproject added after it — fortnum included. Under that flag,lhs = rhson an unallocated allocatable writes through a null descriptor and segfaults rather than allocating. It reproduced in a 20-line standalone program built with SIMPLE's flags but passed underfo, which does not set the flag. Fixed upstream in fortnum by allocating explicitly (this affected the pre-existingode_solve,ode_solve_dopandode_solve_vodetoo, not just the new code). Worth considering whether libneo should scope that flag to its own targets.A bridge between the quasi and symplectic paths. The explicit steppers live on the quasi path and read its module-level state, while the driver calls through the symplectic stepper signature.
orbit_timestep_sympl_radau15/_gbs16copy across. The quasi state is threadprivate so this is thread safe, and the copies are the price of reusingf_odeunchanged — which is what keeps the comparison honest.Known limitation, deliberately not addressed here
The macro-step granularity bounds how large a step the adaptive methods may take, so at the default
npoiper2they pay a high per-step cost for an order they cannot use. Gauss-Radau spends 4.0M evaluations against Bulirsch-Stoer's 0.84M on the same short run largely for this reason. The work-precision study has to sweep that parameter rather than hold it fixed — a fixed-npoiper2comparison would misrepresent both methods.Also worth noting for anyone benchmarking: with the default start mode, initial conditions are regenerated randomly per run, so confined fractions vary run to run (0.625 vs 0.875 on repeat runs of the same input). Benchmark runs need a pinned
start.dat.Verification