| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module orbit_symplectic_base | ||
| 2 | use field_can_mod, only: eval_field => evaluate, field_can_t, get_val, get_derivatives, & | ||
| 3 | get_derivatives2 | ||
| 4 | |||
| 5 | implicit none | ||
| 6 | |||
| 7 | ! Define real(dp) kind parameter | ||
| 8 | integer, parameter :: dp = kind(1.0d0) | ||
| 9 | |||
| 10 | logical, parameter :: extrap_field = .True. ! do extrapolation after final iteration | ||
| 11 | |||
| 12 | ! Upper radial bound for Newton trials: normalized-flux edge s = 1 by | ||
| 13 | ! default; SPECTRE uses Mvol + 1, with physical loss localized at Mvol. | ||
| 14 | real(dp) :: sympl_rmax = 1.0d0 | ||
| 15 | !$acc declare copyin(sympl_rmax) | ||
| 16 | |||
| 17 | ! Integration methods | ||
| 18 | integer, parameter :: RK45 = 0, EXPL_IMPL_EULER = 1, IMPL_EXPL_EULER = 2, & | ||
| 19 | MIDPOINT = 3, GAUSS1 = 4, GAUSS2 = 5, GAUSS3 = 6, GAUSS4 = 7, & | ||
| 20 | CASH_KARP = 8, DORMAND_PRINCE = 9, LOBATTO3 = 15 | ||
| 21 | |||
| 22 | ! High-order explicit methods on the 4D canonical chart, for benchmarking | ||
| 23 | ! against the symplectic schemes above. They integrate the same canonical | ||
| 24 | ! Hamiltonian through the same f_ode right-hand side, so a comparison per | ||
| 25 | ! field evaluation is like for like -- unlike integmode = 0, which runs the | ||
| 26 | ! 5D drift-kinetic form with a different right-hand side and coordinates. | ||
| 27 | ! | ||
| 28 | ! RADAU15 is the first-order-system formulation of IAS15 (Rein & Spiegel | ||
| 29 | ! 2015), the strongest published claim that a non-symplectic method beats | ||
| 30 | ! symplectic ones on long integrations. GBS16 is Gragg-Bulirsch-Stoer | ||
| 31 | ! extrapolation, the classical high-accuracy explicit method of celestial | ||
| 32 | ! mechanics. | ||
| 33 | ! TDRK24 is the Runge-Kutta-Nystrom transplant: an RKN tableau applied to | ||
| 34 | ! the first-order guiding-centre system through zdd = F'(z)F(z), which is a | ||
| 35 | ! special two-derivative Runge-Kutta method (Chan & Tsai 2010). Fixed step, | ||
| 36 | ! so npoiper2 sets its resolution as it does for the symplectic schemes; | ||
| 37 | ! TDRK24A is the same method under tolerance control, which is what makes it | ||
| 38 | ! comparable to the adaptive methods rather than to the symplectic ones. | ||
| 39 | ! | ||
| 40 | ! CASHKARP45 is the classical adaptive workhorse on this same chart. | ||
| 41 | ! integmode = 0 also runs a Cash-Karp, but on the 5D drift-kinetic form with | ||
| 42 | ! a different right-hand side and coordinates, so its accuracy and cost are | ||
| 43 | ! not commensurable with anything here. CASH_KARP = 8 and DORMAND_PRINCE = 9 | ||
| 44 | ! above are a different thing again: they select the RK5(4) kernels of the | ||
| 45 | ! GPU tracing paths (simple_gpu, simple_cuda_native). CASHKARP45 is the CPU | ||
| 46 | ! method on the canonical chart and does not dispatch to them. | ||
| 47 | integer, parameter :: RADAU15 = 20, GBS16 = 21, CASHKARP45 = 22, & | ||
| 48 | TDRK24 = 24, TDRK24A = 25 | ||
| 49 | integer, parameter :: SYMPLECTIC_STEP_OK = 0 | ||
| 50 | integer, parameter :: SYMPLECTIC_STEP_OUTSIDE_DOMAIN = 1 | ||
| 51 | integer, parameter :: SYMPLECTIC_STEP_MAXITER = 2 | ||
| 52 | integer, parameter :: SYMPLECTIC_STEP_LINEAR_SOLVE = 3 | ||
| 53 | integer, parameter :: SYMPLECTIC_STEP_BOUNDARY = 4 | ||
| 54 | integer, parameter :: SYMPLECTIC_STEP_EVENT_NOT_CONVERGED = 5 | ||
| 55 | integer, parameter :: SYMPLECTIC_STEP_BOUNDARY_LIMITED = 6 | ||
| 56 | logical :: symplectic_newton_warning_mode = .true. | ||
| 57 | real(dp), parameter :: symplectic_newton_warning_factor = 100.0_dp | ||
| 58 | real(dp), parameter :: symplectic_spectre_newton_warning_factor = 10.0_dp | ||
| 59 | real(dp) :: boundary_event_fraction_tolerance = -1d0 | ||
| 60 | real(dp) :: boundary_event_radial_tolerance = -1d0 | ||
| 61 | |||
| 62 | type :: symplectic_integrator_t | ||
| 63 | real(dp) :: atol | ||
| 64 | real(dp) :: rtol | ||
| 65 | real(dp) :: warning_factor = symplectic_newton_warning_factor | ||
| 66 | |||
| 67 | ! Current phase-space coordinates z and old pth | ||
| 68 | real(dp), dimension(4) :: z ! z = (r, th, ph, pphi) | ||
| 69 | real(dp) :: pthold | ||
| 70 | |||
| 71 | ! Timestep and variables from z0 | ||
| 72 | integer :: ntau | ||
| 73 | real(dp) :: dt | ||
| 74 | real(dp) :: pabs | ||
| 75 | real(dp) :: last_step_fraction = 1d0 | ||
| 76 | real(dp) :: last_event_radial_residual = 0d0 | ||
| 77 | real(dp) :: last_event_fraction_width = 0d0 | ||
| 78 | end type symplectic_integrator_t | ||
| 79 | |||
| 80 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 81 | ! | ||
| 82 | ! Composition method with 2s internal stages according to Hairer, 2002 V.3.1 | ||
| 83 | ! | ||
| 84 | integer, parameter :: S_MAX = 32 | ||
| 85 | type :: multistage_integrator_t | ||
| 86 | integer :: s | ||
| 87 | real(dp) :: alpha(S_MAX), beta(S_MAX) | ||
| 88 | type(symplectic_integrator_t) stages(2*S_MAX) | ||
| 89 | end type multistage_integrator_t | ||
| 90 | |||
| 91 | abstract interface | ||
| 92 | subroutine orbit_timestep_sympl_i(si, f, ierr) | ||
| 93 | import :: symplectic_integrator_t, field_can_t | ||
| 94 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 95 | type(field_can_t), intent(inout) :: f | ||
| 96 | integer, intent(out) :: ierr | ||
| 97 | end subroutine orbit_timestep_sympl_i | ||
| 98 | end interface | ||
| 99 | |||
| 100 | abstract interface | ||
| 101 | subroutine orbit_timestep_quasi_i(ierr) | ||
| 102 | integer, intent(out) :: ierr | ||
| 103 | end subroutine orbit_timestep_quasi_i | ||
| 104 | end interface | ||
| 105 | |||
| 106 | contains | ||
| 107 | |||
| 108 | 20528 | subroutine coeff_rk_gauss(n, a, b, c) | |
| 109 | !$acc routine seq | ||
| 110 | integer, intent(in) :: n | ||
| 111 | real(dp), intent(inout) :: a(n, n), b(n), c(n) | ||
| 112 | |||
| 113 |
2/2✓ Branch 0 taken 3699 times.
✓ Branch 1 taken 16829 times.
|
20528 | if (n == 1) then |
| 114 | 3699 | a(1, 1) = 0.5d0 | |
| 115 | 3699 | b(1) = 1.0d0 | |
| 116 | 3699 | c(1) = 0.5d0 | |
| 117 |
2/2✓ Branch 0 taken 5979 times.
✓ Branch 1 taken 10850 times.
|
16829 | elseif (n == 2) then |
| 118 | 5979 | a(1, 1) = 0.25d0 | |
| 119 | 5979 | a(1, 2) = -0.038675134594812d0 | |
| 120 | 5979 | a(2, 1) = 0.538675134594812d0 | |
| 121 | 5979 | a(2, 2) = 0.25d0 | |
| 122 | |||
| 123 | 5979 | b(1) = 0.5d0 | |
| 124 | 5979 | b(2) = 0.5d0 | |
| 125 | |||
| 126 | 5979 | c(1) = 0.211324865405187d0 | |
| 127 | 5979 | c(2) = 0.788675134594812d0 | |
| 128 |
2/2✓ Branch 0 taken 8259 times.
✓ Branch 1 taken 2591 times.
|
10850 | elseif (n == 3) then |
| 129 | 8259 | a(1, 1) = 0.1388888888888889d0 | |
| 130 | 8259 | a(1, 2) = -0.03597666752493894d0 | |
| 131 | 8259 | a(1, 3) = 0.009789444015308318d0 | |
| 132 | 8259 | a(2, 1) = 0.3002631949808646d0 | |
| 133 | 8259 | a(2, 2) = 0.2222222222222222d0 | |
| 134 | 8259 | a(2, 3) = -0.022485417203086805d0 | |
| 135 | 8259 | a(3, 1) = 0.26798833376246944d0 | |
| 136 | 8259 | a(3, 2) = 0.48042111196938336d0 | |
| 137 | 8259 | a(3, 3) = 0.1388888888888889d0 | |
| 138 | |||
| 139 | 8259 | b(1) = 0.2777777777777778d0 | |
| 140 | 8259 | b(2) = 0.4444444444444444d0 | |
| 141 | 8259 | b(3) = 0.2777777777777778d0 | |
| 142 | |||
| 143 | 8259 | c(1) = 0.1127016653792583d0 | |
| 144 | 8259 | c(2) = 0.5d0 | |
| 145 | 8259 | c(3) = 0.8872983346207417d0 | |
| 146 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2590 times.
|
2591 | elseif (n == 4) then ! with help of coefficients from GeometricIntegrators.jl of Michael Kraus |
| 147 | 2590 | a(1, 1) = 0.086963711284363462428182d0 | |
| 148 | 2590 | a(1, 2) = -0.026604180084998794303397d0 | |
| 149 | 2590 | a(1, 3) = 0.012627462689404725035280d0 | |
| 150 | 2590 | a(1, 4) = -0.003555149685795683332096d0 | |
| 151 | |||
| 152 | 2590 | a(2, 1) = 0.188118117499868064967927d0 | |
| 153 | 2590 | a(2, 2) = 0.163036288715636523694030d0 | |
| 154 | 2590 | a(2, 3) = -0.027880428602470894855481d0 | |
| 155 | 2590 | a(2, 4) = 0.006735500594538155853808d0 | |
| 156 | |||
| 157 | 2590 | a(3, 1) = 0.167191921974188778543535d0 | |
| 158 | 2590 | a(3, 2) = 0.353953006033743966529670d0 | |
| 159 | 2590 | a(3, 3) = 0.163036288715636523694030d0 | |
| 160 | 2590 | a(3, 4) = -0.014190694931141143581010d0 | |
| 161 | |||
| 162 | 2590 | a(4, 1) = 0.177482572254522602550608d0 | |
| 163 | 2590 | a(4, 2) = 0.313445114741868369190314d0 | |
| 164 | 2590 | a(4, 3) = 0.352676757516271865977586d0 | |
| 165 | 2590 | a(4, 4) = 0.086963711284363462428182d0 | |
| 166 | |||
| 167 | 2590 | b(1) = 0.173927422568726924856364d0 | |
| 168 | 2590 | b(2) = 0.326072577431273047388061d0 | |
| 169 | 2590 | b(3) = 0.326072577431273047388061d0 | |
| 170 | 2590 | b(4) = 0.173927422568726924856364d0 | |
| 171 | |||
| 172 | 2590 | c(1) = 0.069431844202973713731097d0 | |
| 173 | 2590 | c(2) = 0.330009478207571871344328d0 | |
| 174 | 2590 | c(3) = 0.669990521792428128655672d0 | |
| 175 | 2590 | c(4) = 0.930568155797026341780054d0 | |
| 176 | else | ||
| 177 | ! not implemented | ||
| 178 |
4/4✓ Branch 0 taken 1 time.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 25 times.
✓ Branch 3 taken 5 times.
|
31 | a = 0d0 |
| 179 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 5 times.
|
6 | b = 0d0 |
| 180 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1 time.
|
6 | c = 0d0 |
| 181 | end if | ||
| 182 | 20528 | end subroutine coeff_rk_gauss | |
| 183 | |||
| 184 | 5881 | subroutine coeff_rk_lobatto(n, a, ahat, b, c) | |
| 185 | integer, intent(in) :: n | ||
| 186 | real(dp), intent(inout) :: a(n, n), ahat(n, n), b(n), c(n) | ||
| 187 | |||
| 188 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 5880 times.
|
5881 | if (n == 3) then |
| 189 | 5880 | a(1, 1) = 0d0 | |
| 190 | 5880 | a(1, 2) = 0d0 | |
| 191 | 5880 | a(1, 3) = 0d0 | |
| 192 | |||
| 193 | 5880 | a(2, 1) = 0.20833333333333334d0 | |
| 194 | 5880 | a(2, 2) = 0.33333333333333333d0 | |
| 195 | 5880 | a(2, 3) = -0.041666666666666664d0 | |
| 196 | |||
| 197 | 5880 | a(3, 1) = 0.16666666666666667d0 | |
| 198 | 5880 | a(3, 2) = 0.66666666666666667d0 | |
| 199 | 5880 | a(3, 3) = 0.16666666666666667d0 | |
| 200 | |||
| 201 | 5880 | ahat(1, 1) = 0.16666666666666667d0 | |
| 202 | 5880 | ahat(1, 2) = -0.16666666666666667d0 | |
| 203 | 5880 | ahat(1, 3) = 0d0 | |
| 204 | |||
| 205 | 5880 | ahat(2, 1) = 0.16666666666666667d0 | |
| 206 | 5880 | ahat(2, 2) = 0.33333333333333333d0 | |
| 207 | 5880 | ahat(2, 3) = 0d0 | |
| 208 | |||
| 209 | 5880 | ahat(3, 1) = 0.16666666666666667d0 | |
| 210 | 5880 | ahat(3, 2) = 0.83333333333333333d0 | |
| 211 | 5880 | ahat(3, 3) = 0d0 | |
| 212 | |||
| 213 | 5880 | b(1) = 0.16666666666666667d0 | |
| 214 | 5880 | b(2) = 0.66666666666666667d0 | |
| 215 | 5880 | b(3) = 0.16666666666666667d0 | |
| 216 | |||
| 217 | 5880 | c(1) = 0d0 | |
| 218 | 5880 | c(2) = 0.5d0 | |
| 219 | 5880 | c(3) = 1.0d0 | |
| 220 | |||
| 221 | else | ||
| 222 | ! not implemented | ||
| 223 |
4/4✓ Branch 0 taken 1 time.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 16 times.
✓ Branch 3 taken 4 times.
|
21 | a = 0d0 |
| 224 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 4 times.
|
5 | b = 0d0 |
| 225 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 time.
|
5 | c = 0d0 |
| 226 | end if | ||
| 227 | 5881 | end subroutine coeff_rk_lobatto | |
| 228 | |||
| 229 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 230 | ! | ||
| 231 | ! Lobatto (IIIA)-(IIIB) Runge-Kutta method with s internal stages (n=4*s variables) | ||
| 232 | ! | ||
| 233 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2858 times.
|
2858 | subroutine f_rk_lobatto(si, fs, s, x, fvec, jactype) |
| 234 | ! | ||
| 235 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 236 | integer, intent(in) :: s | ||
| 237 | type(field_can_t), intent(inout) :: fs(:) | ||
| 238 | real(dp), intent(in) :: x(4*s) ! = (rend, thend, phend, pphend) | ||
| 239 | real(dp), intent(out) :: fvec(4*s) | ||
| 240 | integer, intent(in) :: jactype ! 0 = no second derivatives, 2 = second derivatives | ||
| 241 | |||
| 242 | 5716 | real(dp) :: a(s, s), ahat(s, s), b(s), c(s), Hprime(s) | |
| 243 | integer :: k, l ! counters | ||
| 244 | |||
| 245 | 2858 | call coeff_rk_lobatto(s, a, ahat, b, c) | |
| 246 | |||
| 247 | 2858 | call eval_field(fs(1), x(1), si%z(2), si%z(3), jactype) | |
| 248 | 2858 | call get_derivatives(fs(1), x(2)) | |
| 249 | |||
| 250 |
2/2✓ Branch 1 taken 2858 times.
✓ Branch 2 taken 5716 times.
|
8574 | do k = 2, s |
| 251 | 5716 | call eval_field(fs(k), x(4*k - 3 - 2), x(4*k - 2 - 2), x(4*k - 1 - 2), jactype) | |
| 252 | 8574 | call get_derivatives(fs(k), x(4*k - 2)) | |
| 253 | end do | ||
| 254 | |||
| 255 |
2/2✓ Branch 0 taken 8574 times.
✓ Branch 1 taken 2858 times.
|
11432 | Hprime = fs%dH(1)/fs%dpth(1) |
| 256 | |||
| 257 | 2858 | fvec(1) = fs(1)%pth - si%pthold | |
| 258 | 2858 | fvec(2) = x(2) - si%z(4) | |
| 259 | |||
| 260 |
2/2✓ Branch 0 taken 2858 times.
✓ Branch 1 taken 8574 times.
|
11432 | do l = 1, s |
| 261 | 8574 | fvec(1) = fvec(1) + si%dt*ahat(1, l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) ! pthdot | |
| 262 | 11432 | fvec(2) = fvec(2) + si%dt*ahat(1, l)*(fs(l)%dH(3) - Hprime(l)*fs(l)%dpth(3)) ! pphdot | |
| 263 | end do | ||
| 264 | |||
| 265 |
2/2✓ Branch 0 taken 2858 times.
✓ Branch 1 taken 5716 times.
|
8574 | do k = 2, s |
| 266 | 5716 | fvec(4*k - 3 - 2) = fs(k)%pth - si%pthold | |
| 267 | 5716 | fvec(4*k - 2 - 2) = x(4*k - 2 - 2) - si%z(2) | |
| 268 | 5716 | fvec(4*k - 1 - 2) = x(4*k - 1 - 2) - si%z(3) | |
| 269 | 8574 | fvec(4*k - 2) = x(4*k - 2) - si%z(4) | |
| 270 | end do | ||
| 271 | |||
| 272 |
2/2✓ Branch 0 taken 8574 times.
✓ Branch 1 taken 2858 times.
|
11432 | do l = 1, s |
| 273 |
2/2✓ Branch 0 taken 17148 times.
✓ Branch 1 taken 8574 times.
|
28580 | do k = 2, s |
| 274 | 17148 | fvec(4*k-3-2) = fvec(4*k-3-2) + si%dt*ahat(k,l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) ! pthdot | |
| 275 | 17148 | fvec(4*k-2-2) = fvec(4*k-2-2) - si%dt*a(k,l)*Hprime(l) ! thdot | |
| 276 | 17148 | fvec(4*k-1-2) = fvec(4*k-1-2) - si%dt*a(k,l)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph ! phdot | |
| 277 | 25722 | fvec(4*k-2) = fvec(4*k-2) + si%dt*ahat(k,l)*(fs(l)%dH(3) - Hprime(l)*fs(l)%dpth(3)) ! pphdot | |
| 278 | end do | ||
| 279 | end do | ||
| 280 | |||
| 281 | 2858 | end subroutine f_rk_lobatto | |
| 282 | |||
| 283 | ✗ | end module orbit_symplectic_base | |
| 284 |