diff --git a/.gitignore b/.gitignore index 5ef4edaa..499def09 100644 --- a/.gitignore +++ b/.gitignore @@ -76,3 +76,4 @@ _skbuild/ artifacts/plots/ artifacts/notebooks/ .venv/ +.sloptools/artifacts/ diff --git a/CMakeLists.txt b/CMakeLists.txt index 5aa3b357..6a671d0c 100644 --- a/CMakeLists.txt +++ b/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() diff --git a/examples/simple_full.in b/examples/simple_full.in index 38a69bc0..599d710a 100644 --- a/examples/simple_full.in +++ b/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! diff --git a/src/orbit_fo_boris.f90 b/src/orbit_fo_boris.f90 index 807d24c7..eeb26f95 100644 --- a/src/orbit_fo_boris.f90 +++ b/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) diff --git a/src/orbit_symplectic.f90 b/src/orbit_symplectic.f90 index 39acacc2..5015c877 100644 --- a/src/orbit_symplectic.f90 +++ b/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 ! diff --git a/src/orbit_symplectic_base.f90 b/src/orbit_symplectic_base.f90 index 49747778..40983146 100644 --- a/src/orbit_symplectic_base.f90 +++ b/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 diff --git a/src/orbit_symplectic_quasi.f90 b/src/orbit_symplectic_quasi.f90 index 601deca1..4552ccf8 100644 --- a/src/orbit_symplectic_quasi.f90 +++ b/src/orbit_symplectic_quasi.f90 @@ -1,9 +1,11 @@ module orbit_symplectic_quasi use util, only: pi -use field_can_mod, only: eval_field => evaluate, field_can_t, get_derivatives +use field_can_mod, only: eval_field => evaluate, field_can_t, get_derivatives, & + get_derivatives2 use orbit_symplectic_base, only: symplectic_integrator_t, multistage_integrator_t, & - orbit_timestep_quasi_i, coeff_rk_gauss, coeff_rk_lobatto, f_rk_lobatto, sympl_rmax + orbit_timestep_quasi_i, coeff_rk_gauss, coeff_rk_lobatto, f_rk_lobatto, sympl_rmax, & + SYMPLECTIC_STEP_OK, SYMPLECTIC_STEP_OUTSIDE_DOMAIN, SYMPLECTIC_STEP_MAXITER use fortnum_multiroot, only: multiroot_hybrids use fortnum_status, only: fortnum_status_t use diag_counters, only: count_event, EVT_R_NEGATIVE @@ -24,8 +26,84 @@ module orbit_symplectic_quasi procedure(orbit_timestep_quasi_i), pointer :: orbit_timestep_quasi => null() +! Step size carried across macro-steps by the adaptive explicit integrators. +! Without this each macro-step restarts the step-size search from a default +! guess and spends most of its evaluations re-discovering a step it already +! knew, which inflates the evaluation count by more than an order of magnitude +! and would make any cost comparison meaningless. +! +! The carry belongs to one orbit, not to the thread: leaving it set when the +! next orbit starts makes that orbit's first macro-step depend on which orbit +! the thread happened to run before it, so the same particle gives different +! evaluation counts (and, through the controller's history, slightly different +! trajectories) depending on scheduling. orbit_sympl_init resets it. +real(dp) :: explicit_h_carry = 0d0 + !$omp threadprivate(explicit_h_carry) + contains +subroutine reset_explicit_step_carry + explicit_h_carry = 0d0 +end subroutine reset_explicit_step_carry + + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! +! Restore the chart invariants the rest of SIMPLE maintains inside a step, and +! report the status the caller's boundary bisection expects. +! +! The symplectic schemes reach both of these through their Newton solve, which +! inspects every iterate. The explicit drivers hand a whole substep to an +! adaptive integrator in one call and see nothing between its endpoints, so +! both have to be re-imposed here, at the same dtaumin granularity the +! symplectic schemes work at. +! +! Leaving the plasma. A marker that crossed the last closed flux surface used +! to be integrated onwards through extrapolated splines and reported as still +! confined: at a 1 s trace that put 44 of 1000 markers at s > 1, one as far out +! as s = 2.8, which is why the explicit modes found 76 losses where the +! symplectic schemes and RK4/5 both independently found about 207. Returning +! OUTSIDE_DOMAIN hands the step to locate_symplectic_boundary, which bisects on +! si%dt exactly as it does for the symplectic schemes. +! +! Crossing the axis. r is a radius, so r < 0 is not a position but a chart +! artefact of passing through the magnetic axis; every other integrator in +! SIMPLE continues on the opposite ray, (r, theta) -> (|r|, theta + pi), rather +! than the pre-#370 teleport to r = 0.01. Without it the next substep evaluates +! the field at a negative radius, outside the spline domain entirely. +! +! h_entry is the step-size carry as it stood before the substep. A step that +! does not complete must not leave its trial step size behind: the boundary +! bisection re-runs the step up to 64 times on ever smaller fractions of dt, +! and without this the carry would come back at bisection scale and throttle +! every later step of the orbit. +subroutine explicit_step_finish(h_entry, ierr) + real(dp), intent(in) :: h_entry + integer, intent(out) :: ierr + + ! A non-finite state has to be caught before the comparisons below, because + ! every comparison against NaN is false: NaN > sympl_rmax does not fire, so a + ! NaN orbit would be waved through as a good step and reported as confined. + ! The adaptive drivers usually fail first -- a NaN error norm never satisfies + ! err <= 1, so 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. + if (any(si%z(1:4) /= si%z(1:4))) then + ierr = SYMPLECTIC_STEP_MAXITER + explicit_h_carry = h_entry + return + end if + + if (si%z(1) < 0d0) then + call count_event(EVT_R_NEGATIVE) + si%z(1) = -si%z(1) + si%z(2) = si%z(2) + pi + end if + + ierr = SYMPLECTIC_STEP_OK + if (si%z(1) > sympl_rmax) ierr = SYMPLECTIC_STEP_OUTSIDE_DOMAIN + if (ierr /= SYMPLECTIC_STEP_OK) explicit_h_carry = h_entry +end subroutine explicit_step_finish + ! ! Wrapper routines for ODEPACK ! @@ -534,6 +612,119 @@ subroutine f_ode(tau, z, zdot) end subroutine f_ode + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! +! Index into the packed symmetric second-derivative arrays d2H, d2pth, d2vpar. +! +! field_can_base documents the order as +! 1:(r,r) 2:(r,th) 3:(r,ph) 4:(th,th) 5:(th,ph) 6:(ph,ph) +! 7:(pph,r) 8:(pph,th) 9:(pph,ph) 10:(pph,pph) +! over the variables (r, th, ph, pph) = (1, 2, 3, 4). +pure function d2idx(i, j) result(k) + integer, intent(in) :: i, j + integer :: k + integer :: lo, hi + + lo = min(i, j) + hi = max(i, j) + + if (hi == 4) then + k = 6 + lo ! (r,pph)=7, (th,pph)=8, (ph,pph)=9, (pph,pph)=10 + else if (lo == 1) then + k = hi ! (r,r)=1, (r,th)=2, (r,ph)=3 + else if (lo == 2) then + k = hi + 2 ! (th,th)=4, (th,ph)=5 + else + k = 6 ! (ph,ph) + end if +end function d2idx + + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! +! Second time derivative of the canonical guiding-centre state, +! +! G(z) = F'(z) F(z), +! +! which is what a two-derivative Runge-Kutta method integrates in place of F. +! +! The derivation and its verification live in the study repo at +! derivation/tdrk_guiding_center_G.wl: Mathematica differentiates f_ode's F +! symbolically and compares against this closed form. The symbolic residual is +! exactly zero and a numeric spot-check with nonlinear, fully coupled test +! functions agrees to 3.6e-15. +! +! The cost argument rests on which derivatives appear. That check is in the same +! file and reports second derivatives of H and pth, first of vpar, hth and hph, +! and no third derivatives anywhere -- so G needs exactly one mode_secders = 2 +! evaluation and no extra F evaluations. Note also that only d2H indices 1 to 9 +! are touched: F depends on H solely through dH(1:3), so differentiating it can +! never raise index 10, which get_derivatives2 does not fill. +! +! hth and hph are geometric and carry no pph dependence, hence the explicit zero +! in the fourth slot of their gradients rather than an out-of-bounds read of a +! three-element array. +subroutine g_ode(tau, z, g) + real(dp), intent(in) :: tau + real(dp), intent(in) :: z(:) + real(dp), intent(out) :: g(:) + + real(dp) :: Hprime, dHprime(4), dhth4(4), dhph4(4) + real(dp) :: zdot(4), dF(4, 4) + integer :: j + + call eval_field(f, z(1), z(2), z(3), 2) + call get_derivatives2(f, z(4)) + + Hprime = f%dH(1)/f%dpth(1) + + dhth4(1:3) = f%dhth + dhth4(4) = 0d0 + dhph4(1:3) = f%dhph + dhph4(4) = 0d0 + + zdot(1) = -(f%dH(2) - f%hth/f%hph*f%dH(3))/f%dpth(1) + zdot(2) = Hprime + zdot(3) = (f%vpar - Hprime*f%hth)/f%hph + zdot(4) = -(f%dH(3) - Hprime*f%dpth(3)) + + ! Hprime is a quotient of two state-dependent quantities, so both numerator + ! and denominator contribute. This is the term a hand derivation usually drops. + do j = 1, 4 + dHprime(j) = (f%d2H(d2idx(1, j))*f%dpth(1) & + - f%dH(1)*f%d2pth(d2idx(1, j)))/f%dpth(1)**2 + end do + + do j = 1, 4 + dF(1, j) = -((f%d2H(d2idx(2, j)) & + - (dhth4(j)/f%hph - f%hth*dhph4(j)/f%hph**2)*f%dH(3) & + - f%hth/f%hph*f%d2H(d2idx(3, j)))*f%dpth(1) & + - (f%dH(2) - f%hth/f%hph*f%dH(3))*f%d2pth(d2idx(1, j))) & + /f%dpth(1)**2 + + dF(2, j) = dHprime(j) + + dF(3, j) = ((f%dvpar(j) - dHprime(j)*f%hth - Hprime*dhth4(j))*f%hph & + - (f%vpar - Hprime*f%hth)*dhph4(j))/f%hph**2 + + dF(4, j) = -(f%d2H(d2idx(3, j)) - dHprime(j)*f%dpth(3) & + - Hprime*f%d2pth(d2idx(3, j))) + end do + + g(1:4) = matmul(dF, zdot) +end subroutine g_ode + + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! +! Adapter from g_ode to fortnum's ode_rhs2_t. +subroutine g_ode_fortnum(t, y, g, ctx) + real(dp), intent(in) :: t + real(dp), intent(in) :: y(:) + real(dp), intent(out) :: g(:) + class(*), intent(in), optional :: ctx + + call g_ode(t, y, g) +end subroutine g_ode_fortnum + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! subroutine orbit_timestep_rk45(ierr) @@ -550,4 +741,276 @@ subroutine orbit_timestep_rk45(ierr) end do end subroutine orbit_timestep_rk45 + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! +! Adapter from f_ode to fortnum's ode_rhs_t, which carries an optional context +! argument that this right-hand side does not need. +subroutine f_ode_fortnum(t, y, dydt, ctx) + real(dp), intent(in) :: t + real(dp), intent(in) :: y(:) + real(dp), intent(out) :: dydt(:) + class(*), intent(in), optional :: ctx + + call f_ode(t, y, dydt) +end subroutine f_ode_fortnum + + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! +! Gauss-Radau (IAS15-style, 15th order) on the same canonical Hamiltonian the +! symplectic schemes use. +subroutine orbit_timestep_radau15(ierr) + ! + use fortnum_ode_gauss_radau, only : ode_integrate_radau + use fortnum_ode, only : ode_problem_t, ode_solution_t + use fortnum_status, only : fortnum_status_t, FORTNUM_OK + integer, intent(out) :: ierr + integer :: ktau, nlast + type(ode_problem_t) :: problem + type(ode_solution_t) :: solution + type(fortnum_status_t) :: status + real(dp) :: h_entry + + ierr = 0 + ktau = 0 + problem%rhs => f_ode_fortnum + problem%rtol = si%rtol + problem%atol = si%atol + allocate(problem%y0(4)) + do while(ktau .lt. si%ntau) + h_entry = explicit_h_carry + problem%t0 = ktau*si%dt + problem%t1 = (ktau+1)*si%dt + problem%y0 = si%z(1:4) + problem%h0 = explicit_h_carry + call ode_integrate_radau(problem, solution, status) + if (status%code /= FORTNUM_OK) then + ierr = SYMPLECTIC_STEP_MAXITER + explicit_h_carry = h_entry + return + end if + nlast = size(solution%t) + si%z(1:4) = solution%y(:, nlast) + if (nlast > 1) explicit_h_carry = solution%h(nlast) + call explicit_step_finish(h_entry, ierr) + if (ierr /= SYMPLECTIC_STEP_OK) return + ktau = ktau+1 + end do + call sync_field_to_state +end subroutine orbit_timestep_radau15 + + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! +! Re-evaluate the field at the converged end-of-step state. +! +! 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) +! carries an error the tolerance cannot remove: measured on the sweep, +! components 1-4 of the final state converged to 5e-9 while z(5) stalled at +! 3.3e-4 and got no better from rtol 1e-8 to 1e-10. +! +! One evaluation per macro-step, against the hundreds a step already costs. +subroutine sync_field_to_state + real(dp) :: zdot_discard(4) + + call f_ode(0d0, si%z(1:4), zdot_discard) +end subroutine sync_field_to_state + + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! +! Gragg-Bulirsch-Stoer extrapolation on the same canonical Hamiltonian. +subroutine orbit_timestep_gbs16(ierr) + ! + use fortnum_ode_extrapolation, only : ode_integrate_gbs + use fortnum_ode, only : ode_problem_t, ode_solution_t + use fortnum_status, only : fortnum_status_t, FORTNUM_OK + integer, intent(out) :: ierr + integer :: ktau, nlast + type(ode_problem_t) :: problem + type(ode_solution_t) :: solution + type(fortnum_status_t) :: status + real(dp) :: h_entry + + ierr = 0 + ktau = 0 + problem%rhs => f_ode_fortnum + problem%rtol = si%rtol + problem%atol = si%atol + allocate(problem%y0(4)) + do while(ktau .lt. si%ntau) + h_entry = explicit_h_carry + problem%t0 = ktau*si%dt + problem%t1 = (ktau+1)*si%dt + problem%y0 = si%z(1:4) + problem%h0 = explicit_h_carry + call ode_integrate_gbs(problem, solution, status) + if (status%code /= FORTNUM_OK) then + ierr = SYMPLECTIC_STEP_MAXITER + explicit_h_carry = h_entry + return + end if + nlast = size(solution%t) + si%z(1:4) = solution%y(:, nlast) + if (nlast > 1) explicit_h_carry = solution%h(nlast) + call explicit_step_finish(h_entry, ierr) + if (ierr /= SYMPLECTIC_STEP_OK) return + ktau = ktau+1 + end do + call sync_field_to_state +end subroutine orbit_timestep_gbs16 + + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! +! Two-derivative Runge-Kutta on the canonical Hamiltonian: the Runge-Kutta- +! Nystrom transplant. +! +! 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), which is in special second-order form because zd is itself determined by +! z, and re-synchronising the velocity as F(z) at the start of every step turns +! the scheme into 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 used directly here. +! +! Because the velocity is recomputed rather than propagated, only the position +! order conditions are needed, so order 4 costs two stages instead of three. +! +! Fixed step, like the symplectic schemes: fortnum's TDRK has no embedded error +! estimate yet, so the accuracy parameter is the step count rather than a +! tolerance. npoiper2 therefore sets resolution here exactly as it does for the +! symplectic integrators, which also makes the two directly comparable. +subroutine orbit_timestep_tdrk24(ierr) + ! + use fortnum_ode_tdrk, only : tdrk_integrate_fixed + use fortnum_status, only : fortnum_status_t, FORTNUM_OK + integer, intent(out) :: ierr + integer :: ktau, nfev_f, nfev_g + real(dp) :: yend(4) + type(fortnum_status_t) :: status + real(dp) :: h_entry + + ierr = 0 + ktau = 0 + do while(ktau .lt. si%ntau) + h_entry = explicit_h_carry + call tdrk_integrate_fixed(f_ode_fortnum, g_ode_fortnum, & + ktau*si%dt, (ktau+1)*si%dt, si%z(1:4), 1, 4, .false., yend, & + nfev_f, nfev_g, status) + if (status%code /= FORTNUM_OK) then + ierr = SYMPLECTIC_STEP_MAXITER + explicit_h_carry = h_entry + return + end if + si%z(1:4) = yend + call explicit_step_finish(h_entry, ierr) + if (ierr /= SYMPLECTIC_STEP_OK) return + ktau = ktau+1 + end do + call sync_field_to_state +end subroutine orbit_timestep_tdrk24 + + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! +! Two-derivative Runge-Kutta with error control: the same method as +! orbit_timestep_tdrk24, driven by a tolerance instead of a step count. +! +! The embedded estimate needs the three-stage order-4 tableau. The two-stage +! one, which orbit_timestep_tdrk24 uses, cannot carry one: with c = (0, 1/2) +! the order-3 conditions already force the order-4 weights, so the two +! solutions coincide and their difference is identically zero. One extra G +! evaluation per step buys the estimate. +! +! Both modes are kept. Fixed step is what makes TDRK comparable to the +! symplectic schemes at matched resolution; tolerance control is what makes it +! comparable to Gauss-Radau, Bulirsch-Stoer and Cash-Karp at matched accuracy. +subroutine orbit_timestep_tdrk24a(ierr) + ! + use fortnum_ode_tdrk, only : tdrk_integrate_adaptive + use fortnum_status, only : fortnum_status_t, FORTNUM_OK + integer, intent(out) :: ierr + integer :: ktau, nfev_f, nfev_g, naccept, nreject + real(dp) :: yend(4), hlast + type(fortnum_status_t) :: status + real(dp) :: h_entry + + ierr = 0 + ktau = 0 + do while(ktau .lt. si%ntau) + h_entry = explicit_h_carry + call tdrk_integrate_adaptive(f_ode_fortnum, g_ode_fortnum, & + ktau*si%dt, (ktau+1)*si%dt, si%z(1:4), si%rtol, si%atol, & + explicit_h_carry, yend, hlast, nfev_f, nfev_g, naccept, nreject, status) + if (status%code /= FORTNUM_OK) then + ierr = SYMPLECTIC_STEP_MAXITER + explicit_h_carry = h_entry + return + end if + si%z(1:4) = yend + explicit_h_carry = hlast + call explicit_step_finish(h_entry, ierr) + if (ierr /= SYMPLECTIC_STEP_OK) return + ktau = ktau+1 + end do + call sync_field_to_state +end subroutine orbit_timestep_tdrk24a + + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! +! 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 error and its +! field-evaluation count are not commensurable with the methods above. This +! mode puts the classical adaptive workhorse on exactly the same f_ode as the +! symplectic schemes, Gauss-Radau, Bulirsch-Stoer and TDRK, which is what makes +! the work-precision comparison a comparison of methods rather than of +! formulations. +subroutine orbit_timestep_cashkarp45(ierr) + ! + use fortnum_ode, only : ode_integrate, ode_problem_t, ode_workspace_t, & + ode_solution_t + use fortnum_status, only : fortnum_status_t, FORTNUM_OK + integer, intent(out) :: ierr + integer :: ktau, nlast + type(ode_problem_t) :: problem + type(ode_workspace_t) :: workspace + type(ode_solution_t) :: solution + type(fortnum_status_t) :: status + real(dp) :: h_entry + + ierr = 0 + ktau = 0 + problem%rhs => f_ode_fortnum + problem%rtol = si%rtol + problem%atol = si%atol + allocate(problem%y0(4)) + do while(ktau .lt. si%ntau) + h_entry = explicit_h_carry + problem%t0 = ktau*si%dt + problem%t1 = (ktau+1)*si%dt + problem%y0 = si%z(1:4) + problem%h0 = explicit_h_carry + call ode_integrate(problem, workspace, solution, status) + if (status%code /= FORTNUM_OK) then + ierr = SYMPLECTIC_STEP_MAXITER + explicit_h_carry = h_entry + return + end if + nlast = solution%nsteps + 1 + si%z(1:4) = solution%y(:, nlast) + ! solution%h is trimmed to nsteps entries, one FEWER than solution%t and + ! solution%y, because a step size belongs to the interval between two + ! recorded points rather than to a point. ode_integrate_radau and + ! ode_integrate_gbs allocate nsteps+1 for both, so the same expression is + ! in bounds there; here nlast would read one past the end and carry + ! whatever that memory held into the next macro-step's initial step size. + if (solution%nsteps > 0) explicit_h_carry = solution%h(solution%nsteps) + call explicit_step_finish(h_entry, ierr) + if (ierr /= SYMPLECTIC_STEP_OK) return + ktau = ktau+1 + end do + call sync_field_to_state +end subroutine orbit_timestep_cashkarp45 + end module orbit_symplectic_quasi diff --git a/src/simple.f90 b/src/simple.f90 index b2e90c85..001540b3 100644 --- a/src/simple.f90 +++ b/src/simple.f90 @@ -30,7 +30,7 @@ module simple real(dp) :: dtau, dtaumin, v0 integer :: n_e, n_d - integer :: integmode = 0 ! 0 = RK, 1 = Euler1, 2 = Euler2, 3 = Midpoint, 4-7 = Gauss1-4, 15 = Lobatto3 + integer :: integmode = 0 ! 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) real(dp) :: relerr type(field_can_t) :: f diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index e2177f12..f98fb434 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -426,6 +426,16 @@ target_link_libraries(test_orbits.x simple) add_executable (test_collis.x ../../src/test_collis.f90) target_link_libraries(test_collis.x simple) +# Explicit high-order integmodes on the 4D canonical chart: energy drift must +# respond to the tolerance, and the four adaptive methods must agree. Uses the +# analytic test field, so it needs no VMEC file. +add_executable(test_explicit_integmodes.x test_explicit_integmodes.f90) +target_link_libraries(test_explicit_integmodes.x simple) +add_test(NAME test_explicit_integmodes COMMAND test_explicit_integmodes.x) +set_tests_properties(test_explicit_integmodes PROPERTIES + LABELS "unit" + TIMEOUT 300) + add_executable (test_profiles.x test_profiles.f90) target_link_libraries(test_profiles.x simple) add_test(NAME test_profiles COMMAND test_profiles.x) @@ -627,6 +637,18 @@ add_test(NAME test_array_utils COMMAND test_array_utils.x) LABELS "system;python" TIMEOUT 1800) + # The explicit integmodes must notice a marker leaving the plasma. They + # hand a whole substep to an adaptive integrator and see nothing between + # its endpoints, so without an explicit radius check a marker that crossed + # s = 1 was integrated onwards and reported as confined. + add_test(NAME test_explicit_loss_detection + COMMAND ${BOOZER_CHARTMAP_PYTHON} + ${CMAKE_CURRENT_SOURCE_DIR}/test_explicit_loss_detection.py) + set_tests_properties(test_explicit_loss_detection PROPERTIES + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + LABELS "system;python" + TIMEOUT 3600) + # E2E test: VMEC-Boozer vs chartmap confined fractions add_test(NAME test_e2e_boozer_chartmap COMMAND ${BOOZER_CHARTMAP_PYTHON} diff --git a/test/tests/test_explicit_integmodes.f90 b/test/tests/test_explicit_integmodes.f90 new file mode 100644 index 00000000..ada7d72b --- /dev/null +++ b/test/tests/test_explicit_integmodes.f90 @@ -0,0 +1,142 @@ +program test_explicit_integmodes + ! The explicit high-order integmodes on the 4D canonical chart -- Gauss-Radau + ! (20), Bulirsch-Stoer (21), Cash-Karp (22), and two-derivative Runge-Kutta + ! fixed (24) and adaptive (25) -- all integrate the same canonical + ! Hamiltonian through the same f_ode as the symplectic schemes. + ! + ! Two oracles, neither of them a recording of this code: + ! + ! 1. The canonical Hamiltonian H is an exact invariant of the flow in a + ! static field. Every part of its drift is truncation error, so an + ! adaptive method's drift must shrink when the tolerance is tightened, + ! and a fixed-step one's must shrink under refinement. This is what + ! catches a broken error estimate: without control, tightening rtol + ! changes nothing. + ! + ! 2. The four adaptive methods are independent -- collocation, + ! extrapolation, an embedded explicit pair, and a two-derivative pair -- + ! so their agreement at tight tolerance is a real cross-check. Four + ! unrelated methods do not share a bug. + ! + ! The analytic test field is used, so the test needs no VMEC file. + + use, intrinsic :: iso_fortran_env, only: dp => real64, error_unit + use orbit_symplectic, only: symplectic_integrator_t, orbit_sympl_init, & + orbit_timestep_sympl + use orbit_symplectic_base, only: RADAU15, GBS16, CASHKARP45, TDRK24, TDRK24A + use field_can_mod, only: eval_field => evaluate, field_can_from_name, & + field_can_t, field_can_init, get_val + + implicit none + + integer, parameter :: NMODE = 4 + integer, parameter :: MODES(NMODE) = [RADAU15, GBS16, CASHKARP45, TDRK24A] + character(len=16), parameter :: NAMES(NMODE) = & + [character(len=16) :: 'Gauss-Radau', 'Bulirsch-Stoer', 'Cash-Karp', 'TDRK adaptive'] + + real(dp), parameter :: Z0(4) = [0.1_dp, 0.7_dp, 0.1_dp, 0.0_dp] + real(dp), parameter :: VPAR0 = 0.8_dp + real(dp), parameter :: DT = 500.0_dp + integer, parameter :: NSTEP = 100 + + real(dp) :: zloose(4, NMODE), ztight(4, NMODE) + real(dp) :: hloose(NMODE), htight(NMODE) + real(dp) :: err(3), rate, spread + integer :: nfail, i, j, k + integer :: ns(3) + + nfail = 0 + + do i = 1, NMODE + call trace(MODES(i), 1.0e-6_dp, 1, zloose(:, i), hloose(i), nfail, NAMES(i)) + call trace(MODES(i), 1.0e-11_dp, 1, ztight(:, i), htight(i), nfail, NAMES(i)) + write (*, '(a,a16,a,es10.3,a,es10.3)') ' ', NAMES(i), & + ' |dH/H| at rtol 1e-6 ', hloose(i), ' at 1e-11 ', htight(i) + ! Gauss-Radau is 15th order, so on this smooth field it is already at + ! round-off at rtol 1e-6 and has nothing left to gain -- that is the + ! method working, not the control failing. Everything above round-off + ! must respond. + if (.not. (htight(i) < 0.2_dp*hloose(i) .or. hloose(i) < 1.0e-13_dp)) then + write (error_unit, '(a,a)') ' tightening the tolerance did not ', & + 'reduce the energy drift for '//trim(NAMES(i)) + nfail = nfail + 1 + end if + end do + + ! Cross-method agreement at tight tolerance, component by component. + do k = 1, 4 + spread = maxval(ztight(k, :)) - minval(ztight(k, :)) + if (.not. (spread < 1.0e-6_dp*max(1.0_dp, maxval(abs(ztight(k, :)))))) then + write (error_unit, '(a,i0,a,es10.3)') & + ' methods disagree on final z(', k, '), spread ', spread + nfail = nfail + 1 + end if + end do + write (*, '(a,es10.3)') ' cross-method spread at rtol 1e-11, max component ', & + maxval([(maxval(ztight(k, :)) - minval(ztight(k, :)), k = 1, 4)]) + + ! Fixed-step TDRK: no tolerance to tighten, so the invariant has to improve + ! under refinement instead. trace holds the macro-step DT fixed and splits it + ! into ntau substeps, so raising ntau refines rather than lengthening. + ns = [16, 32, 64] + do j = 1, 3 + call trace(TDRK24, 1.0e-6_dp, ns(j), ztight(:, 1), err(j), nfail, 'TDRK fixed') + end do + rate = -1.0_dp + if (err(3) > 0.0_dp .and. err(1) > 0.0_dp) rate = log(err(1)/err(3))/log(4.0_dp) + write (*, '(a,es10.3,a,es10.3,a,es10.3,a,f6.2)') & + ' TDRK fixed |dH/H| at h, h/2, h/4 = ', err(1), ', ', err(2), ', ', & + err(3), ' observed order ', rate + if (.not. (rate > 3.0_dp)) then + write (error_unit, '(a,f6.2)') & + ' fixed-step TDRK energy drift decays below order 3: ', rate + nfail = nfail + 1 + end if + + if (nfail > 0) then + write (error_unit, '(i0,a)') nfail, ' check(s) failed' + stop 1 + end if + write (*, '(a)') 'test_explicit_integmodes: all checks passed' + +contains + + ! Trace NSTEP macro-steps and return the final state and the largest + ! relative excursion of H along the way. + subroutine trace(mode, rtol, ntau, zend, hdrift, nfail, tag) + integer, intent(in) :: mode, ntau + real(dp), intent(in) :: rtol + real(dp), intent(out) :: zend(4) + real(dp), intent(out) :: hdrift + integer, intent(inout) :: nfail + character(*), intent(in) :: tag + + type(symplectic_integrator_t) :: si + type(field_can_t) :: f + real(dp) :: h0 + integer :: it, ierr + + call field_can_init(f, 1.0e-5_dp, 1.0_dp, VPAR0) + call field_can_from_name('test') + call eval_field(f, Z0(1), Z0(2), Z0(3), 0) + call get_val(f, Z0(4)) + + ! si%dt is the SUBSTEP and si%ntau the number of them, so the macro-step + ! is ntau*dt. Dividing here keeps the interval fixed under refinement. + call orbit_sympl_init(si, f, Z0, DT/real(ntau, dp), ntau, rtol, mode) + h0 = f%H + hdrift = 0.0_dp + do it = 1, NSTEP + call orbit_timestep_sympl(si, f, ierr) + if (ierr /= 0) then + write (error_unit, '(a,a,a,i0)') ' ', trim(tag), & + ' step failed with ierr = ', ierr + nfail = nfail + 1 + exit + end if + hdrift = max(hdrift, abs((f%H - h0)/h0)) + end do + zend = si%z(1:4) + end subroutine trace + +end program test_explicit_integmodes diff --git a/test/tests/test_explicit_loss_detection.py b/test/tests/test_explicit_loss_detection.py new file mode 100644 index 00000000..83fd5ad6 --- /dev/null +++ b/test/tests/test_explicit_loss_detection.py @@ -0,0 +1,220 @@ +#!/usr/bin/env python3 +"""The explicit integmodes must notice when a marker leaves the plasma. + +The symplectic schemes test the radius of every Newton iterate, so a crossing +of the last closed flux surface is caught inside the step. The explicit modes +(20 Gauss-Radau, 21 Bulirsch-Stoer, 22 Cash-Karp, 24/25 TDRK) hand a whole +substep to an adaptive integrator in one call and see nothing between its +endpoints. Until they reported a domain exit, a marker that crossed s = 1 was +integrated onwards through extrapolated splines and recorded as confined: at +paper scale that put 44 of 1000 markers outside the plasma, one as far out as +s = 2.8, and the explicit modes reported 76 losses where the symplectic +schemes and RK4/5 both reported about 207. + +The oracle is a physical invariant, not a comparison against a recorded run: + + a marker reported as confined is inside the last closed flux surface. + +s > 1 is outside the plasma by definition, so a confined marker sitting there +is a contradiction no tolerance can excuse. This needs no reference +trajectory, and one escaped marker is enough to fail it -- which matters, +because a short trace produces few losses and a test that needed many would be +insensitive. + +A cross-method check on the loss count is included as a second, weaker signal: +it depends on the symplectic run being right, but it catches systematic +under-detection that the invariant alone would miss if markers happened to +drift back inside before the trace ended. +""" + +import shutil +import subprocess +import sys +from pathlib import Path + +import numpy as np + +SCRIPT_DIR = Path(__file__).resolve().parent +BUILD_DIR = SCRIPT_DIR.parent.parent / "build" +SIMPLE_X = BUILD_DIR / "simple.x" +WOUT = SCRIPT_DIR.parent / "test_data" / "wout.nc" +WORKDIR = BUILD_DIR / "explicit_loss_detection" + +NPART = 64 +TRACE = 1e-2 + +# Two launch surfaces, because the explicit drivers have to re-impose two +# separate chart invariants that the symplectic schemes get from their Newton +# solve. +# +# 0.75 the paper's launch surface, far enough out that a short trace loses +# markers, which is what exercises the boundary check. +# 0.05 near the axis, which exercises the axis chart switch +# (r < 0 -> (|r|, theta + pi), issue #370). Integmode 21 records two +# r_negative events on this configuration, so the path is live; without +# the switch those markers carry a negative radius into the next +# substep and the field is evaluated outside the spline domain. +# +# Be clear about what this surface does and does not test. It runs the +# chart-switch path and asserts the run completes without faults or +# escapes, which is a smoke check. It does NOT discriminate whether the +# switch is applied: removing it changes one marker of 64 by 4.8 in the +# final state, while two CORRECT methods already differ by up to 5.6 on +# the same markers, because any trace long enough for an orbit to reach +# the axis is long enough for the orbits to have decorrelated. No +# cross-method tolerance separates those two numbers. The switch itself +# was verified by comparing per-marker final states with it on and off +# against an otherwise identical build; that evidence lives in the +# commit, not here. A sensitive regression would need an axis-crossing +# orbit in the analytic test field, where the trajectory is not chaotic +# -- worth adding, not attempted here. +# +# contr_pp is forced far negative so no particle is skipped as +# certainly-confined without ever being integrated. +LAUNCH_SURFACES = [0.75, 0.05] +SIMPLE_IN = """\ +&config +trace_time = {trace}d0 +ntimstep = 100 +ntestpart = {npart} +sbeg = {sbeg}d0 +contr_pp = -1d10 +netcdffile = 'wout.nc' +isw_field_type = 2 +deterministic = .True. +startmode = {startmode} +integmode = {integmode} +relerr = {relerr} +npoiper2 = 32 +facE_al = 1.0d0 +/ +""" + +# Reference plus every explicit mode. The reference generates the start file; +# the rest reuse it, so all of them trace the same markers. +# +# rtol 1e-6 for the adaptive modes. What is under test is whether a crossing of +# the last closed flux surface is noticed, not how accurately the orbit is +# integrated, and at 1e-6 all five explicit modes already agree with the +# symplectic reference on the loss count to within one marker out of 1000. A +# tighter tolerance buys nothing here and costs a great deal: Gauss-Radau at +# 1e-8 alone took about half an hour of the suite's time. Integmode 24 is +# fixed-step, so its tolerance is inert. +REFERENCE = ("sympl", 3, "1d-13") +EXPLICIT = [ + ("radau", 20, "1d-6"), + ("gbs", 21, "1d-6"), + ("cashkarp", 22, "1d-6"), + ("tdrk", 24, "1d-13"), + ("tdrk_adaptive", 25, "1d-6"), +] + + +def run(case, integmode, relerr, startmode, start_src, sbeg): + d = WORKDIR / case + if d.exists(): + shutil.rmtree(d) + d.mkdir(parents=True) + (d / "wout.nc").symlink_to(WOUT) + (d / "simple.in").write_text(SIMPLE_IN.format( + trace=TRACE, npart=NPART, integmode=integmode, relerr=relerr, + startmode=startmode, sbeg=sbeg)) + if start_src is not None: + shutil.copy2(start_src, d / "start.dat") + res = subprocess.run([str(SIMPLE_X)], cwd=d, capture_output=True, + text=True, timeout=1800) + if res.returncode != 0: + print(res.stdout[-2000:]) + print(res.stderr[-2000:]) + sys.exit(f"{case}: simple.x failed") + return d + + +def read(d): + """(loss time, final normalised flux) per marker.""" + data = np.loadtxt(d / "times_lost.dat") + return data[:, 1], data[:, 5] + + +def classify(tlost): + """Confined markers reached the end of the trace; lost ones did not. + + SIMPLE writes times_lost = trace_time for a marker that survives rather + than a sentinel, so the comparison has to be against the trace time. NaN + marks an integrator fault, which is neither. + """ + finite = ~np.isnan(tlost) + confined = finite & (tlost >= 0.999999 * TRACE) + lost = finite & (tlost > 0) & (tlost < 0.999999 * TRACE) + return confined, lost, int(np.sum(~finite)) + + +def check_surface(sbeg, failures): + """Run the reference and every explicit mode from one launch surface.""" + tag = f"s{sbeg:g}".replace(".", "") + name, mode, relerr = REFERENCE + d_ref = run(f"{tag}_{name}", mode, relerr, startmode=1, start_src=None, + sbeg=sbeg) + tl_ref, s_ref = read(d_ref) + _, lost_ref, _ = classify(tl_ref) + n_lost_ref = int(np.sum(lost_ref)) + print(f" sbeg={sbeg:<5g} {name:14s} lost {n_lost_ref:3d}/{NPART} " + f"max final s {np.nanmax(s_ref):.4f}") + + # Monte-Carlo tolerance on a loss count of N markers: loose enough that a + # genuine resolution difference between correct methods passes, far tighter + # than the 3x gap the missing domain check produced. + tol = max(4, int(round(3.0 * np.sqrt(max(n_lost_ref, 1))))) + + for name, mode, relerr in EXPLICIT: + d = run(f"{tag}_{name}", mode, relerr, startmode=2, + start_src=d_ref / "start.dat", sbeg=sbeg) + tlost, s = read(d) + confined, lost, faults = classify(tlost) + n_lost = int(np.sum(lost)) + n_escaped = int(np.sum(confined & (s > 1.0 + 1e-6))) + worst = float(np.nanmax(s[confined])) if np.any(confined) else 0.0 + + print(f" sbeg={sbeg:<5g} {name:14s} lost {n_lost:3d}/{NPART} " + f"max final s of confined {worst:.4f} " + f"escaped {n_escaped} faulted {faults}") + + if n_escaped > 0: + failures.append( + f"sbeg={sbeg:g} {name} (integmode {mode}): {n_escaped} " + f"marker(s) reported confined while outside the plasma, worst " + f"at s = {worst:.4f}") + if abs(n_lost - n_lost_ref) > tol: + failures.append( + f"sbeg={sbeg:g} {name} (integmode {mode}): {n_lost} losses " + f"against the symplectic reference's {n_lost_ref}, outside the " + f"tolerance of {tol}") + if faults > 0: + failures.append( + f"sbeg={sbeg:g} {name} (integmode {mode}): {faults} orbit(s) " + f"faulted") + return n_lost_ref + + +def main(): + if not SIMPLE_X.exists() or not WOUT.exists(): + sys.exit("missing simple.x or wout.nc") + + failures = [] + losses = [check_surface(sbeg, failures) for sbeg in LAUNCH_SURFACES] + + # The outer surface must actually lose markers, or the boundary check has + # nothing to detect and a pass would mean nothing. The near-axis surface is + # there for the chart switch and is not required to lose anything. + if losses[0] == 0: + sys.exit("no losses from the outer surface: the test is vacuous") + + if failures: + for line in failures: + print(f"FAIL: {line}", file=sys.stderr) + sys.exit(1) + print("test_explicit_loss_detection: all checks passed") + + +if __name__ == "__main__": + main() diff --git a/test/tests/test_fo_boris.f90 b/test/tests/test_fo_boris.f90 index 6efdc85d..9916e857 100644 --- a/test/tests/test_fo_boris.f90 +++ b/test/tests/test_fo_boris.f90 @@ -15,7 +15,8 @@ program test_fo_boris use parmot_mod, only: ro0 use simple, only: init_params, orbit_timestep_fo_bridge, tracer_t use simple_main, only: init_field - use orbit_fo_boris, only: fo_state_t, fo_init, fo_step, & + use orbit_fo_boris, only: fo_state_t, fo_init, fo_step, fo_step_rkng, & + fo_step_rkng_adaptive, & fo_energy, fo_mu, fo_to_gc, accept_or_fail, FO_OK, FO_LOCATE_FAIL use orbit_fo_field, only: fo_eval_field use reference_coordinates, only: ref_coords @@ -58,6 +59,16 @@ program test_fo_boris call run_fo([0.5_dp, 0.5_dp, 0.2_dp, 1.0_dp, 0.2_dp], ro0_bar, 'trapped', nfail) call run_fo([0.04_dp, 0.5_dp, 0.2_dp, 1.0_dp, 0.7_dp], ro0_bar, 'near-axis', nfail) + ! Runge-Kutta-Nystrom on the same orbits. Full orbit is genuinely + ! y'' = f(t,y,y'), so RKNG applies directly with no reformulation -- the one + ! place a Nystrom-family method fits SIMPLE without a transplant. + call run_fo_rkng([0.5_dp, 0.5_dp, 0.2_dp, 1.0_dp, 0.9_dp], ro0_bar, 'passing', nfail) + call run_fo_rkng([0.5_dp, 0.5_dp, 0.2_dp, 1.0_dp, 0.2_dp], ro0_bar, 'trapped', nfail) + call run_fo_rkng_adaptive([0.5_dp, 0.5_dp, 0.2_dp, 1.0_dp, 0.9_dp], ro0_bar, & + 'passing', nfail) + call run_fo_rkng_adaptive([0.5_dp, 0.5_dp, 0.2_dp, 1.0_dp, 0.2_dp], ro0_bar, & + 'trapped', nfail) + ! A marker exiting the boundary must be located (so the guiding-centre loss test ! runs), never turned into a confined fault. call test_accept_classification(nfail) @@ -165,6 +176,146 @@ subroutine expect_status(got, want, tag, nfail) end if end subroutine expect_status + ! RKNG on the same orbits, gated on ORDER rather than on an energy bound. + ! + ! Comparing RKNG's energy drift against Boris's would be a category error. + ! Boris conserves energy structurally -- its magnetic rotation is exact for + ! constant B, so the drift is round-off at any step size and says nothing about + ! accuracy. For a non-structure-preserving method the energy error IS the + ! truncation error, so the meaningful check is that it falls at the design rate + ! when the step is refined. In a static magnetic field the exact flow conserves + ! energy identically, which makes this an independent oracle: no reference + ! trajectory and no recorded output are involved. + ! + ! The tableau is order 3, so halving the step should cut the energy error by + ! about 8. The gate is a rate above 2.5 across the refinement, loose enough for + ! the pre-asymptotic wobble and tight enough to fail an order defect. + subroutine run_fo_rkng(z0, ro0_bar, tag, nfail) + real(dp), intent(in) :: z0(5), ro0_bar + character(*), intent(in) :: tag + integer, intent(inout) :: nfail + real(dp) :: err(3), rate + integer :: k, nsub(3) + + nsub = [1, 2, 4] + do k = 1, 3 + call rkng_energy_error(z0, ro0_bar, nsub(k), err(k), nfail, tag) + end do + + rate = -1.0_dp + if (err(3) > 0.0_dp .and. err(1) > 0.0_dp) & + rate = log(err(1)/err(3))/log(4.0_dp) + + print '(a,a,a,es10.2,a,es10.2,a,es10.2,a,f6.2)', ' rkng ', tag, & + ' |dE/E0| at h, h/2, h/4 = ', err(1), ', ', err(2), ', ', err(3), & + ' observed order = ', rate + call check('rkng '//tag//' energy error converges at order >= 2.5', & + rate > 2.5_dp, nfail) + end subroutine run_fo_rkng + + subroutine rkng_energy_error(z0, ro0_bar, nsub, emax, nfail, tag) + real(dp), intent(in) :: z0(5), ro0_bar + integer, intent(in) :: nsub + real(dp), intent(out) :: emax + integer, intent(inout) :: nfail + character(*), intent(in) :: tag + type(fo_state_t) :: st + real(dp) :: bmod, mu, vpar_bar, vperp0, E0, E + real(dp) :: s, th, ph, vpar + real(dp) :: Acov(3), dA(3,3), dBmod(3), hcov(3) + integer :: it, ierr, nstep, lost + + call fo_eval_field([sqrt(z0(1)), z0(2), z0(3)], Acov, dA, bmod, dBmod, hcov) + mu = 0.5_dp*z0(4)**2*(1.0_dp - z0(5)**2)/bmod*2.0_dp + vpar_bar = z0(4)*z0(5)*sqrt(2.0_dp) + vperp0 = sqrt(max(2.0_dp*mu*bmod, 0.0_dp)) + + nstep = 2000 + call fo_init(st, z0(1:3), vpar_bar, vperp0, mu, 1.0_dp, & + 1.0_dp, dtaumin/sqrt(2.0_dp), ro0_bar, z0(4)) + E0 = fo_energy(st); emax = 0.0_dp; lost = 0 + do it = 1, nstep + call fo_step_rkng(st, ierr, nsub) + if (ierr /= 0) then; lost = 1; exit; end if + call fo_to_gc(st, s, th, ph, vpar, ierr) + if (ierr /= 0) then; lost = 1; exit; end if + if (s <= 0.0_dp .or. s >= 1.0_dp) exit + E = fo_energy(st) + emax = max(emax, abs((E - E0)/E0)) + end do + call check('rkng '//tag//' step never failed', lost == 0, nfail) + end subroutine rkng_energy_error + + ! Error control has to make the achieved error track the REQUESTED tolerance, + ! which a fixed step cannot do and a broken embedded estimate would not do + ! either -- it would still give small errors at small rtol. So the gate is the + ! log-log slope of energy drift against rtol, which must be near 1, plus the + ! requirement that a tighter tolerance actually costs more force evaluations. + ! + ! The oracle is exact: a static magnetic force does no work, so every part of + ! the energy drift is truncation error. + subroutine run_fo_rkng_adaptive(z0, ro0_bar, tag, nfail) + real(dp), intent(in) :: z0(5), ro0_bar + character(*), intent(in) :: tag + integer, intent(inout) :: nfail + real(dp) :: tols(3), err(3), slope + integer :: k, nfev(3) + + tols = [1.0e-5_dp, 1.0e-7_dp, 1.0e-9_dp] + do k = 1, 3 + call rkng_adaptive_energy_error(z0, ro0_bar, tols(k), err(k), nfev(k), & + nfail, tag) + end do + + slope = -1.0_dp + if (err(3) > 0.0_dp .and. err(1) > 0.0_dp) & + slope = log10(err(1)/err(3))/log10(tols(1)/tols(3)) + + print '(a,a,a,es10.2,a,es10.2,a,es10.2,a,f6.2,a,i0)', ' rkng-adaptive ', & + tag, ' |dE/E0| at rtol 1e-5, 1e-7, 1e-9 = ', err(1), ', ', err(2), ', ', & + err(3), ' err/tol slope = ', slope, ' nfev = ', nfev(3) + call check('rkng-adaptive '//tag//' error tracks the requested tolerance', & + slope > 0.6_dp .and. slope < 1.4_dp, nfail) + call check('rkng-adaptive '//tag//' tighter tolerance costs more work', & + nfev(3) > nfev(1), nfail) + end subroutine run_fo_rkng_adaptive + + subroutine rkng_adaptive_energy_error(z0, ro0_bar, rtol, emax, nfev_total, & + nfail, tag) + real(dp), intent(in) :: z0(5), ro0_bar, rtol + real(dp), intent(out) :: emax + integer, intent(out) :: nfev_total + integer, intent(inout) :: nfail + character(*), intent(in) :: tag + type(fo_state_t) :: st + real(dp) :: bmod, mu, vpar_bar, vperp0, E0, E, h_carry + real(dp) :: s, th, ph, vpar + real(dp) :: Acov(3), dA(3,3), dBmod(3), hcov(3) + integer :: it, ierr, nstep, lost, nfev + + call fo_eval_field([sqrt(z0(1)), z0(2), z0(3)], Acov, dA, bmod, dBmod, hcov) + mu = 0.5_dp*z0(4)**2*(1.0_dp - z0(5)**2)/bmod*2.0_dp + vpar_bar = z0(4)*z0(5)*sqrt(2.0_dp) + vperp0 = sqrt(max(2.0_dp*mu*bmod, 0.0_dp)) + + nstep = 2000 + call fo_init(st, z0(1:3), vpar_bar, vperp0, mu, 1.0_dp, & + 1.0_dp, dtaumin/sqrt(2.0_dp), ro0_bar, z0(4)) + E0 = fo_energy(st); emax = 0.0_dp; lost = 0 + h_carry = 0.0_dp; nfev_total = 0 + do it = 1, nstep + call fo_step_rkng_adaptive(st, rtol, rtol*1.0e-3_dp, h_carry, ierr, nfev) + nfev_total = nfev_total + nfev + if (ierr /= 0) then; lost = 1; exit; end if + call fo_to_gc(st, s, th, ph, vpar, ierr) + if (ierr /= 0) then; lost = 1; exit; end if + if (s <= 0.0_dp .or. s >= 1.0_dp) exit + E = fo_energy(st) + emax = max(emax, abs((E - E0)/E0)) + end do + call check('rkng-adaptive '//tag//' step never failed', lost == 0, nfail) + end subroutine rkng_adaptive_energy_error + subroutine run_fo(z0, ro0_bar, tag, nfail) real(dp), intent(in) :: z0(5), ro0_bar character(*), intent(in) :: tag