Add integmode 20 (Gauss-Radau) and 21 (Bulirsch-Stoer)

.gitignore

@@ -76,3 +76,4 @@ _skbuild/
 artifacts/plots/
 artifacts/notebooks/
 .venv/
+.sloptools/artifacts/

CMakeLists.txt

@@ -182,6 +182,32 @@ if(SIMPLE_DETERMINISTIC_FP)
     set(LIBNEO_DETERMINISTIC_FP ON CACHE BOOL
         "Use deterministic floating-point options in libneo" FORCE)
 endif()
+
+# fortnum provides the multidimensional root finder (multiroot_hybrids) that
+# replaced the bundled MINPACK hybrd1 in orbit_symplectic_quasi, and the ODE
+# integrators behind the explicit integmodes.
+#
+# This must be declared BEFORE libneo. libneo declares fortnum too, and
+# FetchContent honours the first declaration it sees and silently discards
+# later ones -- so with libneo first, the pin here has no effect and the build
+# gets whatever revision libneo happens to track. The failure is invisible in a
+# working-tree build, where find_or_fetch(libneo) resolves to a local checkout
+# that is usually already up to date, and appears only in a clean clone, as a
+# missing .mod file for whichever module SIMPLE needs and libneo's fortnum
+# lacks.
+#
+# The guard stays for the case where an enclosing project already provides
+# fortnum: there, that one wins deliberately.
+if(NOT TARGET fortnum)
+    include(FetchContent)
+    FetchContent_Declare(
+        fortnum
+        GIT_REPOSITORY https://github.com/lazy-fortran/fortnum.git
+        GIT_TAG c481d8e494d0c2ea8e9897bb77b39a1161b41c71
+    )
+    FetchContent_MakeAvailable(fortnum)
+endif()
+
 find_or_fetch(libneo)
 
 # Consume the scientific-I/O target selected by libneo. The variable fallback
@@ -239,18 +265,6 @@ else()
     message(STATUS "Fortplot disabled for NVHPC compiler (compatibility issues)")
 endif()
 
-# fortnum provides the multidimensional root finder (multiroot_hybrids) that
-# replaced the bundled MINPACK hybrd1 in orbit_symplectic_quasi.
-if(NOT TARGET fortnum)
-    include(FetchContent)
-    FetchContent_Declare(
-        fortnum
-        GIT_REPOSITORY https://github.com/lazy-fortran/fortnum.git
-        GIT_TAG 92de6e949a772cfffc73bb5295fe5e2b056b9c18
-    )
-    FetchContent_MakeAvailable(fortnum)
-endif()
-
 if (SIMPLE_TESTING)
     message(STATUS "Unit Tests enabled!")
 endif()

examples/simple_full.in

@@ -32,7 +32,7 @@ startmode = 1            !  mode for initial conditions:
 !                         5=distribute in volume ("global")
 grid_density = 0d0       ! for startmode 1 only, between 0.0 to 0.99, when 0.0 then no grid is made.
 special_ants_file = .False. ! if .True., a different start file is read (defined in samplers.f90), .False. uses standard filename (defined in samplers.f90)
-integmode = 1            ! mode for integrator: -1 = RK VMEC, 0 = RK, 1 = Euler1, 2 = Euler2, 3 = Midpoint, 4-7 = Gauss1-4, 15 = Lobatto3
+integmode = 1            ! integrator: -1 = RK VMEC, 0 = RK, 1 = Euler1, 2 = Euler2, 3 = Midpoint, 4-7 = Gauss1-4, 15 = Lobatto3, 20 = Gauss-Radau, 21 = Bulirsch-Stoer, 22 = Cash-Karp, 24 = TDRK (fixed step), 25 = TDRK (adaptive)
 relerr = 1d-13           ! tolerance for integrator. Set to 1d-13 for symplectic.
 tcut = -1d0              ! time when to do cut for classification, usually 1d-1, or -1 if no cuts desired
 debug = .False.          ! produce debugging output (.True./.False.). Use only in non-parallel mode!

