| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module orbit_symplectic | ||
| 2 | |||
| 3 | use, intrinsic :: iso_fortran_env, only: dp => real64, int64 | ||
| 4 | use util, only: pi, twopi | ||
| 5 | use field_can_mod, only: field_can_t, get_val, get_derivatives, get_derivatives2, & | ||
| 6 | eval_field => evaluate | ||
| 7 | use orbit_symplectic_base, only: symplectic_integrator_t, multistage_integrator_t, & | ||
| 8 | RK45, EXPL_IMPL_EULER, IMPL_EXPL_EULER, MIDPOINT, GAUSS1, GAUSS2, GAUSS3, GAUSS4, & | ||
| 9 | LOBATTO3, RADAU15, GBS16, CASHKARP45, TDRK24, TDRK24A, S_MAX, orbit_timestep_sympl_i, extrap_field, sympl_rmax, & | ||
| 10 | coeff_rk_gauss, coeff_rk_lobatto, f_rk_lobatto, & | ||
| 11 | SYMPLECTIC_STEP_OK, SYMPLECTIC_STEP_OUTSIDE_DOMAIN, & | ||
| 12 | SYMPLECTIC_STEP_MAXITER, SYMPLECTIC_STEP_LINEAR_SOLVE, & | ||
| 13 | SYMPLECTIC_STEP_BOUNDARY, SYMPLECTIC_STEP_EVENT_NOT_CONVERGED, & | ||
| 14 | SYMPLECTIC_STEP_BOUNDARY_LIMITED, boundary_event_fraction_tolerance, & | ||
| 15 | boundary_event_radial_tolerance, symplectic_newton_warning_mode | ||
| 16 | use orbit_symplectic_quasi, only: orbit_timestep_quasi, timestep_expl_impl_euler_quasi, & | ||
| 17 | timestep_impl_expl_euler_quasi, timestep_midpoint_quasi, orbit_timestep_rk45, & | ||
| 18 | orbit_timestep_radau15, orbit_timestep_gbs16, orbit_timestep_tdrk24, & | ||
| 19 | orbit_timestep_tdrk24a, orbit_timestep_cashkarp45, & | ||
| 20 | reset_explicit_step_carry, & | ||
| 21 | si_quasi => si, f_quasi => f, & | ||
| 22 | timestep_rk_gauss_quasi, timestep_rk_lobatto_quasi | ||
| 23 | use orbit_symplectic_euler1, only: sympl_euler1_residual, sympl_euler1_jacobian, & | ||
| 24 | sympl_euler1_newton_iter, sympl_euler1_extrapolate_field, & | ||
| 25 | sympl_euler1_advance_angles | ||
| 26 | use orbit_rk_core, only: rk_gauss_residual, rk_gauss_jacobian, rk_solve, MODEL_GC | ||
| 27 | use vector_potentail_mod, only: torflux | ||
| 28 | use diag_counters, only: count_event, EVT_NEWTON1_MAXIT, EVT_NEWTON2_MAXIT, & | ||
| 29 | EVT_RK_GAUSS_MAXIT, EVT_RK_LOBATTO_MAXIT, EVT_FIXPOINT_MAXIT, EVT_R_NEGATIVE, & | ||
| 30 | EVT_MIDPOINT_MAXIT, EVT_WARNING_MAXIT_ACCEPT, & | ||
| 31 | EVT_WARNING_MAXIT_REJECT_100, EVT_WARNING_MAXIT_REJECT_1E4, & | ||
| 32 | EVT_WARNING_MAXIT_REJECT_LARGE, EVT_WARNING_MAXIT_NONFINITE, & | ||
| 33 | EVT_RETRY_EXHAUSTED, EVT_WARNING_AXIS_CROSSING_ACCEPT | ||
| 34 | |||
| 35 | implicit none | ||
| 36 | |||
| 37 | procedure(orbit_timestep_sympl_i), pointer :: orbit_timestep_sympl => null() | ||
| 38 | procedure(orbit_timestep_sympl_i), pointer, private :: raw_timestep_sympl => null() | ||
| 39 | |||
| 40 | contains | ||
| 41 | |||
| 42 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 563 times.
|
563 | pure logical function finite_newton_iterate(x) |
| 43 | real(dp), intent(in) :: x(:) | ||
| 44 | |||
| 45 | integer(int64), parameter :: exponent_mask = & | ||
| 46 | int(z'7FF0000000000000', int64) | ||
| 47 | integer(int64) :: bits | ||
| 48 | integer :: i | ||
| 49 | |||
| 50 | 563 | finite_newton_iterate = .false. | |
| 51 |
2/2✓ Branch 0 taken 635 times.
✓ Branch 1 taken 559 times.
|
1194 | do i = 1, size(x) |
| 52 | 635 | bits = transfer(x(i), bits) | |
| 53 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 631 times.
|
1194 | if (iand(bits, exponent_mask) == exponent_mask) return |
| 54 | end do | ||
| 55 | 559 | finite_newton_iterate = .true. | |
| 56 | 559 | end function finite_newton_iterate | |
| 57 | |||
| 58 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 17 times.
|
17 | logical function accept_warning_maxiter(x, xlast, tolref, rtol, warning_factor) |
| 59 | ! Continue only when the final Newton correction is finite and inside this | ||
| 60 | ! integrator's tolerance-scaled guard. Other iterates enter retry recovery. | ||
| 61 | real(dp), intent(in) :: x(:), xlast(:), tolref(:), rtol, warning_factor | ||
| 62 | real(dp) :: correction_ratio | ||
| 63 | |||
| 64 | 17 | accept_warning_maxiter = .false. | |
| 65 |
3/4✓ Branch 0 taken 7 times.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
|
17 | if (.not. symplectic_newton_warning_mode .or. rtol <= 0.0_dp) return |
| 66 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 4 times.
|
7 | if (.not. finite_newton_iterate(x)) then |
| 67 | 3 | call count_event(EVT_WARNING_MAXIT_NONFINITE) | |
| 68 | 3 | return | |
| 69 | end if | ||
| 70 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (.not. finite_newton_iterate(xlast)) then |
| 71 | ✗ | call count_event(EVT_WARNING_MAXIT_NONFINITE) | |
| 72 | ✗ | return | |
| 73 | end if | ||
| 74 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (.not. finite_newton_iterate(tolref)) then |
| 75 | ✗ | call count_event(EVT_WARNING_MAXIT_NONFINITE) | |
| 76 | ✗ | return | |
| 77 | end if | ||
| 78 | correction_ratio = maxval(abs(x - xlast)/ & | ||
| 79 |
6/10✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 8 times.
✓ Branch 7 taken 4 times.
✓ Branch 8 taken 4 times.
✓ Branch 9 taken 4 times.
|
16 | max(rtol*abs(tolref), tiny(1.0_dp))) |
| 80 | 4 | accept_warning_maxiter = correction_ratio <= warning_factor | |
| 81 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
|
4 | if (accept_warning_maxiter) then |
| 82 | 2 | call count_event(EVT_WARNING_MAXIT_ACCEPT) | |
| 83 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | else if (correction_ratio <= 100.0_dp) then |
| 84 | 1 | call count_event(EVT_WARNING_MAXIT_REJECT_100) | |
| 85 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | else if (correction_ratio <= 1.0e4_dp) then |
| 86 | ✗ | call count_event(EVT_WARNING_MAXIT_REJECT_1E4) | |
| 87 | else | ||
| 88 | 1 | call count_event(EVT_WARNING_MAXIT_REJECT_LARGE) | |
| 89 | end if | ||
| 90 | end function accept_warning_maxiter | ||
| 91 | |||
| 92 | 10 | logical function accept_warning_axis_crossing(si, x, xlast, tolref, rtol) | |
| 93 | ! Near the polar-coordinate axis, the endpoint and midpoint radii can change | ||
| 94 | ! sign while the physical step and its angular/momentum components have | ||
| 95 | ! converged. Accept only that local chart crossing; the caller immediately | ||
| 96 | ! maps the negative endpoint to the equivalent positive-radius chart. | ||
| 97 | type(symplectic_integrator_t), intent(in) :: si | ||
| 98 | real(dp), intent(in) :: x(5), xlast(5), tolref(5), rtol | ||
| 99 | real(dp), parameter :: axis_radius = 1.0e-6_dp | ||
| 100 | |||
| 101 | 10 | accept_warning_axis_crossing = .false. | |
| 102 |
3/4✓ Branch 0 taken 7 times.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
|
10 | if (.not. symplectic_newton_warning_mode .or. rtol <= 0.0_dp) return |
| 103 |
3/4✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 6 times.
|
7 | if (si%z(1) < 0.0_dp .or. si%z(1) > axis_radius) return |
| 104 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 5 times.
|
6 | if (.not. finite_newton_iterate(x)) return |
| 105 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if (.not. finite_newton_iterate(xlast)) return |
| 106 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if (.not. finite_newton_iterate(tolref)) return |
| 107 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 4 times.
|
5 | if (x(1) >= 0.0_dp) return |
| 108 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (max(abs(x(1)), abs(x(5))) > axis_radius) return |
| 109 | accept_warning_axis_crossing = all(abs(x(2:4) - xlast(2:4)) <= & | ||
| 110 |
4/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 1 time.
|
13 | si%warning_factor*rtol*abs(tolref(2:4))) |
| 111 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
|
4 | if (accept_warning_axis_crossing) & |
| 112 | 3 | call count_event(EVT_WARNING_AXIS_CROSSING_ACCEPT) | |
| 113 | end function accept_warning_axis_crossing | ||
| 114 | |||
| 115 | 3962177 | subroutine advance_symplectic_with_retry(si, f, stepper, status, & | |
| 116 | accepted_fraction) | ||
| 117 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 118 | type(field_can_t), intent(inout) :: f | ||
| 119 | procedure(orbit_timestep_sympl_i) :: stepper | ||
| 120 | integer, intent(out) :: status | ||
| 121 | real(dp), intent(out), optional :: accepted_fraction | ||
| 122 | real(dp) :: accepted_fraction_local | ||
| 123 | |||
| 124 | call advance_retry_interval(si, f, stepper, si%dt, 0, status, & | ||
| 125 | 3962177 | accepted_fraction_local) | |
| 126 |
2/2✓ Branch 0 taken 17 times.
✓ Branch 1 taken 3962160 times.
|
3962177 | if (present(accepted_fraction)) accepted_fraction = accepted_fraction_local |
| 127 | 3962177 | end subroutine advance_symplectic_with_retry | |
| 128 | |||
| 129 | 3962190 | recursive subroutine advance_retry_interval(si, f, stepper, interval, depth, & | |
| 130 | status, accepted_fraction) | ||
| 131 | integer, parameter :: max_retry_depth = 8 | ||
| 132 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 133 | type(field_can_t), intent(inout) :: f | ||
| 134 | procedure(orbit_timestep_sympl_i) :: stepper | ||
| 135 | real(dp), intent(in) :: interval | ||
| 136 | integer, intent(in) :: depth | ||
| 137 | integer, intent(out) :: status | ||
| 138 | real(dp), intent(out) :: accepted_fraction | ||
| 139 | type(symplectic_integrator_t) :: initial_integrator | ||
| 140 | type(field_can_t) :: initial_field | ||
| 141 | real(dp) :: child_fraction | ||
| 142 | logical :: recoverable | ||
| 143 | |||
| 144 | 3962190 | accepted_fraction = 0.0_dp | |
| 145 | 3962190 | initial_integrator = si | |
| 146 | 3962190 | initial_field = f | |
| 147 | 3962190 | si%dt = interval | |
| 148 | 3962190 | call stepper(si, f, status) | |
| 149 |
2/2✓ Branch 0 taken 3962168 times.
✓ Branch 1 taken 22 times.
|
3962190 | if (status == SYMPLECTIC_STEP_OK .or. & |
| 150 | status == SYMPLECTIC_STEP_BOUNDARY) then | ||
| 151 |
2/2✓ Branch 0 taken 3962101 times.
✓ Branch 1 taken 67 times.
|
3962168 | if (status == SYMPLECTIC_STEP_OK) then |
| 152 | 3962101 | accepted_fraction = 1.0_dp | |
| 153 | else | ||
| 154 | 67 | accepted_fraction = si%last_step_fraction | |
| 155 | end if | ||
| 156 | 3962168 | si%dt = interval | |
| 157 | 3962182 | return | |
| 158 | end if | ||
| 159 | |||
| 160 | 22 | si = initial_integrator | |
| 161 | 22 | f = initial_field | |
| 162 | 22 | si%dt = interval | |
| 163 | recoverable = status == SYMPLECTIC_STEP_MAXITER .or. & | ||
| 164 | 22 | status == SYMPLECTIC_STEP_LINEAR_SOLVE | |
| 165 |
6/6✓ Branch 0 taken 19 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 11 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 10 times.
|
22 | if (.not. symplectic_newton_warning_mode .or. .not. recoverable .or. & |
| 166 | depth >= max_retry_depth) then | ||
| 167 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 1 time.
|
4 | if (recoverable .and. depth >= max_retry_depth) & |
| 168 | 1 | call count_event(EVT_RETRY_EXHAUSTED) | |
| 169 | 12 | return | |
| 170 | end if | ||
| 171 | |||
| 172 | ! For a fixed subdivision, the two accepted maps remain symplectic and cover | ||
| 173 | ! the original interval exactly. State-dependent subdivision is a warning-mode | ||
| 174 | ! recovery path, not a globally fixed symplectic map. | ||
| 175 | call advance_retry_interval(si, f, stepper, 0.5_dp*interval, depth + 1, & | ||
| 176 | 10 | status, child_fraction) | |
| 177 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
10 | if (status == SYMPLECTIC_STEP_BOUNDARY) then |
| 178 | ✗ | si%last_step_fraction = 0.5_dp*si%last_step_fraction | |
| 179 | ✗ | si%last_event_fraction_width = 0.5_dp*si%last_event_fraction_width | |
| 180 | ✗ | accepted_fraction = 0.5_dp*child_fraction | |
| 181 | ✗ | si%dt = interval | |
| 182 | ✗ | return | |
| 183 | end if | ||
| 184 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 7 times.
|
10 | if (status == SYMPLECTIC_STEP_OK) then |
| 185 | call advance_retry_interval(si, f, stepper, 0.5_dp*interval, depth + 1, & | ||
| 186 | 3 | status, child_fraction) | |
| 187 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2 times.
|
3 | if (status == SYMPLECTIC_STEP_BOUNDARY) then |
| 188 | 1 | si%last_step_fraction = 0.5_dp + 0.5_dp*si%last_step_fraction | |
| 189 | 1 | si%last_event_fraction_width = 0.5_dp*si%last_event_fraction_width | |
| 190 | 1 | accepted_fraction = 0.5_dp + 0.5_dp*child_fraction | |
| 191 | 1 | si%dt = interval | |
| 192 | 1 | return | |
| 193 | end if | ||
| 194 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | if (status == SYMPLECTIC_STEP_OK) then |
| 195 | 1 | accepted_fraction = 0.5_dp + 0.5_dp*child_fraction | |
| 196 | 1 | si%dt = interval | |
| 197 | 1 | return | |
| 198 | end if | ||
| 199 | end if | ||
| 200 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
8 | if (status /= SYMPLECTIC_STEP_OK) then |
| 201 | 8 | si = initial_integrator | |
| 202 | 8 | f = initial_field | |
| 203 | 8 | accepted_fraction = 0.0_dp | |
| 204 | end if | ||
| 205 | 8 | si%dt = interval | |
| 206 | 3962190 | end subroutine advance_retry_interval | |
| 207 | |||
| 208 | 5555121 | subroutine apply_axis_chart_switch(radius, theta, crossed, status) | |
| 209 | real(dp), intent(inout) :: radius, theta | ||
| 210 | logical, intent(out) :: crossed | ||
| 211 | integer, intent(out) :: status | ||
| 212 | |||
| 213 | 5555121 | crossed = radius < 0d0 | |
| 214 | 5555121 | status = SYMPLECTIC_STEP_OK | |
| 215 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 5555117 times.
|
5555121 | if (.not. crossed) return |
| 216 | |||
| 217 | 4 | radius = -radius | |
| 218 | 4 | theta = theta + pi | |
| 219 | 4 | call count_event(EVT_R_NEGATIVE) | |
| 220 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
|
4 | if (radius > sympl_rmax) status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN |
| 221 | end subroutine apply_axis_chart_switch | ||
| 222 | |||
| 223 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 224 | ! | ||
| 225 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 784 times.
|
784 | recursive subroutine orbit_sympl_init(si, f, z, dt, ntau, rtol_init, mode_init) |
| 226 | ! | ||
| 227 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 228 | type(field_can_t), intent(inout) :: f | ||
| 229 | real(dp), intent(in) :: z(:) | ||
| 230 | real(dp), intent(in) :: dt | ||
| 231 | integer, intent(in) :: ntau | ||
| 232 | real(dp), intent(in) :: rtol_init | ||
| 233 | integer, intent(in) :: mode_init | ||
| 234 | |||
| 235 | 784 | si%atol = 1d-15 | |
| 236 | 784 | si%rtol = rtol_init | |
| 237 | 784 | call reset_explicit_step_carry | |
| 238 | |||
| 239 | 784 | si%ntau = ntau | |
| 240 | 784 | si%dt = dt | |
| 241 | 784 | si%last_step_fraction = 1d0 | |
| 242 | 784 | si%last_event_radial_residual = 0d0 | |
| 243 | 784 | si%last_event_fraction_width = 0d0 | |
| 244 | |||
| 245 |
2/2✓ Branch 0 taken 3136 times.
✓ Branch 1 taken 784 times.
|
3920 | si%z = z |
| 246 | |||
| 247 | 784 | call eval_field(f, z(1), z(2), z(3), 0) | |
| 248 | 784 | call get_val(f, si%z(4)) ! for pth | |
| 249 | 784 | si%pthold = f%pth | |
| 250 | |||
| 251 | 784 | select case (mode_init) | |
| 252 | case (RK45) | ||
| 253 | ✗ | nullify(raw_timestep_sympl) | |
| 254 | ✗ | nullify(orbit_timestep_sympl) | |
| 255 | ✗ | orbit_timestep_quasi => orbit_timestep_rk45 | |
| 256 | case (EXPL_IMPL_EULER) | ||
| 257 | 147 | raw_timestep_sympl => orbit_timestep_sympl_expl_impl_euler | |
| 258 | 147 | orbit_timestep_quasi => timestep_expl_impl_euler_quasi | |
| 259 | case (IMPL_EXPL_EULER) | ||
| 260 | 49 | raw_timestep_sympl => orbit_timestep_sympl_impl_expl_euler | |
| 261 | 49 | orbit_timestep_quasi => timestep_impl_expl_euler_quasi | |
| 262 | case (MIDPOINT) | ||
| 263 | 468 | raw_timestep_sympl => orbit_timestep_sympl_midpoint | |
| 264 | 468 | orbit_timestep_quasi => timestep_midpoint_quasi | |
| 265 | case (GAUSS1) | ||
| 266 | 5 | raw_timestep_sympl => orbit_timestep_sympl_gauss1 | |
| 267 | 5 | orbit_timestep_quasi => orbit_timestep_quasi_gauss1 | |
| 268 | case (GAUSS2) | ||
| 269 | 5 | raw_timestep_sympl => orbit_timestep_sympl_gauss2 | |
| 270 | 5 | orbit_timestep_quasi => orbit_timestep_quasi_gauss2 | |
| 271 | case (GAUSS3) | ||
| 272 | 5 | raw_timestep_sympl => orbit_timestep_sympl_gauss3 | |
| 273 | 5 | orbit_timestep_quasi => orbit_timestep_quasi_gauss3 | |
| 274 | case (GAUSS4) | ||
| 275 | 4 | raw_timestep_sympl => orbit_timestep_sympl_gauss4 | |
| 276 | 4 | orbit_timestep_quasi => orbit_timestep_quasi_gauss4 | |
| 277 | case (LOBATTO3) | ||
| 278 | 5 | raw_timestep_sympl => orbit_timestep_sympl_lobatto3 | |
| 279 | 5 | orbit_timestep_quasi => orbit_timestep_quasi_lobatto3 | |
| 280 | case (RADAU15) | ||
| 281 | 19 | raw_timestep_sympl => orbit_timestep_sympl_radau15 | |
| 282 | 19 | orbit_timestep_quasi => orbit_timestep_radau15 | |
| 283 | case (GBS16) | ||
| 284 | 19 | raw_timestep_sympl => orbit_timestep_sympl_gbs16 | |
| 285 | 19 | orbit_timestep_quasi => orbit_timestep_gbs16 | |
| 286 | case (TDRK24) | ||
| 287 | 20 | raw_timestep_sympl => orbit_timestep_sympl_tdrk24 | |
| 288 | 20 | orbit_timestep_quasi => orbit_timestep_tdrk24 | |
| 289 | case (TDRK24A) | ||
| 290 | 19 | raw_timestep_sympl => orbit_timestep_sympl_tdrk24a | |
| 291 | 19 | orbit_timestep_quasi => orbit_timestep_tdrk24a | |
| 292 | case (CASHKARP45) | ||
| 293 | 19 | raw_timestep_sympl => orbit_timestep_sympl_cashkarp45 | |
| 294 | 19 | orbit_timestep_quasi => orbit_timestep_cashkarp45 | |
| 295 | case default | ||
| 296 | ✗ | print *, 'invalid mode for orbit_timestep_sympl: ', mode_init | |
| 297 |
13/15✗ Branch 0 not taken.
✓ Branch 1 taken 147 times.
✓ Branch 2 taken 49 times.
✓ Branch 3 taken 468 times.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 5 times.
✓ Branch 6 taken 5 times.
✓ Branch 7 taken 4 times.
✓ Branch 8 taken 5 times.
✓ Branch 9 taken 19 times.
✓ Branch 10 taken 19 times.
✓ Branch 11 taken 20 times.
✓ Branch 12 taken 19 times.
✓ Branch 13 taken 19 times.
✗ Branch 14 not taken.
|
784 | error stop |
| 298 | end select | ||
| 299 |
3/4✓ Branch 0 taken 784 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 456 times.
✓ Branch 3 taken 328 times.
|
784 | if (mode_init /= RK45 .and. sympl_rmax <= 1d0) then |
| 300 | 456 | orbit_timestep_sympl => orbit_timestep_sympl_with_boundary | |
| 301 |
1/2✓ Branch 0 taken 328 times.
✗ Branch 1 not taken.
|
328 | else if (mode_init /= RK45) then |
| 302 | 328 | orbit_timestep_sympl => raw_timestep_sympl | |
| 303 | end if | ||
| 304 | 784 | end subroutine orbit_sympl_init | |
| 305 | |||
| 306 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 307 | ! | ||
| 308 | ! Bridges the explicit high-order steppers, which live on the quasi path and | ||
| 309 | ! read the module-level state there, into the symplectic stepper signature the | ||
| 310 | ! driver calls. The quasi state is threadprivate, so this is thread safe. | ||
| 311 | ! | ||
| 312 | ! The copies are the price of reusing f_ode unchanged, which is what keeps the | ||
| 313 | ! comparison honest: these methods integrate exactly the same right-hand side | ||
| 314 | ! as the symplectic schemes, through the same field evaluations. | ||
| 315 | 11812 | recursive subroutine orbit_timestep_sympl_radau15(si, f, ierr) | |
| 316 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 317 | type(field_can_t), intent(inout) :: f | ||
| 318 | integer, intent(out) :: ierr | ||
| 319 | |||
| 320 | 11812 | si_quasi = si | |
| 321 | 11812 | f_quasi = f | |
| 322 | 11812 | call orbit_timestep_radau15(ierr) | |
| 323 | 11812 | si = si_quasi | |
| 324 | 11812 | f = f_quasi | |
| 325 | 11812 | end subroutine orbit_timestep_sympl_radau15 | |
| 326 | |||
| 327 | 11806 | recursive subroutine orbit_timestep_sympl_gbs16(si, f, ierr) | |
| 328 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 329 | type(field_can_t), intent(inout) :: f | ||
| 330 | integer, intent(out) :: ierr | ||
| 331 | |||
| 332 | 11806 | si_quasi = si | |
| 333 | 11806 | f_quasi = f | |
| 334 | 11806 | call orbit_timestep_gbs16(ierr) | |
| 335 | 11806 | si = si_quasi | |
| 336 | 11806 | f = f_quasi | |
| 337 | 11806 | end subroutine orbit_timestep_sympl_gbs16 | |
| 338 | |||
| 339 | 22811 | recursive subroutine orbit_timestep_sympl_tdrk24(si, f, ierr) | |
| 340 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 341 | type(field_can_t), intent(inout) :: f | ||
| 342 | integer, intent(out) :: ierr | ||
| 343 | |||
| 344 | 22811 | si_quasi = si | |
| 345 | 22811 | f_quasi = f | |
| 346 | 22811 | call orbit_timestep_tdrk24(ierr) | |
| 347 | 22811 | si = si_quasi | |
| 348 | 22811 | f = f_quasi | |
| 349 | 22811 | end subroutine orbit_timestep_sympl_tdrk24 | |
| 350 | |||
| 351 | 11844 | recursive subroutine orbit_timestep_sympl_tdrk24a(si, f, ierr) | |
| 352 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 353 | type(field_can_t), intent(inout) :: f | ||
| 354 | integer, intent(out) :: ierr | ||
| 355 | |||
| 356 | 11844 | si_quasi = si | |
| 357 | 11844 | f_quasi = f | |
| 358 | 11844 | call orbit_timestep_tdrk24a(ierr) | |
| 359 | 11844 | si = si_quasi | |
| 360 | 11844 | f = f_quasi | |
| 361 | 11844 | end subroutine orbit_timestep_sympl_tdrk24a | |
| 362 | |||
| 363 | 11809 | recursive subroutine orbit_timestep_sympl_cashkarp45(si, f, ierr) | |
| 364 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 365 | type(field_can_t), intent(inout) :: f | ||
| 366 | integer, intent(out) :: ierr | ||
| 367 | |||
| 368 | 11809 | si_quasi = si | |
| 369 | 11809 | f_quasi = f | |
| 370 | 11809 | call orbit_timestep_cashkarp45(ierr) | |
| 371 | 11809 | si = si_quasi | |
| 372 | 11809 | f = f_quasi | |
| 373 | 11809 | end subroutine orbit_timestep_sympl_cashkarp45 | |
| 374 | |||
| 375 | |||
| 376 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 377 | ! | ||
| 378 | 96 | recursive subroutine orbit_timestep_sympl_gauss1(si, f, ierr) | |
| 379 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 380 | type(field_can_t), intent(inout) :: f | ||
| 381 | integer, intent(out) :: ierr | ||
| 382 | |||
| 383 | 96 | call orbit_timestep_sympl_rk_gauss(si, f, 1, ierr) | |
| 384 | 96 | end subroutine orbit_timestep_sympl_gauss1 | |
| 385 | |||
| 386 | 96 | recursive subroutine orbit_timestep_sympl_gauss2(si, f, ierr) | |
| 387 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 388 | type(field_can_t), intent(inout) :: f | ||
| 389 | integer, intent(out) :: ierr | ||
| 390 | |||
| 391 | 96 | call orbit_timestep_sympl_rk_gauss(si, f, 2, ierr) | |
| 392 | 96 | end subroutine orbit_timestep_sympl_gauss2 | |
| 393 | |||
| 394 | 96 | recursive subroutine orbit_timestep_sympl_gauss3(si, f, ierr) | |
| 395 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 396 | type(field_can_t), intent(inout) :: f | ||
| 397 | integer, intent(out) :: ierr | ||
| 398 | |||
| 399 | 96 | call orbit_timestep_sympl_rk_gauss(si, f, 3, ierr) | |
| 400 | 96 | end subroutine orbit_timestep_sympl_gauss3 | |
| 401 | |||
| 402 | 240 | recursive subroutine orbit_timestep_sympl_gauss4(si, f, ierr) | |
| 403 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 404 | type(field_can_t), intent(inout) :: f | ||
| 405 | integer, intent(out) :: ierr | ||
| 406 | |||
| 407 | 240 | call orbit_timestep_sympl_rk_gauss(si, f, 4, ierr) | |
| 408 | 240 | end subroutine orbit_timestep_sympl_gauss4 | |
| 409 | |||
| 410 | 212 | recursive subroutine orbit_timestep_sympl_lobatto3(si, f, ierr) | |
| 411 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 412 | type(field_can_t), intent(inout) :: f | ||
| 413 | integer, intent(out) :: ierr | ||
| 414 | |||
| 415 | 212 | call orbit_timestep_sympl_rk_lobatto(si, f, 3, ierr) | |
| 416 | 212 | end subroutine orbit_timestep_sympl_lobatto3 | |
| 417 | |||
| 418 | 95 | recursive subroutine orbit_timestep_quasi_gauss1(ierr) | |
| 419 | integer, intent(out) :: ierr | ||
| 420 | 95 | call timestep_rk_gauss_quasi(1, ierr) | |
| 421 | 95 | end subroutine orbit_timestep_quasi_gauss1 | |
| 422 | |||
| 423 | 95 | recursive subroutine orbit_timestep_quasi_gauss2(ierr) | |
| 424 | integer, intent(out) :: ierr | ||
| 425 | 95 | call timestep_rk_gauss_quasi(2, ierr) | |
| 426 | 95 | end subroutine orbit_timestep_quasi_gauss2 | |
| 427 | |||
| 428 | 95 | recursive subroutine orbit_timestep_quasi_gauss3(ierr) | |
| 429 | integer, intent(out) :: ierr | ||
| 430 | 95 | call timestep_rk_gauss_quasi(3, ierr) | |
| 431 | 95 | end subroutine orbit_timestep_quasi_gauss3 | |
| 432 | |||
| 433 | ✗ | recursive subroutine orbit_timestep_quasi_gauss4(ierr) | |
| 434 | integer, intent(out) :: ierr | ||
| 435 | ✗ | call timestep_rk_gauss_quasi(4, ierr) | |
| 436 | ✗ | end subroutine orbit_timestep_quasi_gauss4 | |
| 437 | |||
| 438 | ✗ | recursive subroutine orbit_timestep_quasi_lobatto3(ierr) | |
| 439 | integer, intent(out) :: ierr | ||
| 440 | ✗ | call timestep_rk_lobatto_quasi(3, ierr) | |
| 441 | ✗ | end subroutine orbit_timestep_quasi_lobatto3 | |
| 442 | |||
| 443 | |||
| 444 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 445 | ! | ||
| 446 | ✗ | recursive subroutine f_sympl_euler1(si, f, n, x, fvec, iflag) | |
| 447 | ! | ||
| 448 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 449 | type(field_can_t), intent(inout) :: f | ||
| 450 | integer, intent(in) :: n | ||
| 451 | real(dp), intent(in) :: x(n) | ||
| 452 | real(dp), intent(out) :: fvec(n) | ||
| 453 | integer, intent(in) :: iflag | ||
| 454 | |||
| 455 | ✗ | call eval_field(f, x(1), si%z(2), si%z(3), 2) | |
| 456 | ✗ | call get_derivatives2(f, x(2)) | |
| 457 | ✗ | call sympl_euler1_residual(si, f, x, fvec) | |
| 458 | |||
| 459 | ✗ | end subroutine f_sympl_euler1 | |
| 460 | |||
| 461 | |||
| 462 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 463 | ! | ||
| 464 | ✗ | recursive subroutine jac_sympl_euler1(si, f, x, jac) | |
| 465 | ! | ||
| 466 | type(symplectic_integrator_t), intent(in) :: si | ||
| 467 | type(field_can_t), intent(inout) :: f | ||
| 468 | |||
| 469 | real(dp), intent(in) :: x(2) | ||
| 470 | real(dp), intent(out) :: jac(2, 2) | ||
| 471 | |||
| 472 | ✗ | call sympl_euler1_jacobian(si, f, x, jac) | |
| 473 | |||
| 474 | ✗ | end subroutine jac_sympl_euler1 | |
| 475 | |||
| 476 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 477 | ! | ||
| 478 | 17289 | recursive subroutine f_sympl_euler2(si, f, n, x, fvec, iflag) | |
| 479 | ! | ||
| 480 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 481 | type(field_can_t), intent(inout) :: f | ||
| 482 | integer, intent(in) :: n | ||
| 483 | real(dp), intent(in) :: x(n) | ||
| 484 | real(dp), intent(out) :: fvec(n) | ||
| 485 | integer, intent(in) :: iflag | ||
| 486 | |||
| 487 | 17289 | call eval_field(f, x(1), x(2), x(3), 2) | |
| 488 | 17289 | call get_derivatives2(f, si%z(4)) | |
| 489 | |||
| 490 | 17289 | fvec(1) = f%pth - si%pthold | |
| 491 | 17289 | fvec(2) = f%dpth(1)*(x(2) - si%z(2)) - si%dt*f%dH(1) | |
| 492 | 17289 | fvec(3) = f%dpth(1)*f%hph*(x(3) - si%z(3)) - si%dt*(f%dpth(1)*f%vpar - f%dH(1)*f%hth) | |
| 493 | |||
| 494 | 17289 | end subroutine f_sympl_euler2 | |
| 495 | |||
| 496 | |||
| 497 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 498 | ! | ||
| 499 | 17289 | recursive subroutine jac_sympl_euler2(si, f, x, jac) | |
| 500 | ! | ||
| 501 | type(symplectic_integrator_t), intent(in) :: si | ||
| 502 | type(field_can_t), intent(inout) :: f | ||
| 503 | real(dp), intent(in) :: x(3) | ||
| 504 | real(dp), intent(out) :: jac(3, 3) | ||
| 505 | |||
| 506 |
2/2✓ Branch 0 taken 51867 times.
✓ Branch 1 taken 17289 times.
|
69156 | jac(1,:) = f%dpth(1:3) |
| 507 |
2/2✓ Branch 0 taken 51867 times.
✓ Branch 1 taken 17289 times.
|
69156 | jac(2,:) = f%d2pth(1:3)*(x(2) - si%z(2)) - si%dt*f%d2H(1:3) |
| 508 | 17289 | jac(2,2) = jac(2,2) + f%dpth(1) | |
| 509 | jac(3,:) = (f%d2pth(1:3)*f%hph + f%dpth(1)*f%dhph)*(x(3) - si%z(3)) & | ||
| 510 |
2/2✓ Branch 0 taken 51867 times.
✓ Branch 1 taken 17289 times.
|
69156 | - si%dt*(f%d2pth(1:3)*f%vpar + f%dpth(1)*f%dvpar(1:3) - f%d2H(1:3)*f%hth - f%dH(1)*f%dhth) |
| 511 | 17289 | jac(3,3) = jac(3,3) + f%dpth(1)*f%hph | |
| 512 | |||
| 513 | 17289 | end subroutine jac_sympl_euler2 | |
| 514 | |||
| 515 | |||
| 516 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 517 | ! | ||
| 518 | 3320989 | recursive subroutine f_midpoint_part1(si, f, n, x, fvec) | |
| 519 | ! | ||
| 520 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 521 | type(field_can_t), intent(inout) :: f | ||
| 522 | integer, intent(in) :: n | ||
| 523 | real(dp), intent(in) :: x(n) ! = (rend, thend, phend, pphend, rmid) | ||
| 524 | real(dp), intent(out) :: fvec(n) | ||
| 525 | |||
| 526 | ! evaluate at midpoint | ||
| 527 | 3320989 | call eval_field(f, x(5), 0.5d0*(x(2) + si%z(2)), 0.5d0*(x(3) + si%z(3)), 2) | |
| 528 | 3320989 | call get_derivatives2(f, 0.5d0*(x(4) + si%z(4))) | |
| 529 | |||
| 530 | 3320989 | fvec(2) = f%dpth(1)*(x(2) - si%z(2)) - si%dt*f%dH(1) | |
| 531 | 3320989 | fvec(3) = f%dpth(1)*f%hph*(x(3) - si%z(3)) - si%dt*(f%dpth(1)*f%vpar - f%dH(1)*f%hth) | |
| 532 | 3320989 | fvec(4) = f%dpth(1)*(x(4) - si%z(4)) + si%dt*(f%dH(3)*f%dpth(1) - f%dH(1)*f%dpth(3)) | |
| 533 | 3320989 | fvec(5) = f%dpth(1)*(f%pth - si%pthold) + 0.5d0*si%dt*(f%dpth(1)*f%dH(2)-f%dpth(2)*f%dH(1)) | |
| 534 | |||
| 535 | 3320989 | end subroutine f_midpoint_part1 | |
| 536 | |||
| 537 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 538 | ! | ||
| 539 | 3320989 | recursive subroutine f_midpoint_part2(si, f, n, x, fvec) | |
| 540 | ! | ||
| 541 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 542 | type(field_can_t), intent(inout) :: f | ||
| 543 | integer, intent(in) :: n | ||
| 544 | real(dp), intent(in) :: x(n) ! = (rend, thend, phend, pphend, rmid) | ||
| 545 | real(dp), intent(out) :: fvec(n) | ||
| 546 | |||
| 547 | real(dp) :: dpthmid, pthdotbar | ||
| 548 | |||
| 549 | ! save evaluation from midpoint | ||
| 550 | 3320989 | dpthmid = f%dpth(1) | |
| 551 | 3320989 | pthdotbar = f%dpth(1)*f%dH(2) - f%dpth(2)*f%dH(1) | |
| 552 | |||
| 553 | ! evaluate at endpoint | ||
| 554 | 3320989 | call eval_field(f, x(1), x(2), x(3), 2) | |
| 555 | 3320989 | call get_derivatives2(f, x(4)) | |
| 556 | 3320989 | fvec(1) = dpthmid*(f%pth - si%pthold) + si%dt*pthdotbar | |
| 557 | |||
| 558 | 3320989 | end subroutine f_midpoint_part2 | |
| 559 | |||
| 560 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 561 | ! | ||
| 562 | 3320989 | recursive subroutine jac_midpoint_part1(si, f, x, jac) | |
| 563 | ! | ||
| 564 | type(symplectic_integrator_t), intent(in) :: si | ||
| 565 | type(field_can_t), intent(inout) :: f | ||
| 566 | real(dp), intent(in) :: x(5) | ||
| 567 | real(dp), intent(out) :: jac(5, 5) | ||
| 568 | |||
| 569 | 3320989 | jac(2,1) = 0d0 | |
| 570 | 3320989 | jac(2,5) = f%d2pth(1)*(x(2) - si%z(2)) - si%dt*f%d2H(1) | |
| 571 |
2/2✓ Branch 0 taken 6641978 times.
✓ Branch 1 taken 3320989 times.
|
9962967 | jac(2,2:3) = 0.5d0*(f%d2pth(2:3)*(x(2) - si%z(2)) - si%dt*f%d2H(2:3)) |
| 572 | 3320989 | jac(2,2) = jac(2,2) + f%dpth(1) | |
| 573 | 3320989 | jac(2,4) = 0.5d0*(f%d2pth(7)*(x(2) - si%z(2)) - si%dt*f%d2H(7)) | |
| 574 | |||
| 575 | 3320989 | jac(3,1) = 0d0 | |
| 576 | jac(3,5) = (f%d2pth(1)*f%hph + f%dpth(1)*f%dhph(1))*(x(3) - si%z(3)) & | ||
| 577 | - si%dt*(f%d2pth(1)*f%vpar + f%dpth(1)*f%dvpar(1) & | ||
| 578 | 3320989 | - f%d2H(1)*f%hth - f%dH(1)*f%dhth(1)) | |
| 579 | jac(3,2:3) = 0.5d0*((f%d2pth(2:3)*f%hph + f%dpth(1)*f%dhph(2:3))*(x(3) - si%z(3)) & | ||
| 580 | - si%dt*(f%d2pth(2:3)*f%vpar + f%dpth(1)*f%dvpar(2:3) & | ||
| 581 |
2/2✓ Branch 0 taken 6641978 times.
✓ Branch 1 taken 3320989 times.
|
9962967 | - f%d2H(2:3)*f%hth - f%dH(1)*f%dhth(2:3))) |
| 582 | 3320989 | jac(3,3) = jac(3,3) + f%dpth(1)*f%hph | |
| 583 | jac(3,4) = 0.5d0*(f%d2pth(7)*f%hph*(x(3) - si%z(3)) & | ||
| 584 | 3320989 | - si%dt*(f%d2pth(7)*f%vpar + f%dpth(1)*f%dvpar(4) - f%d2H(7)*f%hth)) | |
| 585 | |||
| 586 | 3320989 | jac(4,1) = 0d0 | |
| 587 | jac(4,5) = f%d2pth(1)*(x(4) - si%z(4)) & | ||
| 588 | + si%dt*(f%d2H(3)*f%dpth(1) + f%dH(3)*f%d2pth(1) & | ||
| 589 | 3320989 | - f%d2H(1)*f%dpth(3) - f%dH(1)*f%d2pth(3)) | |
| 590 | jac(4,2) = 0.5d0*(f%d2pth(2)*(x(4) - si%z(4)) & | ||
| 591 | + si%dt*(f%d2H(5)*f%dpth(1) + f%dH(3)*f%d2pth(2) & | ||
| 592 | 3320989 | - f%d2H(2)*f%dpth(3) - f%dH(1)*f%d2pth(5))) | |
| 593 | jac(4,3) = 0.5d0*(f%d2pth(3)*(x(4) - si%z(4)) & | ||
| 594 | + si%dt*(f%d2H(6)*f%dpth(1) + f%dH(3)*f%d2pth(3) & | ||
| 595 | 3320989 | - f%d2H(3)*f%dpth(3) - f%dH(1)*f%d2pth(6))) | |
| 596 | jac(4,4) = 0.5d0*(f%d2pth(7)*(x(4) - si%z(4)) & | ||
| 597 | + si%dt*(f%dH(3)*f%d2pth(7) + f%d2H(9)*f%dpth(1)& | ||
| 598 | 3320989 | - f%d2H(7)*f%dpth(3) - f%dH(1)*f%d2pth(9))) | |
| 599 | 3320989 | jac(4,4) = jac(4,4) + f%dpth(1) | |
| 600 | |||
| 601 | 3320989 | jac(5,1) = 0d0 | |
| 602 | jac(5,5) = f%d2pth(1)*(f%pth - si%pthold) + f%dpth(1)*f%dpth(1)& | ||
| 603 | + si%dt/2.0d0*(f%d2pth(1)*f%dH(2) + f%dpth(1)*f%d2H(2) & | ||
| 604 | 3320989 | - f%d2pth(2)*f%dH(1) - f%dpth(2)*f%d2H(1)) | |
| 605 | jac(5,2) = 0.5d0*(f%d2pth(2)*(f%pth - si%pthold) + f%dpth(1)*f%dpth(2) & | ||
| 606 | + si%dt/2.0d0*(f%d2pth(2)*f%dH(2) + f%dpth(1)*f%d2H(4) & | ||
| 607 | 3320989 | - f%d2pth(4)*f%dH(1) - f%dpth(2)*f%d2H(2))) | |
| 608 | jac(5,3) = 0.5d0*(f%d2pth(3)*(f%pth - si%pthold) + f%dpth(1)*f%dpth(3) & | ||
| 609 | + si%dt/2.0d0*(f%d2pth(3)*f%dH(2) + f%dpth(1)*f%d2H(5) & | ||
| 610 | 3320989 | - f%d2pth(5)*f%dH(1) - f%dpth(2)*f%d2H(3))) | |
| 611 | jac(5,4) = 0.5d0*(f%d2pth(7)*(f%pth - si%pthold) + f%dpth(1)*f%dpth(4) & | ||
| 612 | + si%dt/2.0d0*(f%d2pth(7)*f%dH(2) + f%dpth(1)*f%d2H(8) & | ||
| 613 | 3320989 | - f%d2pth(8)*f%dH(1) - f%dpth(2)*f%d2H(7))) | |
| 614 | |||
| 615 | 3320989 | end subroutine jac_midpoint_part1 | |
| 616 | |||
| 617 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 618 | ! | ||
| 619 | 3320989 | recursive subroutine jac_midpoint_part2(si, f, fmid, x, jac) | |
| 620 | ! | ||
| 621 | type(symplectic_integrator_t), intent(in) :: si | ||
| 622 | type(field_can_t), intent(inout) :: f | ||
| 623 | type(field_can_t), intent(inout) :: fmid | ||
| 624 | real(dp), intent(in) :: x(5) | ||
| 625 | real(dp), intent(out) :: jac(5, 5) | ||
| 626 | |||
| 627 | ! fmid%dpth(1)*(f%pth - si%pthold) + si%dt*(fmid%dpth(1)*fmid%dH(2)-fmid%dpth(2)*fmid%dH(1)) | ||
| 628 | |||
| 629 | 3320989 | jac(1,1) = fmid%dpth(1)*f%dpth(1) | |
| 630 | jac(1,2) = 0.5d0*(fmid%d2pth(2)*(f%pth - si%pthold) & | ||
| 631 | + si%dt*(fmid%d2pth(2)*fmid%dH(2) + fmid%dpth(1)*fmid%d2H(4) & | ||
| 632 | 3320989 | - fmid%dpth(2)*fmid%d2H(2) - fmid%d2pth(4)*fmid%dH(1))) + fmid%dpth(1)*f%dpth(2) | |
| 633 | jac(1,3) = 0.5d0*(fmid%d2pth(3)*(f%pth - si%pthold) & | ||
| 634 | + si%dt*(fmid%d2pth(3)*fmid%dH(2) + fmid%dpth(1)*fmid%d2H(5) & | ||
| 635 | 3320989 | - fmid%dpth(2)*fmid%d2H(3) - fmid%d2pth(5)*fmid%dH(1))) + fmid%dpth(1)*f%dpth(3) | |
| 636 | jac(1,4) = 0.5d0*(fmid%d2pth(7)*(f%pth - si%pthold) & | ||
| 637 | + si%dt*(fmid%d2pth(7)*fmid%dH(2) + fmid%dpth(1)*fmid%d2H(8) & | ||
| 638 | 3320989 | - fmid%dpth(2)*fmid%d2H(7) - fmid%d2pth(8)*fmid%dH(1))) + fmid%dpth(1)*f%dpth(4) | |
| 639 | jac(1,5) = fmid%d2pth(1)*(f%pth - si%pthold) & | ||
| 640 | + si%dt*(fmid%d2pth(1)*fmid%dH(2) + fmid%dpth(1)*fmid%d2H(2) & | ||
| 641 | 3320989 | - fmid%dpth(2)*fmid%d2H(1) - fmid%d2pth(2)*fmid%dH(1)) | |
| 642 | |||
| 643 | 3320989 | end subroutine jac_midpoint_part2 | |
| 644 | |||
| 645 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 646 | ! | ||
| 647 | 4403570 | recursive subroutine newton1(si, f, x, maxit, xlast, status) | |
| 648 | ! | ||
| 649 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 650 | type(field_can_t), intent(inout) :: f | ||
| 651 | integer, parameter :: n = 2 | ||
| 652 | integer, parameter :: radial_indices(1) = [1] | ||
| 653 | |||
| 654 | real(dp), intent(inout) :: x(n) | ||
| 655 | integer, intent(in) :: maxit | ||
| 656 | real(dp), intent(out) :: xlast(n) | ||
| 657 | integer, intent(out) :: status | ||
| 658 | |||
| 659 | real(dp) :: tolref(n) | ||
| 660 | integer :: kit | ||
| 661 | logical :: converged, linear_failed, boundary_limited, step_boundary_limited | ||
| 662 | |||
| 663 | 4403570 | status = SYMPLECTIC_STEP_MAXITER | |
| 664 | 4403570 | boundary_limited = .false. | |
| 665 | |||
| 666 | 4403570 | tolref(1) = 1d0 | |
| 667 | 4403570 | tolref(2) = dabs(1d1*torflux/f%ro0) | |
| 668 | |||
| 669 |
2/2✓ Branch 0 taken 12947664 times.
✓ Branch 1 taken 130 times.
|
12947794 | do kit = 1, maxit |
| 670 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12947664 times.
|
12947664 | if (x(1) > sympl_rmax) then |
| 671 | ✗ | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 672 | 4403440 | return | |
| 673 | end if | ||
| 674 | ! Transient guard: in s = rho^2 coordinates the Hamiltonian is not | ||
| 675 | ! smooth at the axis (sqrt(s) behavior), so there is no consistent | ||
| 676 | ! field extension to s < 0 for the solver itself. Intermediate | ||
| 677 | ! negative iterates are floored as before; a converged negative | ||
| 678 | ! solution is committed as an axis crossing by the caller (#370). | ||
| 679 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12947664 times.
|
12947664 | if (x(1) < 0d0) x(1) = 0.01d0 |
| 680 | |||
| 681 | 12947664 | call eval_field(f, x(1), si%z(2), si%z(3), 2) | |
| 682 | 12947664 | call get_derivatives2(f, x(2)) | |
| 683 | call sympl_euler1_newton_iter(si, f, x, tolref, xlast, converged, & | ||
| 684 | 12947664 | linear_failed, step_boundary_limited) | |
| 685 |
4/4✓ Branch 0 taken 12938782 times.
✓ Branch 1 taken 8882 times.
✓ Branch 2 taken 578 times.
✓ Branch 3 taken 12938204 times.
|
12947664 | boundary_limited = boundary_limited .or. step_boundary_limited |
| 686 | |||
| 687 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12947664 times.
|
12947664 | if (linear_failed) then |
| 688 | ✗ | status = SYMPLECTIC_STEP_LINEAR_SOLVE | |
| 689 | ✗ | return | |
| 690 | end if | ||
| 691 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12947664 times.
|
12947664 | if (x(1) > sympl_rmax) then |
| 692 | ✗ | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 693 | ✗ | return | |
| 694 | end if | ||
| 695 |
2/2✓ Branch 1 taken 4403440 times.
✓ Branch 2 taken 8544224 times.
|
25895458 | if (converged) then |
| 696 | 4403440 | status = SYMPLECTIC_STEP_OK | |
| 697 | 4403440 | return | |
| 698 | end if | ||
| 699 | enddo | ||
| 700 |
2/2✓ Branch 0 taken 128 times.
✓ Branch 1 taken 2 times.
|
130 | if (boundary_limited) then |
| 701 |
2/4✓ Branch 0 taken 128 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 128 times.
✗ Branch 3 not taken.
|
128 | if (step_boundary_limited .and. & |
| 702 | radial_boundary_reached(x, radial_indices, si%rtol)) then | ||
| 703 | 128 | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 704 | else | ||
| 705 | ✗ | status = SYMPLECTIC_STEP_BOUNDARY_LIMITED | |
| 706 | end if | ||
| 707 | else | ||
| 708 | 2 | call count_event(EVT_NEWTON1_MAXIT) | |
| 709 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
|
2 | if (accept_warning_maxiter(x, xlast, tolref, si%rtol, & |
| 710 | si%warning_factor)) & | ||
| 711 | ✗ | status = SYMPLECTIC_STEP_OK | |
| 712 | end if | ||
| 713 | end subroutine | ||
| 714 | |||
| 715 | 4363 | recursive subroutine newton2(si, f, x, atol, rtol, maxit, xlast, status) | |
| 716 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 717 | type(field_can_t), intent(inout) :: f | ||
| 718 | |||
| 719 | integer, parameter :: n = 3 | ||
| 720 | integer, parameter :: radial_indices(1) = [1] | ||
| 721 | integer :: kit | ||
| 722 | |||
| 723 | real(dp), intent(inout) :: x(n) | ||
| 724 | real(dp), intent(in) :: atol, rtol | ||
| 725 | integer, intent(in) :: maxit | ||
| 726 | real(dp), intent(out) :: xlast(n) | ||
| 727 | integer, intent(out) :: status | ||
| 728 | |||
| 729 | real(dp) :: fvec(n), fjac(n,n), jinv(n,n) | ||
| 730 | |||
| 731 | real(dp) :: xabs(n), tolref(n), fabs(n), correction(n) | ||
| 732 | real(dp) :: det, step_scale | ||
| 733 | logical :: boundary_limited, step_limited | ||
| 734 | |||
| 735 | 4363 | status = SYMPLECTIC_STEP_MAXITER | |
| 736 | 4363 | boundary_limited = .false. | |
| 737 | |||
| 738 |
2/2✓ Branch 0 taken 17289 times.
✓ Branch 1 taken 1 time.
|
17290 | do kit = 1, maxit |
| 739 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 17289 times.
|
17289 | if(x(1) > sympl_rmax) then |
| 740 | ✗ | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 741 | 4362 | return | |
| 742 | end if | ||
| 743 | ! Transient guard for intermediate iterates; the converged-negative | ||
| 744 | ! case is handled by the caller via a chart switch (#370). | ||
| 745 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 17289 times.
|
17289 | if(x(1) < 0.0) x(1) = 0.01 |
| 746 | 17289 | call f_sympl_euler2(si, f, n, x, fvec, 1) | |
| 747 |
2/2✓ Branch 0 taken 51867 times.
✓ Branch 1 taken 17289 times.
|
69156 | fabs = dabs(fvec) |
| 748 | 17289 | call jac_sympl_euler2(si, f, x, fjac) | |
| 749 | 17289 | xlast = x | |
| 750 | |||
| 751 | det = fjac(1,1)*fjac(2,2)*fjac(3,3) + fjac(1,2)*fjac(2,3)*fjac(3,1) + fjac(1,3)*fjac(2,1)*fjac(3,2) & | ||
| 752 | 17289 | - fjac(1,3)*fjac(2,2)*fjac(3,1) - fjac(1,1)*fjac(2,3)*fjac(3,2) - fjac(1,2)*fjac(2,1)*fjac(3,3) | |
| 753 | |||
| 754 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 17289 times.
|
17289 | if (matrix3_near_singular(fjac, det)) then |
| 755 | ✗ | status = SYMPLECTIC_STEP_LINEAR_SOLVE | |
| 756 | ✗ | return | |
| 757 | end if | ||
| 758 | |||
| 759 | 17289 | jinv(1,1) = fjac(2,2)*fjac(3,3) - fjac(2,3)*fjac(3,2) | |
| 760 | 17289 | jinv(1,2) = fjac(1,3)*fjac(3,2) - fjac(1,2)*fjac(3,3) | |
| 761 | 17289 | jinv(1,3) = fjac(1,2)*fjac(2,3) - fjac(1,3)*fjac(2,2) | |
| 762 | |||
| 763 | 17289 | jinv(2,1) = fjac(2,3)*fjac(3,1) - fjac(2,1)*fjac(3,3) | |
| 764 | 17289 | jinv(2,2) = fjac(1,1)*fjac(3,3) - fjac(1,3)*fjac(3,1) | |
| 765 | 17289 | jinv(2,3) = fjac(1,3)*fjac(2,1) - fjac(1,1)*fjac(2,3) | |
| 766 | |||
| 767 | 17289 | jinv(3,1) = fjac(2,1)*fjac(3,2) - fjac(3,1)*fjac(2,2) | |
| 768 | 17289 | jinv(3,2) = fjac(1,2)*fjac(3,1) - fjac(3,2)*fjac(1,1) | |
| 769 | 17289 | jinv(3,3) = fjac(1,1)*fjac(2,2) - fjac(1,2)*fjac(2,1) | |
| 770 | |||
| 771 |
6/6✓ Branch 0 taken 51867 times.
✓ Branch 1 taken 17289 times.
✓ Branch 2 taken 155601 times.
✓ Branch 3 taken 51867 times.
✓ Branch 4 taken 51867 times.
✓ Branch 5 taken 17289 times.
|
276624 | correction = matmul(jinv, fvec)/det |
| 772 | call limit_radial_newton_step(x, correction, radial_indices, step_scale, & | ||
| 773 | 17289 | step_limited) | |
| 774 |
4/4✓ Branch 0 taken 17227 times.
✓ Branch 1 taken 62 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 17224 times.
|
17289 | boundary_limited = boundary_limited .or. step_limited |
| 775 |
2/2✓ Branch 0 taken 51867 times.
✓ Branch 1 taken 17289 times.
|
69156 | x = x - step_scale*correction |
| 776 | |||
| 777 | !call dgesv(n, 1, fjac, n, pivot, fvec, n, info) | ||
| 778 | ! after solution: fvec = (xold-xnew)_Newton | ||
| 779 | !x = x - fvec | ||
| 780 | |||
| 781 |
2/2✓ Branch 0 taken 51867 times.
✓ Branch 1 taken 17289 times.
|
69156 | xabs = dabs(x-xlast) |
| 782 | |||
| 783 | 17289 | tolref(1) = 1d0 | |
| 784 | 17289 | tolref(2) = twopi | |
| 785 | 17289 | tolref(3) = twopi | |
| 786 | |||
| 787 |
8/8✓ Branch 0 taken 19447 times.
✓ Branch 1 taken 175 times.
✓ Branch 2 taken 2333 times.
✓ Branch 3 taken 17114 times.
✓ Branch 4 taken 17227 times.
✓ Branch 5 taken 62 times.
✓ Branch 6 taken 170 times.
✓ Branch 7 taken 17057 times.
|
19622 | if (.not. step_limited .and. all(fabs < atol)) then |
| 788 | 170 | status = SYMPLECTIC_STEP_OK | |
| 789 | 170 | return | |
| 790 | end if | ||
| 791 |
8/8✓ Branch 0 taken 28290 times.
✓ Branch 1 taken 4210 times.
✓ Branch 2 taken 15381 times.
✓ Branch 3 taken 12909 times.
✓ Branch 4 taken 17057 times.
✓ Branch 5 taken 62 times.
✓ Branch 6 taken 4192 times.
✓ Branch 7 taken 12865 times.
|
32501 | if (.not. step_limited .and. all(xabs < rtol*tolref)) then |
| 792 | 4192 | status = SYMPLECTIC_STEP_OK | |
| 793 | 4192 | return | |
| 794 | end if | ||
| 795 | enddo | ||
| 796 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (boundary_limited) then |
| 797 | ✗ | if (step_limited .and. & | |
| 798 | radial_boundary_reached(x, radial_indices, rtol)) then | ||
| 799 | ✗ | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 800 | else | ||
| 801 | ✗ | status = SYMPLECTIC_STEP_BOUNDARY_LIMITED | |
| 802 | end if | ||
| 803 | else | ||
| 804 | 1 | call count_event(EVT_NEWTON2_MAXIT) | |
| 805 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
|
1 | if (accept_warning_maxiter(x, xlast, tolref, rtol, si%warning_factor)) & |
| 806 | ✗ | status = SYMPLECTIC_STEP_OK | |
| 807 | end if | ||
| 808 | end subroutine | ||
| 809 | |||
| 810 | 17290 | pure logical function matrix3_near_singular(matrix, determinant) | |
| 811 | real(dp), intent(in) :: matrix(3, 3), determinant | ||
| 812 | real(dp) :: matrix_scale | ||
| 813 | |||
| 814 |
10/12✓ Branch 0 taken 51870 times.
✓ Branch 1 taken 17290 times.
✓ Branch 2 taken 155610 times.
✓ Branch 3 taken 51870 times.
✓ Branch 4 taken 138320 times.
✓ Branch 5 taken 17290 times.
✓ Branch 6 taken 34022 times.
✓ Branch 7 taken 104298 times.
✓ Branch 8 taken 17290 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 17290 times.
|
224770 | matrix_scale = maxval(abs(matrix)) |
| 815 | matrix3_near_singular = matrix_scale == 0d0 .or. & | ||
| 816 |
2/4✓ Branch 0 taken 17290 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 17290 times.
|
17290 | abs(determinant) <= epsilon(1d0)*matrix_scale**3 |
| 817 | 17290 | end function matrix3_near_singular | |
| 818 | |||
| 819 | ! Small dense Newton solve (n = 2 .. 4*s) by in-house LU with partial pivoting. | ||
| 820 | ! No library call and no allocation: LAPACK/OpenBLAS dgesv per step exhausted | ||
| 821 | ! OpenBLAS memory regions and scaled badly with thread count. | ||
| 822 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 3326289 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3326289 times.
|
3326289 | subroutine solve_newton_system(matrix, residual, status) |
| 823 | real(dp), intent(inout) :: matrix(:, :), residual(:) | ||
| 824 | integer, intent(out) :: status | ||
| 825 | |||
| 826 | integer :: info, n | ||
| 827 | |||
| 828 | 3326289 | n = size(residual) | |
| 829 |
4/24✗ Branch 0 not taken.
✓ Branch 1 taken 3326289 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 3326289 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 3326289 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 3326289 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
3326289 | call rk_solve(n, matrix, residual, info) |
| 830 |
2/2✓ Branch 0 taken 3326288 times.
✓ Branch 1 taken 1 time.
|
3326289 | if (info == 0) then |
| 831 | 3326288 | status = SYMPLECTIC_STEP_OK | |
| 832 | else | ||
| 833 | 1 | status = SYMPLECTIC_STEP_LINEAR_SOLVE | |
| 834 | end if | ||
| 835 | 3326289 | end subroutine solve_newton_system | |
| 836 | |||
| 837 | 1216 | pure logical function boundary_event_converged(fraction_width, radial_residual, & | |
| 838 | fraction_tolerance, radial_tolerance) | ||
| 839 | real(dp), intent(in) :: fraction_width, radial_residual | ||
| 840 | real(dp), intent(in) :: fraction_tolerance, radial_tolerance | ||
| 841 | |||
| 842 | boundary_event_converged = fraction_width <= fraction_tolerance .and. & | ||
| 843 |
4/4✓ Branch 0 taken 74 times.
✓ Branch 1 taken 1142 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 72 times.
|
1216 | radial_residual <= radial_tolerance |
| 844 | 1216 | end function boundary_event_converged | |
| 845 | |||
| 846 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 1078905 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1078905 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1078905 times.
|
1078905 | pure subroutine limit_radial_newton_step(x, correction, radial_indices, & |
| 847 | step_scale, limited) | ||
| 848 | real(dp), intent(in) :: x(:), correction(:) | ||
| 849 | integer, intent(in) :: radial_indices(:) | ||
| 850 | real(dp), intent(out) :: step_scale | ||
| 851 | logical, intent(out) :: limited | ||
| 852 | |||
| 853 | real(dp) :: outward_step | ||
| 854 | integer :: k, radius_index | ||
| 855 | |||
| 856 | 1078905 | step_scale = 1d0 | |
| 857 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1078905 times.
|
1078905 | if (sympl_rmax > 1d0) then |
| 858 | ✗ | limited = .false. | |
| 859 | ✗ | return | |
| 860 | end if | ||
| 861 |
2/2✓ Branch 0 taken 2145724 times.
✓ Branch 1 taken 1078905 times.
|
3224629 | do k = 1, size(radial_indices) |
| 862 | 2145724 | radius_index = radial_indices(k) | |
| 863 | 2145724 | outward_step = -correction(radius_index) | |
| 864 |
2/2✓ Branch 0 taken 1084637 times.
✓ Branch 1 taken 1061087 times.
|
3224629 | if (outward_step > 0d0) then |
| 865 | step_scale = min(step_scale, & | ||
| 866 | 1084637 | 0.8d0*max(0d0, sympl_rmax - x(radius_index))/outward_step) | |
| 867 | end if | ||
| 868 | end do | ||
| 869 | 1078905 | limited = step_scale < 1d0 | |
| 870 | end subroutine limit_radial_newton_step | ||
| 871 | |||
| 872 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 261 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 261 times.
|
261 | pure logical function radial_boundary_reached(x, radial_indices, rtol) |
| 873 | real(dp), intent(in) :: x(:), rtol | ||
| 874 | integer, intent(in) :: radial_indices(:) | ||
| 875 | |||
| 876 | radial_boundary_reached = any(abs(sympl_rmax - x(radial_indices)) <= & | ||
| 877 |
3/4✓ Branch 0 taken 265 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 261 times.
|
265 | max(1d-10, 10d0*rtol)) |
| 878 | 261 | end function radial_boundary_reached | |
| 879 | |||
| 880 | 1147410 | recursive subroutine newton_midpoint(si, f, x, atol, rtol, maxit, xlast, status) | |
| 881 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 882 | type(field_can_t), intent(inout) :: f | ||
| 883 | |||
| 884 | type(field_can_t) :: fmid | ||
| 885 | |||
| 886 | integer, parameter :: n = 5 | ||
| 887 | integer :: kit | ||
| 888 | integer, parameter :: radial_indices(2) = [1, 5] | ||
| 889 | |||
| 890 | real(dp), intent(inout) :: x(n) ! = (rend, thend, phend, pphend, rmid) | ||
| 891 | real(dp), intent(in) :: atol, rtol | ||
| 892 | integer, intent(in) :: maxit | ||
| 893 | real(dp), intent(out) :: xlast(n) | ||
| 894 | integer, intent(out) :: status | ||
| 895 | |||
| 896 | real(dp) :: fvec(n), fjac(n,n) | ||
| 897 | |||
| 898 | real(dp) :: xabs(n), tolref(n), fabs(n), step_scale | ||
| 899 | logical :: boundary_limited, step_limited | ||
| 900 | |||
| 901 | 1147410 | xlast = x | |
| 902 | 1147410 | status = SYMPLECTIC_STEP_MAXITER | |
| 903 | 1147410 | boundary_limited = .false. | |
| 904 | |||
| 905 | 1147410 | tolref(1) = 1d0 | |
| 906 | 1147410 | tolref(2) = twopi | |
| 907 | 1147410 | tolref(3) = twopi | |
| 908 | 1147410 | tolref(4) = max(dabs(f%Aph), dabs(1d1*torflux/f%ro0)) | |
| 909 | 1147410 | tolref(5) = 1d0 | |
| 910 | |||
| 911 |
2/2✓ Branch 0 taken 3320990 times.
✓ Branch 1 taken 90 times.
|
3321080 | do kit = 1, maxit |
| 912 |
3/4✓ Branch 0 taken 3320989 times.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 3320989 times.
|
3320990 | if(x(1) > sympl_rmax .or. x(5) > sympl_rmax) then |
| 913 | 1 | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 914 | 1147320 | return | |
| 915 | end if | ||
| 916 | ! Transient guards for intermediate iterates; the converged-negative | ||
| 917 | ! case is handled by the caller via a chart switch (#370). | ||
| 918 |
2/2✓ Branch 0 taken 62 times.
✓ Branch 1 taken 3320927 times.
|
3320989 | if(x(1) < 0.0) x(1) = 0.01 |
| 919 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3320989 times.
|
3320989 | if(x(5) < 0.0) x(5) = 0.01 |
| 920 | 3320989 | call f_midpoint_part1(si, f, n, x, fvec) | |
| 921 | 3320989 | call jac_midpoint_part1(si, f, x, fjac) | |
| 922 | 3320989 | fmid = f | |
| 923 | 3320989 | call f_midpoint_part2(si, f, n, x, fvec) | |
| 924 | 3320989 | call jac_midpoint_part2(si, f, fmid, x, fjac) | |
| 925 |
2/2✓ Branch 0 taken 16604945 times.
✓ Branch 1 taken 3320989 times.
|
19925934 | fabs = dabs(fvec) |
| 926 | 3320989 | xlast = x | |
| 927 | 3320989 | call solve_newton_system(fjac, fvec, status) | |
| 928 |
1/2✓ Branch 0 taken 3320989 times.
✗ Branch 1 not taken.
|
3320989 | if (status /= SYMPLECTIC_STEP_OK) return |
| 929 | 3320989 | status = SYMPLECTIC_STEP_MAXITER | |
| 930 |
2/2✓ Branch 0 taken 2264671 times.
✓ Branch 1 taken 1056318 times.
|
3320989 | if (sympl_rmax > 1d0) then |
| 931 | 2264671 | step_limited = .false. | |
| 932 |
2/2✓ Branch 0 taken 11323355 times.
✓ Branch 1 taken 2264671 times.
|
13588026 | x = x - fvec |
| 933 | else | ||
| 934 | call limit_radial_newton_step(x, fvec, radial_indices, step_scale, & | ||
| 935 | 1056318 | step_limited) | |
| 936 |
4/4✓ Branch 0 taken 1049985 times.
✓ Branch 1 taken 6333 times.
✓ Branch 2 taken 443 times.
✓ Branch 3 taken 1049542 times.
|
1056318 | boundary_limited = boundary_limited .or. step_limited |
| 937 |
2/2✓ Branch 0 taken 5281590 times.
✓ Branch 1 taken 1056318 times.
|
6337908 | x = x - step_scale*fvec |
| 938 | end if | ||
| 939 |
4/6✓ Branch 0 taken 6641978 times.
✓ Branch 1 taken 3320989 times.
✓ Branch 2 taken 6641978 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3320989 times.
|
9962967 | if (any(x(radial_indices) > sympl_rmax)) then |
| 940 | ✗ | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 941 | return | ||
| 942 | end if | ||
| 943 |
2/2✓ Branch 0 taken 16604945 times.
✓ Branch 1 taken 3320989 times.
|
19925934 | xabs = dabs(x - xlast) |
| 944 | |||
| 945 | ! Don't take too small values in pphi as tolerance reference | ||
| 946 | 3320989 | tolref(4) = max(dabs(x(4)), tolref(4)) | |
| 947 | |||
| 948 |
8/8✓ Branch 0 taken 3359668 times.
✓ Branch 1 taken 318 times.
✓ Branch 2 taken 38997 times.
✓ Branch 3 taken 3320671 times.
✓ Branch 4 taken 3314920 times.
✓ Branch 5 taken 6069 times.
✓ Branch 6 taken 280 times.
✓ Branch 7 taken 3314640 times.
|
3359986 | if (.not. step_limited .and. all(fabs < atol)) then |
| 949 | 280 | status = SYMPLECTIC_STEP_OK | |
| 950 | 280 | return | |
| 951 | end if | ||
| 952 |
8/8✓ Branch 0 taken 8056729 times.
✓ Branch 1 taken 1148608 times.
✓ Branch 2 taken 5884628 times.
✓ Branch 3 taken 2172101 times.
✓ Branch 4 taken 3314640 times.
✓ Branch 5 taken 6069 times.
✓ Branch 6 taken 1147039 times.
✓ Branch 7 taken 2167601 times.
|
9205427 | if (.not. step_limited .and. all(xabs < rtol*tolref)) then |
| 953 | 1147039 | status = SYMPLECTIC_STEP_OK | |
| 954 | 1147039 | return | |
| 955 | end if | ||
| 956 | enddo | ||
| 957 |
2/2✓ Branch 0 taken 85 times.
✓ Branch 1 taken 5 times.
|
90 | if (boundary_limited) then |
| 958 |
2/4✓ Branch 0 taken 85 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 85 times.
✗ Branch 3 not taken.
|
85 | if (step_limited .and. & |
| 959 | radial_boundary_reached(x, radial_indices, rtol)) then | ||
| 960 | 85 | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 961 | else | ||
| 962 | ✗ | status = SYMPLECTIC_STEP_BOUNDARY_LIMITED | |
| 963 | end if | ||
| 964 | else | ||
| 965 | 5 | call count_event(EVT_MIDPOINT_MAXIT) | |
| 966 |
2/2✓ Branch 1 taken 2 times.
✓ Branch 2 taken 3 times.
|
5 | if (accept_warning_axis_crossing(si, x, xlast, tolref, rtol)) then |
| 967 | 2 | status = SYMPLECTIC_STEP_OK | |
| 968 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
|
3 | else if (accept_warning_maxiter(x, xlast, tolref, rtol, & |
| 969 | si%warning_factor)) then | ||
| 970 | ✗ | status = SYMPLECTIC_STEP_OK | |
| 971 | end if | ||
| 972 | end if | ||
| 973 |
4/4✓ Branch 0 taken 3442230 times.
✓ Branch 1 taken 1147410 times.
✓ Branch 2 taken 6884460 times.
✓ Branch 3 taken 1147410 times.
|
11474100 | end subroutine |
| 974 | |||
| 975 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 976 | ! | ||
| 977 | ! Gauss-Legendre Runge-Kutta method with s internal stages (n=4*s variables) | ||
| 978 | ! | ||
| 979 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2439 times.
|
2439 | recursive subroutine f_rk_gauss(si, fs, s, x, fvec) |
| 980 | ! GC Gauss residual: thin wrapper over the shared core (MODEL_GC). The | ||
| 981 | ! arithmetic lives in orbit_rk_core::gauss_canfield_residual. | ||
| 982 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 983 | integer, intent(in) :: s | ||
| 984 | type(field_can_t), intent(inout) :: fs(:) | ||
| 985 | real(dp), intent(in) :: x(4*s) ! = (rend, thend, phend, pphend) | ||
| 986 | real(dp), intent(out) :: fvec(4*s) | ||
| 987 | |||
| 988 | 2439 | call rk_gauss_residual(MODEL_GC, si, fs, s, x, fvec) | |
| 989 | 2439 | end subroutine f_rk_gauss | |
| 990 | |||
| 991 | |||
| 992 | |||
| 993 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 994 | ! | ||
| 995 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2439 times.
|
2439 | recursive subroutine jac_rk_gauss(si, fs, s, jac) |
| 996 | ! GC Gauss Jacobian: thin wrapper over the shared core (MODEL_GC). The | ||
| 997 | ! arithmetic lives in orbit_rk_core::gauss_canfield_jacobian. | ||
| 998 | type(symplectic_integrator_t), intent(in) :: si | ||
| 999 | integer, intent(in) :: s | ||
| 1000 | type(field_can_t), intent(in) :: fs(:) | ||
| 1001 | real(dp), intent(out) :: jac(4*s, 4*s) | ||
| 1002 | |||
| 1003 | 2439 | call rk_gauss_jacobian(MODEL_GC, si, fs, s, jac) | |
| 1004 | |||
| 1005 | 2439 | end subroutine jac_rk_gauss | |
| 1006 | |||
| 1007 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 544 times.
|
544 | recursive subroutine newton_rk_gauss(si, fs, s, x, atol, rtol, maxit, xlast, status) |
| 1008 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 1009 | |||
| 1010 | integer, intent(in) :: s | ||
| 1011 | type(field_can_t), intent(inout) :: fs(:) | ||
| 1012 | 1088 | integer :: kit, ks, radial_indices(s) | |
| 1013 | |||
| 1014 | real(dp), intent(inout) :: x(4*s) | ||
| 1015 | real(dp), intent(in) :: atol, rtol | ||
| 1016 | integer, intent(in) :: maxit | ||
| 1017 | real(dp), intent(out) :: xlast(4*s) | ||
| 1018 | integer, intent(out) :: status | ||
| 1019 | |||
| 1020 | 1088 | real(dp) :: fvec(4*s), fjac(4*s, 4*s) | |
| 1021 | |||
| 1022 | 1088 | real(dp) :: xabs(4*s), tolref(4*s), fabs(4*s), step_scale | |
| 1023 | logical :: boundary_limited, step_limited | ||
| 1024 | |||
| 1025 |
2/2✓ Branch 0 taken 6304 times.
✓ Branch 1 taken 544 times.
|
6848 | xlast = x |
| 1026 | 544 | status = SYMPLECTIC_STEP_MAXITER | |
| 1027 | 544 | boundary_limited = .false. | |
| 1028 |
2/2✓ Branch 0 taken 1576 times.
✓ Branch 1 taken 544 times.
|
2120 | do ks = 1, s |
| 1029 | 2120 | radial_indices(ks) = 4*ks - 3 | |
| 1030 | end do | ||
| 1031 | |||
| 1032 |
2/2✓ Branch 0 taken 2439 times.
✓ Branch 1 taken 4 times.
|
2443 | do kit = 1, maxit |
| 1033 | |||
| 1034 | ! Check if radius left the boundary | ||
| 1035 |
2/2✓ Branch 0 taken 7224 times.
✓ Branch 1 taken 2439 times.
|
9663 | do ks = 1, s |
| 1036 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7224 times.
|
7224 | if (x(4*ks-3) > sympl_rmax) then |
| 1037 | ✗ | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 1038 | ✗ | return | |
| 1039 | end if | ||
| 1040 | ! Transient guard for intermediate iterates; the converged-negative | ||
| 1041 | ! case is handled by the caller via a chart switch (#370). | ||
| 1042 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7224 times.
|
9663 | if (x(4*ks-3) < 0.0) x(4*ks-3) = 0.01d0 |
| 1043 | end do | ||
| 1044 | |||
| 1045 | 2439 | call f_rk_gauss(si, fs, s, x, fvec) | |
| 1046 | 2439 | call jac_rk_gauss(si, fs, s, fjac) | |
| 1047 |
2/2✓ Branch 0 taken 28896 times.
✓ Branch 1 taken 2439 times.
|
31335 | fabs = dabs(fvec) |
| 1048 |
2/2✓ Branch 0 taken 28896 times.
✓ Branch 1 taken 2439 times.
|
31335 | xlast = x |
| 1049 | 2439 | call solve_newton_system(fjac, fvec, status) | |
| 1050 |
1/2✓ Branch 0 taken 2439 times.
✗ Branch 1 not taken.
|
2439 | if (status /= SYMPLECTIC_STEP_OK) return |
| 1051 | 2439 | status = SYMPLECTIC_STEP_MAXITER | |
| 1052 | ! after solution: fvec = (xold-xnew)_Newton | ||
| 1053 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2439 times.
|
2439 | if (sympl_rmax > 1d0) then |
| 1054 | ✗ | step_limited = .false. | |
| 1055 | ✗ | x = x - fvec | |
| 1056 | else | ||
| 1057 | call limit_radial_newton_step(x, fvec, radial_indices, step_scale, & | ||
| 1058 | 2439 | step_limited) | |
| 1059 |
2/4✓ Branch 0 taken 2439 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2439 times.
|
2439 | boundary_limited = boundary_limited .or. step_limited |
| 1060 |
2/2✓ Branch 0 taken 28896 times.
✓ Branch 1 taken 2439 times.
|
31335 | x = x - step_scale*fvec |
| 1061 | end if | ||
| 1062 |
4/6✓ Branch 0 taken 7224 times.
✓ Branch 1 taken 2439 times.
✓ Branch 2 taken 7224 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2439 times.
|
9663 | if (any(x(radial_indices) > sympl_rmax)) then |
| 1063 | ✗ | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 1064 | return | ||
| 1065 | end if | ||
| 1066 |
2/2✓ Branch 0 taken 28896 times.
✓ Branch 1 taken 2439 times.
|
31335 | xabs = dabs(x - xlast) |
| 1067 | |||
| 1068 |
2/2✓ Branch 0 taken 7224 times.
✓ Branch 1 taken 2439 times.
|
9663 | do ks = 1, s |
| 1069 | 7224 | tolref(4*ks-3) = 1d0 | |
| 1070 | 7224 | tolref(4*ks-2) = twopi | |
| 1071 | 7224 | tolref(4*ks-1) = twopi | |
| 1072 | 9663 | tolref(4*ks) = dabs(xlast(4*ks)) | |
| 1073 | end do | ||
| 1074 | |||
| 1075 |
10/12✓ Branch 0 taken 7224 times.
✓ Branch 1 taken 2439 times.
✓ Branch 2 taken 7224 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 5190 times.
✓ Branch 5 taken 150 times.
✓ Branch 6 taken 2901 times.
✓ Branch 7 taken 2289 times.
✓ Branch 8 taken 2439 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 150 times.
✓ Branch 11 taken 2289 times.
|
12564 | if (.not. step_limited .and. all(x(radial_indices) >= 0d0) .and. & |
| 1076 | all(fabs < atol)) then | ||
| 1077 | 150 | status = SYMPLECTIC_STEP_OK | |
| 1078 | return | ||
| 1079 | end if | ||
| 1080 |
10/12✓ Branch 0 taken 6624 times.
✓ Branch 1 taken 2289 times.
✓ Branch 2 taken 6624 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6905 times.
✓ Branch 5 taken 390 times.
✓ Branch 6 taken 5006 times.
✓ Branch 7 taken 1899 times.
✓ Branch 8 taken 2289 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 390 times.
✓ Branch 11 taken 1899 times.
|
13919 | if (.not. step_limited .and. all(x(radial_indices) >= 0d0) .and. & |
| 1081 | 4 | all(xabs < rtol*tolref)) then | |
| 1082 | 390 | status = SYMPLECTIC_STEP_OK | |
| 1083 | return | ||
| 1084 | end if | ||
| 1085 | enddo | ||
| 1086 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (boundary_limited) then |
| 1087 | ✗ | if (step_limited .and. & | |
| 1088 | radial_boundary_reached(x, radial_indices, rtol)) then | ||
| 1089 | ✗ | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 1090 | else | ||
| 1091 | ✗ | status = SYMPLECTIC_STEP_BOUNDARY_LIMITED | |
| 1092 | end if | ||
| 1093 | else | ||
| 1094 | 4 | call count_event(EVT_RK_GAUSS_MAXIT) | |
| 1095 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
|
4 | if (accept_warning_maxiter(x, xlast, tolref, rtol, si%warning_factor)) & |
| 1096 | ✗ | status = SYMPLECTIC_STEP_OK | |
| 1097 | end if | ||
| 1098 | end subroutine newton_rk_gauss | ||
| 1099 | |||
| 1100 | |||
| 1101 | ✗ | recursive subroutine fixpoint_rk_gauss(si, fs, s, x, atol, rtol, maxit, xlast) | |
| 1102 | ! TODO: this doesn't work well yet | ||
| 1103 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 1104 | |||
| 1105 | integer, intent(in) :: s | ||
| 1106 | type(field_can_t), intent(inout) :: fs(:) | ||
| 1107 | integer :: kit, ks | ||
| 1108 | |||
| 1109 | real(dp), intent(inout) :: x(4*s) | ||
| 1110 | real(dp), intent(in) :: atol, rtol | ||
| 1111 | integer, intent(in) :: maxit | ||
| 1112 | real(dp), intent(out) :: xlast(4*s) | ||
| 1113 | |||
| 1114 | ✗ | real(dp) :: fvec(4*s) | |
| 1115 | |||
| 1116 | ✗ | real(dp) :: xabs(4*s), tolref(4*s), fabs(4*s) | |
| 1117 | |||
| 1118 | ✗ | real(dp) :: a(s,s), b(s), c(s), Hprime(s), dHprimedr(s) | |
| 1119 | integer :: k, l | ||
| 1120 | |||
| 1121 | real(dp) :: pthnew, dpthnewdr, damp | ||
| 1122 | |||
| 1123 | ✗ | call coeff_rk_gauss(s, a, b, c) ! TODO: move this to preprocessing | |
| 1124 | |||
| 1125 | ✗ | do kit = 1, 4096 | |
| 1126 | |||
| 1127 | ! Check if radius left the boundary | ||
| 1128 | ✗ | do ks = 1, s | |
| 1129 | ✗ | if (x(4*ks-3) > sympl_rmax) return | |
| 1130 | ! Transient guard for intermediate iterates; the converged-negative | ||
| 1131 | ! case is handled by the caller via a chart switch (#370). | ||
| 1132 | ✗ | if (x(4*ks-3) < 0.0) x(4*ks-3) = 0.01d0 | |
| 1133 | end do | ||
| 1134 | |||
| 1135 | ✗ | call f_rk_gauss(si, fs, s, x, fvec) | |
| 1136 | ✗ | fabs = dabs(fvec) | |
| 1137 | ✗ | xlast = x | |
| 1138 | |||
| 1139 | ✗ | do k = 1, s | |
| 1140 | ✗ | call eval_field(fs(k), x(4*k-3), x(4*k-2), x(4*k-1), 2) | |
| 1141 | ✗ | call get_derivatives2(fs(k), x(4*k)) | |
| 1142 | ✗ | Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) | |
| 1143 | ✗ | dHprimedr(k) = (fs(k)%d2H(1) - Hprime(k)*fs(k)%d2pth(1))/fs(k)%dpth(1) | |
| 1144 | end do | ||
| 1145 | |||
| 1146 | ✗ | k = 1 | |
| 1147 | ✗ | l = 1 | |
| 1148 | |||
| 1149 | ✗ | pthnew = si%pthold - si%dt*a(k,l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) | |
| 1150 | dpthnewdr = -si%dt*a(k,l)*( & | ||
| 1151 | ✗ | fs(l)%d2H(2) - Hprime(l)*fs(l)%d2pth(2) - dHprimedr(l)*fs(l)%dpth(2)) | |
| 1152 | !print *, pthnew, dpthnewdr, pthnew/dpthnewdr | ||
| 1153 | |||
| 1154 | ✗ | damp = 0.99d0 | |
| 1155 | |||
| 1156 | ✗ | x(1) = x(1) - (1d0-damp)*pthnew/dpthnewdr | |
| 1157 | ✗ | x(2) = damp*x(2) + (1d0-damp)*(si%z(2) + si%dt*a(k,l)*si%dt*a(k,l)*Hprime(l)) | |
| 1158 | x(3) = damp*x(3) + (1d0-damp)*(si%z(3) + si%dt*a(k,l)*(fs(l)%vpar & | ||
| 1159 | ✗ | - Hprime(l)*fs(l)%hth)/fs(l)%hph) | |
| 1160 | x(4) = damp*x(4) + (1d0-damp)*(si%z(4) - si%dt*a(k,l)*(fs(l)%dH(3) & | ||
| 1161 | ✗ | - Hprime(l)*fs(l)%dpth(3))) | |
| 1162 | |||
| 1163 | ✗ | print *, x | |
| 1164 | |||
| 1165 | ✗ | xabs = dabs(x - xlast) | |
| 1166 | |||
| 1167 | ✗ | do ks = 1, s | |
| 1168 | ✗ | xabs(4*ks-2) = modulo(xabs(4*ks-2), pi) | |
| 1169 | ✗ | xabs(4*ks-1) = modulo(xabs(4*ks-1), pi) | |
| 1170 | end do | ||
| 1171 | |||
| 1172 | ✗ | do ks = 1, s | |
| 1173 | ✗ | tolref(4*ks-3) = 1d0 | |
| 1174 | ✗ | tolref(4*ks-2) = twopi | |
| 1175 | ✗ | tolref(4*ks-1) = twopi | |
| 1176 | ✗ | tolref(4*ks) = dabs(x(4*ks)) | |
| 1177 | end do | ||
| 1178 | |||
| 1179 | ✗ | if (all(fabs < atol)) return | |
| 1180 | ✗ | if (all(xabs < rtol*tolref)) return | |
| 1181 | enddo | ||
| 1182 | ✗ | call count_event(EVT_FIXPOINT_MAXIT) | |
| 1183 | ✗ | end subroutine fixpoint_rk_gauss | |
| 1184 | |||
| 1185 | |||
| 1186 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1187 | ! | ||
| 1188 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2858 times.
|
2858 | recursive subroutine jac_rk_lobatto(si, fs, s, jac) |
| 1189 | ! | ||
| 1190 | type(symplectic_integrator_t), intent(in) :: si | ||
| 1191 | integer, intent(in) :: s | ||
| 1192 | type(field_can_t), intent(in) :: fs(:) | ||
| 1193 | real(dp), intent(out) :: jac(4*s-2, 4*s-2) | ||
| 1194 | |||
| 1195 | 5716 | real(dp) :: a(s,s), ahat(s,s), b(s), c(s), Hprime(s), dHprime(4*s-2) | |
| 1196 | integer :: k,l,m,n ! counters | ||
| 1197 | |||
| 1198 | 2858 | call coeff_rk_lobatto(s, a, ahat, b, c) | |
| 1199 |
4/4✓ Branch 0 taken 2858 times.
✓ Branch 1 taken 28580 times.
✓ Branch 2 taken 285800 times.
✓ Branch 3 taken 28580 times.
|
317238 | jac = 0d0 |
| 1200 |
2/2✓ Branch 0 taken 2858 times.
✓ Branch 1 taken 28580 times.
|
31438 | dHprime = 0.0d0 |
| 1201 | |||
| 1202 |
2/2✓ Branch 0 taken 8574 times.
✓ Branch 1 taken 2858 times.
|
11432 | Hprime = 0.0d0 |
| 1203 | 2858 | Hprime(1) = fs(1)%dH(1)/fs(1)%dpth(1) | |
| 1204 |
2/2✓ Branch 0 taken 5716 times.
✓ Branch 1 taken 2858 times.
|
8574 | do k = 2, s |
| 1205 | 8574 | Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) | |
| 1206 | end do | ||
| 1207 | 2858 | dHprime(1) = (fs(1)%d2H(1)-Hprime(1)*fs(1)%d2pth(1))/fs(1)%dpth(1) ! d/dr | |
| 1208 | 2858 | dHprime(2) = (fs(1)%d2H(7)-Hprime(1)*fs(1)%d2pth(7))/fs(1)%dpth(1) ! d/dpph | |
| 1209 |
2/2✓ Branch 0 taken 5716 times.
✓ Branch 1 taken 2858 times.
|
8574 | do k = 2, s |
| 1210 | 5716 | m = 4*k-2 | |
| 1211 | 5716 | dHprime(m-3)=(fs(k)%d2H(1)-Hprime(k)*fs(k)%d2pth(1))/fs(k)%dpth(1) ! d/dr | |
| 1212 | 5716 | dHprime(m-2)=(fs(k)%d2H(2)-Hprime(k)*fs(k)%d2pth(2))/fs(k)%dpth(1) ! d/dth | |
| 1213 | 5716 | dHprime(m-1)=(fs(k)%d2H(3)-Hprime(k)*fs(k)%d2pth(3))/fs(k)%dpth(1) ! d/dph | |
| 1214 | 8574 | dHprime(m) =(fs(k)%d2H(7)-Hprime(k)*fs(k)%d2pth(7))/fs(k)%dpth(1) ! d/dpph | |
| 1215 | end do | ||
| 1216 | |||
| 1217 | ! evaluate first stage with only r and pph as unknowns | ||
| 1218 | 2858 | jac(1, 1) = fs(1)%dpth(1) | |
| 1219 | 2858 | jac(1, 2) = fs(1)%dpth(4) | |
| 1220 | 2858 | jac(2, 2) = 1d0 | |
| 1221 | |||
| 1222 | jac(1, 1) = jac(1, 1) + si%dt*ahat(1,1) * & | ||
| 1223 | 2858 | (fs(1)%d2H(2) - fs(1)%d2pth(2)*Hprime(1) - fs(1)%dpth(2)*dHprime(1)) | |
| 1224 | jac(1, 2) = jac(1, 2) + si%dt*ahat(1,1)* & | ||
| 1225 | 2858 | (fs(1)%d2H(8) - fs(1)%d2pth(8)*Hprime(1) - fs(1)%dpth(2)*dHprime(2)) | |
| 1226 | jac(2, 1) = jac(2, 1) + si%dt*ahat(1,1) * & | ||
| 1227 | 2858 | (fs(1)%d2H(3) - fs(1)%d2pth(3)*Hprime(1) - fs(1)%dpth(3)*dHprime(1)) | |
| 1228 | jac(2, 2) = jac(2, 2) + si%dt*ahat(1,1) * & | ||
| 1229 | 2858 | (fs(1)%d2H(9) - fs(1)%d2pth(9)*Hprime(1) - fs(1)%dpth(3)*dHprime(2)) | |
| 1230 | |||
| 1231 |
2/2✓ Branch 0 taken 2858 times.
✓ Branch 1 taken 5716 times.
|
8574 | do l = 2, s |
| 1232 | 5716 | n = 4*l-2 | |
| 1233 | jac(1, n-3) = jac(1, n-3) + si%dt*ahat(1,l) * & | ||
| 1234 | 5716 | (fs(l)%d2H(2) - fs(l)%d2pth(2)*Hprime(l) - fs(l)%dpth(2)*dHprime(n-3)) | |
| 1235 | jac(1, n-2) = jac(1, n-2) + si%dt*ahat(1,l) * & | ||
| 1236 | 5716 | (fs(l)%d2H(4) - fs(l)%d2pth(4)*Hprime(l) - fs(l)%dpth(2)*dHprime(n-2)) | |
| 1237 | jac(1, n-1) = jac(1, n-1) + si%dt*ahat(1,l)* & | ||
| 1238 | 5716 | (fs(l)%d2H(5) - fs(l)%d2pth(5)*Hprime(l) - fs(l)%dpth(2)*dHprime(n-1)) | |
| 1239 | jac(1, n) = jac(1, n) + si%dt*ahat(1,l)* & | ||
| 1240 | 5716 | (fs(l)%d2H(8) - fs(l)%d2pth(8)*Hprime(l) - fs(l)%dpth(2)*dHprime(n)) | |
| 1241 | |||
| 1242 | jac(2, n-3) = jac(2, n-3) + si%dt*ahat(1,l) * & | ||
| 1243 | 5716 | (fs(l)%d2H(3) - fs(l)%d2pth(3)*Hprime(l) - fs(l)%dpth(3)*dHprime(n-3)) | |
| 1244 | jac(2, n-2) = jac(2, n-2) + si%dt*ahat(1,l) * & | ||
| 1245 | 5716 | (fs(l)%d2H(5) - fs(l)%d2pth(5)*Hprime(l) - fs(l)%dpth(3)*dHprime(n-2)) | |
| 1246 | jac(2, n-1) = jac(2, n-1) + si%dt*ahat(1,l) * & | ||
| 1247 | 5716 | (fs(l)%d2H(6) - fs(l)%d2pth(6)*Hprime(l) - fs(l)%dpth(3)*dHprime(n-1)) | |
| 1248 | jac(2, n) = jac(2, n) + si%dt*ahat(1,l) * & | ||
| 1249 | 8574 | (fs(l)%d2H(9) - fs(l)%d2pth(9)*Hprime(l) - fs(l)%dpth(3)*dHprime(n)) | |
| 1250 | end do | ||
| 1251 | |||
| 1252 | ! evaluate remaining stages with r, th, ph, pph as unknowns | ||
| 1253 | |||
| 1254 |
2/2✓ Branch 0 taken 2858 times.
✓ Branch 1 taken 5716 times.
|
8574 | do k = 2, s |
| 1255 | 5716 | m = 4*k-2 | |
| 1256 | 5716 | jac(m-3, m-3) = fs(k)%dpth(1) | |
| 1257 | 5716 | jac(m-3, m-2) = fs(k)%dpth(2) | |
| 1258 | 5716 | jac(m-3, m-1) = fs(k)%dpth(3) | |
| 1259 | 5716 | jac(m-3, m) = fs(k)%dpth(4) | |
| 1260 | 5716 | jac(m-2, m-2) = 1d0 | |
| 1261 | 5716 | jac(m-1, m-1) = 1d0 | |
| 1262 | 8574 | jac(m, m) = 1d0 | |
| 1263 | end do | ||
| 1264 | |||
| 1265 | 8574 | l = 1 | |
| 1266 |
2/2✓ Branch 0 taken 2858 times.
✓ Branch 1 taken 5716 times.
|
8574 | do k = 2, s |
| 1267 | 5716 | m = 4*k-2 | |
| 1268 | jac(m-3, 1) = jac(m-3, 1) + si%dt*ahat(k,l) * & ! d/dr | ||
| 1269 | 5716 | (fs(l)%d2H(2) - fs(l)%d2pth(2)*Hprime(l) - fs(l)%dpth(2)*dHprime(1)) | |
| 1270 | jac(m-3, 2) = jac(m-3, 2) + si%dt*ahat(k,l)* & ! d/dpph | ||
| 1271 | 5716 | (fs(l)%d2H(8) - fs(l)%d2pth(8)*Hprime(l) - fs(l)%dpth(2)*dHprime(2)) | |
| 1272 | |||
| 1273 | 5716 | jac(m-2, 1) = jac(m-2, 1) - si%dt*a(k,l)*dHprime(1) ! d/dr | |
| 1274 | 5716 | jac(m-2, 2) = jac(m-2, 2) - si%dt*a(k,l)*dHprime(2) ! d/dpph | |
| 1275 | |||
| 1276 | jac(m-1, 1) = jac(m-1, 1) - si%dt*a(k,l) * & ! d/dr | ||
| 1277 | (-fs(l)%dhph(1)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & | ||
| 1278 | 5716 | +(fs(l)%dvpar(1)-dHprime(1)*fs(l)%hth-Hprime(l)*fs(l)%dhth(1))/fs(l)%hph) | |
| 1279 | jac(m-1, 2) = jac(m-1, 2) - si%dt*a(k,l) * & ! d/dpph | ||
| 1280 | 5716 | (fs(l)%dvpar(4)-dHprime(2)*fs(l)%hth)/fs(l)%hph | |
| 1281 | |||
| 1282 | jac(m, 1) = jac(m, 1) + si%dt*ahat(k,l) * & ! d/dr | ||
| 1283 | 5716 | (fs(l)%d2H(3) - fs(l)%d2pth(3)*Hprime(l) - fs(l)%dpth(3)*dHprime(1)) | |
| 1284 | jac(m, 2) = jac(m, 2) + si%dt*ahat(k,l) * & ! d/dpph | ||
| 1285 | 8574 | (fs(l)%d2H(9) - fs(l)%d2pth(9)*Hprime(l) - fs(l)%dpth(3)*dHprime(2)) | |
| 1286 | end do | ||
| 1287 | |||
| 1288 |
2/2✓ Branch 0 taken 5716 times.
✓ Branch 1 taken 2858 times.
|
8574 | do l = 2, s |
| 1289 |
2/2✓ Branch 0 taken 11432 times.
✓ Branch 1 taken 5716 times.
|
20006 | do k = 2, s |
| 1290 | 11432 | m = 4*k-2 | |
| 1291 | 11432 | n = 4*l-2 | |
| 1292 | jac(m-3, n-3) = jac(m-3, n-3) + si%dt*ahat(k,l) * & ! d/dr | ||
| 1293 | 11432 | (fs(l)%d2H(2) - fs(l)%d2pth(2)*Hprime(l) - fs(l)%dpth(2)*dHprime(n-3)) | |
| 1294 | jac(m-3, n-2) = jac(m-3, n-2) + si%dt*ahat(k,l) * & ! d/dth | ||
| 1295 | 11432 | (fs(l)%d2H(4) - fs(l)%d2pth(4)*Hprime(l) - fs(l)%dpth(2)*dHprime(n-2)) | |
| 1296 | jac(m-3, n-1) = jac(m-3, n-1) + si%dt*ahat(k,l)* & ! d/dph | ||
| 1297 | 11432 | (fs(l)%d2H(5) - fs(l)%d2pth(5)*Hprime(l) - fs(l)%dpth(2)*dHprime(n-1)) | |
| 1298 | jac(m-3, n) = jac(m-3, n) + si%dt*ahat(k,l)* & ! d/dpph | ||
| 1299 | 11432 | (fs(l)%d2H(8) - fs(l)%d2pth(8)*Hprime(l) - fs(l)%dpth(2)*dHprime(n)) | |
| 1300 | |||
| 1301 | 11432 | jac(m-2, n-3) = jac(m-2, n-3) - si%dt*a(k,l)*dHprime(n-3) | |
| 1302 | 11432 | jac(m-2, n-2) = jac(m-2, n-2) - si%dt*a(k,l)*dHprime(n-2) | |
| 1303 | 11432 | jac(m-2, n-1) = jac(m-2, n-1) - si%dt*a(k,l)*dHprime(n-1) | |
| 1304 | 11432 | jac(m-2, n) = jac(m-2, n) - si%dt*a(k,l)*dHprime(n) | |
| 1305 | |||
| 1306 | jac(m-1, n-3) = jac(m-1, n-3) - si%dt*a(k,l) * & ! d/dr | ||
| 1307 | (-fs(l)%dhph(1)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & | ||
| 1308 | 11432 | +(fs(l)%dvpar(1)-dHprime(n-3)*fs(l)%hth-Hprime(l)*fs(l)%dhth(1))/fs(l)%hph) | |
| 1309 | jac(m-1, n-2) = jac(m-1, n-2) - si%dt*a(k,l) * & ! d/dth | ||
| 1310 | (-fs(l)%dhph(2)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & | ||
| 1311 | 11432 | +(fs(l)%dvpar(2)-dHprime(n-2)*fs(l)%hth-Hprime(l)*fs(l)%dhth(2))/fs(l)%hph) | |
| 1312 | jac(m-1, n-1) = jac(m-1, n-1) - si%dt*a(k,l) * & ! d/dph | ||
| 1313 | (-fs(l)%dhph(3)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & | ||
| 1314 | 11432 | +(fs(l)%dvpar(3)-dHprime(n-1)*fs(l)%hth-Hprime(l)*fs(l)%dhth(3))/fs(l)%hph) | |
| 1315 | jac(m-1, n) = jac(m-1, n) - si%dt*a(k,l) * & ! d/dpph | ||
| 1316 | 11432 | (fs(l)%dvpar(4)-dHprime(n)*fs(l)%hth)/fs(l)%hph | |
| 1317 | |||
| 1318 | jac(m, n-3) = jac(m, n-3) + si%dt*ahat(k,l) * & ! d/dr | ||
| 1319 | 11432 | (fs(l)%d2H(3) - fs(l)%d2pth(3)*Hprime(l) - fs(l)%dpth(3)*dHprime(n-3)) | |
| 1320 | jac(m, n-2) = jac(m, n-2) + si%dt*ahat(k,l) * & ! d/dth | ||
| 1321 | 11432 | (fs(l)%d2H(5) - fs(l)%d2pth(5)*Hprime(l) - fs(l)%dpth(3)*dHprime(n-2)) | |
| 1322 | jac(m, n-1) = jac(m, n-1) + si%dt*ahat(k,l) * & ! d/dph | ||
| 1323 | 11432 | (fs(l)%d2H(6) - fs(l)%d2pth(6)*Hprime(l) - fs(l)%dpth(3)*dHprime(n-1)) | |
| 1324 | jac(m, n) = jac(m, n) + si%dt*ahat(k,l) * & ! d/dpph | ||
| 1325 | 17148 | (fs(l)%d2H(9) - fs(l)%d2pth(9)*Hprime(l) - fs(l)%dpth(3)*dHprime(n)) | |
| 1326 | end do | ||
| 1327 | end do | ||
| 1328 | 2858 | end subroutine jac_rk_lobatto | |
| 1329 | |||
| 1330 | 2860 | subroutine guard_lobatto_stage_radii(x, s, status) | |
| 1331 | integer, intent(in) :: s | ||
| 1332 | real(dp), intent(inout) :: x(4*s-2) | ||
| 1333 | integer, intent(out) :: status | ||
| 1334 | |||
| 1335 | integer :: ks, radius_index | ||
| 1336 | |||
| 1337 | 2860 | status = SYMPLECTIC_STEP_OK | |
| 1338 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2860 times.
|
2860 | if (x(1) > sympl_rmax) then |
| 1339 | ✗ | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 1340 | ✗ | return | |
| 1341 | end if | ||
| 1342 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2859 times.
|
2860 | if (x(1) < 0d0) x(1) = 0.01d0 |
| 1343 | |||
| 1344 |
2/2✓ Branch 0 taken 5720 times.
✓ Branch 1 taken 2859 times.
|
8579 | do ks = 2, s |
| 1345 | 5720 | radius_index = 4*ks - 5 | |
| 1346 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 5719 times.
|
5720 | if (x(radius_index) > sympl_rmax) then |
| 1347 | 1 | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 1348 | 1 | return | |
| 1349 | end if | ||
| 1350 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 5718 times.
|
8578 | if (x(radius_index) < 0d0) x(radius_index) = 0.01d0 |
| 1351 | end do | ||
| 1352 | end subroutine guard_lobatto_stage_radii | ||
| 1353 | |||
| 1354 | |||
| 1355 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 212 times.
|
212 | recursive subroutine newton_rk_lobatto(si, fs, s, x, atol, rtol, maxit, xlast, status) |
| 1356 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 1357 | |||
| 1358 | integer, intent(in) :: s | ||
| 1359 | type(field_can_t), intent(inout) :: fs(:) | ||
| 1360 | 424 | integer :: kit, ks, radial_indices(s) | |
| 1361 | |||
| 1362 | real(dp), intent(inout) :: x(4*s-2) | ||
| 1363 | real(dp), intent(in) :: atol, rtol | ||
| 1364 | integer, intent(in) :: maxit | ||
| 1365 | real(dp), intent(out) :: xlast(4*s-2) | ||
| 1366 | integer, intent(out) :: status | ||
| 1367 | |||
| 1368 | 424 | real(dp) :: fvec(4*s-2), fjac(4*s-2, 4*s-2) | |
| 1369 | |||
| 1370 | 424 | real(dp) :: xabs(4*s-2), tolref(4*s-2), fabs(4*s-2), step_scale | |
| 1371 | logical :: boundary_limited, step_limited | ||
| 1372 | |||
| 1373 |
2/2✓ Branch 0 taken 2120 times.
✓ Branch 1 taken 212 times.
|
2332 | xlast = x |
| 1374 | 212 | status = SYMPLECTIC_STEP_MAXITER | |
| 1375 | 212 | boundary_limited = .false. | |
| 1376 | 212 | radial_indices(1) = 1 | |
| 1377 |
2/2✓ Branch 0 taken 424 times.
✓ Branch 1 taken 212 times.
|
636 | do ks = 2, s |
| 1378 | 636 | radial_indices(ks) = 4*ks - 5 | |
| 1379 | end do | ||
| 1380 | |||
| 1381 |
2/2✓ Branch 0 taken 2858 times.
✓ Branch 1 taken 49 times.
|
2907 | do kit = 1, maxit |
| 1382 | |||
| 1383 | 2858 | call guard_lobatto_stage_radii(x, s, status) | |
| 1384 |
1/2✓ Branch 0 taken 2858 times.
✗ Branch 1 not taken.
|
2858 | if (status /= SYMPLECTIC_STEP_OK) return |
| 1385 | 2858 | status = SYMPLECTIC_STEP_MAXITER | |
| 1386 | |||
| 1387 | 2858 | call f_rk_lobatto(si, fs, s, x, fvec, 2) | |
| 1388 | 2858 | call jac_rk_lobatto(si, fs, s, fjac) | |
| 1389 |
2/2✓ Branch 0 taken 28580 times.
✓ Branch 1 taken 2858 times.
|
31438 | fabs = dabs(fvec) |
| 1390 |
2/2✓ Branch 0 taken 28580 times.
✓ Branch 1 taken 2858 times.
|
31438 | xlast = x |
| 1391 | 2858 | call solve_newton_system(fjac, fvec, status) | |
| 1392 |
1/2✓ Branch 0 taken 2858 times.
✗ Branch 1 not taken.
|
2858 | if (status /= SYMPLECTIC_STEP_OK) return |
| 1393 | 2858 | status = SYMPLECTIC_STEP_MAXITER | |
| 1394 | ! after solution: fvec = (xold-xnew)_Newton | ||
| 1395 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2858 times.
|
2858 | if (sympl_rmax > 1d0) then |
| 1396 | ✗ | step_limited = .false. | |
| 1397 | ✗ | x = x - fvec | |
| 1398 | else | ||
| 1399 | call limit_radial_newton_step(x, fvec, radial_indices, step_scale, & | ||
| 1400 | 2858 | step_limited) | |
| 1401 |
4/4✓ Branch 0 taken 990 times.
✓ Branch 1 taken 1868 times.
✓ Branch 2 taken 81 times.
✓ Branch 3 taken 909 times.
|
2858 | boundary_limited = boundary_limited .or. step_limited |
| 1402 |
2/2✓ Branch 0 taken 28580 times.
✓ Branch 1 taken 2858 times.
|
31438 | x = x - step_scale*fvec |
| 1403 | end if | ||
| 1404 |
4/6✓ Branch 0 taken 8574 times.
✓ Branch 1 taken 2858 times.
✓ Branch 2 taken 8574 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2858 times.
|
11432 | if (any(x(radial_indices) > sympl_rmax)) then |
| 1405 | ✗ | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 1406 | return | ||
| 1407 | end if | ||
| 1408 |
2/2✓ Branch 0 taken 28580 times.
✓ Branch 1 taken 2858 times.
|
31438 | xabs = dabs(x - xlast) |
| 1409 | |||
| 1410 | 2858 | tolref(1) = 1d0 | |
| 1411 | 2858 | tolref(2) = dabs(xlast(2)) | |
| 1412 |
2/2✓ Branch 0 taken 5716 times.
✓ Branch 1 taken 2858 times.
|
8574 | do ks = 2, s |
| 1413 | 5716 | tolref(4*ks-2-3) = 1d0 | |
| 1414 | 5716 | tolref(4*ks-2-2) = twopi | |
| 1415 | 5716 | tolref(4*ks-2-1) = twopi | |
| 1416 | 8574 | tolref(4*ks-2) = dabs(xlast(4*ks-2)) | |
| 1417 | end do | ||
| 1418 | |||
| 1419 |
8/8✓ Branch 0 taken 7386 times.
✓ Branch 1 taken 144 times.
✓ Branch 2 taken 4672 times.
✓ Branch 3 taken 2714 times.
✓ Branch 4 taken 1031 times.
✓ Branch 5 taken 1827 times.
✓ Branch 6 taken 68 times.
✓ Branch 7 taken 963 times.
|
7530 | if (.not. step_limited .and. all(fabs < atol)) then |
| 1420 | 68 | status = SYMPLECTIC_STEP_OK | |
| 1421 | 68 | return | |
| 1422 | end if | ||
| 1423 |
8/8✓ Branch 0 taken 16868 times.
✓ Branch 1 taken 1498 times.
✓ Branch 2 taken 15576 times.
✓ Branch 3 taken 1292 times.
✓ Branch 4 taken 963 times.
✓ Branch 5 taken 1827 times.
✓ Branch 6 taken 95 times.
✓ Branch 7 taken 868 times.
|
18415 | if (.not. step_limited .and. all(xabs < rtol*tolref)) then |
| 1424 | 95 | status = SYMPLECTIC_STEP_OK | |
| 1425 | 95 | return | |
| 1426 | end if | ||
| 1427 | enddo | ||
| 1428 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 1 time.
|
49 | if (boundary_limited) then |
| 1429 |
2/4✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
|
48 | if (step_limited .and. & |
| 1430 | radial_boundary_reached(x, radial_indices, rtol)) then | ||
| 1431 | 48 | status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 1432 | else | ||
| 1433 | ✗ | status = SYMPLECTIC_STEP_BOUNDARY_LIMITED | |
| 1434 | end if | ||
| 1435 | else | ||
| 1436 | 1 | call count_event(EVT_RK_LOBATTO_MAXIT) | |
| 1437 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
|
1 | if (accept_warning_maxiter(x, xlast, tolref, rtol, si%warning_factor)) & |
| 1438 | ✗ | status = SYMPLECTIC_STEP_OK | |
| 1439 | end if | ||
| 1440 | end subroutine newton_rk_lobatto | ||
| 1441 | |||
| 1442 | |||
| 1443 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1444 | ! | ||
| 1445 | 380 | recursive subroutine orbit_timestep_sympl_multi(mi, f, ierr) | |
| 1446 | ! | ||
| 1447 | type(multistage_integrator_t), intent(inout) :: mi | ||
| 1448 | type(field_can_t), intent(inout) :: f | ||
| 1449 | |||
| 1450 | integer, intent(out) :: ierr | ||
| 1451 | |||
| 1452 | integer :: kstage | ||
| 1453 | |||
| 1454 | 380 | call orbit_timestep_sympl(mi%stages(1), f, ierr) | |
| 1455 | |||
| 1456 |
2/2✓ Branch 0 taken 3610 times.
✓ Branch 1 taken 380 times.
|
3990 | do kstage = 2, 2*mi%s |
| 1457 |
2/2✓ Branch 0 taken 14440 times.
✓ Branch 1 taken 3610 times.
|
18050 | mi%stages(kstage)%z = mi%stages(kstage-1)%z |
| 1458 | 3610 | mi%stages(kstage)%pthold = f%pth | |
| 1459 | 3990 | call orbit_timestep_sympl(mi%stages(kstage), f, ierr) | |
| 1460 | end do | ||
| 1461 | |||
| 1462 |
2/2✓ Branch 0 taken 1520 times.
✓ Branch 1 taken 380 times.
|
1900 | mi%stages(1)%z = mi%stages(2*mi%s)%z |
| 1463 | 380 | mi%stages(1)%pthold = f%pth | |
| 1464 | |||
| 1465 | 380 | end subroutine orbit_timestep_sympl_multi | |
| 1466 | |||
| 1467 | |||
| 1468 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1469 | ! | ||
| 1470 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 8 times.
|
8 | recursive subroutine orbit_sympl_init_multi(mi, f, z, dtau, ntau, rtol_init, alpha, beta) |
| 1471 | ! | ||
| 1472 | type(multistage_integrator_t), intent(inout) :: mi | ||
| 1473 | type(field_can_t), intent(inout) :: f | ||
| 1474 | |||
| 1475 | real(dp), intent(in) :: z(:) | ||
| 1476 | real(dp), intent(in) :: dtau | ||
| 1477 | integer, intent(in) :: ntau | ||
| 1478 | real(dp), intent(in) :: rtol_init | ||
| 1479 | |||
| 1480 | real(dp), intent(in) :: alpha(:), beta(:) | ||
| 1481 | |||
| 1482 | integer :: ks | ||
| 1483 | |||
| 1484 | 8 | mi%s = size(alpha) | |
| 1485 | |||
| 1486 |
2/2✓ Branch 0 taken 42 times.
✓ Branch 1 taken 8 times.
|
50 | do ks = 1, mi%s |
| 1487 | call orbit_sympl_init(mi%stages(2*ks-1), f, z, & | ||
| 1488 | 42 | alpha(ks)*dtau, ntau, rtol_init, 1) | |
| 1489 | call orbit_sympl_init(mi%stages(2*ks), f, z, & | ||
| 1490 | 50 | beta(ks)*dtau, ntau, rtol_init, 2) | |
| 1491 | end do | ||
| 1492 | 8 | end subroutine orbit_sympl_init_multi | |
| 1493 | |||
| 1494 | |||
| 1495 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1496 | ! | ||
| 1497 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | recursive subroutine orbit_sympl_init_verlet(mi, f, z, dtau, ntau, rtol_init) |
| 1498 | ! | ||
| 1499 | type(multistage_integrator_t), intent(inout) :: mi | ||
| 1500 | |||
| 1501 | type(field_can_t), intent(inout) :: f | ||
| 1502 | |||
| 1503 | real(dp), intent(in) :: z(:) | ||
| 1504 | real(dp), intent(in) :: dtau | ||
| 1505 | integer, intent(in) :: ntau | ||
| 1506 | real(dp), intent(in) :: rtol_init | ||
| 1507 | |||
| 1508 | real(dp) :: alpha(1), beta(1) | ||
| 1509 | |||
| 1510 | 2 | alpha(1) = 0.5d0 | |
| 1511 | 2 | beta(1) = 0.5d0 | |
| 1512 | |||
| 1513 | 2 | call orbit_sympl_init_multi(mi, f, z, dtau, ntau, rtol_init, alpha, beta) | |
| 1514 | 2 | end subroutine orbit_sympl_init_verlet | |
| 1515 | |||
| 1516 | |||
| 1517 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1518 | ! | ||
| 1519 | ✗ | recursive subroutine orbit_sympl_init_order4(mi, f, z, dtau, ntau, rtol_init) | |
| 1520 | ! | ||
| 1521 | ! Composition method of order 4 with s=3 | ||
| 1522 | ! | ||
| 1523 | ! | ||
| 1524 | type(multistage_integrator_t), intent(inout) :: mi | ||
| 1525 | type(field_can_t), intent(inout) :: f | ||
| 1526 | |||
| 1527 | real(dp), intent(in) :: z(:) | ||
| 1528 | real(dp), intent(in) :: dtau | ||
| 1529 | integer, intent(in) :: ntau | ||
| 1530 | real(dp), intent(in) :: rtol_init | ||
| 1531 | |||
| 1532 | real(dp) :: alpha(5), beta(5) | ||
| 1533 | |||
| 1534 | ✗ | alpha(1) = 1d0/(2d0*(2d0 - 2d0**(1d0/3d0))) | |
| 1535 | ✗ | alpha(2) = 2d0**(1d0/3d0)/(2d0*(2d0 - 2d0**(1d0/3d0))) | |
| 1536 | ✗ | alpha(3) = alpha(1) | |
| 1537 | |||
| 1538 | ✗ | beta = alpha | |
| 1539 | |||
| 1540 | ✗ | call orbit_sympl_init_multi(mi, f, z, dtau, ntau, rtol_init, alpha, beta) | |
| 1541 | ✗ | end subroutine orbit_sympl_init_order4 | |
| 1542 | |||
| 1543 | |||
| 1544 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1545 | ! | ||
| 1546 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | recursive subroutine orbit_sympl_init_mclachlan4(mi, f, z, dtau, ntau, rtol_init) |
| 1547 | ! | ||
| 1548 | ! Composition method of order 4 with s=5 by McLachlan (1995) | ||
| 1549 | ! | ||
| 1550 | ! | ||
| 1551 | type(multistage_integrator_t), intent(inout) :: mi | ||
| 1552 | type(field_can_t), intent(inout) :: f | ||
| 1553 | |||
| 1554 | real(dp), intent(in) :: z(:) | ||
| 1555 | real(dp), intent(in) :: dtau | ||
| 1556 | integer, intent(in) :: ntau | ||
| 1557 | real(dp), intent(in) :: rtol_init | ||
| 1558 | |||
| 1559 | real(dp) :: alpha(5), beta(5) | ||
| 1560 | |||
| 1561 | 2 | alpha(1) = (146d0 + 5d0*dsqrt(19d0))/540d0 | |
| 1562 | 2 | alpha(2) = (-2d0 + 10d0*dsqrt(19d0))/135d0 | |
| 1563 | 2 | alpha(3) = 1d0/5d0 | |
| 1564 | 2 | alpha(4) = (-23d0 + 20d0*dsqrt(19d0))/270d0 | |
| 1565 | 2 | alpha(5) = (14d0 - dsqrt(19d0))/108d0 | |
| 1566 | |||
| 1567 | 2 | beta(5) = alpha(1) | |
| 1568 | 2 | beta(4) = alpha(2) | |
| 1569 | 2 | beta(3) = alpha(3) | |
| 1570 | 2 | beta(2) = alpha(4) | |
| 1571 | 2 | beta(1) = alpha(5) | |
| 1572 | |||
| 1573 | 2 | call orbit_sympl_init_multi(mi, f, z, dtau, ntau, rtol_init, alpha, beta) | |
| 1574 | 2 | end subroutine orbit_sympl_init_mclachlan4 | |
| 1575 | |||
| 1576 | |||
| 1577 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1578 | ! | ||
| 1579 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | recursive subroutine orbit_sympl_init_blanes4(mi, f, z, dtau, ntau, rtol_init) |
| 1580 | ! | ||
| 1581 | ! Composition method of order 4 with s=6 by Blanes&Moan (2002) | ||
| 1582 | ! with coefficients in the form of Hairer (2002) | ||
| 1583 | ! | ||
| 1584 | ! | ||
| 1585 | type(multistage_integrator_t), intent(inout) :: mi | ||
| 1586 | type(field_can_t), intent(inout) :: f | ||
| 1587 | |||
| 1588 | real(dp), intent(in) :: z(:) | ||
| 1589 | real(dp), intent(in) :: dtau | ||
| 1590 | integer, intent(in) :: ntau | ||
| 1591 | real(dp), intent(in) :: rtol_init | ||
| 1592 | |||
| 1593 | real(dp) :: alpha(6), beta(6) | ||
| 1594 | |||
| 1595 | 2 | alpha(1) = 0.16231455076687d0 | |
| 1596 | 2 | alpha(2) = 0.37087741497958d0 | |
| 1597 | 2 | alpha(3) = 0.059762097006575d0 | |
| 1598 | 2 | alpha(4) = 0.40993371990193d0 | |
| 1599 | 2 | alpha(5) = 0.23399525073150d0 | |
| 1600 | 2 | alpha(6) = 0.082984406417405d0 | |
| 1601 | |||
| 1602 | 2 | beta(6) = alpha(1) | |
| 1603 | 2 | beta(5) = alpha(2) | |
| 1604 | 2 | beta(4) = alpha(3) | |
| 1605 | 2 | beta(3) = alpha(4) | |
| 1606 | 2 | beta(2) = alpha(5) | |
| 1607 | 2 | beta(1) = alpha(6) | |
| 1608 | |||
| 1609 | 2 | call orbit_sympl_init_multi(mi, f, z, dtau, ntau, rtol_init, alpha, beta) | |
| 1610 | 2 | end subroutine orbit_sympl_init_blanes4 | |
| 1611 | |||
| 1612 | |||
| 1613 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1614 | ! | ||
| 1615 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | recursive subroutine orbit_sympl_init_kahan6(mi, f, z, dtau, ntau, rtol_init) |
| 1616 | ! | ||
| 1617 | ! Composition method of order 6 with s=9 by Kahan&Li (1995) | ||
| 1618 | ! with coefficients in the form of Hairer (2002) | ||
| 1619 | ! | ||
| 1620 | ! | ||
| 1621 | type(multistage_integrator_t), intent(inout) :: mi | ||
| 1622 | type(field_can_t), intent(inout) :: f | ||
| 1623 | |||
| 1624 | real(dp), intent(in) :: z(:) | ||
| 1625 | real(dp), intent(in) :: dtau | ||
| 1626 | integer, intent(in) :: ntau | ||
| 1627 | real(dp), intent(in) :: rtol_init | ||
| 1628 | |||
| 1629 | real(dp) :: gam(9) | ||
| 1630 | |||
| 1631 | 2 | gam(1) = 0.39216144400731413927925056d0 | |
| 1632 | 2 | gam(2) = 0.33259913678935943859974864d0 | |
| 1633 | 2 | gam(3) = -0.70624617255763935980996482d0 | |
| 1634 | 2 | gam(4) = 0.08221359629355080023149045d0 | |
| 1635 | 2 | gam(5) = 0.79854399093482996339895035d0 | |
| 1636 | 2 | gam(6) = gam(4) | |
| 1637 | 2 | gam(7) = gam(3) | |
| 1638 | 2 | gam(8) = gam(2) | |
| 1639 | 2 | gam(9) = gam(1) | |
| 1640 | |||
| 1641 |
4/4✓ Branch 0 taken 18 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 18 times.
✓ Branch 3 taken 2 times.
|
38 | call orbit_sympl_init_multi(mi,f,z,dtau,ntau,rtol_init,gam/2d0,gam/2d0) |
| 1642 | 2 | end subroutine orbit_sympl_init_kahan6 | |
| 1643 | |||
| 1644 | |||
| 1645 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1646 | ! | ||
| 1647 | ✗ | recursive subroutine orbit_sympl_init_kahan8(mi, f, z, dtau, ntau, rtol_init) | |
| 1648 | ! | ||
| 1649 | ! Composition method of order 8 with s=17 by Kahan&Li (1995) | ||
| 1650 | ! with coefficients in the form of Hairer (2002) | ||
| 1651 | ! | ||
| 1652 | ! | ||
| 1653 | type(multistage_integrator_t), intent(inout) :: mi | ||
| 1654 | type(field_can_t), intent(inout) :: f | ||
| 1655 | |||
| 1656 | real(dp), intent(in) :: z(:) | ||
| 1657 | real(dp), intent(in) :: dtau | ||
| 1658 | integer, intent(in) :: ntau | ||
| 1659 | real(dp), intent(in) :: rtol_init | ||
| 1660 | |||
| 1661 | real(dp) :: gam(17) | ||
| 1662 | |||
| 1663 | ✗ | gam(1) = 0.13020248308889008087881763d0 | |
| 1664 | ✗ | gam(2) = 0.56116298177510838456196441d0 | |
| 1665 | ✗ | gam(3) = -0.38947496264484728640807860d0 | |
| 1666 | ✗ | gam(4) = 0.15884190655515560089621075d0 | |
| 1667 | ✗ | gam(5) = -0.39590389413323757733623154d0 | |
| 1668 | ✗ | gam(6) = 0.18453964097831570709183254d0 | |
| 1669 | ✗ | gam(7) = 0.25837438768632204729397911d0 | |
| 1670 | ✗ | gam(8) = 0.29501172360931029887096624d0 | |
| 1671 | ✗ | gam(9) = -0.60550853383003451169892108d0 | |
| 1672 | ✗ | gam(10) = gam(8) | |
| 1673 | ✗ | gam(11) = gam(7) | |
| 1674 | ✗ | gam(12) = gam(6) | |
| 1675 | ✗ | gam(13) = gam(5) | |
| 1676 | ✗ | gam(14) = gam(4) | |
| 1677 | ✗ | gam(15) = gam(3) | |
| 1678 | ✗ | gam(16) = gam(2) | |
| 1679 | ✗ | gam(17) = gam(1) | |
| 1680 | |||
| 1681 | ✗ | call orbit_sympl_init_multi(mi,f,z,dtau,ntau,rtol_init,gam/2d0,gam/2d0) | |
| 1682 | ✗ | end subroutine orbit_sympl_init_kahan8 | |
| 1683 | |||
| 1684 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1685 | ! | ||
| 1686 | 4403570 | recursive subroutine orbit_timestep_sympl_expl_impl_euler(si, f, ierr) | |
| 1687 | ! | ||
| 1688 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 1689 | type(field_can_t), intent(inout) :: f | ||
| 1690 | integer, intent(out) :: ierr | ||
| 1691 | |||
| 1692 | integer, parameter :: n = 2 | ||
| 1693 | integer, parameter :: maxit = 32 | ||
| 1694 | |||
| 1695 | real(dp), dimension(n) :: x, xlast | ||
| 1696 | integer :: ktau, newton_status | ||
| 1697 | logical :: crossed | ||
| 1698 | type(symplectic_integrator_t) :: accepted_integrator | ||
| 1699 | type(field_can_t) :: accepted_field | ||
| 1700 | |||
| 1701 | 4403570 | ierr = 0 | |
| 1702 | 4403570 | ktau = 0 | |
| 1703 |
2/2✓ Branch 0 taken 4403570 times.
✓ Branch 1 taken 4403440 times.
|
8807010 | do while(ktau .lt. si%ntau) |
| 1704 | 4403570 | accepted_integrator = si | |
| 1705 | 4403570 | accepted_field = f | |
| 1706 | 4403570 | si%pthold = f%pth | |
| 1707 | |||
| 1708 | 4403570 | x(1)=si%z(1) | |
| 1709 | 4403570 | x(2)=si%z(4) | |
| 1710 | |||
| 1711 | 4403570 | call newton1(si, f, x, maxit, xlast, newton_status) | |
| 1712 |
2/2✓ Branch 0 taken 130 times.
✓ Branch 1 taken 4403440 times.
|
4403570 | if (newton_status /= SYMPLECTIC_STEP_OK) then |
| 1713 | 130 | si = accepted_integrator | |
| 1714 | 130 | f = accepted_field | |
| 1715 | 130 | ierr = newton_status | |
| 1716 | 130 | return | |
| 1717 | end if | ||
| 1718 | |||
| 1719 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4403440 times.
|
4403440 | if (x(1) > sympl_rmax) then |
| 1720 | ✗ | si = accepted_integrator | |
| 1721 | ✗ | f = accepted_field | |
| 1722 | ✗ | ierr = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 1723 | ✗ | return | |
| 1724 | end if | ||
| 1725 | |||
| 1726 | 4403440 | call apply_axis_chart_switch(x(1), si%z(2), crossed, ierr) | |
| 1727 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4403440 times.
|
4403440 | if (ierr /= SYMPLECTIC_STEP_OK) then |
| 1728 | ✗ | si = accepted_integrator | |
| 1729 | ✗ | f = accepted_field | |
| 1730 | ✗ | return | |
| 1731 | end if | ||
| 1732 | |||
| 1733 | 4403440 | si%z(1) = x(1) | |
| 1734 | 4403440 | si%z(4) = x(2) | |
| 1735 | |||
| 1736 |
1/2✓ Branch 0 taken 4403440 times.
✗ Branch 1 not taken.
|
4403440 | if (extrap_field .and. .not. crossed) then |
| 1737 | 4403440 | call sympl_euler1_extrapolate_field(si, f, x, xlast) | |
| 1738 | else | ||
| 1739 | ! After a chart switch xlast lives in the other chart; extrapolation | ||
| 1740 | ! across the flip is invalid, evaluate the field fresh instead. | ||
| 1741 | ✗ | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | |
| 1742 | ✗ | call get_derivatives(f, si%z(4)) | |
| 1743 | endif | ||
| 1744 | |||
| 1745 | 4403440 | call sympl_euler1_advance_angles(si, f) | |
| 1746 | |||
| 1747 | 4403440 | ktau = ktau+1 | |
| 1748 | enddo | ||
| 1749 | |||
| 1750 | 4403570 | end subroutine orbit_timestep_sympl_expl_impl_euler | |
| 1751 | |||
| 1752 | |||
| 1753 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1754 | ! | ||
| 1755 | 4363 | recursive subroutine orbit_timestep_sympl_impl_expl_euler(si, f, ierr) | |
| 1756 | ! | ||
| 1757 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 1758 | type(field_can_t), intent(inout) :: f | ||
| 1759 | |||
| 1760 | integer, intent(out) :: ierr | ||
| 1761 | |||
| 1762 | integer, parameter :: n = 3 | ||
| 1763 | integer, parameter :: maxit = 32 | ||
| 1764 | |||
| 1765 | real(dp), dimension(n) :: x, xlast, dz | ||
| 1766 | integer :: ktau, newton_status | ||
| 1767 | logical :: crossed | ||
| 1768 | type(symplectic_integrator_t) :: accepted_integrator | ||
| 1769 | type(field_can_t) :: accepted_field | ||
| 1770 | |||
| 1771 | 4363 | ierr = 0 | |
| 1772 | 4363 | ktau = 0 | |
| 1773 |
2/2✓ Branch 0 taken 4363 times.
✓ Branch 1 taken 4362 times.
|
8725 | do while(ktau .lt. si%ntau) |
| 1774 | 4363 | accepted_integrator = si | |
| 1775 | 4363 | accepted_field = f | |
| 1776 | 4363 | si%pthold = f%pth | |
| 1777 | |||
| 1778 |
2/2✓ Branch 0 taken 13089 times.
✓ Branch 1 taken 4363 times.
|
17452 | x = si%z(1:3) |
| 1779 | |||
| 1780 | 4363 | call newton2(si, f, x, si%atol, si%rtol, maxit, xlast, newton_status) | |
| 1781 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 4362 times.
|
4363 | if (newton_status /= SYMPLECTIC_STEP_OK) then |
| 1782 | 1 | si = accepted_integrator | |
| 1783 | 1 | f = accepted_field | |
| 1784 | 1 | ierr = newton_status | |
| 1785 | 1 | return | |
| 1786 | end if | ||
| 1787 | |||
| 1788 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4362 times.
|
4362 | if (x(1) > sympl_rmax) then |
| 1789 | ✗ | si = accepted_integrator | |
| 1790 | ✗ | f = accepted_field | |
| 1791 | ✗ | ierr = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 1792 | ✗ | return | |
| 1793 | end if | ||
| 1794 | |||
| 1795 | 4362 | call apply_axis_chart_switch(x(1), x(2), crossed, ierr) | |
| 1796 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4362 times.
|
4362 | if (ierr /= SYMPLECTIC_STEP_OK) then |
| 1797 | ✗ | si = accepted_integrator | |
| 1798 | ✗ | f = accepted_field | |
| 1799 | ✗ | return | |
| 1800 | end if | ||
| 1801 | |||
| 1802 |
2/2✓ Branch 0 taken 13086 times.
✓ Branch 1 taken 4362 times.
|
17448 | si%z(1:3) = x |
| 1803 | |||
| 1804 |
1/2✓ Branch 0 taken 4362 times.
✗ Branch 1 not taken.
|
4362 | if (extrap_field .and. .not. crossed) then |
| 1805 | 4362 | dz(1) = x(1)-xlast(1) | |
| 1806 | 4362 | dz(2) = x(2)-xlast(2) | |
| 1807 | 4362 | dz(3) = x(3)-xlast(3) | |
| 1808 | |||
| 1809 | 4362 | f%dH(1) = f%dH(1) + f%d2H(1)*dz(1) + f%d2H(2)*dz(2) + f%d2H(3)*dz(3) | |
| 1810 | 4362 | f%dH(2) = f%dH(2) + f%d2H(2)*dz(1) + f%d2H(4)*dz(2) + f%d2H(5)*dz(3) | |
| 1811 | 4362 | f%dH(3) = f%dH(3) + f%d2H(3)*dz(1) + f%d2H(5)*dz(2) + f%d2H(6)*dz(3) | |
| 1812 | |||
| 1813 | 4362 | f%dpth(1) = f%dpth(1)+f%d2pth(1)*dz(1)+f%d2pth(2)*dz(2)+f%d2pth(3)*dz(3) | |
| 1814 | 4362 | f%dpth(2) = f%dpth(2)+f%d2pth(2)*dz(1)+f%d2pth(4)*dz(2)+f%d2pth(5)*dz(3) | |
| 1815 | 4362 | f%dpth(3) = f%dpth(3)+f%d2pth(3)*dz(1)+f%d2pth(5)*dz(2)+f%d2pth(6)*dz(3) | |
| 1816 | else | ||
| 1817 | ✗ | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | |
| 1818 | ✗ | call get_derivatives(f, si%z(4)) | |
| 1819 | endif | ||
| 1820 | |||
| 1821 | 4362 | f%pth = si%pthold - si%dt*(f%dH(2) - f%dH(1)*f%dpth(2)/f%dpth(1)) | |
| 1822 | 4362 | si%z(4) = si%z(4) - si%dt*(f%dH(3) - f%dH(1)*f%dpth(3)/f%dpth(1)) | |
| 1823 | |||
| 1824 | 4362 | ktau = ktau+1 | |
| 1825 | enddo | ||
| 1826 | |||
| 1827 | 4363 | end subroutine orbit_timestep_sympl_impl_expl_euler | |
| 1828 | |||
| 1829 | |||
| 1830 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1831 | ! | ||
| 1832 | 642594 | subroutine attempt_midpoint_step(si, f, status) | |
| 1833 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 1834 | type(field_can_t), intent(inout) :: f | ||
| 1835 | integer, intent(out) :: status | ||
| 1836 | |||
| 1837 | integer, parameter :: n = 5, maxit = 32 | ||
| 1838 | real(dp) :: x(n), xlast(n) | ||
| 1839 | logical :: crossed | ||
| 1840 | |||
| 1841 | ! 32 as in the other schemes: the atol/rtol exits keep converged steps cheap, | ||
| 1842 | ! and the exact-landing substep solve (#441) needs fully converged trial | ||
| 1843 | ! steps, which 8 iterations do not guarantee near interfaces. | ||
| 1844 | 321340 | si%pthold = f%pth | |
| 1845 |
2/2✓ Branch 0 taken 1285360 times.
✓ Branch 1 taken 321340 times.
|
1606700 | x(1:4) = si%z |
| 1846 | 321340 | x(5) = si%z(1) | |
| 1847 | |||
| 1848 | 321340 | call newton_midpoint(si, f, x, si%atol, si%rtol, maxit, xlast, status) | |
| 1849 |
2/2✓ Branch 0 taken 321254 times.
✓ Branch 1 taken 86 times.
|
321340 | if (status /= SYMPLECTIC_STEP_OK) return |
| 1850 | |||
| 1851 | 321254 | call apply_axis_chart_switch(x(1), x(2), crossed, status) | |
| 1852 |
1/2✓ Branch 0 taken 321254 times.
✗ Branch 1 not taken.
|
321254 | if (status /= SYMPLECTIC_STEP_OK) return |
| 1853 | |||
| 1854 |
2/2✓ Branch 0 taken 1285016 times.
✓ Branch 1 taken 321254 times.
|
1606270 | si%z = x(1:4) |
| 1855 |
2/2✓ Branch 0 taken 321253 times.
✓ Branch 1 taken 1 time.
|
321254 | if (extrap_field .and. .not. crossed) then |
| 1856 | ! The carried momentum uses endpoint coordinates, not x(5). | ||
| 1857 | ! Old-chart increments cannot extrapolate p_theta after a switch. | ||
| 1858 | f%pth = f%pth + f%dpth(1)*(x(1)-xlast(1)) & | ||
| 1859 | + f%dpth(2)*(x(2)-xlast(2)) + f%dpth(3)*(x(3)-xlast(3)) & | ||
| 1860 | 321253 | + f%dpth(4)*(x(4)-xlast(4)) | |
| 1861 | else | ||
| 1862 | 1 | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | |
| 1863 | 1 | call get_val(f, si%z(4)) | |
| 1864 | end if | ||
| 1865 | end subroutine attempt_midpoint_step | ||
| 1866 | |||
| 1867 | 65 | subroutine get_boundary_event_tolerances(rtol, fraction_tolerance, radial_tolerance) | |
| 1868 | real(dp), intent(in) :: rtol | ||
| 1869 | real(dp), intent(out) :: fraction_tolerance, radial_tolerance | ||
| 1870 | |||
| 1871 | 65 | fraction_tolerance = max(1d-12, 10d0*rtol) | |
| 1872 | 65 | radial_tolerance = max(1d-10, 10d0*rtol) | |
| 1873 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 62 times.
|
65 | if (boundary_event_fraction_tolerance > 0d0) then |
| 1874 | 3 | fraction_tolerance = boundary_event_fraction_tolerance | |
| 1875 | end if | ||
| 1876 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 62 times.
|
65 | if (boundary_event_radial_tolerance > 0d0) then |
| 1877 | 3 | radial_tolerance = boundary_event_radial_tolerance | |
| 1878 | end if | ||
| 1879 | 65 | end subroutine get_boundary_event_tolerances | |
| 1880 | |||
| 1881 |
4/4✓ Branch 0 taken 192 times.
✓ Branch 1 taken 64 times.
✓ Branch 2 taken 384 times.
✓ Branch 3 taken 64 times.
|
640 | subroutine locate_symplectic_boundary(accepted_integrator, accepted_field, & |
| 1882 | raw_step, si, f, status, event_fraction, radial_residual, fraction_width) | ||
| 1883 | type(symplectic_integrator_t), intent(in) :: accepted_integrator | ||
| 1884 | type(field_can_t), intent(in) :: accepted_field | ||
| 1885 | procedure(orbit_timestep_sympl_i) :: raw_step | ||
| 1886 | type(symplectic_integrator_t), intent(out) :: si | ||
| 1887 | type(field_can_t), intent(out) :: f | ||
| 1888 | integer, intent(out) :: status | ||
| 1889 | real(dp), intent(out) :: event_fraction, radial_residual, fraction_width | ||
| 1890 | |||
| 1891 | integer, parameter :: max_event_iterations = 64 | ||
| 1892 | type(symplectic_integrator_t) :: trial_integrator, best_integrator | ||
| 1893 | type(field_can_t) :: trial_field, best_field | ||
| 1894 | real(dp) :: fraction_low, fraction_ceiling, trial_fraction | ||
| 1895 | real(dp) :: previous_fraction, previous_radius | ||
| 1896 | real(dp) :: root_estimate | ||
| 1897 | real(dp) :: fraction_tolerance, radial_tolerance | ||
| 1898 | integer :: iteration, trial_status | ||
| 1899 | logical :: root_estimate_available | ||
| 1900 | |||
| 1901 | 64 | fraction_low = 0d0 | |
| 1902 | 64 | fraction_ceiling = 1d0 | |
| 1903 | 64 | trial_fraction = 0.5d0 | |
| 1904 | 64 | root_estimate_available = .false. | |
| 1905 | 64 | best_integrator = accepted_integrator | |
| 1906 | 64 | best_field = accepted_field | |
| 1907 | call get_boundary_event_tolerances(accepted_integrator%rtol, & | ||
| 1908 | 64 | fraction_tolerance, radial_tolerance) | |
| 1909 | |||
| 1910 |
2/2✓ Branch 1 taken 1552 times.
✓ Branch 2 taken 1 time.
|
1553 | do iteration = 1, max_event_iterations |
| 1911 | 1552 | trial_integrator = accepted_integrator | |
| 1912 | 1552 | trial_field = accepted_field | |
| 1913 | 1552 | trial_integrator%dt = accepted_integrator%dt*trial_fraction | |
| 1914 | 1552 | trial_integrator%ntau = 1 | |
| 1915 | 1552 | call raw_step(trial_integrator, trial_field, trial_status) | |
| 1916 | |||
| 1917 | 1238 | select case (trial_status) | |
| 1918 | case (SYMPLECTIC_STEP_OK) | ||
| 1919 | 1238 | previous_fraction = fraction_low | |
| 1920 | 1238 | previous_radius = best_integrator%z(1) | |
| 1921 | 1238 | fraction_low = trial_fraction | |
| 1922 | 1238 | best_integrator = trial_integrator | |
| 1923 | 1238 | best_field = trial_field | |
| 1924 | 1238 | radial_residual = abs(sympl_rmax - best_integrator%z(1)) | |
| 1925 | |||
| 1926 |
2/2✓ Branch 0 taken 1206 times.
✓ Branch 1 taken 32 times.
|
1238 | if (radial_residual == 0d0) then |
| 1927 | fraction_ceiling = fraction_low | ||
| 1928 | root_estimate = fraction_low | ||
| 1929 | root_estimate_available = .true. | ||
| 1930 |
3/4✓ Branch 0 taken 1169 times.
✓ Branch 1 taken 37 times.
✓ Branch 2 taken 1169 times.
✗ Branch 3 not taken.
|
1206 | else if (best_integrator%z(1) > previous_radius .and. & |
| 1931 | fraction_low > previous_fraction) then | ||
| 1932 | root_estimate = fraction_low + radial_residual* & | ||
| 1933 | (fraction_low - previous_fraction)/ & | ||
| 1934 | 1169 | (best_integrator%z(1) - previous_radius) | |
| 1935 | 1169 | root_estimate = min(fraction_ceiling, max(fraction_low, root_estimate)) | |
| 1936 | 1169 | root_estimate_available = .true. | |
| 1937 | else | ||
| 1938 | root_estimate_available = .false. | ||
| 1939 | end if | ||
| 1940 | 1238 | fraction_width = fraction_ceiling - fraction_low | |
| 1941 | |||
| 1942 |
4/4✓ Branch 0 taken 1201 times.
✓ Branch 1 taken 37 times.
✓ Branch 2 taken 59 times.
✓ Branch 3 taken 1142 times.
|
1238 | if (root_estimate_available .and. & |
| 1943 | boundary_event_converged(fraction_width, radial_residual, & | ||
| 1944 | fraction_tolerance, radial_tolerance)) then | ||
| 1945 | 59 | si = best_integrator | |
| 1946 | 59 | f = best_field | |
| 1947 | 59 | si%dt = accepted_integrator%dt | |
| 1948 | 59 | si%ntau = accepted_integrator%ntau | |
| 1949 | 59 | si%z(1) = sympl_rmax | |
| 1950 | 59 | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | |
| 1951 | 59 | call get_val(f, si%z(4)) | |
| 1952 | 59 | event_fraction = fraction_low | |
| 1953 | 59 | status = SYMPLECTIC_STEP_BOUNDARY | |
| 1954 | 60 | return | |
| 1955 | end if | ||
| 1956 | |||
| 1957 |
1/2✓ Branch 0 taken 1142 times.
✗ Branch 1 not taken.
|
1142 | if (root_estimate_available .and. root_estimate > fraction_low) then |
| 1958 | trial_fraction = fraction_low + 0.8d0* & | ||
| 1959 | 1142 | (root_estimate - fraction_low) | |
| 1960 | else | ||
| 1961 | 37 | trial_fraction = 0.5d0*(fraction_low + fraction_ceiling) | |
| 1962 | end if | ||
| 1963 | case (SYMPLECTIC_STEP_OUTSIDE_DOMAIN) | ||
| 1964 | 313 | fraction_ceiling = trial_fraction | |
| 1965 | 313 | trial_fraction = 0.5d0*(fraction_low + fraction_ceiling) | |
| 1966 | case (SYMPLECTIC_STEP_BOUNDARY_LIMITED) | ||
| 1967 | 1 | si = accepted_integrator | |
| 1968 | 1 | f = accepted_field | |
| 1969 | 1 | status = SYMPLECTIC_STEP_EVENT_NOT_CONVERGED | |
| 1970 | 1 | event_fraction = 0d0 | |
| 1971 | 1 | radial_residual = abs(sympl_rmax - best_integrator%z(1)) | |
| 1972 | 1 | fraction_width = fraction_ceiling - fraction_low | |
| 1973 | 1 | return | |
| 1974 | case default | ||
| 1975 | ✗ | si = accepted_integrator | |
| 1976 | ✗ | f = accepted_field | |
| 1977 | ✗ | status = trial_status | |
| 1978 | ✗ | event_fraction = 0d0 | |
| 1979 | ✗ | radial_residual = abs(sympl_rmax - best_integrator%z(1)) | |
| 1980 | ✗ | fraction_width = fraction_ceiling - fraction_low | |
| 1981 |
3/4✓ Branch 0 taken 1238 times.
✓ Branch 1 taken 313 times.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
|
1552 | return |
| 1982 | end select | ||
| 1983 | |||
| 1984 |
3/4✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1489 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1489 times.
|
1492 | if (trial_fraction <= fraction_low .or. & |
| 1985 | 1556 | trial_fraction >= fraction_ceiling) then | |
| 1986 | exit | ||
| 1987 | end if | ||
| 1988 | end do | ||
| 1989 | |||
| 1990 | 4 | si = accepted_integrator | |
| 1991 | 4 | f = accepted_field | |
| 1992 | 4 | status = SYMPLECTIC_STEP_EVENT_NOT_CONVERGED | |
| 1993 | 4 | event_fraction = 0d0 | |
| 1994 | 4 | radial_residual = abs(sympl_rmax - best_integrator%z(1)) | |
| 1995 | 4 | fraction_width = fraction_ceiling - fraction_low | |
| 1996 |
4/4✓ Branch 0 taken 192 times.
✓ Branch 1 taken 64 times.
✓ Branch 2 taken 384 times.
✓ Branch 3 taken 64 times.
|
640 | end subroutine locate_symplectic_boundary |
| 1997 | |||
| 1998 | 1147401 | recursive subroutine orbit_timestep_sympl_midpoint(si, f, ierr) | |
| 1999 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 2000 | type(field_can_t), intent(inout) :: f | ||
| 2001 | integer, intent(out) :: ierr | ||
| 2002 | |||
| 2003 | integer, parameter :: n = 5, maxit = 32 | ||
| 2004 | real(dp) :: x(n), xlast(n) | ||
| 2005 | integer :: ktau, step_status | ||
| 2006 | type(symplectic_integrator_t) :: accepted_integrator | ||
| 2007 | type(field_can_t) :: accepted_field | ||
| 2008 | logical :: crossed | ||
| 2009 | |||
| 2010 | 1147401 | if (sympl_rmax > 1d0) then | |
| 2011 | 826063 | ierr = SYMPLECTIC_STEP_OK | |
| 2012 |
2/2✓ Branch 0 taken 826063 times.
✓ Branch 1 taken 826063 times.
|
1652126 | do ktau = 1, si%ntau |
| 2013 | 826063 | accepted_integrator = si | |
| 2014 | 826063 | accepted_field = f | |
| 2015 | 826063 | si%pthold = f%pth | |
| 2016 |
2/2✓ Branch 0 taken 3304252 times.
✓ Branch 1 taken 826063 times.
|
4130315 | x(1:4) = si%z |
| 2017 | 826063 | x(5) = si%z(1) | |
| 2018 | call newton_midpoint(si, f, x, si%atol, si%rtol, maxit, xlast, & | ||
| 2019 | 826063 | step_status) | |
| 2020 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 826063 times.
|
826063 | if (step_status /= SYMPLECTIC_STEP_OK) then |
| 2021 | ✗ | si = accepted_integrator | |
| 2022 | ✗ | f = accepted_field | |
| 2023 | ✗ | ierr = step_status | |
| 2024 | 826063 | return | |
| 2025 | end if | ||
| 2026 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 826063 times.
|
826063 | if (x(1) > sympl_rmax) then |
| 2027 | ✗ | ierr = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 2028 | ✗ | return | |
| 2029 | end if | ||
| 2030 | 826063 | call apply_axis_chart_switch(x(1), x(2), crossed, ierr) | |
| 2031 |
1/2✓ Branch 0 taken 826063 times.
✗ Branch 1 not taken.
|
826063 | if (ierr /= SYMPLECTIC_STEP_OK) return |
| 2032 |
2/2✓ Branch 0 taken 3304252 times.
✓ Branch 1 taken 826063 times.
|
4130315 | si%z = x(1:4) |
| 2033 |
2/2✓ Branch 1 taken 826062 times.
✓ Branch 2 taken 1 time.
|
2478189 | if (extrap_field .and. .not. crossed) then |
| 2034 | ! Only the endpoint coordinates enter the carried momentum. | ||
| 2035 | ! Across a chart switch, evaluate p_theta at the new endpoint. | ||
| 2036 | f%pth = f%pth + f%dpth(1)*(x(1)-xlast(1)) & | ||
| 2037 | + f%dpth(2)*(x(2)-xlast(2)) + f%dpth(3)*(x(3)-xlast(3)) & | ||
| 2038 | 826062 | + f%dpth(4)*(x(4)-xlast(4)) | |
| 2039 | else | ||
| 2040 | 1 | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | |
| 2041 | 1 | call get_val(f, si%z(4)) | |
| 2042 | end if | ||
| 2043 | end do | ||
| 2044 | return | ||
| 2045 | end if | ||
| 2046 | |||
| 2047 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 321338 times.
|
321338 | if (si%ntau /= 1) error stop & |
| 2048 | ✗ | 'orbit_timestep_sympl_midpoint raw step requires ntau=1' | |
| 2049 | 321338 | call attempt_midpoint_step(si, f, ierr) | |
| 2050 |
2/2✓ Branch 0 taken 826063 times.
✓ Branch 1 taken 321338 times.
|
1147401 | end subroutine orbit_timestep_sympl_midpoint |
| 2051 | |||
| 2052 | 4391354 | recursive subroutine orbit_timestep_sympl_with_boundary(si, f, ierr) | |
| 2053 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 2054 | type(field_can_t), intent(inout) :: f | ||
| 2055 | integer, intent(out) :: ierr | ||
| 2056 | |||
| 2057 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4391354 times.
|
4391354 | if (.not. associated(raw_timestep_sympl)) then |
| 2058 | ✗ | error stop 'symplectic timestep implementation is not initialized' | |
| 2059 | end if | ||
| 2060 | 4391354 | call advance_symplectic_with_boundary(si, f, raw_timestep_sympl, ierr) | |
| 2061 | 4391354 | end subroutine orbit_timestep_sympl_with_boundary | |
| 2062 | |||
| 2063 | 4391357 | subroutine advance_symplectic_with_boundary(si, f, raw_step, ierr) | |
| 2064 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 2065 | type(field_can_t), intent(inout) :: f | ||
| 2066 | procedure(orbit_timestep_sympl_i) :: raw_step | ||
| 2067 | integer, intent(out) :: ierr | ||
| 2068 | |||
| 2069 | type(symplectic_integrator_t) :: entry_integrator, accepted_integrator | ||
| 2070 | type(field_can_t) :: entry_field, accepted_field | ||
| 2071 | real(dp) :: event_fraction, radial_residual, fraction_width | ||
| 2072 | integer :: ktau, step_status, original_ntau | ||
| 2073 | |||
| 2074 | 4391357 | entry_integrator = si | |
| 2075 | 4391357 | entry_field = f | |
| 2076 | 4391357 | original_ntau = si%ntau | |
| 2077 | 4391357 | si%last_step_fraction = 0d0 | |
| 2078 | 4391357 | si%last_event_radial_residual = 0d0 | |
| 2079 | 4391357 | si%last_event_fraction_width = 0d0 | |
| 2080 | 4391357 | radial_residual = 0d0 | |
| 2081 | 4391357 | fraction_width = 0d0 | |
| 2082 | 4391357 | ierr = SYMPLECTIC_STEP_OK | |
| 2083 | |||
| 2084 | 4391357 | if (sympl_rmax > 1d0) then | |
| 2085 | ✗ | do ktau = 1, original_ntau | |
| 2086 | ✗ | accepted_integrator = si | |
| 2087 | ✗ | accepted_field = f | |
| 2088 | ✗ | si%ntau = 1 | |
| 2089 | ✗ | call raw_step(si, f, step_status) | |
| 2090 | ✗ | si%ntau = original_ntau | |
| 2091 | ✗ | if (step_status /= SYMPLECTIC_STEP_OK) then | |
| 2092 | ✗ | si = accepted_integrator | |
| 2093 | ✗ | f = accepted_field | |
| 2094 | ✗ | si%last_step_fraction = real(ktau - 1, dp)/real(original_ntau, dp) | |
| 2095 | ✗ | ierr = step_status | |
| 2096 | 76 | return | |
| 2097 | end if | ||
| 2098 | end do | ||
| 2099 | ✗ | si%last_step_fraction = 1d0 | |
| 2100 | ✗ | return | |
| 2101 | end if | ||
| 2102 | |||
| 2103 |
2/2✓ Branch 0 taken 4402257 times.
✓ Branch 1 taken 4391281 times.
|
8793538 | do ktau = 1, original_ntau |
| 2104 | 4402257 | accepted_integrator = si | |
| 2105 | 4402257 | accepted_field = f | |
| 2106 | 4402257 | si%ntau = 1 | |
| 2107 | 4402257 | call raw_step(si, f, step_status) | |
| 2108 | 4402257 | si%ntau = original_ntau | |
| 2109 |
2/2✓ Branch 0 taken 4402184 times.
✓ Branch 1 taken 73 times.
|
4402257 | if (step_status == SYMPLECTIC_STEP_OK) then |
| 2110 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 4402181 times.
|
4402184 | if (si%z(1) == sympl_rmax) then |
| 2111 | 3 | si%last_step_fraction = real(ktau, dp)/real(original_ntau, dp) | |
| 2112 | 3 | si%last_event_radial_residual = 0d0 | |
| 2113 | 3 | si%last_event_fraction_width = 0d0 | |
| 2114 | 3 | ierr = SYMPLECTIC_STEP_BOUNDARY | |
| 2115 | 3 | return | |
| 2116 | end if | ||
| 2117 | cycle | ||
| 2118 | end if | ||
| 2119 | |||
| 2120 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 9 times.
|
73 | if (step_status == SYMPLECTIC_STEP_OUTSIDE_DOMAIN) then |
| 2121 | 64 | accepted_integrator%ntau = original_ntau | |
| 2122 | call locate_symplectic_boundary(accepted_integrator, accepted_field, & | ||
| 2123 | raw_step, si, f, step_status, event_fraction, radial_residual, & | ||
| 2124 | 64 | fraction_width) | |
| 2125 |
2/2✓ Branch 0 taken 59 times.
✓ Branch 1 taken 5 times.
|
64 | if (step_status == SYMPLECTIC_STEP_BOUNDARY) then |
| 2126 | si%last_step_fraction = & | ||
| 2127 | 59 | (real(ktau - 1, dp) + event_fraction)/real(original_ntau, dp) | |
| 2128 | 59 | si%last_event_radial_residual = radial_residual | |
| 2129 | si%last_event_fraction_width = & | ||
| 2130 | 59 | fraction_width/real(original_ntau, dp) | |
| 2131 | 59 | ierr = step_status | |
| 2132 | 59 | return | |
| 2133 | end if | ||
| 2134 | end if | ||
| 2135 | |||
| 2136 | 14 | si = entry_integrator | |
| 2137 | 14 | f = entry_field | |
| 2138 | 14 | si%last_step_fraction = 0d0 | |
| 2139 | 14 | si%last_event_radial_residual = radial_residual | |
| 2140 | 14 | si%last_event_fraction_width = fraction_width | |
| 2141 | 14 | ierr = step_status | |
| 2142 | 8793538 | return | |
| 2143 | end do | ||
| 2144 | |||
| 2145 | 4391281 | si%last_step_fraction = 1d0 | |
| 2146 | 4391281 | si%ntau = original_ntau | |
| 2147 |
5/6✓ Branch 0 taken 13174071 times.
✓ Branch 1 taken 4391357 times.
✓ Branch 2 taken 26348142 times.
✓ Branch 3 taken 4391357 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4391357 times.
|
43913570 | end subroutine advance_symplectic_with_boundary |
| 2148 | |||
| 2149 | |||
| 2150 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 2151 | ! | ||
| 2152 | 532 | pure subroutine hermite_endpoint(s, nq, c, value, rate, value_end) | |
| 2153 | ! Extrapolate smooth functions of the Gauss collocation polynomial u(tau) | ||
| 2154 | ! to the step end tau = 1 from their values and tau-derivatives at the s | ||
| 2155 | ! stages. By collocation u'(c_k) = f(u(c_k)), so the stage rates follow | ||
| 2156 | ! from the stage field caches, and the 2s Hermite conditions give the | ||
| 2157 | ! endpoint to O(dt**(2s)) per step, the order of the method, with no | ||
| 2158 | ! field evaluation. Columns of value and rate are the nq quantities. | ||
| 2159 | integer, intent(in) :: s, nq | ||
| 2160 | real(dp), intent(in) :: c(s), value(s, nq), rate(s, nq) | ||
| 2161 | real(dp), intent(out) :: value_end(nq) | ||
| 2162 | 1064 | real(dp) :: node(2*s), table(2*s) | |
| 2163 | integer :: i, j, n, q | ||
| 2164 | |||
| 2165 | 532 | n = 2*s | |
| 2166 |
2/2✓ Branch 0 taken 532 times.
✓ Branch 1 taken 1546 times.
|
2078 | do i = 1, s |
| 2167 | 1546 | node(2*i - 1) = c(i) | |
| 2168 | 2078 | node(2*i) = c(i) | |
| 2169 | end do | ||
| 2170 |
2/2✓ Branch 0 taken 3192 times.
✓ Branch 1 taken 532 times.
|
3724 | do q = 1, nq |
| 2171 |
2/2✓ Branch 0 taken 3192 times.
✓ Branch 1 taken 9276 times.
|
12468 | do i = 1, s |
| 2172 | 9276 | table(2*i - 1) = value(i, q) | |
| 2173 | 12468 | table(2*i) = value(i, q) | |
| 2174 | end do | ||
| 2175 | ! Divided differences in place; repeated nodes take the derivative. | ||
| 2176 |
2/2✓ Branch 0 taken 15360 times.
✓ Branch 1 taken 3192 times.
|
18552 | do j = 1, n - 1 |
| 2177 |
2/2✓ Branch 0 taken 53292 times.
✓ Branch 1 taken 15360 times.
|
71844 | do i = n, j + 1, -1 |
| 2178 |
4/4✓ Branch 0 taken 15360 times.
✓ Branch 1 taken 37932 times.
✓ Branch 2 taken 9276 times.
✓ Branch 3 taken 6084 times.
|
68652 | if (j == 1 .and. mod(i, 2) == 0) then |
| 2179 | 9276 | table(i) = rate(i/2, q) | |
| 2180 | else | ||
| 2181 | 44016 | table(i) = (table(i) - table(i - 1))/(node(i) - node(i - j)) | |
| 2182 | end if | ||
| 2183 | end do | ||
| 2184 | end do | ||
| 2185 | 3192 | value_end(q) = table(n) | |
| 2186 |
2/2✓ Branch 0 taken 15360 times.
✓ Branch 1 taken 3192 times.
|
19084 | do i = n - 1, 1, -1 |
| 2187 | 18552 | value_end(q) = table(i) + (1d0 - node(i))*value_end(q) | |
| 2188 | end do | ||
| 2189 | end do | ||
| 2190 | 532 | end subroutine hermite_endpoint | |
| 2191 | |||
| 2192 | 536 | recursive subroutine orbit_timestep_sympl_rk_gauss(si, f, s, ierr) | |
| 2193 | ! | ||
| 2194 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 2195 | type(field_can_t), intent(inout) :: f | ||
| 2196 | |||
| 2197 | integer, intent(out) :: ierr | ||
| 2198 | |||
| 2199 | integer, parameter :: maxit = 32 | ||
| 2200 | |||
| 2201 | integer, intent(in) :: s | ||
| 2202 | 1072 | real(dp), dimension(4*s) :: x, xlast | |
| 2203 | integer :: k, l, ktau, newton_status | ||
| 2204 | |||
| 2205 |
6/6✓ Branch 0 taken 1608 times.
✓ Branch 1 taken 536 times.
✓ Branch 2 taken 3216 times.
✓ Branch 3 taken 536 times.
✓ Branch 4 taken 1556 times.
✓ Branch 5 taken 536 times.
|
6916 | type(field_can_t) :: fs(s) |
| 2206 | type(field_can_t) :: accepted_field | ||
| 2207 | type(symplectic_integrator_t) :: accepted_integrator | ||
| 2208 | 1072 | real(dp) :: a(s,s), b(s), c(s), Hprime(s), dz(4*s) | |
| 2209 | 536 | real(dp) :: stage_value(s, 6), stage_rate(s, 6), endpoint_value(6) | |
| 2210 | real(dp) :: radius_end, r_rate, theta_rate, phi_rate, pphi_rate, canonical_pth | ||
| 2211 | |||
| 2212 |
2/2✓ Branch 0 taken 1556 times.
✓ Branch 1 taken 536 times.
|
2092 | do k = 1,s |
| 2213 | 1556 | fs(k) = f | |
| 2214 | 2092 | fs(k) = f | |
| 2215 | end do | ||
| 2216 | |||
| 2217 | 536 | ierr = 0 | |
| 2218 | 536 | ktau = 0 | |
| 2219 |
2/2✓ Branch 0 taken 536 times.
✓ Branch 1 taken 532 times.
|
1604 | do while(ktau .lt. si%ntau) |
| 2220 | 536 | accepted_integrator = si | |
| 2221 | 536 | accepted_field = f | |
| 2222 | 536 | si%pthold = f%pth | |
| 2223 | |||
| 2224 |
2/2✓ Branch 0 taken 1556 times.
✓ Branch 1 taken 536 times.
|
2092 | do k = 1,s |
| 2225 |
2/2✓ Branch 0 taken 6224 times.
✓ Branch 1 taken 1556 times.
|
8316 | x((4*k-3):(4*k)) = si%z |
| 2226 | end do | ||
| 2227 | |||
| 2228 | call newton_rk_gauss(si, fs, s, x, si%atol, si%rtol, maxit, xlast, & | ||
| 2229 | 536 | newton_status) | |
| 2230 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 532 times.
|
536 | if (newton_status /= SYMPLECTIC_STEP_OK) then |
| 2231 | 4 | si = accepted_integrator | |
| 2232 | 4 | f = accepted_field | |
| 2233 | 4 | ierr = newton_status | |
| 2234 | 4 | return | |
| 2235 | end if | ||
| 2236 | !optionally try fixed point iterations, doesn't work yet | ||
| 2237 | !call fixpoint_rk_gauss(si, fs, s, x, si%atol, si%rtol, maxit, xlast) | ||
| 2238 | |||
| 2239 |
4/6✓ Branch 0 taken 1546 times.
✓ Branch 1 taken 532 times.
✓ Branch 2 taken 1546 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 532 times.
|
2078 | if (any(x(1:4*s:4) > sympl_rmax)) then |
| 2240 | ✗ | si = accepted_integrator | |
| 2241 | ✗ | f = accepted_field | |
| 2242 | ✗ | ierr = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 2243 | ✗ | return | |
| 2244 | end if | ||
| 2245 | |||
| 2246 | 532 | call coeff_rk_gauss(s, a, b, c) ! TODO: move this to preprocessing | |
| 2247 | |||
| 2248 | if (extrap_field) then | ||
| 2249 |
2/2✓ Branch 1 taken 1546 times.
✓ Branch 2 taken 532 times.
|
2078 | do k = 1, s |
| 2250 | 1546 | dz(1) = x(4*k-3)-xlast(4*k-3) | |
| 2251 | 1546 | dz(2) = x(4*k-2)-xlast(4*k-2) | |
| 2252 | 1546 | dz(3) = x(4*k-1)-xlast(4*k-1) | |
| 2253 | 1546 | dz(4) = x(4*k)-xlast(4*k) | |
| 2254 | |||
| 2255 | fs(k)%pth = fs(k)%pth + fs(k)%dpth(1)*dz(1) & ! d/dr | ||
| 2256 | + fs(k)%dpth(2)*dz(2) & ! d/dth | ||
| 2257 | + fs(k)%dpth(3)*dz(3) & ! d/dph | ||
| 2258 | 1546 | + fs(k)%dpth(4)*dz(4) ! d/dpph | |
| 2259 | |||
| 2260 | fs(k)%dpth(1) = fs(k)%dpth(1) + fs(k)%d2pth(1)*dz(1) & | ||
| 2261 | + fs(k)%d2pth(2)*dz(2) & | ||
| 2262 | + fs(k)%d2pth(3)*dz(3) & | ||
| 2263 | 1546 | + fs(k)%d2pth(7)*dz(4) | |
| 2264 | |||
| 2265 | fs(k)%dpth(2) = fs(k)%dpth(2) + fs(k)%d2pth(2)*dz(1)& | ||
| 2266 | + fs(k)%d2pth(4)*dz(2) & | ||
| 2267 | + fs(k)%d2pth(5)*dz(3) & | ||
| 2268 | 1546 | + fs(k)%d2pth(8)*dz(4) | |
| 2269 | |||
| 2270 | fs(k)%dpth(3) = fs(k)%dpth(3) + fs(k)%d2pth(3)*dz(1) & | ||
| 2271 | + fs(k)%d2pth(5)*dz(2) & | ||
| 2272 | + fs(k)%d2pth(6)*dz(3) & | ||
| 2273 | 1546 | + fs(k)%d2pth(9)*dz(4) | |
| 2274 | |||
| 2275 | fs(k)%dH(1) = fs(k)%dH(1) + fs(k)%d2H(1)*dz(1) & | ||
| 2276 | + fs(k)%d2H(2)*dz(2) & | ||
| 2277 | + fs(k)%d2H(3)*dz(3) & | ||
| 2278 | 1546 | + fs(k)%d2H(7)*dz(4) | |
| 2279 | |||
| 2280 | fs(k)%dH(2) = fs(k)%dH(2) + fs(k)%d2H(2)*dz(1) & | ||
| 2281 | + fs(k)%d2H(4)*dz(2) & | ||
| 2282 | + fs(k)%d2H(5)*dz(3) & | ||
| 2283 | 1546 | + fs(k)%d2H(8)*dz(4) | |
| 2284 | |||
| 2285 | fs(k)%dH(3) = fs(k)%dH(3) + fs(k)%d2H(3)*dz(1) & | ||
| 2286 | + fs(k)%d2H(5)*dz(2) & | ||
| 2287 | + fs(k)%d2H(6)*dz(3) & | ||
| 2288 | 1546 | + fs(k)%d2H(9)*dz(4) | |
| 2289 | |||
| 2290 | fs(k)%vpar = fs(k)%vpar + fs(k)%dvpar(1)*dz(1) & ! d/dr | ||
| 2291 | + fs(k)%dvpar(2)*dz(2) & | ||
| 2292 | + fs(k)%dvpar(3)*dz(3) & | ||
| 2293 | 1546 | + fs(k)%dvpar(4)*dz(4) | |
| 2294 | |||
| 2295 | fs(k)%hth = fs(k)%hth + fs(k)%dhth(1)*dz(1) & ! d/dr | ||
| 2296 | + fs(k)%dhth(2)*dz(2) & | ||
| 2297 | 1546 | + fs(k)%dhth(3)*dz(3) | |
| 2298 | |||
| 2299 | fs(k)%hph = fs(k)%hph + fs(k)%dhph(1)*dz(1) & ! d/dr | ||
| 2300 | + fs(k)%dhph(2)*dz(2) & | ||
| 2301 | 1546 | + fs(k)%dhph(3)*dz(3) | |
| 2302 | |||
| 2303 | ! Not used by the update, but extrapolated to the endpoint cache. | ||
| 2304 | fs(k)%Ath = fs(k)%Ath + fs(k)%dAth(1)*dz(1) + fs(k)%dAth(2)*dz(2) & | ||
| 2305 | 1546 | + fs(k)%dAth(3)*dz(3) | |
| 2306 | fs(k)%Aph = fs(k)%Aph + fs(k)%dAph(1)*dz(1) + fs(k)%dAph(2)*dz(2) & | ||
| 2307 | 1546 | + fs(k)%dAph(3)*dz(3) | |
| 2308 | fs(k)%Bmod = fs(k)%Bmod + fs(k)%dBmod(1)*dz(1) & | ||
| 2309 | 2078 | + fs(k)%dBmod(2)*dz(2) + fs(k)%dBmod(3)*dz(3) | |
| 2310 | end do | ||
| 2311 | else | ||
| 2312 | do k = 1, s | ||
| 2313 | call eval_field(fs(k), x(4*k-3), x(4*k-2), x(4*k-1), 0) | ||
| 2314 | call get_derivatives(fs(k), x(4*k)) | ||
| 2315 | end do | ||
| 2316 | end if | ||
| 2317 | |||
| 2318 |
2/2✓ Branch 0 taken 1546 times.
✓ Branch 1 taken 532 times.
|
2078 | do k = 1, s |
| 2319 | 2078 | Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) | |
| 2320 | end do | ||
| 2321 | |||
| 2322 | 532 | f = fs(s) | |
| 2323 | 532 | f%pth = si%pthold | |
| 2324 | 532 | si%z(1) = x(4*s-3) | |
| 2325 | |||
| 2326 |
2/2✓ Branch 0 taken 1546 times.
✓ Branch 1 taken 532 times.
|
2078 | do l = 1, s |
| 2327 | 1546 | f%pth = f%pth - si%dt*b(l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) | |
| 2328 | 1546 | si%z(2) = si%z(2) + si%dt*b(l)*Hprime(l) | |
| 2329 | 1546 | si%z(3) = si%z(3)+si%dt*b(l)*(fs(l)%vpar-Hprime(l)*fs(l)%hth)/fs(l)%hph | |
| 2330 | 2078 | si%z(4) = si%z(4) - si%dt*b(l)*(fs(l)%dH(3) - Hprime(l)*fs(l)%dpth(3)) | |
| 2331 | end do | ||
| 2332 | |||
| 2333 | ! The update advances (theta, phi, p_theta, p_phi); the radius is | ||
| 2334 | ! defined only implicitly by p_theta. Report the endpoint radius and | ||
| 2335 | ! field cache by Hermite extrapolation of the stages. The error does not | ||
| 2336 | ! accumulate: each step solves its stage radii from the constraint | ||
| 2337 | ! again, and the endpoint only seeds the next Newton guess. | ||
| 2338 | ! Stage values and rates in tau = t/dt of r and of the cached field | ||
| 2339 | ! values callers read after the step (vpar and H follow from p_phi). | ||
| 2340 | ! dr/dt follows from d/dt p_theta(r, theta, phi, p_phi) = dp_theta/dt | ||
| 2341 | ! along the stage; the theta terms cancel against p_theta's own rate. | ||
| 2342 |
2/2✓ Branch 0 taken 1546 times.
✓ Branch 1 taken 532 times.
|
2078 | do k = 1, s |
| 2343 | 1546 | theta_rate = si%dt*Hprime(k) | |
| 2344 | 1546 | phi_rate = si%dt*(fs(k)%vpar - Hprime(k)*fs(k)%hth)/fs(k)%hph | |
| 2345 | 1546 | pphi_rate = -si%dt*(fs(k)%dH(3) - Hprime(k)*fs(k)%dpth(3)) | |
| 2346 | r_rate = (-si%dt*fs(k)%dH(2) - fs(k)%dpth(3)*phi_rate & | ||
| 2347 | 1546 | - fs(k)%hth/fs(k)%hph*pphi_rate)/fs(k)%dpth(1) | |
| 2348 | 1546 | stage_value(k, 1) = x(4*k - 3) | |
| 2349 | 1546 | stage_value(k, 2) = fs(k)%Ath | |
| 2350 | 1546 | stage_value(k, 3) = fs(k)%Aph | |
| 2351 | 1546 | stage_value(k, 4) = fs(k)%hth | |
| 2352 | 1546 | stage_value(k, 5) = fs(k)%hph | |
| 2353 | 1546 | stage_value(k, 6) = fs(k)%Bmod | |
| 2354 | 1546 | stage_rate(k, 1) = r_rate | |
| 2355 | 1546 | stage_rate(k, 2) = chain_rate(fs(k)%dAth) | |
| 2356 | 1546 | stage_rate(k, 3) = chain_rate(fs(k)%dAph) | |
| 2357 | 1546 | stage_rate(k, 4) = chain_rate(fs(k)%dhth) | |
| 2358 | 1546 | stage_rate(k, 5) = chain_rate(fs(k)%dhph) | |
| 2359 | 532 | stage_rate(k, 6) = chain_rate(fs(k)%dBmod) | |
| 2360 | end do | ||
| 2361 | 532 | call hermite_endpoint(s, 6, c, stage_value, stage_rate, endpoint_value) | |
| 2362 | 532 | radius_end = endpoint_value(1) | |
| 2363 |
3/4✓ Branch 0 taken 532 times.
✓ Branch 1 taken 532 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 532 times.
|
1064 | if (.not. finite_newton_iterate([radius_end])) then |
| 2364 | ✗ | si = accepted_integrator | |
| 2365 | ✗ | f = accepted_field | |
| 2366 | ✗ | ierr = SYMPLECTIC_STEP_MAXITER | |
| 2367 | ✗ | return | |
| 2368 | end if | ||
| 2369 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 532 times.
|
532 | if (radius_end > sympl_rmax) then |
| 2370 | ! A crossing inside the last stage interval; the boundary event solve | ||
| 2371 | ! locates it on the step fraction. | ||
| 2372 | ✗ | si = accepted_integrator | |
| 2373 | ✗ | f = accepted_field | |
| 2374 | ✗ | ierr = SYMPLECTIC_STEP_OUTSIDE_DOMAIN | |
| 2375 | ✗ | return | |
| 2376 | end if | ||
| 2377 | ! The polar chart cannot represent a negative radius. Near the axis keep | ||
| 2378 | ! the last stage radius and cache, as before; the next step solves its | ||
| 2379 | ! stages anew. | ||
| 2380 |
1/2✓ Branch 0 taken 532 times.
✗ Branch 1 not taken.
|
532 | if (radius_end >= 0d0) then |
| 2381 | 532 | si%z(1) = radius_end | |
| 2382 | 532 | canonical_pth = f%pth | |
| 2383 | 532 | f%Ath = endpoint_value(2) | |
| 2384 | 532 | f%Aph = endpoint_value(3) | |
| 2385 | 532 | f%hth = endpoint_value(4) | |
| 2386 | 532 | f%hph = endpoint_value(5) | |
| 2387 | 532 | f%Bmod = endpoint_value(6) | |
| 2388 | 532 | call get_val(f, si%z(4)) | |
| 2389 | 532 | f%pth = canonical_pth | |
| 2390 | end if | ||
| 2391 | |||
| 2392 | 532 | ktau = ktau+1 | |
| 2393 | enddo | ||
| 2394 | |||
| 2395 | contains | ||
| 2396 | |||
| 2397 | 1546 | pure real(dp) function chain_rate(gradient) | |
| 2398 | ! d/dtau of a field value along the stage from its (r, theta, phi) gradient. | ||
| 2399 | real(dp), intent(in) :: gradient(3) | ||
| 2400 | |||
| 2401 | chain_rate = gradient(1)*r_rate + gradient(2)*theta_rate & | ||
| 2402 | 1546 | + gradient(3)*phi_rate | |
| 2403 | end function chain_rate | ||
| 2404 | |||
| 2405 | end subroutine orbit_timestep_sympl_rk_gauss | ||
| 2406 | |||
| 2407 | |||
| 2408 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 2409 | ! | ||
| 2410 | 212 | recursive subroutine orbit_timestep_sympl_rk_lobatto(si, f, s, ierr) | |
| 2411 | ! | ||
| 2412 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 2413 | type(field_can_t), intent(inout) :: f | ||
| 2414 | integer, intent(out) :: ierr | ||
| 2415 | |||
| 2416 | integer, parameter :: maxit = 32 | ||
| 2417 | |||
| 2418 | integer, intent(in) :: s | ||
| 2419 | 424 | real(dp), dimension(4*s-2) :: x, xlast | |
| 2420 | integer :: ktau, k, l, newton_status | ||
| 2421 | |||
| 2422 |
6/6✓ Branch 0 taken 636 times.
✓ Branch 1 taken 212 times.
✓ Branch 2 taken 1272 times.
✓ Branch 3 taken 212 times.
✓ Branch 4 taken 636 times.
✓ Branch 5 taken 212 times.
|
2756 | type(field_can_t) :: fs(s) |
| 2423 | type(field_can_t) :: accepted_field | ||
| 2424 | type(symplectic_integrator_t) :: accepted_integrator | ||
| 2425 | 163 | real(dp) :: a(s,s), ahat(s,s), b(s), c(s), Hprime(s) | |
| 2426 | |||
| 2427 |
2/2✓ Branch 0 taken 636 times.
✓ Branch 1 taken 212 times.
|
848 | do k = 1,s |
| 2428 | 848 | fs(k) = f | |
| 2429 | end do | ||
| 2430 | |||
| 2431 | 212 | ierr = 0 | |
| 2432 | 212 | ktau = 0 | |
| 2433 |
2/2✓ Branch 0 taken 212 times.
✓ Branch 1 taken 163 times.
|
375 | do while(ktau .lt. si%ntau) |
| 2434 | 212 | accepted_integrator = si | |
| 2435 | 212 | accepted_field = f | |
| 2436 | 212 | si%pthold = f%pth | |
| 2437 | |||
| 2438 | 212 | x(1) = si%z(1) | |
| 2439 | 212 | x(2) = si%z(4) | |
| 2440 |
2/2✓ Branch 0 taken 424 times.
✓ Branch 1 taken 212 times.
|
636 | do k = 2,s |
| 2441 |
2/2✓ Branch 0 taken 1696 times.
✓ Branch 1 taken 424 times.
|
2332 | x((4*k-3-2):(4*k-2)) = si%z |
| 2442 | end do | ||
| 2443 | |||
| 2444 | call newton_rk_lobatto(si, fs, s, x, si%atol, si%rtol, maxit, xlast, & | ||
| 2445 | 212 | newton_status) | |
| 2446 |
2/2✓ Branch 0 taken 49 times.
✓ Branch 1 taken 163 times.
|
212 | if (newton_status /= SYMPLECTIC_STEP_OK) then |
| 2447 | 49 | si = accepted_integrator | |
| 2448 | 49 | f = accepted_field | |
| 2449 | 49 | ierr = newton_status | |
| 2450 | return | ||
| 2451 | end if | ||
| 2452 | |||
| 2453 | 163 | call coeff_rk_lobatto(s, a, ahat, b, c) | |
| 2454 | |||
| 2455 | 163 | call eval_field(fs(1), x(1), si%z(2), si%z(3), 0) | |
| 2456 | 163 | call get_derivatives(fs(1), x(2)) | |
| 2457 | 163 | Hprime(1) = fs(1)%dH(1)/fs(1)%dpth(1) | |
| 2458 | |||
| 2459 |
2/2✓ Branch 0 taken 326 times.
✓ Branch 1 taken 163 times.
|
489 | do k = 2, s |
| 2460 | 326 | call eval_field(fs(k), x(4*k-3-2), x(4*k-2-2), x(4*k-1-2), 0) | |
| 2461 | 326 | call get_derivatives(fs(k), x(4*k-2)) | |
| 2462 | 489 | Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) | |
| 2463 | end do | ||
| 2464 | |||
| 2465 | 163 | f = fs(s) | |
| 2466 | 163 | f%pth = si%pthold | |
| 2467 | 163 | si%z(1) = x(4*s-3-2) | |
| 2468 | |||
| 2469 |
2/2✓ Branch 0 taken 489 times.
✓ Branch 1 taken 163 times.
|
652 | do l = 1, s |
| 2470 | 489 | f%pth = f%pth - si%dt*b(l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) | |
| 2471 | 489 | si%z(2) = si%z(2) + si%dt*b(l)*Hprime(l) | |
| 2472 | 489 | si%z(3) = si%z(3) + si%dt*b(l)*(fs(l)%vpar-Hprime(l)*fs(l)%hth)/fs(l)%hph | |
| 2473 | 652 | si%z(4) = si%z(4) - si%dt*b(l)*(fs(l)%dH(3) - Hprime(l)*fs(l)%dpth(3)) | |
| 2474 | end do | ||
| 2475 | |||
| 2476 | 163 | ktau = ktau+1 | |
| 2477 | enddo | ||
| 2478 | |||
| 2479 | 212 | end subroutine orbit_timestep_sympl_rk_lobatto | |
| 2480 | |||
| 2481 | |||
| 2482 | ✗ | recursive subroutine debug_root(si, f, x0) | |
| 2483 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 2484 | type(field_can_t), pointer :: f | ||
| 2485 | |||
| 2486 | real(dp) :: x0(2), x(2) | ||
| 2487 | integer :: k, l, iflag | ||
| 2488 | integer, parameter :: n = 100 | ||
| 2489 | real(dp), parameter :: eps = 1d-15 | ||
| 2490 | |||
| 2491 | real(dp) :: fvec(2) | ||
| 2492 | |||
| 2493 | ✗ | do k = -n,n | |
| 2494 | ✗ | do l = -n,n | |
| 2495 | ✗ | x = x0 + l*eps/n*(/x0(1),0d0/) + k*eps/n*(/0d0,x0(2)/) | |
| 2496 | ✗ | call f_sympl_euler1(si, f, 2, x, fvec, iflag) | |
| 2497 | ✗ | write(5001,*) x(1), x(2), fvec | |
| 2498 | end do | ||
| 2499 | end do | ||
| 2500 | |||
| 2501 | ✗ | end subroutine debug_root | |
| 2502 | |||
| 2503 | end module orbit_symplectic | ||
| 2504 |