| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module diag_orbit | ||
| 2 | !> Diagnostic routines for plotting single orbit trajectories | ||
| 3 | !> Provides trajectory plotting functionality for the Nth particle using | ||
| 4 | !> full SIMPLE initialization and real orbit integration | ||
| 5 | |||
| 6 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 7 | use params, only: dtaumin, ntestpart, ntimstep, ntau, zstart, & | ||
| 8 | special_ants_file, reuse_batch, num_surf, sbeg, integmode, relerr, reset_seed_if_deterministic | ||
| 9 | use field_can_mod, only: field_can_t, get_val, eval_field => evaluate, ref_to_integ | ||
| 10 | use orbit_symplectic_base, only: symplectic_integrator_t, extrap_field | ||
| 11 | use orbit_symplectic, only: f_midpoint_part1, f_midpoint_part2, & | ||
| 12 | jac_midpoint_part1, jac_midpoint_part2 | ||
| 13 | use simple, only: init_sympl | ||
| 14 | use vector_potentail_mod, only: torflux | ||
| 15 | use lapack_interfaces, only: dgesv | ||
| 16 | use util, only: twopi | ||
| 17 | |||
| 18 | implicit none | ||
| 19 | private | ||
| 20 | |||
| 21 | public :: integrate_orbit_with_trajectory_debug | ||
| 22 | |||
| 23 | contains | ||
| 24 | |||
| 25 | !> Newton midpoint solver that returns iteration count (no debug output) | ||
| 26 | ✗ | function newton_midpoint_count_iterations(si, f, x, atol, rtol, maxit, xlast, field_evals) result(iterations) | |
| 27 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 28 | type(field_can_t), intent(inout) :: f | ||
| 29 | type(field_can_t) :: fmid | ||
| 30 | integer, parameter :: n = 5 | ||
| 31 | integer :: kit, iterations | ||
| 32 | real(dp), intent(inout) :: x(n) ! = (rend, thend, phend, pphend, rmid) | ||
| 33 | real(dp), intent(in) :: atol, rtol | ||
| 34 | integer, intent(in) :: maxit | ||
| 35 | real(dp), intent(out) :: xlast(n) | ||
| 36 | integer, intent(inout) :: field_evals | ||
| 37 | real(dp) :: fvec(n), fjac(n, n) | ||
| 38 | integer :: pivot(n), info | ||
| 39 | real(dp) :: xabs(n), tolref(n), fabs(n) | ||
| 40 | |||
| 41 | ! Buffers to store all iteration data (only printed if max iterations reached) | ||
| 42 | ✗ | real(dp) :: x_buffer(n, maxit), fabs_buffer(n, maxit), xabs_buffer(n, maxit) | |
| 43 | real(dp) :: x_initial(n) | ||
| 44 | integer :: k | ||
| 45 | |||
| 46 | ✗ | tolref(1) = 1d0 | |
| 47 | ✗ | tolref(2) = twopi | |
| 48 | ✗ | tolref(3) = twopi | |
| 49 | ✗ | tolref(4) = max(dabs(f%Aph), dabs(1d1*torflux/f%ro0)) ! Use actual Aph from field | |
| 50 | ✗ | tolref(5) = 1d0 | |
| 51 | |||
| 52 | ! Store initial conditions | ||
| 53 | ✗ | x_initial = x | |
| 54 | |||
| 55 | ✗ | do kit = 1, maxit | |
| 56 | ✗ | if (x(1) > 1.0) then | |
| 57 | ✗ | iterations = kit - 1 | |
| 58 | ✗ | return | |
| 59 | end if | ||
| 60 | ✗ | if (x(1) < 0.0) x(1) = 0.01 | |
| 61 | ✗ | if (x(5) < 0.0) x(5) = 0.01 | |
| 62 | ✗ | call f_midpoint_part1(si, f, n, x, fvec) | |
| 63 | ✗ | call jac_midpoint_part1(si, f, x, fjac) | |
| 64 | ✗ | fmid = f | |
| 65 | ✗ | call f_midpoint_part2(si, f, n, x, fvec) | |
| 66 | ✗ | call jac_midpoint_part2(si, f, fmid, x, fjac) | |
| 67 | ! Each Newton iteration involves multiple field evaluations | ||
| 68 | ! f_midpoint_part1 and f_midpoint_part2 each do field evaluations | ||
| 69 | ✗ | field_evals = field_evals + 2 | |
| 70 | ✗ | fabs = dabs(fvec) | |
| 71 | ✗ | xlast = x | |
| 72 | ✗ | call dgesv(n, 1, fjac, n, pivot, fvec, n, info) | |
| 73 | ✗ | x = x - fvec | |
| 74 | ✗ | xabs = dabs(x - xlast) | |
| 75 | |||
| 76 | ! Store iteration data in buffers | ||
| 77 | ✗ | x_buffer(:, kit) = x | |
| 78 | ✗ | fabs_buffer(:, kit) = fabs | |
| 79 | ✗ | xabs_buffer(:, kit) = xabs | |
| 80 | |||
| 81 | ! Use reasonable absolute tolerance instead of machine epsilon | ||
| 82 | ✗ | if (all(fabs < atol)) then | |
| 83 | ✗ | iterations = kit | |
| 84 | ✗ | return | |
| 85 | end if | ||
| 86 | ✗ | if (all(xabs < rtol*tolref)) then | |
| 87 | ✗ | iterations = kit | |
| 88 | ✗ | return | |
| 89 | end if | ||
| 90 | end do | ||
| 91 | |||
| 92 | ! Maximum iterations reached - print complete iteration history | ||
| 93 | ✗ | write (*, '(A)') '=== NEWTON SOLVER FAILURE: MAXIMUM ITERATIONS REACHED ===' | |
| 94 | ✗ | write (*, '(A,I0)') 'Maximum iterations: ', maxit | |
| 95 | ✗ | write (*, '(A,5ES12.5)') 'Initial x = [', x_initial, ']' | |
| 96 | ✗ | write (*, *) | |
| 97 | ✗ | write (*, '(A)') 'Complete iteration history:' | |
| 98 | ✗ | write (*, '(A)') 'Iter | max(fabs) | max(xabs) | Result' | |
| 99 | ✗ | write (*, '(A)') '-----|----------------|----------------|-------' | |
| 100 | |||
| 101 | ✗ | do k = 1, maxit | |
| 102 | ✗ | write(*,'(I4,A,ES12.5,A,ES12.5,A)',advance='no') k, ' | ', maxval(fabs_buffer(:,k)), & | |
| 103 | ✗ | ' | ', maxval(xabs_buffer(:, k)), ' | ' | |
| 104 | |||
| 105 | ✗ | if (all(fabs_buffer(:, k) < atol)) then | |
| 106 | ✗ | write (*, '(A)') 'fabs < atol' | |
| 107 | ✗ | elseif (all(xabs_buffer(:, k) < rtol*tolref)) then | |
| 108 | ✗ | write (*, '(A)') 'xabs < rtol*tolref' | |
| 109 | else | ||
| 110 | ✗ | write (*, '(A)') 'continuing...' | |
| 111 | end if | ||
| 112 | end do | ||
| 113 | |||
| 114 | ✗ | write (*, *) | |
| 115 | ✗ | write (*, '(A,5ES12.5)') 'Final fabs = [', fabs_buffer(:, maxit), ']' | |
| 116 | ✗ | write (*, '(A,5ES12.5)') 'Final xabs = [', xabs_buffer(:, maxit), ']' | |
| 117 | ✗ | write (*, '(A,5ES12.5)') 'rtol*tolref= [', rtol*tolref, ']' | |
| 118 | ✗ | write (*, '(A,5ES12.5)') 'Final x = [', x_buffer(:, maxit), ']' | |
| 119 | ✗ | write (*, *) | |
| 120 | |||
| 121 | ✗ | iterations = maxit | |
| 122 | ✗ | error stop 'Newton solver failed to converge within maximum iterations' | |
| 123 | ✗ | end function newton_midpoint_count_iterations | |
| 124 | |||
| 125 | !> Integration wrapper that plots the trajectory of the Nth particle | ||
| 126 | ✗ | subroutine integrate_orbit_with_trajectory_debug(si, f, particle_number) | |
| 127 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 128 | type(field_can_t), intent(inout) :: f | ||
| 129 | integer, intent(in) :: particle_number | ||
| 130 | |||
| 131 | ✗ | real(dp), allocatable :: s_traj(:), theta_traj(:), phi_traj(:), time_traj(:) | |
| 132 | ✗ | real(dp), allocatable :: pphi_traj(:) | |
| 133 | ✗ | integer, allocatable :: newton_iter_traj(:) | |
| 134 | real(dp), dimension(5) :: z, xlast | ||
| 135 | integer :: it, ktau, point_idx, newton_iters | ||
| 136 | integer, parameter :: maxit = 32 | ||
| 137 | real(dp) :: current_time | ||
| 138 | integer :: total_points, field_eval_count | ||
| 139 | |||
| 140 | ! Validate particle number | ||
| 141 | ✗ | if (particle_number < 1 .or. particle_number > ntestpart) then | |
| 142 | ✗ | print '(A,I0,A,I0)', 'ERROR: Invalid particle number ', particle_number, & | |
| 143 | ✗ | '. Must be between 1 and ', ntestpart | |
| 144 | return | ||
| 145 | end if | ||
| 146 | |||
| 147 | ! CRITICAL: Follow simple_main.f90 trace_orbit EXACTLY | ||
| 148 | ! 1. Reset random seed if deterministic | ||
| 149 | ✗ | call reset_seed_if_deterministic | |
| 150 | |||
| 151 | ! 2. Get particle coordinates and transform ref -> integ (CRITICAL STEP MISSING!) | ||
| 152 | ✗ | call ref_to_integ(zstart(1:3, particle_number), z(1:3)) | |
| 153 | ✗ | z(4:5) = zstart(4:5, particle_number) | |
| 154 | |||
| 155 | ! 3. Initialize symplectic integrator with TRANSFORMED coordinates | ||
| 156 | ✗ | if (integmode > 0) then | |
| 157 | ✗ | call init_sympl(si, f, z, dtaumin, dtaumin, relerr, integmode) | |
| 158 | end if | ||
| 159 | |||
| 160 | ✗ | current_time = 0.0_dp | |
| 161 | |||
| 162 | ! Calculate total number of timesteps (macrosteps * substeps + initial) | ||
| 163 | ✗ | total_points = ntimstep*ntau + 1 | |
| 164 | field_eval_count = 0 | ||
| 165 | |||
| 166 | ! Allocate trajectory arrays | ||
| 167 | ✗ | allocate (s_traj(total_points)) | |
| 168 | ✗ | allocate (theta_traj(total_points)) | |
| 169 | ✗ | allocate (phi_traj(total_points)) | |
| 170 | ✗ | allocate (pphi_traj(total_points)) | |
| 171 | ✗ | allocate (newton_iter_traj(total_points)) | |
| 172 | ✗ | allocate (time_traj(total_points)) | |
| 173 | |||
| 174 | ! Store initial conditions (0 Newton iterations for initial state) | ||
| 175 | ✗ | s_traj(1) = si%z(1) | |
| 176 | ✗ | theta_traj(1) = si%z(2) | |
| 177 | ✗ | phi_traj(1) = si%z(3) | |
| 178 | ✗ | pphi_traj(1) = si%z(4) | |
| 179 | ✗ | newton_iter_traj(1) = 0 | |
| 180 | ✗ | time_traj(1) = current_time | |
| 181 | |||
| 182 | ! Initialize field at starting position | ||
| 183 | ✗ | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | |
| 184 | ✗ | call get_val(f, si%z(4)) | |
| 185 | ✗ | field_eval_count = field_eval_count + 1 | |
| 186 | |||
| 187 | ✗ | point_idx = 1 | |
| 188 | |||
| 189 | ! Initialize xlast for field extrapolation to current coordinates | ||
| 190 | ✗ | xlast(1:4) = si%z | |
| 191 | ✗ | xlast(5) = si%z(1) | |
| 192 | |||
| 193 | ! Use our custom Newton solver to get iteration counts (but with proper initialization now) | ||
| 194 | ✗ | do it = 1, ntimstep | |
| 195 | ✗ | do ktau = 1, ntau | |
| 196 | ✗ | si%pthold = f%pth | |
| 197 | |||
| 198 | ! Set up for midpoint integration (like diag_newton) | ||
| 199 | ✗ | z(1:4) = si%z | |
| 200 | ✗ | z(5) = si%z(1) | |
| 201 | |||
| 202 | ! Use custom Newton midpoint solver to get iteration count | ||
| 203 | ✗ | newton_iters = newton_midpoint_count_iterations(si, f, z, si%atol, si%rtol, maxit, xlast, field_eval_count) | |
| 204 | |||
| 205 | ✗ | current_time = current_time + dtaumin | |
| 206 | |||
| 207 | ✗ | if (z(1) > 1.0_dp) then | |
| 208 | exit | ||
| 209 | end if | ||
| 210 | |||
| 211 | ! Update integrator state | ||
| 212 | ✗ | si%z = z(1:4) | |
| 213 | |||
| 214 | ! Store trajectory point | ||
| 215 | ✗ | point_idx = point_idx + 1 | |
| 216 | ✗ | s_traj(point_idx) = si%z(1) | |
| 217 | ✗ | theta_traj(point_idx) = si%z(2) | |
| 218 | ✗ | phi_traj(point_idx) = si%z(3) | |
| 219 | ✗ | pphi_traj(point_idx) = si%z(4) | |
| 220 | ✗ | newton_iter_traj(point_idx) = newton_iters | |
| 221 | ✗ | time_traj(point_idx) = current_time | |
| 222 | |||
| 223 | ! Update field with extrapolation like production integrator | ||
| 224 | ✗ | if (extrap_field) then | |
| 225 | f%pth = f%pth + f%dpth(1)*(z(1) - xlast(1) + z(5) - xlast(5)) & ! d/dr | ||
| 226 | + f%dpth(2)*(z(2) - xlast(2)) & ! d/dth | ||
| 227 | + f%dpth(3)*(z(3) - xlast(3)) & ! d/dph | ||
| 228 | ✗ | + f%dpth(4)*(z(4) - xlast(4)) ! d/dpph | |
| 229 | else | ||
| 230 | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | ||
| 231 | call get_val(f, si%z(4)) | ||
| 232 | field_eval_count = field_eval_count + 1 | ||
| 233 | end if | ||
| 234 | end do | ||
| 235 | ✗ | if (z(1) > 1.0_dp) exit | |
| 236 | end do | ||
| 237 | |||
| 238 | ! Write trajectory data to files for external plotting | ||
| 239 | call write_trajectory_data(time_traj(1:point_idx), s_traj(1:point_idx), & | ||
| 240 | theta_traj(1:point_idx), phi_traj(1:point_idx), pphi_traj(1:point_idx), & | ||
| 241 | ✗ | newton_iter_traj(1:point_idx), point_idx, particle_number) | |
| 242 | |||
| 243 | ! Output field evaluation count | ||
| 244 | ✗ | print '(A,I0)', 'Total field evaluations: ', field_eval_count | |
| 245 | |||
| 246 | ! Cleanup | ||
| 247 | ✗ | deallocate (s_traj, theta_traj, phi_traj, pphi_traj, newton_iter_traj, time_traj) | |
| 248 | |||
| 249 | end subroutine integrate_orbit_with_trajectory_debug | ||
| 250 | |||
| 251 | ✗ | subroutine write_trajectory_data(time_traj, s_traj, theta_traj, phi_traj, pphi_traj, & | |
| 252 | ✗ | newton_iter_traj, npoints, particle_number) | |
| 253 | integer, intent(in) :: npoints, particle_number | ||
| 254 | real(dp), dimension(npoints), intent(in) :: time_traj, s_traj, theta_traj, phi_traj, pphi_traj | ||
| 255 | integer, dimension(npoints), intent(in) :: newton_iter_traj | ||
| 256 | |||
| 257 | integer :: i | ||
| 258 | character(len=100) :: filename | ||
| 259 | |||
| 260 | ✗ | write (filename, '(A,I0,A)') 'orbit_trajectory_particle_', particle_number, '.dat' | |
| 261 | |||
| 262 | ✗ | open (unit=20, file=filename, status='replace') | |
| 263 | ✗ | write (20, '(A)') '# Time s theta phi pphi newton_iters' | |
| 264 | ✗ | do i = 1, npoints | |
| 265 | ✗ | write(20, '(5ES16.8,I8)') time_traj(i), s_traj(i), theta_traj(i), phi_traj(i), pphi_traj(i), newton_iter_traj(i) | |
| 266 | end do | ||
| 267 | ✗ | close (20) | |
| 268 | |||
| 269 | ✗ | end subroutine write_trajectory_data | |
| 270 | |||
| 271 | end module diag_orbit | ||
| 272 |