| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module orbit_rk_core | ||
| 2 | ! Device-callable shared core for the implicit Gauss-collocation step used by | ||
| 3 | ! the guiding-center (GC) integrator. No procedure pointers and no class() | ||
| 4 | ! polymorphism in the hot path: the residual/Jacobian/Newton kernels are | ||
| 5 | ! inlinable !$acc routine seq routines, the way orbit_symplectic_euler1 | ||
| 6 | ! already shares the symplectic-Euler residual/Jacobian/Newton between the | ||
| 7 | ! CPU integrator and the OpenACC GPU kernel. | ||
| 8 | ! | ||
| 9 | ! The arithmetic is lifted byte-identically from the former f_rk_gauss / | ||
| 10 | ! jac_rk_gauss bodies in orbit_symplectic. coeff_rk_gauss stays in | ||
| 11 | ! orbit_symplectic_base (plain data) and is called inside the kernels. The | ||
| 12 | ! Newton shell here uses the LU solve rk_solve, which the CPU symplectic | ||
| 13 | ! integrators in orbit_symplectic also use (no LAPACK call per step). | ||
| 14 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 15 | use util, only: twopi | ||
| 16 | use field_can_mod, only: field_can_t, get_derivatives2, eval_field => evaluate | ||
| 17 | use orbit_symplectic_base, only: symplectic_integrator_t, coeff_rk_gauss, sympl_rmax | ||
| 18 | |||
| 19 | implicit none | ||
| 20 | private | ||
| 21 | |||
| 22 | ! Model code kept as a named constant for the residual/Jacobian/Newton | ||
| 23 | ! signatures so the GC caller in orbit_symplectic reads explicitly. Only the | ||
| 24 | ! GC kernel exists here; the CPP and symplectic full-orbit models were shelved. | ||
| 25 | integer, parameter, public :: MODEL_GC = 0 | ||
| 26 | |||
| 27 | public :: rk_gauss_residual, rk_gauss_jacobian, rk_gauss_newton, rk_solve | ||
| 28 | |||
| 29 | contains | ||
| 30 | |||
| 31 | ! GC Gauss residual: the field-canonical degenerate-Lagrangian residual. | ||
| 32 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2439 times.
|
2439 | subroutine rk_gauss_residual(model, si, fs, s, x, fvec) |
| 33 | !$acc routine seq | ||
| 34 | integer, intent(in) :: model | ||
| 35 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 36 | integer, intent(in) :: s | ||
| 37 | type(field_can_t), intent(inout) :: fs(:) | ||
| 38 | real(dp), intent(in) :: x(4*s) | ||
| 39 | real(dp), intent(out) :: fvec(4*s) | ||
| 40 | |||
| 41 | 2439 | call gauss_canfield_residual(si, fs, s, x, fvec) | |
| 42 | 2439 | end subroutine rk_gauss_residual | |
| 43 | |||
| 44 | ! GC Gauss Jacobian, analytic from the 2nd derivatives of field_can_t. | ||
| 45 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2439 times.
|
2439 | subroutine rk_gauss_jacobian(model, si, fs, s, jac) |
| 46 | !$acc routine seq | ||
| 47 | integer, intent(in) :: model | ||
| 48 | type(symplectic_integrator_t), intent(in) :: si | ||
| 49 | integer, intent(in) :: s | ||
| 50 | type(field_can_t), intent(in) :: fs(:) | ||
| 51 | real(dp), intent(out) :: jac(4*s, 4*s) | ||
| 52 | |||
| 53 | 2439 | call gauss_canfield_jacobian(si, fs, s, jac) | |
| 54 | 2439 | end subroutine rk_gauss_jacobian | |
| 55 | |||
| 56 | ! Concrete residual kernel: layout x(4*k-3:4*k) = (r,theta,phi,pphi) per stage. | ||
| 57 | ! Body lifted verbatim from f_rk_gauss. | ||
| 58 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2439 times.
|
2439 | subroutine gauss_canfield_residual(si, fs, s, x, fvec) |
| 59 | !$acc routine seq | ||
| 60 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 61 | integer, intent(in) :: s | ||
| 62 | type(field_can_t), intent(inout) :: fs(:) | ||
| 63 | real(dp), intent(in) :: x(4*s) | ||
| 64 | real(dp), intent(out) :: fvec(4*s) | ||
| 65 | |||
| 66 | 4878 | real(dp) :: a(s,s), b(s), c(s), Hprime(s) | |
| 67 | integer :: k, l | ||
| 68 | |||
| 69 | 2439 | call coeff_rk_gauss(s, a, b, c) | |
| 70 | |||
| 71 |
2/2✓ Branch 1 taken 2439 times.
✓ Branch 2 taken 7224 times.
|
9663 | do k = 1, s |
| 72 | 7224 | call eval_field(fs(k), x(4*k-3), x(4*k-2), x(4*k-1), 2) | |
| 73 | 7224 | call get_derivatives2(fs(k), x(4*k)) | |
| 74 | 9663 | Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) | |
| 75 | end do | ||
| 76 | |||
| 77 |
2/2✓ Branch 0 taken 7224 times.
✓ Branch 1 taken 2439 times.
|
9663 | do k = 1, s |
| 78 | 7224 | fvec(4*k-3) = fs(k)%pth - si%pthold | |
| 79 | 7224 | fvec(4*k-2) = x(4*k-2) - si%z(2) | |
| 80 | 7224 | fvec(4*k-1) = x(4*k-1) - si%z(3) | |
| 81 | 7224 | fvec(4*k) = x(4*k) - si%z(4) | |
| 82 |
2/2✓ Branch 0 taken 24676 times.
✓ Branch 1 taken 7224 times.
|
34339 | do l = 1, s |
| 83 | 24676 | fvec(4*k-3) = fvec(4*k-3) + si%dt*a(k,l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) ! pthdot | |
| 84 | 24676 | fvec(4*k-2) = fvec(4*k-2) - si%dt*a(k,l)*Hprime(l) ! thdot | |
| 85 | 24676 | fvec(4*k-1) = fvec(4*k-1) - si%dt*a(k,l)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph ! phdot | |
| 86 | 31900 | fvec(4*k) = fvec(4*k) + si%dt*a(k,l)*(fs(l)%dH(3) - Hprime(l)*fs(l)%dpth(3)) ! pphdot | |
| 87 | end do | ||
| 88 | end do | ||
| 89 | 2439 | end subroutine gauss_canfield_residual | |
| 90 | |||
| 91 | ! Concrete Jacobian kernel. Body lifted verbatim from jac_rk_gauss. | ||
| 92 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2439 times.
|
2439 | subroutine gauss_canfield_jacobian(si, fs, s, jac) |
| 93 | !$acc routine seq | ||
| 94 | type(symplectic_integrator_t), intent(in) :: si | ||
| 95 | integer, intent(in) :: s | ||
| 96 | type(field_can_t), intent(in) :: fs(:) | ||
| 97 | real(dp), intent(out) :: jac(4*s, 4*s) | ||
| 98 | |||
| 99 | 4878 | real(dp) :: a(s,s), b(s), c(s), Hprime(s), dHprime(4*s) | |
| 100 | integer :: k, l, m | ||
| 101 | |||
| 102 | 2439 | call coeff_rk_gauss(s, a, b, c) | |
| 103 | |||
| 104 |
2/2✓ Branch 1 taken 2439 times.
✓ Branch 2 taken 7224 times.
|
9663 | do k = 1, s |
| 105 | 7224 | m = 4*k | |
| 106 | 7224 | Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) | |
| 107 | 7224 | dHprime(m-3) = (fs(k)%d2H(1) - Hprime(k)*fs(k)%d2pth(1))/fs(k)%dpth(1) ! d/dr | |
| 108 | 7224 | dHprime(m-2) = (fs(k)%d2H(2) - Hprime(k)*fs(k)%d2pth(2))/fs(k)%dpth(1) ! d/dth | |
| 109 | 7224 | dHprime(m-1) = (fs(k)%d2H(3) - Hprime(k)*fs(k)%d2pth(3))/fs(k)%dpth(1) ! d/dph | |
| 110 | 9663 | dHprime(m) = (fs(k)%d2H(7) - Hprime(k)*fs(k)%d2pth(7))/fs(k)%dpth(1) ! d/dpph | |
| 111 | end do | ||
| 112 | |||
| 113 |
4/4✓ Branch 0 taken 2439 times.
✓ Branch 1 taken 28896 times.
✓ Branch 2 taken 394816 times.
✓ Branch 3 taken 28896 times.
|
426151 | jac = 0d0 |
| 114 | |||
| 115 |
2/2✓ Branch 0 taken 2439 times.
✓ Branch 1 taken 7224 times.
|
9663 | do k = 1, s |
| 116 | 7224 | m = 4*k | |
| 117 | 7224 | jac(m-3, m-3) = fs(k)%dpth(1) | |
| 118 | 7224 | jac(m-3, m-2) = fs(k)%dpth(2) | |
| 119 | 7224 | jac(m-3, m-1) = fs(k)%dpth(3) | |
| 120 | 7224 | jac(m-3, m) = fs(k)%dpth(4) | |
| 121 | 7224 | jac(m-2, m-2) = 1d0 | |
| 122 | 7224 | jac(m-1, m-1) = 1d0 | |
| 123 | 9663 | jac(m, m) = 1d0 | |
| 124 | end do | ||
| 125 | |||
| 126 |
2/2✓ Branch 0 taken 7224 times.
✓ Branch 1 taken 2439 times.
|
9663 | do l = 1, s |
| 127 |
2/2✓ Branch 0 taken 24676 times.
✓ Branch 1 taken 7224 times.
|
34339 | do k = 1, s |
| 128 | 24676 | m = 4*k | |
| 129 | jac(m-3, 4*l-3) = jac(m-3, 4*l-3) & ! d/dr | ||
| 130 | 24676 | + si%dt*a(k,l)*(fs(l)%d2H(2) - fs(l)%d2pth(2)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l-3)) | |
| 131 | jac(m-3, 4*l-2) = jac(m-3, 4*l-2) & ! d/dth | ||
| 132 | 24676 | + si%dt*a(k,l)*(fs(l)%d2H(4) - fs(l)%d2pth(4)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l-2)) | |
| 133 | jac(m-3, 4*l-1) = jac(m-3, 4*l-1) & ! d/dph | ||
| 134 | 24676 | + si%dt*a(k,l)*(fs(l)%d2H(5) - fs(l)%d2pth(5)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l-1)) | |
| 135 | jac(m-3, 4*l) = jac(m-3, 4*l) & ! d/dpph | ||
| 136 | 24676 | + si%dt*a(k,l)*(fs(l)%d2H(8) - fs(l)%d2pth(8)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l)) | |
| 137 | |||
| 138 | 24676 | jac(m-2, 4*l-3) = jac(m-2, 4*l-3) - si%dt*a(k,l)*dHprime(4*l-3) ! d/dr | |
| 139 | 24676 | jac(m-2, 4*l-2) = jac(m-2, 4*l-2) - si%dt*a(k,l)*dHprime(4*l-2) ! d/dth | |
| 140 | 24676 | jac(m-2, 4*l-1) = jac(m-2, 4*l-1) - si%dt*a(k,l)*dHprime(4*l-1) ! d/dph | |
| 141 | 24676 | jac(m-2, 4*l) = jac(m-2, 4*l) - si%dt*a(k,l)*dHprime(4*l) ! d/dpph | |
| 142 | |||
| 143 | jac(m-1, 4*l-3) = jac(m-1, 4*l-3) & ! d/dr | ||
| 144 | - si%dt*a(k,l)*(-fs(l)%dhph(1)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & | ||
| 145 | 24676 | + (fs(l)%dvpar(1) - dHprime(4*l-3)*fs(l)%hth - Hprime(l)*fs(l)%dhth(1))/fs(l)%hph) | |
| 146 | jac(m-1, 4*l-2) = jac(m-1, 4*l-2) & ! d/dth | ||
| 147 | - si%dt*a(k,l)*(-fs(l)%dhph(2)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & | ||
| 148 | 24676 | + (fs(l)%dvpar(2) - dHprime(4*l-2)*fs(l)%hth - Hprime(l)*fs(l)%dhth(2))/fs(l)%hph) | |
| 149 | jac(m-1, 4*l-1) = jac(m-1, 4*l-1) & ! d/dph | ||
| 150 | - si%dt*a(k,l)*(-fs(l)%dhph(3)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & | ||
| 151 | 24676 | + (fs(l)%dvpar(3) - dHprime(4*l-1)*fs(l)%hth - Hprime(l)*fs(l)%dhth(3))/fs(l)%hph) | |
| 152 | jac(m-1, 4*l) = jac(m-1, 4*l) & ! d/dpph | ||
| 153 | 24676 | - si%dt*a(k,l)*((fs(l)%dvpar(4) - dHprime(4*l)*fs(l)%hth)/fs(l)%hph) | |
| 154 | |||
| 155 | jac(m, 4*l-3) = jac(m, 4*l-3) & ! d/dr | ||
| 156 | 24676 | + si%dt*a(k,l)*(fs(l)%d2H(3) - fs(l)%d2pth(3)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l-3)) | |
| 157 | jac(m, 4*l-2) = jac(m, 4*l-2) & ! d/dth | ||
| 158 | 24676 | + si%dt*a(k,l)*(fs(l)%d2H(5) - fs(l)%d2pth(5)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l-2)) | |
| 159 | jac(m, 4*l-1) = jac(m, 4*l-1) & ! d/dph | ||
| 160 | 24676 | + si%dt*a(k,l)*(fs(l)%d2H(6) - fs(l)%d2pth(6)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l-1)) | |
| 161 | jac(m, 4*l) = jac(m, 4*l) & ! d/dpph | ||
| 162 | 31900 | + si%dt*a(k,l)*(fs(l)%d2H(9) - fs(l)%d2pth(9)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l)) | |
| 163 | end do | ||
| 164 | end do | ||
| 165 | 2439 | end subroutine gauss_canfield_jacobian | |
| 166 | |||
| 167 | ! Device-portable dense LU solve A x = rhs with partial pivoting, in place on | ||
| 168 | ! rhs. n <= 4*S_MAX (n <= 16 for s <= 4). info = 0 on success, else the failing | ||
| 169 | ! pivot column. Replaces dgesv on CPU and device: thread-safe, no allocation. | ||
| 170 | 3421227 | pure subroutine rk_solve(n, A, rhs, info) | |
| 171 | !$acc routine seq | ||
| 172 | integer, intent(in) :: n | ||
| 173 | real(dp), intent(inout) :: A(n,n), rhs(n) | ||
| 174 | integer, intent(out) :: info | ||
| 175 | integer :: i, j, k, ipiv | ||
| 176 | real(dp) :: piv, amax, factor, tmp | ||
| 177 | |||
| 178 | 3421227 | info = 0 | |
| 179 |
2/2✓ Branch 0 taken 17137126 times.
✓ Branch 1 taken 3421226 times.
|
20558352 | do k = 1, n |
| 180 | 17137126 | ipiv = k | |
| 181 | 17137126 | amax = abs(A(k,k)) | |
| 182 |
2/2✓ Branch 0 taken 34470908 times.
✓ Branch 1 taken 17137126 times.
|
51608034 | do i = k+1, n |
| 183 |
2/2✓ Branch 0 taken 2915978 times.
✓ Branch 1 taken 31554930 times.
|
51608034 | if (abs(A(i,k)) > amax) then |
| 184 | 2915978 | amax = abs(A(i,k)) | |
| 185 | 2915978 | ipiv = i | |
| 186 | end if | ||
| 187 | end do | ||
| 188 | ! Singular column: amax = max|A(i,k)| >= 0, so .not.(amax > 0) means the | ||
| 189 | ! whole pivot column is zero. Phrased this way to avoid a real-equality | ||
| 190 | ! comparison (-Wcompare-reals). | ||
| 191 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 17137125 times.
|
17137126 | if (.not. (amax > 0d0)) then |
| 192 | 1 | info = k | |
| 193 | 1 | return | |
| 194 | end if | ||
| 195 |
2/2✓ Branch 0 taken 1762567 times.
✓ Branch 1 taken 15374558 times.
|
17137125 | if (ipiv /= k) then |
| 196 |
2/2✓ Branch 0 taken 8948657 times.
✓ Branch 1 taken 1762567 times.
|
10711224 | do j = 1, n |
| 197 | 10711224 | tmp = A(k,j); A(k,j) = A(ipiv,j); A(ipiv,j) = tmp | |
| 198 | end do | ||
| 199 | 1762567 | tmp = rhs(k); rhs(k) = rhs(ipiv); rhs(ipiv) = tmp | |
| 200 | end if | ||
| 201 | 17137125 | piv = A(k,k) | |
| 202 |
2/2✓ Branch 0 taken 17137125 times.
✓ Branch 1 taken 34470907 times.
|
55029258 | do i = k+1, n |
| 203 | 34470907 | factor = A(i,k)/piv | |
| 204 | 34470907 | A(i,k) = factor | |
| 205 |
2/2✓ Branch 0 taken 105025887 times.
✓ Branch 1 taken 34470907 times.
|
139496794 | do j = k+1, n |
| 206 | 139496794 | A(i,j) = A(i,j) - factor*A(k,j) | |
| 207 | end do | ||
| 208 | 51608032 | rhs(i) = rhs(i) - factor*rhs(k) | |
| 209 | end do | ||
| 210 | end do | ||
| 211 | |||
| 212 |
2/2✓ Branch 0 taken 17137125 times.
✓ Branch 1 taken 3421226 times.
|
20558351 | do i = n, 1, -1 |
| 213 | 17137125 | tmp = rhs(i) | |
| 214 |
2/2✓ Branch 0 taken 34470907 times.
✓ Branch 1 taken 17137125 times.
|
51608032 | do j = i+1, n |
| 215 | 51608032 | tmp = tmp - A(i,j)*rhs(j) | |
| 216 | end do | ||
| 217 | 20558351 | rhs(i) = tmp/A(i,i) | |
| 218 | end do | ||
| 219 | end subroutine rk_solve | ||
| 220 | |||
| 221 | ! Device-callable Newton iteration for the Gauss step. Mirrors the | ||
| 222 | ! newton_rk_gauss control flow (atol/rtol/tolref, boundary guards, maxit) with | ||
| 223 | ! the device LU solver rk_solve in place of dgesv. No event counters here: | ||
| 224 | ! hit_maxit is returned so the CPU caller can record EVT_RK_GAUSS_MAXIT, | ||
| 225 | ! exactly as the GPU newton1 path omits fort.6601. Retained for the deferred | ||
| 226 | ! fortnum multiroot wiring; not yet called by the CPU GC Newton step. | ||
| 227 | ✗ | subroutine rk_gauss_newton(model, si, fs, s, x, atol, rtol, maxit, xlast, & | |
| 228 | hit_maxit) | ||
| 229 | !$acc routine seq | ||
| 230 | integer, intent(in) :: model | ||
| 231 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 232 | integer, intent(in) :: s | ||
| 233 | type(field_can_t), intent(inout) :: fs(:) | ||
| 234 | real(dp), intent(inout) :: x(4*s) | ||
| 235 | real(dp), intent(in) :: atol, rtol | ||
| 236 | integer, intent(in) :: maxit | ||
| 237 | real(dp), intent(out) :: xlast(4*s) | ||
| 238 | logical, intent(out) :: hit_maxit | ||
| 239 | |||
| 240 | ✗ | real(dp) :: fvec(4*s), fjac(4*s, 4*s) | |
| 241 | ✗ | real(dp) :: xabs(4*s), tolref(4*s), fabs(4*s) | |
| 242 | integer :: kit, ks, info | ||
| 243 | logical :: conv | ||
| 244 | |||
| 245 | ✗ | hit_maxit = .false. | |
| 246 | ✗ | do kit = 1, maxit | |
| 247 | |||
| 248 | ! Check if radius left the boundary | ||
| 249 | ✗ | do ks = 1, s | |
| 250 | ✗ | if (x(4*ks-3) > sympl_rmax) return | |
| 251 | ! Transient guard for intermediate iterates; the converged-negative | ||
| 252 | ! case is handled by the caller via a chart switch (#370). | ||
| 253 | ✗ | if (x(4*ks-3) < 0.0d0) x(4*ks-3) = 0.01d0 | |
| 254 | end do | ||
| 255 | |||
| 256 | ✗ | call rk_gauss_residual(model, si, fs, s, x, fvec) | |
| 257 | ✗ | call rk_gauss_jacobian(model, si, fs, s, fjac) | |
| 258 | ✗ | fabs = abs(fvec) | |
| 259 | ✗ | xlast = x | |
| 260 | ✗ | call rk_solve(4*s, fjac, fvec, info) | |
| 261 | ✗ | if (info /= 0) return | |
| 262 | ! after solution: fvec = (xold-xnew)_Newton | ||
| 263 | ✗ | x = x - fvec | |
| 264 | ✗ | xabs = abs(x - xlast) | |
| 265 | |||
| 266 | ✗ | do ks = 1, s | |
| 267 | ✗ | tolref(4*ks-3) = 1d0 | |
| 268 | ✗ | tolref(4*ks-2) = twopi | |
| 269 | ✗ | tolref(4*ks-1) = twopi | |
| 270 | ✗ | tolref(4*ks) = abs(xlast(4*ks)) | |
| 271 | end do | ||
| 272 | |||
| 273 | conv = .true. | ||
| 274 | ✗ | do ks = 1, 4*s | |
| 275 | ✗ | if (fabs(ks) >= atol) conv = .false. | |
| 276 | end do | ||
| 277 | ✗ | if (conv) return | |
| 278 | conv = .true. | ||
| 279 | ✗ | do ks = 1, 4*s | |
| 280 | ✗ | if (xabs(ks) >= rtol*tolref(ks)) conv = .false. | |
| 281 | end do | ||
| 282 | ✗ | if (conv) return | |
| 283 | end do | ||
| 284 | ✗ | hit_maxit = .true. | |
| 285 | end subroutine rk_gauss_newton | ||
| 286 | |||
| 287 | end module orbit_rk_core | ||
| 288 |