| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module orbit_symplectic_euler1 | ||
| 2 | use field_can_mod, only: field_can_t | ||
| 3 | use orbit_symplectic_base, only: symplectic_integrator_t, sympl_rmax | ||
| 4 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 5 | |||
| 6 | implicit none | ||
| 7 | |||
| 8 | private | ||
| 9 | |||
| 10 | public :: sympl_euler1_residual, sympl_euler1_jacobian | ||
| 11 | public :: sympl_euler1_newton_iter, sympl_euler1_extrapolate_field | ||
| 12 | public :: sympl_euler1_advance_angles | ||
| 13 | |||
| 14 | contains | ||
| 15 | |||
| 16 | 13030884 | subroutine sympl_euler1_residual(si, f, x, fvec) | |
| 17 | !$acc routine seq | ||
| 18 | type(symplectic_integrator_t), intent(in) :: si | ||
| 19 | type(field_can_t), intent(in) :: f | ||
| 20 | real(dp), intent(in) :: x(2) | ||
| 21 | real(dp), intent(out) :: fvec(2) | ||
| 22 | |||
| 23 | fvec(1) = f%dpth(1)*(f%pth - si%pthold) & | ||
| 24 | 13030884 | + si%dt*(f%dH(2)*f%dpth(1) - f%dH(1)*f%dpth(2)) | |
| 25 | fvec(2) = f%dpth(1)*(x(2) - si%z(4)) & | ||
| 26 | 13030884 | + si%dt*(f%dH(3)*f%dpth(1) - f%dH(1)*f%dpth(3)) | |
| 27 | 13030884 | end subroutine sympl_euler1_residual | |
| 28 | |||
| 29 | 13030884 | subroutine sympl_euler1_jacobian(si, f, x, jac) | |
| 30 | !$acc routine seq | ||
| 31 | type(symplectic_integrator_t), intent(in) :: si | ||
| 32 | type(field_can_t), intent(in) :: f | ||
| 33 | real(dp), intent(in) :: x(2) | ||
| 34 | real(dp), intent(out) :: jac(2, 2) | ||
| 35 | |||
| 36 | jac(1, 1) = f%d2pth(1)*(f%pth - si%pthold) + f%dpth(1)**2 & | ||
| 37 | + si%dt*(f%d2H(2)*f%dpth(1) + f%dH(2)*f%d2pth(1) & | ||
| 38 | 13030884 | - f%d2H(1)*f%dpth(2) - f%dH(1)*f%d2pth(2)) | |
| 39 | jac(1, 2) = f%d2pth(7)*(f%pth - si%pthold) + f%dpth(1)*f%dpth(4) & | ||
| 40 | + si%dt*(f%d2H(8)*f%dpth(1) + f%dH(2)*f%d2pth(7) & | ||
| 41 | 13030884 | - f%d2H(7)*f%dpth(2) - f%dH(1)*f%d2pth(8)) | |
| 42 | jac(2, 1) = f%d2pth(1)*(x(2) - si%z(4)) & | ||
| 43 | + si%dt*(f%d2H(3)*f%dpth(1) + f%dH(3)*f%d2pth(1) & | ||
| 44 | 13030884 | - f%d2H(1)*f%dpth(3) - f%dH(1)*f%d2pth(3)) | |
| 45 | jac(2, 2) = f%d2pth(7)*(x(2) - si%z(4)) + f%dpth(1) & | ||
| 46 | + si%dt*(f%d2H(9)*f%dpth(1) + f%dH(3)*f%d2pth(7) & | ||
| 47 | 13030884 | - f%d2H(7)*f%dpth(3) - f%dH(1)*f%d2pth(9)) | |
| 48 | 13030884 | end subroutine sympl_euler1_jacobian | |
| 49 | |||
| 50 | 13030884 | subroutine sympl_euler1_newton_iter(si, f, x, tolref, xlast, converged, & | |
| 51 | linear_failed, boundary_limited) | ||
| 52 | !$acc routine seq | ||
| 53 | type(symplectic_integrator_t), intent(in) :: si | ||
| 54 | type(field_can_t), intent(in) :: f | ||
| 55 | real(dp), intent(inout) :: x(2) | ||
| 56 | real(dp), intent(inout) :: tolref(2) | ||
| 57 | real(dp), intent(out) :: xlast(2) | ||
| 58 | logical, intent(out) :: converged | ||
| 59 | logical, intent(out) :: linear_failed | ||
| 60 | logical, intent(out) :: boundary_limited | ||
| 61 | |||
| 62 | real(dp) :: fvec(2), fjac(2, 2), ijac(2, 2) | ||
| 63 | real(dp) :: correction(2), determinant, matrix_scale, step_scale | ||
| 64 | |||
| 65 | 13030884 | call sympl_euler1_residual(si, f, x, fvec) | |
| 66 | 13030884 | call sympl_euler1_jacobian(si, f, x, fjac) | |
| 67 | |||
| 68 | 13030884 | determinant = fjac(1, 1)*fjac(2, 2) - fjac(1, 2)*fjac(2, 1) | |
| 69 |
12/12✓ Branch 0 taken 26061768 times.
✓ Branch 1 taken 13030884 times.
✓ Branch 2 taken 52123536 times.
✓ Branch 3 taken 26061768 times.
✓ Branch 4 taken 39092556 times.
✓ Branch 5 taken 13030980 times.
✓ Branch 6 taken 2126 times.
✓ Branch 7 taken 39090430 times.
✓ Branch 8 taken 13030852 times.
✓ Branch 9 taken 128 times.
✓ Branch 10 taken 32 times.
✓ Branch 11 taken 13030852 times.
|
91216188 | matrix_scale = maxval(abs(fjac)) |
| 70 | linear_failed = matrix_scale == 0d0 .or. & | ||
| 71 |
2/4✓ Branch 0 taken 13030884 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 13030884 times.
|
13030884 | abs(determinant) <= epsilon(1d0)*matrix_scale**2 |
| 72 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 13030884 times.
|
13030884 | if (linear_failed) then |
| 73 | ✗ | converged = .false. | |
| 74 | ✗ | boundary_limited = .false. | |
| 75 | ✗ | xlast = x | |
| 76 | ✗ | return | |
| 77 | end if | ||
| 78 | |||
| 79 |
8/8✓ Branch 0 taken 26059896 times.
✓ Branch 1 taken 13029012 times.
✓ Branch 2 taken 52119792 times.
✓ Branch 3 taken 26058024 times.
✓ Branch 4 taken 52117920 times.
✓ Branch 5 taken 1872 times.
✓ Branch 6 taken 13029012 times.
✓ Branch 7 taken 1872 times.
|
91206828 | if (all(fjac /= 0d0)) then |
| 80 | ijac(1, 1) = 1d0/(fjac(1, 1) - & | ||
| 81 | 13029012 | fjac(1, 2)*fjac(2, 1)/fjac(2, 2)) | |
| 82 | 13029012 | ijac(1, 2) = -1d0/(fjac(1, 1)*fjac(2, 2)/fjac(1, 2) - fjac(2, 1)) | |
| 83 | 13029012 | ijac(2, 1) = -1d0/(fjac(1, 1)*fjac(2, 2)/fjac(2, 1) - fjac(1, 2)) | |
| 84 | ijac(2, 2) = 1d0/(fjac(2, 2) - & | ||
| 85 | 13029012 | fjac(1, 2)*fjac(2, 1)/fjac(1, 1)) | |
| 86 | else | ||
| 87 | 1872 | ijac(1, 1) = fjac(2, 2)/determinant | |
| 88 | 1872 | ijac(1, 2) = -fjac(1, 2)/determinant | |
| 89 | 1872 | ijac(2, 1) = -fjac(2, 1)/determinant | |
| 90 | 1872 | ijac(2, 2) = fjac(1, 1)/determinant | |
| 91 | end if | ||
| 92 | |||
| 93 | 13030884 | xlast = x | |
| 94 | 13030884 | correction(1) = ijac(1, 1)*fvec(1) + ijac(1, 2)*fvec(2) | |
| 95 | 13030884 | correction(2) = ijac(2, 1)*fvec(1) + ijac(2, 2)*fvec(2) | |
| 96 | 13030884 | step_scale = 1d0 | |
| 97 |
4/4✓ Branch 0 taken 11990168 times.
✓ Branch 1 taken 1040716 times.
✓ Branch 2 taken 5947273 times.
✓ Branch 3 taken 6042895 times.
|
13030884 | if (sympl_rmax <= 1d0 .and. correction(1) < 0d0) then |
| 98 | step_scale = min(1d0, & | ||
| 99 | 5947273 | 0.8d0*max(0d0, sympl_rmax - x(1))/(-correction(1))) | |
| 100 | end if | ||
| 101 | 13030884 | boundary_limited = step_scale < 1d0 | |
| 102 |
2/2✓ Branch 0 taken 26061768 times.
✓ Branch 1 taken 13030884 times.
|
39092652 | x = x - step_scale*correction |
| 103 | |||
| 104 | 13030884 | tolref(2) = max(dabs(x(2)), tolref(2)) | |
| 105 | |||
| 106 | converged = .not. boundary_limited .and. & | ||
| 107 | ((dabs(fvec(1)) < si%atol .and. dabs(fvec(2)) < si%atol) .or. & | ||
| 108 | ( & | ||
| 109 | dabs(x(1) - xlast(1)) < si%rtol*tolref(1) .and. & | ||
| 110 |
10/10✓ Branch 0 taken 13022134 times.
✓ Branch 1 taken 8750 times.
✓ Branch 2 taken 40254 times.
✓ Branch 3 taken 12981880 times.
✓ Branch 4 taken 39820 times.
✓ Branch 5 taken 434 times.
✓ Branch 6 taken 4471146 times.
✓ Branch 7 taken 8550554 times.
✓ Branch 8 taken 4444606 times.
✓ Branch 9 taken 26540 times.
|
17475924 | dabs(x(2) - xlast(2)) < si%rtol*tolref(2))) |
| 111 | end subroutine sympl_euler1_newton_iter | ||
| 112 | |||
| 113 | 4445040 | subroutine sympl_euler1_extrapolate_field(si, f, x, xlast) | |
| 114 | !$acc routine seq | ||
| 115 | type(symplectic_integrator_t), intent(in) :: si | ||
| 116 | type(field_can_t), intent(inout) :: f | ||
| 117 | real(dp), intent(in) :: x(2) | ||
| 118 | real(dp), intent(in) :: xlast(2) | ||
| 119 | |||
| 120 | 4445040 | f%pth = f%pth + f%dpth(1)*(x(1) - xlast(1)) + f%dpth(4)*(x(2) - xlast(2)) | |
| 121 | 4445040 | f%dH(1) = f%dH(1) + f%d2H(1)*(x(1) - xlast(1)) + f%d2H(7)*(x(2) - xlast(2)) | |
| 122 | f%dpth(1) = f%dpth(1) + f%d2pth(1)*(x(1) - xlast(1)) & | ||
| 123 | 4445040 | + f%d2pth(7)*(x(2) - xlast(2)) | |
| 124 | 4445040 | f%vpar = f%vpar + f%dvpar(1)*(x(1) - xlast(1)) + f%dvpar(4)*(x(2) - xlast(2)) | |
| 125 | 4445040 | f%hth = f%hth + f%dhth(1)*(x(1) - xlast(1)) | |
| 126 | 4445040 | f%hph = f%hph + f%dhph(1)*(x(1) - xlast(1)) | |
| 127 | 4445040 | end subroutine sympl_euler1_extrapolate_field | |
| 128 | |||
| 129 | 4445040 | subroutine sympl_euler1_advance_angles(si, f) | |
| 130 | !$acc routine seq | ||
| 131 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 132 | type(field_can_t), intent(in) :: f | ||
| 133 | |||
| 134 | 4445040 | si%z(2) = si%z(2) + si%dt*f%dH(1)/f%dpth(1) | |
| 135 | 4445040 | si%z(3) = si%z(3) + si%dt*(f%vpar - f%dH(1)/f%dpth(1)*f%hth)/f%hph | |
| 136 | 4445040 | end subroutine sympl_euler1_advance_angles | |
| 137 | |||
| 138 | end module orbit_symplectic_euler1 | ||
| 139 |