| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module orbit_symplectic_quasi | ||
| 2 | |||
| 3 | use util, only: pi | ||
| 4 | use field_can_mod, only: eval_field => evaluate, field_can_t, get_derivatives, & | ||
| 5 | get_derivatives2 | ||
| 6 | use orbit_symplectic_base, only: symplectic_integrator_t, multistage_integrator_t, & | ||
| 7 | orbit_timestep_quasi_i, coeff_rk_gauss, coeff_rk_lobatto, f_rk_lobatto, sympl_rmax, & | ||
| 8 | SYMPLECTIC_STEP_OK, SYMPLECTIC_STEP_OUTSIDE_DOMAIN, SYMPLECTIC_STEP_MAXITER | ||
| 9 | use fortnum_multiroot, only: multiroot_hybrids | ||
| 10 | use fortnum_status, only: fortnum_status_t | ||
| 11 | use diag_counters, only: count_event, EVT_R_NEGATIVE | ||
| 12 | |||
| 13 | implicit none | ||
| 14 | |||
| 15 | ! Define real(dp) kind parameter | ||
| 16 | integer, parameter :: dp = kind(1.0d0) | ||
| 17 | save | ||
| 18 | |||
| 19 | logical, parameter :: exact_steps = .False. | ||
| 20 | integer, parameter :: S_MAX_GAUSS = 3 | ||
| 21 | |||
| 22 | type(field_can_t) :: f | ||
| 23 | type(field_can_t) :: fs(S_MAX_GAUSS) | ||
| 24 | type(symplectic_integrator_t) :: si | ||
| 25 | !$omp threadprivate(f, si) | ||
| 26 | |||
| 27 | procedure(orbit_timestep_quasi_i), pointer :: orbit_timestep_quasi => null() | ||
| 28 | |||
| 29 | ! Step size carried across macro-steps by the adaptive explicit integrators. | ||
| 30 | ! Without this each macro-step restarts the step-size search from a default | ||
| 31 | ! guess and spends most of its evaluations re-discovering a step it already | ||
| 32 | ! knew, which inflates the evaluation count by more than an order of magnitude | ||
| 33 | ! and would make any cost comparison meaningless. | ||
| 34 | ! | ||
| 35 | ! The carry belongs to one orbit, not to the thread: leaving it set when the | ||
| 36 | ! next orbit starts makes that orbit's first macro-step depend on which orbit | ||
| 37 | ! the thread happened to run before it, so the same particle gives different | ||
| 38 | ! evaluation counts (and, through the controller's history, slightly different | ||
| 39 | ! trajectories) depending on scheduling. orbit_sympl_init resets it. | ||
| 40 | real(dp) :: explicit_h_carry = 0d0 | ||
| 41 | !$omp threadprivate(explicit_h_carry) | ||
| 42 | |||
| 43 | contains | ||
| 44 | |||
| 45 | 784 | subroutine reset_explicit_step_carry | |
| 46 | 784 | explicit_h_carry = 0d0 | |
| 47 | 784 | end subroutine reset_explicit_step_carry | |
| 48 | |||
| 49 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 50 | ! | ||
| 51 | ! Restore the chart invariants the rest of SIMPLE maintains inside a step, and | ||
| 52 | ! report the status the caller's boundary bisection expects. | ||
| 53 | ! | ||
| 54 | ! The symplectic schemes reach both of these through their Newton solve, which | ||
| 55 | ! inspects every iterate. The explicit drivers hand a whole substep to an | ||
| 56 | ! adaptive integrator in one call and see nothing between its endpoints, so | ||
| 57 | ! both have to be re-imposed here, at the same dtaumin granularity the | ||
| 58 | ! symplectic schemes work at. | ||
| 59 | ! | ||
| 60 | ! Leaving the plasma. A marker that crossed the last closed flux surface used | ||
| 61 | ! to be integrated onwards through extrapolated splines and reported as still | ||
| 62 | ! confined: at a 1 s trace that put 44 of 1000 markers at s > 1, one as far out | ||
| 63 | ! as s = 2.8, which is why the explicit modes found 76 losses where the | ||
| 64 | ! symplectic schemes and RK4/5 both independently found about 207. Returning | ||
| 65 | ! OUTSIDE_DOMAIN hands the step to locate_symplectic_boundary, which bisects on | ||
| 66 | ! si%dt exactly as it does for the symplectic schemes. | ||
| 67 | ! | ||
| 68 | ! Crossing the axis. r is a radius, so r < 0 is not a position but a chart | ||
| 69 | ! artefact of passing through the magnetic axis; every other integrator in | ||
| 70 | ! SIMPLE continues on the opposite ray, (r, theta) -> (|r|, theta + pi), rather | ||
| 71 | ! than the pre-#370 teleport to r = 0.01. Without it the next substep evaluates | ||
| 72 | ! the field at a negative radius, outside the spline domain entirely. | ||
| 73 | ! | ||
| 74 | ! h_entry is the step-size carry as it stood before the substep. A step that | ||
| 75 | ! does not complete must not leave its trial step size behind: the boundary | ||
| 76 | ! bisection re-runs the step up to 64 times on ever smaller fractions of dt, | ||
| 77 | ! and without this the carry would come back at bisection scale and throttle | ||
| 78 | ! every later step of the orbit. | ||
| 79 | 70082 | subroutine explicit_step_finish(h_entry, ierr) | |
| 80 | real(dp), intent(in) :: h_entry | ||
| 81 | integer, intent(out) :: ierr | ||
| 82 | |||
| 83 | ! A non-finite state has to be caught before the comparisons below, because | ||
| 84 | ! every comparison against NaN is false: NaN > sympl_rmax does not fire, so a | ||
| 85 | ! NaN orbit would be waved through as a good step and reported as confined. | ||
| 86 | ! The adaptive drivers usually fail first -- a NaN error norm never satisfies | ||
| 87 | ! err <= 1, so the step is rejected until the size underflows -- but the | ||
| 88 | ! fixed-step TDRK has no error estimate to fail, and nothing else on this path | ||
| 89 | ! inspects the state. | ||
| 90 |
4/6✓ Branch 0 taken 280328 times.
✓ Branch 1 taken 70082 times.
✓ Branch 2 taken 280328 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 70082 times.
|
350410 | if (any(si%z(1:4) /= si%z(1:4))) then |
| 91 | ✗ | ierr = SYMPLECTIC_STEP_MAXITER | |
| 92 | ✗ | explicit_h_carry = h_entry | |
| 93 | ✗ | return | |
| 94 | end if | ||
| 95 | |||
| 96 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 70082 times.
|
70082 | if (si%z(1) < 0d0) then |
| 97 | ✗ | call count_event(EVT_R_NEGATIVE) | |
| 98 | ✗ | si%z(1) = -si%z(1) | |
| 99 | ✗ | si%z(2) = si%z(2) + pi | |
| 100 | end if | ||
| 101 | |||
| 102 | 70082 | ierr = SYMPLECTIC_STEP_OK | |
| 103 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 70001 times.
|
70082 | if (si%z(1) > sympl_rmax) ierr = SYMPLECTIC_STEP_OUTSIDE_DOMAIN |
| 104 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 70001 times.
|
70082 | if (ierr /= SYMPLECTIC_STEP_OK) explicit_h_carry = h_entry |
| 105 | end subroutine explicit_step_finish | ||
| 106 | |||
| 107 | ! | ||
| 108 | ! Wrapper routines for ODEPACK | ||
| 109 | ! | ||
| 110 | |||
| 111 | ✗ | subroutine f_exact_quasi(x, fvec, ctx) | |
| 112 | real(dp), intent(in) :: x(:) | ||
| 113 | real(dp), intent(out) :: fvec(:) | ||
| 114 | class(*), intent(in), optional :: ctx | ||
| 115 | |||
| 116 | type(field_can_t) :: f2 | ||
| 117 | |||
| 118 | ✗ | f2 = f | |
| 119 | |||
| 120 | ✗ | call eval_field(f2, x(1), si%z(2), si%z(3), 0) | |
| 121 | ✗ | call get_derivatives(f2, si%z(4)) | |
| 122 | |||
| 123 | ✗ | fvec(1) = f%pth - f2%pth | |
| 124 | |||
| 125 | ✗ | f%H = f2%H | |
| 126 | |||
| 127 | ✗ | end subroutine f_exact_quasi | |
| 128 | |||
| 129 | |||
| 130 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 3009 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3009 times.
|
3009 | subroutine f_euler1_quasi(x, fvec, ctx) |
| 131 | real(dp), intent(in) :: x(:) | ||
| 132 | real(dp), intent(out) :: fvec(:) | ||
| 133 | class(*), intent(in), optional :: ctx | ||
| 134 | |||
| 135 | 3009 | call eval_field(f, x(1), si%z(2), si%z(3), 0) | |
| 136 | 3009 | call get_derivatives(f, x(2)) | |
| 137 | |||
| 138 | fvec(1) = f%dpth(1)*(f%pth - si%pthold) & | ||
| 139 | 3009 | + si%dt*(f%dH(2)*f%dpth(1) - f%dH(1)*f%dpth(2)) | |
| 140 | fvec(2) = f%dpth(1)*(x(2) - si%z(4)) & | ||
| 141 | 3009 | + si%dt*(f%dH(3)*f%dpth(1) - f%dH(1)*f%dpth(3)) | |
| 142 | |||
| 143 | 3009 | end subroutine f_euler1_quasi | |
| 144 | |||
| 145 | |||
| 146 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 184416 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 184416 times.
|
184416 | subroutine f_euler2_quasi(x, fvec, ctx) |
| 147 | real(dp), intent(in) :: x(:) | ||
| 148 | real(dp), intent(out) :: fvec(:) | ||
| 149 | class(*), intent(in), optional :: ctx | ||
| 150 | |||
| 151 | 184416 | call eval_field(f, x(1), x(2), x(3), 0) | |
| 152 | 184416 | call get_derivatives(f, si%z(4)) | |
| 153 | |||
| 154 | 184416 | fvec(1) = f%pth - si%pthold | |
| 155 | 184416 | fvec(2) = f%dpth(1)*(x(2) - si%z(2)) - si%dt*f%dH(1) | |
| 156 | fvec(3) = f%dpth(1)*f%hph*(x(3) - si%z(3)) & | ||
| 157 | 184416 | - si%dt*(f%dpth(1)*f%vpar - f%dH(1)*f%hth) | |
| 158 | |||
| 159 | 184416 | end subroutine f_euler2_quasi | |
| 160 | |||
| 161 | |||
| 162 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 4565 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4565 times.
|
4565 | subroutine f_midpoint_quasi(x, fvec, ctx) |
| 163 | real(dp), intent(in) :: x(:) | ||
| 164 | real(dp), intent(out) :: fvec(:) | ||
| 165 | class(*), intent(in), optional :: ctx | ||
| 166 | |||
| 167 | real(dp) :: dpthmid, pthdotbar | ||
| 168 | |||
| 169 | 4565 | call eval_field(f, x(5), 0.5*(x(2) + si%z(2)), 0.5*(x(3) + si%z(3)), 0) | |
| 170 | 4565 | call get_derivatives(f, 0.5*(x(4) + si%z(4))) | |
| 171 | |||
| 172 | 4565 | fvec(2) = f%dpth(1)*(x(2) - si%z(2)) - si%dt*f%dH(1) | |
| 173 | fvec(3) = f%dpth(1)*f%hph*(x(3) - si%z(3)) & | ||
| 174 | 4565 | - si%dt*(f%dpth(1)*f%vpar - f%dH(1)*f%hth) | |
| 175 | fvec(4) = f%dpth(1)*(x(4) - si%z(4)) & | ||
| 176 | 4565 | + si%dt*(f%dH(3)*f%dpth(1) - f%dH(1)*f%dpth(3)) | |
| 177 | fvec(5) = f%dpth(1)*(f%pth - si%pthold) & | ||
| 178 | 4565 | + si%dt/2.0d0*(f%dpth(1)*f%dH(2)-f%dpth(2)*f%dH(1)) | |
| 179 | |||
| 180 | ! save evaluation from midpoint | ||
| 181 | 4565 | dpthmid = f%dpth(1) | |
| 182 | 4565 | pthdotbar = f%dpth(1)*f%dH(2) - f%dpth(2)*f%dH(1) | |
| 183 | |||
| 184 | ! evaluate at endpoint | ||
| 185 | 4565 | call eval_field(f, x(1), x(2), x(3), 0) | |
| 186 | 4565 | call get_derivatives(f, x(4)) | |
| 187 | 4565 | fvec(1) = dpthmid*(f%pth - si%pthold) + si%dt*pthdotbar | |
| 188 | |||
| 189 | 4565 | end subroutine f_midpoint_quasi | |
| 190 | |||
| 191 | |||
| 192 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 14820 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 14820 times.
|
14820 | subroutine f_rk_gauss_quasi(x, fvec, ctx) |
| 193 | |||
| 194 | real(dp), intent(in) :: x(:) ! = (rend, thend, phend, pphend, rmid) | ||
| 195 | real(dp), intent(out) :: fvec(:) | ||
| 196 | class(*), intent(in), optional :: ctx | ||
| 197 | |||
| 198 | integer :: n | ||
| 199 | 29640 | real(dp) :: a(size(x)/4,size(x)/4), b(size(x)/4), c(size(x)/4), Hprime(size(x)/4) | |
| 200 | integer :: k,l ! counters | ||
| 201 | |||
| 202 | 14820 | n = size(x) | |
| 203 | 14820 | call coeff_rk_gauss(n/4, a, b, c) ! TODO: move this to preprocessing | |
| 204 | |||
| 205 | ! evaluate stages | ||
| 206 |
2/2✓ Branch 0 taken 14820 times.
✓ Branch 1 taken 34200 times.
|
49020 | do k = 1, n/4 |
| 207 | 34200 | call eval_field(fs(k), x(4*k-3), x(4*k-2), x(4*k-1), 0) | |
| 208 | 34200 | call get_derivatives(fs(k), x(4*k)) | |
| 209 | 49020 | Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) | |
| 210 | end do | ||
| 211 | |||
| 212 |
2/2✓ Branch 0 taken 34200 times.
✓ Branch 1 taken 14820 times.
|
49020 | do k = 1, n/4 |
| 213 | 34200 | fvec(4*k-3) = fs(k)%pth - si%pthold | |
| 214 | 34200 | fvec(4*k-2) = x(4*k-2) - si%z(2) | |
| 215 | 34200 | fvec(4*k-1) = x(4*k-1) - si%z(3) | |
| 216 | 34200 | fvec(4*k) = x(4*k) - si%z(4) | |
| 217 |
2/2✓ Branch 0 taken 87400 times.
✓ Branch 1 taken 34200 times.
|
136420 | do l = 1, n/4 |
| 218 | fvec(4*k-3) = fvec(4*k-3) & | ||
| 219 | 87400 | + si%dt*a(k,l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) ! pthdot | |
| 220 | 87400 | fvec(4*k-2) = fvec(4*k-2) - si%dt*a(k,l)*Hprime(l) ! thdot | |
| 221 | fvec(4*k-1) = fvec(4*k-1) & | ||
| 222 | 87400 | - si%dt*a(k,l)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph ! phdot | |
| 223 | fvec(4*k) = fvec(4*k) & | ||
| 224 | 121600 | + si%dt*a(k,l)*(fs(l)%dH(3) - Hprime(l)*fs(l)%dpth(3)) ! pphdot | |
| 225 | end do | ||
| 226 | end do | ||
| 227 | |||
| 228 | 14820 | end subroutine f_rk_gauss_quasi | |
| 229 | |||
| 230 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 231 | ! | ||
| 232 | ! Lobatto (IIIA)-(IIIB) Runge-Kutta method, s stages (n=4*s-2 variables) | ||
| 233 | ! | ||
| 234 | ✗ | subroutine f_rk_lobatto_quasi(x, fvec, ctx) | |
| 235 | ! | ||
| 236 | real(dp), intent(in) :: x(:) | ||
| 237 | real(dp), intent(out) :: fvec(:) | ||
| 238 | class(*), intent(in), optional :: ctx | ||
| 239 | |||
| 240 | ✗ | call f_rk_lobatto(si, fs, (size(x)+2)/4, x, fvec, 0) | |
| 241 | |||
| 242 | ✗ | end subroutine f_rk_lobatto_quasi | |
| 243 | |||
| 244 | |||
| 245 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 246 | ! | ||
| 247 | 380 | subroutine orbit_timestep_multi_quasi(mi, ierr) | |
| 248 | ! | ||
| 249 | type(multistage_integrator_t), intent(inout) :: mi | ||
| 250 | |||
| 251 | integer, intent(out) :: ierr | ||
| 252 | |||
| 253 | integer :: kstage | ||
| 254 | |||
| 255 | 380 | si = mi%stages(1) | |
| 256 | 380 | call orbit_timestep_quasi(ierr) | |
| 257 | |||
| 258 |
2/2✓ Branch 0 taken 380 times.
✓ Branch 1 taken 3610 times.
|
3990 | do kstage = 2, 2*mi%s |
| 259 |
2/2✓ Branch 0 taken 14440 times.
✓ Branch 1 taken 3610 times.
|
18050 | mi%stages(kstage)%z = si%z |
| 260 | 3610 | mi%stages(kstage)%pthold = f%pth | |
| 261 | 3610 | si = mi%stages(kstage) | |
| 262 | 3990 | call orbit_timestep_quasi(ierr) | |
| 263 | end do | ||
| 264 |
2/2✓ Branch 0 taken 1520 times.
✓ Branch 1 taken 380 times.
|
1900 | mi%stages(1)%z = si%z |
| 265 | 380 | mi%stages(1)%pthold = f%pth | |
| 266 | |||
| 267 | 380 | end subroutine orbit_timestep_multi_quasi | |
| 268 | |||
| 269 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 270 | ! | ||
| 271 | 95 | subroutine timestep_midpoint_quasi(ierr) | |
| 272 | ! | ||
| 273 | integer, intent(out) :: ierr | ||
| 274 | |||
| 275 | integer, parameter :: n = 5 | ||
| 276 | |||
| 277 | real(dp), dimension(n) :: x, x0 | ||
| 278 | type(fortnum_status_t) :: status | ||
| 279 | integer :: ktau | ||
| 280 | |||
| 281 | 95 | ierr = 0 | |
| 282 | 95 | ktau = 0 | |
| 283 |
2/2✓ Branch 0 taken 95 times.
✓ Branch 1 taken 95 times.
|
190 | do while(ktau .lt. si%ntau) |
| 284 | 95 | si%pthold = f%pth | |
| 285 | |||
| 286 |
2/2✓ Branch 0 taken 380 times.
✓ Branch 1 taken 95 times.
|
475 | x0(1:4) = si%z |
| 287 | 95 | x0(5) = si%z(1) | |
| 288 | |||
| 289 | call multiroot_hybrids(f_midpoint_quasi, n, x0, x, status, & | ||
| 290 | 95 | xtol=si%rtol, ftol=si%rtol) | |
| 291 | |||
| 292 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 95 times.
|
95 | if (x(1) > sympl_rmax) then |
| 293 | ✗ | ierr = 1 | |
| 294 | ✗ | return | |
| 295 | end if | ||
| 296 | |||
| 297 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 95 times.
|
95 | if (x(1) < 0.0) then |
| 298 | ! Axis crossing: continue on the opposite ray, (r, theta) -> | ||
| 299 | ! (|r|, theta + pi), instead of the former teleport to r = 0.01 (#370). | ||
| 300 | ✗ | call count_event(EVT_R_NEGATIVE) | |
| 301 | ✗ | x(1) = -x(1) | |
| 302 | ✗ | x(2) = x(2) + pi | |
| 303 | ✗ | if (x(1) > sympl_rmax) then | |
| 304 | ✗ | ierr = 1 | |
| 305 | ✗ | return | |
| 306 | end if | ||
| 307 | end if | ||
| 308 | |||
| 309 |
2/2✓ Branch 0 taken 380 times.
✓ Branch 1 taken 95 times.
|
475 | si%z = x(1:4) |
| 310 | |||
| 311 | 95 | ktau = ktau+1 | |
| 312 | enddo | ||
| 313 | |||
| 314 | 95 | end subroutine timestep_midpoint_quasi | |
| 315 | |||
| 316 | |||
| 317 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 318 | ! | ||
| 319 | 95 | subroutine timestep_expl_impl_euler_quasi(ierr) | |
| 320 | ! | ||
| 321 | integer, intent(out) :: ierr | ||
| 322 | |||
| 323 | integer, parameter :: n = 2 | ||
| 324 | |||
| 325 | real(dp), dimension(n) :: x, x0 | ||
| 326 | real(dp), dimension(1) :: xe0, xe | ||
| 327 | type(fortnum_status_t) :: status | ||
| 328 | integer :: ktau | ||
| 329 | |||
| 330 | 95 | ierr = 0 | |
| 331 | 95 | ktau = 0 | |
| 332 |
2/2✓ Branch 0 taken 95 times.
✓ Branch 1 taken 95 times.
|
190 | do while(ktau .lt. si%ntau) |
| 333 | 95 | si%pthold = f%pth | |
| 334 | |||
| 335 | 95 | x0(1)=si%z(1) | |
| 336 | 95 | x0(2)=si%z(4) | |
| 337 | |||
| 338 | call multiroot_hybrids(f_euler1_quasi, n, x0, x, status, & | ||
| 339 | 95 | xtol=si%rtol, ftol=si%rtol) | |
| 340 | |||
| 341 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 95 times.
|
95 | if (x(1) > sympl_rmax) then |
| 342 | ✗ | ierr = 1 | |
| 343 | ✗ | return | |
| 344 | end if | ||
| 345 | |||
| 346 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 95 times.
|
95 | if (x(1) < 0.0) then |
| 347 | ! Axis crossing: continue on the opposite ray, (r, theta) -> | ||
| 348 | ! (|r|, theta + pi), instead of the former teleport to r = 0.01 (#370). | ||
| 349 | ✗ | call count_event(EVT_R_NEGATIVE) | |
| 350 | ✗ | x(1) = -x(1) | |
| 351 | ✗ | x(2) = x(2) + pi | |
| 352 | ✗ | if (x(1) > sympl_rmax) then | |
| 353 | ✗ | ierr = 1 | |
| 354 | ✗ | return | |
| 355 | end if | ||
| 356 | end if | ||
| 357 | |||
| 358 | 95 | si%z(1) = x(1) | |
| 359 | 95 | si%z(4) = x(2) | |
| 360 | |||
| 361 | 95 | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | |
| 362 | 95 | call get_derivatives(f, si%z(4)) | |
| 363 | |||
| 364 | 95 | si%z(2) = si%z(2) + si%dt*f%dH(1)/f%dpth(1) | |
| 365 | 95 | si%z(3) = si%z(3) + si%dt*(f%vpar - f%dH(1)/f%dpth(1)*f%hth)/f%hph | |
| 366 | |||
| 367 | if (exact_steps) then | ||
| 368 | xe0(1) = si%z(1) | ||
| 369 | call multiroot_hybrids(f_exact_quasi, 1, xe0, xe, status, & | ||
| 370 | xtol=si%rtol, ftol=si%rtol) | ||
| 371 | si%z(1) = xe(1) | ||
| 372 | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | ||
| 373 | call get_derivatives(f, si%z(4)) | ||
| 374 | end if | ||
| 375 | |||
| 376 | 95 | ktau = ktau+1 | |
| 377 | enddo | ||
| 378 | |||
| 379 | 95 | end subroutine timestep_expl_impl_euler_quasi | |
| 380 | |||
| 381 | |||
| 382 | |||
| 383 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 384 | ! | ||
| 385 | 4085 | subroutine timestep_impl_expl_euler_quasi(ierr) | |
| 386 | ! | ||
| 387 | integer, intent(out) :: ierr | ||
| 388 | |||
| 389 | integer, parameter :: n = 3 | ||
| 390 | |||
| 391 | real(dp), dimension(n) :: x, x0 | ||
| 392 | real(dp), dimension(1) :: xe0, xe | ||
| 393 | type(fortnum_status_t) :: status | ||
| 394 | integer :: ktau | ||
| 395 | |||
| 396 | 4085 | ierr = 0 | |
| 397 | 4085 | ktau = 0 | |
| 398 |
2/2✓ Branch 0 taken 4085 times.
✓ Branch 1 taken 4085 times.
|
8170 | do while(ktau .lt. si%ntau) |
| 399 | 4085 | si%pthold = f%pth | |
| 400 | |||
| 401 |
2/2✓ Branch 0 taken 12255 times.
✓ Branch 1 taken 4085 times.
|
16340 | x0 = si%z(1:3) |
| 402 | |||
| 403 | call multiroot_hybrids(f_euler2_quasi, n, x0, x, status, & | ||
| 404 | 4085 | xtol=si%rtol, ftol=si%rtol) | |
| 405 | |||
| 406 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4085 times.
|
4085 | if (x(1) > sympl_rmax) then |
| 407 | ✗ | ierr = 1 | |
| 408 | ✗ | return | |
| 409 | end if | ||
| 410 | |||
| 411 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4085 times.
|
4085 | if (x(1) < 0.0) then |
| 412 | ! Axis crossing: continue on the opposite ray, (r, theta) -> | ||
| 413 | ! (|r|, theta + pi), instead of the former teleport to r = 0.01 (#370). | ||
| 414 | ✗ | call count_event(EVT_R_NEGATIVE) | |
| 415 | ✗ | x(1) = -x(1) | |
| 416 | ✗ | x(2) = x(2) + pi | |
| 417 | ✗ | if (x(1) > sympl_rmax) then | |
| 418 | ✗ | ierr = 1 | |
| 419 | ✗ | return | |
| 420 | end if | ||
| 421 | end if | ||
| 422 | |||
| 423 |
2/2✓ Branch 0 taken 12255 times.
✓ Branch 1 taken 4085 times.
|
16340 | si%z(1:3) = x |
| 424 | |||
| 425 | 4085 | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | |
| 426 | 4085 | call get_derivatives(f, si%z(4)) | |
| 427 | |||
| 428 | 4085 | f%pth = si%pthold - si%dt*(f%dH(2) - f%dH(1)*f%dpth(2)/f%dpth(1)) | |
| 429 | 4085 | si%z(4) = si%z(4) - si%dt*(f%dH(3) - f%dH(1)*f%dpth(3)/f%dpth(1)) | |
| 430 | |||
| 431 | if (exact_steps) then | ||
| 432 | xe0(1) = si%z(1) | ||
| 433 | call multiroot_hybrids(f_exact_quasi, 1, xe0, xe, status, & | ||
| 434 | xtol=si%rtol, ftol=si%rtol) | ||
| 435 | si%z(1) = xe(1) | ||
| 436 | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | ||
| 437 | call get_derivatives(f, si%z(4)) | ||
| 438 | end if | ||
| 439 | |||
| 440 | 4085 | ktau = ktau+1 | |
| 441 | enddo | ||
| 442 | |||
| 443 | 4085 | end subroutine timestep_impl_expl_euler_quasi | |
| 444 | |||
| 445 | 285 | subroutine timestep_rk_gauss_quasi(s, ierr) | |
| 446 | ! | ||
| 447 | integer, intent(out) :: ierr | ||
| 448 | |||
| 449 | |||
| 450 | integer, intent(in) :: s | ||
| 451 | 570 | real(dp), dimension(4*s) :: x, x0 | |
| 452 | type(fortnum_status_t) :: status | ||
| 453 | integer :: ktau, k, l | ||
| 454 | |||
| 455 | 285 | real(dp) :: a(s,s), b(s), c(s), Hprime(s) | |
| 456 | |||
| 457 |
2/2✓ Branch 0 taken 570 times.
✓ Branch 1 taken 285 times.
|
855 | do k = 1,s |
| 458 | 855 | fs(k) = f | |
| 459 | end do | ||
| 460 | |||
| 461 | 285 | ierr = 0 | |
| 462 | 285 | ktau = 0 | |
| 463 |
2/2✓ Branch 0 taken 285 times.
✓ Branch 1 taken 285 times.
|
570 | do while(ktau .lt. si%ntau) |
| 464 | 285 | si%pthold = f%pth | |
| 465 | |||
| 466 |
2/2✓ Branch 0 taken 570 times.
✓ Branch 1 taken 285 times.
|
855 | do k = 1,s |
| 467 |
2/2✓ Branch 0 taken 2280 times.
✓ Branch 1 taken 570 times.
|
3135 | x0((4*k-3):(4*k)) = si%z |
| 468 | end do | ||
| 469 | |||
| 470 | call multiroot_hybrids(f_rk_gauss_quasi, 4*s, x0, x, status, & | ||
| 471 | 285 | xtol=si%rtol, ftol=si%rtol) | |
| 472 | |||
| 473 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 285 times.
|
285 | if (x(1) > sympl_rmax) then |
| 474 | ✗ | ierr = 1 | |
| 475 | ✗ | return | |
| 476 | end if | ||
| 477 | |||
| 478 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 285 times.
|
285 | if (x(1) < 0.0) then |
| 479 | ! Axis crossing: continue on the opposite ray, (r, theta) -> | ||
| 480 | ! (|r|, theta + pi), instead of the former teleport to r = 0.01 (#370). | ||
| 481 | ✗ | call count_event(EVT_R_NEGATIVE) | |
| 482 | ✗ | x(1) = -x(1) | |
| 483 | ✗ | x(2) = x(2) + pi | |
| 484 | ✗ | if (x(1) > sympl_rmax) then | |
| 485 | ✗ | ierr = 1 | |
| 486 | ✗ | return | |
| 487 | end if | ||
| 488 | end if | ||
| 489 | |||
| 490 | 285 | call coeff_rk_gauss(s, a, b, c) | |
| 491 | |||
| 492 | 285 | f = fs(s) | |
| 493 | 285 | f%pth = si%pthold | |
| 494 | 285 | si%z(1) = x(4*s-3) | |
| 495 | |||
| 496 |
2/2✓ Branch 0 taken 570 times.
✓ Branch 1 taken 285 times.
|
855 | do l = 1, s |
| 497 | 570 | Hprime(l) = fs(l)%dH(1)/fs(l)%dpth(1) | |
| 498 | 570 | f%pth = f%pth - si%dt*b(l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) | |
| 499 | 570 | si%z(2) = si%z(2) + si%dt*b(l)*Hprime(l) | |
| 500 | 570 | si%z(3) = si%z(3) + si%dt*b(l)*(fs(l)%vpar-Hprime(l)*fs(l)%hth)/fs(l)%hph | |
| 501 | 855 | si%z(4) = si%z(4) - si%dt*b(l)*(fs(l)%dH(3) - Hprime(l)*fs(l)%dpth(3)) | |
| 502 | end do | ||
| 503 | |||
| 504 | 285 | ktau = ktau+1 | |
| 505 | enddo | ||
| 506 | |||
| 507 | 285 | end subroutine timestep_rk_gauss_quasi | |
| 508 | |||
| 509 | |||
| 510 | ✗ | subroutine timestep_rk_lobatto_quasi(s, ierr) | |
| 511 | ! | ||
| 512 | integer, intent(out) :: ierr | ||
| 513 | |||
| 514 | |||
| 515 | integer, intent(in) :: s | ||
| 516 | ✗ | real(dp), dimension(4*s-2) :: x, x0 | |
| 517 | type(fortnum_status_t) :: status | ||
| 518 | integer :: ktau, k, l | ||
| 519 | |||
| 520 | ✗ | real(dp) :: a(s,s), ahat(s,s), b(s), c(s), Hprime(s) | |
| 521 | |||
| 522 | ✗ | do k = 1,s | |
| 523 | ✗ | fs(k) = f | |
| 524 | end do | ||
| 525 | |||
| 526 | ✗ | ierr = 0 | |
| 527 | ✗ | ktau = 0 | |
| 528 | ✗ | do while(ktau .lt. si%ntau) | |
| 529 | ✗ | si%pthold = f%pth | |
| 530 | |||
| 531 | ✗ | x0(1) = si%z(1) | |
| 532 | ✗ | x0(2) = si%z(4) | |
| 533 | ✗ | do k = 2,s | |
| 534 | ✗ | x0((4*k-3-2):(4*k-2)) = si%z | |
| 535 | end do | ||
| 536 | |||
| 537 | ! | ||
| 538 | ! q1 = q0 | ||
| 539 | ! q2 = q0 + h*(a21*f(w1) + a22*f(z2) + a23*f(z3) | ||
| 540 | ! q3 = q0 + h*(a31*f(w1) + a32*f(z2) + a33*f(z3)) | ||
| 541 | ! | ||
| 542 | ! p1(w1) = p0 + h*(o11*f(z1) + o12*f(z2)) | ||
| 543 | ! p2(z2) = p0 + h*(o21*f(z1) + o22*f(z2)) | ||
| 544 | ! p3(z3) = p0 + h*(o31*f(z1) + o32*f(z2)) | ||
| 545 | ! | ||
| 546 | ! 0 0 0 0 0 0 | ||
| 547 | ! x x x x x x | ||
| 548 | ! x x x x x x | ||
| 549 | ! x x x x 0 0 | ||
| 550 | ! x x x x 0 0 | ||
| 551 | ! x x x x 0 0 | ||
| 552 | ! | ||
| 553 | ! effectively | ||
| 554 | ! | ||
| 555 | ! x x x x x | ||
| 556 | ! x x x x x | ||
| 557 | ! x x x 0 0 | ||
| 558 | ! x x x 0 0 | ||
| 559 | ! x x x 0 0 | ||
| 560 | |||
| 561 | call multiroot_hybrids(f_rk_lobatto_quasi, 4*s-2, x0, x, status, & | ||
| 562 | ✗ | xtol=si%rtol, ftol=si%rtol) | |
| 563 | |||
| 564 | ✗ | call coeff_rk_lobatto(s, a, ahat, b, c) | |
| 565 | |||
| 566 | ✗ | call eval_field(fs(1), x(1), si%z(2), si%z(3), 0) | |
| 567 | ✗ | call get_derivatives(fs(1), x(2)) | |
| 568 | ✗ | Hprime(1) = fs(1)%dH(1)/fs(1)%dpth(1) | |
| 569 | |||
| 570 | ✗ | do k = 2, s | |
| 571 | ✗ | call eval_field(fs(k), x(4*k-3-2), x(4*k-2-2), x(4*k-1-2), 0) | |
| 572 | ✗ | call get_derivatives(fs(k), x(4*k-2)) | |
| 573 | ✗ | Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) | |
| 574 | end do | ||
| 575 | |||
| 576 | ✗ | f = fs(s) | |
| 577 | ✗ | f%pth = si%pthold | |
| 578 | ✗ | si%z(1) = x(4*s-3-2) | |
| 579 | |||
| 580 | ✗ | do l = 1, s | |
| 581 | ✗ | f%pth = f%pth - si%dt*b(l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) | |
| 582 | ✗ | si%z(2) = si%z(2) + si%dt*b(l)*Hprime(l) | |
| 583 | ✗ | si%z(3) = si%z(3) + si%dt*b(l)*(fs(l)%vpar-Hprime(l)*fs(l)%hth)/fs(l)%hph | |
| 584 | ✗ | si%z(4) = si%z(4) - si%dt*b(l)*(fs(l)%dH(3) - Hprime(l)*fs(l)%dpth(3)) | |
| 585 | end do | ||
| 586 | |||
| 587 | ✗ | ktau = ktau+1 | |
| 588 | enddo | ||
| 589 | |||
| 590 | ✗ | end subroutine timestep_rk_lobatto_quasi | |
| 591 | |||
| 592 | |||
| 593 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 594 | ! | ||
| 595 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2129351 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2129351 times.
|
2129351 | subroutine f_ode(tau, z, zdot) |
| 596 | ! | ||
| 597 | |||
| 598 | real(dp), intent(in) :: tau | ||
| 599 | real(dp), intent(in) :: z(:) | ||
| 600 | real(dp), intent(out) :: zdot(:) | ||
| 601 | real(dp) :: Hprime | ||
| 602 | |||
| 603 | 2129351 | call eval_field(f, z(1), z(2), z(3), 0) | |
| 604 | 2129351 | call get_derivatives(f, z(4)) | |
| 605 | |||
| 606 | 2129351 | Hprime = f%dH(1)/f%dpth(1) | |
| 607 | |||
| 608 | 2129351 | zdot(1) = -(f%dH(2) - f%hth/f%hph*f%dH(3))/f%dpth(1) | |
| 609 | 2129351 | zdot(2) = Hprime | |
| 610 | 2129351 | zdot(3) = (f%vpar-Hprime*f%hth)/f%hph | |
| 611 | 2129351 | zdot(4) = -(f%dH(3) - Hprime*f%dpth(3)) | |
| 612 | |||
| 613 | 2129351 | end subroutine f_ode | |
| 614 | |||
| 615 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 616 | ! | ||
| 617 | ! Index into the packed symmetric second-derivative arrays d2H, d2pth, d2vpar. | ||
| 618 | ! | ||
| 619 | ! field_can_base documents the order as | ||
| 620 | ! 1:(r,r) 2:(r,th) 3:(r,ph) 4:(th,th) 5:(th,ph) 6:(ph,ph) | ||
| 621 | ! 7:(pph,r) 8:(pph,th) 9:(pph,ph) 10:(pph,pph) | ||
| 622 | ! over the variables (r, th, ph, pph) = (1, 2, 3, 4). | ||
| 623 | 4546760 | pure function d2idx(i, j) result(k) | |
| 624 | integer, intent(in) :: i, j | ||
| 625 | integer :: k | ||
| 626 | integer :: lo, hi | ||
| 627 | |||
| 628 | 4546760 | lo = min(i, j) | |
| 629 | 4546760 | hi = max(i, j) | |
| 630 | |||
| 631 |
2/2✓ Branch 0 taken 1136690 times.
✓ Branch 1 taken 3410070 times.
|
4546760 | if (hi == 4) then |
| 632 | 1136690 | k = 6 + lo ! (r,pph)=7, (th,pph)=8, (ph,pph)=9, (pph,pph)=10 | |
| 633 |
2/2✓ Branch 0 taken 1364028 times.
✓ Branch 1 taken 2046042 times.
|
3410070 | else if (lo == 1) then |
| 634 | k = hi ! (r,r)=1, (r,th)=2, (r,ph)=3 | ||
| 635 |
2/2✓ Branch 0 taken 909352 times.
✓ Branch 1 taken 454676 times.
|
1364028 | else if (lo == 2) then |
| 636 | 909352 | k = hi + 2 ! (th,th)=4, (th,ph)=5 | |
| 637 | else | ||
| 638 | k = 6 ! (ph,ph) | ||
| 639 | end if | ||
| 640 | 4546760 | end function d2idx | |
| 641 | |||
| 642 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 643 | ! | ||
| 644 | ! Second time derivative of the canonical guiding-centre state, | ||
| 645 | ! | ||
| 646 | ! G(z) = F'(z) F(z), | ||
| 647 | ! | ||
| 648 | ! which is what a two-derivative Runge-Kutta method integrates in place of F. | ||
| 649 | ! | ||
| 650 | ! The derivation and its verification live in the study repo at | ||
| 651 | ! derivation/tdrk_guiding_center_G.wl: Mathematica differentiates f_ode's F | ||
| 652 | ! symbolically and compares against this closed form. The symbolic residual is | ||
| 653 | ! exactly zero and a numeric spot-check with nonlinear, fully coupled test | ||
| 654 | ! functions agrees to 3.6e-15. | ||
| 655 | ! | ||
| 656 | ! The cost argument rests on which derivatives appear. That check is in the same | ||
| 657 | ! file and reports second derivatives of H and pth, first of vpar, hth and hph, | ||
| 658 | ! and no third derivatives anywhere -- so G needs exactly one mode_secders = 2 | ||
| 659 | ! evaluation and no extra F evaluations. Note also that only d2H indices 1 to 9 | ||
| 660 | ! are touched: F depends on H solely through dH(1:3), so differentiating it can | ||
| 661 | ! never raise index 10, which get_derivatives2 does not fill. | ||
| 662 | ! | ||
| 663 | ! hth and hph are geometric and carry no pph dependence, hence the explicit zero | ||
| 664 | ! in the fourth slot of their gradients rather than an out-of-bounds read of a | ||
| 665 | ! three-element array. | ||
| 666 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 227338 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 227338 times.
|
227338 | subroutine g_ode(tau, z, g) |
| 667 | real(dp), intent(in) :: tau | ||
| 668 | real(dp), intent(in) :: z(:) | ||
| 669 | real(dp), intent(out) :: g(:) | ||
| 670 | |||
| 671 | real(dp) :: Hprime, dHprime(4), dhth4(4), dhph4(4) | ||
| 672 | real(dp) :: zdot(4), dF(4, 4) | ||
| 673 | integer :: j | ||
| 674 | |||
| 675 | 227338 | call eval_field(f, z(1), z(2), z(3), 2) | |
| 676 | 227338 | call get_derivatives2(f, z(4)) | |
| 677 | |||
| 678 | 227338 | Hprime = f%dH(1)/f%dpth(1) | |
| 679 | |||
| 680 |
2/2✓ Branch 0 taken 682014 times.
✓ Branch 1 taken 227338 times.
|
909352 | dhth4(1:3) = f%dhth |
| 681 | 227338 | dhth4(4) = 0d0 | |
| 682 |
2/2✓ Branch 0 taken 682014 times.
✓ Branch 1 taken 227338 times.
|
909352 | dhph4(1:3) = f%dhph |
| 683 | 227338 | dhph4(4) = 0d0 | |
| 684 | |||
| 685 | 227338 | zdot(1) = -(f%dH(2) - f%hth/f%hph*f%dH(3))/f%dpth(1) | |
| 686 | 227338 | zdot(2) = Hprime | |
| 687 | 227338 | zdot(3) = (f%vpar - Hprime*f%hth)/f%hph | |
| 688 | 227338 | zdot(4) = -(f%dH(3) - Hprime*f%dpth(3)) | |
| 689 | |||
| 690 | ! Hprime is a quotient of two state-dependent quantities, so both numerator | ||
| 691 | ! and denominator contribute. This is the term a hand derivation usually drops. | ||
| 692 |
2/2✓ Branch 0 taken 909352 times.
✓ Branch 1 taken 227338 times.
|
1136690 | do j = 1, 4 |
| 693 | dHprime(j) = (f%d2H(d2idx(1, j))*f%dpth(1) & | ||
| 694 | 1136690 | - f%dH(1)*f%d2pth(d2idx(1, j)))/f%dpth(1)**2 | |
| 695 | end do | ||
| 696 | |||
| 697 |
2/2✓ Branch 0 taken 909352 times.
✓ Branch 1 taken 227338 times.
|
1136690 | do j = 1, 4 |
| 698 | dF(1, j) = -((f%d2H(d2idx(2, j)) & | ||
| 699 | - (dhth4(j)/f%hph - f%hth*dhph4(j)/f%hph**2)*f%dH(3) & | ||
| 700 | - f%hth/f%hph*f%d2H(d2idx(3, j)))*f%dpth(1) & | ||
| 701 | - (f%dH(2) - f%hth/f%hph*f%dH(3))*f%d2pth(d2idx(1, j))) & | ||
| 702 | 909352 | /f%dpth(1)**2 | |
| 703 | |||
| 704 | 909352 | dF(2, j) = dHprime(j) | |
| 705 | |||
| 706 | dF(3, j) = ((f%dvpar(j) - dHprime(j)*f%hth - Hprime*dhth4(j))*f%hph & | ||
| 707 | 909352 | - (f%vpar - Hprime*f%hth)*dhph4(j))/f%hph**2 | |
| 708 | |||
| 709 | dF(4, j) = -(f%d2H(d2idx(3, j)) - dHprime(j)*f%dpth(3) & | ||
| 710 | 1136690 | - Hprime*f%d2pth(d2idx(3, j))) | |
| 711 | end do | ||
| 712 | |||
| 713 |
6/6✓ Branch 0 taken 909352 times.
✓ Branch 1 taken 227338 times.
✓ Branch 2 taken 909352 times.
✓ Branch 3 taken 227338 times.
✓ Branch 4 taken 3637408 times.
✓ Branch 5 taken 909352 times.
|
5683450 | g(1:4) = matmul(dF, zdot) |
| 714 | 227338 | end subroutine g_ode | |
| 715 | |||
| 716 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 717 | ! | ||
| 718 | ! Adapter from g_ode to fortnum's ode_rhs2_t. | ||
| 719 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 227338 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 227338 times.
|
227338 | subroutine g_ode_fortnum(t, y, g, ctx) |
| 720 | real(dp), intent(in) :: t | ||
| 721 | real(dp), intent(in) :: y(:) | ||
| 722 | real(dp), intent(out) :: g(:) | ||
| 723 | class(*), intent(in), optional :: ctx | ||
| 724 | |||
| 725 | 227338 | call g_ode(t, y, g) | |
| 726 | 227338 | end subroutine g_ode_fortnum | |
| 727 | |||
| 728 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 729 | ! | ||
| 730 | ✗ | subroutine orbit_timestep_rk45(ierr) | |
| 731 | ! | ||
| 732 | use odeint_allroutines_sub, only : odeint_allroutines | ||
| 733 | integer, intent(out) :: ierr | ||
| 734 | integer :: ktau | ||
| 735 | |||
| 736 | ✗ | ierr = 0 | |
| 737 | ✗ | ktau = 0 | |
| 738 | ✗ | do while(ktau .lt. si%ntau) | |
| 739 | ✗ | call odeint_allroutines(si%z, 4, ktau*si%dt, (ktau+1)*si%dt, si%rtol, f_ode) | |
| 740 | ✗ | ktau = ktau+1 | |
| 741 | end do | ||
| 742 | ✗ | end subroutine orbit_timestep_rk45 | |
| 743 | |||
| 744 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 745 | ! | ||
| 746 | ! Adapter from f_ode to fortnum's ode_rhs_t, which carries an optional context | ||
| 747 | ! argument that this right-hand side does not need. | ||
| 748 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2059350 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2059350 times.
|
2059350 | subroutine f_ode_fortnum(t, y, dydt, ctx) |
| 749 | real(dp), intent(in) :: t | ||
| 750 | real(dp), intent(in) :: y(:) | ||
| 751 | real(dp), intent(out) :: dydt(:) | ||
| 752 | class(*), intent(in), optional :: ctx | ||
| 753 | |||
| 754 | 2059350 | call f_ode(t, y, dydt) | |
| 755 | 2059350 | end subroutine f_ode_fortnum | |
| 756 | |||
| 757 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 758 | ! | ||
| 759 | ! Gauss-Radau (IAS15-style, 15th order) on the same canonical Hamiltonian the | ||
| 760 | ! symplectic schemes use. | ||
| 761 | 11812 | subroutine orbit_timestep_radau15(ierr) | |
| 762 | ! | ||
| 763 | use fortnum_ode_gauss_radau, only: ode_integrate_radau | ||
| 764 | use fortnum_ode, only: ode_problem_t, ode_solution_t | ||
| 765 | use fortnum_status, only: fortnum_status_t, FORTNUM_OK | ||
| 766 | integer, intent(out) :: ierr | ||
| 767 | integer :: ktau, nlast | ||
| 768 | 11812 | type(ode_problem_t) :: problem | |
| 769 | 11812 | type(ode_solution_t) :: solution | |
| 770 | type(fortnum_status_t) :: status | ||
| 771 | real(dp) :: h_entry | ||
| 772 | |||
| 773 | 11812 | ierr = 0 | |
| 774 | 11812 | ktau = 0 | |
| 775 | 11812 | problem%rhs => f_ode_fortnum | |
| 776 | 11812 | problem%rtol = si%rtol | |
| 777 | 11812 | problem%atol = si%atol | |
| 778 | 11812 | allocate (problem%y0(4)) | |
| 779 |
2/2✓ Branch 0 taken 11812 times.
✓ Branch 1 taken 11805 times.
|
23617 | do while (ktau .lt. si%ntau) |
| 780 | 11812 | h_entry = explicit_h_carry | |
| 781 | 11812 | problem%t0 = ktau*si%dt | |
| 782 | 11812 | problem%t1 = (ktau + 1)*si%dt | |
| 783 |
4/10✓ Branch 0 taken 11812 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 11812 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 47248 times.
✓ Branch 9 taken 11812 times.
|
70872 | problem%y0 = si%z(1:4) |
| 784 | 11812 | problem%h0 = explicit_h_carry | |
| 785 | 11812 | call ode_integrate_radau(problem, solution, status) | |
| 786 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11812 times.
|
11812 | if (status%code /= FORTNUM_OK) then |
| 787 | ✗ | ierr = SYMPLECTIC_STEP_MAXITER | |
| 788 | ✗ | explicit_h_carry = h_entry | |
| 789 | 42 | return | |
| 790 | end if | ||
| 791 | 11812 | nlast = size(solution%t) | |
| 792 |
2/2✓ Branch 0 taken 47248 times.
✓ Branch 1 taken 11812 times.
|
59060 | si%z(1:4) = solution%y(:, nlast) |
| 793 |
1/2✓ Branch 0 taken 11812 times.
✗ Branch 1 not taken.
|
11812 | if (nlast > 1) explicit_h_carry = solution%h(nlast) |
| 794 | 11812 | call explicit_step_finish(h_entry, ierr) | |
| 795 |
2/2✓ Branch 0 taken 11805 times.
✓ Branch 1 taken 7 times.
|
11812 | if (ierr /= SYMPLECTIC_STEP_OK) return |
| 796 | 11805 | ktau = ktau + 1 | |
| 797 | end do | ||
| 798 | 11805 | call sync_field_to_state | |
| 799 |
13/26✗ Branch 0 not taken.
✓ Branch 1 taken 11812 times.
✓ Branch 2 taken 11805 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 11805 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 11805 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 11805 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 11805 times.
✓ Branch 12 taken 11805 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 7 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 7 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 7 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 7 times.
✓ Branch 22 taken 7 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 7 times.
✗ Branch 25 not taken.
|
23666 | end subroutine orbit_timestep_radau15 |
| 800 | |||
| 801 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 802 | ! | ||
| 803 | ! Re-evaluate the field at the converged end-of-step state. | ||
| 804 | ! | ||
| 805 | ! Only z(1:4) are integrated; the caller reconstructs the parallel velocity as | ||
| 806 | ! z(5) = f%vpar/(pabs*sqrt(2)), reading f%vpar off the module-global field_can_t | ||
| 807 | ! that f_ode fills as a side effect. After a multi-stage step that side effect | ||
| 808 | ! belongs to the last interior stage node, not to the end of the step, so z(5) | ||
| 809 | ! carries an error the tolerance cannot remove: measured on the sweep, | ||
| 810 | ! components 1-4 of the final state converged to 5e-9 while z(5) stalled at | ||
| 811 | ! 3.3e-4 and got no better from rtol 1e-8 to 1e-10. | ||
| 812 | ! | ||
| 813 | ! One evaluation per macro-step, against the hundreds a step already costs. | ||
| 814 | 70001 | subroutine sync_field_to_state | |
| 815 | real(dp) :: zdot_discard(4) | ||
| 816 | |||
| 817 | 70001 | call f_ode(0d0, si%z(1:4), zdot_discard) | |
| 818 | 70001 | end subroutine sync_field_to_state | |
| 819 | |||
| 820 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 821 | ! | ||
| 822 | ! Gragg-Bulirsch-Stoer extrapolation on the same canonical Hamiltonian. | ||
| 823 | 11806 | subroutine orbit_timestep_gbs16(ierr) | |
| 824 | ! | ||
| 825 | use fortnum_ode_extrapolation, only: ode_integrate_gbs | ||
| 826 | use fortnum_ode, only: ode_problem_t, ode_solution_t | ||
| 827 | use fortnum_status, only: fortnum_status_t, FORTNUM_OK | ||
| 828 | integer, intent(out) :: ierr | ||
| 829 | integer :: ktau, nlast | ||
| 830 | 11806 | type(ode_problem_t) :: problem | |
| 831 | 11806 | type(ode_solution_t) :: solution | |
| 832 | type(fortnum_status_t) :: status | ||
| 833 | real(dp) :: h_entry | ||
| 834 | |||
| 835 | 11806 | ierr = 0 | |
| 836 | 11806 | ktau = 0 | |
| 837 | 11806 | problem%rhs => f_ode_fortnum | |
| 838 | 11806 | problem%rtol = si%rtol | |
| 839 | 11806 | problem%atol = si%atol | |
| 840 | 11806 | allocate (problem%y0(4)) | |
| 841 |
2/2✓ Branch 0 taken 11806 times.
✓ Branch 1 taken 11797 times.
|
23603 | do while (ktau .lt. si%ntau) |
| 842 | 11806 | h_entry = explicit_h_carry | |
| 843 | 11806 | problem%t0 = ktau*si%dt | |
| 844 | 11806 | problem%t1 = (ktau + 1)*si%dt | |
| 845 |
4/10✓ Branch 0 taken 11806 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 11806 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 47224 times.
✓ Branch 9 taken 11806 times.
|
70836 | problem%y0 = si%z(1:4) |
| 846 | 11806 | problem%h0 = explicit_h_carry | |
| 847 | 11806 | call ode_integrate_gbs(problem, solution, status) | |
| 848 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11806 times.
|
11806 | if (status%code /= FORTNUM_OK) then |
| 849 | ✗ | ierr = SYMPLECTIC_STEP_MAXITER | |
| 850 | ✗ | explicit_h_carry = h_entry | |
| 851 | 54 | return | |
| 852 | end if | ||
| 853 | 11806 | nlast = size(solution%t) | |
| 854 |
2/2✓ Branch 0 taken 47224 times.
✓ Branch 1 taken 11806 times.
|
59030 | si%z(1:4) = solution%y(:, nlast) |
| 855 |
1/2✓ Branch 0 taken 11806 times.
✗ Branch 1 not taken.
|
11806 | if (nlast > 1) explicit_h_carry = solution%h(nlast) |
| 856 | 11806 | call explicit_step_finish(h_entry, ierr) | |
| 857 |
2/2✓ Branch 0 taken 11797 times.
✓ Branch 1 taken 9 times.
|
11806 | if (ierr /= SYMPLECTIC_STEP_OK) return |
| 858 | 11797 | ktau = ktau + 1 | |
| 859 | end do | ||
| 860 | 11797 | call sync_field_to_state | |
| 861 |
13/26✗ Branch 0 not taken.
✓ Branch 1 taken 11806 times.
✓ Branch 2 taken 11797 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 11797 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 11797 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 11797 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 11797 times.
✓ Branch 12 taken 11797 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 9 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 9 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 9 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 9 times.
✓ Branch 22 taken 9 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 9 times.
✗ Branch 25 not taken.
|
23666 | end subroutine orbit_timestep_gbs16 |
| 862 | |||
| 863 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 864 | ! | ||
| 865 | ! Two-derivative Runge-Kutta on the canonical Hamiltonian: the Runge-Kutta- | ||
| 866 | ! Nystrom transplant. | ||
| 867 | ! | ||
| 868 | ! RKN integrates y'' = f(t,y) and does not apply to the guiding-centre system as | ||
| 869 | ! written, which is first order. Differentiating once gives zdd = F'(z)F(z) = | ||
| 870 | ! G(z), which is in special second-order form because zd is itself determined by | ||
| 871 | ! z, and re-synchronising the velocity as F(z) at the start of every step turns | ||
| 872 | ! the scheme into a special two-derivative Runge-Kutta method (Chan & Tsai, | ||
| 873 | ! Numer. Algorithms 53 (2010) 171). Its order conditions coincide with the | ||
| 874 | ! Nystrom ones, which is what lets an RKN tableau be used directly here. | ||
| 875 | ! | ||
| 876 | ! Because the velocity is recomputed rather than propagated, only the position | ||
| 877 | ! order conditions are needed, so order 4 costs two stages instead of three. | ||
| 878 | ! | ||
| 879 | ! Fixed step, like the symplectic schemes: fortnum's TDRK has no embedded error | ||
| 880 | ! estimate yet, so the accuracy parameter is the step count rather than a | ||
| 881 | ! tolerance. npoiper2 therefore sets resolution here exactly as it does for the | ||
| 882 | ! symplectic integrators, which also makes the two directly comparable. | ||
| 883 | 22811 | subroutine orbit_timestep_tdrk24(ierr) | |
| 884 | ! | ||
| 885 | use fortnum_ode_tdrk, only: tdrk_integrate_fixed | ||
| 886 | use fortnum_status, only: fortnum_status_t, FORTNUM_OK | ||
| 887 | integer, intent(out) :: ierr | ||
| 888 | integer :: ktau, nfev_f, nfev_g | ||
| 889 | real(dp) :: yend(4) | ||
| 890 | type(fortnum_status_t) :: status | ||
| 891 | real(dp) :: h_entry | ||
| 892 | |||
| 893 | 22811 | ierr = 0 | |
| 894 | 22811 | ktau = 0 | |
| 895 |
2/2✓ Branch 0 taken 22811 times.
✓ Branch 1 taken 22805 times.
|
45616 | do while (ktau .lt. si%ntau) |
| 896 | 22811 | h_entry = explicit_h_carry | |
| 897 | call tdrk_integrate_fixed(f_ode_fortnum, g_ode_fortnum, & | ||
| 898 | ktau*si%dt, (ktau + 1)*si%dt, si%z(1:4), 1, 4, .false., yend, & | ||
| 899 | 22811 | nfev_f, nfev_g, status) | |
| 900 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 22811 times.
|
22811 | if (status%code /= FORTNUM_OK) then |
| 901 | ✗ | ierr = SYMPLECTIC_STEP_MAXITER | |
| 902 | ✗ | explicit_h_carry = h_entry | |
| 903 | 6 | return | |
| 904 | end if | ||
| 905 |
2/2✓ Branch 0 taken 91244 times.
✓ Branch 1 taken 22811 times.
|
114055 | si%z(1:4) = yend |
| 906 | 22811 | call explicit_step_finish(h_entry, ierr) | |
| 907 |
2/2✓ Branch 0 taken 22805 times.
✓ Branch 1 taken 6 times.
|
22811 | if (ierr /= SYMPLECTIC_STEP_OK) return |
| 908 | 22805 | ktau = ktau + 1 | |
| 909 | end do | ||
| 910 | 22805 | call sync_field_to_state | |
| 911 | 22811 | end subroutine orbit_timestep_tdrk24 | |
| 912 | |||
| 913 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 914 | ! | ||
| 915 | ! Two-derivative Runge-Kutta with error control: the same method as | ||
| 916 | ! orbit_timestep_tdrk24, driven by a tolerance instead of a step count. | ||
| 917 | ! | ||
| 918 | ! The embedded estimate needs the three-stage order-4 tableau. The two-stage | ||
| 919 | ! one, which orbit_timestep_tdrk24 uses, cannot carry one: with c = (0, 1/2) | ||
| 920 | ! the order-3 conditions already force the order-4 weights, so the two | ||
| 921 | ! solutions coincide and their difference is identically zero. One extra G | ||
| 922 | ! evaluation per step buys the estimate. | ||
| 923 | ! | ||
| 924 | ! Both modes are kept. Fixed step is what makes TDRK comparable to the | ||
| 925 | ! symplectic schemes at matched resolution; tolerance control is what makes it | ||
| 926 | ! comparable to Gauss-Radau, Bulirsch-Stoer and Cash-Karp at matched accuracy. | ||
| 927 | 11844 | subroutine orbit_timestep_tdrk24a(ierr) | |
| 928 | ! | ||
| 929 | use fortnum_ode_tdrk, only: tdrk_integrate_adaptive | ||
| 930 | use fortnum_status, only: fortnum_status_t, FORTNUM_OK | ||
| 931 | integer, intent(out) :: ierr | ||
| 932 | integer :: ktau, nfev_f, nfev_g, naccept, nreject | ||
| 933 | real(dp) :: yend(4), hlast | ||
| 934 | type(fortnum_status_t) :: status | ||
| 935 | real(dp) :: h_entry | ||
| 936 | |||
| 937 | 11844 | ierr = 0 | |
| 938 | 11844 | ktau = 0 | |
| 939 |
2/2✓ Branch 0 taken 11844 times.
✓ Branch 1 taken 11794 times.
|
23638 | do while (ktau .lt. si%ntau) |
| 940 | 11844 | h_entry = explicit_h_carry | |
| 941 | call tdrk_integrate_adaptive(f_ode_fortnum, g_ode_fortnum, & | ||
| 942 | ktau*si%dt, (ktau + 1)*si%dt, si%z(1:4), si%rtol, si%atol, & | ||
| 943 | 11844 | explicit_h_carry, yend, hlast, nfev_f, nfev_g, naccept, nreject, status) | |
| 944 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11844 times.
|
11844 | if (status%code /= FORTNUM_OK) then |
| 945 | ✗ | ierr = SYMPLECTIC_STEP_MAXITER | |
| 946 | ✗ | explicit_h_carry = h_entry | |
| 947 | 50 | return | |
| 948 | end if | ||
| 949 |
2/2✓ Branch 0 taken 47376 times.
✓ Branch 1 taken 11844 times.
|
59220 | si%z(1:4) = yend |
| 950 | 11844 | explicit_h_carry = hlast | |
| 951 | 11844 | call explicit_step_finish(h_entry, ierr) | |
| 952 |
2/2✓ Branch 0 taken 11794 times.
✓ Branch 1 taken 50 times.
|
11844 | if (ierr /= SYMPLECTIC_STEP_OK) return |
| 953 | 11794 | ktau = ktau + 1 | |
| 954 | end do | ||
| 955 | 11794 | call sync_field_to_state | |
| 956 | 11844 | end subroutine orbit_timestep_tdrk24a | |
| 957 | |||
| 958 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 959 | ! | ||
| 960 | ! Cash-Karp RK5(4) on the 4D canonical chart. | ||
| 961 | ! | ||
| 962 | ! integmode = 0 already runs a Cash-Karp, but on the 5D drift-kinetic form with | ||
| 963 | ! a different right-hand side and different coordinates, so its error and its | ||
| 964 | ! field-evaluation count are not commensurable with the methods above. This | ||
| 965 | ! mode puts the classical adaptive workhorse on exactly the same f_ode as the | ||
| 966 | ! symplectic schemes, Gauss-Radau, Bulirsch-Stoer and TDRK, which is what makes | ||
| 967 | ! the work-precision comparison a comparison of methods rather than of | ||
| 968 | ! formulations. | ||
| 969 | 11809 | subroutine orbit_timestep_cashkarp45(ierr) | |
| 970 | ! | ||
| 971 | use fortnum_ode, only: ode_integrate, ode_problem_t, ode_workspace_t, & | ||
| 972 | ode_solution_t | ||
| 973 | use fortnum_status, only: fortnum_status_t, FORTNUM_OK | ||
| 974 | integer, intent(out) :: ierr | ||
| 975 | integer :: ktau, nlast | ||
| 976 | 11809 | type(ode_problem_t) :: problem | |
| 977 | 11809 | type(ode_workspace_t) :: workspace | |
| 978 | 11809 | type(ode_solution_t) :: solution | |
| 979 | type(fortnum_status_t) :: status | ||
| 980 | real(dp) :: h_entry | ||
| 981 | |||
| 982 | 11809 | ierr = 0 | |
| 983 | 11809 | ktau = 0 | |
| 984 | 11809 | problem%rhs => f_ode_fortnum | |
| 985 | 11809 | problem%rtol = si%rtol | |
| 986 | 11809 | problem%atol = si%atol | |
| 987 | 11809 | allocate (problem%y0(4)) | |
| 988 |
2/2✓ Branch 0 taken 11809 times.
✓ Branch 1 taken 11800 times.
|
23609 | do while (ktau .lt. si%ntau) |
| 989 | 11809 | h_entry = explicit_h_carry | |
| 990 | 11809 | problem%t0 = ktau*si%dt | |
| 991 | 11809 | problem%t1 = (ktau + 1)*si%dt | |
| 992 |
4/10✓ Branch 0 taken 11809 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 11809 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 47236 times.
✓ Branch 9 taken 11809 times.
|
70854 | problem%y0 = si%z(1:4) |
| 993 | 11809 | problem%h0 = explicit_h_carry | |
| 994 | 11809 | call ode_integrate(problem, workspace, solution, status) | |
| 995 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11809 times.
|
11809 | if (status%code /= FORTNUM_OK) then |
| 996 | ✗ | ierr = SYMPLECTIC_STEP_MAXITER | |
| 997 | ✗ | explicit_h_carry = h_entry | |
| 998 | ✗ | return | |
| 999 | end if | ||
| 1000 | 11809 | nlast = solution%nsteps + 1 | |
| 1001 |
2/2✓ Branch 0 taken 47236 times.
✓ Branch 1 taken 11809 times.
|
59045 | si%z(1:4) = solution%y(:, nlast) |
| 1002 | ! solution%h is trimmed to nsteps entries, one FEWER than solution%t and | ||
| 1003 | ! solution%y, because a step size belongs to the interval between two | ||
| 1004 | ! recorded points rather than to a point. ode_integrate_radau and | ||
| 1005 | ! ode_integrate_gbs allocate nsteps+1 for both, so the same expression is | ||
| 1006 | ! in bounds there; here nlast would read one past the end and carry | ||
| 1007 | ! whatever that memory held into the next macro-step's initial step size. | ||
| 1008 |
1/2✓ Branch 0 taken 11809 times.
✗ Branch 1 not taken.
|
11809 | if (solution%nsteps > 0) explicit_h_carry = solution%h(solution%nsteps) |
| 1009 | 11809 | call explicit_step_finish(h_entry, ierr) | |
| 1010 |
2/2✓ Branch 0 taken 11800 times.
✓ Branch 1 taken 9 times.
|
11809 | if (ierr /= SYMPLECTIC_STEP_OK) return |
| 1011 | 11800 | ktau = ktau + 1 | |
| 1012 | end do | ||
| 1013 | 11800 | call sync_field_to_state | |
| 1014 |
26/52✗ Branch 0 not taken.
✓ Branch 1 taken 11809 times.
✓ Branch 2 taken 11809 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 11809 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 11809 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 11809 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 11809 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 11809 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 11809 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 11809 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 11809 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 11809 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 11809 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 11809 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 11809 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 11809 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 11809 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 11809 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 11809 times.
✗ Branch 36 not taken.
✓ Branch 37 taken 11809 times.
✗ Branch 38 not taken.
✓ Branch 39 taken 11809 times.
✓ Branch 40 taken 11809 times.
✗ Branch 41 not taken.
✓ Branch 42 taken 11809 times.
✗ Branch 43 not taken.
✓ Branch 44 taken 11809 times.
✗ Branch 45 not taken.
✓ Branch 46 taken 11809 times.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✓ Branch 49 taken 11809 times.
✓ Branch 50 taken 11809 times.
✗ Branch 51 not taken.
|
23618 | end subroutine orbit_timestep_cashkarp45 |
| 1015 | |||
| 1016 | end module orbit_symplectic_quasi | ||
| 1017 |