src/orbit_fo_boris.f90

@@ -48,7 +48,9 @@ module orbit_fo_boris
   integer, parameter, public :: FO_OK = 0, FO_LOSS = 1, FO_LOCATE_FAIL = 2
 
   public :: fo_state_t, fo_init, fo_init_reference, fo_step, fo_energy, &
-    fo_mu, fo_to_gc, fo_to_reference_gc, accept_or_fail
+    fo_mu, fo_to_gc, fo_to_reference_gc, accept_or_fail, fo_step_rkng, &
+    fo_step_rkng_adaptive
+
 
   type :: fo_state_t
     real(dp) :: x(3)   = 0.0_dp    ! Cartesian position (scaled cm)
@@ -63,6 +65,16 @@ module orbit_fo_boris
     logical :: reference_field = .false.
   end type fo_state_t
 
+  ! State for the RKNG right-hand side. fortnum's rkng interface passes only
+  ! (t, y, yp, ypp) plus an unlimited-polymorphic ctx, while the force needs the
+  ! warm-start logical coordinate and the reference-field flag, and it updates
+  ! that warm start as it goes. Carrying them here mirrors how the guiding-centre
+  ! quasi module holds si and f, and keeps the interface fortnum expects.
+  type(fo_state_t) :: rkng_st
+  integer :: rkng_status = FO_OK
+  !$omp threadprivate(rkng_st, rkng_status)
+
+
 contains
 
   ! Cartesian (wedge) -> logical chart (rho, theta_B, phi_B). Warm damped Newton on
@@ -356,6 +368,127 @@ end subroutine field_at_logical
   ! fault). On fault Bvec etc. are undefined and the caller must not push. Loss is
   ! not decided here -- the field is defined through the clamped edge, and only the
   ! guiding-centre crossing rho>=1 in fo_to_gc is a confinement loss.
