| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module params | ||
| 2 | use, intrinsic :: iso_fortran_env, only: int8, int64 | ||
| 3 | use util, only: pi, c, e_charge, p_mass, ev | ||
| 4 | use parmot_mod, only: ro0, rmu | ||
| 5 | use new_vmec_stuff_mod, only: old_axis_healing, old_axis_healing_boundary, & | ||
| 6 | axis_healing_power_law, rho_axis_heal, & | ||
| 7 | axis_healing, s_axis_heal, & | ||
| 8 | axis_healing_polyfit_degree, & | ||
| 9 | netcdffile, ns_s, ns_tp, multharm, vmec_B_scale, & | ||
| 10 | vmec_RZ_scale | ||
| 11 | use velo_mod, only: isw_field_type | ||
| 12 | use magfie_sub, only: TEST, BOOZER | ||
| 13 | use field_can_mod, only: eval_field => evaluate, field_can_t | ||
| 14 | use field_can_base, only: set_radial_electric_potential_slope | ||
| 15 | use orbit_symplectic_base, only: symplectic_integrator_t, multistage_integrator_t, & | ||
| 16 | EXPL_IMPL_EULER, & | ||
| 17 | CASH_KARP, DORMAND_PRINCE, & | ||
| 18 | boundary_event_fraction_tolerance, & | ||
| 19 | boundary_event_radial_tolerance, & | ||
| 20 | symplectic_newton_warning_mode | ||
| 21 | use vmecin_sub, only: stevvo | ||
| 22 | use callback, only: output_error, output_orbits_macrostep | ||
| 23 | |||
| 24 | implicit none | ||
| 25 | |||
| 26 | private :: config_value_is_finite | ||
| 27 | |||
| 28 | ! Define real(dp) kind parameter | ||
| 29 | integer, parameter :: dp = kind(1.0d0) | ||
| 30 | integer :: nper = 1000, ntestpart = 1024 | ||
| 31 | integer :: zstart_dim1 = 5 | ||
| 32 | real(dp) :: dphi, phibeg = 0d0, bmod00, rlarm, bmax, bmin | ||
| 33 | real(dp) :: tau, dtau, dtaumin, xi | ||
| 34 | real(dp) :: RT0, R0i, cbfi, bz0i, bf0, rbig | ||
| 35 | real(dp), dimension(1024) :: sbeg = 0.5d0 | ||
| 36 | real(dp) :: thetabeg = 0.0d0 | ||
| 37 | real(dp), dimension(:), allocatable :: bstart, volstart | ||
| 38 | !where npoi is used, add a dimension at the end for sbeg | ||
| 39 | real(dp), dimension(:, :), allocatable :: xstart | ||
| 40 | real(dp), dimension(:, :), allocatable :: zstart, zend | ||
| 41 | real(dp), dimension(:), allocatable :: confpart_trap, confpart_pass | ||
| 42 | integer, dimension(:), allocatable :: unresolved_orbits | ||
| 43 | real(dp), dimension(:), allocatable :: times_lost | ||
| 44 | real(dp), dimension(:), allocatable :: boundary_event_radial_residual | ||
| 45 | real(dp), dimension(:), allocatable :: boundary_event_time_width | ||
| 46 | integer, dimension(:), allocatable :: orbit_exit_code | ||
| 47 | integer, parameter :: ORBIT_EXIT_COMPLETED = 0 | ||
| 48 | integer, parameter :: ORBIT_EXIT_LCFS = 1 | ||
| 49 | integer, parameter :: ORBIT_EXIT_WALL = 2 | ||
| 50 | integer, parameter :: ORBIT_EXIT_SKIPPED = 3 | ||
| 51 | ! A collisionless warning-mode marker whose validated last state is in the | ||
| 52 | ! axis-local core after all recovery methods are exhausted. It is not a | ||
| 53 | ! physical loss and is conservatively counted confined, with the distinct | ||
| 54 | ! code retained for auditability. | ||
| 55 | integer, parameter :: ORBIT_EXIT_NUMERICAL_CONFINED = 4 | ||
| 56 | integer, parameter :: ORBIT_EXIT_NUMERICAL_DOMAIN = 101 | ||
| 57 | integer, parameter :: ORBIT_EXIT_NUMERICAL_MAXITER = 102 | ||
| 58 | integer, parameter :: ORBIT_EXIT_NUMERICAL_LINEAR = 103 | ||
| 59 | integer, parameter :: ORBIT_EXIT_NUMERICAL_EVENT = 104 | ||
| 60 | integer, parameter :: ORBIT_EXIT_NUMERICAL_FULL_ORBIT = 105 | ||
| 61 | real(dp) :: contr_pp = -1d0 | ||
| 62 | integer :: ibins | ||
| 63 | logical :: generate_start_only = .False. | ||
| 64 | integer :: startmode = 1 | ||
| 65 | integer :: fusion_grid(3) = (/128, 32, 32/) | ||
| 66 | real(dp) :: grid_density = 0d0 | ||
| 67 | logical :: special_ants_file = .False. | ||
| 68 | |||
| 69 | integer :: ntau ! number of dtaumin in dtau | ||
| 70 | integer(8) :: n_microsteps_total = 0_8 | ||
| 71 | character(16) :: macrostep_time_grid = 'linear' | ||
| 72 | integer, allocatable :: ntau_macro(:) | ||
| 73 | integer(8), allocatable :: kt_macro(:) | ||
| 74 | |||
| 75 | integer :: integmode = EXPL_IMPL_EULER | ||
| 76 | |||
| 77 | ! GPU execution is configured in the same namelist as the CPU run. Native | ||
| 78 | ! CUDA with tuned DOPRI is the default production path. The legacy | ||
| 79 | ! OpenACC path remains available only when selected explicitly. | ||
| 80 | character(16) :: gpu_mode = 'off' | ||
| 81 | character(16) :: gpu_backend = 'cuda-native' | ||
| 82 | character(16) :: gpu_method = 'dopri' | ||
| 83 | character(16) :: gpu_start_coordinates = 'reference' | ||
| 84 | logical :: gpu_landreman = .true. | ||
| 85 | logical :: gpu_compact_init = .false. | ||
| 86 | logical :: gpu_work_order = .true. | ||
| 87 | integer :: gpu_num_devices = 1 | ||
| 88 | real(dp) :: gpu_t_block = 1.0d-3 | ||
| 89 | real(dp) :: gpu_loss_tau = 0.1d0 | ||
| 90 | real(dp) :: gpu_maxloss = 0.02d0 | ||
| 91 | real(dp) :: gpu_min_timestep = 1.0d-10 | ||
| 92 | real(dp) :: gpu_pilot_fraction = 0.04d0 | ||
| 93 | character(1000) :: gpu_particle_profile = '' | ||
| 94 | |||
| 95 | ! Orbit model selector. 0 = guiding-center (GC), the default symplectic | ||
| 96 | ! gyro-averaged path. 7 = full orbit (FO), the gyro-resolved Boris pusher in | ||
| 97 | ! Cartesian on the Boozer/chartmap chart, the ASCOT-style counterpart to GC. | ||
| 98 | ! Values 1-6 are reserved for other models. | ||
| 99 | integer, parameter :: ORBIT_GC = 0 | ||
| 100 | integer, parameter :: ORBIT_FULL_ORBIT = 7 | ||
| 101 | integer :: orbit_model = ORBIT_GC | ||
| 102 | |||
| 103 | ! Chart for the full-orbit field+geometry: full orbit currently supports only | ||
| 104 | ! orbit_coord = 1 (Boozer/chartmap), which shares the production GC field. | ||
| 105 | integer :: orbit_coord = 0 | ||
| 106 | |||
| 107 | integer :: kpart = 0 ! progress counter for particles | ||
| 108 | |||
| 109 | real(dp) :: relerr = 1d-13 | ||
| 110 | |||
| 111 | ! Legacy progress/checkpoint interval, retained for input compatibility. | ||
| 112 | ! Tracing performs no timed I/O; results and event totals are emitted after | ||
| 113 | ! the OpenMP team joins. | ||
| 114 | real(dp) :: checkpoint_interval = 10.0d0 | ||
| 115 | integer :: canonical_grid_nr = 62 | ||
| 116 | integer :: canonical_grid_ntheta = 63 | ||
| 117 | integer :: canonical_grid_nphi = 64 | ||
| 118 | real(dp) :: canonical_ode_relerr = 1d-11 | ||
| 119 | |||
| 120 | real(dp), allocatable :: trap_par(:), perp_inv(:) | ||
| 121 | !> Continuous classifier values, one row per particle: | ||
| 122 | !> 1 = normalized J_parallel RMS drift rate per precession turn, | ||
| 123 | !> 2 = first-half/second-half rotation-number drift, | ||
| 124 | !> 3 = achieved precession turns, 4 = J_parallel samples, | ||
| 125 | !> 5 = rotation samples per half, 6 = number of tips. | ||
| 126 | !> Columns 7:10 retain the legacy spread, reference, topology margin and | ||
| 127 | !> status for published-classifier validation; objectives use columns 1:2. | ||
| 128 | real(dp), allocatable :: class_scores(:, :) | ||
| 129 | integer, parameter :: n_class_scores = 10 | ||
| 130 | integer, allocatable :: iclass(:, :) | ||
| 131 | logical, allocatable :: class_passing(:), class_lost(:) | ||
| 132 | |||
| 133 | integer, parameter :: n_tip_vars = 6 | ||
| 134 | ! variables to evaluate at tip: z(1..5), par_inv | ||
| 135 | integer :: nplagr, nder, npl_half | ||
| 136 | integer :: norbper, nfp | ||
| 137 | real(dp) :: fper, zerolam = 0d0 | ||
| 138 | |||
| 139 | real(dp) :: tcut = -1d0 | ||
| 140 | integer(8) :: ntcut | ||
| 141 | integer :: nturns = 8 | ||
| 142 | logical :: class_plot = .False. !<=AAA | ||
| 143 | real(dp) :: cut_in_per = 0d0 !<=AAA | ||
| 144 | |||
| 145 | logical :: fast_class = .False. | ||
| 146 | !if .true. quit immeadiately after fast classification | ||
| 147 | |||
| 148 | ! Colliding with D-T reactor plasma. TODO: Make configurable | ||
| 149 | logical :: swcoll = .False. | ||
| 150 | real(dp) :: am1 = 2.0d0, am2 = 3.0d0, Z1 = 1.0d0, Z2 = 1.0d0, & | ||
| 151 | densi1 = 0.5d14, densi2 = 0.5d14, tempi1 = 1.0d4, tempi2 = 1.0d4, & | ||
| 152 | tempe = 1.0d4 | ||
| 153 | real(dp) :: dchichi, slowrate, dchichi_norm, slowrate_norm | ||
| 154 | logical :: deterministic = .False. | ||
| 155 | |||
| 156 | ! Further configuration parameters | ||
| 157 | integer :: notrace_passing = 0 | ||
| 158 | real(dp) :: facE_al = 1d0 | ||
| 159 | ! > particle energy in eV; <= 0 selects 3.5 MeV/facE_al | ||
| 160 | real(dp) :: particle_energy_eV = -1.0_dp | ||
| 161 | real(dp) :: particle_energy_eV_effective = 3.5d6 | ||
| 162 | real(dp) :: trace_time = 1d-1 | ||
| 163 | integer :: ntimstep = 10000, npoiper = 100, npoiper2 = 256, n_e = 2 | ||
| 164 | real(dp) :: n_d = 4 | ||
| 165 | |||
| 166 | real(dp) :: v0 | ||
| 167 | |||
| 168 | ! > 'full' (Coulomb operator) or 'lorentz_nustar' (monoenergetic pitch scattering) | ||
| 169 | character(24) :: collision_model = 'full' | ||
| 170 | ! > Lorentz mode: nu* = R0*nu/(abs(iota)*v0), must be >= 0 | ||
| 171 | real(dp) :: nu_star_standard = -1.0_dp | ||
| 172 | ! > Lorentz mode: local major radius R0 in cm, must be > 0 | ||
| 173 | real(dp) :: lorentz_major_radius_cm = -1.0_dp | ||
| 174 | ! > Lorentz mode: local rotational transform, must be nonzero | ||
| 175 | real(dp) :: lorentz_iota = 0.0_dp | ||
| 176 | real(dp) :: lorentz_nu = 0.0_dp | ||
| 177 | ! > slope of Phi(s) = slope*s, added to the Boozer Hamiltonian H (normalized units) | ||
| 178 | real(dp) :: radial_electric_potential_slope = 0.0_dp | ||
| 179 | |||
| 180 | logical :: debug = .False. | ||
| 181 | logical :: output_results_netcdf = .False. | ||
| 182 | logical :: restart = .False. | ||
| 183 | integer :: ierr | ||
| 184 | |||
| 185 | integer :: batch_size = 2000000000 ! Initialize large so batch mode is not default | ||
| 186 | integer :: particle_index_offset = 0 | ||
| 187 | integer :: ran_seed = 42 | ||
| 188 | integer :: num_surf = 1 | ||
| 189 | logical :: reuse_batch = .False. | ||
| 190 | integer, dimension(:), allocatable :: idx | ||
| 191 | |||
| 192 | !> Warning mode may consume one isolated, fully exhausted numerical | ||
| 193 | !> microstep. A second consecutive failure is terminal rather than freezing | ||
| 194 | !> the remainder of the orbit at one phase-space point. | ||
| 195 | integer, parameter :: max_consecutive_warning_holds = 1 | ||
| 196 | |||
| 197 | character(1000) :: field_input = '' | ||
| 198 | character(1000) :: coord_input = '' | ||
| 199 | character(1000) :: wall_input = '' | ||
| 200 | character(16) :: wall_units = 'm' | ||
| 201 | !> Radial threshold for querying an external STL wall. A negative value | ||
| 202 | !> means that no external-wall query gate is active for this run. | ||
| 203 | real(dp) :: wall_query_rho_lcfs = -1.0_dp | ||
| 204 | character(16) :: chart_boundary_kind = 'auto' | ||
| 205 | character(16) :: chart_boundary_kind_effective = 'lcfs' | ||
| 206 | integer :: integ_coords = -1000 ! Sentinel: -1000 means user did not set it | ||
| 207 | !> SPECTRE interface treatment: only full orbit (3) is supported. | ||
| 208 | integer :: crossing_level = 3 | ||
| 209 | !> Full orbit near SPECTRE interfaces: entry distance in radial gyroradii | ||
| 210 | !> and the relative tolerance of RK45 on the sharp per-volume field. | ||
| 211 | !> 1.5 keeps the switch's p_zeta invariant as at 3 (entry 1 breaks it) at | ||
| 212 | !> about half the FO time (concept/spectre_loss_study/results/fo-residence). | ||
| 213 | real(dp) :: spectre_fo_entry_larmor = 1.5_dp | ||
| 214 | real(dp) :: spectre_fo_tolerance = 1.0e-9_dp | ||
| 215 | !> Per-volume FO spline points: rho, theta, zeta per field period. | ||
| 216 | integer :: spectre_fo_grid(3) = (/32, 64, 32/) | ||
| 217 | logical :: spectre_sbeg_is_toroidal_flux = .false. | ||
| 218 | |||
| 219 | !> SPECTRE per-volume Meiss construction grid (n_r, n_th, n_phi). Each volume | ||
| 220 | !> allocates rank-3 arrays plus batch splines of this size, so the | ||
| 221 | !> total peak memory scales as Mvol*n_r*n_th*n_phi. spectre_ncon_phi = -1 means | ||
| 222 | !> automatic: 32 for fields with toroidal harmonics and a minimal phi grid | ||
| 223 | !> for axisymmetric fields (all n = 0), | ||
| 224 | !> the dominant memory saver for tokamak cases. A positive value forces that | ||
| 225 | !> phi count verbatim, bypassing the axisymmetric clamp. | ||
| 226 | integer :: spectre_ncon_r = 48 | ||
| 227 | integer :: spectre_ncon_th = 48 | ||
| 228 | integer :: spectre_ncon_phi = -1 | ||
| 229 | integer :: spectre_ncon_order = 5 | ||
| 230 | integer :: spectre_ncon_ode_max_steps = 1000000 | ||
| 231 | real(dp) :: spectre_ncon_ode_relerr = 1.0e-2_dp | ||
| 232 | |||
| 233 | ! Removed &config keys. They are still read so that old input files run; | ||
| 234 | ! a value other than the historical default prints a warning. | ||
| 235 | logical :: local = .true. | ||
| 236 | integer :: loopskip = 0 | ||
| 237 | real(dp) :: bmod_ref = 5.0d4 | ||
| 238 | real(dp) :: macrostep_log_ratio = 1.0d3 | ||
| 239 | |||
| 240 | namelist /config/ notrace_passing, nper, npoiper, ntimstep, ntestpart, & | ||
| 241 | trace_time, num_surf, sbeg, phibeg, thetabeg, contr_pp, & | ||
| 242 | facE_al, particle_energy_eV, npoiper2, n_e, n_d, netcdffile, & | ||
| 243 | ns_s, ns_tp, multharm, & | ||
| 244 | isw_field_type, generate_start_only, startmode, fusion_grid, grid_density, & | ||
| 245 | special_ants_file, integmode, orbit_model, orbit_coord, relerr, & | ||
| 246 | gpu_mode, gpu_backend, gpu_method, gpu_start_coordinates, & | ||
| 247 | gpu_landreman, gpu_compact_init, gpu_work_order, gpu_num_devices, & | ||
| 248 | gpu_t_block, gpu_loss_tau, gpu_maxloss, gpu_min_timestep, & | ||
| 249 | gpu_pilot_fraction, gpu_particle_profile, & | ||
| 250 | symplectic_newton_warning_mode, & | ||
| 251 | boundary_event_fraction_tolerance, boundary_event_radial_tolerance, & | ||
| 252 | canonical_grid_nr, canonical_grid_ntheta, canonical_grid_nphi, & | ||
| 253 | canonical_ode_relerr, & | ||
| 254 | tcut, nturns, debug, & | ||
| 255 | class_plot, cut_in_per, fast_class, vmec_B_scale, & | ||
| 256 | vmec_RZ_scale, swcoll, collision_model, nu_star_standard, & | ||
| 257 | lorentz_major_radius_cm, lorentz_iota, radial_electric_potential_slope, & | ||
| 258 | deterministic, old_axis_healing, & | ||
| 259 | old_axis_healing_boundary, axis_healing_power_law, rho_axis_heal, & | ||
| 260 | axis_healing, s_axis_heal, axis_healing_polyfit_degree, & | ||
| 261 | am1, am2, Z1, Z2, & | ||
| 262 | densi1, densi2, tempi1, tempi2, tempe, & | ||
| 263 | batch_size, ran_seed, particle_index_offset, reuse_batch, & | ||
| 264 | field_input, coord_input, & | ||
| 265 | wall_input, wall_units, chart_boundary_kind, integ_coords, crossing_level, & | ||
| 266 | spectre_fo_entry_larmor, spectre_fo_tolerance, spectre_fo_grid, & | ||
| 267 | spectre_sbeg_is_toroidal_flux, spectre_ncon_r, spectre_ncon_th, & | ||
| 268 | spectre_ncon_phi, spectre_ncon_order, spectre_ncon_ode_max_steps, & | ||
| 269 | spectre_ncon_ode_relerr, & | ||
| 270 | output_results_netcdf, & | ||
| 271 | output_error, output_orbits_macrostep, & ! callback | ||
| 272 | macrostep_time_grid, checkpoint_interval, restart, & | ||
| 273 | local, loopskip, bmod_ref, macrostep_log_ratio ! removed, ignored | ||
| 274 | |||
| 275 | integer(int8), allocatable :: wall_hit(:) | ||
| 276 | integer, allocatable :: wall_hit_triangle(:) | ||
| 277 | real(dp), allocatable :: wall_hit_cart(:, :) | ||
| 278 | real(dp), allocatable :: wall_hit_normal_cart(:, :) | ||
| 279 | real(dp), allocatable :: wall_hit_cos_incidence(:) | ||
| 280 | real(dp), allocatable :: wall_hit_angle_rad(:) | ||
| 281 | |||
| 282 | contains | ||
| 283 | |||
| 284 | 73 | subroutine read_config(config_file) | |
| 285 | character(*), intent(in) :: config_file | ||
| 286 | |||
| 287 | 73 | call read_config_namelist(config_file) | |
| 288 | 69 | call finalize_config | |
| 289 | 66 | end subroutine read_config | |
| 290 | |||
| 291 | 74 | subroutine read_config_namelist(config_file) | |
| 292 | ! Read the &config group without resolving aliases or validating. | ||
| 293 | character(*), intent(in) :: config_file | ||
| 294 | integer :: unit | ||
| 295 | |||
| 296 | 74 | open (newunit=unit, file=config_file, status='old', action='read') | |
| 297 | 74 | read (unit, nml=config) | |
| 298 | 70 | close (unit) | |
| 299 | 70 | end subroutine read_config_namelist | |
| 300 | |||
| 301 | 12 | subroutine reset_config | |
| 302 | ! Restore every &config value to its compiled default. The first call | ||
| 303 | ! records the defaults, so call it before changing any value. | ||
| 304 | integer, save :: defaults_unit | ||
| 305 | logical, save :: recorded = .false. | ||
| 306 | |||
| 307 |
2/2✓ Branch 0 taken 11 times.
✓ Branch 1 taken 1 time.
|
12 | if (.not. recorded) then |
| 308 | 11 | open (newunit=defaults_unit, status='scratch', delim='apostrophe') | |
| 309 | 11 | write (defaults_unit, nml=config) | |
| 310 | 11 | recorded = .true. | |
| 311 | end if | ||
| 312 | 12 | rewind (defaults_unit) | |
| 313 | 12 | read (defaults_unit, nml=config) | |
| 314 | 12 | end subroutine reset_config | |
| 315 | |||
| 316 | 81 | subroutine finalize_config | |
| 317 | ! Resolve aliases, normalize and validate the current &config values. | ||
| 318 | 81 | call warn_removed_config_keys | |
| 319 | 81 | call apply_config_aliases | |
| 320 | |||
| 321 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | gpu_mode = to_lower(trim(gpu_mode)) |
| 322 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | gpu_backend = to_lower(trim(gpu_backend)) |
| 323 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | gpu_method = to_lower(trim(gpu_method)) |
| 324 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | gpu_start_coordinates = to_lower(trim(gpu_start_coordinates)) |
| 325 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | collision_model = to_lower(trim(collision_model)) |
| 326 | |||
| 327 | 81 | call reset_seed_if_deterministic | |
| 328 | |||
| 329 |
3/4✓ Branch 0 taken 59 times.
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 59 times.
✗ Branch 3 not taken.
|
81 | if (integmode > 0 .and. symplectic_newton_warning_mode) then |
| 330 | 59 | print *, 'WARNING: symplectic integrators accept bounded Newton ', & | |
| 331 | 59 | 'corrections after maxit; failed steps retry from the last ', & | |
| 332 | 118 | 'valid state; see maxit/step-skip diagnostics' | |
| 333 | end if | ||
| 334 | |||
| 335 | 81 | call validate_boundary_event_tolerances | |
| 336 | 81 | call validate_chart_boundary_kind | |
| 337 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
|
81 | if (min(canonical_grid_nr, canonical_grid_ntheta, & |
| 338 | canonical_grid_nphi) < 6) then | ||
| 339 | ✗ | error stop 'canonical map grid dimensions must be at least 6' | |
| 340 | end if | ||
| 341 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
|
81 | if (max(canonical_grid_nr, canonical_grid_ntheta, & |
| 342 | canonical_grid_nphi) > 512) then | ||
| 343 | ✗ | error stop 'canonical map grid dimensions must not exceed 512' | |
| 344 | end if | ||
| 345 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
|
81 | if (int(canonical_grid_nr, int64)*int(canonical_grid_ntheta, int64)* & |
| 346 | int(canonical_grid_nphi, int64) > 2097152_int64) then | ||
| 347 | ✗ | error stop 'canonical map grid exceeds the 2097152-point limit' | |
| 348 | end if | ||
| 349 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
|
81 | if (.not. config_value_is_finite(canonical_ode_relerr) .or. & |
| 350 | canonical_ode_relerr <= 0d0) then | ||
| 351 | ✗ | error stop 'canonical_ode_relerr must be finite and positive' | |
| 352 | end if | ||
| 353 | |||
| 354 |
1/8✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
81 | if (swcoll .and. (tcut > 0.0d0 .or. class_plot .or. fast_class)) then |
| 355 | ✗ | error stop 'Collisions are incompatible with classification' | |
| 356 | end if | ||
| 357 | |||
| 358 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 78 times.
|
81 | if (crossing_level /= 3) then |
| 359 | 3 | error stop 'crossing_level must be 3; fixed-mu GC maps were removed' | |
| 360 | end if | ||
| 361 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | if (spectre_fo_entry_larmor <= 0.0_dp) & |
| 362 | ✗ | error stop 'spectre_fo_entry_larmor must be positive' | |
| 363 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | if (spectre_fo_tolerance <= 0.0_dp) & |
| 364 | ✗ | error stop 'spectre_fo_tolerance must be positive' | |
| 365 | 78 | call validate_collision_config | |
| 366 | 78 | call validate_radial_electric_potential | |
| 367 | 78 | call set_radial_electric_potential_slope(radial_electric_potential_slope) | |
| 368 | |||
| 369 | 78 | call validate_gpu_config | |
| 370 | 78 | end subroutine finalize_config | |
| 371 | |||
| 372 | 81 | subroutine warn_removed_config_keys | |
| 373 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 80 times.
|
82 | if (.not. local) print *, 'WARNING: removed &config key local is ignored; ', & |
| 374 | 2 | 'use num_surf=0 for volume sampling' | |
| 375 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 80 times.
|
81 | if (loopskip /= 0) print *, 'WARNING: removed &config key loopskip is ignored' |
| 376 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 80 times.
|
82 | if (bmod_ref /= 5.0d4) print *, 'WARNING: removed &config key bmod_ref ', & |
| 377 | 2 | 'is ignored; use vmec_B_scale to scale the field' | |
| 378 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 80 times.
|
82 | if (macrostep_log_ratio /= 1.0d3) print *, 'WARNING: removed &config key ', & |
| 379 | 2 | 'macrostep_log_ratio is ignored' | |
| 380 | 81 | end subroutine warn_removed_config_keys | |
| 381 | |||
| 382 | 84 | subroutine validate_radial_electric_potential | |
| 383 | ! Only field_can_boozer sets Phie; every other field, the RK path, full | ||
| 384 | ! orbits and the GPU kernels would silently ignore the potential. | ||
| 385 | 84 | if (.not. config_value_is_finite(radial_electric_potential_slope)) & | |
| 386 | ✗ | error stop 'radial_electric_potential_slope must be finite' | |
| 387 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 79 times.
|
84 | if (radial_electric_potential_slope == 0.0_dp) return |
| 388 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 4 times.
|
5 | if (isw_field_type /= BOOZER) error stop & |
| 389 | 1 | 'radial_electric_potential_slope requires Boozer integ_coords' | |
| 390 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
|
4 | if (integmode <= 0) error stop & |
| 391 | 1 | 'radial_electric_potential_slope requires a symplectic integmode' | |
| 392 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2 times.
|
3 | if (orbit_model /= ORBIT_GC) error stop & |
| 393 | 1 | 'radial_electric_potential_slope requires the guiding-center model' | |
| 394 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | if (trim(gpu_mode) /= 'off') error stop & |
| 395 | 1 | 'radial_electric_potential_slope is not supported with gpu_mode' | |
| 396 | end subroutine validate_radial_electric_potential | ||
| 397 | |||
| 398 | 78 | subroutine validate_collision_config | |
| 399 | real(dp) :: energy | ||
| 400 | integer :: energy_status | ||
| 401 | |||
| 402 | 78 | call resolve_particle_energy(energy, energy_status) | |
| 403 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | if (energy_status /= 0) & |
| 404 | ✗ | error stop 'particle_energy_eV and facE_al specify inconsistent energies' | |
| 405 | |||
| 406 | 78 | select case (trim(collision_model)) | |
| 407 | case ('full') | ||
| 408 | case ('lorentz_nustar') | ||
| 409 | ✗ | if (nu_star_standard < 0.0_dp) & | |
| 410 | ✗ | error stop 'nu_star_standard must be nonnegative in Lorentz mode' | |
| 411 | ✗ | if (lorentz_major_radius_cm <= 0.0_dp) & | |
| 412 | ✗ | error stop 'lorentz_major_radius_cm must be positive in Lorentz mode' | |
| 413 | ✗ | if (abs(lorentz_iota) <= tiny(1.0_dp)) & | |
| 414 | ✗ | error stop 'lorentz_iota must be nonzero in Lorentz mode' | |
| 415 | case default | ||
| 416 |
1/3✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 78 times.
|
78 | error stop 'collision_model must be full or lorentz_nustar' |
| 417 | end select | ||
| 418 | 78 | end subroutine validate_collision_config | |
| 419 | |||
| 420 | 149 | subroutine resolve_particle_energy(energy_eV, status) | |
| 421 | real(dp), intent(out) :: energy_eV | ||
| 422 | integer, intent(out) :: status | ||
| 423 | |||
| 424 | real(dp) :: legacy_energy | ||
| 425 | |||
| 426 | 149 | status = 0 | |
| 427 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 149 times.
|
149 | if (facE_al <= 0.0_dp) then |
| 428 | ✗ | status = 1 | |
| 429 | ✗ | energy_eV = 0.0_dp | |
| 430 | ✗ | return | |
| 431 | end if | ||
| 432 | 149 | legacy_energy = 3.5d6/facE_al | |
| 433 |
2/2✓ Branch 0 taken 147 times.
✓ Branch 1 taken 2 times.
|
149 | if (particle_energy_eV <= 0.0_dp) then |
| 434 | 147 | energy_eV = legacy_energy | |
| 435 | 147 | return | |
| 436 | end if | ||
| 437 | |||
| 438 | 2 | energy_eV = particle_energy_eV | |
| 439 | ! facE_al=1 is the historical default and is treated as unspecified when | ||
| 440 | ! the explicit input is present. Any nondefault value must agree. | ||
| 441 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | if (abs(facE_al - 1.0_dp) > epsilon(1.0_dp)) then |
| 442 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (abs(legacy_energy - energy_eV) > & |
| 443 | 1 | 1.0e-12_dp*max(legacy_energy, energy_eV)) status = 2 | |
| 444 | end if | ||
| 445 | end subroutine resolve_particle_energy | ||
| 446 | |||
| 447 | 78 | subroutine validate_gpu_config | |
| 448 | 78 | select case (trim(gpu_mode)) | |
| 449 | case ('off', 'production', 'benchmark') | ||
| 450 | case default | ||
| 451 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | error stop 'gpu_mode must be off, production, or benchmark' |
| 452 | end select | ||
| 453 | |||
| 454 | 78 | select case (trim(gpu_backend)) | |
| 455 | case ('openacc', 'cuda-native', 'auto') | ||
| 456 | case default | ||
| 457 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | error stop 'gpu_backend must be openacc, cuda-native, or auto' |
| 458 | end select | ||
| 459 | |||
| 460 | 78 | select case (trim(gpu_method)) | |
| 461 | case ('auto', 'euler', 'midpoint', 'cash-karp', 'dopri') | ||
| 462 | case default | ||
| 463 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | error stop 'gpu_method must be auto, euler, midpoint, cash-karp, or dopri' |
| 464 | end select | ||
| 465 | |||
| 466 | 78 | select case (trim(gpu_start_coordinates)) | |
| 467 | case ('reference', 'boozer') | ||
| 468 | case default | ||
| 469 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | error stop 'gpu_start_coordinates must be reference or boozer' |
| 470 | end select | ||
| 471 | |||
| 472 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | if (gpu_num_devices < 1) error stop 'gpu_num_devices must be positive' |
| 473 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | if (gpu_t_block <= 0.0d0) error stop 'gpu_t_block must be positive' |
| 474 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | if (gpu_loss_tau <= 0.0d0) error stop 'gpu_loss_tau must be positive' |
| 475 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | if (gpu_maxloss < 0.0d0) error stop 'gpu_maxloss must be nonnegative' |
| 476 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | if (gpu_min_timestep < 0.0d0) & |
| 477 | ✗ | error stop 'gpu_min_timestep must be nonnegative' | |
| 478 |
2/4✓ Branch 0 taken 78 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 78 times.
|
78 | if (gpu_pilot_fraction < 0.0d0 .or. gpu_pilot_fraction > 0.25d0) & |
| 479 | ✗ | error stop 'gpu_pilot_fraction must be in [0, 0.25]' | |
| 480 | 78 | end subroutine validate_gpu_config | |
| 481 | |||
| 482 | 83 | subroutine validate_chart_boundary_kind | |
| 483 | 84 | select case (trim(chart_boundary_kind)) | |
| 484 | case ('auto', 'lcfs', 'wall', 'domain') | ||
| 485 | case default | ||
| 486 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 82 times.
|
83 | error stop "chart_boundary_kind must be auto, lcfs, wall, or domain" |
| 487 | end select | ||
| 488 | 82 | end subroutine validate_chart_boundary_kind | |
| 489 | |||
| 490 | 85 | subroutine validate_boundary_event_tolerances | |
| 491 |
3/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
83 | if (.not. config_value_is_finite(boundary_event_fraction_tolerance) .or. & |
| 492 | (boundary_event_fraction_tolerance /= -1d0 .and. & | ||
| 493 | boundary_event_fraction_tolerance <= 0d0)) then | ||
| 494 | 2 | error stop 'boundary_event_fraction_tolerance must be finite and positive or -1' | |
| 495 | end if | ||
| 496 |
3/4✓ Branch 0 taken 1 time.
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
|
82 | if (.not. config_value_is_finite(boundary_event_radial_tolerance) .or. & |
| 497 | (boundary_event_radial_tolerance /= -1d0 .and. & | ||
| 498 | boundary_event_radial_tolerance <= 0d0)) then | ||
| 499 | 1 | error stop 'boundary_event_radial_tolerance must be finite and positive or -1' | |
| 500 | end if | ||
| 501 | 82 | end subroutine validate_boundary_event_tolerances | |
| 502 | |||
| 503 | 333 | pure logical function config_value_is_finite(value) | |
| 504 | real(dp), intent(in) :: value | ||
| 505 | integer(int64), parameter :: exponent_mask = & | ||
| 506 | int(z'7ff0000000000000', int64) | ||
| 507 | integer(int64) :: bits | ||
| 508 | |||
| 509 | 333 | bits = transfer(value, bits) | |
| 510 |
6/8✓ Branch 0 taken 81 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 84 times.
✓ Branch 4 taken 83 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 82 times.
✓ Branch 7 taken 1 time.
|
333 | config_value_is_finite = iand(bits, exponent_mask) /= exponent_mask |
| 511 | end function config_value_is_finite | ||
| 512 | |||
| 513 | 68 | subroutine params_init | |
| 514 | real(dp) :: E_alpha | ||
| 515 | integer :: energy_status | ||
| 516 | integer :: L1i | ||
| 517 | real(dp) :: weight_sum, cumul_weight, w | ||
| 518 | integer :: i, nintv | ||
| 519 | integer(8) :: kt_target, kt_prev | ||
| 520 | character(16) :: grid_kind | ||
| 521 | 68 | call resolve_particle_energy(E_alpha, energy_status) | |
| 522 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
|
68 | if (energy_status /= 0) error stop 'Invalid particle energy configuration' |
| 523 | 68 | particle_energy_eV_effective = E_alpha | |
| 524 | |||
| 525 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 65 times.
|
68 | if (isw_field_type == TEST) then |
| 526 | ! TEST field uses normalized units: B0=1, R0=1, a=0.5 | ||
| 527 | ! Use normalized parameters for orbit integration | ||
| 528 | 3 | v0 = 1.0d0 ! Normalized velocity | |
| 529 | 3 | rlarm = 1.0d0 ! Normalized Larmor radius (will be scaled by ro0) | |
| 530 | 3 | ro0 = 1.0d0 ! Normalized ro0 = mc/e * v0 | |
| 531 | 3 | rmu = 1d8 ! Large inverse relativistic temperature | |
| 532 | 3 | tau = trace_time ! Already in normalized units | |
| 533 | 3 | dtau = tau/dble(ntimstep - 1) | |
| 534 | 3 | L1i = 1 ! One field period (full torus) | |
| 535 | 3 | rbig = 1.0d0 ! Major radius R0=1 | |
| 536 | 3 | dphi = 2.d0*pi/(L1i*npoiper) | |
| 537 | 3 | dtaumin = 2.d0*pi*rbig/npoiper2 | |
| 538 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | ntau = ceiling(dtau/dtaumin) |
| 539 | 3 | dtaumin = dtau/ntau | |
| 540 | 3 | fper = 2d0*pi ! Full torus | |
| 541 | else | ||
| 542 | ! set alpha energy, velocity, and Larmor radius | ||
| 543 | 65 | v0 = sqrt(2.d0*E_alpha*ev/(n_d*p_mass)) | |
| 544 | 65 | rlarm = v0*n_d*p_mass*c/(n_e*e_charge) | |
| 545 | 65 | ro0 = rlarm | |
| 546 | ! Neglect relativistic effects by large inverse relativistic temperature | ||
| 547 | 65 | rmu = 1d8 | |
| 548 | ! normalized slowing down time: | ||
| 549 | 65 | tau = trace_time*v0 | |
| 550 | ! normalized time step: | ||
| 551 | 65 | dtau = tau/dble(ntimstep - 1) | |
| 552 | ! parameters for the vacuum chamber: | ||
| 553 | 65 | call stevvo(RT0, R0i, L1i, cbfi, bz0i, bf0) | |
| 554 | 65 | rbig = rt0 | |
| 555 | ! field line integration step step over phi (to check chamber wall crossing) | ||
| 556 | 65 | dphi = 2.d0*pi/(L1i*npoiper) | |
| 557 | ! orbit integration time step (to check chamber wall crossing) | ||
| 558 | 65 | dtaumin = 2.d0*pi*rbig/npoiper2 | |
| 559 |
1/2✓ Branch 0 taken 65 times.
✗ Branch 1 not taken.
|
65 | ntau = ceiling(dtau/dtaumin) |
| 560 | 65 | dtaumin = dtau/ntau | |
| 561 | 65 | fper = 2d0*pi/dble(L1i) | |
| 562 | end if | ||
| 563 | |||
| 564 | ! Macrostep schedule (number of microsteps per macrostep). | ||
| 565 | ! - Default is linear: constant ntau per macrostep. | ||
| 566 | ! - Log schedule: macrosteps are distributed logarithmically in time | ||
| 567 | ! while keeping the microstep resolution dtaumin. The total number | ||
| 568 | ! of microsteps is preserved as (ntimstep-1)*ntau. | ||
| 569 | 68 | nintv = max(1, ntimstep - 1) | |
| 570 | 68 | n_microsteps_total = int(nintv, kind=8) * int(ntau, kind=8) | |
| 571 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
|
68 | if (allocated(ntau_macro)) deallocate(ntau_macro) |
| 572 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
|
68 | if (allocated(kt_macro)) deallocate(kt_macro) |
| 573 |
7/14✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
✓ Branch 2 taken 68 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 68 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 68 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 68 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 68 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 68 times.
|
204 | allocate (ntau_macro(ntimstep)) |
| 574 |
7/14✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
✓ Branch 2 taken 68 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 68 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 68 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 68 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 68 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 68 times.
|
204 | allocate (kt_macro(ntimstep)) |
| 575 |
2/2✓ Branch 0 taken 33471 times.
✓ Branch 1 taken 68 times.
|
33539 | ntau_macro = 0 |
| 576 |
2/2✓ Branch 0 taken 33471 times.
✓ Branch 1 taken 68 times.
|
33539 | kt_macro = 0_8 |
| 577 | |||
| 578 | 68 | grid_kind = to_lower(macrostep_time_grid) | |
| 579 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
|
68 | if (trim(grid_kind) == 'log') then |
| 580 | ! Logarithmic macrostep grid: timesteps appear equally spaced | ||
| 581 | ! on a log10 plot of time. Weight for interval i is 10^((i-1)/(N-1)), | ||
| 582 | ! giving a 10:1 ratio between last and first intervals. | ||
| 583 | ! Uses cumulative rounding to guarantee exact total microsteps. | ||
| 584 | |||
| 585 | ! First pass: compute total weight | ||
| 586 | weight_sum = 0.0d0 | ||
| 587 | ✗ | do i = 1, nintv | |
| 588 | ✗ | weight_sum = weight_sum + 10.0d0 ** (dble(i - 1) / dble(max(1, nintv - 1))) | |
| 589 | end do | ||
| 590 | |||
| 591 | ! Second pass: distribute microsteps via cumulative rounding | ||
| 592 | ! Guarantee at least 1 microstep per interval to avoid empty loops | ||
| 593 | cumul_weight = 0.0d0 | ||
| 594 | kt_prev = 0_8 | ||
| 595 | ✗ | do i = 1, nintv | |
| 596 | ✗ | w = 10.0d0 ** (dble(i - 1) / dble(max(1, nintv - 1))) | |
| 597 | ✗ | cumul_weight = cumul_weight + w | |
| 598 | ✗ | kt_target = nint(cumul_weight / weight_sum * dble(n_microsteps_total), kind=8) | |
| 599 | ✗ | ntau_macro(i + 1) = max(1, int(kt_target - kt_prev)) | |
| 600 | ✗ | kt_macro(i + 1) = kt_prev + int(ntau_macro(i + 1), kind=8) | |
| 601 | ✗ | kt_prev = kt_macro(i + 1) | |
| 602 | end do | ||
| 603 | else | ||
| 604 |
2/2✓ Branch 0 taken 68 times.
✓ Branch 1 taken 33403 times.
|
33471 | do i = 2, ntimstep |
| 605 | 33403 | ntau_macro(i) = ntau | |
| 606 | 33471 | kt_macro(i) = kt_macro(i - 1) + int(ntau, kind=8) | |
| 607 | end do | ||
| 608 | end if | ||
| 609 | |||
| 610 | 68 | ntcut = microstep_cut_index(ntimstep, ntau, tcut, trace_time) | |
| 611 |
1/2✓ Branch 0 taken 68 times.
✗ Branch 1 not taken.
|
68 | norbper = ceiling(1d0*ntau*ntimstep/(L1i*npoiper2)) |
| 612 | 68 | nfp = L1i*norbper | |
| 613 | |||
| 614 | 68 | zerolam = 0.d0 | |
| 615 | 68 | nplagr = 4 | |
| 616 | 68 | nder = 0 | |
| 617 | 68 | npl_half = nplagr/2 | |
| 618 | |||
| 619 | 68 | call init_batch | |
| 620 | 68 | call reallocate_arrays | |
| 621 | 68 | end subroutine params_init | |
| 622 | |||
| 623 | 72 | pure function microstep_cut_index(ntimstep_in, ntau_in, tcut_in, & | |
| 624 | trace_time_in) result(ntcut_out) | ||
| 625 | ! Microstep index of the classification cut time tcut; <=0 disables it. | ||
| 626 | ! ntimstep*ntau reaches ~1e10 for second-scale traces and overflows a | ||
| 627 | ! 32-bit product, flipping the sign and spuriously enabling classification, | ||
| 628 | ! so evaluate in real(dp) and return int64. | ||
| 629 | integer, intent(in) :: ntimstep_in, ntau_in | ||
| 630 | real(dp), intent(in) :: tcut_in, trace_time_in | ||
| 631 | integer(8) :: ntcut_out | ||
| 632 | |||
| 633 |
3/4✓ Branch 0 taken 3 times.
✓ Branch 1 taken 69 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
72 | if (tcut_in > 0d0 .and. trace_time_in > 0d0) then |
| 634 | ntcut_out = ceiling( & | ||
| 635 | real(ntimstep_in, dp)*real(ntau_in, dp)*tcut_in/trace_time_in, & | ||
| 636 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | kind=8) |
| 637 | else | ||
| 638 | ntcut_out = -1_8 | ||
| 639 | end if | ||
| 640 | 72 | end function microstep_cut_index | |
| 641 | |||
| 642 | 473 | pure function to_lower(s) result(out) | |
| 643 | character(*), intent(in) :: s | ||
| 644 | character(len(s)) :: out | ||
| 645 | integer :: i, c | ||
| 646 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 473 times.
|
473 | out = s |
| 647 |
2/2✓ Branch 0 taken 3699 times.
✓ Branch 1 taken 473 times.
|
4172 | do i = 1, len(s) |
| 648 | 3699 | c = iachar(out(i:i)) | |
| 649 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3699 times.
|
4172 | if (c >= iachar('A') .and. c <= iachar('Z')) then |
| 650 | ✗ | out(i:i) = achar(c + 32) | |
| 651 | end if | ||
| 652 | end do | ||
| 653 |
1/2✓ Branch 0 taken 473 times.
✗ Branch 1 not taken.
|
473 | end function to_lower |
| 654 | |||
| 655 | 71 | subroutine init_batch | |
| 656 | integer :: unit, status, source_index, selected, source_count, seedsize | ||
| 657 | 71 | integer, allocatable :: batch_seed(:) | |
| 658 | real(dp) :: draw | ||
| 659 | logical :: old_batch | ||
| 660 | |||
| 661 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 69 times.
|
71 | if (allocated(idx)) deallocate (idx) |
| 662 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 71 times.
|
71 | if (particle_index_offset < 0) & |
| 663 | ✗ | error stop 'particle_index_offset must be nonnegative' | |
| 664 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 71 times.
|
71 | if (batch_size < 1) error stop 'batch_size must be positive' |
| 665 |
2/2✓ Branch 0 taken 68 times.
✓ Branch 1 taken 3 times.
|
71 | if (ntestpart <= batch_size) return |
| 666 | 3 | source_count = ntestpart | |
| 667 |
6/12✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 3 times.
|
9 | allocate (idx(batch_size)) |
| 668 | 3 | inquire (file='batch.dat', exist=old_batch) | |
| 669 |
3/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
3 | if (reuse_batch .and. old_batch) then |
| 670 | 2 | open (newunit=unit, file='batch.dat', status='old', action='read') | |
| 671 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 2 times.
|
12 | do selected = 1, batch_size |
| 672 | 10 | read (unit, *, iostat=status) idx(selected) | |
| 673 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
10 | if (status /= 0) error stop 'Invalid batch.dat index list' |
| 674 |
2/4✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
|
10 | if (idx(selected) < 1 .or. idx(selected) > source_count) & |
| 675 | ✗ | error stop 'Batch index outside start.dat source range' | |
| 676 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
|
12 | if (selected > 1) then |
| 677 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | if (idx(selected) <= idx(selected - 1)) & |
| 678 | ✗ | error stop 'Batch indices must be strictly increasing' | |
| 679 | end if | ||
| 680 | end do | ||
| 681 | 2 | close (unit) | |
| 682 | else | ||
| 683 | 1 | reuse_batch = .false. | |
| 684 | 1 | call random_seed(size=seedsize) | |
| 685 |
7/14✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✓ 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 not taken.
✓ Branch 9 taken 1 time.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 time.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 time.
|
3 | allocate (batch_seed(seedsize)) |
| 686 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 time.
|
9 | batch_seed = ran_seed |
| 687 | 1 | call random_seed(put=batch_seed) | |
| 688 | 1 | selected = 0 | |
| 689 | ! Sequential sampling without replacement gives sorted, one-based | ||
| 690 | ! source rows without allocating the complete source population. | ||
| 691 |
1/2✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
|
10 | do source_index = 1, source_count |
| 692 | 10 | call random_number(draw) | |
| 693 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
|
10 | if (draw*real(source_count - source_index + 1, dp) < & |
| 694 | 1 | real(batch_size - selected, dp)) then | |
| 695 | 5 | selected = selected + 1 | |
| 696 | 5 | idx(selected) = source_index | |
| 697 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 4 times.
|
5 | if (selected == batch_size) exit |
| 698 | end if | ||
| 699 | end do | ||
| 700 | 1 | open (newunit=unit, file='batch.dat', status='replace', action='write') | |
| 701 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1 time.
|
6 | do selected = 1, batch_size |
| 702 | 6 | write (unit, *) idx(selected) | |
| 703 | end do | ||
| 704 | 1 | close (unit) | |
| 705 | end if | ||
| 706 | 3 | ntestpart = batch_size | |
| 707 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2 times.
|
3 | end subroutine init_batch |
| 708 | |||
| 709 | 50 | integer function particle_stream_id(local_index) result(global_index) | |
| 710 | integer, intent(in) :: local_index | ||
| 711 | integer :: source_index | ||
| 712 | |||
| 713 | 50 | source_index = local_index | |
| 714 |
3/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
|
50 | if (startmode == 4 .and. allocated(idx)) then |
| 715 |
2/4✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
|
10 | if (local_index < 1 .or. local_index > size(idx)) & |
| 716 | ✗ | error stop 'Particle index outside selected batch' | |
| 717 | 10 | source_index = idx(local_index) | |
| 718 | end if | ||
| 719 |
2/4✓ Branch 0 taken 50 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 50 times.
|
50 | if (source_index < 1 .or. particle_index_offset < 0) & |
| 720 | ✗ | error stop 'Particle stream identity must be positive' | |
| 721 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 50 times.
|
50 | if (int(particle_index_offset, int64) > & |
| 722 | 2147483646_int64 - int(source_index, int64)) & | ||
| 723 | ✗ | error stop 'Particle stream identity exceeds integer range' | |
| 724 | 50 | global_index = particle_index_offset + source_index | |
| 725 | 50 | end function particle_stream_id | |
| 726 | |||
| 727 | 81 | subroutine reallocate_arrays | |
| 728 | integer :: npoi | ||
| 729 | 81 | npoi = nper*npoiper ! total number of starting points | |
| 730 | |||
| 731 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(zstart)) deallocate (zstart) |
| 732 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(zend)) deallocate (zend) |
| 733 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(times_lost)) deallocate (times_lost) |
| 734 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(boundary_event_radial_residual)) & |
| 735 | 13 | deallocate (boundary_event_radial_residual) | |
| 736 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(boundary_event_time_width)) deallocate (boundary_event_time_width) |
| 737 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(orbit_exit_code)) deallocate (orbit_exit_code) |
| 738 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(trap_par)) deallocate (trap_par) |
| 739 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(perp_inv)) deallocate (perp_inv) |
| 740 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(iclass)) deallocate (iclass) |
| 741 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(class_scores)) deallocate (class_scores) |
| 742 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(class_passing)) deallocate (class_passing) |
| 743 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(class_lost)) deallocate (class_lost) |
| 744 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(xstart)) deallocate (xstart) |
| 745 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(bstart)) deallocate (bstart) |
| 746 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(volstart)) deallocate (volstart) |
| 747 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(confpart_trap)) deallocate (confpart_trap) |
| 748 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(confpart_pass)) deallocate (confpart_pass) |
| 749 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(unresolved_orbits)) deallocate (unresolved_orbits) |
| 750 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(wall_hit_triangle)) deallocate (wall_hit_triangle) |
| 751 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(wall_hit)) deallocate (wall_hit) |
| 752 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(wall_hit_cart)) deallocate (wall_hit_cart) |
| 753 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(wall_hit_normal_cart)) deallocate (wall_hit_normal_cart) |
| 754 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(wall_hit_cos_incidence)) deallocate (wall_hit_cos_incidence) |
| 755 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 68 times.
|
81 | if (allocated(wall_hit_angle_rad)) deallocate (wall_hit_angle_rad) |
| 756 | |||
| 757 |
16/32✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 81 times.
✓ Branch 6 taken 81 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 81 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 81 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 81 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 81 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 81 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 81 times.
✓ Branch 23 taken 81 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 81 times.
✓ Branch 27 taken 81 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 30 taken 81 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 81 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 81 times.
|
486 | allocate (zstart(zstart_dim1, ntestpart), zend(zstart_dim1, ntestpart)) |
| 758 |
17/34✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 81 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 81 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 81 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 81 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 81 times.
✓ Branch 19 taken 81 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 81 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 81 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 81 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 81 times.
✓ Branch 32 taken 81 times.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✓ Branch 35 taken 81 times.
✗ Branch 37 not taken.
✓ Branch 38 taken 81 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 81 times.
|
405 | allocate (times_lost(ntestpart), trap_par(ntestpart), perp_inv(ntestpart)) |
| 759 | ✗ | allocate (boundary_event_radial_residual(ntestpart), & | |
| 760 |
10/20✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 81 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 81 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 81 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 81 times.
✓ Branch 14 taken 81 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 81 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 81 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 81 times.
|
243 | boundary_event_time_width(ntestpart)) |
| 761 |
7/14✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 81 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 81 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 81 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 81 times.
|
243 | allocate (orbit_exit_code(ntestpart)) |
| 762 |
19/38✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 81 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 81 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 81 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 81 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 81 times.
✓ Branch 19 taken 81 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 81 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 81 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 81 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 81 times.
✗ Branch 31 not taken.
✓ Branch 32 taken 81 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 81 times.
✓ Branch 36 taken 81 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 81 times.
✗ Branch 41 not taken.
✓ Branch 42 taken 81 times.
✗ Branch 44 not taken.
✓ Branch 45 taken 81 times.
|
486 | allocate (xstart(3, npoi), bstart(npoi), volstart(npoi)) |
| 763 |
12/24✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 81 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 81 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 81 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 81 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 81 times.
✓ Branch 19 taken 81 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 81 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 81 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 81 times.
|
324 | allocate (confpart_trap(ntimstep), confpart_pass(ntimstep)) |
| 764 |
7/14✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 81 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 81 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 81 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 81 times.
|
243 | allocate (unresolved_orbits(ntimstep)) |
| 765 |
7/14✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 81 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 81 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 81 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 81 times.
|
243 | allocate (iclass(3, ntestpart)) |
| 766 |
7/14✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 81 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 81 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 81 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 81 times.
|
243 | allocate (class_scores(n_class_scores, ntestpart)) |
| 767 |
10/20✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 81 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 81 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 81 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 81 times.
✓ Branch 15 taken 81 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 81 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 81 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 81 times.
|
243 | allocate (class_passing(ntestpart), class_lost(ntestpart)) |
| 768 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 81 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 81 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 81 times.
|
162 | allocate (wall_hit_triangle(ntestpart)) |
| 769 |
6/12✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 81 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 81 times.
|
243 | allocate (wall_hit(ntestpart)) |
| 770 |
7/14✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 81 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 81 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 81 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 81 times.
|
243 | allocate (wall_hit_cart(3, ntestpart)) |
| 771 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 81 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 81 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 81 times.
|
162 | allocate (wall_hit_normal_cart(3, ntestpart)) |
| 772 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 81 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 81 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 81 times.
|
162 | allocate (wall_hit_cos_incidence(ntestpart)) |
| 773 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 81 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 81 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 81 times.
|
162 | allocate (wall_hit_angle_rad(ntestpart)) |
| 774 | |||
| 775 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
|
811200 | times_lost = 0.0d0 |
| 776 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
|
811200 | boundary_event_radial_residual = -1d0 |
| 777 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
|
811200 | boundary_event_time_width = -1d0 |
| 778 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
|
811200 | orbit_exit_code = 0 |
| 779 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
|
811200 | trap_par = 0.0d0 |
| 780 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
|
811200 | perp_inv = 0.0d0 |
| 781 |
4/4✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
✓ Branch 2 taken 2433357 times.
✓ Branch 3 taken 811119 times.
|
3244557 | iclass = 0 |
| 782 |
4/4✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
✓ Branch 2 taken 8111190 times.
✓ Branch 3 taken 811119 times.
|
8922390 | class_scores = 0d0 |
| 783 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
|
811200 | class_passing = .false. |
| 784 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
|
811200 | class_lost = .false. |
| 785 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
|
811200 | wall_hit_triangle = -1 |
| 786 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
|
811200 | wall_hit = 0_int8 |
| 787 |
4/4✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
✓ Branch 2 taken 2433357 times.
✓ Branch 3 taken 811119 times.
|
3244557 | wall_hit_cart = 0.0d0 |
| 788 |
4/4✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
✓ Branch 2 taken 2433357 times.
✓ Branch 3 taken 811119 times.
|
3244557 | wall_hit_normal_cart = 0.0d0 |
| 789 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 811119 times.
|
811200 | wall_hit_cos_incidence = 0.0d0 |
| 790 |
2/2✓ Branch 0 taken 811119 times.
✓ Branch 1 taken 81 times.
|
811200 | wall_hit_angle_rad = 0.0d0 |
| 791 | 81 | end subroutine reallocate_arrays | |
| 792 | |||
| 793 | ✗ | subroutine sort_idx(idx_arr, N) | |
| 794 | ! sort particle indices. | ||
| 795 | integer :: N | ||
| 796 | integer, dimension(N) :: idx_arr | ||
| 797 | integer :: i, j, temp, o, r, num_removed, p | ||
| 798 | logical :: swapped | ||
| 799 | |||
| 800 | ✗ | do j = n - 1, 1, -1 | |
| 801 | swapped = .false. | ||
| 802 | ✗ | do i = 1, j | |
| 803 | ✗ | if (idx_arr(i) > idx_arr(i + 1)) then | |
| 804 | ✗ | temp = idx_arr(i) | |
| 805 | ✗ | idx_arr(i) = idx_arr(i + 1) | |
| 806 | ✗ | idx_arr(i + 1) = temp | |
| 807 | ✗ | swapped = .true. | |
| 808 | end if | ||
| 809 | end do | ||
| 810 | ✗ | if (.not. swapped) exit | |
| 811 | end do | ||
| 812 | |||
| 813 | ! eliminate the duplicates, replace them by either appending from higher | ||
| 814 | ! indices, dependent on available indices. | ||
| 815 | num_removed = 0 | ||
| 816 | ✗ | do o = 1, N | |
| 817 | ✗ | if ((N - o) < (num_removed + 1)) exit | |
| 818 | |||
| 819 | ✗ | if (idx_arr(o) == idx_arr(o + 1)) then | |
| 820 | ✗ | num_removed = num_removed + 1 | |
| 821 | ✗ | do r = o + 1, N - 1 | |
| 822 | ✗ | idx_arr(r) = idx_arr(r + 1) | |
| 823 | end do | ||
| 824 | end if | ||
| 825 | end do | ||
| 826 | |||
| 827 | ✗ | do p = N - num_removed + 1, N | |
| 828 | ✗ | idx_arr(p) = idx_arr(p - 1) + 1 | |
| 829 | end do | ||
| 830 | |||
| 831 | ✗ | if (idx_arr(N) > ntestpart) then | |
| 832 | ✗ | print *, 'ERROR - Invalid indices (Out of Range)!' | |
| 833 | ✗ | error stop | |
| 834 | end if | ||
| 835 | |||
| 836 | ✗ | end subroutine sort_idx | |
| 837 | |||
| 838 | 247 | function should_skip(ipart) | |
| 839 | ! notrace_passing: no tracing of passing particles, assume that all are confined | ||
| 840 | ! or skip strongly passing particles that are certainly confined | ||
| 841 | logical :: should_skip | ||
| 842 | integer, intent(in) :: ipart | ||
| 843 | |||
| 844 |
4/4✓ Branch 0 taken 207 times.
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 62 times.
✓ Branch 3 taken 145 times.
|
247 | should_skip = (notrace_passing .eq. 1) .or. (trap_par(ipart) .le. contr_pp) |
| 845 | 247 | end function should_skip | |
| 846 | |||
| 847 | 81 | subroutine apply_config_aliases | |
| 848 | ! Handle deprecated aliases and apply fallback logic for config parameters. | ||
| 849 | ! Priority: new name > deprecated alias > default | ||
| 850 | ! | ||
| 851 | ! Aliases: | ||
| 852 | ! netcdffile -> fallback for field_input and coord_input (deprecated) | ||
| 853 | ! isw_field_type -> integ_coords (deprecated) | ||
| 854 | ! | ||
| 855 | ! Fallback chain: | ||
| 856 | ! field_input: explicit > netcdffile > '' | ||
| 857 | ! coord_input: explicit > netcdffile > field_input > '' | ||
| 858 | |||
| 859 | ! netcdffile serves as fallback for both field_input and coord_input | ||
| 860 |
4/4✓ Branch 0 taken 38 times.
✓ Branch 1 taken 43 times.
✓ Branch 2 taken 31 times.
✓ Branch 3 taken 7 times.
|
81 | if (field_input == '' .and. len_trim(netcdffile) > 0) then |
| 861 | 31 | field_input = netcdffile | |
| 862 | end if | ||
| 863 | |||
| 864 |
4/4✓ Branch 0 taken 67 times.
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 31 times.
✓ Branch 3 taken 36 times.
|
81 | if (coord_input == '' .and. len_trim(netcdffile) > 0) then |
| 865 | 31 | coord_input = netcdffile | |
| 866 | end if | ||
| 867 | |||
| 868 | ! coord_input falls back to field_input if still not set | ||
| 869 |
2/2✓ Branch 0 taken 36 times.
✓ Branch 1 taken 45 times.
|
81 | if (coord_input == '') then |
| 870 | 36 | coord_input = field_input | |
| 871 | end if | ||
| 872 | |||
| 873 | ! Do not force netcdffile to follow coord_input: | ||
| 874 | ! - coord_input may point to a chartmap file (reference coordinates) | ||
| 875 | ! - netcdffile is still used as the VMEC equilibrium file by legacy code paths | ||
| 876 | ! The VMEC file selection is handled in simple_main.init_field. | ||
| 877 | |||
| 878 | ! isw_field_type is deprecated alias for integ_coords | ||
| 879 | ! integ_coords == -1000 means user did not set it | ||
| 880 |
2/2✓ Branch 0 taken 50 times.
✓ Branch 1 taken 31 times.
|
81 | if (integ_coords == -1000) then |
| 881 | 50 | integ_coords = isw_field_type | |
| 882 | else | ||
| 883 | 31 | isw_field_type = integ_coords | |
| 884 | end if | ||
| 885 | 81 | end subroutine apply_config_aliases | |
| 886 | |||
| 887 | 824 | pure logical function classification_enabled() | |
| 888 | !> Whether a run collects orbit classifications at all. Three call sites | ||
| 889 | !> must agree: the bmin/bmax cache the classifier needs to separate | ||
| 890 | !> trapped from passing, the dispatch to trace_orbit_with_classifiers, | ||
| 891 | !> and the class_parts.dat writer. | ||
| 892 | !> | ||
| 893 | !> fast_class belongs here. Without it, fast_class alone (tcut <= 0 and | ||
| 894 | !> class_plot = .False.) silently traced ordinary orbits and wrote no | ||
| 895 | !> classification output, even though the flag documents the opposite. | ||
| 896 | !> That combination is the only one giving fast classification without | ||
| 897 | !> the Minkowski fractal cut, since the cut fires at kt == ntcut and the | ||
| 898 | !> early exit in classification.f90 requires .not. class_plot. | ||
| 899 |
6/6✓ Branch 0 taken 815 times.
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 63 times.
✓ Branch 3 taken 752 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 742 times.
|
824 | classification_enabled = (ntcut > 0) .or. class_plot .or. fast_class |
| 900 | 824 | end function classification_enabled | |
| 901 | |||
| 902 | 139 | subroutine reset_seed_if_deterministic(stream_id) | |
| 903 | ! for run with fixed random seed | ||
| 904 | integer, intent(in), optional :: stream_id | ||
| 905 | integer :: seedsize, i, stream | ||
| 906 | 139 | integer, allocatable :: seed(:) | |
| 907 | integer(int64) :: state | ||
| 908 | integer(int64), parameter :: modulus = 2147483647_int64 | ||
| 909 | |||
| 910 |
2/2✓ Branch 0 taken 128 times.
✓ Branch 1 taken 11 times.
|
139 | if (deterministic) then |
| 911 | 128 | call random_seed(size=seedsize) | |
| 912 |
8/16✓ Branch 0 taken 128 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 128 times.
✓ Branch 4 taken 128 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 128 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 128 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 128 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 128 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 128 times.
|
384 | if (.not. allocated(seed)) allocate (seed(seedsize)) |
| 913 |
2/2✓ Branch 0 taken 75 times.
✓ Branch 1 taken 53 times.
|
128 | if (.not. present(stream_id)) then |
| 914 |
2/2✓ Branch 0 taken 600 times.
✓ Branch 1 taken 75 times.
|
675 | seed = ran_seed |
| 915 | 75 | call random_seed(put=seed) | |
| 916 | 75 | return | |
| 917 | end if | ||
| 918 | 53 | stream = stream_id | |
| 919 | state = modulo(int(ran_seed, int64) + 104729_int64*int(stream, int64), & | ||
| 920 | 53 | modulus - 1_int64) + 1_int64 | |
| 921 |
2/2✓ Branch 0 taken 424 times.
✓ Branch 1 taken 53 times.
|
477 | do i = 1, seedsize |
| 922 | state = modulo(48271_int64*state + 7919_int64*int(i, int64), & | ||
| 923 | 424 | modulus - 1_int64) + 1_int64 | |
| 924 | 477 | seed(i) = int(state) | |
| 925 | end do | ||
| 926 | 53 | call random_seed(put=seed) | |
| 927 | end if | ||
| 928 |
3/4✓ Branch 0 taken 75 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 53 times.
✓ Branch 3 taken 11 times.
|
214 | end subroutine reset_seed_if_deterministic |
| 929 | |||
| 930 | end module params | ||
| 931 |