| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module spectre_fo_hybrid | ||
| 2 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 3 | use magfie_sub, only: spectre_field | ||
| 4 | use spectre_fo_field, only: spectre_fo_field_build, spectre_fo_field_local, & | ||
| 5 | spectre_fo_field_locate, spectre_fo_local_t, fo_field_axisymmetric, & | ||
| 6 | FO_AXIS_RADIUS | ||
| 7 | implicit none | ||
| 8 | private | ||
| 9 | integer, parameter, public :: SPECTRE_FO_OK = 0 | ||
| 10 | integer, parameter, public :: SPECTRE_FO_UNINITIALIZED = 1 | ||
| 11 | integer, parameter, public :: SPECTRE_FO_LOCATE_FAIL = 2 | ||
| 12 | integer, parameter, public :: SPECTRE_FO_GYROPHASE_FAIL = 3 | ||
| 13 | integer, parameter, public :: SPECTRE_FO_INVALID_STATE = 4 | ||
| 14 | integer, parameter, public :: SPECTRE_FO_STEP_LIMIT = 5 | ||
| 15 | integer, parameter, public :: SPECTRE_FO_LOSS = 6 | ||
| 16 | real(dp), parameter :: OWNER_EPS = 1.0e-12_dp | ||
| 17 | integer, parameter :: STAGE_IN = 0, STAGE_OUT = 1 | ||
| 18 | !> RK45 relative tolerance, configured before parallel tracing. | ||
| 19 | real(dp) :: fo_tolerance = 1.0e-9_dp | ||
| 20 | ! At most 0.1 rad (1/63 of a gyro-period). A cross-and-return hidden | ||
| 21 | ! between stage samples has circular-orbit depth <= rho_L*(1-cos(0.05)), | ||
| 22 | ! about 1.25e-3 rho_L; tolerance alone cannot bound a missed sheet force. | ||
| 23 | real(dp), parameter :: RK_GYRO_ANGLE = 0.1_dp | ||
| 24 | !> Gyroangle between exit checks. The guiding-centre map costs about as much | ||
| 25 | !> as an RK45 substep, so it is not evaluated after every substep; an exit | ||
| 26 | !> is detected at most this angle late, and always at the end of a call. | ||
| 27 | real(dp), parameter :: EXIT_CHECK_ANGLE = 1.0_dp | ||
| 28 | real(dp), parameter :: ENTRY_MOMENTUM_TOL = 1.0e-3_dp | ||
| 29 | integer, parameter :: MAX_SUBSTEPS = 100000 | ||
| 30 | integer, parameter :: NPHASE = 64 | ||
| 31 | type, public :: spectre_fo_state_t | ||
| 32 | real(dp) :: x(3) = 0.0_dp | ||
| 33 | real(dp) :: v(3) = 0.0_dp | ||
| 34 | real(dp) :: u(3) = 0.0_dp | ||
| 35 | real(dp) :: last_y(5) = 0.0_dp | ||
| 36 | real(dp) :: entry_rho = 0.0_dp | ||
| 37 | real(dp) :: entry_pzeta = 0.0_dp | ||
| 38 | ! Actual FO minus target GC momentum, expressed in the entry gauge. | ||
| 39 | real(dp) :: entry_pzeta_defect = 0.0_dp | ||
| 40 | !> Next RK45 step (adapted between substeps). | ||
| 41 | real(dp) :: h_rk = 0.0_dp | ||
| 42 | !> Volume whose gauge entry_pzeta is expressed in. | ||
| 43 | integer :: entry_owner = 0 | ||
| 44 | !> Exit when the guiding centre is farther than this from the interface | ||
| 45 | !> (rho_g); 0 keeps the gyrophase-based legacy rule. | ||
| 46 | real(dp) :: exit_rho = 0.0_dp | ||
| 47 | integer :: owner = 0 | ||
| 48 | integer :: iface = 0 | ||
| 49 | logical :: has_y = .false. | ||
| 50 | logical :: active = .false. | ||
| 51 | end type spectre_fo_state_t | ||
| 52 | public :: spectre_fo_enter, spectre_fo_to_gc, spectre_fo_configure | ||
| 53 | public :: spectre_fo_radial_larmor | ||
| 54 | public :: spectre_fo_advance_until_exit | ||
| 55 | public :: spectre_fo_canonical_pzeta | ||
| 56 | contains | ||
| 57 | 31 | subroutine spectre_fo_configure(tolerance, grid) | |
| 58 | real(dp), intent(in) :: tolerance | ||
| 59 | integer, intent(in), optional :: grid(3) | ||
| 60 | integer :: n(3) | ||
| 61 | |||
| 62 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 31 times.
|
31 | if (tolerance <= 0.0_dp) & |
| 63 | ✗ | error stop 'spectre_fo_configure: tolerance must be positive' | |
| 64 | 31 | fo_tolerance = tolerance | |
| 65 | 31 | n = [32, 64, 32] | |
| 66 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 7 times.
|
31 | if (present(grid)) n = grid |
| 67 |
6/6✓ Branch 0 taken 119 times.
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 107 times.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 19 times.
✓ Branch 5 taken 12 times.
|
138 | if (all(spectre_field%data%in == 0)) n(3) = 7 |
| 68 | 31 | call spectre_fo_field_build(n(1), n(2), n(3)) | |
| 69 | 31 | end subroutine spectre_fo_configure | |
| 70 | 254 | subroutine spectre_fo_enter(state, y, owner, iface, ro0_bar, status) | |
| 71 | type(spectre_fo_state_t), intent(out) :: state | ||
| 72 | real(dp), intent(in) :: y(5), ro0_bar | ||
| 73 | integer, intent(in) :: owner, iface | ||
| 74 | integer, intent(out) :: status | ||
| 75 | real(dp) :: x_gc(3), b(3), e1(3), e2(3), eperp(3), e_zeta(3) | ||
| 76 | real(dp) :: acov(3), bmod, vpar, vperp, target, alpha | ||
| 77 | integer :: locate_status | ||
| 78 |
8/8✓ Branch 0 taken 762 times.
✓ Branch 1 taken 254 times.
✓ Branch 2 taken 762 times.
✓ Branch 3 taken 254 times.
✓ Branch 4 taken 762 times.
✓ Branch 5 taken 254 times.
✓ Branch 6 taken 1270 times.
✓ Branch 7 taken 254 times.
|
3810 | state = spectre_fo_state_t() |
| 79 |
2/4✓ Branch 0 taken 254 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 254 times.
|
254 | if (ro0_bar <= 0.0_dp .or. y(4) <= 0.0_dp) then |
| 80 | ✗ | status = SPECTRE_FO_INVALID_STATE | |
| 81 | 15 | return | |
| 82 | end if | ||
| 83 |
2/2✓ Branch 0 taken 762 times.
✓ Branch 1 taken 254 times.
|
1016 | state%u = y(1:3) |
| 84 |
2/2✓ Branch 0 taken 1270 times.
✓ Branch 1 taken 254 times.
|
1524 | state%last_y = y |
| 85 | 254 | state%entry_rho = y(1) | |
| 86 | 254 | state%has_y = .true. | |
| 87 | 254 | state%owner = owner | |
| 88 | 254 | state%iface = iface | |
| 89 | 254 | call logical_sample(state%u, owner, x_gc, b, bmod, acov, e_zeta, status) | |
| 90 |
1/2✓ Branch 0 taken 254 times.
✗ Branch 1 not taken.
|
254 | if (status /= SPECTRE_FO_OK) return |
| 91 | 254 | vpar = sqrt(2.0_dp)*y(4)*y(5) | |
| 92 | 254 | vperp = sqrt(max(2.0_dp*y(4)**2 - vpar**2, 0.0_dp)) | |
| 93 | 254 | call perpendicular_frame(b, e1, e2) | |
| 94 |
2/2✓ Branch 0 taken 762 times.
✓ Branch 1 taken 254 times.
|
1016 | target = vpar*dot_product(b, e_zeta) + acov(3)/ro0_bar |
| 95 | 254 | state%entry_pzeta = target | |
| 96 | 254 | state%entry_owner = owner | |
| 97 | 254 | alpha = 0.0_dp | |
| 98 | 507 | if (axisymmetric()) then | |
| 99 | call solve_gyrophase(x_gc, state%u, owner, b, e1, e2, vpar, & | ||
| 100 | 42 | vperp, target, ro0_bar, alpha, status) | |
| 101 |
2/2✓ Branch 0 taken 41 times.
✓ Branch 1 taken 1 time.
|
42 | if (status /= SPECTRE_FO_OK) return |
| 102 | end if | ||
| 103 |
2/2✓ Branch 0 taken 759 times.
✓ Branch 1 taken 253 times.
|
1012 | eperp = cos(alpha)*e1 + sin(alpha)*e2 |
| 104 |
2/2✓ Branch 0 taken 759 times.
✓ Branch 1 taken 253 times.
|
1012 | state%v = vpar*b + vperp*eperp |
| 105 |
6/6✓ Branch 0 taken 759 times.
✓ Branch 1 taken 253 times.
✓ Branch 3 taken 759 times.
✓ Branch 4 taken 253 times.
✓ Branch 5 taken 759 times.
✓ Branch 6 taken 253 times.
|
2530 | state%x = x_gc + (ro0_bar/bmod)*cross(b, vperp*eperp) |
| 106 | 253 | call spectre_fo_field_locate(state%x, state%u, state%owner, locate_status) | |
| 107 |
2/2✓ Branch 0 taken 14 times.
✓ Branch 1 taken 239 times.
|
253 | if (locate_status /= 0) then |
| 108 | 14 | status = SPECTRE_FO_LOCATE_FAIL | |
| 109 |
2/4✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 14 times.
✗ Branch 3 not taken.
|
14 | if (locate_status == 1 .and. & |
| 110 | state%u(1) > real(spectre_field%data%Mvol, dp)) then | ||
| 111 | 14 | state%active = .true. | |
| 112 | 14 | status = SPECTRE_FO_LOSS | |
| 113 | end if | ||
| 114 | 14 | return | |
| 115 | end if | ||
| 116 | 478 | if (axisymmetric()) then | |
| 117 | call phase_residual(alpha, x_gc, y(1:3), owner, b, e1, e2, vpar, & | ||
| 118 | 41 | vperp, target, ro0_bar, state%entry_pzeta_defect, status) | |
| 119 |
1/2✓ Branch 0 taken 41 times.
✗ Branch 1 not taken.
|
41 | if (status /= SPECTRE_FO_OK) return |
| 120 | end if | ||
| 121 | 239 | state%active = .true. | |
| 122 | end subroutine spectre_fo_enter | ||
| 123 | 10356 | subroutine spectre_fo_radial_larmor(y, owner, ro0_bar, rho_r, status) | |
| 124 | !> Gyroradius of the guiding-centre state y in rho_g units. | ||
| 125 | real(dp), intent(in) :: y(5), ro0_bar | ||
| 126 | integer, intent(in) :: owner | ||
| 127 | real(dp), intent(out) :: rho_r | ||
| 128 | integer, intent(out) :: status | ||
| 129 | real(dp) :: x(3), b(3), bmod, acov(3), e_zeta(3), vpar, vperp | ||
| 130 | |||
| 131 | 10356 | rho_r = 0.0_dp | |
| 132 | 10356 | call logical_sample(y(1:3), owner, x, b, bmod, acov, e_zeta, status) | |
| 133 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10356 times.
|
10356 | if (status /= SPECTRE_FO_OK) return |
| 134 | 10356 | vpar = sqrt(2.0_dp)*y(4)*y(5) | |
| 135 | 10356 | vperp = sqrt(max(2.0_dp*y(4)**2 - vpar**2, 0.0_dp)) | |
| 136 | 10356 | rho_r = ro0_bar*vperp/bmod*radial_gradient(y(1:3)) | |
| 137 | end subroutine spectre_fo_radial_larmor | ||
| 138 | 10356 | real(dp) function radial_gradient(u) result(grad) | |
| 139 | real(dp), intent(in) :: u(3) | ||
| 140 | type(spectre_fo_local_t) :: loc | ||
| 141 | |||
| 142 | call spectre_fo_field_local(u, min(int(u(1)) + 1, & | ||
| 143 | 10356 | spectre_field%data%Mvol), loc) | |
| 144 |
4/6✓ Branch 0 taken 31068 times.
✓ Branch 1 taken 10356 times.
✓ Branch 2 taken 31068 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 31068 times.
|
41424 | grad = norm2(loc%dual(:, 1)) |
| 145 |
14/14✓ Branch 0 taken 31068 times.
✓ Branch 1 taken 10356 times.
✓ Branch 2 taken 31068 times.
✓ Branch 3 taken 10356 times.
✓ Branch 4 taken 93204 times.
✓ Branch 5 taken 31068 times.
✓ Branch 6 taken 31068 times.
✓ Branch 7 taken 10356 times.
✓ Branch 8 taken 93204 times.
✓ Branch 9 taken 31068 times.
✓ Branch 10 taken 31068 times.
✓ Branch 11 taken 10356 times.
✓ Branch 12 taken 31068 times.
✓ Branch 13 taken 10356 times.
|
352104 | end function radial_gradient |
| 146 | 1821 | subroutine spectre_fo_advance_until_exit(state, dt, ro0_bar, dt_used, y, & | |
| 147 | owner, exited, status, h_scale) | ||
| 148 | type(spectre_fo_state_t), intent(inout) :: state | ||
| 149 | real(dp), intent(in) :: dt, ro0_bar | ||
| 150 | real(dp), intent(out) :: dt_used, y(5) | ||
| 151 | integer, intent(out) :: owner, status | ||
| 152 | logical, intent(out) :: exited | ||
| 153 | !> Gyrostep reduction factor in (0, 1] for convergence studies of the | ||
| 154 | !> recorded canonical-momentum defect and the measured return mu; the | ||
| 155 | !> production path omits it. | ||
| 156 | real(dp), intent(in), optional :: h_scale | ||
| 157 | real(dp) :: remaining, h, hfac, angle | ||
| 158 | type(spectre_fo_local_t) :: particle | ||
| 159 | integer :: nstep, map_status | ||
| 160 | 1821 | dt_used = 0.0_dp | |
| 161 | 1821 | y = 0.0_dp | |
| 162 | 1821 | owner = state%owner | |
| 163 | 1821 | exited = .false. | |
| 164 |
2/4✓ Branch 0 taken 1821 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1821 times.
|
1821 | if (.not. state%active .or. dt < 0.0_dp .or. ro0_bar <= 0.0_dp) then |
| 165 | ✗ | status = SPECTRE_FO_INVALID_STATE | |
| 166 | 1821 | return | |
| 167 | end if | ||
| 168 | 1821 | hfac = 1.0_dp | |
| 169 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1821 times.
|
1821 | if (present(h_scale)) then |
| 170 | ✗ | if (h_scale <= 0.0_dp .or. h_scale > 1.0_dp) then | |
| 171 | ✗ | status = SPECTRE_FO_INVALID_STATE | |
| 172 | ✗ | return | |
| 173 | end if | ||
| 174 | hfac = h_scale | ||
| 175 | end if | ||
| 176 | ! Public states may be initialized from Cartesian particle data. | ||
| 177 | ! Synchronize once per call; accepted logical steps cache x exactly. | ||
| 178 | 1821 | call synchronize_position(state, status) | |
| 179 |
1/2✓ Branch 0 taken 1821 times.
✗ Branch 1 not taken.
|
1821 | if (status /= SPECTRE_FO_OK) return |
| 180 | 1821 | remaining = dt | |
| 181 | 1821 | angle = 0.0_dp | |
| 182 | 1821 | call spectre_fo_field_local(state%u, state%owner, particle) | |
| 183 |
1/2✓ Branch 1 taken 15899 times.
✗ Branch 2 not taken.
|
15899 | do nstep = 1, MAX_SUBSTEPS |
| 184 |
2/2✓ Branch 0 taken 1817 times.
✓ Branch 1 taken 14082 times.
|
15899 | if (remaining <= epsilon(1.0_dp)*max(dt, 1.0_dp)) then |
| 185 | 1817 | dt_used = dt | |
| 186 | 1817 | status = SPECTRE_FO_OK | |
| 187 | 1817 | return | |
| 188 | end if | ||
| 189 | 14082 | h = remaining | |
| 190 |
2/2✓ Branch 0 taken 14017 times.
✓ Branch 1 taken 65 times.
|
14082 | if (state%h_rk > 0.0_dp) h = min(remaining, hfac*state%h_rk) |
| 191 | 14082 | h = min(h, hfac*RK_GYRO_ANGLE*ro0_bar/particle%bmod) | |
| 192 | call rk45_substep(state, h, ro0_bar, status, hfac*RK_GYRO_ANGLE, & | ||
| 193 | 14082 | particle) | |
| 194 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 14078 times.
|
14082 | if (status /= SPECTRE_FO_OK) then |
| 195 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (status == SPECTRE_FO_LOSS) dt_used = dt_used + h |
| 196 | 4 | return | |
| 197 | end if | ||
| 198 | 14078 | dt_used = dt_used + h | |
| 199 | 14078 | remaining = remaining - h | |
| 200 | 14078 | angle = angle + h*particle%bmod/ro0_bar | |
| 201 |
4/4✓ Branch 0 taken 14011 times.
✓ Branch 1 taken 67 times.
✓ Branch 2 taken 1817 times.
✓ Branch 3 taken 12194 times.
|
14078 | if (angle < EXIT_CHECK_ANGLE .and. & |
| 202 | remaining > epsilon(1.0_dp)*max(dt, 1.0_dp)) cycle | ||
| 203 | 1884 | angle = 0.0_dp | |
| 204 | call gc_map(state, particle, ro0_bar, y, owner, map_status, & | ||
| 205 | 1884 | exited=exited) | |
| 206 |
2/2✓ Branch 1 taken 1882 times.
✓ Branch 2 taken 2 times.
|
1884 | if (map_status == SPECTRE_FO_OK) then |
| 207 |
2/2✓ Branch 0 taken 9410 times.
✓ Branch 1 taken 1882 times.
|
11292 | state%last_y = y |
| 208 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1882 times.
|
1882 | if (exited) then |
| 209 | ✗ | state%active = .false. | |
| 210 | ✗ | status = SPECTRE_FO_OK | |
| 211 | ✗ | return | |
| 212 | end if | ||
| 213 | end if | ||
| 214 | end do | ||
| 215 | ✗ | status = SPECTRE_FO_STEP_LIMIT | |
| 216 |
15/16✓ Branch 0 taken 5463 times.
✓ Branch 1 taken 1821 times.
✓ Branch 2 taken 5463 times.
✓ Branch 3 taken 1821 times.
✓ Branch 4 taken 16389 times.
✓ Branch 5 taken 5463 times.
✓ Branch 6 taken 5463 times.
✓ Branch 7 taken 1821 times.
✓ Branch 8 taken 16389 times.
✓ Branch 9 taken 5463 times.
✓ Branch 10 taken 5463 times.
✓ Branch 11 taken 1821 times.
✓ Branch 12 taken 5463 times.
✓ Branch 13 taken 1821 times.
✓ Branch 14 taken 1821 times.
✗ Branch 15 not taken.
|
61914 | end subroutine spectre_fo_advance_until_exit |
| 217 | 3888 | subroutine spectre_fo_to_gc(state, ro0_bar, y, owner, status, mu_bar) | |
| 218 | type(spectre_fo_state_t), intent(inout) :: state | ||
| 219 | real(dp), intent(in) :: ro0_bar | ||
| 220 | real(dp), intent(out) :: y(5) | ||
| 221 | integer, intent(out) :: owner, status | ||
| 222 | real(dp), intent(out), optional :: mu_bar | ||
| 223 | type(spectre_fo_local_t) :: particle | ||
| 224 | |||
| 225 | 3888 | y = 0.0_dp | |
| 226 | 3888 | owner = state%owner | |
| 227 | 3888 | if (present(mu_bar)) mu_bar = 0.0_dp | |
| 228 |
2/4✓ Branch 0 taken 3888 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3888 times.
|
3888 | if (.not. state%active .or. ro0_bar <= 0.0_dp) then |
| 229 | ✗ | status = SPECTRE_FO_INVALID_STATE | |
| 230 | ✗ | return | |
| 231 | end if | ||
| 232 | 3888 | call synchronize_position(state, status) | |
| 233 |
1/2✓ Branch 0 taken 3888 times.
✗ Branch 1 not taken.
|
3888 | if (status /= SPECTRE_FO_OK) return |
| 234 | 3888 | call spectre_fo_field_local(state%u, state%owner, particle) | |
| 235 | 3888 | call gc_map(state, particle, ro0_bar, y, owner, status, mu_bar) | |
| 236 |
15/16✓ Branch 0 taken 11664 times.
✓ Branch 1 taken 3888 times.
✓ Branch 2 taken 11664 times.
✓ Branch 3 taken 3888 times.
✓ Branch 4 taken 34992 times.
✓ Branch 5 taken 11664 times.
✓ Branch 6 taken 11664 times.
✓ Branch 7 taken 3888 times.
✓ Branch 8 taken 34992 times.
✓ Branch 9 taken 11664 times.
✓ Branch 10 taken 11664 times.
✓ Branch 11 taken 3888 times.
✓ Branch 12 taken 11664 times.
✓ Branch 13 taken 3888 times.
✓ Branch 14 taken 3888 times.
✗ Branch 15 not taken.
|
132192 | end subroutine spectre_fo_to_gc |
| 237 | |||
| 238 | 5772 | subroutine gc_map(state, particle, ro0_bar, y, owner, status, mu_bar, exited) | |
| 239 | type(spectre_fo_state_t), intent(in) :: state | ||
| 240 | type(spectre_fo_local_t), intent(in) :: particle | ||
| 241 | real(dp), intent(in) :: ro0_bar | ||
| 242 | real(dp), intent(out) :: y(5) | ||
| 243 | integer, intent(out) :: owner, status | ||
| 244 | real(dp), intent(out), optional :: mu_bar | ||
| 245 | logical, intent(out), optional :: exited | ||
| 246 | type(spectre_fo_local_t) :: gc | ||
| 247 | real(dp) :: b(3), b_gc(3), acov(3), e_zeta(3), bmod | ||
| 248 | real(dp) :: vpar, vperp(3), rho_l(3), x_gc(3), speed, acov_gc(3) | ||
| 249 | real(dp) :: e_zeta_gc(3), pzeta, h_zeta, perpendicular_energy, u_gc(3) | ||
| 250 | real(dp) :: radial_larmor | ||
| 251 | integer :: owner_gc | ||
| 252 | |||
| 253 | 1884 | if (present(exited)) exited = .false. | |
| 254 | 5772 | y = 0.0_dp | |
| 255 | 5772 | owner = state%owner | |
| 256 |
2/2✓ Branch 0 taken 3888 times.
✓ Branch 1 taken 1884 times.
|
5772 | if (present(mu_bar)) mu_bar = 0.0_dp |
| 257 |
2/4✓ Branch 0 taken 5772 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 5772 times.
|
5772 | if (.not. state%active .or. ro0_bar <= 0.0_dp) then |
| 258 | ✗ | status = SPECTRE_FO_INVALID_STATE | |
| 259 | 3890 | return | |
| 260 | end if | ||
| 261 |
2/2✓ Branch 0 taken 17316 times.
✓ Branch 1 taken 5772 times.
|
23088 | b = particle%b |
| 262 | 5772 | bmod = particle%bmod | |
| 263 |
2/2✓ Branch 0 taken 17316 times.
✓ Branch 1 taken 5772 times.
|
23088 | acov = particle%acov |
| 264 |
2/2✓ Branch 0 taken 17316 times.
✓ Branch 1 taken 5772 times.
|
23088 | e_zeta = particle%e(:, 3) |
| 265 |
2/2✓ Branch 0 taken 5772 times.
✓ Branch 1 taken 17316 times.
|
23088 | vpar = dot_product(state%v, b) |
| 266 |
2/2✓ Branch 0 taken 17316 times.
✓ Branch 1 taken 5772 times.
|
23088 | pzeta = dot_product(state%v, e_zeta) + acov(3)/ro0_bar |
| 267 | ! In an axisymmetric field toroidal canonical momentum is an exact | ||
| 268 | ! physical invariant, including across a | ||
| 269 | ! discontinuous radial sheet. Use the entry value when mapping back to | ||
| 270 | ! GC variables so a short recovery excursion cannot create a hidden | ||
| 271 | ! p_zeta kick. | ||
| 272 | 5772 | if (axisymmetric()) pzeta = state%entry_pzeta | |
| 273 |
2/2✓ Branch 0 taken 17316 times.
✓ Branch 1 taken 5772 times.
|
23088 | vperp = state%v - vpar*b |
| 274 |
2/2✓ Branch 1 taken 17316 times.
✓ Branch 2 taken 5772 times.
|
23088 | rho_l = (ro0_bar/bmod)*cross(b, vperp) |
| 275 |
2/2✓ Branch 0 taken 17316 times.
✓ Branch 1 taken 5772 times.
|
23088 | x_gc = particle%x - rho_l |
| 276 |
6/6✓ Branch 0 taken 17316 times.
✓ Branch 1 taken 5772 times.
✓ Branch 2 taken 51948 times.
✓ Branch 3 taken 17316 times.
✓ Branch 4 taken 17316 times.
✓ Branch 5 taken 5772 times.
|
92352 | u_gc = state%u - matmul(rho_l, particle%dual) |
| 277 | 5772 | owner_gc = state%owner | |
| 278 | 5772 | call locate_cart(x_gc, u_gc, owner_gc, status) | |
| 279 |
2/2✓ Branch 0 taken 5770 times.
✓ Branch 1 taken 2 times.
|
5772 | if (status /= SPECTRE_FO_OK) return |
| 280 | 5770 | call spectre_fo_field_local(u_gc, owner_gc, gc) | |
| 281 |
2/2✓ Branch 0 taken 17310 times.
✓ Branch 1 taken 5770 times.
|
23080 | b_gc = gc%b |
| 282 | 5770 | bmod = gc%bmod | |
| 283 |
2/2✓ Branch 0 taken 17310 times.
✓ Branch 1 taken 5770 times.
|
23080 | acov_gc = gc%acov |
| 284 |
2/2✓ Branch 0 taken 17310 times.
✓ Branch 1 taken 5770 times.
|
23080 | e_zeta_gc = gc%e(:, 3) |
| 285 |
5/6✓ Branch 0 taken 17310 times.
✓ Branch 1 taken 5770 times.
✓ Branch 2 taken 17310 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4996 times.
✓ Branch 5 taken 12314 times.
|
23080 | speed = norm2(state%v) |
| 286 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5770 times.
|
5770 | if (speed <= 0.0_dp) then |
| 287 | ✗ | status = SPECTRE_FO_INVALID_STATE | |
| 288 | ✗ | return | |
| 289 | end if | ||
| 290 |
2/2✓ Branch 0 taken 5770 times.
✓ Branch 1 taken 17310 times.
|
23080 | h_zeta = dot_product(b_gc, e_zeta_gc) |
| 291 |
6/8✓ Branch 0 taken 17310 times.
✓ Branch 1 taken 5770 times.
✓ Branch 2 taken 17310 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 8416 times.
✓ Branch 5 taken 8894 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 154 times.
|
23234 | if (axisymmetric() .and. abs(h_zeta) <= 100.0_dp*epsilon(1.0_dp)* & |
| 292 | max(norm2(e_zeta_gc), 1.0_dp)) then | ||
| 293 | ✗ | status = SPECTRE_FO_INVALID_STATE | |
| 294 | ✗ | return | |
| 295 | end if | ||
| 296 | 28234 | if (axisymmetric()) then | |
| 297 | ! Each volume has its own gauge (A vanishes on its inner boundary); | ||
| 298 | ! in axisymmetry A_zeta is constant on an interface, so the entry | ||
| 299 | ! momentum moves to the exit volume's gauge by the jumps of A_zeta | ||
| 300 | ! across the interfaces in between. | ||
| 301 | 154 | pzeta = pzeta + gauge_shift(state%entry_owner, owner_gc)/ro0_bar | |
| 302 | 154 | vpar = (pzeta - acov_gc(3)/ro0_bar)/h_zeta | |
| 303 | else | ||
| 304 |
2/2✓ Branch 0 taken 16848 times.
✓ Branch 1 taken 5616 times.
|
22464 | vpar = dot_product(state%v, b_gc) |
| 305 | end if | ||
| 306 | 5770 | perpendicular_energy = speed**2 - vpar**2 | |
| 307 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5770 times.
|
5770 | if (perpendicular_energy < -100.0_dp*epsilon(1.0_dp)*speed**2) then |
| 308 | ✗ | status = SPECTRE_FO_INVALID_STATE | |
| 309 | ✗ | return | |
| 310 | end if | ||
| 311 | 5770 | perpendicular_energy = max(perpendicular_energy, 0.0_dp) | |
| 312 |
2/2✓ Branch 0 taken 17310 times.
✓ Branch 1 taken 5770 times.
|
23080 | y(1:3) = u_gc |
| 313 | 5770 | y(4) = speed/sqrt(2.0_dp) | |
| 314 | 5770 | y(5) = vpar/speed | |
| 315 |
2/2✓ Branch 0 taken 3888 times.
✓ Branch 1 taken 1882 times.
|
5770 | if (present(mu_bar)) mu_bar = perpendicular_energy/(2.0_dp*bmod) |
| 316 | 5770 | owner = owner_gc | |
| 317 |
2/2✓ Branch 0 taken 1882 times.
✓ Branch 1 taken 3888 times.
|
5770 | if (.not. present(exited)) return |
| 318 | 1882 | radial_larmor = max(abs(state%u(1) - y(1)), 1.0e-8_dp) | |
| 319 |
1/2✓ Branch 0 taken 1882 times.
✗ Branch 1 not taken.
|
1882 | if (state%exit_rho > 0.0_dp) then |
| 320 | 1882 | exited = abs(y(1) - real(state%iface, dp)) > state%exit_rho | |
| 321 | ✗ | else if (state%iface == 0) then | |
| 322 | ✗ | exited = abs(y(1) - state%entry_rho) > 4.0_dp*radial_larmor | |
| 323 | else | ||
| 324 | ✗ | exited = abs(y(1) - real(state%iface, dp)) > 4.0_dp*radial_larmor | |
| 325 | end if | ||
| 326 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 1882 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1882 | exited = exited .and. abs(particle%bmod - gc%bmod)/gc%bmod <= 0.05_dp |
| 327 |
5/10✓ Branch 0 taken 5646 times.
✓ Branch 1 taken 1882 times.
✓ Branch 2 taken 5646 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 5646 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1882 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
9410 | exited = exited .and. norm2(particle%b - gc%b) <= 0.05_dp |
| 328 |
16/16✓ Branch 0 taken 17316 times.
✓ Branch 1 taken 5772 times.
✓ Branch 2 taken 17316 times.
✓ Branch 3 taken 5772 times.
✓ Branch 4 taken 51948 times.
✓ Branch 5 taken 17316 times.
✓ Branch 6 taken 17316 times.
✓ Branch 7 taken 5772 times.
✓ Branch 8 taken 51948 times.
✓ Branch 9 taken 17316 times.
✓ Branch 10 taken 17316 times.
✓ Branch 11 taken 5772 times.
✓ Branch 12 taken 17316 times.
✓ Branch 13 taken 5772 times.
✓ Branch 14 taken 1884 times.
✓ Branch 15 taken 3888 times.
|
196248 | end subroutine gc_map |
| 329 | 3491 | real(dp) function gauge_shift(from, to) result(shift) | |
| 330 | !> Sum of the A_zeta jumps (in the field's acov units) across the | ||
| 331 | !> interfaces from volume from to volume to, at theta = zeta = 0. | ||
| 332 | integer, intent(in) :: from, to | ||
| 333 | real(dp) :: x(3), b(3), bmod, a_low(3), a_high(3), e_zeta(3) | ||
| 334 | integer :: k, step, status | ||
| 335 | |||
| 336 | 3491 | shift = 0.0_dp | |
| 337 |
5/6✓ Branch 0 taken 1142 times.
✓ Branch 1 taken 2349 times.
✓ Branch 2 taken 1141 times.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 1141 times.
✗ Branch 5 not taken.
|
3491 | if (from == to .or. from < 1 .or. to < 1) return |
| 338 |
2/2✓ Branch 0 taken 271 times.
✓ Branch 1 taken 870 times.
|
1141 | step = merge(1, -1, to > from) |
| 339 |
2/2✓ Branch 0 taken 1141 times.
✓ Branch 1 taken 1141 times.
|
2282 | do k = min(from, to), max(from, to) - 1 |
| 340 | call logical_sample([real(k, dp), 0.0_dp, 0.0_dp], k, x, b, bmod, & | ||
| 341 |
2/2✓ Branch 0 taken 3423 times.
✓ Branch 1 taken 1141 times.
|
4564 | a_low, e_zeta, status) |
| 342 | call logical_sample([real(k, dp), 0.0_dp, 0.0_dp], k + 1, x, b, & | ||
| 343 |
2/2✓ Branch 0 taken 3423 times.
✓ Branch 1 taken 1141 times.
|
4564 | bmod, a_high, e_zeta, status) |
| 344 | 2282 | shift = shift + real(step, dp)*(a_high(3) - a_low(3)) | |
| 345 | end do | ||
| 346 | end function gauge_shift | ||
| 347 | 17805 | logical function axisymmetric() | |
| 348 | 17805 | axisymmetric = allocated(spectre_field) | |
| 349 |
15/20✓ Branch 0 taken 254 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 42 times.
✓ Branch 3 taken 212 times.
✓ Branch 4 taken 239 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 41 times.
✓ Branch 7 taken 198 times.
✓ Branch 8 taken 5772 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 156 times.
✓ Branch 11 taken 5616 times.
✓ Branch 12 taken 5770 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 154 times.
✓ Branch 15 taken 5616 times.
✓ Branch 16 taken 5770 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 5616 times.
✓ Branch 19 taken 154 times.
|
12189 | if (axisymmetric) axisymmetric = fo_field_axisymmetric |
| 350 | end function axisymmetric | ||
| 351 | 36 | subroutine spectre_fo_canonical_pzeta(state, ro0_bar, pzeta, status) | |
| 352 | type(spectre_fo_state_t), intent(inout) :: state | ||
| 353 | real(dp), intent(in) :: ro0_bar | ||
| 354 | real(dp), intent(out) :: pzeta | ||
| 355 | integer, intent(out) :: status | ||
| 356 | real(dp) :: b(3), bmod, acov(3), e_zeta(3) | ||
| 357 | 36 | pzeta = 0.0_dp | |
| 358 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
|
36 | if (ro0_bar <= 0.0_dp) then |
| 359 | ✗ | status = SPECTRE_FO_INVALID_STATE | |
| 360 | ✗ | return | |
| 361 | end if | ||
| 362 | call cart_sample(state%x, state%u, state%owner, b, bmod, acov, & | ||
| 363 | 36 | e_zeta, status) | |
| 364 |
1/2✓ Branch 0 taken 36 times.
✗ Branch 1 not taken.
|
36 | if (status /= SPECTRE_FO_OK) return |
| 365 |
2/2✓ Branch 0 taken 108 times.
✓ Branch 1 taken 36 times.
|
144 | pzeta = dot_product(state%v, e_zeta) + acov(3)/ro0_bar |
| 366 | end subroutine spectre_fo_canonical_pzeta | ||
| 367 | 14082 | subroutine rk45_substep(state, h, ro0_bar, status, angle_limit, particle) | |
| 368 | !> Adaptive Dormand-Prince 5(4) step of x' = v, v' = (|B|/ro0) v x b. On | ||
| 369 | !> the sharp field the stages continue the owner's spline polynomial | ||
| 370 | !> beyond the volume edge. Stage cross-and-return | ||
| 371 | !> trials are rejected; changed endpoints land on the first interface, | ||
| 372 | !> where the owner switches. On input h is the trial step, on output the step | ||
| 373 | !> taken. | ||
| 374 | type(spectre_fo_state_t), intent(inout) :: state | ||
| 375 | real(dp), intent(inout) :: h | ||
| 376 | real(dp), intent(in) :: ro0_bar, angle_limit | ||
| 377 | integer, intent(out) :: status | ||
| 378 | type(spectre_fo_local_t), intent(inout) :: particle | ||
| 379 | type(spectre_fo_local_t) :: final | ||
| 380 | real(dp) :: x1(3), v1(3), u1(3), err, a, b, fa, fb, c, fc, iface_rho | ||
| 381 | real(dp) :: trial_angle | ||
| 382 | integer :: attempt, owner1, it, side, grazes, stage_sides, opposite_stage | ||
| 383 | logical :: root_consistent | ||
| 384 | |||
| 385 | 14082 | grazes = 0 | |
| 386 |
1/2✓ Branch 0 taken 14157 times.
✗ Branch 1 not taken.
|
14157 | do attempt = 1, 60 |
| 387 | call dp_trial(state, h, ro0_bar, x1, v1, u1, owner1, err, status, & | ||
| 388 | 14157 | particle, final, stage_sides, trial_angle) | |
| 389 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 14150 times.
|
14157 | if (status /= SPECTRE_FO_OK) then |
| 390 | ! A failed stage has no embedded error estimate. | ||
| 391 | 7 | h = 0.5_dp*h | |
| 392 | 7 | cycle | |
| 393 | end if | ||
| 394 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 14149 times.
|
14150 | if (trial_angle > angle_limit*(1.0_dp + 4.0_dp*epsilon(1.0_dp))) then |
| 395 | ! A zero-time owner change can increase |B|. Recheck the cap | ||
| 396 | ! with this trial's initial field, including h_scale. | ||
| 397 | 1 | h = h*angle_limit/trial_angle | |
| 398 | 1 | cycle | |
| 399 | end if | ||
| 400 |
2/2✓ Branch 0 taken 65 times.
✓ Branch 1 taken 14084 times.
|
14149 | if (err > 1.0_dp) then |
| 401 | 65 | h = h*max(0.2_dp, min(0.5_dp, 0.9_dp*err**(-0.2_dp))) | |
| 402 | 65 | cycle | |
| 403 | end if | ||
| 404 |
2/2✓ Branch 0 taken 51 times.
✓ Branch 1 taken 14033 times.
|
14084 | if (stage_sides /= 0) then |
| 405 | if (owner1 == state%owner .or. & | ||
| 406 |
8/10✓ Branch 0 taken 50 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 26 times.
✓ Branch 3 taken 24 times.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 24 times.
✓ Branch 7 taken 26 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 24 times.
|
51 | (owner1 > state%owner .and. btest(stage_sides, STAGE_IN)) .or. & |
| 407 | (owner1 < state%owner .and. btest(stage_sides, STAGE_OUT))) then | ||
| 408 | ! A stage detects a cross-and-return, or an earlier sheet | ||
| 409 | ! in the opposite direction. Never commit the continuation. | ||
| 410 | 1 | h = 0.5_dp*h | |
| 411 | 1 | cycle | |
| 412 | end if | ||
| 413 | end if | ||
| 414 | 14083 | state%h_rk = h*min(5.0_dp, 0.9_dp*max(err, 1.0e-10_dp)**(-0.2_dp)) | |
| 415 |
2/2✓ Branch 0 taken 50 times.
✓ Branch 1 taken 14033 times.
|
14083 | if (owner1 /= state%owner) then |
| 416 | ! The first sheet of the initial owner, even if the endpoint | ||
| 417 | ! lies beyond a thin intermediate volume. | ||
| 418 | iface_rho = real(merge(state%owner, state%owner - 1, & | ||
| 419 |
2/2✓ Branch 0 taken 26 times.
✓ Branch 1 taken 24 times.
|
50 | owner1 > state%owner), dp) |
| 420 | 50 | opposite_stage = merge(STAGE_IN, STAGE_OUT, owner1 > state%owner) | |
| 421 | 50 | a = 0.0_dp | |
| 422 | 50 | fa = state%u(1) - iface_rho | |
| 423 | 50 | b = h | |
| 424 | 50 | fb = u1(1) - iface_rho | |
| 425 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 47 times.
|
50 | if (abs(fa) <= OWNER_EPS) then |
| 426 |
2/2✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 time.
|
3 | if (real(owner1 - state%owner, dp) & |
| 427 | *radial_velocity(state%u, state%v) > 0.0_dp) then | ||
| 428 | ! The outgoing normal velocity leaves the home side. | ||
| 429 | state%owner = merge(state%owner + 1, state%owner - 1, & | ||
| 430 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | owner1 > state%owner) |
| 431 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | if (state%owner > spectre_field%data%Mvol) then |
| 432 | 1 | state%owner = spectre_field%data%Mvol | |
| 433 | 1 | h = 0.0_dp | |
| 434 | 1 | status = SPECTRE_FO_LOSS | |
| 435 | return | ||
| 436 | end if | ||
| 437 | 1 | grazes = grazes + 1 | |
| 438 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (grazes >= 2) h = 0.5_dp*h |
| 439 | 1 | call spectre_fo_field_local(state%u, state%owner, particle) | |
| 440 | cycle | ||
| 441 | end if | ||
| 442 | ! Starts on the sheet moving into home: exclude the root | ||
| 443 | ! at t=0 and bracket the later return at positive time. | ||
| 444 | 1 | a = h/32.0_dp | |
| 445 | call dp_trial(state, a, ro0_bar, x1, v1, u1, owner1, err, & | ||
| 446 | 1 | status, particle, final, stage_sides) | |
| 447 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (status /= SPECTRE_FO_OK) then |
| 448 | ✗ | h = 0.5_dp*h | |
| 449 | ✗ | cycle | |
| 450 | end if | ||
| 451 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (btest(stage_sides, opposite_stage)) then |
| 452 | ✗ | h = 0.5_dp*h | |
| 453 | ✗ | cycle | |
| 454 | end if | ||
| 455 | 1 | fa = u1(1) - iface_rho | |
| 456 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (fa*fb >= 0.0_dp) then |
| 457 | ✗ | h = 0.5_dp*h | |
| 458 | ✗ | cycle | |
| 459 | end if | ||
| 460 | end if | ||
| 461 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
|
48 | if (fb == fa) then |
| 462 | ✗ | h = 0.5_dp*h | |
| 463 | ✗ | cycle | |
| 464 | end if | ||
| 465 | 211 | root_consistent = .true. | |
| 466 | side = 0 | ||
| 467 |
1/2✓ Branch 0 taken 211 times.
✗ Branch 1 not taken.
|
211 | do it = 1, 100 |
| 468 | 211 | c = (a*fb - b*fa)/(fb - fa) | |
| 469 | call dp_trial(state, c, ro0_bar, x1, v1, u1, owner1, err, & | ||
| 470 | 211 | status, particle, final, stage_sides) | |
| 471 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 211 times.
|
211 | if (status /= SPECTRE_FO_OK) exit |
| 472 |
1/2✓ Branch 0 taken 211 times.
✗ Branch 1 not taken.
|
211 | if (btest(stage_sides, opposite_stage)) then |
| 473 | ! Shortened trials can expose an earlier opposite sheet | ||
| 474 | ! missed by the nominal stage grid. Retry from the start. | ||
| 475 | root_consistent = .false. | ||
| 476 | exit | ||
| 477 | end if | ||
| 478 | 211 | fc = u1(1) - iface_rho | |
| 479 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 163 times.
|
211 | if (abs(fc) <= 1.0e-12_dp) exit |
| 480 |
2/2✓ Branch 1 taken 81 times.
✓ Branch 2 taken 82 times.
|
422 | if (fc*fb > 0.0_dp) then |
| 481 | 81 | b = c | |
| 482 | 81 | fb = fc | |
| 483 |
2/2✓ Branch 0 taken 25 times.
✓ Branch 1 taken 56 times.
|
81 | if (side == -1) fa = 0.5_dp*fa |
| 484 | side = -1 | ||
| 485 | else | ||
| 486 | 82 | a = c | |
| 487 | 82 | fa = fc | |
| 488 |
2/2✓ Branch 0 taken 26 times.
✓ Branch 1 taken 56 times.
|
82 | if (side == 1) fb = 0.5_dp*fb |
| 489 | side = 1 | ||
| 490 | end if | ||
| 491 | end do | ||
| 492 |
2/4✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 48 times.
|
48 | if (status /= SPECTRE_FO_OK .or. .not. root_consistent) then |
| 493 | ✗ | h = 0.5_dp*h | |
| 494 | ✗ | cycle | |
| 495 | end if | ||
| 496 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
|
48 | if (c <= 1.0e-10_dp*h) then |
| 497 | ! The orbit grazes the interface at the step start: it | ||
| 498 | ! belongs to the other volume. Switch and redo the step. | ||
| 499 | ! Repeated grazing must reduce/retry, never bypass landing. | ||
| 500 | state%owner = merge(state%owner + 1, state%owner - 1, & | ||
| 501 | ✗ | iface_rho >= real(state%owner, dp)) | |
| 502 | ✗ | state%u(1) = iface_rho | |
| 503 | ✗ | if (state%owner > spectre_field%data%Mvol) then | |
| 504 | ✗ | state%owner = spectre_field%data%Mvol | |
| 505 | ✗ | h = 0.0_dp | |
| 506 | ✗ | status = SPECTRE_FO_LOSS | |
| 507 | ✗ | return | |
| 508 | end if | ||
| 509 | ✗ | grazes = grazes + 1 | |
| 510 | ✗ | if (grazes >= 2) h = 0.5_dp*h | |
| 511 | ✗ | call spectre_fo_field_local(state%u, state%owner, particle) | |
| 512 | ✗ | cycle | |
| 513 | end if | ||
| 514 |
2/4✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 48 times.
|
48 | if (it > 100 .or. err > 1.0_dp) then |
| 515 | ✗ | h = 0.5_dp*h | |
| 516 | ✗ | cycle | |
| 517 | end if | ||
| 518 | 48 | h = c | |
| 519 | 48 | u1(1) = iface_rho | |
| 520 | owner1 = merge(state%owner + 1, state%owner - 1, & | ||
| 521 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 24 times.
|
48 | iface_rho >= real(state%owner, dp)) |
| 522 | end if | ||
| 523 |
2/2✓ Branch 0 taken 14033 times.
✓ Branch 1 taken 48 times.
|
14081 | if (owner1 == state%owner) then |
| 524 | ! The seventh RK stage is the accepted endpoint. Reuse its | ||
| 525 | ! sample for GC reconstruction and the next first stage. | ||
| 526 | 14033 | particle = final | |
| 527 |
2/2✓ Branch 0 taken 45 times.
✓ Branch 1 taken 3 times.
|
48 | else if (owner1 <= spectre_field%data%Mvol) then |
| 528 | ! A landed sheet changes both the exact coordinate and owner. | ||
| 529 | 45 | call spectre_fo_field_local(u1, owner1, particle) | |
| 530 | end if | ||
| 531 | 14081 | owner1 = min(owner1, spectre_field%data%Mvol) | |
| 532 |
2/2✓ Branch 0 taken 42243 times.
✓ Branch 1 taken 14081 times.
|
56324 | state%x = x1 |
| 533 |
2/2✓ Branch 0 taken 42243 times.
✓ Branch 1 taken 14081 times.
|
56324 | state%v = v1 |
| 534 |
2/2✓ Branch 0 taken 42243 times.
✓ Branch 1 taken 14081 times.
|
56324 | state%u = u1 |
| 535 | 14081 | state%owner = owner1 | |
| 536 | 14081 | status = SPECTRE_FO_OK | |
| 537 |
4/4✓ Branch 0 taken 605 times.
✓ Branch 1 taken 13476 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 602 times.
|
14081 | if (owner1 == spectre_field%data%Mvol .and. & |
| 538 | 3 | u1(1) >= real(owner1, dp) - OWNER_EPS) status = SPECTRE_FO_LOSS | |
| 539 | 14157 | return | |
| 540 | end do | ||
| 541 | ✗ | status = SPECTRE_FO_STEP_LIMIT | |
| 542 |
14/14✓ Branch 0 taken 42246 times.
✓ Branch 1 taken 14082 times.
✓ Branch 2 taken 42246 times.
✓ Branch 3 taken 14082 times.
✓ Branch 4 taken 126738 times.
✓ Branch 5 taken 42246 times.
✓ Branch 6 taken 42246 times.
✓ Branch 7 taken 14082 times.
✓ Branch 8 taken 126738 times.
✓ Branch 9 taken 42246 times.
✓ Branch 10 taken 42246 times.
✓ Branch 11 taken 14082 times.
✓ Branch 12 taken 42246 times.
✓ Branch 13 taken 14082 times.
|
478788 | end subroutine rk45_substep |
| 543 |
14/14✓ Branch 0 taken 43107 times.
✓ Branch 1 taken 14369 times.
✓ Branch 2 taken 43107 times.
✓ Branch 3 taken 14369 times.
✓ Branch 4 taken 129321 times.
✓ Branch 5 taken 43107 times.
✓ Branch 6 taken 43107 times.
✓ Branch 7 taken 14369 times.
✓ Branch 8 taken 129321 times.
✓ Branch 9 taken 43107 times.
✓ Branch 10 taken 43107 times.
✓ Branch 11 taken 14369 times.
✓ Branch 12 taken 43107 times.
✓ Branch 13 taken 14369 times.
|
502915 | subroutine dp_trial(state, h, ro0_bar, x1, v1, u1, owner1, err, status, & |
| 544 | initial, final, stage_sides, gyro_angle) | ||
| 545 | !> Dormand-Prince on logical position and Cartesian velocity. The | ||
| 546 | !> physical Lorentz equation is unchanged; du/dt = E^{-1} v replaces | ||
| 547 | !> the iterative Cartesian-to-logical inversion at every RK stage. | ||
| 548 | type(spectre_fo_state_t), intent(in) :: state | ||
| 549 | real(dp), intent(in) :: h, ro0_bar | ||
| 550 | real(dp), intent(out) :: x1(3), v1(3), u1(3), err | ||
| 551 | integer, intent(out) :: owner1, status | ||
| 552 | type(spectre_fo_local_t), intent(in) :: initial | ||
| 553 | type(spectre_fo_local_t), intent(out) :: final | ||
| 554 | integer, intent(out), optional :: stage_sides | ||
| 555 | real(dp), intent(out), optional :: gyro_angle | ||
| 556 | type(spectre_fo_local_t) :: loc | ||
| 557 | real(dp) :: ku(3, 7), kv(3, 7), us(3), vs(3), qs(3), q0(3), w(7, 7) | ||
| 558 | real(dp) :: rho_l, eu(3), ev(3), ex(3) | ||
| 559 | integer :: i, j, owner_geo | ||
| 560 | logical :: cartesian | ||
| 561 | |||
| 562 | ! Cartesian position remains a regular chart near the polar axis. | ||
| 563 | ! That small recovery region uses the same spline provider, with an | ||
| 564 | ! inverse per stage; the regular logical path never performs one. | ||
| 565 |
4/4✓ Branch 0 taken 13626 times.
✓ Branch 1 taken 743 times.
✓ Branch 2 taken 13624 times.
✓ Branch 3 taken 2 times.
|
28736 | cartesian = state%owner == 1 .and. state%u(1) < FO_AXIS_RADIUS |
| 566 |
2/2✓ Branch 0 taken 43107 times.
✓ Branch 1 taken 14369 times.
|
57476 | q0 = state%u |
| 567 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 14367 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
|
14375 | if (cartesian) q0 = state%x |
| 568 |
2/2✓ Branch 0 taken 43107 times.
✓ Branch 1 taken 14369 times.
|
57476 | u1 = state%u |
| 569 |
1/2✓ Branch 0 taken 14369 times.
✗ Branch 1 not taken.
|
14369 | if (present(stage_sides)) stage_sides = 0 |
| 570 | 14369 | status = SPECTRE_FO_OK | |
| 571 | 14369 | w = 0.0_dp | |
| 572 | 14369 | w(2, 1) = 1.0_dp/5 | |
| 573 |
2/2✓ Branch 0 taken 28738 times.
✓ Branch 1 taken 14369 times.
|
43107 | w(3, 1:2) = [3.0_dp/40, 9.0_dp/40] |
| 574 |
2/2✓ Branch 0 taken 43107 times.
✓ Branch 1 taken 14369 times.
|
57476 | w(4, 1:3) = [44.0_dp/45, -56.0_dp/15, 32.0_dp/9] |
| 575 | w(5, 1:4) = [19372.0_dp/6561, -25360.0_dp/2187, 64448.0_dp/6561, & | ||
| 576 |
2/2✓ Branch 0 taken 57476 times.
✓ Branch 1 taken 14369 times.
|
71845 | -212.0_dp/729] |
| 577 | w(6, 1:5) = [9017.0_dp/3168, -355.0_dp/33, 46732.0_dp/5247, & | ||
| 578 |
2/2✓ Branch 0 taken 71845 times.
✓ Branch 1 taken 14369 times.
|
86214 | 49.0_dp/176, -5103.0_dp/18656] |
| 579 | w(7, 1:6) = [35.0_dp/384, 0.0_dp, 500.0_dp/1113, 125.0_dp/192, & | ||
| 580 |
2/2✓ Branch 0 taken 86214 times.
✓ Branch 1 taken 14369 times.
|
100583 | -2187.0_dp/6784, 11.0_dp/84] |
| 581 |
2/2✓ Branch 0 taken 100552 times.
✓ Branch 1 taken 14362 times.
|
114914 | do i = 1, 7 |
| 582 | 100552 | qs = q0 | |
| 583 |
2/2✓ Branch 0 taken 301656 times.
✓ Branch 1 taken 100552 times.
|
402208 | vs = state%v |
| 584 |
2/2✓ Branch 0 taken 301619 times.
✓ Branch 1 taken 100552 times.
|
402171 | do j = 1, i - 1 |
| 585 |
2/2✓ Branch 0 taken 904857 times.
✓ Branch 1 taken 301619 times.
|
1206476 | qs = qs + h*w(i, j)*ku(:, j) |
| 586 |
2/2✓ Branch 0 taken 904857 times.
✓ Branch 1 taken 301619 times.
|
1307028 | vs = vs + h*w(i, j)*kv(:, j) |
| 587 | end do | ||
| 588 | 100552 | us = qs | |
| 589 | 100552 | owner_geo = state%owner | |
| 590 |
2/2✓ Branch 0 taken 14 times.
✓ Branch 1 taken 100538 times.
|
100552 | if (cartesian) then |
| 591 | 14 | call locate_cart(qs, u1, owner_geo, status) | |
| 592 |
1/2✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
|
21 | if (status /= SPECTRE_FO_OK) return |
| 593 | 14 | us = u1 | |
| 594 | else | ||
| 595 | 100538 | owner_geo = logical_owner(us(1), state%owner) | |
| 596 | end if | ||
| 597 |
1/2✓ Branch 0 taken 100552 times.
✗ Branch 1 not taken.
|
100552 | if (present(stage_sides)) then |
| 598 |
2/2✓ Branch 0 taken 170 times.
✓ Branch 1 taken 100382 times.
|
100552 | if (owner_geo < state%owner) & |
| 599 | 170 | stage_sides = ibset(stage_sides, STAGE_IN) | |
| 600 |
2/2✓ Branch 0 taken 197 times.
✓ Branch 1 taken 100355 times.
|
100552 | if (owner_geo > state%owner) & |
| 601 | 197 | stage_sides = ibset(stage_sides, STAGE_OUT) | |
| 602 | end if | ||
| 603 |
3/4✓ Branch 0 taken 7 times.
✓ Branch 1 taken 100545 times.
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
|
100552 | if (us(1) <= 0.0_dp .and. .not. cartesian) then |
| 604 | 7 | status = SPECTRE_FO_LOCATE_FAIL | |
| 605 | 7 | return | |
| 606 | end if | ||
| 607 | ! Continue only the initial owner's polynomial for trial stages; | ||
| 608 | ! accepted sheet crossings are landed and then change owner. | ||
| 609 |
2/2✓ Branch 0 taken 14367 times.
✓ Branch 1 taken 86178 times.
|
100545 | if (i == 1 .and. .not. cartesian) then |
| 610 | 14367 | loc = initial | |
| 611 | else | ||
| 612 | 86178 | call spectre_fo_field_local(us, state%owner, loc) | |
| 613 | end if | ||
| 614 |
2/2✓ Branch 0 taken 14369 times.
✓ Branch 1 taken 86176 times.
|
100545 | if (i == 1) then |
| 615 |
6/6✓ Branch 0 taken 43107 times.
✓ Branch 1 taken 14369 times.
✓ Branch 2 taken 43083 times.
✓ Branch 3 taken 24 times.
✓ Branch 4 taken 12574 times.
✓ Branch 5 taken 30509 times.
|
57476 | rho_l = ro0_bar*norm2(vs)/loc%bmod |
| 616 |
2/2✓ Branch 0 taken 14157 times.
✓ Branch 1 taken 212 times.
|
14369 | if (present(gyro_angle)) gyro_angle = h*loc%bmod/ro0_bar |
| 617 | end if | ||
| 618 |
2/2✓ Branch 0 taken 301635 times.
✓ Branch 1 taken 100545 times.
|
402180 | ku(:, i) = vs |
| 619 |
8/8✓ Branch 0 taken 100531 times.
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 301593 times.
✓ Branch 3 taken 100531 times.
✓ Branch 4 taken 301593 times.
✓ Branch 5 taken 100531 times.
✓ Branch 6 taken 904779 times.
✓ Branch 7 taken 301593 times.
|
1608510 | if (.not. cartesian) ku(:, i) = matmul(vs, loc%dual) |
| 620 |
2/2✓ Branch 1 taken 301635 times.
✓ Branch 2 taken 100545 times.
|
416542 | kv(:, i) = (loc%bmod/ro0_bar)*cross(vs, loc%b) |
| 621 | end do | ||
| 622 | 14362 | u1 = us | |
| 623 | 14362 | v1 = vs | |
| 624 |
2/2✓ Branch 0 taken 43086 times.
✓ Branch 1 taken 14362 times.
|
57448 | x1 = loc%x |
| 625 | 14362 | final = loc | |
| 626 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 14360 times.
|
14362 | if (cartesian) x1 = qs |
| 627 | 14362 | owner1 = logical_owner(u1(1), state%owner) | |
| 628 | eu = h*(71.0_dp/57600*ku(:, 1) - 71.0_dp/16695*ku(:, 3) + & | ||
| 629 | 71.0_dp/1920*ku(:, 4) - 17253.0_dp/339200*ku(:, 5) + & | ||
| 630 |
2/2✓ Branch 0 taken 43086 times.
✓ Branch 1 taken 14362 times.
|
57448 | 22.0_dp/525*ku(:, 6) - 1.0_dp/40*ku(:, 7)) |
| 631 | ev = h*(71.0_dp/57600*kv(:, 1) - 71.0_dp/16695*kv(:, 3) + & | ||
| 632 | 71.0_dp/1920*kv(:, 4) - 17253.0_dp/339200*kv(:, 5) + & | ||
| 633 |
2/2✓ Branch 0 taken 43086 times.
✓ Branch 1 taken 14362 times.
|
57448 | 22.0_dp/525*kv(:, 6) - 1.0_dp/40*kv(:, 7)) |
| 634 | 14362 | ex = eu | |
| 635 |
6/6✓ Branch 0 taken 14360 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 43080 times.
✓ Branch 3 taken 14360 times.
✓ Branch 4 taken 129240 times.
✓ Branch 5 taken 43080 times.
|
186682 | if (.not. cartesian) ex = matmul(loc%e, eu) |
| 636 |
16/18✓ Branch 0 taken 14362 times.
✓ Branch 1 taken 43086 times.
✓ Branch 2 taken 43083 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 43083 times.
✓ Branch 6 taken 14362 times.
✓ Branch 7 taken 43086 times.
✓ Branch 8 taken 43077 times.
✓ Branch 9 taken 9 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 43077 times.
✓ Branch 12 taken 43086 times.
✓ Branch 13 taken 14362 times.
✓ Branch 14 taken 43062 times.
✓ Branch 15 taken 24 times.
✓ Branch 16 taken 12574 times.
✓ Branch 17 taken 30488 times.
|
172344 | err = max(norm2(ex)/rho_l, norm2(ev)/norm2(state%v))/fo_tolerance |
| 637 |
14/14✓ Branch 0 taken 43107 times.
✓ Branch 1 taken 14369 times.
✓ Branch 2 taken 43107 times.
✓ Branch 3 taken 14369 times.
✓ Branch 4 taken 129321 times.
✓ Branch 5 taken 43107 times.
✓ Branch 6 taken 43107 times.
✓ Branch 7 taken 14369 times.
✓ Branch 8 taken 129321 times.
✓ Branch 9 taken 43107 times.
✓ Branch 10 taken 43107 times.
✓ Branch 11 taken 14369 times.
✓ Branch 12 taken 43107 times.
✓ Branch 13 taken 14369 times.
|
488546 | end subroutine dp_trial |
| 638 | |||
| 639 | 114900 | integer function logical_owner(rho, home) result(owner) | |
| 640 | real(dp), intent(in) :: rho | ||
| 641 | integer, intent(in) :: home | ||
| 642 | |||
| 643 | 114900 | owner = home | |
| 644 |
2/2✓ Branch 0 taken 223 times.
✓ Branch 1 taken 114677 times.
|
114900 | if (rho < real(home - 1, dp) - OWNER_EPS) & |
| 645 | 223 | owner = max(1, int(rho) + 1) | |
| 646 |
2/2✓ Branch 0 taken 277 times.
✓ Branch 1 taken 114623 times.
|
114900 | if (rho > real(home, dp) + OWNER_EPS) & |
| 647 | 277 | owner = min(spectre_field%data%Mvol + 1, int(rho) + 1) | |
| 648 | 114900 | end function logical_owner | |
| 649 | |||
| 650 | 5709 | subroutine synchronize_position(state, status) | |
| 651 | type(spectre_fo_state_t), intent(inout) :: state | ||
| 652 | integer, intent(out) :: status | ||
| 653 | type(spectre_fo_local_t) :: loc | ||
| 654 | |||
| 655 | 5709 | call spectre_fo_field_local(state%u, state%owner, loc) | |
| 656 | 5709 | status = SPECTRE_FO_OK | |
| 657 |
14/14✓ Branch 0 taken 5709 times.
✓ Branch 1 taken 17127 times.
✓ Branch 2 taken 6189 times.
✓ Branch 3 taken 10938 times.
✓ Branch 4 taken 109 times.
✓ Branch 5 taken 6080 times.
✓ Branch 6 taken 17127 times.
✓ Branch 7 taken 5709 times.
✓ Branch 8 taken 16161 times.
✓ Branch 9 taken 966 times.
✓ Branch 10 taken 5749 times.
✓ Branch 11 taken 10412 times.
✓ Branch 12 taken 2026 times.
✓ Branch 13 taken 3683 times.
|
45672 | if (norm2(state%x - loc%x) > 1.0e-12_dp*max(norm2(loc%x), 1.0_dp)) & |
| 658 | 2026 | call locate_cart(state%x, state%u, state%owner, status) | |
| 659 |
14/14✓ Branch 0 taken 17127 times.
✓ Branch 1 taken 5709 times.
✓ Branch 2 taken 17127 times.
✓ Branch 3 taken 5709 times.
✓ Branch 4 taken 51381 times.
✓ Branch 5 taken 17127 times.
✓ Branch 6 taken 17127 times.
✓ Branch 7 taken 5709 times.
✓ Branch 8 taken 51381 times.
✓ Branch 9 taken 17127 times.
✓ Branch 10 taken 17127 times.
✓ Branch 11 taken 5709 times.
✓ Branch 12 taken 17127 times.
✓ Branch 13 taken 5709 times.
|
194106 | end subroutine synchronize_position |
| 660 | |||
| 661 | 42 | subroutine solve_gyrophase(x_gc, u_gc, owner, b, e1, e2, vpar, vperp, & | |
| 662 | target, ro0_bar, alpha, status) | ||
| 663 | real(dp), intent(in) :: x_gc(3), u_gc(3), b(3), e1(3), e2(3) | ||
| 664 | real(dp), intent(in) :: vpar, vperp, target, ro0_bar | ||
| 665 | integer, intent(in) :: owner | ||
| 666 | real(dp), intent(out) :: alpha | ||
| 667 | integer, intent(out) :: status | ||
| 668 | real(dp) :: a0, a1, f0, f1, amid, fmid, scale, best_residual, best_alpha | ||
| 669 | integer :: i, ierr | ||
| 670 | logical :: have_previous | ||
| 671 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
|
42 | if (vperp <= sqrt(epsilon(1.0_dp))) then |
| 672 | ✗ | alpha = 0.0_dp | |
| 673 | call phase_residual(alpha, x_gc, u_gc, owner, b, e1, e2, vpar, & | ||
| 674 | ✗ | vperp, target, ro0_bar, f0, ierr) | |
| 675 | ✗ | scale = max(abs(target), 1.0_dp) | |
| 676 | ✗ | if (ierr /= SPECTRE_FO_OK) then | |
| 677 | ✗ | status = ierr | |
| 678 | ✗ | else if (abs(f0) <= 1.0e-11_dp*scale) then | |
| 679 | ✗ | status = SPECTRE_FO_OK | |
| 680 | else | ||
| 681 | ✗ | status = SPECTRE_FO_GYROPHASE_FAIL | |
| 682 | end if | ||
| 683 | 4 | return | |
| 684 | end if | ||
| 685 | have_previous = .false. | ||
| 686 | best_residual = huge(1.0_dp) | ||
| 687 | best_alpha = 0.0_dp | ||
| 688 |
2/2✓ Branch 0 taken 1018 times.
✓ Branch 1 taken 4 times.
|
1022 | do i = 0, NPHASE |
| 689 | 1018 | a1 = 8.0_dp*atan(1.0_dp)*real(i, dp)/real(NPHASE, dp) | |
| 690 | call phase_residual(a1, x_gc, u_gc, owner, b, e1, e2, vpar, & | ||
| 691 | 1018 | vperp, target, ro0_bar, f1, ierr) | |
| 692 |
2/2✓ Branch 0 taken 1016 times.
✓ Branch 1 taken 2 times.
|
1018 | if (ierr /= SPECTRE_FO_OK) then |
| 693 | have_previous = .false. | ||
| 694 | cycle | ||
| 695 | end if | ||
| 696 |
2/2✓ Branch 0 taken 541 times.
✓ Branch 1 taken 475 times.
|
1016 | if (abs(f1) < best_residual) then |
| 697 | 541 | best_residual = abs(f1) | |
| 698 | 541 | best_alpha = a1 | |
| 699 | end if | ||
| 700 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1016 times.
|
1016 | if (f1 == 0.0_dp) then |
| 701 | ✗ | alpha = a1 | |
| 702 | ✗ | status = SPECTRE_FO_OK | |
| 703 | ✗ | return | |
| 704 | end if | ||
| 705 |
2/2✓ Branch 0 taken 974 times.
✓ Branch 1 taken 42 times.
|
1016 | if (have_previous) then |
| 706 |
2/2✓ Branch 0 taken 38 times.
✓ Branch 1 taken 936 times.
|
974 | if (f0*f1 <= 0.0_dp) exit |
| 707 | end if | ||
| 708 | 978 | a0 = a1 | |
| 709 | 978 | f0 = f1 | |
| 710 | 2038 | have_previous = .true. | |
| 711 | end do | ||
| 712 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 38 times.
|
42 | if (i > NPHASE) then |
| 713 | ! The gyrophase variation of p_zeta, proportional to v_perp, does | ||
| 714 | ! not reach the guiding-centre value (strongly passing markers): | ||
| 715 | ! take the closest gyrophase only with a bounded momentum defect. | ||
| 716 | ! Entry records the actual signed defect; no order claim is made. | ||
| 717 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
|
4 | if (best_residual > ENTRY_MOMENTUM_TOL*max(abs(target), 1.0_dp)) then |
| 718 | 1 | status = SPECTRE_FO_GYROPHASE_FAIL | |
| 719 | 1 | return | |
| 720 | end if | ||
| 721 | 3 | alpha = best_alpha | |
| 722 | 3 | status = SPECTRE_FO_OK | |
| 723 | 3 | return | |
| 724 | end if | ||
| 725 |
2/2✓ Branch 0 taken 2280 times.
✓ Branch 1 taken 38 times.
|
2318 | do i = 1, 60 |
| 726 | 2280 | amid = 0.5_dp*(a0 + a1) | |
| 727 | call phase_residual(amid, x_gc, u_gc, owner, b, e1, e2, vpar, & | ||
| 728 | 2280 | vperp, target, ro0_bar, fmid, ierr) | |
| 729 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2280 times.
|
2280 | if (ierr /= SPECTRE_FO_OK) then |
| 730 | ✗ | status = ierr | |
| 731 | ✗ | return | |
| 732 | end if | ||
| 733 |
2/2✓ Branch 1 taken 1170 times.
✓ Branch 2 taken 1110 times.
|
4598 | if (f0*fmid <= 0.0_dp) then |
| 734 | 1170 | a1 = amid | |
| 735 | else | ||
| 736 | 1110 | a0 = amid | |
| 737 | 1110 | f0 = fmid | |
| 738 | end if | ||
| 739 | end do | ||
| 740 | 38 | alpha = 0.5_dp*(a0 + a1) | |
| 741 | 38 | status = SPECTRE_FO_OK | |
| 742 | end subroutine solve_gyrophase | ||
| 743 | 10017 | subroutine phase_residual(alpha, x_gc, u_gc, owner, b_gc, e1, e2, vpar, & | |
| 744 | vperp, target, ro0_bar, residual, status) | ||
| 745 | real(dp), intent(in) :: alpha, x_gc(3), u_gc(3), b_gc(3), e1(3), e2(3) | ||
| 746 | real(dp), intent(in) :: vpar, vperp, target, ro0_bar | ||
| 747 | integer, intent(in) :: owner | ||
| 748 | real(dp), intent(out) :: residual | ||
| 749 | integer, intent(out) :: status | ||
| 750 | real(dp) :: eperp(3), x(3), u(3), b(3), bmod, acov(3), e_zeta(3), v(3) | ||
| 751 | integer :: trial_owner | ||
| 752 |
2/2✓ Branch 0 taken 10017 times.
✓ Branch 1 taken 3339 times.
|
13356 | eperp = cos(alpha)*e1 + sin(alpha)*e2 |
| 753 | 3339 | call logical_sample(u_gc, owner, x, b, bmod, acov, e_zeta, status) | |
| 754 |
1/2✓ Branch 0 taken 3339 times.
✗ Branch 1 not taken.
|
3341 | if (status /= SPECTRE_FO_OK) return |
| 755 |
4/4✓ Branch 0 taken 10017 times.
✓ Branch 1 taken 3339 times.
✓ Branch 3 taken 10017 times.
✓ Branch 4 taken 3339 times.
|
23373 | x = x_gc + (ro0_bar/bmod)*cross(b_gc, vperp*eperp) |
| 756 | 3339 | u = u_gc | |
| 757 | 3339 | trial_owner = owner | |
| 758 | 3339 | call cart_sample(x, u, trial_owner, b, bmod, acov, e_zeta, status) | |
| 759 |
2/2✓ Branch 0 taken 3337 times.
✓ Branch 1 taken 2 times.
|
3339 | if (status /= SPECTRE_FO_OK) return |
| 760 |
2/2✓ Branch 0 taken 10011 times.
✓ Branch 1 taken 3337 times.
|
13348 | v = vpar*b_gc + vperp*eperp |
| 761 | residual = dot_product(v, e_zeta) & | ||
| 762 |
2/2✓ Branch 0 taken 10011 times.
✓ Branch 1 taken 3337 times.
|
13348 | + (acov(3) - gauge_shift(owner, trial_owner))/ro0_bar - target |
| 763 | end subroutine phase_residual | ||
| 764 | 3375 | subroutine cart_sample(x, u, owner, b, bmod, acov, e_zeta, status) | |
| 765 | real(dp), intent(in) :: x(3) | ||
| 766 | real(dp), intent(inout) :: u(3) | ||
| 767 | integer, intent(inout) :: owner | ||
| 768 | real(dp), intent(out) :: b(3), bmod, acov(3), e_zeta(3) | ||
| 769 | integer, intent(out) :: status | ||
| 770 | real(dp) :: x_at_u(3) | ||
| 771 | 3375 | call locate_cart(x, u, owner, status) | |
| 772 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 3373 times.
|
3375 | if (status /= SPECTRE_FO_OK) return |
| 773 | 3373 | call logical_sample(u, owner, x_at_u, b, bmod, acov, e_zeta, status) | |
| 774 | end subroutine cart_sample | ||
| 775 | 19604 | subroutine logical_sample(u, owner, x, b, bmod, acov, e_zeta, status) | |
| 776 | real(dp), intent(in) :: u(3) | ||
| 777 | integer, intent(in) :: owner | ||
| 778 | real(dp), intent(out) :: x(3), b(3), bmod, acov(3), e_zeta(3) | ||
| 779 | integer, intent(out) :: status | ||
| 780 | type(spectre_fo_local_t) :: loc | ||
| 781 | |||
| 782 | 19604 | call spectre_fo_field_local(u, owner, loc) | |
| 783 |
2/2✓ Branch 0 taken 19604 times.
✓ Branch 1 taken 58812 times.
|
78416 | x = loc%x |
| 784 |
2/2✓ Branch 0 taken 58812 times.
✓ Branch 1 taken 19604 times.
|
78416 | b = loc%b |
| 785 | 19604 | bmod = loc%bmod | |
| 786 |
2/2✓ Branch 0 taken 19604 times.
✓ Branch 1 taken 58812 times.
|
78416 | acov = loc%acov |
| 787 |
2/2✓ Branch 0 taken 58812 times.
✓ Branch 1 taken 19604 times.
|
78416 | e_zeta = loc%e(:, 3) |
| 788 | 19604 | status = SPECTRE_FO_OK | |
| 789 |
14/14✓ Branch 0 taken 19604 times.
✓ Branch 1 taken 58812 times.
✓ Branch 2 taken 19604 times.
✓ Branch 3 taken 58812 times.
✓ Branch 4 taken 176436 times.
✓ Branch 5 taken 58812 times.
✓ Branch 6 taken 19604 times.
✓ Branch 7 taken 58812 times.
✓ Branch 8 taken 176436 times.
✓ Branch 9 taken 58812 times.
✓ Branch 10 taken 19604 times.
✓ Branch 11 taken 58812 times.
✓ Branch 12 taken 58812 times.
✓ Branch 13 taken 19604 times.
|
744952 | end subroutine logical_sample |
| 790 | |||
| 791 | 11187 | subroutine locate_cart(x, u, owner, status) | |
| 792 | real(dp), intent(in) :: x(3) | ||
| 793 | real(dp), intent(inout) :: u(3) | ||
| 794 | integer, intent(inout) :: owner | ||
| 795 | integer, intent(out) :: status | ||
| 796 | integer :: ierr | ||
| 797 | |||
| 798 | 11187 | call spectre_fo_field_locate(x, u, owner, ierr) | |
| 799 | 11187 | status = SPECTRE_FO_OK | |
| 800 |
6/8✓ Branch 0 taken 2 times.
✓ Branch 1 taken 5770 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2026 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 3373 times.
|
11187 | if (ierr /= 0) status = SPECTRE_FO_LOCATE_FAIL |
| 801 | end subroutine locate_cart | ||
| 802 | |||
| 803 | 3 | real(dp) function radial_velocity(u, v) result(vrho) | |
| 804 | real(dp), intent(in) :: u(3), v(3) | ||
| 805 | type(spectre_fo_local_t) :: loc | ||
| 806 | |||
| 807 | call spectre_fo_field_local(u, min(int(u(1)) + 1, & | ||
| 808 | 3 | spectre_field%data%Mvol), loc) | |
| 809 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
|
12 | vrho = dot_product(v, loc%dual(:, 1)) |
| 810 |
14/14✓ Branch 0 taken 3 times.
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 9 times.
✓ Branch 4 taken 27 times.
✓ Branch 5 taken 9 times.
✓ Branch 6 taken 3 times.
✓ Branch 7 taken 9 times.
✓ Branch 8 taken 27 times.
✓ Branch 9 taken 9 times.
✓ Branch 10 taken 3 times.
✓ Branch 11 taken 9 times.
✓ Branch 12 taken 9 times.
✓ Branch 13 taken 3 times.
|
114 | end function radial_velocity |
| 811 | |||
| 812 | 254 | pure subroutine perpendicular_frame(b, e1, e2) | |
| 813 | real(dp), intent(in) :: b(3) | ||
| 814 | real(dp), intent(out) :: e1(3), e2(3) | ||
| 815 | real(dp) :: axis(3) | ||
| 816 | |||
| 817 |
1/2✓ Branch 0 taken 254 times.
✗ Branch 1 not taken.
|
254 | if (abs(b(3)) < 0.9_dp) then |
| 818 | 254 | axis = [0.0_dp, 0.0_dp, 1.0_dp] | |
| 819 | else | ||
| 820 | ✗ | axis = [1.0_dp, 0.0_dp, 0.0_dp] | |
| 821 | end if | ||
| 822 |
4/4✓ Branch 0 taken 254 times.
✓ Branch 1 taken 762 times.
✓ Branch 2 taken 762 times.
✓ Branch 3 taken 254 times.
|
2032 | e1 = axis - dot_product(axis, b)*b |
| 823 |
6/8✓ Branch 0 taken 762 times.
✓ Branch 1 taken 254 times.
✓ Branch 2 taken 762 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 762 times.
✓ Branch 6 taken 762 times.
✓ Branch 7 taken 254 times.
|
1778 | e1 = e1/norm2(e1) |
| 824 | 254 | e2 = cross(b, e1) | |
| 825 | 254 | end subroutine perpendicular_frame | |
| 826 | |||
| 827 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 110163 times.
|
110163 | pure function cross(a, b) result(c) |
| 828 | real(dp), intent(in) :: a(3), b(3) | ||
| 829 | real(dp) :: c(3) | ||
| 830 | c = [a(2)*b(3) - a(3)*b(2), a(3)*b(1) - a(1)*b(3), & | ||
| 831 |
2/2✓ Branch 0 taken 330489 times.
✓ Branch 1 taken 110163 times.
|
440652 | a(1)*b(2) - a(2)*b(1)] |
| 832 | 110163 | end function cross | |
| 833 | ✗ | end module spectre_fo_hybrid | |
| 834 |