+  ! Lorentz acceleration in general Runge-Kutta-Nystrom form.
+  !
+  ! 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: it is the one
+  ! place a Nystrom-family method applies to SIMPLE directly, and so the most
+  ! faithful test of the original suggestion.
+  !
+  ! The warm-start logical coordinate is advanced on every force evaluation, as
+  ! in fo_step, so the inversion stays close to its previous solution. A failed
+  ! inversion is recorded in rkng_status and the acceleration is returned as
+  ! zero: fortnum has no way to abort a fixed-step integration mid-flight, so the
+  ! failure is reported after the step rather than propagated through it.
+  subroutine fo_rkng_force(t, y, yp, ypp, ctx)
+    real(dp), intent(in)  :: t
+    real(dp), intent(in)  :: y(:), yp(:)
+    real(dp), intent(out) :: ypp(:)
+    class(*), intent(in), optional :: ctx
+    real(dp) :: Bvec(3), Bmod, gradB(3), u(3), qcm
+    integer :: status
+
+    ypp = 0.0_dp
+    if (rkng_status /= FO_OK) return
+
+    call cart_field(y(1:3), rkng_st%u, Bvec, Bmod, gradB, u, status, &
+      rkng_st%reference_field)
+    if (status /= FO_OK) then
+      rkng_status = status
+      return
+    end if
+    rkng_st%u = u
+
+    qcm = rkng_st%charge/(c*rkng_st%ro0*rkng_st%mass)
+    ypp(1:3) = qcm*cross(yp(1:3), Bvec)
+  end subroutine fo_rkng_force
+
+  ! One full-orbit step with a general Runge-Kutta-Nystrom method, as an
+  ! alternative to the Boris pusher over the same dt and the same field path.
+  !
+  ! Boris is second order and volume preserving, and conserves energy to
+  ! round-off in a static magnetic field because its rotation is exact. RKNG is
+  ! higher order but not structure preserving, so the comparison to make is
+  ! energy drift at matched cost, not energy drift alone.
+  subroutine fo_step_rkng(st, status, nsub)
+    use fortnum_ode_tdrk, only: rkng_integrate_fixed
+    use fortnum_status, only: fortnum_status_t, FORTNUM_OK
+    type(fo_state_t), intent(inout) :: st
+    integer, intent(out) :: status
+    integer, intent(in), optional :: nsub
+    real(dp) :: yend(3), ypend(3)
+    integer :: nfev, nsteps
+    type(fortnum_status_t) :: fstat
+
+    nsteps = 1
+    if (present(nsub)) nsteps = max(1, nsub)
+
+    rkng_st = st
+    rkng_status = FO_OK
+
+    call rkng_integrate_fixed(fo_rkng_force, 0.0_dp, st%dt, st%x, st%v, &
+      nsteps, yend, ypend, nfev, fstat)
+
+    if (rkng_status /= FO_OK) then
+      status = rkng_status        ! leave st at the last resolved state
+      return
+    end if
+    if (fstat%code /= FORTNUM_OK) then
+      status = FO_LOCATE_FAIL
+      return
+    end if
+
+    st%x = yend
+    st%v = ypend
+    st%u = rkng_st%u
+    status = FO_OK
+  end subroutine fo_step_rkng
+
+  ! The same RKNG method under tolerance control rather than a fixed substep
+  ! count. rtol/atol replace nsub as the accuracy knob, which is what lets the
+  ! full-orbit arena be swept the same way the guiding-centre one is.
+  !
+  ! h_carry passes the accepted step size from one macro-step to the next: a
+  ! gyro-orbit's timescale barely changes between steps, so restarting the
+  ! controller from scratch every dt would waste the rejections it already paid
+  ! for. Pass 0 on the first call and feed the returned value back afterwards.
+  subroutine fo_step_rkng_adaptive(st, rtol, atol, h_carry, status, nfev)
+    use fortnum_ode_tdrk, only: rkng_integrate_adaptive
+    use fortnum_status, only: fortnum_status_t, FORTNUM_OK
+    type(fo_state_t), intent(inout) :: st
+    real(dp), intent(in) :: rtol, atol
+    real(dp), intent(inout) :: h_carry
+    integer, intent(out) :: status
+    integer, intent(out), optional :: nfev
+    real(dp) :: yend(3), ypend(3), hlast
+    integer :: nf, naccept, nreject
+    type(fortnum_status_t) :: fstat
+
+    rkng_st = st
+    rkng_status = FO_OK
+
+    call rkng_integrate_adaptive(fo_rkng_force, 0.0_dp, st%dt, st%x, st%v, &
+      rtol, atol, h_carry, yend, ypend, hlast, nf, naccept, nreject, fstat)
+
+    if (present(nfev)) nfev = nf
+
+    if (rkng_status /= FO_OK) then
+      status = rkng_status        ! leave st at the last resolved state
+      return
+    end if
+    if (fstat%code /= FORTNUM_OK) then
+      status = FO_LOCATE_FAIL
+      return
+    end if
+
+    st%x = yend
+    st%v = ypend
+    st%u = rkng_st%u
+    h_carry = hlast
+    status = FO_OK
+  end subroutine fo_step_rkng_adaptive
+
   subroutine cart_field(x, u_guess, Bvec, Bmod, gradB, u_out, status, &
       reference_field)
     real(dp), intent(in) :: x(3), u_guess(3)

src/orbit_symplectic.f90

@@ -6,7 +6,7 @@ module orbit_symplectic
   eval_field => evaluate
 use orbit_symplectic_base, only: symplectic_integrator_t, multistage_integrator_t, &
   RK45, EXPL_IMPL_EULER, IMPL_EXPL_EULER, MIDPOINT, GAUSS1, GAUSS2, GAUSS3, GAUSS4, &
