| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module spectre_sympl_orbit | ||
| 2 | !> Per-volume guiding-centre stepping with full-orbit interface excursions. | ||
| 3 | !> Interior sheets are crossed only by the full orbit; no fixed-mu map. | ||
| 4 | |||
| 5 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 6 | use finite_check, only: is_finite | ||
| 7 | !$ use omp_lib, only: omp_get_max_threads, omp_get_thread_num | ||
| 8 | use util, only: twopi, sqrt2 | ||
| 9 | use parmot_mod, only: ro0 | ||
| 10 | use field_can_mod, only: field_can_t, eval_field => evaluate, integ_to_ref, & | ||
| 11 | ref_to_integ | ||
| 12 | use field_can_spectre, only: set_spectre_volume_lock | ||
| 13 | use orbit_symplectic_base, only: symplectic_integrator_t, & | ||
| 14 | symplectic_newton_warning_mode, symplectic_spectre_newton_warning_factor, & | ||
| 15 | SYMPLECTIC_STEP_OUTSIDE_DOMAIN | ||
| 16 | use orbit_symplectic, only: orbit_timestep_sympl, orbit_sympl_init | ||
| 17 | use spectre_fo_hybrid, only: spectre_fo_state_t, spectre_fo_enter, & | ||
| 18 | spectre_fo_advance_until_exit, spectre_fo_radial_larmor, SPECTRE_FO_OK, & | ||
| 19 | SPECTRE_FO_LOSS | ||
| 20 | use params, only: spectre_fo_entry_larmor | ||
| 21 | use diag_counters, only: count_event, EVT_WARNING_STEP_SKIP | ||
| 22 | |||
| 23 | implicit none | ||
| 24 | private | ||
| 25 | |||
| 26 | public :: sympl_spectre_state_t, sympl_spectre_reset, recanon_pphi, & | ||
| 27 | orbit_microstep_sympl_spectre, sympl_landing_stats_reset, & | ||
| 28 | sympl_landing_stats, sympl_landing_eval_stats | ||
| 29 | public :: sympl_fo_stats, sympl_bulk_recovery_stats | ||
| 30 | public :: SYMPL_SPECTRE_OK, SYMPL_SPECTRE_LOSS, SYMPL_SPECTRE_STOP | ||
| 31 | |||
| 32 | integer, parameter :: SYMPL_SPECTRE_OK = 0 | ||
| 33 | integer, parameter :: SYMPL_SPECTRE_LOSS = 88 | ||
| 34 | integer, parameter :: SYMPL_SPECTRE_STOP = 89 | ||
| 35 | integer, parameter :: STOP_SHEET = 1, STOP_TELEPORT = 2, STOP_STEP = 3 | ||
| 36 | !> Landing target and acceptance for |rho_g - k| at the interface. The | ||
| 37 | !> h*-solve aims at landing_tol; the per-volume splines are clamped at the | ||
| 38 | !> volume edge, which kinks the landed rho_g just outside the interface, so a | ||
| 39 | !> best-effort landing below landing_accept is still committed. | ||
| 40 | real(dp), parameter :: landing_tol = 1.0d-10 | ||
| 41 | real(dp), parameter :: landing_accept = 1.0d-8 | ||
| 42 | integer, parameter :: max_locate = 64 | ||
| 43 | |||
| 44 | !> Near-tangent crossings can put the landing root where the implicit solve | ||
| 45 | !> of a long substep is multi-branched (Newton basin boundaries), making | ||
| 46 | !> g(h) discontinuous. The substep is then halved and re-attempted: committed | ||
| 47 | !> half-steps are fixed-length symplectic maps, and the landing solve re-runs | ||
| 48 | !> in a smaller-h regime where the step map is single-valued. | ||
| 49 | integer, parameter :: max_iters = 200 | ||
| 50 | real(dp), parameter :: h_min_frac = 1.0d-8 | ||
| 51 | real(dp), parameter :: budget_eps = 1.0d-12 | ||
| 52 | !> Loose state-sanity bounds for a committed substep. The Hamiltonian is not | ||
| 53 | !> required to stay constant from step to step: symplectic discretizations | ||
| 54 | !> normally oscillate around a nearby modified Hamiltonian. The shell bound | ||
| 55 | !> only rejects a state that has left the tested bounded-error band around | ||
| 56 | !> the fixed physical speed pabs. GC/FO conversion preserves the momentum shell. | ||
| 57 | real(dp), parameter :: max_step_dr = 0.5d0 | ||
| 58 | real(dp), parameter :: max_relative_shell_defect = 5.0d-4 | ||
| 59 | |||
| 60 | !> Home-side offset for the per-volume angle-transform evaluation at a landed | ||
| 61 | !> point: the volume dispatch keys on int(rho_g), so exactly at rho_g = k it | ||
| 62 | !> would pick the outer volume, whose zeta gauge differs by O(1). | ||
| 63 | real(dp), parameter :: transform_bias = 1.0d-8 | ||
| 64 | |||
| 65 | type :: sympl_spectre_state_t | ||
| 66 | real(dp) :: home_lo = 0.0_dp | ||
| 67 | real(dp) :: home_hi = 1.0_dp | ||
| 68 | integer :: mvol = 1 | ||
| 69 | integer :: mode = 0 | ||
| 70 | integer :: level = 3 | ||
| 71 | real(dp) :: dt_std = 0.0_dp | ||
| 72 | !> Accepted FO clock, in the same units as dt_std; study diagnostic. | ||
| 73 | real(dp) :: fo_elapsed = 0.0_dp | ||
| 74 | !> An exhausted warning-mode recovery holds one unresolved interval. | ||
| 75 | !> The next microstep retries guiding-centre advancement from the last | ||
| 76 | !> valid state. A later identical failure consumes one unresolved | ||
| 77 | !> interval directly instead of repeating the full recovery cascade; | ||
| 78 | !> default warning mode never turns that retry into a marker stop. | ||
| 79 | logical :: warning_hold_latched = .false. | ||
| 80 | integer :: warning_hold_count = 0 | ||
| 81 | integer :: warning_hold_reason = STOP_STEP | ||
| 82 | integer :: warning_hold_iface = 0 | ||
| 83 | integer :: warning_hold_direction = 0 | ||
| 84 | real(dp) :: warning_hold_z(4) = 0.0_dp | ||
| 85 | type(spectre_fo_state_t) :: fo | ||
| 86 | end type sympl_spectre_state_t | ||
| 87 | |||
| 88 | !> Hot-path statistics are per-thread. Each OpenMP worker writes only its | ||
| 89 | !> own slot; totals are reduced after the parallel trace has joined. | ||
| 90 | integer, allocatable :: n_landings(:), n_landing_evals(:), n_stops(:) | ||
| 91 | integer, allocatable :: n_fo_entries(:), n_fo_exits(:), n_fo_losses(:) | ||
| 92 | integer, allocatable :: n_fo_failures(:), n_fo_failure_status(:, :) | ||
| 93 | !> Bulk GC failures handed to FO recovery: iteration cap, persistent | ||
| 94 | !> teleport, persistent Newton failure at the minimum substep. | ||
| 95 | integer, parameter :: BULK_ITERS = 1, BULK_TELEPORT = 2, BULK_NEWTON = 3 | ||
| 96 | integer, allocatable :: n_bulk_recover(:, :) | ||
| 97 | real(dp), allocatable :: max_landing_resid(:) | ||
| 98 | |||
| 99 | contains | ||
| 100 | |||
| 101 |
10/10✓ Branch 0 taken 321 times.
✓ Branch 1 taken 1284 times.
✓ Branch 2 taken 321 times.
✓ Branch 3 taken 963 times.
✓ Branch 4 taken 321 times.
✓ Branch 5 taken 963 times.
✓ Branch 6 taken 321 times.
✓ Branch 7 taken 963 times.
✓ Branch 8 taken 1605 times.
✓ Branch 9 taken 321 times.
|
7383 | subroutine sympl_spectre_reset(state, si, mvol, mode, level) |
| 102 | type(sympl_spectre_state_t), intent(out) :: state | ||
| 103 | type(symplectic_integrator_t), intent(in) :: si | ||
| 104 | integer, intent(in) :: mvol, mode, level | ||
| 105 | |||
| 106 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 321 times.
|
321 | if (level /= 3) error stop 'crossing_level must be 3; GC maps were removed' |
| 107 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 321 times.
|
321 | if (mode <= 0) then |
| 108 | error stop 'spectre_sympl_orbit: crossing pipeline requires a '// & | ||
| 109 | ✗ | 'symplectic scheme (integmode > 0)' | |
| 110 | end if | ||
| 111 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 321 times.
|
321 | if (si%ntau /= 1) then |
| 112 | ! Re-canonicalization restarts the scheme from a single point, so no | ||
| 113 | ! multistep history may carry over; all schemes here are one-step and | ||
| 114 | ! must be driven one step per call. | ||
| 115 | ✗ | error stop 'spectre_sympl_orbit: requires one-step driving (ntau = 1)' | |
| 116 | end if | ||
| 117 | |||
| 118 | 321 | state%mvol = mvol | |
| 119 | 321 | state%mode = mode | |
| 120 | 321 | state%level = level | |
| 121 | 321 | state%dt_std = si%dt | |
| 122 | 321 | call set_home(state, si%z(1)) | |
| 123 | 321 | end subroutine sympl_spectre_reset | |
| 124 | |||
| 125 | 321 | subroutine set_home(state, rho) | |
| 126 | !> Fix the home volume [home_lo, home_hi] bracketing rho and pin the | ||
| 127 | !> canonical field evaluation to it, so Newton iterates probing past an | ||
| 128 | !> interior interface see the home volume's clamped splines instead of | ||
| 129 | !> the neighbour's discontinuous field. | ||
| 130 | type(sympl_spectre_state_t), intent(inout) :: state | ||
| 131 | real(dp), intent(in) :: rho | ||
| 132 | |||
| 133 | real(dp) :: lo | ||
| 134 | |||
| 135 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 321 times.
|
321 | lo = real(floor(rho), dp) |
| 136 | 321 | lo = max(0.0_dp, min(lo, real(state%mvol - 1, dp))) | |
| 137 | 321 | state%home_lo = lo | |
| 138 | 321 | state%home_hi = lo + 1.0_dp | |
| 139 | 321 | call set_spectre_volume_lock(nint(state%home_hi)) | |
| 140 | 321 | end subroutine set_home | |
| 141 | |||
| 142 | ✗ | subroutine set_home_volume(state, lvol) | |
| 143 | type(sympl_spectre_state_t), intent(inout) :: state | ||
| 144 | integer, intent(in) :: lvol | ||
| 145 | |||
| 146 | ✗ | state%home_lo = real(lvol - 1, dp) | |
| 147 | ✗ | state%home_hi = real(lvol, dp) | |
| 148 | ✗ | call set_spectre_volume_lock(lvol) | |
| 149 | end subroutine set_home_volume | ||
| 150 | |||
| 151 | 821862 | subroutine orbit_microstep_sympl_spectre(state, si, f, ierr, t_frac) | |
| 152 | !> Advance one microstep of length state%dt_std, resolving interface | ||
| 153 | !> events by exact-landing substeps. On SYMPL_SPECTRE_LOSS or | ||
| 154 | !> SYMPL_SPECTRE_STOP the orbit terminates at fraction t_frac of the | ||
| 155 | !> microstep, with (si, f) holding the terminal (landed) state. | ||
| 156 | type(sympl_spectre_state_t), intent(inout) :: state | ||
| 157 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 158 | type(field_can_t), intent(inout) :: f | ||
| 159 | integer, intent(out) :: ierr | ||
| 160 | real(dp), intent(out) :: t_frac | ||
| 161 | |||
| 162 | type(symplectic_integrator_t) :: si0 | ||
| 163 | type(field_can_t) :: f0 | ||
| 164 | real(dp) :: budget, used, h_try, h_land, resid, y_iface(5), y_out(5) | ||
| 165 | real(dp) :: exit_rho | ||
| 166 | integer :: iters, iface, direction, ierr_step, ierr_fo | ||
| 167 | integer :: fo_owner | ||
| 168 | logical :: boundary, solved, fo_exited, resume_gc | ||
| 169 | |||
| 170 | 821862 | ierr = SYMPL_SPECTRE_OK | |
| 171 | 821862 | t_frac = 1.0_dp | |
| 172 | 821862 | si%warning_factor = symplectic_spectre_newton_warning_factor | |
| 173 | 821862 | budget = state%dt_std | |
| 174 | 821862 | used = 0.0_dp | |
| 175 | 821862 | h_try = budget | |
| 176 | |||
| 177 | 821918 | if (state%fo%active) then | |
| 178 | call continue_fo(state, si, f, budget, used, y_out, fo_owner, & | ||
| 179 | 56 | fo_exited, ierr_fo) | |
| 180 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
|
56 | if (ierr_fo /= SPECTRE_FO_OK) then |
| 181 | ✗ | if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & | |
| 182 | budget, used, ierr, t_frac)) then | ||
| 183 | ✗ | si%dt = state%dt_std | |
| 184 | 56 | return | |
| 185 | end if | ||
| 186 | ✗ | call finish_fo_error(ierr_fo, used, state%dt_std, ierr, t_frac) | |
| 187 | ✗ | si%dt = state%dt_std | |
| 188 | ✗ | return | |
| 189 | end if | ||
| 190 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
56 | if (.not. fo_exited .or. budget <= budget_eps*state%dt_std) then |
| 191 | 56 | si%dt = state%dt_std | |
| 192 | 56 | return | |
| 193 | end if | ||
| 194 | ✗ | h_try = budget | |
| 195 | end if | ||
| 196 | |||
| 197 | iters = 0 | ||
| 198 | 1643580 | do | |
| 199 | 821806 | iters = iters + 1 | |
| 200 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 821806 times.
|
821806 | if (iters > max_iters) then |
| 201 | ✗ | call count_bulk_recovery(BULK_ITERS) | |
| 202 | call recover_local_error(state, si, f, STOP_STEP, budget, used, & | ||
| 203 | ✗ | ierr, t_frac, resume_gc) | |
| 204 | ✗ | if (resume_gc) then | |
| 205 | ✗ | iters = 0 | |
| 206 | ✗ | h_try = budget | |
| 207 | ✗ | cycle | |
| 208 | end if | ||
| 209 | exit | ||
| 210 | end if | ||
| 211 |
1/2✓ Branch 0 taken 821806 times.
✗ Branch 1 not taken.
|
821806 | if (state%level == 3) then |
| 212 | 821806 | call interface_proximity(state, si, f, y_iface, iface, exit_rho) | |
| 213 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 821790 times.
|
821806 | if (iface > 0) then |
| 214 | call start_fo(state, si, f, y_iface, nint(state%home_hi), & | ||
| 215 | iface, budget, used, y_out, fo_owner, fo_exited, ierr_fo, & | ||
| 216 | 16 | exit_rho) | |
| 217 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
16 | if (ierr_fo /= SPECTRE_FO_OK) then |
| 218 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & |
| 219 | budget, used, ierr, t_frac)) exit | ||
| 220 | call finish_fo_error(ierr_fo, used, state%dt_std, ierr, & | ||
| 221 | 821814 | t_frac) | |
| 222 | exit | ||
| 223 | end if | ||
| 224 |
1/4✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
8 | if (.not. fo_exited .or. & |
| 225 | budget <= budget_eps*state%dt_std) exit | ||
| 226 | ✗ | h_try = budget | |
| 227 | ✗ | cycle | |
| 228 | end if | ||
| 229 | end if | ||
| 230 | 821790 | si0 = si | |
| 231 | 821790 | f0 = f | |
| 232 | 821790 | si%dt = h_try | |
| 233 | 821790 | call orbit_timestep_sympl(si, f, ierr_step) | |
| 234 | |||
| 235 |
1/2✓ Branch 0 taken 821790 times.
✗ Branch 1 not taken.
|
821790 | if (ierr_step == 0) then |
| 236 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 821790 times.
|
821790 | if (step_teleported(si0, si, f)) then |
| 237 | ✗ | si = si0 | |
| 238 | ✗ | f = f0 | |
| 239 | ✗ | if (same_warning_hold(state, si, STOP_TELEPORT, 0, 0)) then | |
| 240 | ✗ | if (recover_warning_skip(state, si, budget, used, ierr, & | |
| 241 | t_frac, STOP_TELEPORT, 0, 0)) exit | ||
| 242 | end if | ||
| 243 | ✗ | h_try = 0.5_dp*h_try | |
| 244 | ✗ | if (h_try >= h_min_frac*state%dt_std) cycle | |
| 245 | ✗ | call count_bulk_recovery(BULK_TELEPORT) | |
| 246 | call recover_local_error(state, si, f, STOP_TELEPORT, budget, & | ||
| 247 | ✗ | used, ierr, t_frac, resume_gc) | |
| 248 | ✗ | if (resume_gc) then | |
| 249 | ✗ | h_try = budget | |
| 250 | ✗ | cycle | |
| 251 | end if | ||
| 252 | exit | ||
| 253 | end if | ||
| 254 | 821790 | state%warning_hold_latched = .false. | |
| 255 | end if | ||
| 256 | |||
| 257 | call classify_step(state, si0%z(1), si, ierr_step, boundary, iface, & | ||
| 258 | 821790 | direction) | |
| 259 | |||
| 260 |
2/2✓ Branch 0 taken 821654 times.
✓ Branch 1 taken 136 times.
|
821790 | if (.not. boundary) then |
| 261 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 821654 times.
|
821654 | if (ierr_step /= 0) then |
| 262 | ! Newton divergence away from the outer boundary: retry with | ||
| 263 | ! a halved substep, like the teleport and stalled-landing | ||
| 264 | ! branches. | ||
| 265 | ✗ | si = si0 | |
| 266 | ✗ | f = f0 | |
| 267 | ✗ | if (same_warning_hold(state, si, STOP_STEP, 0, 0)) then | |
| 268 | ✗ | if (recover_warning_skip(state, si, budget, used, ierr, & | |
| 269 | t_frac, STOP_STEP, 0, 0)) exit | ||
| 270 | end if | ||
| 271 | ✗ | h_try = 0.5_dp*h_try | |
| 272 | ✗ | if (h_try >= h_min_frac*state%dt_std) cycle | |
| 273 | ✗ | call count_bulk_recovery(BULK_NEWTON) | |
| 274 | call recover_local_error(state, si, f, STOP_STEP, budget, & | ||
| 275 | ✗ | used, ierr, t_frac, resume_gc) | |
| 276 | ✗ | if (resume_gc) then | |
| 277 | ✗ | h_try = budget | |
| 278 | ✗ | cycle | |
| 279 | end if | ||
| 280 | exit | ||
| 281 | end if | ||
| 282 | 821654 | used = used + h_try | |
| 283 | 821654 | budget = state%dt_std - used | |
| 284 |
1/2✓ Branch 0 taken 821654 times.
✗ Branch 1 not taken.
|
821654 | if (budget <= budget_eps*state%dt_std) exit |
| 285 | ! Grow the working substep gently after halvings: resetting to | ||
| 286 | ! the full remaining budget would re-trigger the same failed | ||
| 287 | ! solve indefinitely. | ||
| 288 | ✗ | h_try = min(2.0_dp*h_try, budget) | |
| 289 | ✗ | cycle | |
| 290 | end if | ||
| 291 | |||
| 292 | call locate_landing(si, f, si0, f0, h_try, real(iface, dp), & | ||
| 293 | 136 | direction, ierr_step == 0, h_land, resid, solved) | |
| 294 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 136 times.
|
136 | if (.not. solved) then |
| 295 | ! The landing root sits where the implicit solve of this substep | ||
| 296 | ! length is multi-branched: halve and re-attempt from the same | ||
| 297 | ! pre-step state. | ||
| 298 | ✗ | si = si0 | |
| 299 | ✗ | f = f0 | |
| 300 | ✗ | h_try = 0.5_dp*h_try | |
| 301 | ✗ | if (h_try >= h_min_frac*state%dt_std) cycle | |
| 302 | ✗ | call integrator_state(si, f, y_iface) | |
| 303 | call start_fo(state, si, f, y_iface, nint(state%home_hi), & | ||
| 304 | ✗ | iface, budget, used, y_out, fo_owner, fo_exited, ierr_fo) | |
| 305 | ✗ | if (ierr_fo /= SPECTRE_FO_OK) then | |
| 306 | ✗ | if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & | |
| 307 | budget, used, ierr, t_frac)) exit | ||
| 308 | 821806 | call finish_fo_error(ierr_fo, used, state%dt_std, ierr, t_frac) | |
| 309 | exit | ||
| 310 | end if | ||
| 311 | ✗ | if (.not. fo_exited .or. & | |
| 312 | budget <= budget_eps*state%dt_std) exit | ||
| 313 | ✗ | h_try = budget | |
| 314 | ✗ | cycle | |
| 315 | end if | ||
| 316 | 136 | call landed_state(si, f, real(iface, dp), direction, y_iface) | |
| 317 | 136 | used = used + h_land | |
| 318 | 136 | call update_landing_stats(resid) | |
| 319 |
1/2✓ Branch 0 taken 136 times.
✗ Branch 1 not taken.
|
136 | if (state%mvol == 1) then |
| 320 | ! The separate single-volume GC model retains its GC edge. | ||
| 321 | 136 | ierr = SYMPL_SPECTRE_LOSS | |
| 322 | 136 | t_frac = used/state%dt_std | |
| 323 | 136 | exit | |
| 324 | end if | ||
| 325 | ! A multi-volume hybrid resolves both sheets and its absorbing | ||
| 326 | ! particle boundary in FO if a GC step missed the proximity band. | ||
| 327 | ✗ | budget = state%dt_std - used | |
| 328 | ✗ | call fo_exit_distance(state, y_iface, exit_rho) | |
| 329 | call start_fo(state, si, f, y_iface, nint(state%home_hi), iface, & | ||
| 330 | ✗ | budget, used, y_out, fo_owner, fo_exited, ierr_fo, exit_rho) | |
| 331 | ✗ | if (ierr_fo /= SPECTRE_FO_OK) then | |
| 332 | ✗ | if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & | |
| 333 | budget, used, ierr, t_frac)) exit | ||
| 334 | 821806 | call finish_fo_error(ierr_fo, used, state%dt_std, ierr, t_frac) | |
| 335 | exit | ||
| 336 | end if | ||
| 337 | ✗ | if (.not. fo_exited .or. budget <= budget_eps*state%dt_std) exit | |
| 338 | ✗ | h_try = budget | |
| 339 | 821662 | cycle | |
| 340 | |||
| 341 | end do | ||
| 342 | |||
| 343 | 821806 | si%dt = state%dt_std | |
| 344 |
6/6✓ Branch 0 taken 2465586 times.
✓ Branch 1 taken 821862 times.
✓ Branch 2 taken 4931172 times.
✓ Branch 3 taken 821862 times.
✓ Branch 4 taken 56 times.
✓ Branch 5 taken 821806 times.
|
8218620 | end subroutine orbit_microstep_sympl_spectre |
| 345 | |||
| 346 | 821790 | subroutine classify_step(state, rho_start, si, ierr_step, boundary, iface, & | |
| 347 | direction) | ||
| 348 | type(sympl_spectre_state_t), intent(in) :: state | ||
| 349 | real(dp), intent(in) :: rho_start | ||
| 350 | type(symplectic_integrator_t), intent(in) :: si | ||
| 351 | integer, intent(in) :: ierr_step | ||
| 352 | logical, intent(out) :: boundary | ||
| 353 | integer, intent(out) :: iface, direction | ||
| 354 | |||
| 355 | 821790 | boundary = .false. | |
| 356 | 821790 | iface = 0 | |
| 357 | 821790 | direction = 0 | |
| 358 | |||
| 359 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 821790 times.
|
821790 | if (ierr_step /= 0) then |
| 360 | ! A trial exceeding the Newton margin (Mvol + 1) from the outer | ||
| 361 | ! volume still brackets its physical edge at Mvol, with si%z | ||
| 362 | ! holding the pre-step state. Successful trials are checked below. | ||
| 363 | ✗ | if (ierr_step == SYMPLECTIC_STEP_OUTSIDE_DOMAIN .and. & | |
| 364 | nint(state%home_hi) == state%mvol) then | ||
| 365 | ✗ | boundary = .true. | |
| 366 | ✗ | iface = state%mvol | |
| 367 | ✗ | direction = 1 | |
| 368 | end if | ||
| 369 | ✗ | return | |
| 370 | end if | ||
| 371 | |||
| 372 |
3/6✓ Branch 0 taken 821654 times.
✓ Branch 1 taken 136 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 821654 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
821790 | if (si%z(1) > state%home_hi .or. & |
| 373 | (si%z(1) == state%home_hi .and. rho_start < state%home_hi)) then | ||
| 374 | 136 | boundary = .true. | |
| 375 | 136 | iface = nint(state%home_hi) | |
| 376 | 136 | direction = 1 | |
| 377 |
2/6✓ Branch 0 taken 821654 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 821654 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
821654 | else if (si%z(1) < state%home_lo .or. & |
| 378 | (si%z(1) == state%home_lo .and. rho_start > state%home_lo)) then | ||
| 379 | ! home_lo = 0 is the coordinate axis, not an interface: the chart | ||
| 380 | ! flip (r, theta) -> (-r, theta + pi) inside the steppers handles | ||
| 381 | ! axis flybys (#370), so only real interfaces raise events here. | ||
| 382 | ✗ | if (nint(state%home_lo) >= 1) then | |
| 383 | ✗ | boundary = .true. | |
| 384 | ✗ | iface = nint(state%home_lo) | |
| 385 | ✗ | direction = -1 | |
| 386 | end if | ||
| 387 | end if | ||
| 388 | end subroutine classify_step | ||
| 389 | |||
| 390 | 136 | subroutine locate_landing(si, f, si0, f0, h_full, rho_k, direction, & | |
| 391 | full_committed, h_land, resid, solved) | ||
| 392 | !> Solve g(h) = dir*(z1(1; h) - rho_k) = 0 for the substep length by | ||
| 393 | !> Illinois false position; every evaluation is one full implicit step of | ||
| 394 | !> length h from the same pre-step state (si0, f0). Spline evaluations | ||
| 395 | !> past the volume edge are clamped (field_can_spectre), which kinks g | ||
| 396 | !> just outside the interface; the shrinking bracket walks through the | ||
| 397 | !> kink, and trials the sympl_rmax guard refuses to commit (outermost | ||
| 398 | !> interface only) tighten the upper bound by bisection. On return | ||
| 399 | !> (si, f) hold the best landed step of length h_land. | ||
| 400 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 401 | type(symplectic_integrator_t), intent(in) :: si0 | ||
| 402 | type(field_can_t), intent(inout) :: f | ||
| 403 | type(field_can_t), intent(in) :: f0 | ||
| 404 | real(dp), intent(in) :: h_full, rho_k | ||
| 405 | integer, intent(in) :: direction | ||
| 406 | logical, intent(in) :: full_committed | ||
| 407 | real(dp), intent(out) :: h_land, resid | ||
| 408 | logical, intent(out) :: solved | ||
| 409 | |||
| 410 | type(symplectic_integrator_t) :: si_best | ||
| 411 | type(field_can_t) :: f_best | ||
| 412 | real(dp) :: dir, t_lo, t_hi, f_lo, f_hi, t_mid, g, best | ||
| 413 | integer :: it, ierr_step, last_side, evals | ||
| 414 | logical :: hi_valid | ||
| 415 | |||
| 416 | 136 | dir = real(direction, dp) | |
| 417 | 136 | t_lo = 0.0_dp | |
| 418 | 136 | f_lo = dir*(si0%z(1) - rho_k) | |
| 419 | 136 | t_hi = h_full | |
| 420 | 136 | f_hi = 0.0_dp | |
| 421 | 136 | hi_valid = .false. | |
| 422 | 136 | last_side = 0 | |
| 423 | 136 | evals = 0 | |
| 424 | |||
| 425 | 136 | best = huge(1.0_dp) | |
| 426 | 136 | h_land = h_full | |
| 427 | 136 | if (full_committed) then | |
| 428 | 136 | f_hi = dir*(si%z(1) - rho_k) | |
| 429 | 136 | hi_valid = .true. | |
| 430 | 136 | best = abs(f_hi) | |
| 431 | 136 | si_best = si | |
| 432 | 136 | f_best = f | |
| 433 | end if | ||
| 434 | |||
| 435 |
1/2✓ Branch 0 taken 660 times.
✗ Branch 1 not taken.
|
660 | do it = 1, max_locate |
| 436 |
2/2✓ Branch 0 taken 136 times.
✓ Branch 1 taken 524 times.
|
660 | if (best < landing_tol) exit |
| 437 |
1/2✓ Branch 0 taken 524 times.
✗ Branch 1 not taken.
|
524 | if (hi_valid) then |
| 438 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 524 times.
|
524 | if (abs(f_hi - f_lo) <= tiny(1.0_dp)) exit |
| 439 | 524 | t_mid = t_hi - f_hi*(t_hi - t_lo)/(f_hi - f_lo) | |
| 440 |
2/4✓ Branch 0 taken 524 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 524 times.
|
524 | if (t_mid <= t_lo .or. t_mid >= t_hi) then |
| 441 | ✗ | t_mid = 0.5_dp*(t_lo + t_hi) | |
| 442 | end if | ||
| 443 | else | ||
| 444 | ✗ | t_mid = 0.5_dp*(t_lo + t_hi) | |
| 445 | end if | ||
| 446 | |||
| 447 | 524 | si = si0 | |
| 448 | 524 | f = f0 | |
| 449 | 524 | si%dt = t_mid | |
| 450 | 524 | call orbit_timestep_sympl(si, f, ierr_step) | |
| 451 | 524 | evals = evals + 1 | |
| 452 | 524 | g = dir*(si%z(1) - rho_k) | |
| 453 |
3/6✓ Branch 0 taken 524 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 524 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 524 times.
✗ Branch 5 not taken.
|
524 | if (ierr_step /= 0 .or. g /= g .or. & |
| 454 | step_teleported(si0, si, f)) then | ||
| 455 | t_hi = t_mid | ||
| 456 | hi_valid = .false. | ||
| 457 | cycle | ||
| 458 | end if | ||
| 459 | |||
| 460 |
1/2✓ Branch 0 taken 524 times.
✗ Branch 1 not taken.
|
524 | if (abs(g) < best) then |
| 461 | 524 | best = abs(g) | |
| 462 | 524 | h_land = t_mid | |
| 463 | 524 | si_best = si | |
| 464 | 524 | f_best = f | |
| 465 | end if | ||
| 466 | |||
| 467 |
2/2✓ Branch 1 taken 324 times.
✓ Branch 2 taken 200 times.
|
1184 | if (g < 0.0_dp) then |
| 468 | 324 | t_lo = t_mid | |
| 469 | 324 | f_lo = g | |
| 470 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 214 times.
|
324 | if (last_side == 1) then |
| 471 |
1/2✓ Branch 0 taken 110 times.
✗ Branch 1 not taken.
|
110 | if (hi_valid) f_hi = 0.5_dp*f_hi |
| 472 | end if | ||
| 473 | last_side = 1 | ||
| 474 | else | ||
| 475 | 200 | t_hi = t_mid | |
| 476 | 200 | f_hi = g | |
| 477 | 200 | hi_valid = .true. | |
| 478 |
2/2✓ Branch 0 taken 52 times.
✓ Branch 1 taken 148 times.
|
200 | if (last_side == -1) f_lo = 0.5_dp*f_lo |
| 479 | last_side = -1 | ||
| 480 | end if | ||
| 481 | end do | ||
| 482 | |||
| 483 | 136 | solved = best < landing_accept | |
| 484 | 136 | resid = best | |
| 485 |
1/2✓ Branch 0 taken 136 times.
✗ Branch 1 not taken.
|
136 | if (solved) then |
| 486 | 136 | si = si_best | |
| 487 | 136 | f = f_best | |
| 488 | end if | ||
| 489 |
1/2✓ Branch 0 taken 136 times.
✗ Branch 1 not taken.
|
136 | if (allocated(n_landing_evals)) then |
| 490 | n_landing_evals(stats_thread()) = & | ||
| 491 | 408 | n_landing_evals(stats_thread()) + evals | |
| 492 | end if | ||
| 493 |
5/6✓ Branch 0 taken 408 times.
✓ Branch 1 taken 136 times.
✓ Branch 2 taken 816 times.
✓ Branch 3 taken 136 times.
✓ Branch 4 taken 136 times.
✗ Branch 5 not taken.
|
1360 | end subroutine locate_landing |
| 494 | |||
| 495 | 136 | subroutine landed_state(si, f, rho_k, direction, y) | |
| 496 | !> Physical state y = (rho_g, theta, zeta, p, lambda) in SPECTRE reference | ||
| 497 | !> coordinates at the landed integrator state. The zeta transform is read | ||
| 498 | !> a hair inside the home volume (transform_bias) and unwrapped so the | ||
| 499 | !> angle winding survives the modulo in integ_to_ref. | ||
| 500 | type(symplectic_integrator_t), intent(in) :: si | ||
| 501 | type(field_can_t), intent(in) :: f | ||
| 502 | real(dp), intent(in) :: rho_k | ||
| 503 | integer, intent(in) :: direction | ||
| 504 | real(dp), intent(out) :: y(5) | ||
| 505 | |||
| 506 | real(dp) :: xinteg(3), xref(3), p | ||
| 507 | |||
| 508 | 136 | xinteg(1) = rho_k - real(direction, dp)*transform_bias | |
| 509 |
2/2✓ Branch 0 taken 272 times.
✓ Branch 1 taken 136 times.
|
408 | xinteg(2:3) = si%z(2:3) |
| 510 | 136 | call integ_to_ref(xinteg, xref) | |
| 511 | |||
| 512 | 136 | y(1) = si%z(1) | |
| 513 | 136 | y(2) = si%z(2) | |
| 514 | 136 | y(3) = unwrap_near(xref(3), si%z(3)) | |
| 515 | 136 | p = si%pabs | |
| 516 | 136 | y(4) = p | |
| 517 | 136 | y(5) = shell_pitch(f, p) | |
| 518 | 136 | end subroutine landed_state | |
| 519 | |||
| 520 | 821806 | subroutine integrator_state(si, f, y) | |
| 521 | type(symplectic_integrator_t), intent(in) :: si | ||
| 522 | type(field_can_t), intent(in) :: f | ||
| 523 | real(dp), intent(out) :: y(5) | ||
| 524 | |||
| 525 | real(dp) :: xref(3), p | ||
| 526 | |||
| 527 | 821806 | call integ_to_ref(si%z(1:3), xref) | |
| 528 | 821806 | y(1) = si%z(1) | |
| 529 | 821806 | y(2) = si%z(2) | |
| 530 | 821806 | y(3) = unwrap_near(xref(3), si%z(3)) | |
| 531 | 821806 | p = si%pabs | |
| 532 | 821806 | y(4) = p | |
| 533 | 821806 | y(5) = shell_pitch(f, p) | |
| 534 | 821806 | end subroutine integrator_state | |
| 535 | |||
| 536 | 821942 | pure function shell_pitch(f, p) result(lambda) | |
| 537 | type(field_can_t), intent(in) :: f | ||
| 538 | real(dp), intent(in) :: p | ||
| 539 | real(dp) :: lambda, lambda2 | ||
| 540 | |||
| 541 | 821942 | lambda2 = max(1.0_dp - f%mu*f%Bmod/p**2, 0.0_dp) | |
| 542 | 821942 | lambda = sign(sqrt(lambda2), f%vpar) | |
| 543 | 821942 | end function shell_pitch | |
| 544 | |||
| 545 | 821806 | subroutine interface_proximity(state, si, f, y, iface, exit_rho) | |
| 546 | !> Nearest interface or absorbing outer edge if the GC is closer to | ||
| 547 | !> it than spectre_fo_entry_larmor radial gyroradii, else iface = 0. y | ||
| 548 | !> is the guiding-centre state and exit_rho the matching exit distance. | ||
| 549 | type(sympl_spectre_state_t), intent(in) :: state | ||
| 550 | type(symplectic_integrator_t), intent(in) :: si | ||
| 551 | type(field_can_t), intent(in) :: f | ||
| 552 | real(dp), intent(out) :: y(5), exit_rho | ||
| 553 | integer, intent(out) :: iface | ||
| 554 | |||
| 555 | real(dp) :: rho_r | ||
| 556 | integer :: k, status | ||
| 557 | |||
| 558 | 821806 | iface = 0 | |
| 559 | 821806 | exit_rho = 0.0_dp | |
| 560 | 821806 | call integrator_state(si, f, y) | |
| 561 |
2/2✓ Branch 0 taken 10016 times.
✓ Branch 1 taken 811790 times.
|
831806 | if (state%mvol == 1) return |
| 562 | 10016 | k = nint(y(1)) | |
| 563 |
2/4✓ Branch 0 taken 10016 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10016 times.
✗ Branch 3 not taken.
|
10016 | if (k < 1 .or. k > state%mvol) return |
| 564 | call spectre_fo_radial_larmor(y, nint(state%home_hi), ro0/sqrt2, rho_r, & | ||
| 565 | 10016 | status) | |
| 566 |
1/2✓ Branch 0 taken 10016 times.
✗ Branch 1 not taken.
|
10016 | if (status /= SPECTRE_FO_OK) return |
| 567 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 10000 times.
|
10016 | if (abs(y(1) - real(k, dp)) >= spectre_fo_entry_larmor*rho_r) return |
| 568 | 16 | iface = k | |
| 569 | 16 | exit_rho = (spectre_fo_entry_larmor + 1.0_dp)*rho_r | |
| 570 | end subroutine interface_proximity | ||
| 571 | |||
| 572 | ✗ | subroutine fo_exit_distance(state, y, exit_rho) | |
| 573 | !> Exit distance (rho_g) for a full orbit started on an interface. | ||
| 574 | type(sympl_spectre_state_t), intent(in) :: state | ||
| 575 | real(dp), intent(in) :: y(5) | ||
| 576 | real(dp), intent(out) :: exit_rho | ||
| 577 | |||
| 578 | real(dp) :: rho_r | ||
| 579 | integer :: status | ||
| 580 | |||
| 581 | call spectre_fo_radial_larmor(y, nint(state%home_hi), ro0/sqrt2, rho_r, & | ||
| 582 | ✗ | status) | |
| 583 | ✗ | exit_rho = 0.0_dp | |
| 584 | ✗ | if (status == SPECTRE_FO_OK) & | |
| 585 | ✗ | exit_rho = (spectre_fo_entry_larmor + 1.0_dp)*rho_r | |
| 586 | ✗ | end subroutine fo_exit_distance | |
| 587 | |||
| 588 | 16 | subroutine start_fo(state, si, f, y, owner, iface, budget, used, y_out, & | |
| 589 | owner_out, exited, ierr, exit_rho) | ||
| 590 | type(sympl_spectre_state_t), intent(inout) :: state | ||
| 591 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 592 | type(field_can_t), intent(inout) :: f | ||
| 593 | real(dp), intent(in) :: y(5) | ||
| 594 | integer, intent(in) :: owner, iface | ||
| 595 | real(dp), intent(inout) :: budget, used | ||
| 596 | real(dp), intent(out) :: y_out(5) | ||
| 597 | integer, intent(out) :: owner_out, ierr | ||
| 598 | logical, intent(out) :: exited | ||
| 599 | !> Leave the full orbit once the guiding centre is this far (rho_g) | ||
| 600 | !> from the interface; default: the gyrophase-based legacy rule. | ||
| 601 | real(dp), intent(in), optional :: exit_rho | ||
| 602 | |||
| 603 | 16 | y_out = y | |
| 604 | 16 | owner_out = owner | |
| 605 | 16 | exited = .false. | |
| 606 | 16 | call spectre_fo_enter(state%fo, y, owner, iface, ro0/sqrt2, ierr) | |
| 607 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 9 times.
|
16 | if (ierr == SPECTRE_FO_LOSS) then |
| 608 | ! A GC initially inside the wall can map to a particle outside it. | ||
| 609 | ! This is an immediate physical loss, with no accepted FO time. | ||
| 610 | 7 | call count_fo_entry | |
| 611 |
1/2✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
|
7 | if (allocated(n_fo_losses)) n_fo_losses(stats_thread()) = & |
| 612 | 21 | n_fo_losses(stats_thread()) + 1 | |
| 613 | 7 | return | |
| 614 | end if | ||
| 615 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
|
9 | if (ierr /= SPECTRE_FO_OK) then |
| 616 | ✗ | call count_fo_failure(ierr) | |
| 617 | ✗ | return | |
| 618 | end if | ||
| 619 |
1/2✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
|
9 | if (present(exit_rho)) state%fo%exit_rho = exit_rho |
| 620 | 9 | call count_fo_entry | |
| 621 | call continue_fo(state, si, f, budget, used, y_out, owner_out, & | ||
| 622 | 9 | exited, ierr) | |
| 623 | end subroutine start_fo | ||
| 624 | |||
| 625 | 65 | subroutine continue_fo(state, si, f, budget, used, y_out, owner_out, & | |
| 626 | exited, ierr) | ||
| 627 | type(sympl_spectre_state_t), intent(inout) :: state | ||
| 628 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 629 | type(field_can_t), intent(inout) :: f | ||
| 630 | real(dp), intent(inout) :: budget, used | ||
| 631 | real(dp), intent(out) :: y_out(5) | ||
| 632 | integer, intent(out) :: owner_out, ierr | ||
| 633 | logical, intent(out) :: exited | ||
| 634 | |||
| 635 | real(dp) :: dt_used | ||
| 636 | |||
| 637 | call spectre_fo_advance_until_exit(state%fo, budget, ro0/sqrt2, & | ||
| 638 | 65 | dt_used, y_out, owner_out, exited, ierr) | |
| 639 | 65 | state%fo_elapsed = state%fo_elapsed + dt_used | |
| 640 | 65 | used = used + dt_used | |
| 641 | 65 | budget = state%dt_std - used | |
| 642 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 64 times.
|
65 | if (ierr /= SPECTRE_FO_OK) then |
| 643 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (ierr == SPECTRE_FO_LOSS) then |
| 644 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (allocated(n_fo_losses)) n_fo_losses(stats_thread()) = & |
| 645 | 3 | n_fo_losses(stats_thread()) + 1 | |
| 646 | else | ||
| 647 | ✗ | call count_fo_failure(ierr) | |
| 648 | end if | ||
| 649 | 65 | return | |
| 650 | end if | ||
| 651 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 64 times.
|
64 | if (.not. exited) return |
| 652 | ✗ | call count_fo_exit | |
| 653 | ✗ | call set_home_volume(state, owner_out) | |
| 654 | ✗ | call recanonicalize(state, si, f, y_out) | |
| 655 | end subroutine continue_fo | ||
| 656 | |||
| 657 | ✗ | subroutine recanonicalize(state, si, f, y) | |
| 658 | !> Re-initialize the integrator from the physical state y in the explicitly | ||
| 659 | !> owned volume, exactly like a fresh orbit start | ||
| 660 | !> (init_sympl): per-volume gauges differ, so canonical momenta are | ||
| 661 | !> rebuilt from the target volume's own Ath, Aph, hth, hph. | ||
| 662 | type(sympl_spectre_state_t), intent(in) :: state | ||
| 663 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 664 | type(field_can_t), intent(inout) :: f | ||
| 665 | real(dp), intent(in) :: y(5) | ||
| 666 | |||
| 667 | real(dp) :: xref(3), xinteg(3), z(4) | ||
| 668 | |||
| 669 | ✗ | xref = y(1:3) | |
| 670 | ✗ | call ref_to_integ(xref, xinteg) | |
| 671 | ✗ | xinteg(2) = y(2) | |
| 672 | ✗ | xinteg(3) = unwrap_near(xinteg(3), y(3)) | |
| 673 | |||
| 674 | ✗ | call eval_field(f, xinteg(1), xinteg(2), xinteg(3), 0) | |
| 675 | ✗ | si%pabs = y(4) | |
| 676 | ✗ | f%mu = y(4)**2*(1.0_dp - y(5)**2)/f%Bmod | |
| 677 | ✗ | f%ro0 = ro0/sqrt2 | |
| 678 | ✗ | f%vpar = y(4)*y(5)*sqrt2 | |
| 679 | ✗ | z(1:3) = xinteg | |
| 680 | ✗ | z(4) = recanon_pphi(f, y(4), y(5)) | |
| 681 | ✗ | call orbit_sympl_init(si, f, z, state%dt_std, 1, si%rtol, state%mode) | |
| 682 | ✗ | si%warning_factor = symplectic_spectre_newton_warning_factor | |
| 683 | ✗ | end subroutine recanonicalize | |
| 684 | |||
| 685 | ✗ | pure function recanon_pphi(f, p, lambda) result(pphi) | |
| 686 | !> Canonical toroidal momentum of the physical state (p, lambda) in the | ||
| 687 | !> gauge of the volume f is evaluated in. Correctness of re-canonicalization | ||
| 688 | !> is that f carries the TARGET volume's Aph/hph: decoding pphi via get_val | ||
| 689 | !> in the same volume returns vpar = p*lambda*sqrt2. A neighbour volume's | ||
| 690 | !> gauge shifts the decoded vpar by [[Aph]]/(f%ro0*f%hph). | ||
| 691 | type(field_can_t), intent(in) :: f | ||
| 692 | real(dp), intent(in) :: p, lambda | ||
| 693 | real(dp) :: pphi | ||
| 694 | |||
| 695 | ✗ | pphi = p*lambda*sqrt2*f%hph + f%Aph/f%ro0 | |
| 696 | ✗ | end function recanon_pphi | |
| 697 | |||
| 698 | 822314 | pure function step_teleported(si0, si, f) result(bad) | |
| 699 | type(symplectic_integrator_t), intent(in) :: si0, si | ||
| 700 | type(field_can_t), intent(in) :: f | ||
| 701 | logical :: bad | ||
| 702 | |||
| 703 | bad = .not. all(is_finite(si%z)) .or. & | ||
| 704 | .not. is_finite(si%pabs) .or. & | ||
| 705 | .not. is_finite(f%Bmod) .or. & | ||
| 706 | .not. is_finite(f%vpar) .or. & | ||
| 707 | .not. is_finite(f%mu) .or. & | ||
| 708 |
9/16✓ Branch 0 taken 3289256 times.
✓ Branch 1 taken 822314 times.
✓ Branch 2 taken 3289256 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 822314 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 822314 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 822314 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 822314 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 822314 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 822314 times.
✗ Branch 15 not taken.
|
4111570 | .not. is_finite(f%H) |
| 709 | if (bad) return | ||
| 710 | 822314 | bad = abs(si%z(1) - si0%z(1)) > max_step_dr | |
| 711 |
1/2✓ Branch 0 taken 822314 times.
✗ Branch 1 not taken.
|
822314 | if (bad) return |
| 712 | bad = abs(f%H - si%pabs**2) > & | ||
| 713 | 822314 | max_relative_shell_defect*max(si%pabs**2, tiny(1.0_dp)) | |
| 714 |
1/2✓ Branch 0 taken 822314 times.
✗ Branch 1 not taken.
|
822314 | if (bad) return |
| 715 | ! The steppers commit a negative-radius Newton root as an axis chart | ||
| 716 | ! flip (r, theta) -> (|r|, theta + pi). Away from the axis a negative | ||
| 717 | ! root is a Newton divergence, not an axis crossing: a committed theta | ||
| 718 | ! jump of order pi at large radius is a teleport. | ||
| 719 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 822314 times.
|
822314 | if (abs(si%z(2) - si0%z(2)) > 1.0_dp) then |
| 720 | ✗ | bad = si%z(1) > 0.5_dp | |
| 721 | if (bad) return | ||
| 722 | end if | ||
| 723 | end function step_teleported | ||
| 724 | |||
| 725 | ✗ | subroutine recover_local_error(state, si, f, reason, budget, used, & | |
| 726 | ierr, t_frac, resume_gc) | ||
| 727 | !> A bulk or axis-local solver failure has no interface metadata. Do not | ||
| 728 | !> invent an interior interface from nint(rho): interface FO recovery is | ||
| 729 | !> valid only after an actual landed interface event. | ||
| 730 | type(sympl_spectre_state_t), intent(inout) :: state | ||
| 731 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 732 | type(field_can_t), intent(inout) :: f | ||
| 733 | integer, intent(in) :: reason | ||
| 734 | real(dp), intent(inout) :: budget, used | ||
| 735 | integer, intent(out) :: ierr | ||
| 736 | real(dp), intent(out) :: t_frac | ||
| 737 | logical, intent(out) :: resume_gc | ||
| 738 | |||
| 739 | real(dp) :: y(5), y_out(5) | ||
| 740 | integer :: ierr_fo, owner, owner_out | ||
| 741 | logical :: fo_exited | ||
| 742 | |||
| 743 | ✗ | resume_gc = .false. | |
| 744 | ✗ | call integrator_state(si, f, y) | |
| 745 | ✗ | owner = nint(state%home_hi) | |
| 746 | call start_fo(state, si, f, y, owner, 0, budget, used, y_out, & | ||
| 747 | ✗ | owner_out, fo_exited, ierr_fo) | |
| 748 | ✗ | if (ierr_fo == SPECTRE_FO_OK) then | |
| 749 | ✗ | resume_gc = fo_exited .and. budget > budget_eps*state%dt_std | |
| 750 | ✗ | return | |
| 751 | end if | ||
| 752 | ✗ | if (ierr_fo == SPECTRE_FO_LOSS) then | |
| 753 | ✗ | call finish_fo_error(ierr_fo, used, state%dt_std, ierr, t_frac) | |
| 754 | ✗ | return | |
| 755 | end if | ||
| 756 | ! Bulk recovery has no interface or sheet. Preserve the last valid FO | ||
| 757 | ! state if one exists, but retain the original local failure reason and | ||
| 758 | ! the local stop statistics. | ||
| 759 | ✗ | if (state%fo%has_y) then | |
| 760 | ✗ | call set_home_volume(state, owner_out) | |
| 761 | ✗ | call recanonicalize(state, si, f, state%fo%last_y) | |
| 762 | end if | ||
| 763 | |||
| 764 | ✗ | if (recover_warning_skip(state, si, budget, used, ierr, t_frac, reason, & | |
| 765 | 0, 0)) return | ||
| 766 | |||
| 767 | ✗ | call record_stop | |
| 768 | ✗ | ierr = SYMPL_SPECTRE_STOP | |
| 769 | ✗ | t_frac = used/state%dt_std | |
| 770 | end subroutine recover_local_error | ||
| 771 | |||
| 772 | 821942 | pure function unwrap_near(a, near) result(b) | |
| 773 | real(dp), intent(in) :: a, near | ||
| 774 | real(dp) :: b | ||
| 775 | |||
| 776 | ✗ | b = a + twopi*nint((near - a)/twopi) | |
| 777 | end function unwrap_near | ||
| 778 | |||
| 779 | 8 | logical function recover_fo_failure(state, si, f, fo_status, owner, budget, & | |
| 780 | used, ierr, t_frac) | ||
| 781 | !> Full orbit is the final active recovery. If it also fails numerically, | ||
| 782 | !> return to its last valid guiding-centre state and use the warning-mode | ||
| 783 | !> hold fallback. A physical full-orbit loss is never intercepted. | ||
| 784 | type(sympl_spectre_state_t), intent(inout) :: state | ||
| 785 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 786 | type(field_can_t), intent(inout) :: f | ||
| 787 | integer, intent(in) :: fo_status, owner | ||
| 788 | real(dp), intent(inout) :: budget, used | ||
| 789 | integer, intent(inout) :: ierr | ||
| 790 | real(dp), intent(inout) :: t_frac | ||
| 791 | |||
| 792 | 8 | recover_fo_failure = .false. | |
| 793 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | if (fo_status == SPECTRE_FO_LOSS) return |
| 794 | ✗ | if (.not. symplectic_newton_warning_mode) return | |
| 795 | ✗ | if (state%fo%has_y) then | |
| 796 | ✗ | call set_home_volume(state, owner) | |
| 797 | ✗ | call recanonicalize(state, si, f, state%fo%last_y) | |
| 798 | end if | ||
| 799 | recover_fo_failure = recover_warning_skip(state, si, budget, used, & | ||
| 800 | ierr, t_frac, STOP_SHEET, state%fo%iface, & | ||
| 801 | ✗ | merge(1, -1, owner == state%fo%iface)) | |
| 802 | ✗ | end function recover_fo_failure | |
| 803 | |||
| 804 | ✗ | logical function recover_warning_skip(state, si, budget, used, ierr, & | |
| 805 | t_frac, reason, iface, direction) | ||
| 806 | !> Last-resort default: after all active recovery has failed, retain the | ||
| 807 | !> last valid state for this unresolved interval and resume the marker on | ||
| 808 | !> the next microstep. Strict mode returns .false. to expose the stop. | ||
| 809 | type(sympl_spectre_state_t), intent(inout) :: state | ||
| 810 | type(symplectic_integrator_t), intent(in) :: si | ||
| 811 | real(dp), intent(inout) :: budget, used | ||
| 812 | integer, intent(inout) :: ierr | ||
| 813 | real(dp), intent(inout) :: t_frac | ||
| 814 | integer, intent(in) :: reason, iface, direction | ||
| 815 | |||
| 816 | ✗ | recover_warning_skip = symplectic_newton_warning_mode | |
| 817 | ✗ | if (.not. recover_warning_skip) return | |
| 818 | |||
| 819 | ✗ | used = used + budget | |
| 820 | ✗ | budget = 0.0_dp | |
| 821 | ✗ | ierr = SYMPL_SPECTRE_OK | |
| 822 | ✗ | t_frac = 1.0_dp | |
| 823 | ✗ | call count_event(EVT_WARNING_STEP_SKIP) | |
| 824 | |||
| 825 | ✗ | state%warning_hold_latched = .true. | |
| 826 | ✗ | state%warning_hold_count = state%warning_hold_count + 1 | |
| 827 | ✗ | state%warning_hold_reason = reason | |
| 828 | ✗ | state%warning_hold_iface = iface | |
| 829 | ✗ | state%warning_hold_direction = direction | |
| 830 | ✗ | state%warning_hold_z = si%z | |
| 831 | |||
| 832 | ✗ | state%fo%active = .false. | |
| 833 | ✗ | end function recover_warning_skip | |
| 834 | |||
| 835 | ✗ | pure logical function same_warning_hold(state, si, reason, iface, direction) | |
| 836 | type(sympl_spectre_state_t), intent(in) :: state | ||
| 837 | type(symplectic_integrator_t), intent(in) :: si | ||
| 838 | integer, intent(in) :: reason, iface, direction | ||
| 839 | |||
| 840 | same_warning_hold = state%warning_hold_latched .and. & | ||
| 841 | reason == state%warning_hold_reason .and. & | ||
| 842 | iface == state%warning_hold_iface .and. & | ||
| 843 | direction == state%warning_hold_direction .and. & | ||
| 844 | maxval(abs(si%z - state%warning_hold_z)) <= & | ||
| 845 | 100.0_dp*epsilon(1.0_dp)*max(1.0_dp, & | ||
| 846 | ✗ | maxval(abs(state%warning_hold_z))) | |
| 847 | ✗ | end function same_warning_hold | |
| 848 | |||
| 849 | ✗ | subroutine record_stop | |
| 850 | ✗ | if (allocated(n_stops)) then | |
| 851 | ✗ | n_stops(stats_thread()) = n_stops(stats_thread()) + 1 | |
| 852 | end if | ||
| 853 | ✗ | end subroutine record_stop | |
| 854 | |||
| 855 | 8 | subroutine finish_fo_error(fo_status, used, dt, ierr, t_frac) | |
| 856 | integer, intent(in) :: fo_status | ||
| 857 | real(dp), intent(in) :: used, dt | ||
| 858 | integer, intent(out) :: ierr | ||
| 859 | real(dp), intent(out) :: t_frac | ||
| 860 | |||
| 861 |
1/8✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
8 | if (fo_status == SPECTRE_FO_LOSS) then |
| 862 | 8 | ierr = SYMPL_SPECTRE_LOSS | |
| 863 | else | ||
| 864 | ✗ | call record_stop | |
| 865 | ✗ | ierr = SYMPL_SPECTRE_STOP | |
| 866 | end if | ||
| 867 | ✗ | t_frac = used/dt | |
| 868 | 8 | end subroutine finish_fo_error | |
| 869 | |||
| 870 | 136 | subroutine update_landing_stats(resid) | |
| 871 | real(dp), intent(in) :: resid | ||
| 872 | |||
| 873 |
1/2✓ Branch 0 taken 136 times.
✗ Branch 1 not taken.
|
136 | if (allocated(n_landings)) then |
| 874 | 408 | n_landings(stats_thread()) = n_landings(stats_thread()) + 1 | |
| 875 | max_landing_resid(stats_thread()) = & | ||
| 876 | 408 | max(max_landing_resid(stats_thread()), resid) | |
| 877 | end if | ||
| 878 | 136 | end subroutine update_landing_stats | |
| 879 | |||
| 880 | 24 | subroutine sympl_landing_stats_reset | |
| 881 | integer :: nthreads | ||
| 882 | |||
| 883 | 24 | nthreads = 1 | |
| 884 | 24 | !$ nthreads = omp_get_max_threads() | |
| 885 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 17 times.
|
24 | if (allocated(n_landings)) then |
| 886 | ✗ | deallocate(n_landings, n_landing_evals, n_stops, & | |
| 887 | ✗ | n_fo_entries, n_fo_exits, & | |
| 888 | ✗ | n_fo_losses, n_fo_failures, n_fo_failure_status, & | |
| 889 |
9/18✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 7 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 7 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 7 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 7 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 7 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 7 times.
|
7 | n_bulk_recover, max_landing_resid) |
| 890 | end if | ||
| 891 | ✗ | allocate(n_landings(0:nthreads - 1), n_landing_evals(0:nthreads - 1), & | |
| 892 | ✗ | n_stops(0:nthreads - 1), n_fo_entries(0:nthreads - 1), & | |
| 893 | ✗ | n_fo_exits(0:nthreads - 1), n_fo_losses(0:nthreads - 1), & | |
| 894 | ✗ | n_fo_failures(0:nthreads - 1), & | |
| 895 | ✗ | n_fo_failure_status(5, 0:nthreads - 1), & | |
| 896 | ✗ | n_bulk_recover(3, 0:nthreads - 1), & | |
| 897 |
58/116✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 24 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 24 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 24 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 24 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 24 times.
✓ Branch 18 taken 24 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 24 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 24 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 24 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 24 times.
✓ Branch 30 taken 24 times.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✓ Branch 33 taken 24 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 24 times.
✗ Branch 37 not taken.
✓ Branch 38 taken 24 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 24 times.
✓ Branch 42 taken 24 times.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✓ Branch 45 taken 24 times.
✗ Branch 47 not taken.
✓ Branch 48 taken 24 times.
✗ Branch 49 not taken.
✓ Branch 50 taken 24 times.
✗ Branch 52 not taken.
✓ Branch 53 taken 24 times.
✓ Branch 54 taken 24 times.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✓ Branch 57 taken 24 times.
✗ Branch 59 not taken.
✓ Branch 60 taken 24 times.
✗ Branch 61 not taken.
✓ Branch 62 taken 24 times.
✗ Branch 64 not taken.
✓ Branch 65 taken 24 times.
✓ Branch 66 taken 24 times.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✓ Branch 69 taken 24 times.
✗ Branch 71 not taken.
✓ Branch 72 taken 24 times.
✗ Branch 73 not taken.
✓ Branch 74 taken 24 times.
✗ Branch 76 not taken.
✓ Branch 77 taken 24 times.
✓ Branch 78 taken 24 times.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✓ Branch 81 taken 24 times.
✗ Branch 83 not taken.
✓ Branch 84 taken 24 times.
✗ Branch 85 not taken.
✓ Branch 86 taken 24 times.
✗ Branch 88 not taken.
✓ Branch 89 taken 24 times.
✓ Branch 90 taken 24 times.
✗ Branch 91 not taken.
✓ Branch 92 taken 24 times.
✗ Branch 93 not taken.
✓ Branch 94 taken 24 times.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✓ Branch 97 taken 24 times.
✗ Branch 99 not taken.
✓ Branch 100 taken 24 times.
✗ Branch 101 not taken.
✓ Branch 102 taken 24 times.
✗ Branch 104 not taken.
✓ Branch 105 taken 24 times.
✓ Branch 106 taken 24 times.
✗ Branch 107 not taken.
✓ Branch 108 taken 24 times.
✗ Branch 109 not taken.
✓ Branch 110 taken 24 times.
✗ Branch 111 not taken.
✗ Branch 112 not taken.
✓ Branch 113 taken 24 times.
✗ Branch 115 not taken.
✓ Branch 116 taken 24 times.
✗ Branch 117 not taken.
✓ Branch 118 taken 24 times.
✗ Branch 120 not taken.
✓ Branch 121 taken 24 times.
✓ Branch 122 taken 24 times.
✗ Branch 123 not taken.
✓ Branch 124 taken 24 times.
✗ Branch 125 not taken.
✓ Branch 126 taken 24 times.
✗ Branch 127 not taken.
✗ Branch 128 not taken.
✓ Branch 129 taken 24 times.
✗ Branch 131 not taken.
✓ Branch 132 taken 24 times.
✗ Branch 134 not taken.
✓ Branch 135 taken 24 times.
|
360 | max_landing_resid(0:nthreads - 1)) |
| 898 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 27 times.
|
51 | n_landings = 0 |
| 899 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 27 times.
|
51 | n_landing_evals = 0 |
| 900 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 27 times.
|
51 | n_stops = 0 |
| 901 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 27 times.
|
51 | n_fo_entries = 0 |
| 902 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 27 times.
|
51 | n_fo_exits = 0 |
| 903 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 27 times.
|
51 | n_fo_losses = 0 |
| 904 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 27 times.
|
51 | n_fo_failures = 0 |
| 905 |
4/4✓ Branch 0 taken 24 times.
✓ Branch 1 taken 27 times.
✓ Branch 2 taken 135 times.
✓ Branch 3 taken 27 times.
|
186 | n_fo_failure_status = 0 |
| 906 |
4/4✓ Branch 0 taken 24 times.
✓ Branch 1 taken 27 times.
✓ Branch 2 taken 81 times.
✓ Branch 3 taken 27 times.
|
132 | n_bulk_recover = 0 |
| 907 |
2/2✓ Branch 0 taken 27 times.
✓ Branch 1 taken 24 times.
|
51 | max_landing_resid = 0.0_dp |
| 908 | 24 | end subroutine sympl_landing_stats_reset | |
| 909 | |||
| 910 | 2 | subroutine sympl_landing_stats(landings, max_resid, stops) | |
| 911 | integer, intent(out) :: landings, stops | ||
| 912 | real(dp), intent(out) :: max_resid | ||
| 913 | |||
| 914 | 2 | landings = 0 | |
| 915 | 2 | max_resid = 0.0_dp | |
| 916 | 2 | stops = 0 | |
| 917 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (allocated(n_landings)) then |
| 918 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
|
4 | landings = sum(n_landings) |
| 919 |
5/10✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
|
6 | max_resid = maxval(max_landing_resid) |
| 920 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
|
4 | stops = sum(n_stops) |
| 921 | end if | ||
| 922 | 2 | end subroutine sympl_landing_stats | |
| 923 | |||
| 924 | ✗ | subroutine sympl_landing_eval_stats(evals) | |
| 925 | !> Total implicit-step evaluations spent inside locate_landing since | ||
| 926 | !> the last reset; evaluations/landing is the root-iteration cost the | ||
| 927 | !> seeded bracket shrink is measured by. | ||
| 928 | integer, intent(out) :: evals | ||
| 929 | |||
| 930 | ✗ | evals = 0 | |
| 931 | ✗ | if (allocated(n_landing_evals)) evals = sum(n_landing_evals) | |
| 932 | ✗ | end subroutine sympl_landing_eval_stats | |
| 933 | |||
| 934 | |||
| 935 | 16 | subroutine count_fo_entry | |
| 936 |
1/2✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
|
16 | if (allocated(n_fo_entries)) n_fo_entries(stats_thread()) = & |
| 937 | 48 | n_fo_entries(stats_thread()) + 1 | |
| 938 | 16 | end subroutine count_fo_entry | |
| 939 | |||
| 940 | ✗ | subroutine count_fo_exit | |
| 941 | ✗ | if (allocated(n_fo_exits)) n_fo_exits(stats_thread()) = & | |
| 942 | ✗ | n_fo_exits(stats_thread()) + 1 | |
| 943 | ✗ | end subroutine count_fo_exit | |
| 944 | |||
| 945 | ✗ | subroutine count_fo_failure(status) | |
| 946 | integer, intent(in) :: status | ||
| 947 | integer :: tid | ||
| 948 | |||
| 949 | ✗ | if (.not. allocated(n_fo_failures)) return | |
| 950 | ✗ | tid = stats_thread() | |
| 951 | ✗ | n_fo_failures(tid) = n_fo_failures(tid) + 1 | |
| 952 | ✗ | if (status >= 1 .and. status <= size(n_fo_failure_status, 1)) & | |
| 953 | n_fo_failure_status(status, tid) = & | ||
| 954 | ✗ | n_fo_failure_status(status, tid) + 1 | |
| 955 | end subroutine count_fo_failure | ||
| 956 | |||
| 957 | 10009 | subroutine sympl_fo_stats(entries, exits, losses, failures, failure_status) | |
| 958 | integer, intent(out) :: entries, exits, losses, failures | ||
| 959 | integer, intent(out) :: failure_status(5) | ||
| 960 | |||
| 961 | 10009 | entries = 0 | |
| 962 | 10009 | exits = 0 | |
| 963 | 10009 | losses = 0 | |
| 964 | 10009 | failures = 0 | |
| 965 | 10009 | failure_status = 0 | |
| 966 |
1/2✓ Branch 0 taken 10009 times.
✗ Branch 1 not taken.
|
10009 | if (allocated(n_fo_entries)) then |
| 967 |
2/2✓ Branch 0 taken 10009 times.
✓ Branch 1 taken 10009 times.
|
20018 | entries = sum(n_fo_entries) |
| 968 |
2/2✓ Branch 0 taken 10009 times.
✓ Branch 1 taken 10009 times.
|
20018 | exits = sum(n_fo_exits) |
| 969 |
2/2✓ Branch 0 taken 10009 times.
✓ Branch 1 taken 10009 times.
|
20018 | losses = sum(n_fo_losses) |
| 970 |
2/2✓ Branch 0 taken 10009 times.
✓ Branch 1 taken 10009 times.
|
20018 | failures = sum(n_fo_failures) |
| 971 |
4/4✓ Branch 0 taken 50045 times.
✓ Branch 1 taken 10009 times.
✓ Branch 2 taken 50045 times.
✓ Branch 3 taken 50045 times.
|
110099 | failure_status = sum(n_fo_failure_status, dim=2) |
| 972 | end if | ||
| 973 | 10009 | end subroutine sympl_fo_stats | |
| 974 | |||
| 975 | ✗ | subroutine count_bulk_recovery(kind) | |
| 976 | integer, intent(in) :: kind | ||
| 977 | |||
| 978 | ✗ | if (allocated(n_bulk_recover)) n_bulk_recover(kind, stats_thread()) = & | |
| 979 | ✗ | n_bulk_recover(kind, stats_thread()) + 1 | |
| 980 | ✗ | end subroutine count_bulk_recovery | |
| 981 | |||
| 982 | ✗ | subroutine sympl_bulk_recovery_stats(counts) | |
| 983 | !> Bulk recoveries by cause: iteration cap, teleport, Newton. | ||
| 984 | integer, intent(out) :: counts(3) | ||
| 985 | |||
| 986 | ✗ | counts = 0 | |
| 987 | ✗ | if (allocated(n_bulk_recover)) counts = sum(n_bulk_recover, dim=2) | |
| 988 | ✗ | end subroutine sympl_bulk_recovery_stats | |
| 989 | |||
| 990 | 864 | integer function stats_thread() | |
| 991 | 864 | stats_thread = 0 | |
| 992 |
0/2✗ Branch 17 not taken.
✗ Branch 18 not taken.
|
296 | !$ stats_thread = omp_get_thread_num() |
| 993 | end function stats_thread | ||
| 994 | |||
| 995 | ✗ | end module spectre_sympl_orbit | |
| 996 |