| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module simple_main | ||
| 2 | use, intrinsic :: iso_fortran_env, only: int8 | ||
| 3 | !$ use omp_lib | ||
| 4 | use timing, only: get_wtime | ||
| 5 | use util, only: sqrt2, twopi | ||
| 6 | use simple, only: init_vmec, init_sympl, init_fo, orbit_timestep_fo, & | ||
| 7 | orbit_timestep_fo_bridge, reseed_sympl, tracer_t, ORBIT_FO_LOSS, & | ||
| 8 | ORBIT_FO_NUMERICAL | ||
| 9 | use diag_mod, only: icounter | ||
| 10 | use collis_alp, only: loacol_alpha, stost, init_collision_profiles, & | ||
| 11 | configure_lorentz_nustar, stost_lorentz | ||
| 12 | use magfie_sub, only: ALBERT, SPECTRE | ||
| 13 | use samplers, only: sample | ||
| 14 | use field_can_mod, only: field_can_t, evaluate_field => evaluate, & | ||
| 15 | integ_to_ref, ref_to_integ, init_field_can | ||
| 16 | use callback, only: output_orbits_macrostep | ||
| 17 | use params, only: swcoll, collision_model, ntestpart, startmode, & | ||
| 18 | special_ants_file, num_surf, & | ||
| 19 | grid_density, dtau, dtaumin, ntau, v0, kpart, confpart_pass, & | ||
| 20 | confpart_trap, unresolved_orbits, times_lost, orbit_exit_code, & | ||
| 21 | boundary_event_radial_residual, & | ||
| 22 | boundary_event_time_width, integmode, relerr, trace_time, class_plot, & | ||
| 23 | fast_class, generate_start_only, ntcut, iclass, bmin, bmax, zstart, & | ||
| 24 | zend, trap_par, perp_inv, class_scores, sbeg, ntimstep, & | ||
| 25 | should_skip, & | ||
| 26 | reset_seed_if_deterministic, particle_stream_id, & | ||
| 27 | classification_enabled, field_input, & | ||
| 28 | isw_field_type, reuse_batch, & | ||
| 29 | max_consecutive_warning_holds, & | ||
| 30 | coord_input, wall_input, wall_units, wall_hit, wall_hit_triangle, & | ||
| 31 | wall_hit_cart, & | ||
| 32 | wall_query_rho_lcfs, & | ||
| 33 | chart_boundary_kind, chart_boundary_kind_effective, & | ||
| 34 | wall_hit_normal_cart, wall_hit_cos_incidence, wall_hit_angle_rad, & | ||
| 35 | ntau_macro, kt_macro, checkpoint_interval, orbit_model, orbit_coord, & | ||
| 36 | gpu_mode, gpu_backend, gpu_method, gpu_start_coordinates, & | ||
| 37 | gpu_landreman, gpu_compact_init, gpu_t_block, gpu_loss_tau, & | ||
| 38 | gpu_maxloss, gpu_min_timestep, gpu_particle_profile, & | ||
| 39 | ORBIT_GC, ORBIT_FULL_ORBIT, ORBIT_EXIT_COMPLETED, ORBIT_EXIT_LCFS, & | ||
| 40 | ORBIT_EXIT_WALL, ORBIT_EXIT_SKIPPED, ORBIT_EXIT_NUMERICAL_CONFINED, & | ||
| 41 | ORBIT_EXIT_NUMERICAL_DOMAIN, & | ||
| 42 | ORBIT_EXIT_NUMERICAL_MAXITER, ORBIT_EXIT_NUMERICAL_LINEAR, & | ||
| 43 | ORBIT_EXIT_NUMERICAL_EVENT, ORBIT_EXIT_NUMERICAL_FULL_ORBIT | ||
| 44 | use params, only: radial_electric_potential_slope | ||
| 45 | use params, only: canonical_grid_nr, canonical_grid_ntheta, & | ||
| 46 | canonical_grid_nphi, canonical_ode_relerr | ||
| 47 | use diag_counters, only: diag_counters_init, count_event, & | ||
| 48 | EVT_WARNING_STEP_SKIP, EVT_SYMPLECTIC_RK_RECOVERY, & | ||
| 49 | EVT_SYMPLECTIC_RESUME, EVT_TOROIDAL_REGULARIZATION, & | ||
| 50 | EVT_AXIS_FO_BRIDGE | ||
| 51 | use progress_monitor, only: progress_init, progress_tick, progress_finalize | ||
| 52 | use restart_mod, only: particle_done, read_restart_data, restore_confined_counts | ||
| 53 | use chartmap_metadata, only: chartmap_metadata_t, read_chartmap_metadata | ||
| 54 | use reference_coordinates, only: ref_coords | ||
| 55 | use stl_wall_intersection, only: stl_wall_t, stl_wall_init, & | ||
| 56 | stl_wall_finalize, & | ||
| 57 | stl_wall_first_hit_segment_with_normal | ||
| 58 | use libneo_coordinates, only: chartmap_coordinate_system_t | ||
| 59 | use orbit_symplectic_base, only: symplectic_integrator_t, & | ||
| 60 | SYMPLECTIC_STEP_BOUNDARY, & | ||
| 61 | SYMPLECTIC_STEP_OUTSIDE_DOMAIN, SYMPLECTIC_STEP_MAXITER, & | ||
| 62 | SYMPLECTIC_STEP_LINEAR_SOLVE, SYMPLECTIC_STEP_EVENT_NOT_CONVERGED, & | ||
| 63 | SYMPLECTIC_STEP_BOUNDARY_LIMITED, symplectic_newton_warning_mode | ||
| 64 | |||
| 65 | implicit none | ||
| 66 | |||
| 67 | ! Define real(dp) kind parameter | ||
| 68 | integer, parameter :: dp = kind(1.0d0) | ||
| 69 | |||
| 70 | logical, save :: wall_enabled = .false. | ||
| 71 | integer, save :: map_boundary_exit_code = ORBIT_EXIT_LCFS | ||
| 72 | real(dp), save :: chartmap_cart_scale_to_m = -1.0d0 | ||
| 73 | real(dp), save :: invariant_edge_radius = 1.0_dp | ||
| 74 | type(stl_wall_t), save :: wall | ||
| 75 | |||
| 76 | integer, parameter :: INPUT_VMEC = 0 | ||
| 77 | integer, parameter :: INPUT_TEST = 1 | ||
| 78 | integer, parameter :: INPUT_SPECTRE = 2 | ||
| 79 | integer, parameter :: INPUT_CHARTMAP = 3 | ||
| 80 | integer, parameter :: INPUT_VMECPP = 4 | ||
| 81 | |||
| 82 | integer, parameter :: RK_RECOVERY_GENERIC = 1 | ||
| 83 | integer, parameter :: RK_RECOVERY_AXIS = 2 | ||
| 84 | integer, parameter :: RK_RECOVERY_EDGE = 3 | ||
| 85 | real(dp), parameter :: RK_AXIS_ENTER = 0.1_dp | ||
| 86 | real(dp), parameter :: RK_AXIS_EXIT = 0.12_dp | ||
| 87 | real(dp), parameter :: FO_AXIS_ENTER = 0.02_dp | ||
| 88 | real(dp), parameter :: RK_EDGE_ENTER = 0.98_dp | ||
| 89 | integer, parameter :: RK_RECOVERY_STEP_LIMIT = 10000 | ||
| 90 | integer, parameter :: TOROIDAL_RECOVERY_STEP_LIMIT = 100000 | ||
| 91 | integer, parameter :: RK_SAFE_STREAK_TO_RESUME = 256 | ||
| 92 | real(dp), parameter :: RK_RECOVERY_AXIS_ENTER = 0.1_dp | ||
| 93 | real(dp), parameter :: NUMERICAL_CONFINED_S_MAX = 0.01_dp | ||
| 94 | ! Retry only a failed legacy axis solve with a smoother, still axis-local | ||
| 95 | ! Cartesian regularization. Successful release-era RK paths stay unchanged. | ||
| 96 | real(dp), parameter :: RK_RECOVERY_AXIS_FLOOR = 1.0e-3_dp | ||
| 97 | real(dp), parameter :: FO_BRIDGE_BSTAR_ENTER = 0.2_dp | ||
| 98 | real(dp), parameter :: RK_RECOVERY_BSTAR_EXIT = 0.3_dp | ||
| 99 | integer, parameter :: RECOVERY_STEP_STANDARD = 0 | ||
| 100 | integer, parameter :: RECOVERY_STEP_AXIS = 1 | ||
| 101 | integer, parameter :: RECOVERY_STEP_TOROIDAL = 2 | ||
| 102 | integer, parameter :: RECOVERY_STEP_FULL_ORBIT = 3 | ||
| 103 | real(dp), parameter :: MOMENTUM_SANITY_FRACTION = 0.005_dp | ||
| 104 | |||
| 105 | type :: rk_recovery_state_t | ||
| 106 | logical :: active = .false. | ||
| 107 | logical :: has_momentum_reference = .false. | ||
| 108 | integer :: reason = RK_RECOVERY_GENERIC | ||
| 109 | integer :: steps = 0 | ||
| 110 | integer :: safe_steps = 0 | ||
| 111 | integer :: failures = 0 | ||
| 112 | real(dp) :: momentum_reference = 0.0_dp | ||
| 113 | end type rk_recovery_state_t | ||
| 114 | |||
| 115 | contains | ||
| 116 | |||
| 117 | 15 | logical function recovery_state_has_valid_phase_space(z) | |
| 118 | use finite_check, only: is_finite | ||
| 119 | |||
| 120 | real(dp), intent(in) :: z(5) | ||
| 121 | |||
| 122 | recovery_state_has_valid_phase_space = all(is_finite(z)) .and. & | ||
| 123 | z(4) > 0.0_dp .and. & | ||
| 124 |
7/10✓ Branch 0 taken 75 times.
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 75 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 15 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✓ Branch 9 taken 13 times.
|
90 | abs(z(5)) <= 1.0_dp + 100.0_dp*epsilon(1.0_dp) |
| 125 | 15 | end function recovery_state_has_valid_phase_space | |
| 126 | |||
| 127 | 10 | logical function recovery_state_is_physical(z, expected_momentum) | |
| 128 | real(dp), intent(in) :: z(5), expected_momentum | ||
| 129 | |||
| 130 | recovery_state_is_physical = recovery_state_has_valid_phase_space(z) .and. & | ||
| 131 | abs(z(4) - expected_momentum) <= & | ||
| 132 |
4/4✓ Branch 0 taken 9 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 7 times.
|
10 | MOMENTUM_SANITY_FRACTION*max(abs(expected_momentum), tiny(1.0_dp)) |
| 133 | 10 | end function recovery_state_is_physical | |
| 134 | |||
| 135 | 7 | subroutine orbit_timestep_recovery(fo, z, interval, tolerance, & | |
| 136 | expected_momentum, ierr, method) | ||
| 137 | use alpha_lifetime_sub, only: orbit_timestep_axis, & | ||
| 138 | bstar_parallel_factor, orbit_timestep_toroidal_regularized | ||
| 139 | use orbit_fo_boris, only: fo_state_t | ||
| 140 | |||
| 141 | type(fo_state_t), intent(inout) :: fo | ||
| 142 | real(dp), intent(inout) :: z(5) | ||
| 143 | real(dp), intent(in) :: interval, tolerance, expected_momentum | ||
| 144 | integer, intent(out) :: ierr | ||
| 145 | integer, intent(out) :: method | ||
| 146 | real(dp) :: z_start(5) | ||
| 147 | real(dp) :: bstar_factor | ||
| 148 | real(dp) :: x_reference(3) | ||
| 149 | |||
| 150 | 7 | method = RECOVERY_STEP_STANDARD | |
| 151 | 7 | z_start = z | |
| 152 | call orbit_timestep_axis(z, interval, interval, tolerance, ierr, & | ||
| 153 | 7 | RK_RECOVERY_STEP_LIMIT) | |
| 154 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
|
7 | if (ierr == 0) then |
| 155 | ✗ | if (.not. recovery_state_is_physical(z, expected_momentum)) ierr = 2 | |
| 156 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 4 times.
|
7 | else if (ierr == 1) then |
| 157 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if (.not. recovery_state_has_valid_phase_space(z)) ierr = 2 |
| 158 | end if | ||
| 159 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 3 times.
|
7 | if (ierr /= 2) return |
| 160 | |||
| 161 | 4 | call bstar_parallel_factor(z_start, bstar_factor) | |
| 162 | 4 | call integ_to_ref(z_start(1:3), x_reference) | |
| 163 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
4 | if (x_reference(1) < RK_RECOVERY_AXIS_ENTER .and. & |
| 164 | abs(bstar_factor) >= FO_BRIDGE_BSTAR_ENTER) then | ||
| 165 | ✗ | z = z_start | |
| 166 | call orbit_timestep_axis(z, interval, interval, tolerance, ierr, & | ||
| 167 | RK_RECOVERY_STEP_LIMIT, RK_RECOVERY_AXIS_FLOOR, & | ||
| 168 | ✗ | RK_RECOVERY_AXIS_ENTER) | |
| 169 | ✗ | if (ierr == 0) then | |
| 170 | ✗ | if (.not. recovery_state_is_physical(z, expected_momentum)) & | |
| 171 | ✗ | ierr = 2 | |
| 172 | ✗ | else if (ierr == 1) then | |
| 173 | ✗ | if (.not. recovery_state_has_valid_phase_space(z)) ierr = 2 | |
| 174 | end if | ||
| 175 | ✗ | if (ierr /= 2) then | |
| 176 | ✗ | method = RECOVERY_STEP_AXIS | |
| 177 | ✗ | return | |
| 178 | end if | ||
| 179 | end if | ||
| 180 | |||
| 181 | 4 | z = z_start | |
| 182 | call orbit_timestep_toroidal_regularized(z, interval, & | ||
| 183 | max(tolerance, 1.0e-10_dp), ierr, & | ||
| 184 | 4 | TOROIDAL_RECOVERY_STEP_LIMIT) | |
| 185 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (ierr == 0) then |
| 186 | ✗ | if (recovery_state_is_physical(z, expected_momentum)) then | |
| 187 | ✗ | method = RECOVERY_STEP_TOROIDAL | |
| 188 | ✗ | call count_event(EVT_TOROIDAL_REGULARIZATION) | |
| 189 | ✗ | return | |
| 190 | end if | ||
| 191 | ✗ | ierr = 2 | |
| 192 | end if | ||
| 193 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (ierr == 1) then |
| 194 | ✗ | if (.not. recovery_state_has_valid_phase_space(z)) ierr = 2 | |
| 195 | end if | ||
| 196 |
2/4✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
|
4 | if (ierr /= 2 .or. x_reference(1) >= FO_AXIS_ENTER) return |
| 197 | |||
| 198 | ✗ | z = z_start | |
| 199 | ✗ | call orbit_timestep_fo_bridge(fo, z, interval, ierr) | |
| 200 | ✗ | if (ierr == 0) then | |
| 201 | ✗ | if (recovery_state_is_physical(z, expected_momentum)) then | |
| 202 | ✗ | method = RECOVERY_STEP_FULL_ORBIT | |
| 203 | ✗ | call count_event(EVT_AXIS_FO_BRIDGE) | |
| 204 | else | ||
| 205 | ✗ | ierr = 2 | |
| 206 | end if | ||
| 207 | end if | ||
| 208 | 7 | end subroutine orbit_timestep_recovery | |
| 209 | |||
| 210 | 4 | subroutine activate_rk_recovery(state, z, momentum_reference) | |
| 211 | type(rk_recovery_state_t), intent(inout) :: state | ||
| 212 | real(dp), intent(in) :: z(5), momentum_reference | ||
| 213 | real(dp) :: x_reference(3) | ||
| 214 | |||
| 215 | 4 | call integ_to_ref(z(1:3), x_reference) | |
| 216 | 4 | state%active = .true. | |
| 217 | 4 | state%has_momentum_reference = .true. | |
| 218 | 4 | state%steps = 0 | |
| 219 | 4 | state%safe_steps = 0 | |
| 220 | 4 | state%momentum_reference = momentum_reference | |
| 221 | 4 | state%failures = min(state%failures + 1, 9) | |
| 222 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
|
4 | if (x_reference(1) < RK_AXIS_ENTER) then |
| 223 | 1 | state%reason = RK_RECOVERY_AXIS | |
| 224 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2 times.
|
3 | else if (x_reference(1) > RK_EDGE_ENTER) then |
| 225 | 1 | state%reason = RK_RECOVERY_EDGE | |
| 226 | else | ||
| 227 | 2 | state%reason = RK_RECOVERY_GENERIC | |
| 228 | end if | ||
| 229 | 7 | end subroutine activate_rk_recovery | |
| 230 | |||
| 231 | 8 | logical function should_resume_symplectic(state, z) | |
| 232 | type(rk_recovery_state_t), intent(in) :: state | ||
| 233 | real(dp), intent(in) :: z(5) | ||
| 234 | integer :: cooldown | ||
| 235 | real(dp) :: x_reference(3) | ||
| 236 | |||
| 237 | 8 | call integ_to_ref(z(1:3), x_reference) | |
| 238 | 10 | select case (state%reason) | |
| 239 | case (RK_RECOVERY_AXIS) | ||
| 240 | should_resume_symplectic = x_reference(1) >= RK_AXIS_EXIT .and. & | ||
| 241 |
3/4✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 1 time.
|
2 | state%safe_steps >= RK_SAFE_STREAK_TO_RESUME |
| 242 | case (RK_RECOVERY_EDGE) | ||
| 243 | should_resume_symplectic = x_reference(1) <= RK_EDGE_ENTER .and. & | ||
| 244 |
3/4✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
|
2 | state%safe_steps >= RK_SAFE_STREAK_TO_RESUME |
| 245 | case default | ||
| 246 |
2/4✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
|
4 | cooldown = 2**(7 + state%failures) |
| 247 |
3/3✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 4 times.
|
12 | should_resume_symplectic = state%safe_steps >= cooldown |
| 248 | end select | ||
| 249 | 8 | end function should_resume_symplectic | |
| 250 | |||
| 251 | ✗ | subroutine update_recovery_streak(state, z, method) | |
| 252 | use alpha_lifetime_sub, only: bstar_parallel_factor | ||
| 253 | |||
| 254 | type(rk_recovery_state_t), intent(inout) :: state | ||
| 255 | real(dp), intent(in) :: z(5) | ||
| 256 | integer, intent(in) :: method | ||
| 257 | real(dp) :: bstar_factor | ||
| 258 | |||
| 259 | ✗ | state%steps = state%steps + 1 | |
| 260 | ✗ | if (method /= RECOVERY_STEP_STANDARD) then | |
| 261 | ✗ | state%safe_steps = 0 | |
| 262 | ✗ | return | |
| 263 | end if | ||
| 264 | ✗ | call bstar_parallel_factor(z, bstar_factor) | |
| 265 | ✗ | if (abs(bstar_factor) < RK_RECOVERY_BSTAR_EXIT) then | |
| 266 | ✗ | state%safe_steps = 0 | |
| 267 | ✗ | return | |
| 268 | end if | ||
| 269 | ✗ | state%safe_steps = state%safe_steps + 1 | |
| 270 | end subroutine update_recovery_streak | ||
| 271 | |||
| 272 | 65 | subroutine main | |
| 273 | use params, only: restart | ||
| 274 | use timing, only: init_timer, print_phase_time | ||
| 275 | use version, only: simple_version | ||
| 276 | |||
| 277 | character(256) :: config_file | ||
| 278 | type(tracer_t) :: norb | ||
| 279 | |||
| 280 | 65 | print '(A,A)', 'SIMPLE version ', simple_version | |
| 281 | 65 | call init_timer() | |
| 282 | |||
| 283 |
2/2✓ Branch 1 taken 33 times.
✓ Branch 2 taken 32 times.
|
65 | if (command_argument_count() == 0) then |
| 284 | 33 | config_file = 'simple.in' | |
| 285 | else | ||
| 286 | 32 | call get_command_argument(1, config_file) | |
| 287 | end if | ||
| 288 | |||
| 289 | 65 | call load_config(config_file) | |
| 290 | 57 | call init_run(norb) | |
| 291 | 55 | call init_sampling_field | |
| 292 | 55 | call sample_start | |
| 293 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 50 times.
|
55 | if (generate_start_only) stop 'stopping after generating start.dat' |
| 294 | 50 | call init_tracing_field | |
| 295 | 50 | call init_counters | |
| 296 | 50 | call print_phase_time('Counter initialization completed') | |
| 297 | |||
| 298 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 49 times.
|
50 | if (restart) then |
| 299 | 1 | call read_restart_data() | |
| 300 | 1 | call restore_confined_counts() | |
| 301 | 1 | call print_phase_time('Restart data loaded') | |
| 302 | end if | ||
| 303 | |||
| 304 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 48 times.
|
50 | if (trim(gpu_mode) == 'benchmark') then |
| 305 | 2 | call trace_compare_gpu(norb) | |
| 306 | return | ||
| 307 | end if | ||
| 308 | |||
| 309 | 48 | call diag_counters_init | |
| 310 |
2/2✓ Branch 0 taken 15 times.
✓ Branch 1 taken 33 times.
|
48 | if (isw_field_type == SPECTRE) then |
| 311 | block | ||
| 312 | use spectre_sympl_orbit, only: sympl_landing_stats_reset | ||
| 313 | 15 | call sympl_landing_stats_reset | |
| 314 | end block | ||
| 315 | end if | ||
| 316 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
|
48 | if (trim(gpu_mode) == 'production') then |
| 317 | ✗ | call trace_gpu_production(norb) | |
| 318 | else | ||
| 319 | 48 | call progress_init(checkpoint_interval, ntestpart) | |
| 320 | 48 | call trace_parallel(norb) | |
| 321 | 48 | call progress_finalize | |
| 322 | end if | ||
| 323 | 48 | call print_phase_time('Particle tracing completed') | |
| 324 | |||
| 325 | 48 | call write_output | |
| 326 | 48 | call print_phase_time('Output writing completed') | |
| 327 | |||
| 328 | 48 | call stl_wall_finalize(wall) | |
| 329 |
12/12✓ Branch 0 taken 195 times.
✓ Branch 1 taken 65 times.
✓ Branch 2 taken 390 times.
✓ Branch 3 taken 65 times.
✓ Branch 4 taken 4160 times.
✓ Branch 5 taken 65 times.
✓ Branch 6 taken 195 times.
✓ Branch 7 taken 65 times.
✓ Branch 8 taken 195 times.
✓ Branch 9 taken 65 times.
✓ Branch 10 taken 195 times.
✓ Branch 11 taken 65 times.
|
5445 | end subroutine main |
| 330 | |||
| 331 | !> Read and validate a simple.in file, including its &profiles group. | ||
| 332 | 65 | subroutine load_config(config_file) | |
| 333 | 65 | use params, only: read_config | |
| 334 | use timing, only: print_phase_time | ||
| 335 | |||
| 336 | character(*), intent(in) :: config_file | ||
| 337 | |||
| 338 | 65 | call read_config(config_file) | |
| 339 | 58 | call validate_orbit_model_config | |
| 340 | 58 | call validate_gpu_production_config | |
| 341 | 57 | call print_phase_time('Configuration reading completed') | |
| 342 | 57 | call read_profiles_config(config_file) | |
| 343 | 57 | call print_phase_time('Profiles configuration reading completed') | |
| 344 | 57 | end subroutine load_config | |
| 345 | |||
| 346 | !> Load the field and derive run parameters from the current &config values. | ||
| 347 | 69 | subroutine init_run(norb) | |
| 348 | 65 | use params, only: netcdffile, ns_s, ns_tp, multharm, params_init | |
| 349 | use timing, only: print_phase_time | ||
| 350 | |||
| 351 | type(tracer_t), intent(inout) :: norb | ||
| 352 | |||
| 353 | 69 | call init_field(norb, netcdffile, ns_s, ns_tp, multharm, integmode) | |
| 354 | 67 | call print_phase_time('Field initialization completed') | |
| 355 | 67 | call params_init | |
| 356 | 67 | call print_phase_time('Parameter initialization completed') | |
| 357 | 67 | call print_parameters | |
| 358 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 67 times.
|
67 | if (swcoll) then |
| 359 | ✗ | call init_collisions | |
| 360 | ✗ | call print_phase_time('Collision initialization completed') | |
| 361 | end if | ||
| 362 | 67 | end subroutine init_run | |
| 363 | |||
| 364 | !> Classify the configured input as TEST, SPECTRE, Boozer chartmap, VMEC or | ||
| 365 | !> direct VMEC++ (a .json input). | ||
| 366 | 325 | integer function run_input_kind() | |
| 367 | 69 | use magfie_sub, only: TEST | |
| 368 | use field_boozer_chartmap, only: is_boozer_chartmap | ||
| 369 | use field, only: is_spectre_file | ||
| 370 | use params, only: netcdffile | ||
| 371 | |||
| 372 | 325 | run_input_kind = INPUT_VMEC | |
| 373 |
2/2✓ Branch 0 taken 315 times.
✓ Branch 1 taken 10 times.
|
325 | if (isw_field_type == TEST) then |
| 374 | run_input_kind = INPUT_TEST | ||
| 375 |
1/2✓ Branch 0 taken 315 times.
✗ Branch 1 not taken.
|
315 | else if (is_vmecpp_input(netcdffile)) then |
| 376 | run_input_kind = INPUT_VMECPP | ||
| 377 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 315 times.
|
315 | else if (len_trim(field_input) > 0) then |
| 378 |
2/2✓ Branch 1 taken 229 times.
✓ Branch 2 taken 86 times.
|
315 | if (is_spectre_file(field_input)) then |
| 379 | run_input_kind = INPUT_SPECTRE | ||
| 380 |
2/2✓ Branch 1 taken 199 times.
✓ Branch 2 taken 30 times.
|
229 | else if (is_boozer_chartmap(field_input)) then |
| 381 | 325 | run_input_kind = INPUT_CHARTMAP | |
| 382 | end if | ||
| 383 | end if | ||
| 384 | 325 | end function run_input_kind | |
| 385 | |||
| 386 | !> Select the field used for sampling and integrate the starting field line. | ||
| 387 | 80 | subroutine init_sampling_field | |
| 388 | 325 | use magfie_sub, only: TEST, VMEC, BOOZER, SPECTRE, init_magfie | |
| 389 | #ifdef SIMPLE_ENABLE_VMECPP | ||
| 390 | use magfie_sub, only: VMECPP | ||
| 391 | #endif | ||
| 392 | use samplers, only: init_starting_surf | ||
| 393 | use timing, only: print_phase_time | ||
| 394 | |||
| 395 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 62 times.
|
80 | if (startmode == 7) then |
| 396 |
2/4✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 18 times.
|
18 | if (run_input_kind() == INPUT_TEST .or. & |
| 397 | run_input_kind() == INPUT_SPECTRE) & | ||
| 398 | ✗ | error stop 'fusion births require VMEC or a Boozer chartmap' | |
| 399 |
3/4✓ Branch 1 taken 3 times.
✓ Branch 2 taken 15 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
|
18 | if (run_input_kind() == INPUT_CHARTMAP .and. isw_field_type /= BOOZER) & |
| 400 | ✗ | error stop 'fusion births on a chartmap require Boozer coordinates' | |
| 401 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
|
18 | if (gpu_compact_init) & |
| 402 | ✗ | error stop 'fusion births require gpu_compact_init=.false.' | |
| 403 | end if | ||
| 404 | 2 | select case (run_input_kind()) | |
| 405 | case (INPUT_TEST) | ||
| 406 | 2 | call init_magfie(TEST) | |
| 407 | case (INPUT_SPECTRE) | ||
| 408 | 21 | call init_magfie(SPECTRE) | |
| 409 | case (INPUT_CHARTMAP) | ||
| 410 |
1/2✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
|
7 | if (.not. gpu_compact_init) then |
| 411 | 7 | call init_magfie(isw_field_type) | |
| 412 | 7 | call init_starting_surf | |
| 413 | end if | ||
| 414 | case (INPUT_VMECPP) | ||
| 415 | #ifdef SIMPLE_ENABLE_VMECPP | ||
| 416 | call init_magfie(VMECPP) | ||
| 417 | #else | ||
| 418 | ✗ | error stop 'SIMPLE was built without SIMPLE_ENABLE_VMECPP' | |
| 419 | #endif | ||
| 420 | 50 | call init_starting_surf | |
| 421 | case default | ||
| 422 | 50 | call init_magfie(VMEC) | |
| 423 |
4/5✓ Branch 1 taken 2 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 7 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 50 times.
|
130 | call init_starting_surf |
| 424 | end select | ||
| 425 | 80 | call print_phase_time('Sampling field initialization completed') | |
| 426 | 80 | end subroutine init_sampling_field | |
| 427 | |||
| 428 | !> Fill zstart according to startmode. Requires init_sampling_field. | ||
| 429 | 68 | subroutine sample_start | |
| 430 | 80 | use samplers, only: sample_spectre_surface, init_spectre_start_bounds | |
| 431 | use timing, only: print_phase_time | ||
| 432 | |||
| 433 | 2 | select case (run_input_kind()) | |
| 434 | case (INPUT_TEST) | ||
| 435 | 2 | call sample_particles_test_field | |
| 436 | case (INPUT_SPECTRE) | ||
| 437 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 1 time.
|
19 | if (startmode == 1) then |
| 438 | 18 | call sample_spectre_surface(zstart) | |
| 439 | else | ||
| 440 | 1 | call sample_particles(.true.) | |
| 441 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (startmode == 2) call init_spectre_start_bounds(zstart) |
| 442 | end if | ||
| 443 | case (INPUT_CHARTMAP) | ||
| 444 | 5 | call sample_particles(.true.) | |
| 445 | case default | ||
| 446 |
4/4✓ Branch 1 taken 2 times.
✓ Branch 2 taken 19 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 42 times.
|
68 | call sample_particles |
| 447 | end select | ||
| 448 | 67 | call print_phase_time('Particle sampling completed') | |
| 449 | 67 | end subroutine sample_start | |
| 450 | |||
| 451 | !> Select the integration field and the bmin/bmax cache for tracing. | ||
| 452 | 52 | subroutine init_tracing_field | |
| 453 | use magfie_sub, only: TEST, init_magfie, spectre_field | ||
| 454 | use timing, only: print_phase_time | ||
| 455 | |||
| 456 | ! Library callers can activate tracing after sampling-only initialization. | ||
| 457 |
2/2✓ Branch 1 taken 16 times.
✓ Branch 2 taken 36 times.
|
52 | if (run_input_kind() == INPUT_SPECTRE) then |
| 458 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 1 time.
|
16 | if (spectre_field%data%Mvol > 1 .and. integmode <= 0) error stop & |
| 459 | 'multi-volume SPECTRE tracing requires integmode > 0 and crossing_level = 3; '// & | ||
| 460 | 1 | 'see DOC/spectre-rejected-interface-models.md' | |
| 461 | end if | ||
| 462 |
3/4✓ Branch 1 taken 20 times.
✓ Branch 2 taken 31 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 20 times.
|
51 | if (run_input_kind() == INPUT_VMEC .or. run_input_kind() == INPUT_VMECPP) then |
| 463 | 31 | call init_magfie(isw_field_type) | |
| 464 | 31 | call print_phase_time('Field type initialization completed') | |
| 465 | end if | ||
| 466 |
4/4✓ Branch 0 taken 49 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 7 times.
✓ Branch 3 taken 42 times.
|
51 | if ((isw_field_type /= TEST) .and. needs_bminmax_cache()) then |
| 467 | 7 | call init_bminmax | |
| 468 | 7 | call print_phase_time('Bmin/Bmax initialization completed') | |
| 469 | end if | ||
| 470 | 51 | end subroutine init_tracing_field | |
| 471 | |||
| 472 | ! Reject orbit_model values this build does not implement, with a clear message, | ||
| 473 | ! before any tracing starts. Only guiding-center (the default symplectic path) | ||
| 474 | ! and full orbit (the gyro-resolved Boris pusher) are available here. | ||
| 475 | 70 | subroutine validate_orbit_model_config | |
| 476 | select case (orbit_model) | ||
| 477 | case (ORBIT_GC) | ||
| 478 | ✗ | continue | |
| 479 | case (ORBIT_FULL_ORBIT) | ||
| 480 | ✗ | if (orbit_coord /= 1) error stop & | |
| 481 | ✗ | 'orbit_model=ORBIT_FULL_ORBIT supports only orbit_coord=1 (Boozer)' | |
| 482 | ✗ | if (class_plot .or. fast_class) error stop & | |
| 483 | ✗ | 'orbit_model=ORBIT_FULL_ORBIT does not support orbit classification' | |
| 484 | case default | ||
| 485 | error stop 'unsupported orbit_model (use 0 = guiding-center or '// & | ||
| 486 |
1/3✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 70 times.
|
70 | '7 = full orbit)' |
| 487 | end select | ||
| 488 | 52 | end subroutine validate_orbit_model_config | |
| 489 | |||
| 490 | !> Reject unsupported gpu_mode='production' runs before any field setup. | ||
| 491 | !> The native CUDA kernel reads compact RK tables that libneo fills only | ||
| 492 | !> from a Boozer chartmap carrying them (#569). | ||
| 493 | 70 | subroutine validate_gpu_production_config | |
| 494 | use field_boozer_chartmap, only: has_boozer_rk_tables | ||
| 495 | use magfie_sub, only: BOOZER | ||
| 496 | use params, only: field_input | ||
| 497 | |||
| 498 |
2/2✓ Branch 0 taken 69 times.
✓ Branch 1 taken 1 time.
|
70 | if (trim(gpu_mode) /= 'production') return |
| 499 | if (isw_field_type /= BOOZER .or. swcoll .or. & | ||
| 500 | len_trim(wall_input) > 0 .or. class_plot .or. fast_class .or. & | ||
| 501 |
8/16✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 time.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 time.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 time.
|
1 | ntcut > 0 .or. output_orbits_macrostep .or. orbit_model /= ORBIT_GC) & |
| 502 | error stop 'gpu_mode=production requires Boozer GC loss tracing '// & | ||
| 503 | ✗ | 'without collisions, wall, orbit output, or classifiers' | |
| 504 |
1/2✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
|
1 | if (.not. has_boozer_rk_tables(field_input)) then |
| 505 | 1 | print *, 'gpu_mode=production: field_input = ', trim(field_input) | |
| 506 | error stop 'gpu_mode=production needs field_input to be a Boozer '// & | ||
| 507 | 'chartmap with RK tables (libneo booz_xform_to_boozer_chartmap); '// & | ||
| 508 | 1 | 'VMEC, coil, SPECTRE and chartmap-coordinate fields are unsupported' | |
| 509 | end if | ||
| 510 | 69 | end subroutine validate_gpu_production_config | |
| 511 | |||
| 512 | 92 | subroutine init_field(self, vmec_file, ans_s, ans_tp, amultharm, aintegmode) | |
| 513 | 70 | use spectre_fo_field, only: spectre_fo_field_clear | |
| 514 | use field_base, only: magnetic_field_t | ||
| 515 | use field, only: field_from_file, is_spectre_file | ||
| 516 | use field_boozer_chartmap, only: boozer_chartmap_field_t, is_boozer_chartmap | ||
| 517 | use boozer_chartmap, only: set_boozer_chartmap_spline_orders, & | ||
| 518 | load_boozer_from_chartmap | ||
| 519 | use timing, only: print_phase_time | ||
| 520 | use magfie_sub, only: TEST, CANFLUX, VMEC, BOOZER, MEISS, ALBERT, & | ||
| 521 | REFCOORDS, SPECTRE, set_magfie_refcoords_field | ||
| 522 | #ifdef SIMPLE_ENABLE_VMECPP | ||
| 523 | use magfie_sub, only: VMECPP | ||
| 524 | #endif | ||
| 525 | use field_splined, only: splined_field_t, create_splined_field | ||
| 526 | use field_vmec, only: vmec_field_t | ||
| 527 | use reference_coordinates, only: init_reference_coordinates, ref_coords | ||
| 528 | use params, only: coord_input, field_input, wall_input, wall_units | ||
| 529 | use field_can_mod, only: init_field_can_boozer_chartmap_compact | ||
| 530 | use new_vmec_stuff_mod, only: nper | ||
| 531 | use orbit_symplectic_base, only: sympl_rmax | ||
| 532 | |||
| 533 | character(*), intent(in) :: vmec_file | ||
| 534 | type(tracer_t), intent(inout) :: self | ||
| 535 | integer, intent(in) :: ans_s, ans_tp, amultharm, aintegmode | ||
| 536 | 182 | class(magnetic_field_t), allocatable :: field_temp | |
| 537 | 92 | character(:), allocatable :: vmec_equilibrium_file | |
| 538 | logical :: use_boozer_chartmap, use_spectre, use_vmecpp | ||
| 539 | |||
| 540 | ! Release the previous FO tables when a library run changes fields. | ||
| 541 | 92 | call spectre_fo_field_clear | |
| 542 | 92 | self%integmode = aintegmode | |
| 543 | ! Reset the default domain when a library run changes field backend. | ||
| 544 | 92 | sympl_rmax = 1.0_dp | |
| 545 | !$acc update device(sympl_rmax) | ||
| 546 | 92 | map_boundary_exit_code = ORBIT_EXIT_LCFS | |
| 547 | 92 | chart_boundary_kind_effective = 'lcfs' | |
| 548 | 92 | invariant_edge_radius = 1.0_dp | |
| 549 | |||
| 550 | ! Check if field_input is a Boozer chartmap or SPECTRE file (no VMEC needed) | ||
| 551 | 92 | use_boozer_chartmap = .false. | |
| 552 | 92 | use_spectre = .false. | |
| 553 | 92 | use_vmecpp = is_vmecpp_input(vmec_file) | |
| 554 |
3/4✓ Branch 0 taken 89 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 89 times.
✗ Branch 3 not taken.
|
92 | if (len_trim(field_input) > 0 .and. .not. use_vmecpp) then |
| 555 | 89 | use_spectre = is_spectre_file(field_input) | |
| 556 |
2/2✓ Branch 0 taken 63 times.
✓ Branch 1 taken 26 times.
|
89 | if (.not. use_spectre) use_boozer_chartmap = is_boozer_chartmap(field_input) |
| 557 | end if | ||
| 558 | |||
| 559 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 92 times.
|
92 | if (gpu_compact_init) then |
| 560 | if (trim(gpu_backend) /= 'cuda-native' .or. & | ||
| 561 | trim(gpu_start_coordinates) /= 'boozer' .or. & | ||
| 562 | ✗ | .not. use_boozer_chartmap .or. len_trim(wall_input) > 0) & | |
| 563 | ✗ | error stop 'compact init requires native CUDA Boozer starts without wall' | |
| 564 | end if | ||
| 565 | |||
| 566 | ! TEST field is analytic - no VMEC or field files needed | ||
| 567 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 88 times.
|
92 | if (isw_field_type == TEST) then |
| 568 | 4 | invariant_edge_radius = 0.5_dp | |
| 569 | 4 | self%fper = twopi ! Full torus for analytic tokamak | |
| 570 | 4 | call print_phase_time('TEST field mode - no input files required') | |
| 571 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 88 times.
|
88 | else if (use_vmecpp) then |
| 572 | #ifdef SIMPLE_ENABLE_VMECPP | ||
| 573 | block | ||
| 574 | use field_vmecpp, only: get_vmecpp_metadata, open_vmecpp_geometry | ||
| 575 | use new_vmec_stuff_mod, only: rmajor | ||
| 576 | integer :: nper_vmecpp | ||
| 577 | call open_vmecpp_geometry(vmec_file) | ||
| 578 | call get_vmecpp_metadata(nper_vmecpp, rmajor) | ||
| 579 | nper = nper_vmecpp | ||
| 580 | self%fper = twopi/real(nper, dp) | ||
| 581 | isw_field_type = VMECPP | ||
| 582 | end block | ||
| 583 | call print_phase_time('VMEC++ direct geometry initialization completed') | ||
| 584 | #else | ||
| 585 | ✗ | error stop 'SIMPLE was built without SIMPLE_ENABLE_VMECPP' | |
| 586 | #endif | ||
| 587 |
2/2✓ Branch 0 taken 26 times.
✓ Branch 1 taken 62 times.
|
88 | else if (use_spectre) then |
| 588 | 26 | call init_spectre_field(self) | |
| 589 | 24 | call print_phase_time('SPECTRE field loading completed') | |
| 590 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 56 times.
|
62 | else if (use_boozer_chartmap) then |
| 591 | 6 | call set_boozer_chartmap_spline_orders(ans_s, ans_tp) | |
| 592 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (gpu_compact_init) then |
| 593 | ✗ | call load_boozer_from_chartmap(field_input, .true.) | |
| 594 | ✗ | call init_field_can_boozer_chartmap_compact | |
| 595 | ✗ | self%fper = twopi/real(nper, dp) | |
| 596 | ✗ | call print_phase_time('Compact Boozer RK table loading completed') | |
| 597 | else | ||
| 598 | ! Boozer chartmap: file-based, no VMEC initialization needed | ||
| 599 | 6 | call init_reference_coordinates(coord_input) | |
| 600 | call print_phase_time('Reference coordinate system '// & | ||
| 601 | 6 | 'initialization completed') | |
| 602 | |||
| 603 | block | ||
| 604 | use libneo_coordinates, only: chartmap_coordinate_system_t | ||
| 605 | select type (cs => ref_coords) | ||
| 606 | class is (chartmap_coordinate_system_t) | ||
| 607 | 6 | self%fper = twopi / real(cs%num_field_periods, dp) | |
| 608 | class default | ||
| 609 | ✗ | self%fper = twopi | |
| 610 | end select | ||
| 611 | end block | ||
| 612 | |||
| 613 | 6 | call init_stl_wall_if_enabled(coord_input) | |
| 614 | 6 | call print_phase_time('STL wall initialization completed') | |
| 615 | |||
| 616 |
1/2✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
|
6 | if (self%integmode >= 0) then |
| 617 | call field_from_file(field_input, field_temp) | ||
| 618 | 6 | call print_phase_time('Boozer chartmap field loading completed') | |
| 619 | end if | ||
| 620 | end if | ||
| 621 | else | ||
| 622 | vmec_equilibrium_file = select_vmec_equilibrium_file(vmec_file, & | ||
| 623 | field_input, & | ||
| 624 |
1/2✓ Branch 1 taken 56 times.
✗ Branch 2 not taken.
|
56 | coord_input) |
| 625 | call init_vmec(vmec_equilibrium_file, ans_s, ans_tp, amultharm, & | ||
| 626 | 56 | self%fper) | |
| 627 | 56 | call print_phase_time('VMEC initialization completed') | |
| 628 | |||
| 629 | 56 | call init_reference_coordinates(coord_input) | |
| 630 | call print_phase_time('Reference coordinate system '// & | ||
| 631 | 56 | 'initialization completed') | |
| 632 | |||
| 633 | 56 | call init_stl_wall_if_enabled(coord_input) | |
| 634 | 56 | call print_phase_time('STL wall initialization completed') | |
| 635 | |||
| 636 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
|
112 | if (self%integmode >= 0) then |
| 637 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
|
56 | if (trim(field_input) == '') then |
| 638 | ✗ | print *, 'simple_main.init_field: field_input must be set (see ', & | |
| 639 | ✗ | 'params.apply_config_aliases)' | |
| 640 | ✗ | error stop | |
| 641 | end if | ||
| 642 | |||
| 643 | call field_from_file(field_input, field_temp) | ||
| 644 | 56 | call print_phase_time('Field from file loading completed') | |
| 645 | |||
| 646 |
2/2✓ Branch 0 taken 54 times.
✓ Branch 1 taken 2 times.
|
56 | if (isw_field_type == REFCOORDS) then |
| 647 | select type (field_temp) | ||
| 648 | type is (splined_field_t) | ||
| 649 | ✗ | call set_magfie_refcoords_field(field_temp) | |
| 650 | type is (vmec_field_t) | ||
| 651 | 4 | block | |
| 652 |
14/14✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 6 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 6 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 6 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 6 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 6 times.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 6 times.
✓ Branch 13 taken 2 times.
|
44 | type(splined_field_t) :: splined_vmec |
| 653 | call create_splined_field(field_temp, ref_coords, & | ||
| 654 | 2 | splined_vmec) | |
| 655 | 2 | call set_magfie_refcoords_field(splined_vmec) | |
| 656 | end block | ||
| 657 | class default | ||
| 658 | print *, & | ||
| 659 | 'simple_main.init_field: REFCOORDS requires '// & | ||
| 660 | ✗ | 'a splined field' | |
| 661 | print *, & | ||
| 662 | ✗ | 'Supported inputs: coils (auto-splined) or VMEC wout', & | |
| 663 | ✗ | ' (splined onto coord_input)' | |
| 664 | ✗ | error stop | |
| 665 | end select | ||
| 666 | end if | ||
| 667 | end if | ||
| 668 | end if | ||
| 669 | |||
| 670 |
2/2✓ Branch 0 taken 60 times.
✓ Branch 1 taken 30 times.
|
90 | if (self%integmode > 0) then |
| 671 | ✗ | select case (isw_field_type) | |
| 672 | case (VMEC) | ||
| 673 | error stop & | ||
| 674 | 'Symplectic guiding-center integrators require '// & | ||
| 675 | 'canonical field representation (set isw_field_type to TEST, '// & | ||
| 676 | ✗ | 'CANFLUX, BOOZER, MEISS, ALBERT, or SPECTRE)' | |
| 677 | case (TEST, CANFLUX, BOOZER, MEISS, ALBERT, SPECTRE) | ||
| 678 | ✗ | continue | |
| 679 | case default | ||
| 680 | error stop & | ||
| 681 | 'Unknown canonical field type for symplectic guiding-center '// & | ||
| 682 |
1/3✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
|
60 | 'integrator' |
| 683 | end select | ||
| 684 | end if | ||
| 685 | |||
| 686 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 86 times.
|
90 | if (isw_field_type == TEST) then |
| 687 | ! TEST field is fully analytic - no field file needed | ||
| 688 | call init_field_can(isw_field_type, n_r=canonical_grid_nr, & | ||
| 689 | n_th=canonical_grid_ntheta, n_phi=canonical_grid_nphi, & | ||
| 690 | 4 | transformation_relerr=canonical_ode_relerr) | |
| 691 | 4 | call print_phase_time('Canonical field initialization completed') | |
| 692 | else if (isw_field_type == CANFLUX .or. isw_field_type == BOOZER .or. & | ||
| 693 |
6/6✓ Branch 0 taken 40 times.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 38 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 35 times.
|
86 | isw_field_type == MEISS .or. isw_field_type == ALBERT) then |
| 694 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 51 times.
|
51 | if (gpu_compact_init) return |
| 695 | call init_field_can(isw_field_type, field_temp, canonical_grid_nr, & | ||
| 696 | canonical_grid_ntheta, canonical_grid_nphi, & | ||
| 697 | 51 | canonical_ode_relerr) | |
| 698 |
2/2✓ Branch 0 taken 46 times.
✓ Branch 1 taken 5 times.
|
51 | if (isw_field_type == BOOZER .or. & |
| 699 | isw_field_type == CANFLUX) then | ||
| 700 | 46 | call set_magfie_refcoords_field(field_temp) | |
| 701 | end if | ||
| 702 | 51 | call print_phase_time('Canonical field initialization completed') | |
| 703 | end if | ||
| 704 |
5/12✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 56 times.
✓ Branch 3 taken 34 times.
✓ Branch 4 taken 62 times.
✓ Branch 5 taken 28 times.
✓ Branch 6 taken 62 times.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
90 | end subroutine init_field |
| 705 | |||
| 706 | 407 | logical function is_vmecpp_input(path) | |
| 707 | character(*), intent(in) :: path | ||
| 708 | integer :: length | ||
| 709 | |||
| 710 | 407 | length = len_trim(path) | |
| 711 | 407 | is_vmecpp_input = length >= 5 | |
| 712 |
2/2✓ Branch 0 taken 248 times.
✓ Branch 1 taken 159 times.
|
407 | if (is_vmecpp_input) is_vmecpp_input = path(length - 4:length) == '.json' |
| 713 | 92 | end function is_vmecpp_input | |
| 714 | |||
| 715 | 26 | subroutine init_spectre_field(self) | |
| 716 | !> VMEC-free SPECTRE setup. Mirrors the Boozer-chartmap precedent: no | ||
| 717 | !> init_vmec; the equilibrium globals nper and rmajor are taken from the | ||
| 718 | !> SPECTRE file so stevvo and params_init produce a consistent | ||
| 719 | !> dphi/dtaumin/fper. rmajor is the m=0,n=0 R harmonic of the outermost | ||
| 720 | !> interface (SI meters); stevvo scales it to cm. For integmode > 0 the | ||
| 721 | !> per-volume Meiss canonical coordinates are built here (#439). | ||
| 722 | use field_spectre, only: spectre_field_t, create_spectre_field | ||
| 723 | use field_can_spectre, only: set_spectre_construction_grid | ||
| 724 | use magfie_sub, only: set_magfie_spectre_field, SPECTRE | ||
| 725 | use new_vmec_stuff_mod, only: nper, rmajor | ||
| 726 | use util, only: twopi | ||
| 727 | use params, only: field_input, spectre_ncon_r, spectre_ncon_th, & | ||
| 728 | spectre_ncon_phi, spectre_ncon_order, & | ||
| 729 | spectre_ncon_ode_max_steps | ||
| 730 | use params, only: spectre_ncon_ode_relerr | ||
| 731 | use timing, only: print_phase_time | ||
| 732 | use orbit_symplectic_base, only: sympl_rmax | ||
| 733 | use spectre_fo_hybrid, only: spectre_fo_configure | ||
| 734 | use params, only: spectre_fo_tolerance, spectre_fo_grid | ||
| 735 | |||
| 736 | type(tracer_t), intent(inout) :: self | ||
| 737 | |||
| 738 | 26 | type(spectre_field_t) :: sf | |
| 739 | integer :: ierr, ii | ||
| 740 | real(dp) :: r00 | ||
| 741 | |||
| 742 | 26 | call create_spectre_field(sf, field_input, ierr) | |
| 743 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
|
26 | if (ierr /= 0) then |
| 744 | ✗ | print *, 'init_spectre_field: create_spectre_field failed for ', & | |
| 745 | ✗ | trim(field_input), ' (ierr = ', ierr, ')' | |
| 746 | ✗ | error stop | |
| 747 | end if | ||
| 748 |
6/6✓ Branch 0 taken 12 times.
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 6 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 4 times.
|
26 | if (sf%data%Mvol > 1 .and. integmode <= 0 .and. & |
| 749 | .not. generate_start_only) error stop & | ||
| 750 | 'multi-volume SPECTRE tracing requires integmode > 0 and crossing_level = 3; '// & | ||
| 751 | 2 | 'see DOC/spectre-rejected-interface-models.md' | |
| 752 | 24 | call set_magfie_spectre_field(sf) | |
| 753 | |||
| 754 | ! Shared by simple.x and simple_api, before any parallel tracing. | ||
| 755 | 24 | call spectre_fo_configure(spectre_fo_tolerance, spectre_fo_grid) | |
| 756 | |||
| 757 | 24 | r00 = 0.0d0 | |
| 758 |
2/2✓ Branch 0 taken 318 times.
✓ Branch 1 taken 24 times.
|
342 | do ii = 1, sf%data%mn |
| 759 |
2/2✓ Branch 0 taken 46 times.
✓ Branch 1 taken 272 times.
|
342 | if (sf%data%im(ii) == 0) then |
| 760 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 22 times.
|
46 | if (sf%data%in(ii) == 0) then |
| 761 | 24 | r00 = sf%data%Rbc(ii, sf%data%Mvol) | |
| 762 | end if | ||
| 763 | end if | ||
| 764 | end do | ||
| 765 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
|
24 | if (r00 <= 0.0d0) error stop & |
| 766 | ✗ | 'init_spectre_field: no positive m=0,n=0 R harmonic on outer interface' | |
| 767 | |||
| 768 | 24 | nper = sf%data%Nfp | |
| 769 | 24 | rmajor = r00 | |
| 770 | 24 | self%fper = twopi/real(sf%data%Nfp, dp) | |
| 771 | |||
| 772 | 24 | print *, 'SPECTRE field: Nfp = ', sf%data%Nfp, ' Mvol = ', sf%data%Mvol, & | |
| 773 | 48 | ' rmajor = ', rmajor, ' m' | |
| 774 | |||
| 775 |
2/2✓ Branch 0 taken 17 times.
✓ Branch 1 taken 7 times.
|
24 | if (self%integmode > 0) then |
| 776 | ! Allow Newton trials beyond every volume, including the outer edge. | ||
| 777 | ! Physical loss is localized at rho_g = Mvol by the event handler. | ||
| 778 | 17 | sympl_rmax = real(sf%data%Mvol + 1, dp) | |
| 779 | !$acc update device(sympl_rmax) | ||
| 780 | call set_spectre_construction_grid(spectre_ncon_r, spectre_ncon_th, & | ||
| 781 | spectre_ncon_phi, spectre_ncon_order, & | ||
| 782 | spectre_ncon_ode_max_steps, & | ||
| 783 | 17 | spectre_ncon_ode_relerr) | |
| 784 | 17 | call init_field_can(SPECTRE, sf) | |
| 785 | 17 | call print_phase_time('SPECTRE per-volume canonical construction completed') | |
| 786 | end if | ||
| 787 |
21/40✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 24 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 24 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 24 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 24 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 24 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 24 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 24 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 24 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 24 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 24 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 24 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 24 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 24 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 40 times.
✓ Branch 32 taken 24 times.
✓ Branch 33 taken 40 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 40 times.
✗ Branch 36 not taken.
✓ Branch 37 taken 40 times.
✗ Branch 38 not taken.
✓ Branch 39 taken 40 times.
✗ Branch 40 not taken.
|
64 | end subroutine init_spectre_field |
| 788 | |||
| 789 | 62 | subroutine init_stl_wall_if_enabled(coord_file) | |
| 790 | character(len=*), intent(in) :: coord_file | ||
| 791 | |||
| 792 | type(chartmap_metadata_t) :: meta | ||
| 793 | character(len=16) :: units | ||
| 794 | |||
| 795 | 62 | wall_enabled = .false. | |
| 796 | 62 | map_boundary_exit_code = ORBIT_EXIT_LCFS | |
| 797 | 62 | chart_boundary_kind_effective = 'lcfs' | |
| 798 | 62 | chartmap_cart_scale_to_m = -1.0d0 | |
| 799 | 62 | wall_query_rho_lcfs = -1.0d0 | |
| 800 | |||
| 801 | 62 | if (.not. allocated(ref_coords)) then | |
| 802 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
62 | if (len_trim(wall_input) == 0) return |
| 803 | ✗ | error stop "wall_input set but reference coordinates not initialized" | |
| 804 | end if | ||
| 805 | |||
| 806 | select type (ref_coords) | ||
| 807 | class is (chartmap_coordinate_system_t) | ||
| 808 | 8 | call read_chartmap_metadata(coord_file, meta) | |
| 809 | 8 | chartmap_cart_scale_to_m = meta%cart_scale_to_m | |
| 810 | 16 | select case (trim(chart_boundary_kind)) | |
| 811 | case ('auto') | ||
| 812 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
|
8 | if (abs(meta%rho_lcfs - 1d0) <= 1d-12) then |
| 813 | 6 | map_boundary_exit_code = ORBIT_EXIT_LCFS | |
| 814 | 6 | chart_boundary_kind_effective = 'lcfs' | |
| 815 | else | ||
| 816 | 2 | map_boundary_exit_code = ORBIT_EXIT_NUMERICAL_DOMAIN | |
| 817 | 2 | chart_boundary_kind_effective = 'domain' | |
| 818 | end if | ||
| 819 | case ('lcfs') | ||
| 820 | ✗ | map_boundary_exit_code = ORBIT_EXIT_LCFS | |
| 821 | ✗ | chart_boundary_kind_effective = 'lcfs' | |
| 822 | case ('wall') | ||
| 823 | ✗ | map_boundary_exit_code = ORBIT_EXIT_WALL | |
| 824 | ✗ | chart_boundary_kind_effective = 'wall' | |
| 825 | case ('domain') | ||
| 826 | ✗ | map_boundary_exit_code = ORBIT_EXIT_NUMERICAL_DOMAIN | |
| 827 |
1/5✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
8 | chart_boundary_kind_effective = 'domain' |
| 828 | end select | ||
| 829 |
3/4✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
|
8 | if (meta%rho_lcfs > 0.0_dp .and. meta%rho_lcfs <= 1.0_dp) then |
| 830 | 6 | invariant_edge_radius = meta%rho_lcfs**2 | |
| 831 | end if | ||
| 832 | class default | ||
| 833 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
|
54 | if (len_trim(wall_input) == 0) return |
| 834 | ✗ | error stop "wall_input requires chartmap reference coordinates" | |
| 835 | end select | ||
| 836 | |||
| 837 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | if (len_trim(wall_input) == 0) return |
| 838 | |||
| 839 | ! The STL is an external physical wall enclosing the LCFS. Querying a | ||
| 840 | ! straight Cartesian chord while the accepted endpoint is still inside | ||
| 841 | ! the LCFS can create a false intersection in a non-convex toroidal | ||
| 842 | ! geometry. Preserve the historical SIMPLE rule: activate wall queries | ||
| 843 | ! only once an accepted microstep endpoint has passed the recorded LCFS. | ||
| 844 | ✗ | wall_query_rho_lcfs = meta%rho_lcfs | |
| 845 | |||
| 846 | ✗ | units = adjustl(wall_units) | |
| 847 | ✗ | select case (trim(units)) | |
| 848 | case ("m", "M") | ||
| 849 | ✗ | call stl_wall_init(wall, trim(wall_input), 1.0d0) | |
| 850 | case ("mm", "MM") | ||
| 851 | ✗ | call stl_wall_init(wall, trim(wall_input), 1.0d-3) | |
| 852 | case default | ||
| 853 | ✗ | print *, "init_stl_wall_if_enabled: invalid wall_units=", trim(wall_units) | |
| 854 | ✗ | error stop "wall_units must be m or mm" | |
| 855 | end select | ||
| 856 | |||
| 857 | ✗ | wall_enabled = .true. | |
| 858 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 62 times.
|
88 | end subroutine init_stl_wall_if_enabled |
| 859 | |||
| 860 | 56 | function select_vmec_equilibrium_file(vmec_file_default, field_file, coord_file) & | |
| 861 | 56 | result(vmec_file) | |
| 862 | use libneo_coordinates, only: detect_refcoords_file_type, & | ||
| 863 | refcoords_file_vmec_wout | ||
| 864 | character(*), intent(in) :: vmec_file_default | ||
| 865 | character(*), intent(in) :: field_file | ||
| 866 | character(*), intent(in) :: coord_file | ||
| 867 | character(:), allocatable :: vmec_file | ||
| 868 | |||
| 869 | integer :: file_type, ierr | ||
| 870 | character(len=2048) :: message | ||
| 871 | |||
| 872 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 8 times.
|
56 | vmec_file = trim(vmec_file_default) |
| 873 | |||
| 874 |
1/2✓ Branch 0 taken 56 times.
✗ Branch 1 not taken.
|
56 | if (len_trim(field_file) > 0) then |
| 875 |
1/2✓ Branch 0 taken 56 times.
✗ Branch 1 not taken.
|
56 | if (ends_with_nc(field_file)) then |
| 876 | call detect_refcoords_file_type(trim(field_file), file_type, & | ||
| 877 |
1/2✓ Branch 2 taken 56 times.
✗ Branch 3 not taken.
|
56 | ierr, message) |
| 878 |
2/4✓ Branch 0 taken 56 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 56 times.
✗ Branch 3 not taken.
|
56 | if (ierr == 0 .and. file_type == refcoords_file_vmec_wout) then |
| 879 |
4/6✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 48 times.
✓ Branch 4 taken 56 times.
✗ Branch 5 not taken.
|
56 | vmec_file = trim(field_file) |
| 880 | 56 | return | |
| 881 | end if | ||
| 882 | end if | ||
| 883 | end if | ||
| 884 | |||
| 885 | ✗ | if (len_trim(coord_file) > 0) then | |
| 886 | ✗ | if (ends_with_nc(coord_file)) then | |
| 887 | call detect_refcoords_file_type(trim(coord_file), file_type, & | ||
| 888 | ✗ | ierr, message) | |
| 889 | ✗ | if (ierr == 0 .and. file_type == refcoords_file_vmec_wout) then | |
| 890 | ✗ | vmec_file = trim(coord_file) | |
| 891 | ✗ | return | |
| 892 | end if | ||
| 893 | end if | ||
| 894 | end if | ||
| 895 | 56 | end function select_vmec_equilibrium_file | |
| 896 | |||
| 897 | 56 | function ends_with_nc(filename) result(is_nc) | |
| 898 | character(*), intent(in) :: filename | ||
| 899 | logical :: is_nc | ||
| 900 | integer :: n | ||
| 901 | |||
| 902 | 56 | n = len_trim(filename) | |
| 903 | 56 | is_nc = .false. | |
| 904 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
|
56 | if (n >= 3) then |
| 905 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
|
56 | if (filename(n - 2:n) == '.nc') then |
| 906 | 56 | is_nc = .true. | |
| 907 | end if | ||
| 908 | end if | ||
| 909 | 56 | end function ends_with_nc | |
| 910 | |||
| 911 | 49 | subroutine trace_parallel(norb) | |
| 912 | use netcdf_orbit_output, only: init_orbit_netcdf, close_orbit_netcdf, & | ||
| 913 | write_orbit_to_netcdf | ||
| 914 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 915 | use params, only: debug | ||
| 916 | #endif | ||
| 917 | |||
| 918 | type(tracer_t), intent(inout) :: norb | ||
| 919 | integer :: i | ||
| 920 | 49 | real(dp), allocatable :: traj(:, :), times(:) | |
| 921 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 922 | integer :: debug_thread | ||
| 923 | #endif | ||
| 924 | |||
| 925 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 42 times.
|
49 | if (output_orbits_macrostep) then |
| 926 | 7 | call init_orbit_netcdf(ntestpart, ntimstep) | |
| 927 | end if | ||
| 928 | |||
| 929 | 49 | !$omp parallel firstprivate(norb) private(traj, times, i) | |
| 930 | allocate (traj(5, ntimstep), times(ntimstep)) | ||
| 931 | |||
| 932 | ! Independent markers have highly nonuniform trace costs. One-marker | ||
| 933 | ! dynamic chunks keep unfinished work available to every thread instead | ||
| 934 | ! of queueing it behind a pathological marker in a static block. | ||
| 935 | !$omp do schedule(dynamic, 1) | ||
| 936 | do i = 1, ntestpart | ||
| 937 | if (allocated(particle_done)) then | ||
| 938 | if (particle_done(i)) then | ||
| 939 | call progress_tick | ||
| 940 | cycle | ||
| 941 | end if | ||
| 942 | end if | ||
| 943 | |||
| 944 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 945 | if (debug) then | ||
| 946 | !$omp critical | ||
| 947 | kpart = kpart + 1 | ||
| 948 | debug_thread = 0 | ||
| 949 | !$ debug_thread = omp_get_thread_num() | ||
| 950 | print *, kpart, ' / ', ntestpart, 'particle: ', i, 'thread: ', & | ||
| 951 | debug_thread | ||
| 952 | !$omp end critical | ||
| 953 | end if | ||
| 954 | #endif | ||
| 955 | |||
| 956 | call trace_orbit(norb, i, traj, times) | ||
| 957 | |||
| 958 | if (output_orbits_macrostep) then | ||
| 959 | !$omp critical | ||
| 960 | call write_orbit_to_netcdf(i, traj, times) | ||
| 961 | !$omp end critical | ||
| 962 | end if | ||
| 963 | |||
| 964 | call progress_tick | ||
| 965 | end do | ||
| 966 | !$omp end do | ||
| 967 | !$omp end parallel | ||
| 968 | |||
| 969 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 42 times.
|
49 | if (output_orbits_macrostep) then |
| 970 | 7 | call close_orbit_netcdf() | |
| 971 | end if | ||
| 972 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 49 times.
|
49 | end subroutine trace_parallel |
| 973 | |||
| 974 | ✗ | subroutine trace_gpu_production(norb) | |
| 975 | ! End-to-end loss-tracing path used by the Landreman comparison. This | ||
| 976 | ! related optimization and loss-threshold protocol is described at | ||
| 977 | ! https://arxiv.org/abs/2606.19523. | ||
| 978 | ! path deliberately excludes the separate continuous fast-classifier | ||
| 979 | ! campaign (J_parallel and rotation only) and all unsupported physics. | ||
| 980 | 49 | use finite_check, only: quiet_nan | |
| 981 | use simple, only: init_sympl | ||
| 982 | use simple_gpu, only: trace_orbits_gpu_landreman, trace_orbits_gpu_method | ||
| 983 | use field_can_boozer, only: eval_field_booz | ||
| 984 | use field_can_mod, only: get_val | ||
| 985 | use boozer_sub, only: boozer_rk_device_supported, sync_boozer_state | ||
| 986 | use magfie_sub, only: BOOZER | ||
| 987 | use orbit_symplectic_base, only: EXPL_IMPL_EULER, MIDPOINT, & | ||
| 988 | CASH_KARP, DORMAND_PRINCE | ||
| 989 | #ifdef SIMPLE_ENABLE_CUDA_NATIVE | ||
| 990 | use simple_cuda_native, only: trace_orbits_cuda_native_landreman, & | ||
| 991 | initialize_cuda_native_particle, evaluate_cuda_native_particle | ||
| 992 | #endif | ||
| 993 | |||
| 994 | type(tracer_t), intent(inout) :: norb | ||
| 995 | ✗ | type(symplectic_integrator_t), allocatable :: si_gpu(:) | |
| 996 | ✗ | type(field_can_t), allocatable :: f_gpu(:) | |
| 997 | ✗ | integer, allocatable :: loss_step(:), nfev(:) | |
| 998 | ✗ | real(dp), allocatable :: loss_time_normalized(:), zcanonical(:, :) | |
| 999 | ✗ | logical, allocatable :: passing(:) | |
| 1000 | character(32) :: backend_name | ||
| 1001 | character(32) :: start_coordinates | ||
| 1002 | integer :: method, init_mode, i, it | ||
| 1003 | integer :: particle_profile_unit | ||
| 1004 | integer(8) :: nfev_total, warp_nfev_slots | ||
| 1005 | real(dp) :: z(5), total_time_normalized, t0, t_init, t_trace, t_finish | ||
| 1006 | real(dp) :: block_time, loss_tau, maxloss, minimum_timestep | ||
| 1007 | real(dp) :: time_scale, observed_duration, energy_loss_fraction | ||
| 1008 | real(dp) :: cuda_profile_ms(5) | ||
| 1009 | logical :: landreman_mode, cuda_native_backend | ||
| 1010 | |||
| 1011 | ✗ | call validate_gpu_production_config | |
| 1012 | |||
| 1013 | ✗ | backend_name = trim(gpu_backend) | |
| 1014 | ✗ | select case (trim(backend_name)) | |
| 1015 | case ('auto') | ||
| 1016 | ✗ | backend_name = 'cuda-native' | |
| 1017 | case ('openacc') | ||
| 1018 | ✗ | error stop 'gpu_backend=openacc is not supported; build with SIMPLE_ENABLE_CUDA' | |
| 1019 | case ('cuda-native') | ||
| 1020 | #ifndef SIMPLE_ENABLE_CUDA_NATIVE | ||
| 1021 | ✗ | error stop 'gpu_backend=cuda-native requires SIMPLE_ENABLE_CUDA' | |
| 1022 | #endif | ||
| 1023 | case default | ||
| 1024 | ✗ | error stop 'gpu_backend must be openacc or cuda-native' | |
| 1025 | end select | ||
| 1026 | ✗ | select case (trim(backend_name)) | |
| 1027 | case ('openacc') | ||
| 1028 | ✗ | cuda_native_backend = .false. | |
| 1029 | case ('cuda-native') | ||
| 1030 | #ifdef SIMPLE_ENABLE_CUDA_NATIVE | ||
| 1031 | cuda_native_backend = .true. | ||
| 1032 | #else | ||
| 1033 | ✗ | error stop 'gpu_backend=cuda-native requires SIMPLE_ENABLE_CUDA' | |
| 1034 | #endif | ||
| 1035 | end select | ||
| 1036 | |||
| 1037 | ✗ | select case (trim(gpu_method)) | |
| 1038 | case ('auto') | ||
| 1039 | ✗ | if (cuda_native_backend) then | |
| 1040 | method = DORMAND_PRINCE | ||
| 1041 | else | ||
| 1042 | ✗ | method = EXPL_IMPL_EULER | |
| 1043 | end if | ||
| 1044 | case default | ||
| 1045 | ✗ | select case (trim(gpu_method)) | |
| 1046 | case ('euler') | ||
| 1047 | ✗ | method = EXPL_IMPL_EULER | |
| 1048 | case ('midpoint') | ||
| 1049 | ✗ | method = MIDPOINT | |
| 1050 | case ('cash-karp', 'cash_karp', 'ck') | ||
| 1051 | ✗ | method = CASH_KARP | |
| 1052 | case ('dormand-prince', 'dormand_prince', 'dopri', 'rk45') | ||
| 1053 | ✗ | method = DORMAND_PRINCE | |
| 1054 | case default | ||
| 1055 | ✗ | error stop 'gpu_method must be auto, euler, midpoint, cash-karp, or dopri' | |
| 1056 | end select | ||
| 1057 | end select | ||
| 1058 | if (method /= EXPL_IMPL_EULER .and. method /= MIDPOINT .and. & | ||
| 1059 | ✗ | method /= CASH_KARP .and. method /= DORMAND_PRINCE) & | |
| 1060 | ✗ | error stop 'gpu_method supports euler, midpoint, cash-karp, or dopri' | |
| 1061 | if (cuda_native_backend .and. method /= CASH_KARP .and. & | ||
| 1062 | method /= DORMAND_PRINCE) & | ||
| 1063 | error stop 'cuda-native backend supports only cash-karp or dopri' | ||
| 1064 | |||
| 1065 | ✗ | init_mode = method | |
| 1066 | ✗ | if (method == CASH_KARP .or. method == DORMAND_PRINCE) & | |
| 1067 | ✗ | init_mode = EXPL_IMPL_EULER | |
| 1068 | allocate (si_gpu(ntestpart), f_gpu(ntestpart), loss_step(ntestpart), & | ||
| 1069 | nfev(ntestpart), loss_time_normalized(ntestpart), & | ||
| 1070 | ✗ | zcanonical(4, ntestpart), passing(ntestpart)) | |
| 1071 | |||
| 1072 | ✗ | start_coordinates = trim(gpu_start_coordinates) | |
| 1073 | ✗ | if (trim(start_coordinates) /= 'reference' .and. & | |
| 1074 | trim(start_coordinates) /= 'boozer') & | ||
| 1075 | ✗ | error stop 'gpu_start_coordinates must be reference or boozer' | |
| 1076 | |||
| 1077 | ✗ | t0 = get_wtime() | |
| 1078 | ✗ | do i = 1, ntestpart | |
| 1079 | ✗ | if (trim(start_coordinates) == 'boozer') then | |
| 1080 | ✗ | z(1:3) = zstart(1:3, i) | |
| 1081 | else | ||
| 1082 | ✗ | call ref_to_integ(zstart(1:3, i), z(1:3)) | |
| 1083 | end if | ||
| 1084 | ✗ | z(4:5) = zstart(4:5, i) | |
| 1085 | ✗ | if (gpu_compact_init .and. cuda_native_backend) then | |
| 1086 | #ifdef SIMPLE_ENABLE_CUDA_NATIVE | ||
| 1087 | call initialize_cuda_native_particle(si_gpu(i), f_gpu(i), z, & | ||
| 1088 | dtaumin, relerr) | ||
| 1089 | passing(i) = .false. | ||
| 1090 | trap_par(i) = 0.0_dp | ||
| 1091 | perp_inv(i) = z(4)**2*(1.0_dp - z(5)**2)/f_gpu(i)%Bmod | ||
| 1092 | #endif | ||
| 1093 | else | ||
| 1094 | call init_sympl(si_gpu(i), f_gpu(i), z, dtaumin, dtaumin, & | ||
| 1095 | ✗ | relerr, init_mode) | |
| 1096 | call compute_pitch_angle_params(z, passing(i), trap_par(i), & | ||
| 1097 | ✗ | perp_inv(i)) | |
| 1098 | end if | ||
| 1099 | end do | ||
| 1100 | ✗ | call sync_boozer_state | |
| 1101 | ✗ | if ((method == CASH_KARP .or. method == DORMAND_PRINCE) .and. & | |
| 1102 | .not. boozer_rk_device_supported()) & | ||
| 1103 | ✗ | error stop 'GPU RK requires scalar cubic or quintic Boozer splines' | |
| 1104 | ✗ | t_init = get_wtime() - t0 | |
| 1105 | |||
| 1106 | ✗ | landreman_mode = gpu_landreman | |
| 1107 | if (cuda_native_backend .and. .not. landreman_mode) & | ||
| 1108 | error stop 'gpu_landreman must be true for the native CUDA backend' | ||
| 1109 | |||
| 1110 | ✗ | t0 = get_wtime() | |
| 1111 | ✗ | warp_nfev_slots = 0_8 | |
| 1112 | cuda_profile_ms = 0.0_dp | ||
| 1113 | ✗ | if (landreman_mode) then | |
| 1114 | ✗ | block_time = gpu_t_block | |
| 1115 | ✗ | loss_tau = gpu_loss_tau | |
| 1116 | ✗ | maxloss = gpu_maxloss | |
| 1117 | ✗ | minimum_timestep = gpu_min_timestep | |
| 1118 | ✗ | if (block_time <= 0.0_dp .or. block_time > trace_time) & | |
| 1119 | ✗ | error stop 'gpu_t_block must be in (0, trace_time]' | |
| 1120 | ✗ | if (loss_tau <= 0.0_dp) & | |
| 1121 | ✗ | error stop 'gpu_loss_tau must be positive' | |
| 1122 | ✗ | if (maxloss < 0.0_dp) & | |
| 1123 | ✗ | error stop 'gpu_maxloss must be nonnegative' | |
| 1124 | ✗ | if (minimum_timestep < 0.0_dp) & | |
| 1125 | ✗ | error stop 'gpu_min_timestep must be nonnegative' | |
| 1126 | ✗ | time_scale = sqrt2/v0 | |
| 1127 | total_time_normalized = & | ||
| 1128 | ✗ | real(sum(ntau_macro(2:ntimstep)), dp)*si_gpu(1)%dt | |
| 1129 | if (cuda_native_backend) then | ||
| 1130 | #ifdef SIMPLE_ENABLE_CUDA_NATIVE | ||
| 1131 | call trace_orbits_cuda_native_landreman(si_gpu, f_gpu, & | ||
| 1132 | ntestpart, method, total_time_normalized, & | ||
| 1133 | block_time/time_scale, time_scale, loss_tau, maxloss, & | ||
| 1134 | minimum_timestep/time_scale, loss_step, & | ||
| 1135 | loss_time_normalized, zcanonical, nfev, observed_duration, & | ||
| 1136 | energy_loss_fraction, warp_nfev_slots, cuda_profile_ms) | ||
| 1137 | #endif | ||
| 1138 | else | ||
| 1139 | call trace_orbits_gpu_landreman(si_gpu, f_gpu, ntestpart, & | ||
| 1140 | ntimstep, ntau_macro, method, block_time/time_scale, & | ||
| 1141 | time_scale, loss_tau, maxloss, & | ||
| 1142 | minimum_timestep/time_scale, loss_step, & | ||
| 1143 | loss_time_normalized, zcanonical, nfev, observed_duration, & | ||
| 1144 | ✗ | energy_loss_fraction, warp_nfev_slots) | |
| 1145 | end if | ||
| 1146 | else | ||
| 1147 | call trace_orbits_gpu_method(si_gpu, f_gpu, ntestpart, ntimstep, & | ||
| 1148 | ✗ | ntau_macro, method, loss_step, loss_time_normalized, zcanonical, nfev) | |
| 1149 | observed_duration = & | ||
| 1150 | ✗ | real(sum(ntau_macro(2:ntimstep)), dp)*si_gpu(1)%dt | |
| 1151 | ✗ | energy_loss_fraction = 0.0_dp | |
| 1152 | end if | ||
| 1153 | ✗ | t_trace = get_wtime() - t0 | |
| 1154 | |||
| 1155 | ✗ | total_time_normalized = real(sum(ntau_macro(2:ntimstep)), dp)*si_gpu(1)%dt | |
| 1156 | ✗ | confpart_pass = 0.0_dp | |
| 1157 | ✗ | confpart_trap = 0.0_dp | |
| 1158 | ✗ | unresolved_orbits = 0 | |
| 1159 | ✗ | nfev_total = 0_8 | |
| 1160 | ✗ | t0 = get_wtime() | |
| 1161 | ✗ | do i = 1, ntestpart | |
| 1162 | ✗ | nfev_total = nfev_total + int(nfev(i), 8) | |
| 1163 | ✗ | si_gpu(i)%z = zcanonical(:, i) | |
| 1164 | if (gpu_compact_init .and. cuda_native_backend) then | ||
| 1165 | #ifdef SIMPLE_ENABLE_CUDA_NATIVE | ||
| 1166 | call evaluate_cuda_native_particle(f_gpu(i), zcanonical(1, i), & | ||
| 1167 | zcanonical(2, i), zcanonical(3, i)) | ||
| 1168 | #endif | ||
| 1169 | else | ||
| 1170 | call eval_field_booz(f_gpu(i), zcanonical(1, i), & | ||
| 1171 | ✗ | zcanonical(2, i), zcanonical(3, i), 0) | |
| 1172 | end if | ||
| 1173 | ✗ | call get_val(f_gpu(i), zcanonical(4, i)) | |
| 1174 | ✗ | call integ_to_ref(zcanonical(1:3, i), zend(1:3, i)) | |
| 1175 | ✗ | zend(4, i) = si_gpu(i)%pabs | |
| 1176 | ✗ | zend(5, i) = f_gpu(i)%vpar/(si_gpu(i)%pabs*sqrt2) | |
| 1177 | |||
| 1178 | ✗ | if (loss_step(i) < 0) then | |
| 1179 | ✗ | times_lost(i) = quiet_nan() | |
| 1180 | ✗ | orbit_exit_code(i) = ORBIT_EXIT_NUMERICAL_EVENT | |
| 1181 | ✗ | else if (loss_time_normalized(i) < total_time_normalized - & | |
| 1182 | 16.0_dp*epsilon(1.0_dp)*max(total_time_normalized, 1.0_dp)) then | ||
| 1183 | ✗ | times_lost(i) = loss_time_normalized(i)*sqrt2/v0 | |
| 1184 | ✗ | orbit_exit_code(i) = ORBIT_EXIT_LCFS | |
| 1185 | else | ||
| 1186 | ✗ | times_lost(i) = trace_time | |
| 1187 | ✗ | orbit_exit_code(i) = ORBIT_EXIT_COMPLETED | |
| 1188 | end if | ||
| 1189 | ✗ | do it = 1, ntimstep | |
| 1190 | ✗ | if (orbit_exit_code(i) == ORBIT_EXIT_COMPLETED .or. & | |
| 1191 | times_lost(i) >= real(kt_macro(it), dp)*dtaumin/v0) & | ||
| 1192 | ✗ | call increase_confined_count(it, passing(i)) | |
| 1193 | end do | ||
| 1194 | end do | ||
| 1195 | ✗ | if (.not. landreman_mode) then | |
| 1196 | ✗ | do i = 1, ntestpart, 32 | |
| 1197 | warp_nfev_slots = warp_nfev_slots + 32_8*int( & | ||
| 1198 | ✗ | maxval(nfev(i:min(i + 31, ntestpart))), 8) | |
| 1199 | end do | ||
| 1200 | end if | ||
| 1201 | ✗ | t_finish = get_wtime() - t0 | |
| 1202 | |||
| 1203 | ✗ | print *, '==================== GPU production trace ===================' | |
| 1204 | ✗ | print '(a,i0)', ' integrator mode = ', method | |
| 1205 | ✗ | print '(a,a)', ' backend = ', trim(backend_name) | |
| 1206 | ✗ | print '(a,i0)', ' particles = ', ntestpart | |
| 1207 | ✗ | print '(a,a)', ' start coordinates = ', trim(start_coordinates) | |
| 1208 | ✗ | print '(a,f12.6,a)', ' initialization = ', t_init, ' s' | |
| 1209 | ✗ | print '(a,f12.6,a)', ' tracing = ', t_trace, ' s' | |
| 1210 | ✗ | print '(a,f12.6,a)', ' finalization = ', t_finish, ' s' | |
| 1211 | ✗ | print '(a,f12.6,a)', ' end-to-end = ', t_init + t_trace + t_finish, ' s' | |
| 1212 | ✗ | print '(a,i0)', ' field evaluations = ', nfev_total | |
| 1213 | ✗ | print '(a,i0)', ' warp RHS slots = ', warp_nfev_slots | |
| 1214 | ✗ | if (warp_nfev_slots > 0_8) print '(a,f8.4)', & | |
| 1215 | ✗ | ' warp RHS utilization = ', real(nfev_total, dp)/ & | |
| 1216 | ✗ | real(warp_nfev_slots, dp) | |
| 1217 | ✗ | print '(a,i0)', ' physical losses = ', count(orbit_exit_code == ORBIT_EXIT_LCFS) | |
| 1218 | ✗ | print '(a,i0)', ' numerical failures = ', & | |
| 1219 | ✗ | count(orbit_exit_code >= ORBIT_EXIT_NUMERICAL_DOMAIN) | |
| 1220 | ✗ | if (landreman_mode) then | |
| 1221 | ✗ | print '(a,f12.6,a)', ' observed until = ', & | |
| 1222 | ✗ | observed_duration*sqrt2/v0, ' s' | |
| 1223 | ✗ | print '(a,es12.4)', ' energy loss fraction = ', energy_loss_fraction | |
| 1224 | end if | ||
| 1225 | if (cuda_native_backend) then | ||
| 1226 | print '(a,f12.6,a)', ' CUDA allocation = ', cuda_profile_ms(1)/1.0e3_dp, ' s' | ||
| 1227 | print '(a,f12.6,a)', ' CUDA upload = ', cuda_profile_ms(2)/1.0e3_dp, ' s' | ||
| 1228 | print '(a,f12.6,a)', ' CUDA kernel = ', cuda_profile_ms(3)/1.0e3_dp, ' s' | ||
| 1229 | print '(a,f12.6,a)', ' CUDA download = ', cuda_profile_ms(4)/1.0e3_dp, ' s' | ||
| 1230 | print '(a,f12.6,a)', ' CUDA C-ABI total= ', cuda_profile_ms(5)/1.0e3_dp, ' s' | ||
| 1231 | end if | ||
| 1232 | ✗ | print *, '============================================================' | |
| 1233 | |||
| 1234 | ✗ | if (len_trim(gpu_particle_profile) > 0) then | |
| 1235 | open (newunit=particle_profile_unit, & | ||
| 1236 | file=trim(gpu_particle_profile), & | ||
| 1237 | ✗ | status='replace', action='write') | |
| 1238 | write (particle_profile_unit, '(a)') & | ||
| 1239 | 'particle,s_start,theta_start,zeta_start,speed_start,'// & | ||
| 1240 | 'pitch_start,s_end,theta_end,zeta_end,speed_end,pitch_end,'// & | ||
| 1241 | ✗ | 'nfev,loss_step,loss_time' | |
| 1242 | ✗ | do i = 1, ntestpart | |
| 1243 | write (particle_profile_unit, & | ||
| 1244 | '(i0,10(",",es24.16),2(",",i0),",",es24.16)') & | ||
| 1245 | ✗ | i, zstart(:, i), zend(:, i), nfev(i), loss_step(i), & | |
| 1246 | ✗ | loss_time_normalized(i)*time_scale | |
| 1247 | end do | ||
| 1248 | ✗ | close (particle_profile_unit) | |
| 1249 | end if | ||
| 1250 | ✗ | end subroutine trace_gpu_production | |
| 1251 | |||
| 1252 | 2 | subroutine trace_compare_gpu(norb) | |
| 1253 | ! Validate and benchmark the OpenACC GPU tracing kernel against the CPU | ||
| 1254 | ! symplectic integrator on identical per-particle initial states. | ||
| 1255 | ! Triggered by gpu_mode='benchmark' and restricted to the Boozer | ||
| 1256 | ! explicit-implicit Euler or midpoint path without wall, collision, or | ||
| 1257 | ! classifier options. | ||
| 1258 | ✗ | use orbit_symplectic, only: orbit_timestep_sympl | |
| 1259 | use orbit_symplectic_base, only: symplectic_integrator_t, & | ||
| 1260 | EXPL_IMPL_EULER, MIDPOINT | ||
| 1261 | use field_can_mod, only: field_can_t, get_derivatives | ||
| 1262 | use field_can_boozer, only: eval_field_booz | ||
| 1263 | use magfie_sub, only: BOOZER | ||
| 1264 | use boozer_sub, only: sync_boozer_state | ||
| 1265 | use simple_gpu, only: evaluate_rhs_gpu, trace_orbits_gpu_method | ||
| 1266 | |||
| 1267 | type(tracer_t), intent(inout) :: norb | ||
| 1268 | |||
| 1269 | 2 | type(symplectic_integrator_t), allocatable :: si_cpu(:), si_gpu(:) | |
| 1270 | 2 | type(field_can_t), allocatable :: f_cpu(:), f_gpu(:) | |
| 1271 | 2 | integer, allocatable :: cpu_loss(:), gpu_loss(:), gpu_nfev(:) | |
| 1272 | 2 | real(dp), allocatable :: gpu_loss_time(:) | |
| 1273 | 2 | real(dp), allocatable :: cpu_zend(:, :), gpu_zend(:, :) | |
| 1274 | 2 | real(dp), allocatable :: rhs_y(:, :), rhs_mu(:), rhs_ro0(:) | |
| 1275 | 2 | real(dp), allocatable :: rhs_host(:, :), rhs_device(:, :) | |
| 1276 | real(dp) :: z(5), hprime, rhs_scale, rhs_error | ||
| 1277 | integer :: i, it, ktau, ierr, loss_mismatch, rhs_location(2) | ||
| 1278 | integer :: cpu_lost, gpu_lost, flip | ||
| 1279 | real(dp) :: t0, t1, t_cpu, t_gpu, maxz | ||
| 1280 | |||
| 1281 | if (isw_field_type /= BOOZER .or. & | ||
| 1282 | (integmode /= EXPL_IMPL_EULER .and. integmode /= MIDPOINT) .or. swcoll .or. & | ||
| 1283 |
7/14✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
|
2 | len_trim(wall_input) > 0 .or. class_plot .or. fast_class .or. generate_start_only) then |
| 1284 | error stop "simple_main.trace_compare_gpu: gpu_mode='benchmark' requires Boozer, " // & | ||
| 1285 | ✗ | "EXPL_IMPL_EULER or MIDPOINT, and no wall/collision/classifier options" | |
| 1286 | end if | ||
| 1287 | |||
| 1288 | 2 | call sync_boozer_state | |
| 1289 | |||
| 1290 |
15/26✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✓ Branch 17 taken 16 times.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 2 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 2 times.
✓ Branch 29 taken 16 times.
✓ Branch 30 taken 2 times.
|
40 | allocate (si_cpu(ntestpart), si_gpu(ntestpart)) |
| 1291 |
23/34✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✓ Branch 17 taken 6 times.
✓ Branch 18 taken 2 times.
✓ Branch 19 taken 12 times.
✓ Branch 20 taken 2 times.
✓ Branch 21 taken 16 times.
✓ Branch 22 taken 2 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 2 times.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✓ Branch 28 taken 2 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 2 times.
✓ Branch 33 taken 6 times.
✓ Branch 34 taken 2 times.
✓ Branch 35 taken 12 times.
✓ Branch 36 taken 2 times.
✓ Branch 37 taken 16 times.
✓ Branch 38 taken 2 times.
|
76 | allocate (f_cpu(ntestpart), f_gpu(ntestpart)) |
| 1292 |
15/30✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 2 times.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 2 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 2 times.
✓ Branch 29 taken 2 times.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✓ Branch 32 taken 2 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 2 times.
|
10 | allocate (cpu_loss(ntestpart), gpu_loss(ntestpart), gpu_nfev(ntestpart)) |
| 1293 |
6/12✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
|
6 | allocate (gpu_loss_time(ntestpart)) |
| 1294 |
11/22✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 2 times.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 2 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
|
8 | allocate (cpu_zend(4, ntestpart), gpu_zend(4, ntestpart)) |
| 1295 |
12/24✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✓ Branch 22 taken 2 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 2 times.
|
8 | allocate (rhs_y(4, ntestpart), rhs_mu(ntestpart), rhs_ro0(ntestpart)) |
| 1296 |
8/16✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 2 times.
|
6 | allocate (rhs_host(4, ntestpart), rhs_device(4, ntestpart)) |
| 1297 | |||
| 1298 | ! Identical per-particle initialisation (host). init_sympl sets the | ||
| 1299 | ! orbit_timestep_sympl procedure pointer for the CPU reference. | ||
| 1300 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
|
18 | do i = 1, ntestpart |
| 1301 | 16 | call ref_to_integ(zstart(1:3, i), z(1:3)) | |
| 1302 |
2/2✓ Branch 0 taken 32 times.
✓ Branch 1 taken 16 times.
|
48 | z(4:5) = zstart(4:5, i) |
| 1303 | 16 | call init_sympl(si_cpu(i), f_cpu(i), z, dtaumin, dtaumin, relerr, integmode) | |
| 1304 | 16 | si_gpu(i) = si_cpu(i) | |
| 1305 | 16 | f_gpu(i) = f_cpu(i) | |
| 1306 | |||
| 1307 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 16 times.
|
80 | rhs_y(:, i) = si_cpu(i)%z |
| 1308 | 16 | rhs_mu(i) = f_cpu(i)%mu | |
| 1309 | 16 | rhs_ro0(i) = f_cpu(i)%ro0 | |
| 1310 | call eval_field_booz(f_cpu(i), rhs_y(1, i), rhs_y(2, i), & | ||
| 1311 | 16 | rhs_y(3, i), 0) | |
| 1312 | 16 | call get_derivatives(f_cpu(i), rhs_y(4, i)) | |
| 1313 | 16 | hprime = f_cpu(i)%dH(1)/f_cpu(i)%dpth(1) | |
| 1314 | rhs_host(1, i) = -(f_cpu(i)%dH(2) - & | ||
| 1315 | 16 | f_cpu(i)%hth/f_cpu(i)%hph*f_cpu(i)%dH(3))/f_cpu(i)%dpth(1) | |
| 1316 | 16 | rhs_host(2, i) = hprime | |
| 1317 | 16 | rhs_host(3, i) = (f_cpu(i)%vpar - hprime*f_cpu(i)%hth)/f_cpu(i)%hph | |
| 1318 | 18 | rhs_host(4, i) = -(f_cpu(i)%dH(3) - hprime*f_cpu(i)%dpth(3)) | |
| 1319 | end do | ||
| 1320 | |||
| 1321 | 2 | call evaluate_rhs_gpu(rhs_y, rhs_mu, rhs_ro0, ntestpart, rhs_device) | |
| 1322 |
10/14✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 64 times.
✓ Branch 3 taken 16 times.
✓ Branch 4 taken 62 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 8 times.
✓ Branch 7 taken 54 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
82 | rhs_scale = max(1.0_dp, maxval(abs(rhs_host))) |
| 1323 |
10/14✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 64 times.
✓ Branch 3 taken 16 times.
✓ Branch 4 taken 62 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 6 times.
✓ Branch 7 taken 56 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
82 | rhs_error = maxval(abs(rhs_device - rhs_host))/rhs_scale |
| 1324 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (rhs_error > 1.0e-11_dp) then |
| 1325 | ✗ | rhs_location = maxloc(abs(rhs_device - rhs_host)) | |
| 1326 | ✗ | print '(a,es12.4)', ' compact RHS relative error = ', rhs_error | |
| 1327 | ✗ | print *, 'maximum RHS mismatch at ', rhs_location | |
| 1328 | ✗ | print *, 'host RHS = ', rhs_host(:, rhs_location(2)) | |
| 1329 | ✗ | print *, 'compact RHS = ', rhs_device(:, rhs_location(2)) | |
| 1330 | ✗ | error stop 'compact GPU canonical RHS differs from host reference' | |
| 1331 | end if | ||
| 1332 | |||
| 1333 | ! CPU reference (OpenMP over particles) | ||
| 1334 | 2 | t0 = get_wtime() | |
| 1335 | 2 | !$omp parallel do private(i, it, ktau, ierr) schedule(dynamic) | |
| 1336 | do i = 1, ntestpart | ||
| 1337 | ierr = 0 | ||
| 1338 | cpu_loss(i) = ntimstep | ||
| 1339 | do it = 2, ntimstep | ||
| 1340 | do ktau = 1, ntau_macro(it) | ||
| 1341 | call orbit_timestep_sympl(si_cpu(i), f_cpu(i), ierr) | ||
| 1342 | if (ierr /= 0) exit | ||
| 1343 | end do | ||
| 1344 | if (ierr /= 0) then | ||
| 1345 | cpu_loss(i) = it | ||
| 1346 | exit | ||
| 1347 | end if | ||
| 1348 | end do | ||
| 1349 | cpu_zend(:, i) = si_cpu(i)%z(1:4) | ||
| 1350 | end do | ||
| 1351 | !$omp end parallel do | ||
| 1352 | 2 | t1 = get_wtime() | |
| 1353 | 2 | t_cpu = t1 - t0 | |
| 1354 | |||
| 1355 |
4/6✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
18 | if (any(cpu_loss < ntimstep)) then |
| 1356 | ✗ | error stop 'gpu_mode=benchmark does not support boundary-reaching traces' | |
| 1357 | end if | ||
| 1358 | |||
| 1359 | ! GPU kernel | ||
| 1360 | 2 | t0 = get_wtime() | |
| 1361 | call trace_orbits_gpu_method(si_gpu, f_gpu, ntestpart, ntimstep, & | ||
| 1362 | 2 | ntau_macro, integmode, gpu_loss, gpu_loss_time, gpu_zend, gpu_nfev) | |
| 1363 | 2 | t1 = get_wtime() | |
| 1364 | 2 | t_gpu = t1 - t0 | |
| 1365 | |||
| 1366 | ! Compare | ||
| 1367 | 2 | maxz = 0d0 | |
| 1368 | 2 | loss_mismatch = 0 | |
| 1369 | 2 | cpu_lost = 0 | |
| 1370 | 2 | gpu_lost = 0 | |
| 1371 | 2 | flip = 0 | |
| 1372 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
|
18 | do i = 1, ntestpart |
| 1373 |
6/10✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 64 times.
✓ Branch 7 taken 16 times.
✓ Branch 8 taken 61 times.
✓ Branch 9 taken 3 times.
|
96 | maxz = max(maxz, maxval(dabs(cpu_zend(:, i) - gpu_zend(:, i)))) |
| 1374 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
|
16 | if (cpu_loss(i) /= gpu_loss(i)) loss_mismatch = loss_mismatch + 1 |
| 1375 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
|
16 | if (cpu_loss(i) < ntimstep) cpu_lost = cpu_lost + 1 |
| 1376 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
|
16 | if (gpu_loss(i) < ntimstep) gpu_lost = gpu_lost + 1 |
| 1377 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
|
18 | if ((cpu_loss(i) < ntimstep) .neqv. (gpu_loss(i) < ntimstep)) flip = flip + 1 |
| 1378 | end do | ||
| 1379 | |||
| 1380 | 2 | print *, '==================== GPU vs CPU tracing ====================' | |
| 1381 | 2 | print '(a,i0,a,i0)', ' particles = ', ntestpart, ' timesteps = ', ntimstep | |
| 1382 | 2 | print '(a,i0)', ' integrator mode = ', integmode | |
| 1383 | 2 | print '(a,es12.4)', ' compact RHS relative error = ', rhs_error | |
| 1384 | 2 | print '(a,es12.4)', ' max |z_cpu - z_gpu| (final state) = ', maxz | |
| 1385 | 2 | print '(a,i0,a,i0)', ' loss-step mismatches = ', loss_mismatch, ' / ', ntestpart | |
| 1386 | 2 | print '(a,i0,a,i0,a,f7.4)', ' CPU lost = ', cpu_lost, ' / ', ntestpart, & | |
| 1387 | 4 | ' confined frac = ', 1d0 - real(cpu_lost, dp)/real(ntestpart, dp) | |
| 1388 | 2 | print '(a,i0,a,i0,a,f7.4)', ' GPU lost = ', gpu_lost, ' / ', ntestpart, & | |
| 1389 | 4 | ' confined frac = ', 1d0 - real(gpu_lost, dp)/real(ntestpart, dp) | |
| 1390 | 2 | print '(a,i0,a,i0)', ' lost<->confined flips = ', flip, ' / ', ntestpart | |
| 1391 | 2 | print '(a,f10.4,a)', ' CPU time (OpenMP) = ', t_cpu, ' s' | |
| 1392 | 2 | print '(a,f10.4,a)', ' GPU time = ', t_gpu, ' s' | |
| 1393 |
2/2✓ Branch 2 taken 16 times.
✓ Branch 3 taken 2 times.
|
18 | print '(a,i0)', ' GPU field evaluations = ', sum(gpu_nfev) |
| 1394 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (t_gpu > 0d0) print '(a,f8.2,a)', ' speedup (CPU/GPU) = ', t_cpu/t_gpu, ' x' |
| 1395 | 2 | print *, '============================================================' | |
| 1396 |
15/30✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 2 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 2 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 2 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 2 times.
✗ Branch 29 not taken.
|
2 | end subroutine trace_compare_gpu |
| 1397 | |||
| 1398 | ✗ | subroutine classify_parallel(norb) | |
| 1399 | 2 | use classification, only: trace_orbit_with_classifiers, classification_result_t | |
| 1400 | use params, only: class_passing, class_lost | ||
| 1401 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 1402 | use params, only: debug | ||
| 1403 | #endif | ||
| 1404 | |||
| 1405 | type(tracer_t), intent(inout) :: norb | ||
| 1406 | integer :: i | ||
| 1407 | type(classification_result_t) :: class_result | ||
| 1408 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 1409 | integer :: debug_thread | ||
| 1410 | #endif | ||
| 1411 | |||
| 1412 | ✗ | !$omp parallel firstprivate(norb) private(class_result, i) | |
| 1413 | !$omp do | ||
| 1414 | do i = 1, ntestpart | ||
| 1415 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 1416 | if (debug) then | ||
| 1417 | !$omp critical | ||
| 1418 | kpart = kpart + 1 | ||
| 1419 | debug_thread = 0 | ||
| 1420 | !$ debug_thread = omp_get_thread_num() | ||
| 1421 | print *, kpart, ' / ', ntestpart, 'particle: ', i, 'thread: ', & | ||
| 1422 | debug_thread | ||
| 1423 | !$omp end critical | ||
| 1424 | end if | ||
| 1425 | #endif | ||
| 1426 | |||
| 1427 | if (swcoll) call reset_seed_if_deterministic(particle_stream_id(i)) | ||
| 1428 | call trace_orbit_with_classifiers(norb, i, class_result) | ||
| 1429 | |||
| 1430 | ! Store classification flags in params arrays | ||
| 1431 | class_passing(i) = class_result%passing | ||
| 1432 | class_lost(i) = class_result%lost | ||
| 1433 | ! iclass already populated by trace_orbit_with_classifiers | ||
| 1434 | ! Other results (zend, times_lost, trap_par, perp_inv) also already stored | ||
| 1435 | end do | ||
| 1436 | !$omp end do | ||
| 1437 | !$omp end parallel | ||
| 1438 | ✗ | end subroutine classify_parallel | |
| 1439 | |||
| 1440 | 67 | subroutine print_parameters | |
| 1441 | 67 | print *, 'tau: ', dtau, dtaumin, min(dabs(mod(dtau, dtaumin)), & | |
| 1442 | dabs(mod(dtau, dtaumin) - & | ||
| 1443 | 134 | dtaumin))/dtaumin, ntau | |
| 1444 | 67 | print *, 'v0 = ', v0 | |
| 1445 | ✗ | end subroutine print_parameters | |
| 1446 | |||
| 1447 | 57 | subroutine read_profiles_config(config_file) | |
| 1448 | use simple_profiles, only: read_profiles | ||
| 1449 | |||
| 1450 | character(*), intent(in) :: config_file | ||
| 1451 | |||
| 1452 | 57 | call read_profiles(config_file) | |
| 1453 | 57 | end subroutine read_profiles_config | |
| 1454 | |||
| 1455 | ✗ | subroutine init_collisions | |
| 1456 | use profile_radius, only: init_profile_radius | ||
| 1457 | use params, only: am1, am2, Z1, Z2, particle_energy_eV_effective, & | ||
| 1458 | collision_model, nu_star_standard, lorentz_major_radius_cm, & | ||
| 1459 | lorentz_iota, lorentz_nu, dchichi, slowrate, & | ||
| 1460 | dchichi_norm, slowrate_norm, v0 | ||
| 1461 | use simple_profiles, only: Te_scale, Ti1_scale, Ti2_scale, & | ||
| 1462 | ni1_scale, ni2_scale | ||
| 1463 | |||
| 1464 | real(dp) :: v0_coll, ealpha | ||
| 1465 | real(dp) :: densi1, densi2, tempi1, tempi2, tempe | ||
| 1466 | |||
| 1467 | ✗ | ealpha = particle_energy_eV_effective | |
| 1468 | ✗ | if (trim(collision_model) == 'lorentz_nustar') then | |
| 1469 | call configure_lorentz_nustar(nu_star_standard, & | ||
| 1470 | ✗ | lorentz_major_radius_cm, lorentz_iota, v0) | |
| 1471 | lorentz_nu = nu_star_standard*abs(lorentz_iota)*v0 & | ||
| 1472 | ✗ | /lorentz_major_radius_cm | |
| 1473 | ✗ | dchichi = lorentz_nu | |
| 1474 | ✗ | slowrate = 0.0_dp | |
| 1475 | ✗ | dchichi_norm = lorentz_nu/v0 | |
| 1476 | ✗ | slowrate_norm = 0.0_dp | |
| 1477 | return | ||
| 1478 | end if | ||
| 1479 | ✗ | densi1 = ni1_scale*1.0d-6 | |
| 1480 | ✗ | densi2 = ni2_scale*1.0d-6 | |
| 1481 | ✗ | tempi1 = Ti1_scale | |
| 1482 | ✗ | tempi2 = Ti2_scale | |
| 1483 | ✗ | tempe = Te_scale | |
| 1484 | |||
| 1485 | call loacol_alpha(am1, am2, Z1, Z2, densi1, densi2, tempi1, tempi2, tempe, & | ||
| 1486 | ealpha, v0_coll, dchichi, slowrate, dchichi_norm, & | ||
| 1487 | ✗ | slowrate_norm) | |
| 1488 | |||
| 1489 | ✗ | if (abs(v0_coll - v0) > 1d-6) then | |
| 1490 | ✗ | error stop 'simple_main.init_collisions: v0_coll != v0' | |
| 1491 | end if | ||
| 1492 | |||
| 1493 | ✗ | call init_collision_profiles(am1, am2, Z1, Z2, ealpha, v0) | |
| 1494 | ✗ | call init_profile_radius | |
| 1495 | ✗ | end subroutine init_collisions | |
| 1496 | |||
| 1497 | 48 | subroutine sample_particles(xstart_is_integ_coords) | |
| 1498 | use samplers, only: sample, START_FILE, sample_grid, & | ||
| 1499 | sample_read_integration, & | ||
| 1500 | ✗ | sample_surface_fieldline, sample_surface_fieldline_from_integ, & | |
| 1501 | sample_fusion_birth | ||
| 1502 | |||
| 1503 | logical, intent(in), optional :: xstart_is_integ_coords | ||
| 1504 | logical :: convert_surface_starts | ||
| 1505 | |||
| 1506 | 48 | convert_surface_starts = .false. | |
| 1507 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 42 times.
|
48 | if (present(xstart_is_integ_coords)) then |
| 1508 | 6 | convert_surface_starts = xstart_is_integ_coords | |
| 1509 | end if | ||
| 1510 | |||
| 1511 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 24 times.
|
48 | if (1 == startmode) then |
| 1512 |
3/4✓ Branch 0 taken 1 time.
✓ Branch 1 taken 23 times.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
|
24 | if ((0d0 < grid_density) .and. (1d0 > grid_density)) then |
| 1513 | 1 | call sample_grid(zstart, grid_density, convert_surface_starts) | |
| 1514 | else | ||
| 1515 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 20 times.
|
23 | if (convert_surface_starts) then |
| 1516 | 3 | call sample_surface_fieldline_from_integ(zstart) | |
| 1517 | else | ||
| 1518 | 20 | call sample_surface_fieldline(zstart) | |
| 1519 | end if | ||
| 1520 | end if | ||
| 1521 | |||
| 1522 | elseif (2 == startmode) then | ||
| 1523 | 14 | call sample(zstart, START_FILE) | |
| 1524 | |||
| 1525 | elseif (3 == startmode) then | ||
| 1526 | ✗ | call sample(special_ants_file) | |
| 1527 | |||
| 1528 | elseif (4 == startmode) then | ||
| 1529 | ✗ | call sample(zstart, reuse_batch) | |
| 1530 | |||
| 1531 | elseif (7 == startmode) then | ||
| 1532 | 5 | select case (num_surf) | |
| 1533 | case (0) | ||
| 1534 | call sample_fusion_birth(zstart, 0.0_dp, 1.0_dp, & | ||
| 1535 | 5 | convert_surface_starts) | |
| 1536 | case (1) | ||
| 1537 | call sample_fusion_birth(zstart, 0.0_dp, sbeg(1), & | ||
| 1538 | 1 | convert_surface_starts) | |
| 1539 | case (2) | ||
| 1540 | call sample_fusion_birth(zstart, sbeg(1), sbeg(2), & | ||
| 1541 | 3 | convert_surface_starts) | |
| 1542 | case default | ||
| 1543 |
3/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
9 | error stop 'fusion births require num_surf=0, 1 or 2' |
| 1544 | end select | ||
| 1545 | |||
| 1546 | elseif (5 == startmode) then | ||
| 1547 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (0 == num_surf) then |
| 1548 | 1 | call sample(zstart, 0.0d0, 1.0d0) | |
| 1549 | ✗ | elseif (1 == num_surf) then | |
| 1550 | ✗ | call sample(zstart, 0.0d0, sbeg(1)) | |
| 1551 | ✗ | elseif (2 == num_surf) then | |
| 1552 | ✗ | call sample(zstart, sbeg(1), sbeg(num_surf)) | |
| 1553 | else | ||
| 1554 | ✗ | print *, 'Invalid surface range for volume sample defined ', & | |
| 1555 | ✗ | '(2 < num_surf), stopping.' | |
| 1556 | ✗ | stop | |
| 1557 | end if | ||
| 1558 | |||
| 1559 | elseif (6 == startmode) then | ||
| 1560 | ✗ | call sample_read_integration(zstart, START_FILE) | |
| 1561 | |||
| 1562 | else | ||
| 1563 | ✗ | print *, 'Unknown startmode: ', startmode | |
| 1564 | end if | ||
| 1565 | 47 | end subroutine sample_particles | |
| 1566 | |||
| 1567 | 2 | subroutine sample_particles_test_field | |
| 1568 | !> Sample particles for the analytic circular tokamak TEST field. | ||
| 1569 | !> TEST field uses (r, theta, phi) coordinates with B0=1, R0=1, a=0.5, iota=1. | ||
| 1570 | !> bmod = B0 * (1 - r/R0 * cos(theta)) | ||
| 1571 | use params, only: ntestpart, sbeg, bmod00, bmin, bmax, zstart, & | ||
| 1572 | reset_seed_if_deterministic | ||
| 1573 | |||
| 1574 | real(dp), parameter :: B0 = 1.0d0, R0 = 1.0d0, a = 0.5d0 | ||
| 1575 | real(dp) :: r_start, tmp_rand | ||
| 1576 | integer :: ipart | ||
| 1577 | |||
| 1578 | ! Use sbeg(1) as the starting minor radius (mapped to r for tokamak) | ||
| 1579 | ! sbeg(1) is input as a flux-like value; for TEST field, interpret as r/a | ||
| 1580 | 2 | r_start = sbeg(1)*a | |
| 1581 | |||
| 1582 | ! Set magnetic field bounds for this r value | ||
| 1583 | ! bmod = B0*(1 - r/R0*cos(theta)) | ||
| 1584 | ! Maximum at theta=pi: bmax = B0*(1 + r/R0) | ||
| 1585 | ! Minimum at theta=0: bmin = B0*(1 - r/R0) | ||
| 1586 | 2 | bmod00 = B0 | |
| 1587 | 2 | bmax = B0*(1.0d0 + r_start/R0) | |
| 1588 | 2 | bmin = B0*(1.0d0 - r_start/R0) | |
| 1589 | |||
| 1590 | 2 | print *, 'TEST field: bmod00 = ', bmod00, 'bmin = ', bmin, 'bmax = ', bmax | |
| 1591 | |||
| 1592 | 2 | call reset_seed_if_deterministic | |
| 1593 | |||
| 1594 | ! Sample particles uniformly in (theta, phi) at fixed r | ||
| 1595 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 2 times.
|
66 | do ipart = 1, ntestpart |
| 1596 | 64 | zstart(1, ipart) = r_start | |
| 1597 | 64 | call random_number(tmp_rand) | |
| 1598 | 64 | zstart(2, ipart) = twopi*tmp_rand | |
| 1599 | 64 | call random_number(tmp_rand) | |
| 1600 | 64 | zstart(3, ipart) = twopi*tmp_rand | |
| 1601 | 64 | zstart(4, ipart) = 1.0d0 | |
| 1602 | 64 | call random_number(tmp_rand) | |
| 1603 | 66 | zstart(5, ipart) = 2.0d0*(tmp_rand - 0.5d0) | |
| 1604 | end do | ||
| 1605 | |||
| 1606 | 2 | call save_starting_points_test(zstart) | |
| 1607 | |||
| 1608 | contains | ||
| 1609 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | subroutine save_starting_points_test(zs) |
| 1610 | real(dp), intent(in) :: zs(:, :) | ||
| 1611 | integer :: i, unit | ||
| 1612 | |||
| 1613 | 2 | open (newunit=unit, file='start.dat', recl=1024) | |
| 1614 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 2 times.
|
66 | do i = 1, size(zs, 2) |
| 1615 | 66 | write (unit, *) zs(:, i) | |
| 1616 | end do | ||
| 1617 | 2 | close (unit) | |
| 1618 | 2 | end subroutine save_starting_points_test | |
| 1619 | end subroutine sample_particles_test_field | ||
| 1620 | |||
| 1621 | 7 | subroutine init_bminmax | |
| 1622 | use find_bminmax_sub, only: init_bminmax_arrays | ||
| 1623 | |||
| 1624 | ! Populate bminmax arrays with the active tracing magfie backend. | ||
| 1625 | ! init_starting_surf supplies bmin/bmax for a single sampled surface; | ||
| 1626 | ! this cache is only needed when particles span multiple surfaces. | ||
| 1627 | 7 | call init_bminmax_arrays | |
| 1628 | 7 | end subroutine init_bminmax | |
| 1629 | |||
| 1630 | 68 | logical function needs_bminmax_cache() | |
| 1631 | ! Match the current readers. | ||
| 1632 |
2/2✓ Branch 0 taken 25 times.
✓ Branch 1 taken 43 times.
|
68 | if (classification_enabled()) then |
| 1633 | 25 | needs_bminmax_cache = num_surf > 1 | |
| 1634 | else | ||
| 1635 | 43 | needs_bminmax_cache = num_surf /= 1 | |
| 1636 | end if | ||
| 1637 | 68 | end function needs_bminmax_cache | |
| 1638 | |||
| 1639 | 51 | subroutine init_counters | |
| 1640 | 51 | icounter = 0 ! evaluation counter | |
| 1641 | 51 | kpart = 0 | |
| 1642 | |||
| 1643 | ! initialize array of confined particle percentage | ||
| 1644 |
2/2✓ Branch 0 taken 22979 times.
✓ Branch 1 taken 51 times.
|
23030 | confpart_trap = 0.d0 |
| 1645 |
2/2✓ Branch 0 taken 22979 times.
✓ Branch 1 taken 51 times.
|
23030 | confpart_pass = 0.d0 |
| 1646 |
2/2✓ Branch 0 taken 22979 times.
✓ Branch 1 taken 51 times.
|
23030 | unresolved_orbits = 0 |
| 1647 | |||
| 1648 | ! initialize lost times when particles get lost | ||
| 1649 |
2/2✓ Branch 0 taken 731 times.
✓ Branch 1 taken 51 times.
|
782 | times_lost = -1.d0 |
| 1650 |
2/2✓ Branch 0 taken 731 times.
✓ Branch 1 taken 51 times.
|
782 | orbit_exit_code = 0 |
| 1651 |
2/2✓ Branch 0 taken 731 times.
✓ Branch 1 taken 51 times.
|
782 | boundary_event_radial_residual = -1d0 |
| 1652 |
2/2✓ Branch 0 taken 731 times.
✓ Branch 1 taken 51 times.
|
782 | boundary_event_time_width = -1d0 |
| 1653 | 51 | end subroutine init_counters | |
| 1654 | |||
| 1655 | ✗ | subroutine compute_canonical_frequencies_flat(anorb, initial_state, & | |
| 1656 | n_periods, max_steps, & | ||
| 1657 | values, metadata) | ||
| 1658 | use orbit_frequencies, only: frequency_options_t, frequency_result_t, & | ||
| 1659 | compute_canonical_frequencies | ||
| 1660 | |||
| 1661 | type(tracer_t), intent(inout) :: anorb | ||
| 1662 | real(dp), intent(in) :: initial_state(5) | ||
| 1663 | integer, intent(in) :: n_periods, max_steps | ||
| 1664 | real(dp), intent(out) :: values(6) | ||
| 1665 | integer, intent(out) :: metadata(5) | ||
| 1666 | |||
| 1667 | type(frequency_options_t) :: options | ||
| 1668 | type(frequency_result_t) :: result | ||
| 1669 | real(dp) :: z(5) | ||
| 1670 | |||
| 1671 | ✗ | call ref_to_integ(initial_state(1:3), z(1:3)) | |
| 1672 | ✗ | z(4:5) = initial_state(4:5) | |
| 1673 | |||
| 1674 | ✗ | anorb%dtaumin = dtaumin | |
| 1675 | ✗ | anorb%dtau = dtaumin | |
| 1676 | ✗ | anorb%v0 = v0 | |
| 1677 | ✗ | anorb%relerr = relerr | |
| 1678 | ✗ | anorb%integmode = integmode | |
| 1679 | ✗ | options%n_periods = n_periods | |
| 1680 | ✗ | options%max_steps = max_steps | |
| 1681 | |||
| 1682 | ✗ | call compute_canonical_frequencies(anorb, z, options, result) | |
| 1683 | |||
| 1684 | ✗ | values(1) = result%period | |
| 1685 | ✗ | values(2) = result%period_std | |
| 1686 | ✗ | values(3) = result%delta_phi | |
| 1687 | ✗ | values(4) = result%delta_phi_std | |
| 1688 | ✗ | values(5) = result%omega_b | |
| 1689 | ✗ | values(6) = result%omega_phi | |
| 1690 | ✗ | metadata(1) = result%status | |
| 1691 | ✗ | metadata(2) = result%orbit_class | |
| 1692 | ✗ | metadata(3) = result%parallel_direction | |
| 1693 | ✗ | metadata(4) = result%n_periods | |
| 1694 | ✗ | metadata(5) = result%n_steps | |
| 1695 | ✗ | end subroutine compute_canonical_frequencies_flat | |
| 1696 | |||
| 1697 | ✗ | subroutine invariants_from_state_flat(anorb, initial_state, values, status) | |
| 1698 | use orbit_invariants, only: guiding_center_invariants_t, & | ||
| 1699 | ✗ | invariants_from_state, invariant_flux_convention | |
| 1700 | |||
| 1701 | type(tracer_t), intent(in) :: anorb | ||
| 1702 | real(dp), intent(in) :: initial_state(5) | ||
| 1703 | real(dp), intent(out) :: values(6) | ||
| 1704 | integer, intent(out) :: status | ||
| 1705 | |||
| 1706 | type(guiding_center_invariants_t) :: invariants | ||
| 1707 | type(tracer_t) :: edge_field | ||
| 1708 | real(dp) :: state(5) | ||
| 1709 | |||
| 1710 | ✗ | call ref_to_integ(initial_state(1:3), state(1:3)) | |
| 1711 | ✗ | state(4:5) = initial_state(4:5) | |
| 1712 | ✗ | call invariants_from_state(anorb, state, invariants, status) | |
| 1713 | ✗ | values(1:3) = [invariants%h0, invariants%j_perp, invariants%p_phi] | |
| 1714 | ✗ | call invariant_flux_convention(anorb, values(4), values(5)) | |
| 1715 | ✗ | edge_field = anorb | |
| 1716 | call evaluate_field(edge_field%f, invariant_edge_radius, 0.0_dp, & | ||
| 1717 | ✗ | 0.0_dp, 0) | |
| 1718 | ✗ | values(6) = edge_field%f%Aph | |
| 1719 | ✗ | end subroutine invariants_from_state_flat | |
| 1720 | |||
| 1721 | ✗ | subroutine potato_invariants_to_simple_flat(potato_values, psi_axis, & | |
| 1722 | psi_edge, v0_ratio, simple_values) | ||
| 1723 | ✗ | use orbit_invariants, only: guiding_center_invariants_t, & | |
| 1724 | potato_to_simple_invariants | ||
| 1725 | |||
| 1726 | real(dp), intent(in) :: potato_values(3) | ||
| 1727 | real(dp), intent(in) :: psi_axis, psi_edge, v0_ratio | ||
| 1728 | real(dp), intent(out) :: simple_values(3) | ||
| 1729 | |||
| 1730 | type(guiding_center_invariants_t) :: potato, converted | ||
| 1731 | |||
| 1732 | ✗ | potato%h0 = potato_values(1) | |
| 1733 | ✗ | potato%j_perp = potato_values(2) | |
| 1734 | ✗ | potato%p_phi = potato_values(3) | |
| 1735 | call potato_to_simple_invariants(potato, psi_axis, psi_edge, v0_ratio, & | ||
| 1736 | ✗ | converted) | |
| 1737 | ✗ | simple_values = [converted%h0, converted%j_perp, converted%p_phi] | |
| 1738 | ✗ | end subroutine potato_invariants_to_simple_flat | |
| 1739 | |||
| 1740 | ✗ | subroutine states_from_invariants_flat(anorb, values, max_solutions, & | |
| 1741 | ✗ | states, metadata, residuals, cylindrical, n_solutions, status) | |
| 1742 | use orbit_invariants, only: guiding_center_invariants_t, & | ||
| 1743 | ✗ | invariant_start_result_t, states_from_invariants, & | |
| 1744 | INVARIANT_CAPACITY | ||
| 1745 | |||
| 1746 | type(tracer_t), intent(in) :: anorb | ||
| 1747 | real(dp), intent(in) :: values(3) | ||
| 1748 | integer, intent(in) :: max_solutions | ||
| 1749 | real(dp), intent(out) :: states(5, max_solutions) | ||
| 1750 | integer, intent(out) :: metadata(3, max_solutions) | ||
| 1751 | real(dp), intent(out) :: residuals(max_solutions) | ||
| 1752 | real(dp), intent(out) :: cylindrical(3, max_solutions) | ||
| 1753 | integer, intent(out) :: n_solutions, status | ||
| 1754 | |||
| 1755 | type(guiding_center_invariants_t) :: invariants | ||
| 1756 | ✗ | type(invariant_start_result_t) :: result | |
| 1757 | integer :: i | ||
| 1758 | |||
| 1759 | ✗ | states = 0.0_dp | |
| 1760 | ✗ | metadata = 0 | |
| 1761 | ✗ | residuals = 0.0_dp | |
| 1762 | ✗ | cylindrical = 0.0_dp | |
| 1763 | ✗ | n_solutions = 0 | |
| 1764 | ✗ | invariants%h0 = values(1) | |
| 1765 | ✗ | invariants%j_perp = values(2) | |
| 1766 | ✗ | invariants%p_phi = values(3) | |
| 1767 | ✗ | call states_from_invariants(anorb, invariants, result) | |
| 1768 | ✗ | status = result%status | |
| 1769 | ✗ | n_solutions = min(size(result%starts), max_solutions) | |
| 1770 | ✗ | do i = 1, n_solutions | |
| 1771 | ✗ | states(:, i) = result%starts(i)%state | |
| 1772 | metadata(:, i) = [result%starts(i)%sigma, & | ||
| 1773 | ✗ | result%starts(i)%section_branch, result%starts(i)%section_kind] | |
| 1774 | ✗ | residuals(i) = result%starts(i)%residual | |
| 1775 | ✗ | if (allocated(ref_coords)) then | |
| 1776 | ✗ | call ref_coords%evaluate_cyl(states(1:3, i), cylindrical(:, i)) | |
| 1777 | end if | ||
| 1778 | end do | ||
| 1779 | call sort_invariant_starts(states, metadata, residuals, cylindrical, & | ||
| 1780 | ✗ | n_solutions) | |
| 1781 | ✗ | if (size(result%starts) > max_solutions) status = INVARIANT_CAPACITY | |
| 1782 | ✗ | end subroutine states_from_invariants_flat | |
| 1783 | |||
| 1784 | ✗ | subroutine sort_invariant_starts(states, metadata, residuals, cylindrical, & | |
| 1785 | count) | ||
| 1786 | real(dp), intent(inout) :: states(:, :), residuals(:), cylindrical(:, :) | ||
| 1787 | integer, intent(inout) :: metadata(:, :) | ||
| 1788 | integer, intent(in) :: count | ||
| 1789 | |||
| 1790 | ✗ | real(dp) :: state_tmp(size(states, 1)), residual_tmp | |
| 1791 | ✗ | real(dp) :: cylindrical_tmp(size(cylindrical, 1)) | |
| 1792 | ✗ | integer :: metadata_tmp(size(metadata, 1)) | |
| 1793 | integer :: i, j | ||
| 1794 | |||
| 1795 | ✗ | if (.not. allocated(ref_coords)) return | |
| 1796 | ✗ | do i = 2, count | |
| 1797 | ✗ | state_tmp = states(:, i) | |
| 1798 | ✗ | metadata_tmp = metadata(:, i) | |
| 1799 | ✗ | residual_tmp = residuals(i) | |
| 1800 | ✗ | cylindrical_tmp = cylindrical(:, i) | |
| 1801 | ✗ | j = i - 1 | |
| 1802 | ✗ | do while (j >= 1) | |
| 1803 | ✗ | if (cylindrical(1, j) <= cylindrical_tmp(1)) exit | |
| 1804 | ✗ | states(:, j + 1) = states(:, j) | |
| 1805 | ✗ | metadata(:, j + 1) = metadata(:, j) | |
| 1806 | ✗ | residuals(j + 1) = residuals(j) | |
| 1807 | ✗ | cylindrical(:, j + 1) = cylindrical(:, j) | |
| 1808 | ✗ | j = j - 1 | |
| 1809 | end do | ||
| 1810 | ✗ | states(:, j + 1) = state_tmp | |
| 1811 | ✗ | metadata(:, j + 1) = metadata_tmp | |
| 1812 | ✗ | residuals(j + 1) = residual_tmp | |
| 1813 | ✗ | cylindrical(:, j + 1) = cylindrical_tmp | |
| 1814 | end do | ||
| 1815 | ✗ | end subroutine sort_invariant_starts | |
| 1816 | |||
| 1817 | ✗ | subroutine trace_to_cut_flat(anorb, initial_state, requested_cut_type, & | |
| 1818 | max_events, cut_state, cut_type, status) | ||
| 1819 | use cut_detector, only: cut_detector_t, init_cut_detector => init, & | ||
| 1820 | trace_to_cut | ||
| 1821 | |||
| 1822 | type(tracer_t), intent(inout) :: anorb | ||
| 1823 | real(dp), intent(in) :: initial_state(5) | ||
| 1824 | integer, intent(in) :: requested_cut_type, max_events | ||
| 1825 | real(dp), intent(out) :: cut_state(6) | ||
| 1826 | integer, intent(out) :: cut_type, status | ||
| 1827 | |||
| 1828 | type(cut_detector_t) :: detector | ||
| 1829 | real(dp) :: z(5), reference_position(3) | ||
| 1830 | integer :: event_index, ierr | ||
| 1831 | |||
| 1832 | ✗ | cut_state = 0.0_dp | |
| 1833 | ✗ | cut_type = -1 | |
| 1834 | ✗ | status = 1 | |
| 1835 | ✗ | if (requested_cut_type < -1 .or. requested_cut_type > 1) return | |
| 1836 | ✗ | if (max_events < 1 .or. integmode <= 0) return | |
| 1837 | |||
| 1838 | ✗ | call ref_to_integ(initial_state(1:3), z(1:3)) | |
| 1839 | ✗ | z(4:5) = initial_state(4:5) | |
| 1840 | ✗ | anorb%dtaumin = dtaumin | |
| 1841 | ✗ | anorb%dtau = dtaumin | |
| 1842 | ✗ | anorb%v0 = v0 | |
| 1843 | ✗ | anorb%relerr = relerr | |
| 1844 | ✗ | anorb%integmode = integmode | |
| 1845 | |||
| 1846 | ✗ | call init_sympl(anorb%si, anorb%f, z, dtaumin, dtaumin, relerr, integmode) | |
| 1847 | ✗ | call init_cut_detector(detector, anorb%fper, z) | |
| 1848 | |||
| 1849 | ✗ | do event_index = 1, max_events | |
| 1850 | call trace_to_cut(detector, anorb%si, anorb%f, z, cut_state, & | ||
| 1851 | ✗ | cut_type, ierr) | |
| 1852 | ✗ | if (ierr /= 0) then | |
| 1853 | ✗ | status = ierr | |
| 1854 | ✗ | return | |
| 1855 | end if | ||
| 1856 | ✗ | if (requested_cut_type == -1) exit | |
| 1857 | ✗ | if (cut_type == requested_cut_type) exit | |
| 1858 | end do | ||
| 1859 | |||
| 1860 | ✗ | if (requested_cut_type /= -1) then | |
| 1861 | ✗ | if (cut_type /= requested_cut_type) then | |
| 1862 | ✗ | status = 4 | |
| 1863 | ✗ | return | |
| 1864 | end if | ||
| 1865 | end if | ||
| 1866 | |||
| 1867 | ✗ | call integ_to_ref(cut_state(1:3), reference_position) | |
| 1868 | ✗ | cut_state(1:3) = reference_position | |
| 1869 | ✗ | status = 0 | |
| 1870 | ✗ | end subroutine trace_to_cut_flat | |
| 1871 | 77 | pure integer function classify_orbit_exit(ierr_orbit, model, mode, & | |
| 1872 | boundary_exit_code) | ||
| 1873 | integer, intent(in) :: ierr_orbit, model, mode | ||
| 1874 | integer, intent(in) :: boundary_exit_code | ||
| 1875 | |||
| 1876 |
1/2✓ Branch 0 taken 77 times.
✗ Branch 1 not taken.
|
77 | if (ierr_orbit == 0) then |
| 1877 | classify_orbit_exit = ORBIT_EXIT_COMPLETED | ||
| 1878 |
2/2✓ Branch 0 taken 75 times.
✓ Branch 1 taken 2 times.
|
77 | else if (ierr_orbit == 77) then |
| 1879 | classify_orbit_exit = ORBIT_EXIT_WALL | ||
| 1880 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 72 times.
|
75 | else if (model == ORBIT_FULL_ORBIT) then |
| 1881 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | if (ierr_orbit == ORBIT_FO_LOSS) then |
| 1882 | 2 | classify_orbit_exit = boundary_exit_code | |
| 1883 | else | ||
| 1884 | classify_orbit_exit = ORBIT_EXIT_NUMERICAL_FULL_ORBIT | ||
| 1885 | end if | ||
| 1886 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 68 times.
|
72 | else if (mode <= 0) then |
| 1887 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
|
4 | if (ierr_orbit == 1) then |
| 1888 | 3 | classify_orbit_exit = boundary_exit_code | |
| 1889 | else | ||
| 1890 | classify_orbit_exit = ORBIT_EXIT_NUMERICAL_EVENT | ||
| 1891 | end if | ||
| 1892 | else | ||
| 1893 | 55 | select case (ierr_orbit) | |
| 1894 | case (SYMPLECTIC_STEP_BOUNDARY) | ||
| 1895 | 55 | classify_orbit_exit = boundary_exit_code | |
| 1896 | case (SYMPLECTIC_STEP_OUTSIDE_DOMAIN) | ||
| 1897 | ✗ | classify_orbit_exit = ORBIT_EXIT_NUMERICAL_DOMAIN | |
| 1898 | case (SYMPLECTIC_STEP_MAXITER) | ||
| 1899 | ✗ | classify_orbit_exit = ORBIT_EXIT_NUMERICAL_MAXITER | |
| 1900 | case (SYMPLECTIC_STEP_LINEAR_SOLVE) | ||
| 1901 | ✗ | classify_orbit_exit = ORBIT_EXIT_NUMERICAL_LINEAR | |
| 1902 | case (SYMPLECTIC_STEP_EVENT_NOT_CONVERGED) | ||
| 1903 | 13 | classify_orbit_exit = ORBIT_EXIT_NUMERICAL_EVENT | |
| 1904 | case (SYMPLECTIC_STEP_BOUNDARY_LIMITED) | ||
| 1905 | 13 | classify_orbit_exit = ORBIT_EXIT_NUMERICAL_EVENT | |
| 1906 | case default | ||
| 1907 |
3/5✓ Branch 0 taken 12 times.
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
|
69 | classify_orbit_exit = ORBIT_EXIT_NUMERICAL_EVENT |
| 1908 | end select | ||
| 1909 | end if | ||
| 1910 | ✗ | end function classify_orbit_exit | |
| 1911 | |||
| 1912 | 55 | pure real(dp) function normalized_flux_from_reference_radial( & | |
| 1913 | reference_radial, radial_is_rho) | ||
| 1914 | real(dp), intent(in) :: reference_radial | ||
| 1915 | logical, intent(in) :: radial_is_rho | ||
| 1916 | |||
| 1917 | 55 | normalized_flux_from_reference_radial = reference_radial | |
| 1918 |
3/4✓ Branch 0 taken 20 times.
✓ Branch 1 taken 35 times.
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
|
55 | if (radial_is_rho .and. reference_radial >= 0.0_dp) & |
| 1919 | 20 | normalized_flux_from_reference_radial = reference_radial**2 | |
| 1920 | 55 | end function normalized_flux_from_reference_radial | |
| 1921 | |||
| 1922 | 53 | real(dp) function reference_normalized_flux_s(reference_radial) | |
| 1923 | real(dp), intent(in) :: reference_radial | ||
| 1924 | logical :: radial_is_rho | ||
| 1925 | |||
| 1926 | 53 | radial_is_rho = .false. | |
| 1927 |
1/2✓ Branch 0 taken 53 times.
✗ Branch 1 not taken.
|
53 | if (allocated(ref_coords)) then |
| 1928 | select type (ref_coords) | ||
| 1929 | class is (chartmap_coordinate_system_t) | ||
| 1930 | 19 | radial_is_rho = .true. | |
| 1931 | class default | ||
| 1932 | continue | ||
| 1933 | end select | ||
| 1934 | end if | ||
| 1935 | reference_normalized_flux_s = normalized_flux_from_reference_radial( & | ||
| 1936 | 53 | reference_radial, radial_is_rho) | |
| 1937 | 53 | end function reference_normalized_flux_s | |
| 1938 | |||
| 1939 | 61 | pure logical function should_classify_numerically_confined(exit_code, & | |
| 1940 | model, mode, warning_mode, collisions_enabled, recovery_active, & | ||
| 1941 | normalized_flux_s) | ||
| 1942 | use finite_check, only: is_finite | ||
| 1943 | |||
| 1944 | integer, intent(in) :: exit_code, model, mode | ||
| 1945 | logical, intent(in) :: warning_mode, collisions_enabled, recovery_active | ||
| 1946 | real(dp), intent(in) :: normalized_flux_s | ||
| 1947 | |||
| 1948 | should_classify_numerically_confined = & | ||
| 1949 | exit_code == ORBIT_EXIT_NUMERICAL_MAXITER .and. & | ||
| 1950 | warning_mode .and. .not. collisions_enabled .and. & | ||
| 1951 | model == ORBIT_GC .and. mode > 0 .and. recovery_active .and. & | ||
| 1952 | is_finite(normalized_flux_s) .and. & | ||
| 1953 | normalized_flux_s >= 0.0_dp .and. & | ||
| 1954 |
16/18✓ Branch 0 taken 7 times.
✓ Branch 1 taken 54 times.
✓ Branch 2 taken 6 times.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 1 time.
✓ Branch 8 taken 3 times.
✓ Branch 9 taken 1 time.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 1 time.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 time.
✓ Branch 17 taken 1 time.
|
61 | normalized_flux_s < NUMERICAL_CONFINED_S_MAX |
| 1955 | 61 | end function should_classify_numerically_confined | |
| 1956 | |||
| 1957 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 704 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 704 times.
|
704 | subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) |
| 1958 | use classification, only: trace_orbit_with_classifiers, & | ||
| 1959 | classification_result_t, & | ||
| 1960 | write_classification_results | ||
| 1961 | use magfie_sub, only: SPECTRE | ||
| 1962 | use orbit_symplectic_base, only: symplectic_spectre_newton_warning_factor | ||
| 1963 | use finite_check, only: quiet_nan | ||
| 1964 | |||
| 1965 | type(tracer_t), intent(inout) :: anorb | ||
| 1966 | integer, intent(in) :: ipart | ||
| 1967 | real(dp), intent(out) :: orbit_traj(:, :) ! (5, ntimstep) | ||
| 1968 | real(dp), intent(out) :: orbit_times(:) ! (ntimstep) | ||
| 1969 | |||
| 1970 | real(dp), dimension(5) :: z | ||
| 1971 | real(dp) :: u_ref_prev(3), x_prev(3), x_prev_m(3), exit_step | ||
| 1972 | integer :: it, ierr_orbit, it_final, it_f | ||
| 1973 | integer :: first_unresolved_it, hold_streak | ||
| 1974 | integer(8) :: kt | ||
| 1975 | logical :: passing, faulted, numerical_hold, physical_exit | ||
| 1976 | logical :: numerically_confined | ||
| 1977 | type(rk_recovery_state_t) :: rk_recovery | ||
| 1978 | type(classification_result_t) :: class_result | ||
| 1979 | |||
| 1980 | 704 | ierr_orbit = 0 | |
| 1981 | 704 | faulted = .false. | |
| 1982 | 704 | physical_exit = .false. | |
| 1983 | 704 | numerically_confined = .false. | |
| 1984 | 704 | first_unresolved_it = 0 | |
| 1985 | 704 | hold_streak = 0 | |
| 1986 | rk_recovery = rk_recovery_state_t() | ||
| 1987 | 704 | exit_step = -1d0 | |
| 1988 |
4/4✓ Branch 0 taken 680262 times.
✓ Branch 1 taken 704 times.
✓ Branch 2 taken 3401310 times.
✓ Branch 3 taken 680262 times.
|
4082276 | orbit_traj = quiet_nan() |
| 1989 |
2/2✓ Branch 0 taken 680262 times.
✓ Branch 1 taken 704 times.
|
680966 | orbit_times = quiet_nan() |
| 1990 | |||
| 1991 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 704 times.
|
704 | if (swcoll) call reset_seed_if_deterministic(particle_stream_id(ipart)) |
| 1992 | |||
| 1993 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 656 times.
|
704 | if (classification_enabled()) then |
| 1994 | 48 | call trace_orbit_with_classifiers(anorb, ipart, class_result) | |
| 1995 |
2/2✓ Branch 0 taken 40 times.
✓ Branch 1 taken 8 times.
|
48 | if (class_plot) then |
| 1996 | 40 | call write_classification_results(ipart, class_result) | |
| 1997 | end if | ||
| 1998 | 48 | return | |
| 1999 | end if | ||
| 2000 | |||
| 2001 | 656 | call ref_to_integ(zstart(1:3, ipart), z(1:3)) | |
| 2002 |
2/2✓ Branch 0 taken 1312 times.
✓ Branch 1 taken 656 times.
|
1968 | z(4:5) = zstart(4:5, ipart) |
| 2003 |
2/2✓ Branch 0 taken 3280 times.
✓ Branch 1 taken 656 times.
|
3936 | orbit_traj(:, 1) = zstart(:, ipart) |
| 2004 | 656 | orbit_times(1) = 0.0_dp | |
| 2005 |
2/2✓ Branch 0 taken 3280 times.
✓ Branch 1 taken 656 times.
|
3936 | zend(:, ipart) = 0d0 |
| 2006 | |||
| 2007 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 656 times.
|
656 | if (wall_enabled) then |
| 2008 | ✗ | u_ref_prev = zstart(1:3, ipart) | |
| 2009 | ✗ | call ref_coords%evaluate_cart(u_ref_prev, x_prev) | |
| 2010 | ✗ | x_prev_m = x_prev*chartmap_cart_scale_to_m | |
| 2011 | end if | ||
| 2012 | |||
| 2013 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 656 times.
|
656 | if (orbit_model == ORBIT_FULL_ORBIT) then |
| 2014 | ✗ | if (wall_enabled) error stop & | |
| 2015 | ✗ | 'orbit_model=ORBIT_FULL_ORBIT does not support wall loss yet' | |
| 2016 | ✗ | if (swcoll) error stop & | |
| 2017 | ✗ | 'orbit_model=ORBIT_FULL_ORBIT does not support collisions yet' | |
| 2018 | ✗ | call init_fo(anorb%fo, z, dtaumin) | |
| 2019 |
2/2✓ Branch 0 taken 519 times.
✓ Branch 1 taken 137 times.
|
656 | else if (integmode > 0) then |
| 2020 | 519 | call init_sympl(anorb%si, anorb%f, z, dtaumin, dtaumin, relerr, integmode) | |
| 2021 |
2/2✓ Branch 0 taken 312 times.
✓ Branch 1 taken 207 times.
|
519 | if (isw_field_type == SPECTRE) & |
| 2022 | 312 | anorb%si%warning_factor = symplectic_spectre_newton_warning_factor | |
| 2023 | end if | ||
| 2024 | |||
| 2025 | 656 | call compute_pitch_angle_params(z, passing, trap_par(ipart), perp_inv(ipart)) | |
| 2026 | |||
| 2027 |
2/2✓ Branch 0 taken 380 times.
✓ Branch 1 taken 276 times.
|
656 | if (isw_field_type == SPECTRE) then |
| 2028 |
2/2✓ Branch 0 taken 312 times.
✓ Branch 1 taken 68 times.
|
380 | if (integmode > 0) then |
| 2029 | call trace_orbit_spectre_sympl(anorb, ipart, z, passing, & | ||
| 2030 | 312 | orbit_traj, orbit_times) | |
| 2031 | else | ||
| 2032 | 68 | call trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times) | |
| 2033 | end if | ||
| 2034 | 380 | return | |
| 2035 | end if | ||
| 2036 | |||
| 2037 |
4/4✓ Branch 0 taken 211 times.
✓ Branch 1 taken 65 times.
✓ Branch 2 taken 66 times.
✓ Branch 3 taken 145 times.
|
276 | if (passing .and. should_skip(ipart) .and. isw_field_type /= SPECTRE) then |
| 2038 |
2/2✓ Branch 0 taken 330 times.
✓ Branch 1 taken 66 times.
|
396 | zend(:, ipart) = zstart(:, ipart) |
| 2039 | 66 | times_lost(ipart) = -1.d0 | |
| 2040 | 66 | orbit_exit_code(ipart) = ORBIT_EXIT_SKIPPED | |
| 2041 |
2/2✓ Branch 0 taken 350961 times.
✓ Branch 1 taken 66 times.
|
351027 | do it = 1, ntimstep |
| 2042 | 351027 | !$omp atomic update | |
| 2043 | confpart_pass(it) = confpart_pass(it) + 1.d0 | ||
| 2044 | end do | ||
| 2045 | return | ||
| 2046 | end if | ||
| 2047 | |||
| 2048 | 210 | kt = 0 | |
| 2049 | 210 | it_final = 0 | |
| 2050 |
2/2✓ Branch 0 taken 157 times.
✓ Branch 1 taken 170444 times.
|
170601 | do it = 1, ntimstep |
| 2051 |
2/2✓ Branch 0 taken 170234 times.
✓ Branch 1 taken 210 times.
|
170444 | if (it >= 2) then |
| 2052 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 170234 times.
|
170234 | if (wall_enabled) then |
| 2053 | ! Use microstep-level wall checking for STL walls | ||
| 2054 | call macrostep_with_wall_check(anorb, z, kt, ierr_orbit, & | ||
| 2055 | ntau_macro(it), ipart, x_prev_m, exit_step, & | ||
| 2056 | hold_streak=hold_streak, & | ||
| 2057 | numerical_hold_any=numerical_hold, & | ||
| 2058 | ✗ | rk_recovery=rk_recovery) | |
| 2059 | else | ||
| 2060 | call macrostep(anorb, z, kt, ierr_orbit, ntau_macro(it), & | ||
| 2061 | exit_step, hold_streak=hold_streak, & | ||
| 2062 | numerical_hold_any=numerical_hold, & | ||
| 2063 | 170234 | rk_recovery=rk_recovery) | |
| 2064 | end if | ||
| 2065 | end if | ||
| 2066 | |||
| 2067 |
2/2✓ Branch 0 taken 53 times.
✓ Branch 1 taken 170391 times.
|
170444 | if (ierr_orbit .ne. 0) then |
| 2068 | orbit_exit_code(ipart) = classify_orbit_exit( & | ||
| 2069 | 53 | ierr_orbit, orbit_model, integmode, map_boundary_exit_code) | |
| 2070 | 53 | it_final = it | |
| 2071 | 53 | call integ_to_ref(z(1:3), u_ref_prev) | |
| 2072 | numerically_confined = should_classify_numerically_confined( & | ||
| 2073 | orbit_exit_code(ipart), orbit_model, integmode, & | ||
| 2074 | symplectic_newton_warning_mode, swcoll, rk_recovery%active, & | ||
| 2075 | 53 | reference_normalized_flux_s(u_ref_prev(1))) | |
| 2076 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 53 times.
|
53 | if (numerically_confined) then |
| 2077 | ! All pushers failed from the same validated state in the | ||
| 2078 | ! axis-local core (VMEC s<0.01, or chartmap rho<0.1). | ||
| 2079 | ! This is neither a wall nor an LCFS event. Preserve it as | ||
| 2080 | ! an explicit audited outcome and apply the established | ||
| 2081 | ! numerical-fault-as- | ||
| 2082 | ! confined convention instead of fabricating a loss. | ||
| 2083 | ✗ | orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_CONFINED | |
| 2084 | ✗ | do it_f = it, ntimstep | |
| 2085 | ✗ | call increase_confined_count(it_f, passing) | |
| 2086 | end do | ||
| 2087 | end if | ||
| 2088 | physical_exit = orbit_exit_code(ipart) == ORBIT_EXIT_LCFS .or. & | ||
| 2089 | 53 | orbit_exit_code(ipart) == ORBIT_EXIT_WALL | |
| 2090 |
1/2✓ Branch 0 taken 53 times.
✗ Branch 1 not taken.
|
53 | if (physical_exit) then |
| 2091 |
2/2✓ Branch 0 taken 265 times.
✓ Branch 1 taken 53 times.
|
318 | orbit_traj(:, it) = z |
| 2092 | 53 | orbit_times(it) = exit_step*dtaumin/v0 | |
| 2093 | if (orbit_exit_code(ipart) == ORBIT_EXIT_LCFS .and. & | ||
| 2094 |
5/8✓ Branch 0 taken 53 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 53 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 52 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 52 times.
✗ Branch 7 not taken.
|
53 | orbit_model == ORBIT_GC .and. integmode > 0 .and. & |
| 2095 | ierr_orbit == SYMPLECTIC_STEP_BOUNDARY) then | ||
| 2096 | boundary_event_radial_residual(ipart) = & | ||
| 2097 | 52 | anorb%si%last_event_radial_residual | |
| 2098 | boundary_event_time_width(ipart) = & | ||
| 2099 | 52 | anorb%si%last_event_fraction_width*dtaumin/v0 | |
| 2100 | end if | ||
| 2101 | ✗ | else if (.not. numerically_confined) then | |
| 2102 | if (first_unresolved_it == 0) first_unresolved_it = it | ||
| 2103 | ✗ | faulted = .true. | |
| 2104 | end if | ||
| 2105 | exit | ||
| 2106 | end if | ||
| 2107 | |||
| 2108 | ! Store trajectory data (after macrostep so time is correct) | ||
| 2109 |
2/2✓ Branch 0 taken 851955 times.
✓ Branch 1 taken 170391 times.
|
1022346 | orbit_traj(:, it) = z |
| 2110 | 170391 | orbit_times(it) = kt*dtaumin/v0 | |
| 2111 | |||
| 2112 | 170391 | if (first_unresolved_it == 0) call increase_confined_count(it, passing) | |
| 2113 | 170548 | it_final = it | |
| 2114 | end do | ||
| 2115 | |||
| 2116 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 210 times.
|
210 | if (first_unresolved_it > 0) then |
| 2117 | ✗ | do it_f = first_unresolved_it, ntimstep | |
| 2118 | ✗ | !$omp atomic update | |
| 2119 | unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 | ||
| 2120 | end do | ||
| 2121 | ✗ | physical_exit = .false. | |
| 2122 | ✗ | faulted = .true. | |
| 2123 | ✗ | if (orbit_exit_code(ipart) < ORBIT_EXIT_NUMERICAL_DOMAIN) & | |
| 2124 | ✗ | orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_EVENT | |
| 2125 | end if | ||
| 2126 | |||
| 2127 | ! Fill remaining timesteps with NaN if particle left domain early | ||
| 2128 |
2/2✓ Branch 0 taken 50 times.
✓ Branch 1 taken 160 times.
|
210 | if (it_final < ntimstep) then |
| 2129 |
2/2✓ Branch 0 taken 134817 times.
✓ Branch 1 taken 50 times.
|
134867 | do it = it_final + 1, ntimstep |
| 2130 |
2/2✓ Branch 0 taken 674085 times.
✓ Branch 1 taken 134817 times.
|
808902 | orbit_traj(:, it) = quiet_nan() |
| 2131 | 134867 | orbit_times(it) = quiet_nan() | |
| 2132 | end do | ||
| 2133 | end if | ||
| 2134 | |||
| 2135 | 420 | !$omp critical | |
| 2136 | 210 | call integ_to_ref(z(1:3), zend(1:3, ipart)) | |
| 2137 |
2/2✓ Branch 0 taken 420 times.
✓ Branch 1 taken 210 times.
|
630 | zend(4:5, ipart) = z(4:5) |
| 2138 |
2/2✓ Branch 0 taken 53 times.
✓ Branch 1 taken 157 times.
|
210 | if (physical_exit) then |
| 2139 | 53 | times_lost(ipart) = exit_step*dtaumin/v0 | |
| 2140 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 157 times.
|
157 | else if (numerically_confined) then |
| 2141 | ✗ | times_lost(ipart) = trace_time | |
| 2142 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 157 times.
|
157 | else if (faulted) then |
| 2143 | ✗ | times_lost(ipart) = quiet_nan() | |
| 2144 | else | ||
| 2145 | 157 | times_lost(ipart) = kt*dtaumin/v0 | |
| 2146 | end if | ||
| 2147 | !$omp end critical | ||
| 2148 | 1408 | end subroutine trace_orbit | |
| 2149 | |||
| 2150 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 312 times.
|
312 | subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & |
| 2151 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 312 times.
|
312 | orbit_times) |
| 2152 | !> Per-volume guiding centres switch to full orbit at interior sheets. | ||
| 2153 | 704 | use spectre_sympl_orbit, only: sympl_spectre_state_t, sympl_spectre_reset, & | |
| 2154 | orbit_microstep_sympl_spectre, & | ||
| 2155 | SYMPL_SPECTRE_OK, SYMPL_SPECTRE_LOSS, & | ||
| 2156 | SYMPL_SPECTRE_STOP | ||
| 2157 | use field_can_spectre, only: spectre_mvol, set_spectre_volume_lock | ||
| 2158 | use params, only: crossing_level | ||
| 2159 | use finite_check, only: quiet_nan | ||
| 2160 | |||
| 2161 | type(tracer_t), intent(inout) :: anorb | ||
| 2162 | integer, intent(in) :: ipart | ||
| 2163 | real(dp), intent(inout) :: z(5) | ||
| 2164 | logical, intent(in) :: passing | ||
| 2165 | real(dp), intent(out) :: orbit_traj(:, :) | ||
| 2166 | real(dp), intent(out) :: orbit_times(:) | ||
| 2167 | |||
| 2168 | type(sympl_spectre_state_t) :: state | ||
| 2169 | integer :: it, it_f, ktau, ierr_orbit, it_final | ||
| 2170 | integer :: first_unresolved_it | ||
| 2171 | integer(8) :: kt | ||
| 2172 | real(dp) :: t_stop, t_frac | ||
| 2173 | |||
| 2174 | call sympl_spectre_reset(state, anorb%si, spectre_mvol, integmode, & | ||
| 2175 | 312 | crossing_level) | |
| 2176 | |||
| 2177 | 312 | kt = 0 | |
| 2178 | 312 | it_final = 0 | |
| 2179 | 312 | ierr_orbit = SYMPL_SPECTRE_OK | |
| 2180 | 312 | t_frac = 1.0d0 | |
| 2181 | 312 | first_unresolved_it = 0 | |
| 2182 | |||
| 2183 |
2/2✓ Branch 0 taken 7958 times.
✓ Branch 1 taken 176 times.
|
8134 | do it = 1, ntimstep |
| 2184 |
2/2✓ Branch 0 taken 7646 times.
✓ Branch 1 taken 312 times.
|
7958 | if (it >= 2) then |
| 2185 |
2/2✓ Branch 0 taken 811854 times.
✓ Branch 1 taken 7510 times.
|
819364 | do ktau = 1, ntau_macro(it) |
| 2186 | call orbit_microstep_sympl_spectre(state, anorb%si, anorb%f, & | ||
| 2187 | 811854 | ierr_orbit, t_frac) | |
| 2188 |
2/2✓ Branch 0 taken 136 times.
✓ Branch 1 taken 811718 times.
|
811854 | if (ierr_orbit /= SYMPL_SPECTRE_OK) exit |
| 2189 | 1631218 | kt = kt + 1 | |
| 2190 | end do | ||
| 2191 | end if | ||
| 2192 | |||
| 2193 |
2/2✓ Branch 0 taken 136 times.
✓ Branch 1 taken 7822 times.
|
7958 | if (ierr_orbit /= SYMPL_SPECTRE_OK) then |
| 2194 | 136 | it_final = it | |
| 2195 | 136 | select case (ierr_orbit) | |
| 2196 | case (SYMPL_SPECTRE_LOSS) | ||
| 2197 | 136 | orbit_exit_code(ipart) = ORBIT_EXIT_LCFS | |
| 2198 | case (SYMPL_SPECTRE_STOP) | ||
| 2199 | ✗ | orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_FULL_ORBIT | |
| 2200 | case default | ||
| 2201 |
1/3✓ Branch 0 taken 136 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
136 | orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_EVENT |
| 2202 | end select | ||
| 2203 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 136 times.
|
136 | if (ierr_orbit /= SYMPL_SPECTRE_LOSS .and. & |
| 2204 | ✗ | first_unresolved_it == 0) first_unresolved_it = it | |
| 2205 | exit | ||
| 2206 | end if | ||
| 2207 | |||
| 2208 |
3/4✓ Branch 0 taken 32 times.
✓ Branch 1 taken 7790 times.
✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
|
7822 | if (state%fo%active .and. state%fo%has_y) then |
| 2209 |
2/2✓ Branch 0 taken 160 times.
✓ Branch 1 taken 32 times.
|
192 | z = state%fo%last_y |
| 2210 | else | ||
| 2211 | 7790 | call to_standard_z_coordinates(anorb, z) | |
| 2212 | end if | ||
| 2213 |
2/2✓ Branch 0 taken 39110 times.
✓ Branch 1 taken 7822 times.
|
46932 | orbit_traj(:, it) = z |
| 2214 | 7822 | orbit_times(it) = kt*dtaumin/v0 | |
| 2215 | 7822 | if (first_unresolved_it == 0) call increase_confined_count(it, passing) | |
| 2216 | 7998 | it_final = it | |
| 2217 | end do | ||
| 2218 | |||
| 2219 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 312 times.
|
312 | if (first_unresolved_it > 0) then |
| 2220 | ✗ | do it_f = first_unresolved_it, ntimstep | |
| 2221 | ✗ | !$omp atomic update | |
| 2222 | unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 | ||
| 2223 | end do | ||
| 2224 | end if | ||
| 2225 | |||
| 2226 |
2/2✓ Branch 0 taken 136 times.
✓ Branch 1 taken 176 times.
|
312 | if (it_final < ntimstep) then |
| 2227 |
2/2✓ Branch 0 taken 8762 times.
✓ Branch 1 taken 136 times.
|
8898 | do it = it_final + 1, ntimstep |
| 2228 |
2/2✓ Branch 0 taken 43810 times.
✓ Branch 1 taken 8762 times.
|
52572 | orbit_traj(:, it) = quiet_nan() |
| 2229 | 8898 | orbit_times(it) = quiet_nan() | |
| 2230 | end do | ||
| 2231 | end if | ||
| 2232 | |||
| 2233 |
2/2✓ Branch 0 taken 176 times.
✓ Branch 1 taken 136 times.
|
312 | if (ierr_orbit == SYMPL_SPECTRE_OK) then |
| 2234 | 176 | t_stop = real(kt, dp)*dtaumin/v0 | |
| 2235 | 176 | orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED | |
| 2236 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 136 times.
|
136 | else if (ierr_orbit /= SYMPL_SPECTRE_LOSS) then |
| 2237 | ✗ | t_stop = quiet_nan() | |
| 2238 | else | ||
| 2239 | 136 | t_stop = (real(kt, dp) + t_frac)*dtaumin/v0 | |
| 2240 | end if | ||
| 2241 | |||
| 2242 |
3/4✓ Branch 0 taken 8 times.
✓ Branch 1 taken 304 times.
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
|
312 | if (state%fo%active .and. state%fo%has_y) then |
| 2243 |
2/2✓ Branch 0 taken 40 times.
✓ Branch 1 taken 8 times.
|
48 | z = state%fo%last_y |
| 2244 | else | ||
| 2245 | 304 | call to_standard_z_coordinates(anorb, z) | |
| 2246 | end if | ||
| 2247 | ! The next marker on this thread starts with an unlocked field dispatch. | ||
| 2248 | 312 | call set_spectre_volume_lock(0) | |
| 2249 | |||
| 2250 | 624 | !$omp critical | |
| 2251 | 312 | call integ_to_ref(z(1:3), zend(1:3, ipart)) | |
| 2252 |
2/2✓ Branch 0 taken 624 times.
✓ Branch 1 taken 312 times.
|
936 | zend(4:5, ipart) = z(4:5) |
| 2253 | 312 | times_lost(ipart) = t_stop | |
| 2254 | !$omp end critical | ||
| 2255 |
10/10✓ Branch 0 taken 1248 times.
✓ Branch 1 taken 312 times.
✓ Branch 2 taken 936 times.
✓ Branch 3 taken 312 times.
✓ Branch 4 taken 936 times.
✓ Branch 5 taken 312 times.
✓ Branch 6 taken 936 times.
✓ Branch 7 taken 312 times.
✓ Branch 8 taken 1560 times.
✓ Branch 9 taken 312 times.
|
6240 | end subroutine trace_orbit_spectre_sympl |
| 2256 | |||
| 2257 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 68 times.
|
68 | subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times) |
| 2258 | !> Single-volume RK45 guiding-center trace for SPECTRE. | ||
| 2259 | use spectre_orbit, only: spectre_orbit_state_t, spectre_event_t, & | ||
| 2260 | 312 | spectre_state_reset, orbit_timestep_spectre, & | |
| 2261 | SPECTRE_OK, SPECTRE_BOUNDARY, SPECTRE_FAULT | ||
| 2262 | use magfie_sub, only: spectre_field | ||
| 2263 | use finite_check, only: quiet_nan | ||
| 2264 | |||
| 2265 | integer, intent(in) :: ipart | ||
| 2266 | real(dp), intent(inout) :: z(5) | ||
| 2267 | logical, intent(in) :: passing | ||
| 2268 | real(dp), intent(out) :: orbit_traj(:, :) | ||
| 2269 | real(dp), intent(out) :: orbit_times(:) | ||
| 2270 | |||
| 2271 | type(spectre_orbit_state_t) :: state | ||
| 2272 | type(spectre_event_t) :: event | ||
| 2273 | integer :: it, it_f, ktau, ierr_orbit, it_final | ||
| 2274 | integer(8) :: kt | ||
| 2275 | real(dp) :: t_stop | ||
| 2276 | |||
| 2277 | 68 | call spectre_state_reset(state, spectre_field%data%Mvol) | |
| 2278 | |||
| 2279 | 68 | kt = 0 | |
| 2280 | 68 | it_final = 0 | |
| 2281 | 68 | ierr_orbit = SPECTRE_OK | |
| 2282 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 930 times.
|
934 | do it = 1, ntimstep |
| 2283 |
2/2✓ Branch 0 taken 862 times.
✓ Branch 1 taken 68 times.
|
930 | if (it >= 2) then |
| 2284 |
2/2✓ Branch 0 taken 397982 times.
✓ Branch 1 taken 798 times.
|
398780 | do ktau = 1, ntau_macro(it) |
| 2285 | call orbit_timestep_spectre(state, z, dtaumin, relerr, & | ||
| 2286 | 397982 | ierr_orbit, event) | |
| 2287 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 397982 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
397982 | if (ierr_orbit == SPECTRE_FAULT .and. & |
| 2288 | symplectic_newton_warning_mode) then | ||
| 2289 | ! orbit_timestep_spectre rolls back to the last accepted | ||
| 2290 | ! position. Consume this held interval and retry the same | ||
| 2291 | ! marker on the next microstep; warning mode never turns an | ||
| 2292 | ! isolated RK failure into a terminal marker. | ||
| 2293 | ✗ | call count_event(EVT_WARNING_STEP_SKIP) | |
| 2294 | ✗ | ierr_orbit = SPECTRE_OK | |
| 2295 | ✗ | kt = kt + 1 | |
| 2296 | ✗ | cycle | |
| 2297 | end if | ||
| 2298 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 397918 times.
|
397982 | if (ierr_orbit /= SPECTRE_OK) exit |
| 2299 | 796762 | kt = kt + 1 | |
| 2300 | end do | ||
| 2301 | end if | ||
| 2302 | |||
| 2303 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 866 times.
|
930 | if (ierr_orbit /= SPECTRE_OK) then |
| 2304 | 64 | it_final = it | |
| 2305 |
1/2✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
|
64 | if (ierr_orbit == SPECTRE_BOUNDARY) then |
| 2306 | 64 | orbit_exit_code(ipart) = ORBIT_EXIT_LCFS | |
| 2307 | else | ||
| 2308 | ✗ | orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_EVENT | |
| 2309 | ✗ | do it_f = it, ntimstep | |
| 2310 | ✗ | !$omp atomic update | |
| 2311 | unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 | ||
| 2312 | end do | ||
| 2313 | end if | ||
| 2314 | exit | ||
| 2315 | end if | ||
| 2316 | |||
| 2317 |
2/2✓ Branch 0 taken 4330 times.
✓ Branch 1 taken 866 times.
|
5196 | orbit_traj(:, it) = z |
| 2318 | 866 | orbit_times(it) = kt*dtaumin/v0 | |
| 2319 | 866 | call increase_confined_count(it, passing) | |
| 2320 | 870 | it_final = it | |
| 2321 | end do | ||
| 2322 | |||
| 2323 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 4 times.
|
68 | if (it_final < ntimstep) then |
| 2324 |
2/2✓ Branch 0 taken 6270 times.
✓ Branch 1 taken 64 times.
|
6334 | do it = it_final + 1, ntimstep |
| 2325 |
2/2✓ Branch 0 taken 31350 times.
✓ Branch 1 taken 6270 times.
|
37620 | orbit_traj(:, it) = quiet_nan() |
| 2326 | 6334 | orbit_times(it) = quiet_nan() | |
| 2327 | end do | ||
| 2328 | end if | ||
| 2329 | |||
| 2330 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 4 times.
|
68 | if (ierr_orbit == SPECTRE_BOUNDARY) then |
| 2331 | 64 | t_stop = (real(kt, dp) + event%t_frac)*dtaumin/v0 | |
| 2332 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | else if (ierr_orbit /= SPECTRE_OK) then |
| 2333 | ✗ | t_stop = quiet_nan() | |
| 2334 | else | ||
| 2335 | 4 | t_stop = real(kt, dp)*dtaumin/v0 | |
| 2336 | 4 | orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED | |
| 2337 | end if | ||
| 2338 | |||
| 2339 | 136 | !$omp critical | |
| 2340 | 68 | call integ_to_ref(z(1:3), zend(1:3, ipart)) | |
| 2341 |
2/2✓ Branch 0 taken 136 times.
✓ Branch 1 taken 68 times.
|
204 | zend(4:5, ipart) = z(4:5) |
| 2342 | 68 | times_lost(ipart) = t_stop | |
| 2343 | !$omp end critical | ||
| 2344 | 136 | end subroutine trace_orbit_spectre | |
| 2345 | |||
| 2346 | 8 | pure subroutine locate_linear_lcfs(z_start, z_end, field_period, z_event, & | |
| 2347 | event_fraction) | ||
| 2348 | real(dp), intent(in) :: z_start(5), z_end(5), field_period | ||
| 2349 | real(dp), intent(out) :: z_event(5), event_fraction | ||
| 2350 | real(dp) :: theta_delta, phi_delta | ||
| 2351 | |||
| 2352 | 8 | event_fraction = 1d0 | |
| 2353 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 1 time.
|
8 | if (z_end(1) > z_start(1)) then |
| 2354 | event_fraction = min(1d0, max(0d0, & | ||
| 2355 | 7 | (1d0 - z_start(1))/(z_end(1) - z_start(1)))) | |
| 2356 | end if | ||
| 2357 | theta_delta = modulo(z_end(2) - z_start(2) + 0.5d0*twopi, twopi) - & | ||
| 2358 |
3/4✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 7 times.
|
8 | 0.5d0*twopi |
| 2359 | phi_delta = modulo(z_end(3) - z_start(3) + 0.5d0*field_period, & | ||
| 2360 |
3/4✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 7 times.
|
8 | field_period) - 0.5d0*field_period |
| 2361 | 8 | z_event(1) = 1d0 | |
| 2362 | 8 | z_event(2) = z_start(2) + event_fraction*theta_delta | |
| 2363 | 8 | z_event(3) = z_start(3) + event_fraction*phi_delta | |
| 2364 | z_event(4:5) = z_start(4:5) + & | ||
| 2365 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 8 times.
|
24 | event_fraction*(z_end(4:5) - z_start(4:5)) |
| 2366 | 68 | end subroutine locate_linear_lcfs | |
| 2367 | |||
| 2368 | 7 | subroutine locate_validated_lcfs(z_start, z_end, field_period, z_event, & | |
| 2369 | event_fraction, ierr) | ||
| 2370 | use finite_check, only: is_finite | ||
| 2371 | |||
| 2372 | real(dp), intent(in) :: z_start(5), z_end(5), field_period | ||
| 2373 | real(dp), intent(out) :: z_event(5), event_fraction | ||
| 2374 | integer, intent(inout) :: ierr | ||
| 2375 | |||
| 2376 |
7/10✓ Branch 0 taken 35 times.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 35 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 35 times.
✓ Branch 5 taken 7 times.
✓ Branch 6 taken 35 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 7 times.
|
84 | if (.not. all(is_finite(z_start)) .or. & |
| 2377 | .not. all(is_finite(z_end))) then | ||
| 2378 | ✗ | z_event = z_start | |
| 2379 | ✗ | event_fraction = 0.0_dp | |
| 2380 | ✗ | ierr = 2 | |
| 2381 | ✗ | return | |
| 2382 | end if | ||
| 2383 | |||
| 2384 | call locate_linear_lcfs(z_start, z_end, field_period, z_event, & | ||
| 2385 | 7 | event_fraction) | |
| 2386 | if (all(is_finite(z_event)) .and. & | ||
| 2387 |
8/12✓ Branch 0 taken 35 times.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 35 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 7 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 7 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 5 times.
|
42 | event_fraction >= 0.0_dp .and. event_fraction <= 1.0_dp .and. & |
| 2388 | recovery_state_is_physical(z_event, z_start(4))) return | ||
| 2389 | |||
| 2390 | ! The legacy chamber is the native integration-coordinate surface | ||
| 2391 | ! z(1)=1, which need not coincide with reference rho=1 for a perturbed | ||
| 2392 | ! canonical map. Accept only a physical event inside the attempted | ||
| 2393 | ! microstep; never extrapolate a chamber flag into a loss. | ||
| 2394 | 2 | z_event = z_start | |
| 2395 | 2 | event_fraction = 0.0_dp | |
| 2396 | 2 | ierr = 2 | |
| 2397 | end subroutine locate_validated_lcfs | ||
| 2398 | |||
| 2399 | ✗ | logical function crosses_reference_boundary(z_start, z_end) | |
| 2400 | use finite_check, only: is_finite | ||
| 2401 | |||
| 2402 | real(dp), intent(in) :: z_start(5), z_end(5) | ||
| 2403 | real(dp) :: u_start(3), u_end(3) | ||
| 2404 | |||
| 2405 | ✗ | crosses_reference_boundary = .false. | |
| 2406 | ✗ | if (.not. all(is_finite(z_start)) .or. & | |
| 2407 | .not. all(is_finite(z_end))) return | ||
| 2408 | ✗ | call integ_to_ref(z_start(1:3), u_start) | |
| 2409 | ✗ | call integ_to_ref(z_end(1:3), u_end) | |
| 2410 | ✗ | if (.not. all(is_finite(u_start)) .or. & | |
| 2411 | .not. all(is_finite(u_end))) return | ||
| 2412 | crosses_reference_boundary = u_start(1) < 1.0_dp .and. & | ||
| 2413 | ✗ | u_end(1) >= 1.0_dp | |
| 2414 | end function crosses_reference_boundary | ||
| 2415 | |||
| 2416 | ✗ | subroutine locate_reference_crossing(z_start, z_end, field_period, & | |
| 2417 | z_event, event_fraction, found) | ||
| 2418 | use finite_check, only: is_finite | ||
| 2419 | |||
| 2420 | real(dp), intent(in) :: z_start(5), z_end(5), field_period | ||
| 2421 | real(dp), intent(out) :: z_event(5), event_fraction | ||
| 2422 | logical, intent(out) :: found | ||
| 2423 | real(dp) :: delta(5), z_trial(5), u_start(3), u_end(3), u_trial(3) | ||
| 2424 | real(dp) :: fraction_low, fraction_high, fraction_mid | ||
| 2425 | integer :: iteration | ||
| 2426 | |||
| 2427 | ✗ | found = .false. | |
| 2428 | ✗ | z_event = z_start | |
| 2429 | ✗ | event_fraction = 0.0_dp | |
| 2430 | ✗ | if (.not. all(is_finite(z_start)) .or. & | |
| 2431 | ✗ | .not. all(is_finite(z_end))) return | |
| 2432 | ✗ | call integ_to_ref(z_start(1:3), u_start) | |
| 2433 | ✗ | call integ_to_ref(z_end(1:3), u_end) | |
| 2434 | ✗ | if (.not. all(is_finite(u_start)) .or. & | |
| 2435 | .not. all(is_finite(u_end))) return | ||
| 2436 | ✗ | if (u_start(1) >= 1.0_dp .or. u_end(1) < 1.0_dp) return | |
| 2437 | |||
| 2438 | ✗ | delta = z_end - z_start | |
| 2439 | ✗ | delta(2) = modulo(delta(2) + 0.5_dp*twopi, twopi) - 0.5_dp*twopi | |
| 2440 | delta(3) = modulo(delta(3) + 0.5_dp*field_period, field_period) - & | ||
| 2441 | ✗ | 0.5_dp*field_period | |
| 2442 | ✗ | fraction_low = 0.0_dp | |
| 2443 | ✗ | fraction_high = 1.0_dp | |
| 2444 | ✗ | do iteration = 1, 64 | |
| 2445 | ✗ | fraction_mid = 0.5_dp*(fraction_low + fraction_high) | |
| 2446 | ✗ | z_trial = z_start + fraction_mid*delta | |
| 2447 | ✗ | call integ_to_ref(z_trial(1:3), u_trial) | |
| 2448 | ✗ | if (.not. all(is_finite(u_trial))) return | |
| 2449 | ✗ | if (u_trial(1) >= 1.0_dp) then | |
| 2450 | fraction_high = fraction_mid | ||
| 2451 | else | ||
| 2452 | ✗ | fraction_low = fraction_mid | |
| 2453 | end if | ||
| 2454 | end do | ||
| 2455 | ✗ | event_fraction = fraction_high | |
| 2456 | ✗ | z_event = z_start + event_fraction*delta | |
| 2457 | ✗ | call integ_to_ref(z_event(1:3), u_trial) | |
| 2458 | ✗ | if (.not. all(is_finite(u_trial))) return | |
| 2459 | ✗ | if (abs(u_trial(1) - 1.0_dp) > 0.05_dp) return | |
| 2460 | ✗ | found = .true. | |
| 2461 | end subroutine locate_reference_crossing | ||
| 2462 | |||
| 2463 | 11677873 | subroutine validate_rk_state(z_start, z_end, ierr, expected_momentum) | |
| 2464 | use finite_check, only: is_finite | ||
| 2465 | |||
| 2466 | real(dp), intent(in) :: z_start(5) | ||
| 2467 | real(dp), intent(inout) :: z_end(5) | ||
| 2468 | integer, intent(inout) :: ierr | ||
| 2469 | real(dp), intent(in), optional :: expected_momentum | ||
| 2470 | real(dp) :: momentum_reference | ||
| 2471 | logical :: momentum_valid | ||
| 2472 | |||
| 2473 |
1/2✓ Branch 0 taken 11677873 times.
✗ Branch 1 not taken.
|
11677873 | if (ierr /= 0) return |
| 2474 | 11677873 | momentum_reference = z_start(4) | |
| 2475 | 11677873 | momentum_valid = .true. | |
| 2476 |
2/2✓ Branch 0 taken 7716487 times.
✓ Branch 1 taken 3961386 times.
|
11677873 | if (present(expected_momentum)) then |
| 2477 | 7716487 | momentum_reference = expected_momentum | |
| 2478 | momentum_valid = abs(z_end(4) - momentum_reference) <= & | ||
| 2479 | MOMENTUM_SANITY_FRACTION* & | ||
| 2480 | 7716487 | max(abs(momentum_reference), tiny(1.0_dp)) | |
| 2481 | end if | ||
| 2482 | if (all(is_finite(z_end)) .and. z_end(4) > 0.0_dp .and. & | ||
| 2483 |
9/12✓ Branch 0 taken 58389365 times.
✓ Branch 1 taken 11677873 times.
✓ Branch 2 taken 58389365 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 11677873 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 11677873 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 11677871 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 11 times.
✓ Branch 11 taken 11677860 times.
|
70067238 | abs(z_end(5)) <= 1.0_dp + 100.0_dp*epsilon(1.0_dp) .and. & |
| 2484 | momentum_valid) return | ||
| 2485 | |||
| 2486 | 13 | z_end = z_start | |
| 2487 | 13 | ierr = 2 | |
| 2488 | end subroutine validate_rk_state | ||
| 2489 | |||
| 2490 | 7716493 | subroutine validate_rk_step_or_reference_boundary(z_start, z_candidate, & | |
| 2491 | field_period, expected_momentum, z, event_fraction, ierr, & | ||
| 2492 | reference_boundary) | ||
| 2493 | real(dp), intent(in) :: z_start(5), z_candidate(5), field_period | ||
| 2494 | real(dp), intent(in), optional :: expected_momentum | ||
| 2495 | real(dp), intent(inout) :: z(5) | ||
| 2496 | real(dp), intent(out) :: event_fraction | ||
| 2497 | integer, intent(inout) :: ierr | ||
| 2498 | logical, intent(out) :: reference_boundary | ||
| 2499 | integer :: step_ierr | ||
| 2500 | |||
| 2501 | 7716493 | reference_boundary = .false. | |
| 2502 | 7716493 | event_fraction = 0.0_dp | |
| 2503 | 7716493 | step_ierr = ierr | |
| 2504 |
2/2✓ Branch 0 taken 7716489 times.
✓ Branch 1 taken 4 times.
|
7716493 | if (step_ierr /= 0 .and. step_ierr /= 1) return |
| 2505 | |||
| 2506 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7716489 times.
|
7716489 | if (isw_field_type == ALBERT) then |
| 2507 | ✗ | if (crosses_reference_boundary(z_start, z_candidate)) then | |
| 2508 | call locate_reference_crossing(z_start, z_candidate, & | ||
| 2509 | ✗ | field_period, z, event_fraction, reference_boundary) | |
| 2510 | ✗ | if (reference_boundary) then | |
| 2511 | ✗ | ierr = 0 | |
| 2512 | ✗ | call validate_rk_state(z_start, z, ierr, expected_momentum) | |
| 2513 | ✗ | if (ierr == 0) then | |
| 2514 | ✗ | ierr = 1 | |
| 2515 | ✗ | return | |
| 2516 | end if | ||
| 2517 | ✗ | reference_boundary = .false. | |
| 2518 | end if | ||
| 2519 | end if | ||
| 2520 | end if | ||
| 2521 | |||
| 2522 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 7716485 times.
|
7716489 | if (step_ierr == 1) then |
| 2523 | 4 | z = z_candidate | |
| 2524 | 4 | ierr = step_ierr | |
| 2525 | 4 | return | |
| 2526 | end if | ||
| 2527 | |||
| 2528 | 7716485 | z = z_candidate | |
| 2529 | 7716485 | ierr = 0 | |
| 2530 | 7716485 | call validate_rk_state(z_start, z, ierr, expected_momentum) | |
| 2531 | end subroutine validate_rk_step_or_reference_boundary | ||
| 2532 | |||
| 2533 | 3961383 | subroutine validate_symplectic_endpoint(si_start, anorb, z_start, & | |
| 2534 | z_candidate, z, event_fraction, ierr, reference_boundary) | ||
| 2535 | use field_can_mod, only: get_val | ||
| 2536 | |||
| 2537 | type(symplectic_integrator_t), intent(in) :: si_start | ||
| 2538 | type(tracer_t), intent(in) :: anorb | ||
| 2539 | real(dp), intent(in) :: z_start(5), z_candidate(5) | ||
| 2540 | real(dp), intent(inout) :: z(5) | ||
| 2541 | real(dp), intent(out) :: event_fraction | ||
| 2542 | integer, intent(out) :: ierr | ||
| 2543 | logical, intent(out) :: reference_boundary | ||
| 2544 | type(field_can_t) :: endpoint_field | ||
| 2545 | real(dp) :: z_physical(5), pphi, momentum_reference | ||
| 2546 | integer :: momentum_ierr | ||
| 2547 | |||
| 2548 | 3961383 | ierr = 0 | |
| 2549 | call validate_rk_step_or_reference_boundary(z_start, z_candidate, & | ||
| 2550 | anorb%fper, z=z, event_fraction=event_fraction, ierr=ierr, & | ||
| 2551 | 3961383 | reference_boundary=reference_boundary) | |
| 2552 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3961383 times.
|
3961383 | if (ierr /= 0) then |
| 2553 | ✗ | if (.not. reference_boundary .or. ierr /= 1) return | |
| 2554 | end if | ||
| 2555 | |||
| 2556 | ! Euler's carried field is extrapolated before the angular update; it | ||
| 2557 | ! is not the field at si%z. Validate a copy at the accepted endpoint, | ||
| 2558 | ! without changing the next map's carried field or public coordinates. | ||
| 2559 | 3961383 | endpoint_field = anorb%f | |
| 2560 | 3961383 | pphi = anorb%si%z(4) | |
| 2561 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3961383 times.
|
3961383 | if (reference_boundary) pphi = si_start%z(4) + & |
| 2562 | ✗ | event_fraction*(pphi - si_start%z(4)) | |
| 2563 | 3961383 | call evaluate_field(endpoint_field, z(1), z(2), z(3), 0) | |
| 2564 | 3961383 | call get_val(endpoint_field, pphi) | |
| 2565 | 3961383 | z_physical = z | |
| 2566 | z_physical(4) = sqrt(endpoint_field%mu*endpoint_field%Bmod + & | ||
| 2567 | 3961383 | 0.5_dp*endpoint_field%vpar**2) | |
| 2568 | 3961383 | z_physical(5) = endpoint_field%vpar/(sqrt2*z_physical(4)) | |
| 2569 | 3961383 | momentum_reference = si_start%pabs | |
| 2570 | ! Kinetic momentum is invariant only without electric work. | ||
| 2571 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 3961379 times.
|
3961383 | if (abs(radial_electric_potential_slope) > 0.0_dp) & |
| 2572 | 4 | momentum_reference = z_physical(4) | |
| 2573 | 3961383 | momentum_ierr = 0 | |
| 2574 | call validate_rk_state(z_start, z_physical, momentum_ierr, & | ||
| 2575 | 3961383 | momentum_reference) | |
| 2576 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 3961375 times.
|
3961383 | if (momentum_ierr /= 0) then |
| 2577 | 8 | z = z_start | |
| 2578 | 8 | ierr = momentum_ierr | |
| 2579 | 8 | reference_boundary = .false. | |
| 2580 | end if | ||
| 2581 | 3961383 | end subroutine validate_symplectic_endpoint | |
| 2582 | |||
| 2583 | 170251 | subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step, & | |
| 2584 | hold_streak, numerical_hold_any, rk_recovery) | ||
| 2585 | 3961383 | use alpha_lifetime_sub, only: orbit_timestep_axis | |
| 2586 | use orbit_fo_boris, only: fo_state_t | ||
| 2587 | use orbit_symplectic, only: advance_symplectic_with_retry, & | ||
| 2588 | orbit_timestep_sympl | ||
| 2589 | |||
| 2590 | type(tracer_t), intent(inout) :: anorb | ||
| 2591 | real(dp), intent(inout) :: z(5) | ||
| 2592 | integer(8), intent(inout) :: kt | ||
| 2593 | integer, intent(out) :: ierr_orbit | ||
| 2594 | integer, intent(in) :: ntau_local | ||
| 2595 | real(dp), intent(out), optional :: exit_step | ||
| 2596 | integer, intent(inout), optional :: hold_streak | ||
| 2597 | logical, intent(out), optional :: numerical_hold_any | ||
| 2598 | type(rk_recovery_state_t), intent(inout), optional :: rk_recovery | ||
| 2599 | |||
| 2600 | integer :: hold_streak_local, ktau, recovery_method, symplectic_status | ||
| 2601 | integer :: ierr_coll, collision_hold_start | ||
| 2602 | real(dp) :: z_step_start(5), z_step_end(5), loss_fraction | ||
| 2603 | real(dp) :: u_event(3) | ||
| 2604 | logical :: numerical_hold, numerical_hold_any_local | ||
| 2605 | logical :: reference_boundary | ||
| 2606 | type(rk_recovery_state_t) :: rk_recovery_local | ||
| 2607 | type(symplectic_integrator_t) :: si_step_start | ||
| 2608 | type(field_can_t) :: f_step_start | ||
| 2609 | type(symplectic_integrator_t) :: collision_si_start | ||
| 2610 | type(field_can_t) :: collision_field_start | ||
| 2611 | type(fo_state_t) :: collision_fo_start | ||
| 2612 | type(rk_recovery_state_t) :: collision_recovery_start | ||
| 2613 | |||
| 2614 | 170246 | if (present(exit_step)) exit_step = real(kt, dp) | |
| 2615 | 170251 | hold_streak_local = 0 | |
| 2616 |
2/2✓ Branch 0 taken 170238 times.
✓ Branch 1 taken 13 times.
|
170251 | if (present(hold_streak)) hold_streak_local = hold_streak |
| 2617 | 170251 | numerical_hold_any_local = .false. | |
| 2618 | rk_recovery_local = rk_recovery_state_t() | ||
| 2619 |
2/2✓ Branch 0 taken 170234 times.
✓ Branch 1 taken 17 times.
|
170251 | if (present(rk_recovery)) rk_recovery_local = rk_recovery |
| 2620 | |||
| 2621 |
2/2✓ Branch 0 taken 7716476 times.
✓ Branch 1 taken 170183 times.
|
7886659 | do ktau = 1, ntau_local |
| 2622 | 7716476 | numerical_hold = .false. | |
| 2623 | 7716476 | reference_boundary = .false. | |
| 2624 | 7716476 | z_step_start = z | |
| 2625 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 7716469 times.
|
7716476 | if (swcoll) then |
| 2626 | 7 | collision_si_start = anorb%si | |
| 2627 | 7 | collision_field_start = anorb%f | |
| 2628 | 7 | collision_fo_start = anorb%fo | |
| 2629 | 7 | collision_recovery_start = rk_recovery_local | |
| 2630 | 7 | collision_hold_start = hold_streak_local | |
| 2631 | end if | ||
| 2632 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7716476 times.
|
7716476 | if (orbit_model == ORBIT_FULL_ORBIT) then |
| 2633 | ✗ | call orbit_timestep_fo(anorb%fo, z, ierr_orbit) | |
| 2634 | ✗ | if (ierr_orbit .ne. 0) then | |
| 2635 | ✗ | if (ierr_orbit == ORBIT_FO_LOSS) then | |
| 2636 | ✗ | z_step_end = z | |
| 2637 | call locate_linear_lcfs(z_step_start, z_step_end, & | ||
| 2638 | anorb%fper, & | ||
| 2639 | ✗ | z, loss_fraction) | |
| 2640 | ✗ | if (present(exit_step)) then | |
| 2641 | ✗ | exit_step = real(kt, dp) + loss_fraction | |
| 2642 | end if | ||
| 2643 | ✗ | else if (symplectic_newton_warning_mode) then | |
| 2644 | call hold_isolated_warning_failure(z, z_step_start, & | ||
| 2645 | ierr_orbit, hold_streak_local, numerical_hold, & | ||
| 2646 | ✗ | numerical_hold_any_local) | |
| 2647 | end if | ||
| 2648 | ✗ | if (ierr_orbit .ne. 0) exit | |
| 2649 | end if | ||
| 2650 | ✗ | if (numerical_hold) then | |
| 2651 | ✗ | kt = kt + 1 | |
| 2652 | ✗ | cycle | |
| 2653 | end if | ||
| 2654 | ✗ | hold_streak_local = 0 | |
| 2655 | ✗ | kt = kt + 1 | |
| 2656 | ✗ | cycle | |
| 2657 | end if | ||
| 2658 |
3/4✓ Branch 0 taken 3961376 times.
✓ Branch 1 taken 3755100 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3961376 times.
|
7716476 | if (integmode <= 0 .or. rk_recovery_local%active) then |
| 2659 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3755100 times.
|
3755100 | if (rk_recovery_local%active) then |
| 2660 | call orbit_timestep_recovery(anorb%fo, z, dtaumin, relerr, & | ||
| 2661 | merge(rk_recovery_local%momentum_reference, & | ||
| 2662 | z_step_start(4), & | ||
| 2663 | rk_recovery_local%has_momentum_reference), ierr_orbit, & | ||
| 2664 | ✗ | recovery_method) | |
| 2665 | else | ||
| 2666 | call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, & | ||
| 2667 | 3755100 | ierr_orbit) | |
| 2668 | end if | ||
| 2669 | 3755100 | z_step_end = z | |
| 2670 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3755100 times.
|
3755100 | if (rk_recovery_local%has_momentum_reference) then |
| 2671 | call validate_rk_step_or_reference_boundary(z_step_start, & | ||
| 2672 | z_step_end, anorb%fper, & | ||
| 2673 | rk_recovery_local%momentum_reference, z, loss_fraction, & | ||
| 2674 | ✗ | ierr_orbit, reference_boundary) | |
| 2675 | else | ||
| 2676 | call validate_rk_step_or_reference_boundary(z_step_start, & | ||
| 2677 | z_step_end, anorb%fper, z_step_start(4), z, & | ||
| 2678 | 3755100 | loss_fraction, ierr_orbit, reference_boundary) | |
| 2679 | end if | ||
| 2680 |
3/4✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3755099 times.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
|
3755100 | if (ierr_orbit == 1 .and. .not. reference_boundary) then |
| 2681 | 1 | z_step_end = z | |
| 2682 | call locate_validated_lcfs(z_step_start, z_step_end, & | ||
| 2683 | 1 | anorb%fper, z, loss_fraction, ierr_orbit) | |
| 2684 | end if | ||
| 2685 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3755099 times.
|
3755100 | if (ierr_orbit == 1) then |
| 2686 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (present(exit_step)) then |
| 2687 | 1 | exit_step = real(kt, dp) + loss_fraction | |
| 2688 | end if | ||
| 2689 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (rk_recovery_local%active) then |
| 2690 | ✗ | anorb%si%last_step_fraction = loss_fraction | |
| 2691 | ✗ | if (reference_boundary) then | |
| 2692 | ✗ | call integ_to_ref(z(1:3), u_event) | |
| 2693 | anorb%si%last_event_radial_residual = & | ||
| 2694 | ✗ | abs(u_event(1) - 1.d0) | |
| 2695 | else | ||
| 2696 | anorb%si%last_event_radial_residual = & | ||
| 2697 | ✗ | abs(z(1) - 1.d0) | |
| 2698 | end if | ||
| 2699 | ✗ | anorb%si%last_event_fraction_width = 1.d0 | |
| 2700 | ✗ | ierr_orbit = SYMPLECTIC_STEP_BOUNDARY | |
| 2701 | end if | ||
| 2702 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3755099 times.
|
3755099 | else if (ierr_orbit /= 0) then |
| 2703 | call hold_isolated_warning_failure(z, z_step_start, & | ||
| 2704 | ierr_orbit, hold_streak_local, numerical_hold, & | ||
| 2705 | ✗ | numerical_hold_any_local) | |
| 2706 | end if | ||
| 2707 | else | ||
| 2708 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 3961369 times.
|
3961376 | if (swcoll) call update_momentum(anorb, z) |
| 2709 | 3961376 | si_step_start = anorb%si | |
| 2710 | 3961376 | f_step_start = anorb%f | |
| 2711 | call advance_symplectic_with_retry(anorb%si, anorb%f, & | ||
| 2712 | 3961376 | orbit_timestep_sympl, ierr_orbit) | |
| 2713 |
2/2✓ Branch 0 taken 3961370 times.
✓ Branch 1 taken 6 times.
|
3961376 | if (ierr_orbit == 0 .or. & |
| 2714 | ierr_orbit == SYMPLECTIC_STEP_BOUNDARY) then | ||
| 2715 | 3961370 | symplectic_status = ierr_orbit | |
| 2716 | 3961370 | ierr_orbit = 0 | |
| 2717 | 3961370 | call to_standard_z_coordinates(anorb, z) | |
| 2718 | 3961370 | z_step_end = z | |
| 2719 | call validate_symplectic_endpoint(si_step_start, anorb, & | ||
| 2720 | z_step_start, z_step_end, z, loss_fraction, ierr_orbit, & | ||
| 2721 | 3961370 | reference_boundary) | |
| 2722 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 3961370 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
3961370 | if (reference_boundary .and. ierr_orbit == 1) then |
| 2723 | ✗ | call integ_to_ref(z(1:3), u_event) | |
| 2724 | ✗ | anorb%si%last_step_fraction = loss_fraction | |
| 2725 | anorb%si%last_event_radial_residual = & | ||
| 2726 | ✗ | abs(u_event(1) - 1.d0) | |
| 2727 | ✗ | anorb%si%last_event_fraction_width = 1.d0 | |
| 2728 | ✗ | ierr_orbit = SYMPLECTIC_STEP_BOUNDARY | |
| 2729 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 3961366 times.
|
3961370 | else if (ierr_orbit /= 0) then |
| 2730 | 4 | ierr_orbit = SYMPLECTIC_STEP_EVENT_NOT_CONVERGED | |
| 2731 | 4 | anorb%si = si_step_start | |
| 2732 | 4 | anorb%f = f_step_start | |
| 2733 |
2/2✓ Branch 0 taken 56 times.
✓ Branch 1 taken 3961310 times.
|
3961366 | else if (symplectic_status == SYMPLECTIC_STEP_BOUNDARY) then |
| 2734 | 56 | ierr_orbit = SYMPLECTIC_STEP_BOUNDARY | |
| 2735 | end if | ||
| 2736 | end if | ||
| 2737 |
2/2✓ Branch 0 taken 56 times.
✓ Branch 1 taken 3961320 times.
|
3961376 | if (ierr_orbit == SYMPLECTIC_STEP_BOUNDARY) then |
| 2738 |
1/2✓ Branch 0 taken 56 times.
✗ Branch 1 not taken.
|
56 | if (.not. reference_boundary) & |
| 2739 | 56 | call to_standard_z_coordinates(anorb, z) | |
| 2740 |
1/2✓ Branch 0 taken 56 times.
✗ Branch 1 not taken.
|
56 | if (present(exit_step)) exit_step = real(kt, dp) + & |
| 2741 | 56 | anorb%si%last_step_fraction | |
| 2742 | exit | ||
| 2743 | end if | ||
| 2744 |
4/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 3961310 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 5 times.
|
3961320 | if (ierr_orbit .ne. 0 .and. & |
| 2745 | symplectic_newton_warning_mode) then | ||
| 2746 | ! A failed symplectic map has restored its last accepted | ||
| 2747 | ! state. Advance that same interval with the established | ||
| 2748 | ! adaptive RK/axis path, then thread-locally reseed the | ||
| 2749 | ! symplectic state. Keep this marker on RK while it remains | ||
| 2750 | ! in the failing axis/edge region; generic failures use an | ||
| 2751 | ! exponentially backed-off probe interval. | ||
| 2752 | ! Recover from the reference momentum, not the H-derived | ||
| 2753 | ! z(4) of the symplectic state, so that RK, validation and | ||
| 2754 | ! LCFS interpolation share one start state. | ||
| 2755 | 5 | z_step_start(4) = si_step_start%pabs | |
| 2756 | 5 | z = z_step_start | |
| 2757 | call orbit_timestep_recovery(anorb%fo, z, dtaumin, relerr, & | ||
| 2758 | 5 | si_step_start%pabs, ierr_orbit, recovery_method) | |
| 2759 | 5 | z_step_end = z | |
| 2760 | call validate_rk_step_or_reference_boundary(z_step_start, & | ||
| 2761 | z_step_end, anorb%fper, si_step_start%pabs, z, & | ||
| 2762 | 5 | loss_fraction, ierr_orbit, reference_boundary) | |
| 2763 |
3/4✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
5 | if (ierr_orbit == 1 .and. .not. reference_boundary) then |
| 2764 | 3 | z_step_end = z | |
| 2765 | call locate_validated_lcfs(z_step_start, z_step_end, & | ||
| 2766 | 3 | anorb%fper, z, loss_fraction, ierr_orbit) | |
| 2767 | end if | ||
| 2768 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if (ierr_orbit == 0) then |
| 2769 | call reseed_sympl(anorb%si, anorb%f, z, & | ||
| 2770 | ✗ | si_step_start%pabs) | |
| 2771 | ✗ | call count_event(EVT_SYMPLECTIC_RK_RECOVERY) | |
| 2772 | call activate_rk_recovery(rk_recovery_local, & | ||
| 2773 | ✗ | z_step_start, si_step_start%pabs) | |
| 2774 | ✗ | hold_streak_local = 0 | |
| 2775 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
|
5 | else if (ierr_orbit == 1) then |
| 2776 | 3 | anorb%si%last_step_fraction = loss_fraction | |
| 2777 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if (reference_boundary) then |
| 2778 | ✗ | call integ_to_ref(z(1:3), u_event) | |
| 2779 | anorb%si%last_event_radial_residual = & | ||
| 2780 | ✗ | abs(u_event(1) - 1.d0) | |
| 2781 | else | ||
| 2782 | anorb%si%last_event_radial_residual = & | ||
| 2783 | 3 | abs(z(1) - 1.d0) | |
| 2784 | end if | ||
| 2785 | 3 | anorb%si%last_event_fraction_width = 1.d0 | |
| 2786 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | if (present(exit_step)) exit_step = real(kt, dp) + & |
| 2787 | 3 | loss_fraction | |
| 2788 | 3 | ierr_orbit = SYMPLECTIC_STEP_BOUNDARY | |
| 2789 | 3 | exit | |
| 2790 | else | ||
| 2791 | 2 | z = z_step_start | |
| 2792 | 2 | anorb%si = si_step_start | |
| 2793 | 2 | anorb%f = f_step_start | |
| 2794 | call hold_isolated_warning_failure(z, z_step_start, & | ||
| 2795 | ierr_orbit, hold_streak_local, numerical_hold, & | ||
| 2796 | 2 | numerical_hold_any_local) | |
| 2797 | end if | ||
| 2798 |
2/2✓ Branch 0 taken 3961310 times.
✓ Branch 1 taken 5 times.
|
3961315 | else if (ierr_orbit == 0) then |
| 2799 | 3961310 | call to_standard_z_coordinates(anorb, z) | |
| 2800 | 3961310 | hold_streak_local = 0 | |
| 2801 | 3961310 | rk_recovery_local = rk_recovery_state_t() | |
| 2802 | end if | ||
| 2803 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 3961311 times.
|
3961317 | if (ierr_orbit .ne. 0) exit |
| 2804 | end if | ||
| 2805 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 7716410 times.
|
7716411 | if (ierr_orbit .ne. 0) exit |
| 2806 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 7716409 times.
|
7716410 | if (numerical_hold) then |
| 2807 | 1 | kt = kt + 1 | |
| 2808 | 1 | cycle | |
| 2809 | end if | ||
| 2810 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7716409 times.
|
7716409 | if (rk_recovery_local%active .and. .not. numerical_hold) then |
| 2811 | ✗ | hold_streak_local = 0 | |
| 2812 | call update_recovery_streak(rk_recovery_local, z, & | ||
| 2813 | ✗ | recovery_method) | |
| 2814 | ✗ | if (should_resume_symplectic(rk_recovery_local, z)) then | |
| 2815 | ✗ | if (rk_recovery_local%has_momentum_reference) then | |
| 2816 | call reseed_sympl(anorb%si, anorb%f, z, & | ||
| 2817 | ✗ | rk_recovery_local%momentum_reference) | |
| 2818 | else | ||
| 2819 | ✗ | call reseed_sympl(anorb%si, anorb%f, z) | |
| 2820 | end if | ||
| 2821 | ✗ | rk_recovery_local%active = .false. | |
| 2822 | ✗ | call count_event(EVT_SYMPLECTIC_RESUME) | |
| 2823 | end if | ||
| 2824 | end if | ||
| 2825 |
3/4✓ Branch 0 taken 3961310 times.
✓ Branch 1 taken 3755099 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3961310 times.
|
7716409 | if (integmode <= 0 .or. rk_recovery_local%active) & |
| 2826 | 3755099 | hold_streak_local = 0 | |
| 2827 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 7716407 times.
|
7716409 | if (swcoll) then |
| 2828 | 2 | call collide(z, dtaumin, ierr_coll) | |
| 2829 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (ierr_coll /= 0) then |
| 2830 | ! Reject the complete drift-plus-collision microstep. | ||
| 2831 | ! A collision fault is numerical, never a boundary or hold. | ||
| 2832 | 2 | z = z_step_start | |
| 2833 | 2 | anorb%si = collision_si_start | |
| 2834 | 2 | anorb%f = collision_field_start | |
| 2835 | 2 | anorb%fo = collision_fo_start | |
| 2836 | 2 | rk_recovery_local = collision_recovery_start | |
| 2837 | 2 | hold_streak_local = collision_hold_start | |
| 2838 | 2 | ierr_orbit = SYMPLECTIC_STEP_EVENT_NOT_CONVERGED | |
| 2839 | 2 | exit | |
| 2840 | end if | ||
| 2841 | ✗ | if (rk_recovery_local%active) then | |
| 2842 | ✗ | rk_recovery_local%momentum_reference = z(4) | |
| 2843 | end if | ||
| 2844 | end if | ||
| 2845 | 7886597 | kt = kt + 1 | |
| 2846 | end do | ||
| 2847 |
2/2✓ Branch 0 taken 170238 times.
✓ Branch 1 taken 13 times.
|
170251 | if (present(hold_streak)) hold_streak = hold_streak_local |
| 2848 |
2/2✓ Branch 0 taken 170238 times.
✓ Branch 1 taken 13 times.
|
170251 | if (present(numerical_hold_any)) & |
| 2849 | 170238 | numerical_hold_any = numerical_hold_any_local | |
| 2850 |
2/2✓ Branch 0 taken 170234 times.
✓ Branch 1 taken 17 times.
|
170251 | if (present(rk_recovery)) rk_recovery = rk_recovery_local |
| 2851 |
12/12✓ Branch 0 taken 170251 times.
✓ Branch 1 taken 510753 times.
✓ Branch 2 taken 1021506 times.
✓ Branch 3 taken 170251 times.
✓ Branch 4 taken 170251 times.
✓ Branch 5 taken 510753 times.
✓ Branch 6 taken 170251 times.
✓ Branch 7 taken 510753 times.
✓ Branch 8 taken 510753 times.
✓ Branch 9 taken 170251 times.
✓ Branch 10 taken 170246 times.
✓ Branch 11 taken 5 times.
|
3915773 | end subroutine macrostep |
| 2852 | |||
| 2853 | 4 | subroutine hold_isolated_warning_failure(z, z_step_start, ierr_orbit, & | |
| 2854 | hold_streak, numerical_hold, numerical_hold_any) | ||
| 2855 | real(dp), intent(inout) :: z(5) | ||
| 2856 | real(dp), intent(in) :: z_step_start(5) | ||
| 2857 | integer, intent(inout) :: ierr_orbit, hold_streak | ||
| 2858 | logical, intent(inout) :: numerical_hold, numerical_hold_any | ||
| 2859 | |||
| 2860 | 4 | z = z_step_start | |
| 2861 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (.not. symplectic_newton_warning_mode) return |
| 2862 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
|
4 | if (hold_streak >= max_consecutive_warning_holds) return |
| 2863 | |||
| 2864 | ! Preserve the historical warning-level behavior for one isolated | ||
| 2865 | ! unresolved microstep. Unlike the former permanent latch, the next | ||
| 2866 | ! microstep retries both pushers; a successful step resets the streak. | ||
| 2867 | 2 | call count_event(EVT_WARNING_STEP_SKIP) | |
| 2868 | 2 | hold_streak = hold_streak + 1 | |
| 2869 | 2 | numerical_hold = .true. | |
| 2870 | 2 | numerical_hold_any = .true. | |
| 2871 | 2 | ierr_orbit = 0 | |
| 2872 | 170251 | end subroutine hold_isolated_warning_failure | |
| 2873 | |||
| 2874 | 9 | subroutine locate_wall_segment(z, z_start, z_end, ipart, x_start_m, & | |
| 2875 | u_start, x_end_m, u_end, field_period, start_step, segment_steps, & | ||
| 2876 | hit, hit_step) | ||
| 2877 | real(dp), intent(inout) :: z(5) | ||
| 2878 | real(dp), intent(in) :: z_start(5), z_end(5) | ||
| 2879 | integer, intent(in) :: ipart | ||
| 2880 | real(dp), intent(in) :: x_start_m(3), u_start(3), x_end_m(3), u_end(3) | ||
| 2881 | real(dp), intent(in) :: field_period | ||
| 2882 | real(dp), intent(in) :: start_step, segment_steps | ||
| 2883 | logical, intent(out) :: hit | ||
| 2884 | real(dp), intent(out) :: hit_step | ||
| 2885 | |||
| 2886 | real(dp) :: x_hit_m(3), x_hit(3), normal_m(3), vhat(3) | ||
| 2887 | real(dp) :: segment_length, hit_distance, hit_fraction, cos_inc | ||
| 2888 | real(dp) :: u_hit(3), theta_delta, phi_delta | ||
| 2889 | integer :: ierr_from_cart, triangle | ||
| 2890 | |||
| 2891 | 9 | hit = .false. | |
| 2892 | 9 | hit_step = start_step + segment_steps | |
| 2893 |
2/4✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
|
9 | if (wall_query_rho_lcfs >= 0.0_dp .and. & |
| 2894 | u_end(1) <= wall_query_rho_lcfs) return | ||
| 2895 | call stl_wall_first_hit_segment_with_normal( & | ||
| 2896 | ✗ | wall, x_start_m, x_end_m, hit, x_hit_m, normal_m, triangle) | |
| 2897 | ✗ | if (.not. hit) return | |
| 2898 | |||
| 2899 | ✗ | x_hit = x_hit_m/chartmap_cart_scale_to_m | |
| 2900 | ✗ | wall_hit_triangle(ipart) = triangle | |
| 2901 | ✗ | wall_hit(ipart) = 1_int8 | |
| 2902 | ✗ | wall_hit_cart(:, ipart) = x_hit | |
| 2903 | ✗ | wall_hit_normal_cart(:, ipart) = normal_m | |
| 2904 | |||
| 2905 | ✗ | vhat = x_end_m - x_start_m | |
| 2906 | ✗ | segment_length = sqrt(sum(vhat*vhat)) | |
| 2907 | ✗ | hit_fraction = 1d0 | |
| 2908 | ✗ | if (segment_length > 0d0) then | |
| 2909 | ✗ | hit_distance = sqrt(sum((x_hit_m - x_start_m)**2)) | |
| 2910 | ✗ | hit_fraction = min(1d0, max(0d0, hit_distance/segment_length)) | |
| 2911 | ✗ | vhat = vhat/segment_length | |
| 2912 | ✗ | cos_inc = min(1d0, max(0d0, abs(sum(vhat*normal_m)))) | |
| 2913 | ✗ | wall_hit_cos_incidence(ipart) = cos_inc | |
| 2914 | ✗ | wall_hit_angle_rad(ipart) = acos(cos_inc) | |
| 2915 | end if | ||
| 2916 | |||
| 2917 | ierr_from_cart = 0 | ||
| 2918 | select type (ccs => ref_coords) | ||
| 2919 | class is (chartmap_coordinate_system_t) | ||
| 2920 | ✗ | call ccs%from_cart(x_hit, u_hit, ierr_from_cart) | |
| 2921 | class default | ||
| 2922 | ✗ | ierr_from_cart = 1 | |
| 2923 | end select | ||
| 2924 | |||
| 2925 | ✗ | if (ierr_from_cart /= 0) then | |
| 2926 | theta_delta = modulo(u_end(2) - u_start(2) + 0.5d0*twopi, & | ||
| 2927 | ✗ | twopi) - 0.5d0*twopi | |
| 2928 | phi_delta = modulo(u_end(3) - u_start(3) + 0.5d0*field_period, & | ||
| 2929 | ✗ | field_period) - 0.5d0*field_period | |
| 2930 | ✗ | u_hit(1) = u_start(1) + hit_fraction*(u_end(1) - u_start(1)) | |
| 2931 | ✗ | u_hit(2) = u_start(2) + hit_fraction*theta_delta | |
| 2932 | ✗ | u_hit(3) = u_start(3) + hit_fraction*phi_delta | |
| 2933 | end if | ||
| 2934 | ✗ | call ref_to_integ(u_hit, z(1:3)) | |
| 2935 | ✗ | z(4:5) = z_start(4:5) + hit_fraction*(z_end(4:5) - z_start(4:5)) | |
| 2936 | ✗ | hit_step = start_step + hit_fraction*segment_steps | |
| 2937 | ✗ | end subroutine locate_wall_segment | |
| 2938 | |||
| 2939 | 16 | subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & | |
| 2940 | ipart, x_prev_m, exit_step, hold_streak, numerical_hold_any, & | ||
| 2941 | rk_recovery) | ||
| 2942 | use alpha_lifetime_sub, only: orbit_timestep_axis | ||
| 2943 | use orbit_fo_boris, only: fo_state_t | ||
| 2944 | use orbit_symplectic, only: advance_symplectic_with_retry, & | ||
| 2945 | orbit_timestep_sympl | ||
| 2946 | |||
| 2947 | type(tracer_t), intent(inout) :: anorb | ||
| 2948 | real(dp), intent(inout) :: z(5) | ||
| 2949 | integer(8), intent(inout) :: kt | ||
| 2950 | integer, intent(out) :: ierr_orbit | ||
| 2951 | integer, intent(in) :: ntau_local | ||
| 2952 | integer, intent(in) :: ipart | ||
| 2953 | real(dp), intent(inout) :: x_prev_m(3) | ||
| 2954 | real(dp), intent(out), optional :: exit_step | ||
| 2955 | integer, intent(inout), optional :: hold_streak | ||
| 2956 | logical, intent(out), optional :: numerical_hold_any | ||
| 2957 | type(rk_recovery_state_t), intent(inout), optional :: rk_recovery | ||
| 2958 | |||
| 2959 | integer :: hold_streak_local, ktau, recovery_method, symplectic_status | ||
| 2960 | integer :: ierr_coll, collision_hold_start | ||
| 2961 | real(dp) :: u_ref_prev(3), u_ref_cur(3), x_cur(3), x_cur_m(3) | ||
| 2962 | real(dp) :: z_step_start(5), z_step_end(5), segment_duration | ||
| 2963 | real(dp) :: boundary_fraction | ||
| 2964 | real(dp) :: wall_exit_step | ||
| 2965 | logical :: hit, numerical_hold, numerical_hold_any_local | ||
| 2966 | logical :: reference_boundary | ||
| 2967 | type(rk_recovery_state_t) :: rk_recovery_local | ||
| 2968 | type(symplectic_integrator_t) :: si_step_start | ||
| 2969 | type(field_can_t) :: f_step_start | ||
| 2970 | type(symplectic_integrator_t) :: collision_si_start | ||
| 2971 | type(field_can_t) :: collision_field_start | ||
| 2972 | type(fo_state_t) :: collision_fo_start | ||
| 2973 | type(rk_recovery_state_t) :: collision_recovery_start | ||
| 2974 | |||
| 2975 | 16 | call integ_to_ref(z(1:3), u_ref_prev) | |
| 2976 |
2/2✓ Branch 0 taken 11 times.
✓ Branch 1 taken 5 times.
|
16 | if (present(exit_step)) exit_step = real(kt, dp) |
| 2977 | 16 | hold_streak_local = 0 | |
| 2978 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 12 times.
|
16 | if (present(hold_streak)) hold_streak_local = hold_streak |
| 2979 | 16 | numerical_hold_any_local = .false. | |
| 2980 | rk_recovery_local = rk_recovery_state_t() | ||
| 2981 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
|
16 | if (present(rk_recovery)) rk_recovery_local = rk_recovery |
| 2982 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
|
18 | do ktau = 1, ntau_local |
| 2983 | 16 | numerical_hold = .false. | |
| 2984 | 16 | reference_boundary = .false. | |
| 2985 | 16 | segment_duration = 1.0_dp | |
| 2986 | 16 | z_step_start = z | |
| 2987 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 9 times.
|
16 | if (swcoll) then |
| 2988 | 7 | collision_si_start = anorb%si | |
| 2989 | 7 | collision_field_start = anorb%f | |
| 2990 | 7 | collision_fo_start = anorb%fo | |
| 2991 | 7 | collision_recovery_start = rk_recovery_local | |
| 2992 | 7 | collision_hold_start = hold_streak_local | |
| 2993 | end if | ||
| 2994 |
2/4✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 16 times.
|
16 | if (integmode <= 0 .or. rk_recovery_local%active) then |
| 2995 | ✗ | if (rk_recovery_local%active) then | |
| 2996 | call orbit_timestep_recovery(anorb%fo, z, dtaumin, relerr, & | ||
| 2997 | merge(rk_recovery_local%momentum_reference, & | ||
| 2998 | z_step_start(4), & | ||
| 2999 | rk_recovery_local%has_momentum_reference), ierr_orbit, & | ||
| 3000 | ✗ | recovery_method) | |
| 3001 | else | ||
| 3002 | call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, & | ||
| 3003 | ✗ | ierr_orbit) | |
| 3004 | end if | ||
| 3005 | ✗ | z_step_end = z | |
| 3006 | ✗ | if (rk_recovery_local%has_momentum_reference) then | |
| 3007 | call validate_rk_step_or_reference_boundary(z_step_start, & | ||
| 3008 | z_step_end, anorb%fper, & | ||
| 3009 | rk_recovery_local%momentum_reference, z, & | ||
| 3010 | ✗ | boundary_fraction, ierr_orbit, reference_boundary) | |
| 3011 | else | ||
| 3012 | call validate_rk_step_or_reference_boundary(z_step_start, & | ||
| 3013 | z_step_end, anorb%fper, z_step_start(4), z, & | ||
| 3014 | ✗ | boundary_fraction, ierr_orbit, reference_boundary) | |
| 3015 | end if | ||
| 3016 | ✗ | if (ierr_orbit == 1 .and. .not. reference_boundary) then | |
| 3017 | ✗ | z_step_end = z | |
| 3018 | call locate_validated_lcfs(z_step_start, z_step_end, & | ||
| 3019 | ✗ | anorb%fper, z, boundary_fraction, ierr_orbit) | |
| 3020 | end if | ||
| 3021 | ✗ | if (ierr_orbit == 1) then | |
| 3022 | ✗ | z_step_end = z | |
| 3023 | ✗ | call integ_to_ref(z(1:3), u_ref_cur) | |
| 3024 | ✗ | call ref_coords%evaluate_cart(u_ref_cur, x_cur) | |
| 3025 | ✗ | x_cur_m = x_cur*chartmap_cart_scale_to_m | |
| 3026 | call locate_wall_segment(z, z_step_start, z_step_end, ipart, & | ||
| 3027 | x_prev_m, u_ref_prev, x_cur_m, u_ref_cur, anorb%fper, & | ||
| 3028 | ✗ | real(kt, dp), boundary_fraction, hit, wall_exit_step) | |
| 3029 | ✗ | if (hit) ierr_orbit = 77 | |
| 3030 | ✗ | if (.not. hit .and. rk_recovery_local%active) then | |
| 3031 | ✗ | ierr_orbit = SYMPLECTIC_STEP_BOUNDARY | |
| 3032 | ✗ | anorb%si%last_step_fraction = boundary_fraction | |
| 3033 | ✗ | anorb%si%last_event_radial_residual = abs(z(1) - 1.d0) | |
| 3034 | ✗ | anorb%si%last_event_fraction_width = 1.d0 | |
| 3035 | end if | ||
| 3036 | ✗ | if (present(exit_step)) then | |
| 3037 | ✗ | if (hit) then | |
| 3038 | ✗ | exit_step = wall_exit_step | |
| 3039 | else | ||
| 3040 | ✗ | exit_step = real(kt, dp) + boundary_fraction | |
| 3041 | end if | ||
| 3042 | end if | ||
| 3043 | exit | ||
| 3044 | ✗ | else if (ierr_orbit /= 0) then | |
| 3045 | call hold_isolated_warning_failure(z, z_step_start, & | ||
| 3046 | ierr_orbit, hold_streak_local, numerical_hold, & | ||
| 3047 | ✗ | numerical_hold_any_local) | |
| 3048 | end if | ||
| 3049 | else | ||
| 3050 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 9 times.
|
16 | if (swcoll) call update_momentum(anorb, z) |
| 3051 | 16 | si_step_start = anorb%si | |
| 3052 | 16 | f_step_start = anorb%f | |
| 3053 | call advance_symplectic_with_retry(anorb%si, anorb%f, & | ||
| 3054 | 16 | orbit_timestep_sympl, ierr_orbit, segment_duration) | |
| 3055 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 3 times.
|
16 | if (ierr_orbit == 0 .or. & |
| 3056 | ierr_orbit == SYMPLECTIC_STEP_BOUNDARY) then | ||
| 3057 | 13 | symplectic_status = ierr_orbit | |
| 3058 | 13 | ierr_orbit = 0 | |
| 3059 | 13 | call to_standard_z_coordinates(anorb, z) | |
| 3060 | 13 | z_step_end = z | |
| 3061 | call validate_symplectic_endpoint(si_step_start, anorb, & | ||
| 3062 | z_step_start, z_step_end, z, boundary_fraction, ierr_orbit, & | ||
| 3063 | 13 | reference_boundary) | |
| 3064 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
13 | if (reference_boundary .and. ierr_orbit == 1) then |
| 3065 | ✗ | anorb%si%last_step_fraction = boundary_fraction | |
| 3066 | ✗ | anorb%si%last_event_fraction_width = 1.d0 | |
| 3067 | ✗ | ierr_orbit = SYMPLECTIC_STEP_BOUNDARY | |
| 3068 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 9 times.
|
13 | else if (ierr_orbit /= 0) then |
| 3069 | 4 | ierr_orbit = SYMPLECTIC_STEP_EVENT_NOT_CONVERGED | |
| 3070 | 4 | anorb%si = si_step_start | |
| 3071 | 4 | anorb%f = f_step_start | |
| 3072 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 3 times.
|
9 | else if (symplectic_status == SYMPLECTIC_STEP_BOUNDARY) then |
| 3073 | 6 | ierr_orbit = SYMPLECTIC_STEP_BOUNDARY | |
| 3074 | end if | ||
| 3075 | end if | ||
| 3076 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 10 times.
|
16 | if (ierr_orbit == SYMPLECTIC_STEP_BOUNDARY) then |
| 3077 |
1/2✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
|
6 | if (.not. reference_boundary) & |
| 3078 | 6 | call to_standard_z_coordinates(anorb, z) | |
| 3079 | 6 | call integ_to_ref(z(1:3), u_ref_cur) | |
| 3080 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (reference_boundary) then |
| 3081 | anorb%si%last_event_radial_residual = & | ||
| 3082 | ✗ | abs(u_ref_cur(1) - 1.d0) | |
| 3083 | end if | ||
| 3084 | 6 | call ref_coords%evaluate_cart(u_ref_cur, x_cur) | |
| 3085 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 6 times.
|
24 | x_cur_m = x_cur*chartmap_cart_scale_to_m |
| 3086 | 6 | z_step_end = z | |
| 3087 | 6 | segment_duration = anorb%si%last_step_fraction | |
| 3088 | call locate_wall_segment(z, z_step_start, z_step_end, ipart, & | ||
| 3089 | x_prev_m, u_ref_prev, x_cur_m, u_ref_cur, anorb%fper, & | ||
| 3090 | 6 | real(kt, dp), segment_duration, hit, wall_exit_step) | |
| 3091 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (hit) ierr_orbit = 77 |
| 3092 |
1/2✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
|
6 | if (present(exit_step)) then |
| 3093 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (hit) then |
| 3094 | ✗ | exit_step = wall_exit_step | |
| 3095 | else | ||
| 3096 | 6 | exit_step = real(kt, dp) + segment_duration | |
| 3097 | end if | ||
| 3098 | end if | ||
| 3099 | exit | ||
| 3100 | end if | ||
| 3101 |
4/4✓ Branch 0 taken 7 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 5 times.
|
10 | if (ierr_orbit .ne. 0 .and. & |
| 3102 | symplectic_newton_warning_mode) then | ||
| 3103 | ! Recover from the reference momentum, not the H-derived | ||
| 3104 | ! z(4) of the symplectic state, so that RK, validation and | ||
| 3105 | ! LCFS interpolation share one start state. | ||
| 3106 | 2 | z_step_start(4) = si_step_start%pabs | |
| 3107 | 2 | z = z_step_start | |
| 3108 | call orbit_timestep_recovery(anorb%fo, z, dtaumin, relerr, & | ||
| 3109 | 2 | si_step_start%pabs, ierr_orbit, recovery_method) | |
| 3110 | 2 | z_step_end = z | |
| 3111 | call validate_rk_step_or_reference_boundary(z_step_start, & | ||
| 3112 | z_step_end, anorb%fper, si_step_start%pabs, z, & | ||
| 3113 | 2 | boundary_fraction, ierr_orbit, reference_boundary) | |
| 3114 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
2 | if (ierr_orbit == 1 .and. .not. reference_boundary) then |
| 3115 | ✗ | z_step_end = z | |
| 3116 | call locate_validated_lcfs(z_step_start, z_step_end, & | ||
| 3117 | ✗ | anorb%fper, z, boundary_fraction, ierr_orbit) | |
| 3118 | end if | ||
| 3119 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (ierr_orbit == 0) then |
| 3120 | call reseed_sympl(anorb%si, anorb%f, z, & | ||
| 3121 | ✗ | si_step_start%pabs) | |
| 3122 | ✗ | call count_event(EVT_SYMPLECTIC_RK_RECOVERY) | |
| 3123 | call activate_rk_recovery(rk_recovery_local, & | ||
| 3124 | ✗ | z_step_start, si_step_start%pabs) | |
| 3125 | ✗ | hold_streak_local = 0 | |
| 3126 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | else if (ierr_orbit == 1) then |
| 3127 | ✗ | call integ_to_ref(z(1:3), u_ref_cur) | |
| 3128 | ✗ | call ref_coords%evaluate_cart(u_ref_cur, x_cur) | |
| 3129 | ✗ | x_cur_m = x_cur*chartmap_cart_scale_to_m | |
| 3130 | ✗ | z_step_end = z | |
| 3131 | call locate_wall_segment(z, z_step_start, z_step_end, & | ||
| 3132 | ipart, & | ||
| 3133 | x_prev_m, u_ref_prev, x_cur_m, u_ref_cur, anorb%fper, & | ||
| 3134 | ✗ | real(kt, dp), boundary_fraction, hit, wall_exit_step) | |
| 3135 | ✗ | if (hit) then | |
| 3136 | ✗ | ierr_orbit = 77 | |
| 3137 | else | ||
| 3138 | ✗ | ierr_orbit = SYMPLECTIC_STEP_BOUNDARY | |
| 3139 | ✗ | anorb%si%last_step_fraction = boundary_fraction | |
| 3140 | ✗ | anorb%si%last_event_radial_residual = abs(z(1) - 1.d0) | |
| 3141 | ✗ | anorb%si%last_event_fraction_width = 1.d0 | |
| 3142 | end if | ||
| 3143 | ✗ | if (present(exit_step)) then | |
| 3144 | ✗ | if (hit) then | |
| 3145 | ✗ | exit_step = wall_exit_step | |
| 3146 | else | ||
| 3147 | ✗ | exit_step = real(kt, dp) + boundary_fraction | |
| 3148 | end if | ||
| 3149 | end if | ||
| 3150 | exit | ||
| 3151 | else | ||
| 3152 | 2 | z = z_step_start | |
| 3153 | 2 | anorb%si = si_step_start | |
| 3154 | 2 | anorb%f = f_step_start | |
| 3155 | call hold_isolated_warning_failure(z, z_step_start, & | ||
| 3156 | ierr_orbit, hold_streak_local, numerical_hold, & | ||
| 3157 | 2 | numerical_hold_any_local) | |
| 3158 | end if | ||
| 3159 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 5 times.
|
8 | else if (ierr_orbit == 0) then |
| 3160 | 3 | call to_standard_z_coordinates(anorb, z) | |
| 3161 | 3 | hold_streak_local = 0 | |
| 3162 | 3 | rk_recovery_local = rk_recovery_state_t() | |
| 3163 | end if | ||
| 3164 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 4 times.
|
10 | if (ierr_orbit .ne. 0) exit |
| 3165 | end if | ||
| 3166 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (ierr_orbit .ne. 0) exit |
| 3167 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
|
4 | if (numerical_hold) then |
| 3168 | 1 | kt = kt + 1 | |
| 3169 | 1 | cycle | |
| 3170 | end if | ||
| 3171 | 3 | z_step_end = z | |
| 3172 | 3 | call integ_to_ref(z(1:3), u_ref_cur) | |
| 3173 | 3 | call ref_coords%evaluate_cart(u_ref_cur, x_cur) | |
| 3174 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
|
12 | x_cur_m = x_cur*chartmap_cart_scale_to_m |
| 3175 | call locate_wall_segment(z, z_step_start, z_step_end, ipart, & | ||
| 3176 | x_prev_m, u_ref_prev, x_cur_m, u_ref_cur, anorb%fper, & | ||
| 3177 | 3 | real(kt, dp), segment_duration, hit, wall_exit_step) | |
| 3178 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if (hit) then |
| 3179 | ✗ | ierr_orbit = 77 | |
| 3180 | ✗ | if (present(exit_step)) exit_step = wall_exit_step | |
| 3181 | exit | ||
| 3182 | end if | ||
| 3183 | |||
| 3184 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | if (swcoll) then |
| 3185 | 2 | call collide(z, dtaumin, ierr_coll) | |
| 3186 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (ierr_coll /= 0) then |
| 3187 | ! Reject the complete drift-plus-collision microstep. | ||
| 3188 | ! A collision fault is numerical, never a boundary or hold. | ||
| 3189 | 2 | z = z_step_start | |
| 3190 | 2 | anorb%si = collision_si_start | |
| 3191 | 2 | anorb%f = collision_field_start | |
| 3192 | 2 | anorb%fo = collision_fo_start | |
| 3193 | 2 | rk_recovery_local = collision_recovery_start | |
| 3194 | 2 | hold_streak_local = collision_hold_start | |
| 3195 | 2 | ierr_orbit = SYMPLECTIC_STEP_EVENT_NOT_CONVERGED | |
| 3196 | 2 | exit | |
| 3197 | end if | ||
| 3198 | ✗ | if (rk_recovery_local%active) then | |
| 3199 | ✗ | rk_recovery_local%momentum_reference = z(4) | |
| 3200 | end if | ||
| 3201 | end if | ||
| 3202 | 1 | kt = kt + 1 | |
| 3203 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (rk_recovery_local%active) then |
| 3204 | ✗ | hold_streak_local = 0 | |
| 3205 | call update_recovery_streak(rk_recovery_local, z, & | ||
| 3206 | ✗ | recovery_method) | |
| 3207 | ✗ | if (should_resume_symplectic(rk_recovery_local, z)) then | |
| 3208 | ✗ | if (rk_recovery_local%has_momentum_reference) then | |
| 3209 | call reseed_sympl(anorb%si, anorb%f, z, & | ||
| 3210 | ✗ | rk_recovery_local%momentum_reference) | |
| 3211 | else | ||
| 3212 | ✗ | call reseed_sympl(anorb%si, anorb%f, z) | |
| 3213 | end if | ||
| 3214 | ✗ | rk_recovery_local%active = .false. | |
| 3215 | ✗ | call count_event(EVT_SYMPLECTIC_RESUME) | |
| 3216 | end if | ||
| 3217 | end if | ||
| 3218 |
2/4✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
|
1 | if (integmode <= 0 .or. rk_recovery_local%active) & |
| 3219 | ✗ | hold_streak_local = 0 | |
| 3220 | 1 | x_prev_m = x_cur_m | |
| 3221 | 9 | u_ref_prev = u_ref_cur | |
| 3222 | end do | ||
| 3223 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 12 times.
|
16 | if (present(hold_streak)) hold_streak = hold_streak_local |
| 3224 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 12 times.
|
16 | if (present(numerical_hold_any)) & |
| 3225 | 4 | numerical_hold_any = numerical_hold_any_local | |
| 3226 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
|
16 | if (present(rk_recovery)) rk_recovery = rk_recovery_local |
| 3227 |
10/10✓ Branch 0 taken 16 times.
✓ Branch 1 taken 48 times.
✓ Branch 2 taken 96 times.
✓ Branch 3 taken 16 times.
✓ Branch 4 taken 16 times.
✓ Branch 5 taken 48 times.
✓ Branch 6 taken 16 times.
✓ Branch 7 taken 48 times.
✓ Branch 8 taken 48 times.
✓ Branch 9 taken 16 times.
|
368 | end subroutine macrostep_with_wall_check |
| 3228 | |||
| 3229 | 7930852 | subroutine to_standard_z_coordinates(anorb, z) | |
| 3230 | type(tracer_t), intent(in) :: anorb | ||
| 3231 | real(dp), intent(inout) :: z(5) | ||
| 3232 | |||
| 3233 |
2/2✓ Branch 0 taken 23792556 times.
✓ Branch 1 taken 7930852 times.
|
31723408 | z(1:3) = anorb%si%z(1:3) |
| 3234 | 7930852 | z(4) = dsqrt(anorb%f%mu*anorb%f%Bmod + 0.5d0*anorb%f%vpar**2) | |
| 3235 | 7930852 | z(5) = anorb%f%vpar/(z(4)*sqrt2) | |
| 3236 | 16 | end subroutine to_standard_z_coordinates | |
| 3237 | |||
| 3238 | 179079 | subroutine increase_confined_count(it, passing) | |
| 3239 | integer, intent(in) :: it | ||
| 3240 | logical, intent(in) :: passing | ||
| 3241 | |||
| 3242 |
2/2✓ Branch 0 taken 138880 times.
✓ Branch 1 taken 40199 times.
|
179079 | if (passing) then |
| 3243 | 138880 | !$omp atomic update | |
| 3244 | confpart_pass(it) = confpart_pass(it) + 1.d0 | ||
| 3245 | else | ||
| 3246 | 40199 | !$omp atomic update | |
| 3247 | confpart_trap(it) = confpart_trap(it) + 1.d0 | ||
| 3248 | end if | ||
| 3249 | 179079 | end subroutine increase_confined_count | |
| 3250 | |||
| 3251 | 656 | subroutine compute_pitch_angle_params(z, passing, trap_par_, perp_inv_) | |
| 3252 | use find_bminmax_sub, only: get_bminmax | ||
| 3253 | |||
| 3254 | real(dp), intent(in) :: z(5) | ||
| 3255 | logical, intent(out) :: passing | ||
| 3256 | real(dp), intent(out) :: trap_par_, perp_inv_ | ||
| 3257 | |||
| 3258 | real(dp) :: bmod | ||
| 3259 | |||
| 3260 | 1312 | !$omp critical | |
| 3261 | 656 | bmod = compute_bmod(z(1:3)) | |
| 3262 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 646 times.
|
656 | if (num_surf /= 1) then |
| 3263 | 10 | call get_bminmax(z(1), bmin, bmax) | |
| 3264 | end if | ||
| 3265 | 656 | passing = z(5)**2 .gt. 1.d0 - bmod/bmax | |
| 3266 | 656 | trap_par_ = ((1.d0 - z(5)**2)*bmax/bmod - 1.d0)*bmin/(bmax - bmin) | |
| 3267 | 656 | perp_inv_ = z(4)**2*(1.d0 - z(5)**2)/bmod | |
| 3268 | !$omp end critical | ||
| 3269 | 656 | end subroutine compute_pitch_angle_params | |
| 3270 | |||
| 3271 | 656 | function compute_bmod(z) result(bmod) | |
| 3272 | use magfie_sub, only: magfie | ||
| 3273 | |||
| 3274 | real(dp) :: bmod | ||
| 3275 | real(dp), intent(in) :: z(3) | ||
| 3276 | |||
| 3277 | real(dp) :: sqrtg, bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 3278 | |||
| 3279 | 656 | call magfie(z(1:3), bmod, sqrtg, bder, hcovar, hctrvr, hcurl) | |
| 3280 | 656 | end function compute_bmod | |
| 3281 | |||
| 3282 | 24 | subroutine update_momentum(anorb, z) | |
| 3283 | 656 | use orbit_symplectic, only: get_val | |
| 3284 | |||
| 3285 | type(tracer_t), intent(inout) :: anorb | ||
| 3286 | real(dp), intent(in) :: z(5) | ||
| 3287 | |||
| 3288 | 24 | anorb%si%pabs = z(4) | |
| 3289 | 24 | anorb%f%vpar = z(4)*z(5)*sqrt2 | |
| 3290 | 24 | anorb%f%mu = z(4)**2*(1.d0 - z(5)**2)/anorb%f%Bmod | |
| 3291 | 24 | anorb%si%z(4) = anorb%f%vpar*anorb%f%hph + anorb%f%Aph/anorb%f%ro0 | |
| 3292 | 24 | call get_val(anorb%f, anorb%si%z(4)) ! for pth | |
| 3293 | 24 | anorb%si%pthold = anorb%f%pth | |
| 3294 | 24 | end subroutine update_momentum | |
| 3295 | |||
| 3296 | 4 | subroutine collide(z, dt, ierr_coll) | |
| 3297 | 24 | use profile_radius, only: profile_s | |
| 3298 | real(dp), intent(inout) :: z(5) | ||
| 3299 | real(dp), intent(in) :: dt | ||
| 3300 | integer, intent(out) :: ierr_coll | ||
| 3301 | |||
| 3302 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (trim(collision_model) == 'lorentz_nustar') then |
| 3303 | ✗ | call stost_lorentz(z, dt, ierr_coll) | |
| 3304 | else | ||
| 3305 | 4 | call stost(z, dt, 1, ierr_coll, s_profile=profile_s(z(1:3))) | |
| 3306 | end if | ||
| 3307 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (ierr_coll /= 0) then |
| 3308 | 4 | print *, 'Error in stost: ', ierr_coll, 'z = ', z, 'dtaumin = ', dtaumin | |
| 3309 | end if | ||
| 3310 | 4 | end subroutine collide | |
| 3311 | |||
| 3312 | 49 | subroutine write_output | |
| 3313 | use field_can_base, only: n_field_evaluations, n_field_evaluations_d1, & | ||
| 3314 | n_field_evaluations_d2 | ||
| 3315 | use params, only: output_results_netcdf | ||
| 3316 | use netcdf_results_output, only: write_results_netcdf | ||
| 3317 | |||
| 3318 | integer(8) :: total_field_evaluations | ||
| 3319 | integer(8) :: total_field_evaluations_d1, total_field_evaluations_d2 | ||
| 3320 | |||
| 3321 | ! Sum field evaluations across all threads | ||
| 3322 | total_field_evaluations = 0 | ||
| 3323 | total_field_evaluations_d1 = 0 | ||
| 3324 | total_field_evaluations_d2 = 0 | ||
| 3325 | !$omp parallel reduction(+:total_field_evaluations, & | ||
| 3326 | !$omp& total_field_evaluations_d1, & | ||
| 3327 | 49 | !$omp& total_field_evaluations_d2) | |
| 3328 | total_field_evaluations = total_field_evaluations + n_field_evaluations | ||
| 3329 | total_field_evaluations_d1 = total_field_evaluations_d1 + & | ||
| 3330 | n_field_evaluations_d1 | ||
| 3331 | total_field_evaluations_d2 = total_field_evaluations_d2 + & | ||
| 3332 | n_field_evaluations_d2 | ||
| 3333 | !$omp end parallel | ||
| 3334 | |||
| 3335 | 49 | print *, "Total field evaluations: ", total_field_evaluations | |
| 3336 | ! Split by derivative depth: first derivatives only (what explicit | ||
| 3337 | ! Runge-Kutta paths need) versus second derivatives (what the implicit | ||
| 3338 | ! symplectic Jacobians need). The two cost different amounts, so a fair | ||
| 3339 | ! cost-per-accuracy comparison has to weight them separately. | ||
| 3340 | 49 | print *, " first-derivative (mode_secders=0): ", total_field_evaluations_d1 | |
| 3341 | 49 | print *, " second-derivative (mode_secders>0): ", total_field_evaluations_d2 | |
| 3342 | |||
| 3343 | 49 | call write_results | |
| 3344 | |||
| 3345 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 48 times.
|
49 | if (output_results_netcdf) then |
| 3346 | 1 | call write_results_netcdf('results.nc', wall%facet_count) | |
| 3347 | end if | ||
| 3348 | 49 | end subroutine write_output | |
| 3349 | |||
| 3350 | 52 | subroutine write_results | |
| 3351 | !> Write the per-particle and confined-fraction result files from the | ||
| 3352 | !> shared result arrays after particle tracing has reached a team-safe | ||
| 3353 | !> point. Confined fractions are conditional on the numerically resolved | ||
| 3354 | !> population at each time. Particles not yet finished (for example in an | ||
| 3355 | !> explicitly requested partial write) keep times_lost = -1. | ||
| 3356 | 49 | use finite_check, only: quiet_nan | |
| 3357 | |||
| 3358 | integer :: i, num_lost, resolved_count, unit | ||
| 3359 | real(dp) :: inverse_times_lost_sum, norm, pass_fraction, trap_fraction | ||
| 3360 | |||
| 3361 | 52 | norm = real(max(ntestpart, 1), dp) | |
| 3362 | |||
| 3363 | open (newunit=unit, file='orbit_exit_code.dat', status='replace', & | ||
| 3364 | 52 | action='write') | |
| 3365 | 52 | close (unit) | |
| 3366 | |||
| 3367 | open (newunit=unit, file='times_lost.dat', status='replace', & | ||
| 3368 | 52 | action='write', recl=1024) | |
| 3369 | 52 | num_lost = 0 | |
| 3370 | 52 | inverse_times_lost_sum = 0.0d0 | |
| 3371 |
2/2✓ Branch 0 taken 724 times.
✓ Branch 1 taken 52 times.
|
776 | do i = 1, ntestpart |
| 3372 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 723 times.
|
724 | if (orbit_exit_code(i) >= ORBIT_EXIT_NUMERICAL_DOMAIN) then |
| 3373 | 1 | times_lost(i) = quiet_nan() | |
| 3374 | end if | ||
| 3375 | 724 | write (unit, *) i, times_lost(i), trap_par(i), zstart(1, i), & | |
| 3376 | 1448 | perp_inv(i), zend(:, i) | |
| 3377 |
2/2✓ Branch 0 taken 256 times.
✓ Branch 1 taken 468 times.
|
724 | if (orbit_exit_code(i) == ORBIT_EXIT_LCFS .or. & |
| 3378 | 52 | orbit_exit_code(i) == ORBIT_EXIT_WALL) then | |
| 3379 |
1/2✓ Branch 0 taken 256 times.
✗ Branch 1 not taken.
|
256 | if (times_lost(i) > 0d0) then |
| 3380 | 256 | num_lost = num_lost + 1 | |
| 3381 | inverse_times_lost_sum = inverse_times_lost_sum + & | ||
| 3382 | 256 | 1.0/times_lost(i) | |
| 3383 | end if | ||
| 3384 | end if | ||
| 3385 | end do | ||
| 3386 | 52 | close (unit) | |
| 3387 | |||
| 3388 | open (newunit=unit, file='orbit_exit_code.dat', status='replace', & | ||
| 3389 | 52 | action='write', recl=1024) | |
| 3390 |
2/2✓ Branch 0 taken 724 times.
✓ Branch 1 taken 52 times.
|
776 | do i = 1, ntestpart |
| 3391 | 724 | write (unit, *) i, orbit_exit_code(i), & | |
| 3392 | 724 | times_lost(i), boundary_event_radial_residual(i), & | |
| 3393 | 1500 | boundary_event_time_width(i) | |
| 3394 | end do | ||
| 3395 | 52 | close (unit) | |
| 3396 | |||
| 3397 | open (newunit=unit, file='avg_inverse_t_lost.dat', status='replace', & | ||
| 3398 | 52 | action='write', recl=1024) | |
| 3399 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 28 times.
|
52 | if (num_lost > 0) then |
| 3400 | 24 | write (unit, *) inverse_times_lost_sum/num_lost | |
| 3401 | else | ||
| 3402 | 28 | write (unit, *) 0d0 | |
| 3403 | end if | ||
| 3404 | 52 | close (unit) | |
| 3405 | |||
| 3406 | open (newunit=unit, file='confined_fraction.dat', status='replace', & | ||
| 3407 | 52 | action='write', recl=1024) | |
| 3408 |
2/2✓ Branch 0 taken 22880 times.
✓ Branch 1 taken 52 times.
|
22932 | do i = 1, ntimstep |
| 3409 | 22880 | resolved_count = ntestpart - unresolved_orbits(i) | |
| 3410 |
2/2✓ Branch 0 taken 22879 times.
✓ Branch 1 taken 1 time.
|
22880 | if (resolved_count > 0) then |
| 3411 | 22879 | pass_fraction = confpart_pass(i)/real(resolved_count, dp) | |
| 3412 | 22879 | trap_fraction = confpart_trap(i)/real(resolved_count, dp) | |
| 3413 | else | ||
| 3414 | 1 | pass_fraction = quiet_nan() | |
| 3415 | 1 | trap_fraction = quiet_nan() | |
| 3416 | end if | ||
| 3417 | 22880 | write (unit, *) dble(kt_macro(i))*dtaumin/v0, pass_fraction, & | |
| 3418 | 45812 | trap_fraction, resolved_count | |
| 3419 | end do | ||
| 3420 | 52 | close (unit) | |
| 3421 | |||
| 3422 | open (newunit=unit, file='unresolved_fraction.dat', status='replace', & | ||
| 3423 | 52 | action='write', recl=1024) | |
| 3424 |
2/2✓ Branch 0 taken 22880 times.
✓ Branch 1 taken 52 times.
|
22932 | do i = 1, ntimstep |
| 3425 | 22880 | write (unit, *) dble(kt_macro(i))*dtaumin/v0, & | |
| 3426 | 45812 | real(unresolved_orbits(i), dp)/norm, ntestpart | |
| 3427 | end do | ||
| 3428 | 52 | close (unit) | |
| 3429 | |||
| 3430 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 43 times.
|
52 | if (classification_enabled()) then |
| 3431 | 9 | open (newunit=unit, file='class_parts.dat', recl=1024) | |
| 3432 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 9 times.
|
57 | do i = 1, ntestpart |
| 3433 | 57 | write (unit, *) i, zstart(1, i), perp_inv(i), iclass(:, i) | |
| 3434 | end do | ||
| 3435 | 9 | close (unit) | |
| 3436 | |||
| 3437 | ! Continuous fast-classifier scores and validation diagnostics, in | ||
| 3438 | ! a separate file so class_parts.dat keeps its column contract. | ||
| 3439 | open (newunit=unit, file='class_scores.dat', status='replace', & | ||
| 3440 | 9 | action='write', recl=1024) | |
| 3441 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 9 times.
|
57 | do i = 1, ntestpart |
| 3442 | 57 | write (unit, *) i, class_scores(:, i), trap_par(i) | |
| 3443 | end do | ||
| 3444 | 9 | close (unit) | |
| 3445 | |||
| 3446 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 4 times.
|
9 | if (needs_bminmax_cache()) then |
| 3447 | block | ||
| 3448 | use bminmax_mod, only: nsbmnx, hsbmnx, bmin_arr, bmax_arr | ||
| 3449 | |||
| 3450 | 5 | open (newunit=unit, file='bminmax.dat', recl=1024) | |
| 3451 |
2/2✓ Branch 0 taken 505 times.
✓ Branch 1 taken 5 times.
|
510 | do i = 0, nsbmnx |
| 3452 | 510 | write (unit, *) hsbmnx*dble(i), bmin_arr(i), bmax_arr(i) | |
| 3453 | end do | ||
| 3454 | 5 | close (unit) | |
| 3455 | end block | ||
| 3456 | end if | ||
| 3457 | end if | ||
| 3458 | 52 | end subroutine write_results | |
| 3459 | |||
| 3460 |
7/12✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 8 times.
✓ Branch 9 taken 54 times.
✓ Branch 11 taken 19 times.
✓ Branch 12 taken 34 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
123 | end module simple_main |
| 3461 |