| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module diag_newton | ||
| 2 | !> Diagnostic routines for newton_midpoint function | ||
| 3 | !> Provides detailed analysis of Newton iteration convergence behavior | ||
| 4 | !> during orbit integration using midpoint rule | ||
| 5 | |||
| 6 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 7 | use util, only: twopi | ||
| 8 | use field_can_mod, only: field_can_t, get_val, eval_field => evaluate | ||
| 9 | use orbit_symplectic_base, only: symplectic_integrator_t | ||
| 10 | use vector_potentail_mod, only: torflux | ||
| 11 | use lapack_interfaces, only: dgesv | ||
| 12 | use params, only: dtau | ||
| 13 | use orbit_symplectic, only: f_midpoint_part1, f_midpoint_part2, & | ||
| 14 | jac_midpoint_part1, jac_midpoint_part2 | ||
| 15 | |||
| 16 | implicit none | ||
| 17 | private | ||
| 18 | |||
| 19 | public :: newton_midpoint_debug, integrate_orbit_with_newton_debug | ||
| 20 | |||
| 21 | contains | ||
| 22 | |||
| 23 | !> EXACT copy of newton_midpoint from orbit_symplectic.f90 with debug output | ||
| 24 | ✗ | subroutine newton_midpoint_debug(si, f, x, atol, rtol, maxit, xlast, step_num) | |
| 25 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 26 | type(field_can_t), intent(inout) :: f | ||
| 27 | type(field_can_t) :: fmid | ||
| 28 | integer, parameter :: n = 5 | ||
| 29 | integer :: kit | ||
| 30 | real(dp), intent(inout) :: x(n) ! = (rend, thend, phend, pphend, rmid) | ||
| 31 | real(dp), intent(in) :: atol, rtol | ||
| 32 | integer, intent(in) :: maxit | ||
| 33 | real(dp), intent(out) :: xlast(n) | ||
| 34 | integer, intent(in) :: step_num | ||
| 35 | real(dp) :: fvec(n), fjac(n, n) | ||
| 36 | integer :: pivot(n), info | ||
| 37 | real(dp) :: xabs(n), tolref(n), fabs(n) | ||
| 38 | |||
| 39 | ! DEBUG OUTPUT | ||
| 40 | ✗ | print '(A,I0)', 'Newton Midpoint Debug - Step ', step_num | |
| 41 | ✗ | print '(A,ES12.5,A,ES12.5)', 'Tolerances: atol = ', atol, ', rtol = ', rtol | |
| 42 | |||
| 43 | ✗ | tolref(1) = 1d0 | |
| 44 | ✗ | tolref(2) = twopi | |
| 45 | ✗ | tolref(3) = twopi | |
| 46 | ✗ | tolref(4) = dabs(1d1*torflux/f%ro0) | |
| 47 | ✗ | tolref(5) = 1d0 | |
| 48 | |||
| 49 | ✗ | print '(A,5ES12.5)', 'Initial Tolref = [', tolref, ']' | |
| 50 | ✗ | print '(A,5ES12.5)', 'Initial x = [', x, ']' | |
| 51 | ✗ | print * | |
| 52 | |||
| 53 | ✗ | do kit = 1, maxit | |
| 54 | ✗ | if (x(1) > 1.0) return | |
| 55 | ✗ | if (x(1) < 0.0) x(1) = 0.01 | |
| 56 | ✗ | if (x(5) < 0.0) x(5) = 0.01 | |
| 57 | ✗ | call f_midpoint_part1(si, f, n, x, fvec) | |
| 58 | ✗ | call jac_midpoint_part1(si, f, x, fjac) | |
| 59 | ✗ | fmid = f | |
| 60 | ✗ | call f_midpoint_part2(si, f, n, x, fvec) | |
| 61 | ✗ | call jac_midpoint_part2(si, f, fmid, x, fjac) | |
| 62 | ✗ | fabs = dabs(fvec) | |
| 63 | ✗ | xlast = x | |
| 64 | ✗ | call dgesv(n, 1, fjac, n, pivot, fvec, n, info) | |
| 65 | ! after solution: fvec = (xold-xnew)_Newton | ||
| 66 | ✗ | x = x - fvec | |
| 67 | ✗ | xabs = dabs(x - xlast) | |
| 68 | ! Don't take too small values in pphi as tolerance reference | ||
| 69 | ✗ | tolref(4) = max(dabs(x(4)), tolref(4)) | |
| 70 | |||
| 71 | ! DEBUG OUTPUT | ||
| 72 | ✗ | print '(A,I0,A,5ES12.5)', 'Iteration ', kit, ': fabs = [', fabs, ']' | |
| 73 | ✗ | print '(A,I0,A,5ES12.5)', 'Iteration ', kit, ': xabs = [', xabs, ']' | |
| 74 | ✗ | print '(A,5ES12.5)', 'Tolref = [', tolref, ']' | |
| 75 | ✗ | print '(A,5ES12.5)', 'rtol*tolref = [', rtol*tolref, ']' | |
| 76 | |||
| 77 | ✗ | if (all(fabs < atol)) then | |
| 78 | ✗ | print '(A,I0,A)', 'Iteration ', kit, ': Convergence achieved (fabs < atol)' | |
| 79 | ✗ | return | |
| 80 | end if | ||
| 81 | ✗ | if (all(xabs < rtol*tolref)) then | |
| 82 | ✗ | print '(A,I0,A)', 'Iteration ', kit, ': Convergence achieved (xabs < rtol*tolref)' | |
| 83 | ✗ | return | |
| 84 | end if | ||
| 85 | |||
| 86 | ✗ | print '(A,I0,A,5ES12.5)', 'Iteration ', kit, ': Updated x = [', x, ']' | |
| 87 | ✗ | print * | |
| 88 | end do | ||
| 89 | ✗ | print '(A,I0)', 'newton_midpoint: maximum iterations reached: ', maxit | |
| 90 | !write(6603,*) x(1), x(2), x(3), x(4), x(5), xabs, fvec | ||
| 91 | ! TODO fix criterion for convergence | ||
| 92 | ✗ | end subroutine newton_midpoint_debug | |
| 93 | |||
| 94 | !> Integration wrapper that calls debug newton_midpoint for specified number of steps | ||
| 95 | ✗ | subroutine integrate_orbit_with_newton_debug(si, f, num_steps) | |
| 96 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 97 | type(field_can_t), intent(inout) :: f | ||
| 98 | integer, intent(in) :: num_steps | ||
| 99 | |||
| 100 | integer, parameter :: n = 5, maxit = 32 | ||
| 101 | real(dp), dimension(n) :: x, xlast | ||
| 102 | integer :: step | ||
| 103 | |||
| 104 | ✗ | print *, 'Starting Newton Midpoint Integration Debug' | |
| 105 | ✗ | print '(A,I0,A)', 'Will integrate for ', num_steps, ' time steps' | |
| 106 | ✗ | print '(A,ES12.5)', 'dtau (large time step): ', dtau | |
| 107 | ✗ | print '(A,ES12.5)', 'dtaumin (integration time step): ', si%dt | |
| 108 | ✗ | print '(A,I0)', 'ntau (substeps per dtau): ', si%ntau | |
| 109 | ✗ | print '(A,ES12.5)', 'Absolute tolerance: ', si%atol | |
| 110 | ✗ | print '(A,ES12.5)', 'Relative tolerance: ', si%rtol | |
| 111 | ✗ | print '(A,4ES12.5)', 'Initial conditions: ', si%z | |
| 112 | ✗ | print * | |
| 113 | |||
| 114 | ✗ | do step = 1, num_steps | |
| 115 | ✗ | si%pthold = f%pth | |
| 116 | |||
| 117 | ✗ | x(1:4) = si%z | |
| 118 | ✗ | x(5) = si%z(1) | |
| 119 | |||
| 120 | ✗ | call newton_midpoint_debug(si, f, x, si%atol, si%rtol, maxit, xlast, step) | |
| 121 | |||
| 122 | ✗ | if (x(1) > 1.0_dp) then | |
| 123 | ✗ | print *, 'Particle lost: s > 1.0 at step ', step | |
| 124 | ✗ | exit | |
| 125 | end if | ||
| 126 | |||
| 127 | ✗ | si%z = x(1:4) | |
| 128 | |||
| 129 | ! Update field | ||
| 130 | ✗ | call eval_field(f, si%z(1), si%z(2), si%z(3), 0) | |
| 131 | ✗ | call get_val(f, si%z(4)) | |
| 132 | |||
| 133 | ✗ | print '(A,I0,A,4ES12.5)', 'Step ', step, ' completed. Final state: ', si%z | |
| 134 | ✗ | print * | |
| 135 | end do | ||
| 136 | |||
| 137 | ✗ | print *, 'Newton Midpoint Integration Debug completed successfully!' | |
| 138 | |||
| 139 | ✗ | end subroutine integrate_orbit_with_newton_debug | |
| 140 | |||
| 141 | end module diag_newton | ||
| 142 |