-  LOBATTO3, S_MAX, orbit_timestep_sympl_i, extrap_field, sympl_rmax, &
+  LOBATTO3, RADAU15, GBS16, CASHKARP45, TDRK24, TDRK24A, S_MAX, orbit_timestep_sympl_i, extrap_field, sympl_rmax, &
   coeff_rk_gauss, coeff_rk_lobatto, f_rk_lobatto, &
   SYMPLECTIC_STEP_OK, SYMPLECTIC_STEP_OUTSIDE_DOMAIN, &
   SYMPLECTIC_STEP_MAXITER, SYMPLECTIC_STEP_LINEAR_SOLVE, &
@@ -15,6 +15,10 @@ module orbit_symplectic
   boundary_event_radial_tolerance, symplectic_newton_warning_mode
 use orbit_symplectic_quasi, only: orbit_timestep_quasi, timestep_expl_impl_euler_quasi, &
   timestep_impl_expl_euler_quasi, timestep_midpoint_quasi, orbit_timestep_rk45, &
+  orbit_timestep_radau15, orbit_timestep_gbs16, orbit_timestep_tdrk24, &
+  orbit_timestep_tdrk24a, orbit_timestep_cashkarp45, &
+  reset_explicit_step_carry, &
+  si_quasi => si, f_quasi => f, &
   timestep_rk_gauss_quasi, timestep_rk_lobatto_quasi
 use orbit_symplectic_euler1, only: sympl_euler1_residual, sympl_euler1_jacobian, &
   sympl_euler1_newton_iter, sympl_euler1_extrapolate_field, &
@@ -231,6 +235,7 @@ recursive subroutine orbit_sympl_init(si, f, z, dt, ntau, rtol_init, mode_init)
 
   si%atol = 1d-15
   si%rtol = rtol_init
+  call reset_explicit_step_carry
 
   si%ntau = ntau
   si%dt = dt
@@ -273,6 +278,21 @@ recursive subroutine orbit_sympl_init(si, f, z, dt, ntau, rtol_init, mode_init)
     case (LOBATTO3)
       raw_timestep_sympl => orbit_timestep_sympl_lobatto3
       orbit_timestep_quasi => orbit_timestep_quasi_lobatto3
+    case (RADAU15)
+      raw_timestep_sympl => orbit_timestep_sympl_radau15
+      orbit_timestep_quasi => orbit_timestep_radau15
+    case (GBS16)
+      raw_timestep_sympl => orbit_timestep_sympl_gbs16
+      orbit_timestep_quasi => orbit_timestep_gbs16
+    case (TDRK24)
+      raw_timestep_sympl => orbit_timestep_sympl_tdrk24
+      orbit_timestep_quasi => orbit_timestep_tdrk24
+    case (TDRK24A)
+      raw_timestep_sympl => orbit_timestep_sympl_tdrk24a
+      orbit_timestep_quasi => orbit_timestep_tdrk24a
+    case (CASHKARP45)
+      raw_timestep_sympl => orbit_timestep_sympl_cashkarp45
+      orbit_timestep_quasi => orbit_timestep_cashkarp45
     case default
       print *, 'invalid mode for orbit_timestep_sympl: ', mode_init
       error stop
@@ -284,6 +304,75 @@ recursive subroutine orbit_sympl_init(si, f, z, dt, ntau, rtol_init, mode_init)
   end if
 end subroutine orbit_sympl_init
 
