Locate the bounce event by chunked integration with the DVODE root convention - #164
Draft
krystophny wants to merge 7 commits into
Draft
Locate the bounce event by chunked integration with the DVODE root convention#164krystophny wants to merge 7 commits into
krystophny wants to merge 7 commits into
Conversation
Wire the fortnum numerical-core library via CMake FetchContent (git@github.com:lazy-fortran/fortnum.git, pinned to a7faa3c) guarded by if(NOT TARGET fortnum) so a parent build that already provides the target is reused. Link the static library into neo_rt; its module directory is exported PUBLIC, so neo_rt and the diagnostics library that links neo_rt both see the fortnum modules. This commit adds the dependency without changing any solver call site; the existing DVODE path is untouched and still in use. Subsequent commits move the orbit integrations onto the fortnum ODE solver and then drop DVODE.
Replace the DVODE calls in the two fixed-interval, no-event integrations with the fortnum DOP853 solver: - orbit.f90 bounce_fast: ode_solve_dop from 0 to taub, take the final recorded state. The fortnum status maps onto the legacy istate the transport caller branches on (2 success, -1 step budget exhausted, 0 other), so dvode_error_context still fires on the exhausted-budget case and the public bounce_fast signature is unchanged. - diag/diag_bounce_debug.f90 probe_bounce: ode_integrate_dop, then read solution%nsteps, solution%nrejected and the last solution%h into the diagnostic columns the scan previously filled from get_stats. - freq.f90: drop the vode_thread_init shim. The fortnum integrator keeps no module-level state, so freq_thread_init has nothing VODE-specific to initialise. The event-based bounce_integral and the vode dependency itself stay until the next commit so this change is reviewable on its own.
Two defects surfaced once bounce_fast ran on fortnum DOP853 instead of DVODE, both caught only by the CI Debug build (gfortran 13, -ffpe-trap=zero,overflow,invalid,underflow): - timestep_transport left ydot(7), the unused abs(B) slot, uninitialised. DVODE ignored the stray derivative, but DOP853 carries every component through its stage combinations, so the leftover denormal tripped the underflow trap inside dop853_step. Zero the trailing integrands the routine does not compute, matching what timestep already does. - The bounceavg(3:4) perturbed-Hamiltonian integrands oscillate at the resonant harmonic but are driven, not fed back into the orbit, so the rtol error estimate barely constrains them. The unconstrained adaptive step grew too coarse and the oscillatory integral came out wrong (D11 off by an order of magnitude). Drive bounce_fast through ode_integrate_dop with hmax = taub/200 to keep the step fine enough; this recovers the DVODE result to four significant figures.
bounce_integral switches from the chunked DVODE root search to the fortnum dop853 integrator plus a post-integration event scan. ode_problem_t carries the bounce event: event => bounce_event, event_direction set per branch (passing crosses a full 2*pi turn in the sign_vpar_htheta direction -> RISING for sign_vpar_htheta>0, FALLING otherwise; trapped returns to th0 rising from below -> RISING, matching the old (yold(1)-th0)<0 accept filter), terminal_event=.true. After ode_integrate_dop, ode_event_scan locates the crossing and yields t_event/y_event; the no-event branch keeps the existing diagnostic dump. The event functions reproduce the original bounceroots zero locations: passing g = (theta-th0) - sign*2*pi, trapped g = theta-th0. CMakeLists drops find_or_fetch(vode) and the vode include/link/dependency. neo_rt no longer pulls the external DVODE solver; fortnum is the sole ODE backend. .gitignore ignores .sloptools/artifacts/.
…nvention The fortnum migration integrated a single fixed [0, n_turns*dt] span and scanned the recorded trace once. Two defects made 9/14 ripple golden cases abort with no bounce event for near-separatrix passing orbits: - The passing event was g = (theta-th0) - sign_vpar_htheta*2*pi, which roots at theta-th0 = sign_vpar_htheta*2*pi. When sign_vpar_htheta < 0 the orbit still advances theta by +2*pi, so g never crossed zero and the scan found nothing. Restore the old bounceroots form sign_vpar_htheta*(2*pi-(theta-th0)) (root always at +2*pi advance) and the matching trapped form, and select the crossing direction from the sign of theta motion through each root. - ode_integrate_dop has no terminal-event support, so a single fixed span both failed to bracket late events and integrated ~200 turns of a near-separatrix orbit, exhausting max_steps. Integrate in dt-sized chunks and scan after each, re-seeding from the chunk endpoint until the crossing is bracketed, mirroring the old DVODE until-event loop. Bump fortnum to 4c560b0.
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.
This draft PR publishes the existing committed state of
ode-verify-eventsso the local branch/worktree can be retired after review.Checks were not run in this cleanup operation; the branch-specific CI should provide validation.