| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | ✗ | program diag_traj_main | |
| 2 | !> Stream the trajectory of one particle through the PRODUCTION integrator | ||
| 3 | !> selected by integmode (RK orbit_timestep_axis for integmode<=0, symplectic | ||
| 4 | !> orbit_timestep_sympl otherwise). Unlike diag_orbit.x it does not hardcode a | ||
| 5 | !> solver and does not preallocate ntimstep*ntau points: it writes s, theta, phi, | ||
| 6 | !> pphi, H every `stride` substeps and reports min(s), loss time, and the number | ||
| 7 | !> of axis crossings (EVT_R_NEGATIVE). Built to localize the symplectic-vs-RK | ||
| 8 | !> non-prompt loss discrepancy on the internal Boozer field (#370). | ||
| 9 | !> | ||
| 10 | !> Usage: ./diag_traj.x [config_file] particle_number [stride] | ||
| 11 | |||
| 12 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 13 | use params, only: read_config, netcdffile, ns_s, ns_tp, multharm, integmode, & | ||
| 14 | params_init, dtaumin, relerr, ntestpart, zstart, startmode, grid_density, & | ||
| 15 | ✗ | special_ants_file, reuse_batch, num_surf, sbeg, trace_time, v0, & | |
| 16 | orbit_model, ORBIT_FULL_ORBIT | ||
| 17 | use simple, only: tracer_t, init_sympl, init_fo | ||
| 18 | use simple_main, only: init_field | ||
| 19 | use magfie_sub, only: init_magfie, VMEC | ||
| 20 | use samplers, only: init_starting_surf, sample, START_FILE | ||
| 21 | use field_can_mod, only: field_can_t, get_val, ref_to_integ | ||
| 22 | use orbit_symplectic, only: orbit_timestep_sympl | ||
| 23 | use orbit_symplectic_base, only: symplectic_integrator_t | ||
| 24 | use alpha_lifetime_sub, only: orbit_timestep_axis | ||
| 25 | use orbit_fo_boris, only: fo_step, fo_to_gc, fo_energy, FO_OK, FO_LOSS | ||
| 26 | use diag_counters, only: diag_counters_init, diag_counters_reset, & | ||
| 27 | diag_counters_total, EVT_R_NEGATIVE | ||
| 28 | use util, only: twopi | ||
| 29 | |||
| 30 | implicit none | ||
| 31 | |||
| 32 | character(256) :: config_file, arg | ||
| 33 | type(tracer_t) :: norb | ||
| 34 | type(symplectic_integrator_t) :: si | ||
| 35 | type(field_can_t) :: f | ||
| 36 | integer :: pnum, stride, nargs | ||
| 37 | real(dp), dimension(5) :: z | ||
| 38 | real(dp) :: t, s, smin, H | ||
| 39 | integer(8) :: kt | ||
| 40 | integer :: ierr, unit | ||
| 41 | character(64) :: fname | ||
| 42 | |||
| 43 | ✗ | config_file = 'simple.in' | |
| 44 | ✗ | stride = 1 | |
| 45 | ✗ | nargs = command_argument_count() | |
| 46 | ✗ | if (nargs == 1) then | |
| 47 | ✗ | call get_command_argument(1, arg); read (arg, *) pnum | |
| 48 | ✗ | elseif (nargs == 2) then | |
| 49 | ✗ | call get_command_argument(1, arg) | |
| 50 | ! second arg is either config (non-numeric) or stride; assume "pnum stride" | ||
| 51 | ✗ | read (arg, *) pnum | |
| 52 | ✗ | call get_command_argument(2, arg); read (arg, *) stride | |
| 53 | ✗ | elseif (nargs == 3) then | |
| 54 | ✗ | call get_command_argument(1, config_file) | |
| 55 | ✗ | call get_command_argument(2, arg); read (arg, *) pnum | |
| 56 | ✗ | call get_command_argument(3, arg); read (arg, *) stride | |
| 57 | else | ||
| 58 | ✗ | print *, 'Usage: ./diag_traj.x [config] particle_number [stride]' | |
| 59 | ✗ | stop | |
| 60 | end if | ||
| 61 | |||
| 62 | ✗ | call read_config(config_file) | |
| 63 | ✗ | call init_field(norb, netcdffile, ns_s, ns_tp, multharm, integmode) | |
| 64 | ✗ | call params_init | |
| 65 | ✗ | call init_magfie(VMEC) | |
| 66 | ✗ | call init_starting_surf | |
| 67 | ✗ | if (startmode == 2) then | |
| 68 | ✗ | call sample(zstart, START_FILE) | |
| 69 | ✗ | elseif (startmode == 5) then | |
| 70 | ✗ | if (num_surf == 1) then | |
| 71 | ✗ | call sample(zstart, 0.0d0, sbeg(1)) | |
| 72 | else | ||
| 73 | ✗ | call sample(zstart, sbeg(1), sbeg(num_surf)) | |
| 74 | end if | ||
| 75 | else | ||
| 76 | ✗ | call sample(zstart, START_FILE) | |
| 77 | end if | ||
| 78 | |||
| 79 | ✗ | if (pnum < 1 .or. pnum > ntestpart) then | |
| 80 | ✗ | print *, 'particle out of range', pnum, ntestpart; error stop | |
| 81 | end if | ||
| 82 | |||
| 83 | ✗ | call diag_counters_init() | |
| 84 | ✗ | call diag_counters_reset() | |
| 85 | |||
| 86 | ✗ | call ref_to_integ(zstart(1:3, pnum), z(1:3)) | |
| 87 | ✗ | z(4:5) = zstart(4:5, pnum) | |
| 88 | |||
| 89 | ✗ | if (orbit_model == ORBIT_FULL_ORBIT) then | |
| 90 | ✗ | call trace_full_orbit(norb, z, pnum, stride) | |
| 91 | ✗ | stop | |
| 92 | end if | ||
| 93 | |||
| 94 | ✗ | if (integmode > 0) call init_sympl(si, f, z, dtaumin, dtaumin, relerr, integmode) | |
| 95 | |||
| 96 | ✗ | write (fname, '(A,I0,A,I0,A)') 'traj_p', pnum, '_im', integmode, '.dat' | |
| 97 | ✗ | open (newunit=unit, file=trim(fname), status='replace') | |
| 98 | ✗ | write (unit, '(A)') '# t[s] s theta phi pphi/z4 H' | |
| 99 | |||
| 100 | ✗ | kt = 0; t = 0.0_dp; smin = z(1); ierr = 0 | |
| 101 | do | ||
| 102 | ✗ | if (integmode <= 0) then | |
| 103 | ✗ | call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, ierr) | |
| 104 | ✗ | s = z(1) | |
| 105 | ✗ | H = -1.0_dp | |
| 106 | else | ||
| 107 | ✗ | call orbit_timestep_sympl(si, f, ierr) | |
| 108 | ✗ | z(1:4) = si%z | |
| 109 | ✗ | s = si%z(1) | |
| 110 | ✗ | call get_val(f, si%z(4)) | |
| 111 | ✗ | H = f%H | |
| 112 | end if | ||
| 113 | ✗ | kt = kt + 1 | |
| 114 | ✗ | t = kt*dtaumin/v0 | |
| 115 | ✗ | smin = min(smin, s) | |
| 116 | ✗ | if (mod(kt, int(stride, 8)) == 0_8 .or. ierr /= 0) then | |
| 117 | ✗ | write (unit, '(6ES16.8)') t, s, mod(z(2), twopi), mod(z(3), twopi), z(4), H | |
| 118 | end if | ||
| 119 | ✗ | if (ierr /= 0) exit | |
| 120 | ✗ | if (t >= trace_time) exit | |
| 121 | end do | ||
| 122 | ✗ | close (unit) | |
| 123 | |||
| 124 | ✗ | print '(A,I0,A,I0)', 'particle ', pnum, ' integmode ', integmode | |
| 125 | ✗ | print '(A,L1,A,ES12.5)', 'lost = ', (ierr /= 0), ' t_end[s] = ', t | |
| 126 | ✗ | print '(A,ES12.5)', 'min(s) reached = ', smin | |
| 127 | ✗ | print '(A,I0)', 'axis crossings (R<0) = ', int(diag_counters_total(EVT_R_NEGATIVE)) | |
| 128 | ✗ | print '(A,A)', 'trajectory written : ', trim(fname) | |
| 129 | |||
| 130 | contains | ||
| 131 | |||
| 132 | ! Stream the full-orbit (orbit_model=7) trajectory of one particle, recording | ||
| 133 | ! the guiding-centre s and the particle s each step. A row's status flags 0 OK, | ||
| 134 | ! 1 guiding-centre loss (s>=1), 2 particle-field inversion fault. The point is | ||
| 135 | ! to see whether a terminal fault happens while the guiding centre is still | ||
| 136 | ! inside the plasma (a confined orbit miscounted) or genuinely at the edge. | ||
| 137 | ✗ | subroutine trace_full_orbit(norb, z, pnum, stride) | |
| 138 | type(tracer_t), intent(inout) :: norb | ||
| 139 | real(dp), intent(in) :: z(5) | ||
| 140 | integer, intent(in) :: pnum, stride | ||
| 141 | real(dp) :: zz(5), tt, sgc, spart, the, phi, vp, E0, sgcmin | ||
| 142 | integer(8) :: k | ||
| 143 | integer :: u, st, gst, code | ||
| 144 | character(64) :: fn | ||
| 145 | |||
| 146 | ✗ | zz = z | |
| 147 | ✗ | call init_fo(norb%fo, zz, dtaumin) | |
| 148 | ✗ | E0 = max(fo_energy(norb%fo), 1.0e-300_dp) | |
| 149 | ✗ | sgc = zz(1); spart = norb%fo%u(1)**2; the = zz(2); phi = zz(3); sgcmin = sgc | |
| 150 | |||
| 151 | ✗ | write (fn, '(A,I0,A)') 'traj_fo_p', pnum, '.dat' | |
| 152 | ✗ | open (newunit=u, file=trim(fn), status='replace') | |
| 153 | ✗ | write (u, '(A)') '# t[s] s_gc s_part theta phi E/E0 status' | |
| 154 | ✗ | write (u, '(6ES16.8,I3)') 0.0_dp, sgc, spart, the, phi, 1.0_dp, 0 | |
| 155 | |||
| 156 | ✗ | k = 0; gst = FO_OK | |
| 157 | ✗ | do | |
| 158 | ✗ | call fo_step(norb%fo, st) | |
| 159 | ✗ | k = k + 1 | |
| 160 | ✗ | tt = k*dtaumin/v0 | |
| 161 | ✗ | if (st /= FO_OK) then ! inversion fault | |
| 162 | ✗ | write (u, '(6ES16.8,I3)') tt, sgc, spart, the, phi, & | |
| 163 | ✗ | fo_energy(norb%fo)/E0, 2 | |
| 164 | ✗ | exit | |
| 165 | end if | ||
| 166 | ✗ | spart = norb%fo%u(1)**2 | |
| 167 | ✗ | call fo_to_gc(norb%fo, sgc, the, phi, vp, gst) | |
| 168 | ✗ | sgcmin = min(sgcmin, sgc) | |
| 169 | ✗ | code = 0 | |
| 170 | ✗ | if (gst == FO_LOSS) code = 1 | |
| 171 | ✗ | if (gst /= FO_OK .and. gst /= FO_LOSS) code = 2 | |
| 172 | ✗ | if (mod(k, int(stride, 8)) == 0_8 .or. gst /= FO_OK) & | |
| 173 | ✗ | write (u, '(6ES16.8,I3)') tt, sgc, spart, mod(the, twopi), & | |
| 174 | ✗ | mod(phi, twopi), fo_energy(norb%fo)/E0, code | |
| 175 | ✗ | if (gst /= FO_OK) exit | |
| 176 | ✗ | if (tt >= trace_time) exit | |
| 177 | end do | ||
| 178 | ✗ | close (u) | |
| 179 | |||
| 180 | ✗ | print '(A,I0,A)', 'particle ', pnum, ' orbit_model = 7 (full orbit)' | |
| 181 | ✗ | print '(A,ES12.5)', 't_end[s] = ', tt | |
| 182 | ✗ | if (st /= FO_OK) then | |
| 183 | ✗ | print '(A)', 'exit: particle inversion FAULT (not a physical loss)' | |
| 184 | ✗ | else if (gst == FO_LOSS) then | |
| 185 | ✗ | print '(A)', 'exit: guiding-centre loss (s>=1)' | |
| 186 | else | ||
| 187 | ✗ | print '(A)', 'exit: reached trace_time (confined)' | |
| 188 | end if | ||
| 189 | ✗ | print '(A,ES12.5,A,ES12.5)', 'min(s_gc) = ', sgcmin, ' last s_part = ', spart | |
| 190 | ✗ | print '(A,A)', 'trajectory written : ', trim(fn) | |
| 191 | ✗ | end subroutine trace_full_orbit | |
| 192 | |||
| 193 | end program diag_traj_main | ||
| 194 |