+  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+  !
+! Bridges the explicit high-order steppers, which live on the quasi path and
+! read the module-level state there, into the symplectic stepper signature the
+! driver calls. 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: these methods integrate exactly the same right-hand side
+! as the symplectic schemes, through the same field evaluations.
+recursive subroutine orbit_timestep_sympl_radau15(si, f, ierr)
+  type(symplectic_integrator_t), intent(inout) :: si
+  type(field_can_t), intent(inout) :: f
+  integer, intent(out) :: ierr
+
+  si_quasi = si
+  f_quasi = f
+  call orbit_timestep_radau15(ierr)
+  si = si_quasi
+  f = f_quasi
+end subroutine orbit_timestep_sympl_radau15
+
+recursive subroutine orbit_timestep_sympl_gbs16(si, f, ierr)
+  type(symplectic_integrator_t), intent(inout) :: si
+  type(field_can_t), intent(inout) :: f
+  integer, intent(out) :: ierr
+
+  si_quasi = si
+  f_quasi = f
+  call orbit_timestep_gbs16(ierr)
+  si = si_quasi
+  f = f_quasi
+end subroutine orbit_timestep_sympl_gbs16
+
+recursive subroutine orbit_timestep_sympl_tdrk24(si, f, ierr)
+  type(symplectic_integrator_t), intent(inout) :: si
+  type(field_can_t), intent(inout) :: f
+  integer, intent(out) :: ierr
+
+  si_quasi = si
+  f_quasi = f
+  call orbit_timestep_tdrk24(ierr)
+  si = si_quasi
+  f = f_quasi
+end subroutine orbit_timestep_sympl_tdrk24
+
+recursive subroutine orbit_timestep_sympl_tdrk24a(si, f, ierr)
+  type(symplectic_integrator_t), intent(inout) :: si
+  type(field_can_t), intent(inout) :: f
+  integer, intent(out) :: ierr
+
+  si_quasi = si
+  f_quasi = f
+  call orbit_timestep_tdrk24a(ierr)
+  si = si_quasi
+  f = f_quasi
+end subroutine orbit_timestep_sympl_tdrk24a
+
+recursive subroutine orbit_timestep_sympl_cashkarp45(si, f, ierr)
+  type(symplectic_integrator_t), intent(inout) :: si
+  type(field_can_t), intent(inout) :: f
+  integer, intent(out) :: ierr
+
+  si_quasi = si
+  f_quasi = f
+  call orbit_timestep_cashkarp45(ierr)
+  si = si_quasi
+  f = f_quasi
+end subroutine orbit_timestep_sympl_cashkarp45
+
 
   !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
   !

src/orbit_symplectic_base.f90

@@ -19,6 +19,31 @@ module orbit_symplectic_base
     ! Integration methods
     integer, parameter :: RK45 = 0, EXPL_IMPL_EULER = 1, IMPL_EXPL_EULER = 2, &
              MIDPOINT = 3, GAUSS1 = 4, GAUSS2 = 5, GAUSS3 = 6, GAUSS4 = 7, LOBATTO3 = 15
+
+    ! High-order explicit methods on the 4D canonical chart, for benchmarking
+    ! against the symplectic schemes above. They integrate the same canonical
+    ! Hamiltonian through the same f_ode right-hand side, so a comparison per
+    ! field evaluation is like for like -- unlike integmode = 0, which runs the
+    ! 5D drift-kinetic form with a different right-hand side and coordinates.
+    !
+    ! RADAU15 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. GBS16 is Gragg-Bulirsch-Stoer
+    ! extrapolation, the classical high-accuracy explicit method of celestial
+    ! mechanics.
+    ! TDRK24 is the Runge-Kutta-Nystrom transplant: an RKN tableau applied to
+    ! the first-order guiding-centre system through zdd = F'(z)F(z), which is a
+    ! special two-derivative Runge-Kutta method (Chan & Tsai 2010). Fixed step,
+    ! so npoiper2 sets its resolution as it does for the symplectic schemes;
+    ! TDRK24A is the same method under tolerance control, which is what makes it
+    ! comparable to the adaptive methods rather than to the symplectic ones.
+    !
+    ! CASHKARP45 is the classical adaptive workhorse on this same chart.
+    ! integmode = 0 also runs a Cash-Karp, but on the 5D drift-kinetic form with
+    ! a different right-hand side and coordinates, so its accuracy and cost are
+    ! not commensurable with anything here.
+    integer, parameter :: RADAU15 = 20, GBS16 = 21, CASHKARP45 = 22, &
+             TDRK24 = 24, TDRK24A = 25
     integer, parameter :: SYMPLECTIC_STEP_OK = 0
     integer, parameter :: SYMPLECTIC_STEP_OUTSIDE_DOMAIN = 1
     integer, parameter :: SYMPLECTIC_STEP_MAXITER = 2