| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module alpha_lifetime_sub | ||
| 2 | |||
| 3 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 4 | |||
| 5 | implicit none | ||
| 6 | |||
| 7 | integer, parameter :: axis_rhs_count = 5 | ||
| 8 | |||
| 9 | type :: toroidal_regularization_context_t | ||
| 10 | real(dp) :: mu = 0.0_dp | ||
| 11 | real(dp) :: major_radius_cm = 1.0_dp | ||
| 12 | end type toroidal_regularization_context_t | ||
| 13 | |||
| 14 | ! The Burby--Ellison position map is retained only through first order in | ||
| 15 | ! gyroradius. Near a mirror point, its neglected second-order terms can put | ||
| 16 | ! the mapped state just below the original constant-momentum energy shell. | ||
| 17 | ! Clamp only that small mismatch; a larger negative radicand is a failed | ||
| 18 | ! regularized step and must fall through to the next recovery method. | ||
| 19 | real(dp), parameter :: toroidal_energy_shell_tolerance = 1.0e-5_dp | ||
| 20 | |||
| 21 | contains | ||
| 22 | |||
| 23 | 4 | subroutine bstar_parallel_factor(z, factor) | |
| 24 | use parmot_mod, only: ro0 | ||
| 25 | use magfie_sub, only: magfie | ||
| 26 | real(dp), intent(in) :: z(5) | ||
| 27 | real(dp), intent(out) :: factor | ||
| 28 | real(dp) :: bmod, sqrtg, bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 29 | |||
| 30 | 4 | call magfie(z(1:3), bmod, sqrtg, bder, hcovar, hctrvr, hcurl) | |
| 31 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 4 times.
|
16 | factor = 1.0_dp + z(4)*z(5)*(ro0/bmod)*dot_product(hcurl, hcovar) |
| 32 | 4 | end subroutine bstar_parallel_factor | |
| 33 | |||
| 34 | 2 | subroutine toroidal_regularized_position_map(x, ppar, direction, x_mapped, & | |
| 35 | ierr) | ||
| 36 | !> Leading near-identity map from Burby and Ellison, Phys. Plasmas 24, | ||
| 37 | !> 110703 (2017), Eq. (13). direction=1 maps the conventional guiding | ||
| 38 | !> centre to the toroidally regularized coordinate; direction=-1 | ||
| 39 | !> numerically inverts that first-order map. The toroidal component is | ||
| 40 | !> unchanged, which preserves field-period closure without an angular | ||
| 41 | !> seam. | ||
| 42 | real(dp), intent(in) :: x(3), ppar | ||
| 43 | integer, intent(in) :: direction | ||
| 44 | real(dp), intent(out) :: x_mapped(3) | ||
| 45 | integer, intent(out) :: ierr | ||
| 46 | real(dp), parameter :: derivative_step = 1.0e-6_dp | ||
| 47 | real(dp) :: shift(3), shift_plus(3), shift_minus(3) | ||
| 48 | real(dp) :: residual(2), residual_trial(2), jacobian(2, 2), delta(2) | ||
| 49 | real(dp) :: x_plus(3), x_minus(3), x_trial(3), determinant, damping | ||
| 50 | real(dp) :: residual_norm, trial_norm | ||
| 51 | integer :: iteration, line_search, shift_error | ||
| 52 | |||
| 53 | 2 | ierr = 2 | |
| 54 | 2 | x_mapped = x | |
| 55 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
4 | if (abs(direction) /= 1) return |
| 56 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | if (direction == 1) then |
| 57 | 1 | call toroidal_position_shift(x, ppar, shift, shift_error) | |
| 58 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (shift_error /= 0) return |
| 59 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
|
4 | x_mapped = x + shift |
| 60 | 1 | ierr = 0 | |
| 61 | 1 | return | |
| 62 | end if | ||
| 63 | |||
| 64 | ! Invert x_regular=x_standard+shift(x_standard) with a damped two- | ||
| 65 | ! dimensional Newton solve. Fixed-point inversion is not contractive | ||
| 66 | ! near the conventional B_parallel_star singularity. | ||
| 67 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | do iteration = 1, 20 |
| 68 | 3 | call toroidal_position_shift(x_mapped, ppar, shift, shift_error) | |
| 69 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if (shift_error /= 0) then |
| 70 | ✗ | x_mapped = x | |
| 71 | ✗ | return | |
| 72 | end if | ||
| 73 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 3 times.
|
9 | residual = x_mapped(1:2) + shift(1:2) - x(1:2) |
| 74 |
6/8✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 2 times.
|
12 | residual_norm = maxval(abs(residual)) |
| 75 |
8/10✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 9 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 6 times.
✓ Branch 7 taken 3 times.
✓ Branch 8 taken 1 time.
✓ Branch 9 taken 2 times.
|
15 | if (residual_norm <= & |
| 76 | 1.0e-12_dp*max(1.0_dp, maxval(abs(x)))) then | ||
| 77 | 1 | ierr = 0 | |
| 78 | 1 | return | |
| 79 | end if | ||
| 80 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
|
6 | do line_search = 1, 2 |
| 81 | 4 | x_plus = x_mapped | |
| 82 | 4 | x_minus = x_mapped | |
| 83 | 4 | x_plus(line_search) = x_plus(line_search) + derivative_step | |
| 84 | 4 | x_minus(line_search) = x_minus(line_search) - derivative_step | |
| 85 | call toroidal_position_shift(x_plus, ppar, shift_plus, & | ||
| 86 | 4 | shift_error) | |
| 87 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (shift_error /= 0) then |
| 88 | ✗ | x_mapped = x | |
| 89 | ✗ | return | |
| 90 | end if | ||
| 91 | call toroidal_position_shift(x_minus, ppar, shift_minus, & | ||
| 92 | 4 | shift_error) | |
| 93 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (shift_error /= 0) then |
| 94 | ✗ | x_mapped = x | |
| 95 | ✗ | return | |
| 96 | end if | ||
| 97 | jacobian(:, line_search) = & | ||
| 98 | (shift_plus(1:2) - shift_minus(1:2))/ & | ||
| 99 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 4 times.
|
12 | (2.0_dp*derivative_step) |
| 100 | jacobian(line_search, line_search) = & | ||
| 101 | 14 | jacobian(line_search, line_search) + 1.0_dp | |
| 102 | end do | ||
| 103 | determinant = jacobian(1, 1)*jacobian(2, 2) - & | ||
| 104 | 2 | jacobian(1, 2)*jacobian(2, 1) | |
| 105 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (abs(determinant) <= 100.0_dp*tiny(1.0_dp)) then |
| 106 | ✗ | x_mapped = x | |
| 107 | ✗ | return | |
| 108 | end if | ||
| 109 | delta(1) = (-jacobian(2, 2)*residual(1) + & | ||
| 110 | 2 | jacobian(1, 2)*residual(2))/determinant | |
| 111 | delta(2) = (jacobian(2, 1)*residual(1) - & | ||
| 112 | 2 | jacobian(1, 1)*residual(2))/determinant | |
| 113 | |||
| 114 | 2 | damping = 1.0_dp | |
| 115 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | do line_search = 1, 12 |
| 116 | 2 | x_trial = x_mapped | |
| 117 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
|
6 | x_trial(1:2) = x_trial(1:2) + damping*delta |
| 118 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (x_trial(1) < 0.0_dp) then |
| 119 | ✗ | damping = 0.5_dp*damping | |
| 120 | ✗ | cycle | |
| 121 | end if | ||
| 122 | call toroidal_position_shift(x_trial, ppar, shift_plus, & | ||
| 123 | 2 | shift_error) | |
| 124 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (shift_error /= 0) then |
| 125 | ✗ | damping = 0.5_dp*damping | |
| 126 | ✗ | cycle | |
| 127 | end if | ||
| 128 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
|
6 | residual_trial = x_trial(1:2) + shift_plus(1:2) - x(1:2) |
| 129 |
6/8✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 3 times.
✓ Branch 7 taken 1 time.
|
8 | trial_norm = maxval(abs(residual_trial)) |
| 130 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (trial_norm < residual_norm) then |
| 131 | 2 | x_mapped = x_trial | |
| 132 | exit | ||
| 133 | end if | ||
| 134 | 2 | damping = 0.5_dp*damping | |
| 135 | end do | ||
| 136 | 3 | if (line_search > 12) then | |
| 137 | ✗ | x_mapped = x | |
| 138 | ✗ | return | |
| 139 | end if | ||
| 140 | end do | ||
| 141 | ✗ | x_mapped = x | |
| 142 | 4 | end subroutine toroidal_regularized_position_map | |
| 143 | |||
| 144 | 14 | subroutine toroidal_position_shift(x, ppar, shift, ierr) | |
| 145 | use finite_check, only: is_finite | ||
| 146 | use parmot_mod, only: ro0 | ||
| 147 | use magfie_sub, only: magfie_refcoords | ||
| 148 | real(dp), intent(in) :: x(3), ppar | ||
| 149 | real(dp), intent(out) :: shift(3) | ||
| 150 | integer, intent(out) :: ierr | ||
| 151 | real(dp) :: bmod, sqrtg, bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 152 | real(dp) :: bphi | ||
| 153 | |||
| 154 | 14 | ierr = 2 | |
| 155 | 14 | shift = 0.0_dp | |
| 156 | 14 | call magfie_refcoords(x, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) | |
| 157 | 14 | bphi = bmod*hctrvr(3) | |
| 158 |
1/2✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
|
14 | if (abs(sqrtg*bphi) <= 100.0_dp*tiny(1.0_dp)) return |
| 159 | 14 | shift(1) = ro0*ppar*hcovar(2)/(sqrtg*bphi) | |
| 160 | 14 | shift(2) = -ro0*ppar*hcovar(1)/(sqrtg*bphi) | |
| 161 |
4/6✓ Branch 0 taken 42 times.
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 42 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 14 times.
|
56 | if (.not. all(is_finite(shift))) then |
| 162 | ✗ | shift = 0.0_dp | |
| 163 | ✗ | return | |
| 164 | end if | ||
| 165 | 14 | ierr = 0 | |
| 166 | 14 | end subroutine toroidal_position_shift | |
| 167 | |||
| 168 | ✗ | subroutine toroidal_regularized_rhs(tau, y, dy, context) | |
| 169 | 14 | use parmot_mod, only: rmu, ro0 | |
| 170 | use magfie_sub, only: magfie_refcoords | ||
| 171 | real(dp), intent(in) :: tau, y(:) | ||
| 172 | real(dp), intent(out) :: dy(:) | ||
| 173 | class(*), intent(in) :: context | ||
| 174 | real(dp), parameter :: derivative_step = 1.0e-5_dp | ||
| 175 | real(dp) :: bmod, sqrtg, bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 176 | real(dp) :: bp, gp, bdp(3), hcp(3), htp(3), hxp(3) | ||
| 177 | real(dp) :: bm, gm, bdm(3), hcm(3), htm(3), hxm(3) | ||
| 178 | real(dp) :: xwork(3), dq(3), dkinetic(3), q, u, ppar, vpa | ||
| 179 | real(dp) :: drift_scale, gamma | ||
| 180 | real(dp) :: h | ||
| 181 | integer :: i | ||
| 182 | |||
| 183 | select type (ctx => context) | ||
| 184 | type is (toroidal_regularization_context_t) | ||
| 185 | call magfie_refcoords(y(1:3), bmod, sqrtg, bder, hcovar, & | ||
| 186 | ✗ | hctrvr, hcurl) | |
| 187 | ✗ | q = ctx%major_radius_cm*hctrvr(3) | |
| 188 | ✗ | if (abs(q) <= 100.0_dp*tiny(1.0_dp)) then | |
| 189 | ✗ | dy = huge(1.0_dp) | |
| 190 | return | ||
| 191 | end if | ||
| 192 | |||
| 193 | ✗ | do i = 1, 3 | |
| 194 | h = derivative_step | ||
| 195 | ✗ | xwork = y(1:3) | |
| 196 | ✗ | xwork(i) = xwork(i) + h | |
| 197 | ✗ | call magfie_refcoords(xwork, bp, gp, bdp, hcp, htp, hxp) | |
| 198 | ✗ | xwork(i) = y(i) - h | |
| 199 | ✗ | call magfie_refcoords(xwork, bm, gm, bdm, hcm, htm, hxm) | |
| 200 | ✗ | dq(i) = ctx%major_radius_cm*(htp(3) - htm(3))/(2.0_dp*h) | |
| 201 | end do | ||
| 202 | |||
| 203 | ✗ | u = y(4) | |
| 204 | ✗ | ppar = q*u | |
| 205 | gamma = sqrt(1.0_dp + & | ||
| 206 | ✗ | 2.0_dp*(ppar**2 + 2.0_dp*ctx%mu*bmod)/rmu) | |
| 207 | ✗ | vpa = ppar/gamma | |
| 208 | ✗ | dkinetic = (u*u*q*dq + ctx%mu*bmod*bder)/gamma | |
| 209 | ✗ | drift_scale = ro0/(sqrtg*bmod*hctrvr(3)) | |
| 210 | ✗ | dy(1) = vpa*hctrvr(1) - drift_scale*dkinetic(2) | |
| 211 | ✗ | dy(2) = vpa*hctrvr(2) + drift_scale*dkinetic(1) | |
| 212 | ✗ | dy(3) = vpa*hctrvr(3) | |
| 213 | ✗ | dy(4) = -dot_product(hctrvr, dkinetic)/q | |
| 214 | class default | ||
| 215 | ✗ | dy = huge(1.0_dp) | |
| 216 | end select | ||
| 217 | ✗ | end subroutine toroidal_regularized_rhs | |
| 218 | |||
| 219 | 4 | subroutine orbit_timestep_toroidal_regularized(z, dtau, relerr, ierr, & | |
| 220 | step_limit) | ||
| 221 | ✗ | use finite_check, only: is_finite | |
| 222 | use odeint_allroutines_sub, only: odeint_allroutines, & | ||
| 223 | odeint_has_failed | ||
| 224 | use chamb_sub, only: chamb_can | ||
| 225 | use magfie_sub, only: magfie_refcoords, has_magfie_refcoords_field | ||
| 226 | use field_can_mod, only: integ_to_ref, ref_to_integ | ||
| 227 | use new_vmec_stuff_mod, only: rmajor | ||
| 228 | real(dp), intent(inout) :: z(5) | ||
| 229 | real(dp), intent(in) :: dtau, relerr | ||
| 230 | integer, intent(out) :: ierr | ||
| 231 | integer, intent(in) :: step_limit | ||
| 232 | type(toroidal_regularization_context_t) :: context | ||
| 233 | real(dp) :: z_initial(5), y(4), chamber_y(2) | ||
| 234 | real(dp) :: x_standard(3), x_regular(3) | ||
| 235 | real(dp) :: bmod, sqrtg, bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 236 | real(dp) :: q, ppar, ppar_raw, ppar_squared | ||
| 237 | integer :: map_error | ||
| 238 | |||
| 239 | 4 | z_initial = z | |
| 240 | 4 | if (.not. has_magfie_refcoords_field()) then | |
| 241 | 4 | ierr = 2 | |
| 242 | 4 | return | |
| 243 | end if | ||
| 244 | ✗ | call integ_to_ref(z(1:3), x_standard) | |
| 245 | call magfie_refcoords(x_standard, bmod, sqrtg, bder, hcovar, hctrvr, & | ||
| 246 | ✗ | hcurl) | |
| 247 | ✗ | context%major_radius_cm = 100.0_dp*rmajor | |
| 248 | ✗ | context%mu = 0.5_dp*z(4)**2*(1.0_dp - z(5)**2)/bmod | |
| 249 | ✗ | ppar = z(4)*z(5) | |
| 250 | call toroidal_regularized_position_map(x_standard, ppar, 1, x_regular, & | ||
| 251 | ✗ | map_error) | |
| 252 | ✗ | if (map_error /= 0) then | |
| 253 | ✗ | ierr = 2 | |
| 254 | ✗ | return | |
| 255 | end if | ||
| 256 | ✗ | y(1:3) = x_regular | |
| 257 | ✗ | call magfie_refcoords(y(1:3), bmod, sqrtg, bder, hcovar, hctrvr, hcurl) | |
| 258 | ✗ | q = context%major_radius_cm*hctrvr(3) | |
| 259 | ✗ | if (abs(q) <= 100.0_dp*tiny(1.0_dp)) then | |
| 260 | ✗ | ierr = 2 | |
| 261 | ✗ | return | |
| 262 | end if | ||
| 263 | ✗ | y(4) = ppar/q | |
| 264 | |||
| 265 | call odeint_allroutines(y, 4, context, 0.0_dp, dtau, relerr, & | ||
| 266 | ✗ | toroidal_regularized_rhs, step_limit=step_limit) | |
| 267 | ✗ | if (odeint_has_failed() .or. .not. all(is_finite(y))) then | |
| 268 | ✗ | z = z_initial | |
| 269 | ✗ | ierr = 2 | |
| 270 | ✗ | return | |
| 271 | end if | ||
| 272 | |||
| 273 | ✗ | call magfie_refcoords(y(1:3), bmod, sqrtg, bder, hcovar, hctrvr, hcurl) | |
| 274 | ✗ | q = context%major_radius_cm*hctrvr(3) | |
| 275 | ✗ | ppar_raw = q*y(4) | |
| 276 | call toroidal_regularized_position_map(y(1:3), ppar_raw, -1, & | ||
| 277 | ✗ | x_standard, map_error) | |
| 278 | ✗ | if (map_error /= 0) then | |
| 279 | ✗ | z = z_initial | |
| 280 | ✗ | ierr = 2 | |
| 281 | ✗ | return | |
| 282 | end if | ||
| 283 | call magfie_refcoords(x_standard, bmod, sqrtg, bder, hcovar, hctrvr, & | ||
| 284 | ✗ | hcurl) | |
| 285 | ✗ | ppar_squared = z_initial(4)**2 - 2.0_dp*context%mu*bmod | |
| 286 | ✗ | if (ppar_squared < & | |
| 287 | -toroidal_energy_shell_tolerance*z_initial(4)**2) then | ||
| 288 | ✗ | z = z_initial | |
| 289 | ✗ | ierr = 2 | |
| 290 | ✗ | return | |
| 291 | end if | ||
| 292 | ✗ | ppar = sign(sqrt(max(ppar_squared, 0.0_dp)), ppar_raw) | |
| 293 | ✗ | call ref_to_integ(x_standard, z(1:3)) | |
| 294 | ✗ | z(4) = z_initial(4) | |
| 295 | ✗ | z(5) = ppar/z(4) | |
| 296 | ✗ | chamber_y = z(1:2) | |
| 297 | ✗ | call chamb_can(chamber_y, z(3), ierr) | |
| 298 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
8 | end subroutine orbit_timestep_toroidal_regularized |
| 299 | |||
| 300 | !ToDo make module from all global things like this one | ||
| 301 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 302 | ! | ||
| 303 | 34974647 | subroutine elefie_can(x,derphi) | |
| 304 | ! | ||
| 305 | ! Computes the derivatives of the electrostatic potential over | ||
| 306 | ! coordinates (covariant electric field). Potential is normalized | ||
| 307 | ! to reference temperature : phinor = e phi / T. | ||
| 308 | ! | ||
| 309 | ! Input parameters: | ||
| 310 | ! formal: x - array of coordinates | ||
| 311 | ! | ||
| 312 | ! Output parameters: | ||
| 313 | ! formal: derphi - array of derivatives | ||
| 314 | ! | ||
| 315 | real(dp), intent(in) :: x(3) | ||
| 316 | real(dp), intent(out) :: derphi(3) | ||
| 317 | ! | ||
| 318 | 34974647 | derphi=0.d0 | |
| 319 | ! | ||
| 320 | 4 | end subroutine elefie_can | |
| 321 | ! | ||
| 322 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 323 | ! | ||
| 324 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 34974647 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 34974647 times.
|
34974647 | subroutine velo_can(tau,z,vz) |
| 325 | ! | ||
| 326 | ! | ||
| 327 | ! Computes the components of the 5-D velocity - vz | ||
| 328 | ! for given set of phase space coordinates - z. | ||
| 329 | ! | ||
| 330 | ! Warning !!! The dimensionless time is used (see below) | ||
| 331 | ! | ||
| 332 | ! Order of the coordinates is the following: | ||
| 333 | ! z(i) = x(i) for i=1,2,3 - spatial coordinates with real | ||
| 334 | ! dimension of the general covariant | ||
| 335 | ! space coordinate system | ||
| 336 | ! z(4) = p - momentum module normalized to | ||
| 337 | ! thermal momentum and sqrt(2); | ||
| 338 | ! z(5) = alambd - cosine of the pitch-angle | ||
| 339 | ! | ||
| 340 | ! Input parameters: | ||
| 341 | ! formal: tau - dimensionless time: tau=sqrt(2*T/m)*t | ||
| 342 | ! z - see above | ||
| 343 | ! common: rmu - inverse relativistic temperature | ||
| 344 | ! ro0 - Larmor radius for the reference | ||
| 345 | ! magnetic field and temperature: | ||
| 346 | ! ro0=sqrt(2*T/m)/(e*B_ref/m*c) | ||
| 347 | ! Output parameters: | ||
| 348 | ! formal: vz - see above | ||
| 349 | ! | ||
| 350 | ! Called routines: magfie, elefie_can | ||
| 351 | ! | ||
| 352 | use parmot_mod, only : rmu,ro0 | ||
| 353 | use magfie_sub, only : magfie | ||
| 354 | use field_can_base, only: count_field_evaluation | ||
| 355 | ! | ||
| 356 | implicit none | ||
| 357 | ! | ||
| 358 | integer :: i | ||
| 359 | ! | ||
| 360 | real(dp), intent(in) :: tau | ||
| 361 | real(dp), intent(in) :: z(:) | ||
| 362 | real(dp), intent(out) :: vz(:) | ||
| 363 | real(dp) x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl | ||
| 364 | real(dp) derphi | ||
| 365 | real(dp) p,alambd,p2,ovmu,gamma2,gamma,ppar,vpa,coala | ||
| 366 | real(dp) rmumag,rovsqg,rosqgb,rovbm | ||
| 367 | real(dp) a_phi,a_b,a_c,hstar | ||
| 368 | real(dp) s_hc,hpstar,phidot,blodot,bra | ||
| 369 | ! | ||
| 370 | dimension x(3),bder(3),hcovar(3),hctrvr(3),hcurl(3) | ||
| 371 | dimension derphi(3) | ||
| 372 | dimension a_phi(3),a_b(3),a_c(3),hstar(3) | ||
| 373 | ! | ||
| 374 |
2/2✓ Branch 0 taken 104923941 times.
✓ Branch 1 taken 34974647 times.
|
139898588 | do 1 i=1,3 |
| 375 | 104923941 | x(i)=z(i) | |
| 376 | 34974647 | 1 continue | |
| 377 | ! | ||
| 378 | ! in magfie: x(i) - set of 3 curvilinear space coordinates (input) | ||
| 379 | ! bmod - dimensionless magnetic field module: bmod=B/B_ref | ||
| 380 | ! sqrtg - Jacobian of space coordinates (square root of | ||
| 381 | ! metric tensor | ||
| 382 | ! bder - derivatives of logarithm of bmod over space coords | ||
| 383 | ! (covariant vector) | ||
| 384 | ! hcovar - covariant components of the unit vector along | ||
| 385 | ! the magnetic field | ||
| 386 | ! hctrvr - contravariant components of the unit vector along | ||
| 387 | ! the magnetic field | ||
| 388 | ! hcurl - contravariant components of the curl of this vector | ||
| 389 | ! | ||
| 390 | 34974647 | call magfie(x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl) | |
| 391 | ! The magfie path carries its own counter. The canonical backends | ||
| 392 | ! bump n_field_evaluations inside field_can, but this drift-kinetic | ||
| 393 | ! right-hand side goes straight to magfie and was previously | ||
| 394 | ! invisible, so any cost comparison between the RK and symplectic | ||
| 395 | ! paths was counting only one of them. magfie supplies first | ||
| 396 | ! derivatives only, so it counts as a d1 evaluation. | ||
| 397 | 34974647 | call count_field_evaluation(0) | |
| 398 | ! | ||
| 399 | ! in elefie: x(i) - space coords (input, see above) | ||
| 400 | ! derphi - derivatives of the dimensionless electric potential | ||
| 401 | ! phihat=e*phi/T over space coords (covar. vector) | ||
| 402 | ! | ||
| 403 | 34974647 | call elefie_can(x,derphi) | |
| 404 | ! | ||
| 405 | 34974647 | p=z(4) | |
| 406 | 34974647 | alambd=z(5) | |
| 407 | ! | ||
| 408 | 34974647 | p2=p*p | |
| 409 | 34974647 | ovmu=2.d0/rmu | |
| 410 | 34974647 | gamma2=p2*ovmu+1.d0 | |
| 411 | 34974647 | gamma=dsqrt(gamma2) | |
| 412 | 34974647 | ppar=p*alambd | |
| 413 | ! vpa - dimensionless parallel velocity: vpa=v_parallel/sqrt(2*T/m) | ||
| 414 | 34974647 | vpa=ppar/gamma | |
| 415 | 34974647 | coala=(1.d0-alambd**2) | |
| 416 | ! rmumag - magnetic moment | ||
| 417 | 34974647 | rmumag=.5d0*p2*coala/bmod | |
| 418 | ! | ||
| 419 | 34974647 | rovsqg=ro0/sqrtg | |
| 420 | 34974647 | rosqgb=.5d0*rovsqg/bmod | |
| 421 | 34974647 | rovbm=ro0/bmod | |
| 422 | ! | ||
| 423 | 34974647 | a_phi(1)=(hcovar(2)*derphi(3)-hcovar(3)*derphi(2))*rosqgb | |
| 424 | 34974647 | a_b(1)=(hcovar(2)*bder(3)-hcovar(3)*bder(2))*rovsqg | |
| 425 | 34974647 | a_phi(2)=(hcovar(3)*derphi(1)-hcovar(1)*derphi(3))*rosqgb | |
| 426 | 34974647 | a_b(2)=(hcovar(3)*bder(1)-hcovar(1)*bder(3))*rovsqg | |
| 427 | 34974647 | a_phi(3)=(hcovar(1)*derphi(2)-hcovar(2)*derphi(1))*rosqgb | |
| 428 | 34974647 | a_b(3)=(hcovar(1)*bder(2)-hcovar(2)*bder(1))*rovsqg | |
| 429 | ! | ||
| 430 | 34974647 | s_hc=0.d0 | |
| 431 |
2/2✓ Branch 0 taken 104923941 times.
✓ Branch 1 taken 34974647 times.
|
139898588 | do i=1,3 |
| 432 | 104923941 | a_c(i)=hcurl(i)*rovbm | |
| 433 | 104923941 | s_hc=s_hc+a_c(i)*hcovar(i) | |
| 434 | 139898588 | hstar(i)=hctrvr(i)+ppar*a_c(i) | |
| 435 | enddo | ||
| 436 | 34974647 | hpstar=1.d0+ppar*s_hc | |
| 437 | ! | ||
| 438 | ! velocities in the coordinate space | ||
| 439 | ! | ||
| 440 | ! phidot - derivative of the dmls el. potential over dmls time | ||
| 441 | ! blodot - derivative of the logarith of the mag. field module over dmls time | ||
| 442 | 34974647 | phidot=0.d0 | |
| 443 | 34974647 | blodot=0.d0 | |
| 444 |
2/2✓ Branch 0 taken 104923941 times.
✓ Branch 1 taken 34974647 times.
|
139898588 | do i=1,3 |
| 445 | 104923941 | bra=vpa*hstar(i)+a_phi(i)+a_b(i)*rmumag/gamma | |
| 446 | 104923941 | vz(i)=bra/hpstar | |
| 447 | 104923941 | phidot=phidot+vz(i)*derphi(i) | |
| 448 | 139898588 | blodot=blodot+vz(i)*bder(i) | |
| 449 | enddo | ||
| 450 | ! | ||
| 451 | ! velocities in the phase space | ||
| 452 | ! | ||
| 453 | 34974647 | vz(4)=-0.5d0*gamma*phidot/p | |
| 454 | ! vz(5)=coala/(alambd+dsign(1.d-32,alambd))*(vz(4)/p-0.5d0*blodot) | ||
| 455 | vz(5)=-(0.5d0*coala/hpstar)*(sum(hstar*derphi)/p & | ||
| 456 |
6/6✓ Branch 0 taken 34974647 times.
✓ Branch 1 taken 104923941 times.
✓ Branch 2 taken 34974647 times.
✓ Branch 3 taken 104923941 times.
✓ Branch 4 taken 104923941 times.
✓ Branch 5 taken 34974647 times.
|
419695764 | + p*sum(hstar*bder)/gamma+alambd*sum(a_phi*bder)) |
| 457 | ! | ||
| 458 | 34974647 | end subroutine velo_can | |
| 459 | ! | ||
| 460 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 461 | ! | ||
| 462 | ✗ | subroutine orbit_timestep_can(z,dtau,dtaumin,relerr,ierr) | |
| 463 | 34974647 | use diag_mod, only : dodiag | |
| 464 | use odeint_allroutines_sub, only : odeint_allroutines | ||
| 465 | use chamb_sub, only : chamb_can | ||
| 466 | ! | ||
| 467 | implicit none | ||
| 468 | ! | ||
| 469 | integer, parameter :: ndim=5, nstepmax=1000000 | ||
| 470 | ! | ||
| 471 | integer :: ierr | ||
| 472 | real(dp) :: dtau,dtaumin,phi,tau1,tau2 | ||
| 473 | ! | ||
| 474 | real(dp), dimension(2) :: y | ||
| 475 | real(dp), dimension(ndim) :: z | ||
| 476 | real(dp) :: relerr | ||
| 477 | ! | ||
| 478 | ✗ | if(abs(dtaumin*nstepmax).le.abs(dtau)) then | |
| 479 | ✗ | ierr=2 | |
| 480 | ✗ | print *,'orbit_timestep: number of steps exceeds nstepmax' | |
| 481 | ✗ | return | |
| 482 | endif | ||
| 483 | ! | ||
| 484 | ✗ | ierr=0 | |
| 485 | ✗ | y(1)=z(1) | |
| 486 | ✗ | y(2)=z(2) | |
| 487 | ✗ | phi=z(3) | |
| 488 | ! | ||
| 489 | ✗ | call chamb_can(y,phi,ierr) | |
| 490 | ! | ||
| 491 | ✗ | if(ierr.eq.1) return | |
| 492 | ✗ | tau1=0.d0 | |
| 493 | ✗ | tau2=dtaumin | |
| 494 | ! | ||
| 495 | ✗ | do while((dtau>0 .and. (tau2 .lt. dtau)) .or. (dtau<0 .and. (tau2 .gt. dtau))) | |
| 496 | ! | ||
| 497 | ✗ | call odeint_allroutines(z,ndim,tau1,tau2,relerr,velo_can) | |
| 498 | ! | ||
| 499 | ✗ | if(dodiag) write (123,*) tau2,z | |
| 500 | ✗ | y(1)=z(1) | |
| 501 | ✗ | y(2)=z(2) | |
| 502 | ✗ | phi=z(3) | |
| 503 | ! | ||
| 504 | ✗ | call chamb_can(y,phi,ierr) | |
| 505 | ! | ||
| 506 | ✗ | if(ierr.eq.1) return | |
| 507 | ✗ | tau1=tau2 | |
| 508 | ✗ | tau2=tau2+dtaumin | |
| 509 | enddo | ||
| 510 | ! | ||
| 511 | ✗ | tau2=dtau | |
| 512 | ! | ||
| 513 | ✗ | call odeint_allroutines(z,ndim,tau1,tau2,relerr,velo_can) | |
| 514 | ! | ||
| 515 | ✗ | y(1)=z(1) | |
| 516 | ✗ | y(2)=z(2) | |
| 517 | ✗ | phi=z(3) | |
| 518 | ! | ||
| 519 | ✗ | call chamb_can(y,phi,ierr) | |
| 520 | ! | ||
| 521 | ✗ | if(ierr.eq.1) return | |
| 522 | ! | ||
| 523 | end subroutine orbit_timestep_can | ||
| 524 | ! | ||
| 525 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 526 | ! | ||
| 527 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 13243681 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13243681 times.
|
13243681 | subroutine rhs_mflint_can(phi,y,dery) |
| 528 | ! | ||
| 529 | ! Computes the right hand side of the magnetic field line equation for | ||
| 530 | ! the integration over the toroidal angle, subintegrand for the flux tube | ||
| 531 | ! volume $1/B^\varphi$ and subintegrants for Boozer $B_{00}$ computation | ||
| 532 | ! $B^2/B^\varphi$ and $B^3/B^\varphi$ - dery | ||
| 533 | ! | ||
| 534 | ! Oder of coordinates in magfie is the following: x(1)=R (big radius), | ||
| 535 | ! x(2)=phi (toroidal angle), x(3)=z (altitude). | ||
| 536 | ! | ||
| 537 | ! Input parameters: | ||
| 538 | ! formal: | ||
| 539 | ! y(1:2) - coordinates in the poloidal plane (phi=const) | ||
| 540 | ! y(3:5) - integrals | ||
| 541 | ! Output parameters: | ||
| 542 | ! formal: | ||
| 543 | ! dery(1:5) - vector of the right-hand side | ||
| 544 | ! Called routines: magfie_vmec, magfie_can, magfie_boozer | ||
| 545 | ! | ||
| 546 | use magfie_sub, only : magfie | ||
| 547 | ! | ||
| 548 | real(dp), intent(in) :: phi | ||
| 549 | real(dp), intent(in) :: y(:) | ||
| 550 | real(dp), intent(out) :: dery(:) | ||
| 551 | real(dp) x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl | ||
| 552 | dimension x(3),bder(3),hcovar(3),hctrvr(3),hcurl(3) | ||
| 553 | ! | ||
| 554 | 13243681 | x(1)=y(1) | |
| 555 | 13243681 | x(2)=y(2) | |
| 556 | 13243681 | x(3)=phi | |
| 557 | ! | ||
| 558 | 13243681 | call magfie(x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl) | |
| 559 | ! | ||
| 560 | 13243681 | dery(1)=hctrvr(1)/hctrvr(3) | |
| 561 | 13243681 | dery(2)=hctrvr(2)/hctrvr(3) | |
| 562 | 13243681 | dery(3)=1.d0/(bmod*hctrvr(3)) | |
| 563 | 13243681 | dery(4)=bmod/hctrvr(3) | |
| 564 | 13243681 | dery(5)=bmod**2/hctrvr(3) | |
| 565 | ! | ||
| 566 | 13243681 | end subroutine rhs_mflint_can | |
| 567 | ! | ||
| 568 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 569 | !Integrates along a magnetic field line to generate equidistant homogeneous volumetric distribution on flux surface | ||
| 570 | 57 | subroutine integrate_mfl_can(npoi,dphi,rbeg,phibeg,zbeg, & | |
| 571 | 57 | xstart,bstart,volstart,bmod00,ierr) | |
| 572 | ! | ||
| 573 | 13243681 | use odeint_allroutines_sub, only : odeint_allroutines | |
| 574 | use chamb_sub, only : chamb_can | ||
| 575 | use magfie_sub, only : magfie | ||
| 576 | ! | ||
| 577 | implicit none | ||
| 578 | ! | ||
| 579 | real(dp), intent(out) :: bmod00 | ||
| 580 | ! | ||
| 581 | integer, parameter :: ndim=5 | ||
| 582 | real(dp), parameter :: relerr=1d-10 | ||
| 583 | integer :: i | ||
| 584 | integer, intent(in) :: npoi, ierr | ||
| 585 | real(dp) :: phi,phiold | ||
| 586 | real(dp) , intent(in):: dphi,rbeg,phibeg,zbeg | ||
| 587 | real(dp), dimension(3,npoi), intent(out) :: xstart | ||
| 588 | real(dp), dimension(npoi), intent(out) :: bstart,volstart | ||
| 589 | real(dp), dimension(ndim) :: y | ||
| 590 | ! | ||
| 591 | real(dp) x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl | ||
| 592 | dimension x(3),bder(3),hcovar(3),hctrvr(3),hcurl(3) | ||
| 593 | ! | ||
| 594 | 57 | phi=phibeg | |
| 595 | 57 | y(1)=rbeg | |
| 596 | 57 | y(2)=zbeg | |
| 597 |
2/2✓ Branch 0 taken 171 times.
✓ Branch 1 taken 57 times.
|
228 | y(3:5)=0.d0 |
| 598 | ! | ||
| 599 | 57 | call chamb_can(y(1:2),phi,ierr) | |
| 600 | ! | ||
| 601 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 57 times.
|
57 | if(ierr.eq.1) return |
| 602 | 57 | x(1)=y(1) | |
| 603 | 57 | x(2)=y(2) | |
| 604 | 57 | x(3)=phi | |
| 605 | ! | ||
| 606 | 57 | call magfie(x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl) | |
| 607 | ! | ||
| 608 |
2/2✓ Branch 0 taken 171 times.
✓ Branch 1 taken 57 times.
|
228 | xstart(:,1)=x |
| 609 | 57 | bstart(1)=bmod | |
| 610 | 57 | volstart(1)=y(3) | |
| 611 | ! | ||
| 612 |
2/2✓ Branch 0 taken 2200223 times.
✓ Branch 1 taken 57 times.
|
2200280 | do i=2,npoi |
| 613 | 2200223 | phiold=phi | |
| 614 | 2200223 | phi=phiold+dphi | |
| 615 | ! | ||
| 616 | 2200223 | call odeint_allroutines(y,ndim,phiold,phi,relerr,rhs_mflint_can) | |
| 617 | ! | ||
| 618 | 2200223 | call chamb_can(y(1:2),phi,ierr) | |
| 619 | ! | ||
| 620 | if(ierr.eq.1) return | ||
| 621 | 2200223 | x(1)=y(1) | |
| 622 | 2200223 | x(2)=y(2) | |
| 623 | 2200223 | x(3)=phi | |
| 624 | ! | ||
| 625 | 2200223 | call magfie(x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl) | |
| 626 | ! | ||
| 627 |
2/2✓ Branch 0 taken 6600669 times.
✓ Branch 1 taken 2200223 times.
|
8800892 | xstart(:,i)=x |
| 628 | 2200223 | bstart(i)=bmod | |
| 629 | 4400503 | volstart(i)=y(3) | |
| 630 | enddo | ||
| 631 | ! | ||
| 632 |
2/2✓ Branch 0 taken 2200280 times.
✓ Branch 1 taken 57 times.
|
2200337 | volstart=volstart/volstart(npoi) |
| 633 | 57 | bmod00=y(5)/y(4) | |
| 634 | ! | ||
| 635 | 57 | end subroutine integrate_mfl_can | |
| 636 | ! | ||
| 637 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 638 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 639 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 640 | ! | ||
| 641 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 642 | ! | ||
| 643 | ✗ | subroutine velo_axis(tau,z_axis,vz_axis) | |
| 644 | real(dp), intent(in) :: tau | ||
| 645 | real(dp), intent(in) :: z_axis(:) | ||
| 646 | real(dp), intent(out) :: vz_axis(:) | ||
| 647 | real(dp), parameter :: legacy_s_floor = 1.d-8 | ||
| 648 | |||
| 649 | ✗ | call velo_axis_regularized(tau, z_axis, vz_axis, legacy_s_floor) | |
| 650 | 57 | end subroutine velo_axis | |
| 651 | |||
| 652 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 534744 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 534744 times.
|
534744 | subroutine velo_axis_regularized(tau, z_axis, vz_axis, s_floor) |
| 653 | ! s_floor uses native radial units; reconstruction uses Cartesian R**2. | ||
| 654 | use magfie_sub, only: axis_radial_coordinate | ||
| 655 | real(dp), intent(in) :: tau, z_axis(:), s_floor | ||
| 656 | real(dp), intent(out) :: vz_axis(:) | ||
| 657 | logical :: radial_is_rho | ||
| 658 | integer :: ierr | ||
| 659 | |||
| 660 | 534744 | call axis_radial_coordinate(radial_is_rho, ierr) | |
| 661 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 534744 times.
|
534744 | if (ierr /= 0) error stop 'unsupported reference axis coordinate' |
| 662 | call evaluate_axis_rhs(tau, z_axis, & | ||
| 663 | axis_native_radius_squared(s_floor, radial_is_rho), vz_axis, & | ||
| 664 | 534744 | radial_is_rho) | |
| 665 | 534744 | end subroutine velo_axis_regularized | |
| 666 | |||
| 667 | 4294056 | pure real(dp) function axis_native_radius_squared(q, radial_is_rho) result(r2) | |
| 668 | real(dp), intent(in) :: q | ||
| 669 | logical, intent(in) :: radial_is_rho | ||
| 670 | |||
| 671 | 4294056 | r2 = q | |
| 672 |
2/2✓ Branch 0 taken 1122 times.
✓ Branch 1 taken 4292934 times.
|
4294056 | if (radial_is_rho) r2 = q*q |
| 673 | 534744 | end function axis_native_radius_squared | |
| 674 | |||
| 675 | 89198 | subroutine axis_polar_to_cartesian(zpolar, zaxis, radial_is_rho) | |
| 676 | use axis_map_generated, only: axis_flux_forward, axis_rho_forward | ||
| 677 | real(dp), intent(in) :: zpolar(5) | ||
| 678 | real(dp), intent(out) :: zaxis(5) | ||
| 679 | logical, intent(in) :: radial_is_rho | ||
| 680 | |||
| 681 |
2/2✓ Branch 0 taken 40 times.
✓ Branch 1 taken 89158 times.
|
89198 | if (radial_is_rho) then |
| 682 | call axis_rho_forward(max(zpolar(1), 0.0_dp), zpolar(2), & | ||
| 683 | 40 | zaxis(1), zaxis(2)) | |
| 684 | else | ||
| 685 | call axis_flux_forward(max(zpolar(1), 0.0_dp), zpolar(2), & | ||
| 686 | 89158 | zaxis(1), zaxis(2)) | |
| 687 | end if | ||
| 688 |
2/2✓ Branch 0 taken 267594 times.
✓ Branch 1 taken 89198 times.
|
356792 | zaxis(3:5) = zpolar(3:5) |
| 689 | 89198 | end subroutine axis_polar_to_cartesian | |
| 690 | |||
| 691 | 624274 | subroutine axis_cartesian_to_polar(zaxis, zpolar, radial_is_rho) | |
| 692 | use axis_map_generated, only: axis_flux_inverse, axis_rho_inverse | ||
| 693 | real(dp), intent(in) :: zaxis(5) | ||
| 694 | real(dp), intent(out) :: zpolar(5) | ||
| 695 | logical, intent(in) :: radial_is_rho | ||
| 696 | |||
| 697 |
2/2✓ Branch 0 taken 208 times.
✓ Branch 1 taken 624066 times.
|
624274 | if (radial_is_rho) then |
| 698 | 208 | call axis_rho_inverse(zaxis(1), zaxis(2), zpolar(1)) | |
| 699 | else | ||
| 700 | 624066 | call axis_flux_inverse(zaxis(1), zaxis(2), zpolar(1)) | |
| 701 | end if | ||
| 702 | 624274 | zpolar(2) = atan2(zaxis(2), zaxis(1)) | |
| 703 |
2/2✓ Branch 0 taken 1872822 times.
✓ Branch 1 taken 624274 times.
|
2497096 | zpolar(3:5) = zaxis(3:5) |
| 704 | 624274 | end subroutine axis_cartesian_to_polar | |
| 705 | |||
| 706 | 535136 | subroutine axis_polar_velocity(zaxis, vpolar, vaxis, radial_is_rho) | |
| 707 | use axis_map_generated, only: axis_flux_velocity, axis_rho_velocity | ||
| 708 | real(dp), intent(in) :: zaxis(5), vpolar(5) | ||
| 709 | real(dp), intent(out) :: vaxis(5) | ||
| 710 | logical, intent(in) :: radial_is_rho | ||
| 711 | |||
| 712 |
2/2✓ Branch 0 taken 198 times.
✓ Branch 1 taken 534938 times.
|
535136 | if (radial_is_rho) then |
| 713 | call axis_rho_velocity(zaxis(1), zaxis(2), vpolar(1), vpolar(2), & | ||
| 714 | 198 | vaxis(1), vaxis(2)) | |
| 715 | else | ||
| 716 | call axis_flux_velocity(zaxis(1), zaxis(2), vpolar(1), vpolar(2), & | ||
| 717 | 534938 | vaxis(1), vaxis(2)) | |
| 718 | end if | ||
| 719 |
2/2✓ Branch 0 taken 1605408 times.
✓ Branch 1 taken 535136 times.
|
2140544 | vaxis(3:5) = vpolar(3:5) |
| 720 | 535136 | end subroutine axis_polar_velocity | |
| 721 | |||
| 722 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 534744 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 534744 times.
|
534744 | subroutine evaluate_axis_rhs(tau, z_axis, sample_s, vz_axis, radial_is_rho) |
| 723 | real(dp), intent(in) :: tau, z_axis(:), sample_s | ||
| 724 | logical, intent(in) :: radial_is_rho | ||
| 725 | real(dp), intent(out) :: vz_axis(:) | ||
| 726 | integer, parameter :: nang = 8 | ||
| 727 | real(dp) :: values(axis_rhs_count, nang, 2) | ||
| 728 | real(dp) :: mean(axis_rhs_count, 2) | ||
| 729 | real(dp) :: m1c(axis_rhs_count, 2), m1s(axis_rhs_count, 2) | ||
| 730 | real(dp) :: m2c(axis_rhs_count, 2), m2s(axis_rhs_count, 2) | ||
| 731 | real(dp) :: coeff(axis_rhs_count), reconstructed(axis_rhs_count) | ||
| 732 | real(dp) :: raw(axis_rhs_count), zsample(axis_rhs_count) | ||
| 733 | real(dp) :: angle, radius, radius0, r2, blend, blend_arg | ||
| 734 | integer :: iangle, iring | ||
| 735 | |||
| 736 | 534744 | r2=z_axis(1)*z_axis(1)+z_axis(2)*z_axis(2) | |
| 737 |
2/2✓ Branch 0 taken 534722 times.
✓ Branch 1 taken 22 times.
|
534744 | if(r2.ge.sample_s) then |
| 738 | 534722 | call evaluate_axis_rhs_raw(tau, z_axis, vz_axis, radial_is_rho) | |
| 739 | 534722 | return | |
| 740 | endif | ||
| 741 | |||
| 742 | ! Reconstruct the physical guiding-centre vector field itself in | ||
| 743 | ! regular Cartesian modes. Unlike independent tensor fits, this | ||
| 744 | ! retains all algebraic identities among the magfie quantities. | ||
| 745 | 22 | radius0=sqrt(sample_s) | |
| 746 |
2/2✓ Branch 0 taken 66 times.
✓ Branch 1 taken 22 times.
|
88 | zsample(3:5)=z_axis(3:5) |
| 747 |
2/2✓ Branch 0 taken 44 times.
✓ Branch 1 taken 22 times.
|
66 | do iring=1,2 |
| 748 | 44 | radius=real(iring,dp)*radius0 | |
| 749 |
2/2✓ Branch 0 taken 352 times.
✓ Branch 1 taken 44 times.
|
418 | do iangle=1,nang |
| 750 | 352 | angle=2.d0*acos(-1.d0)*real(iangle-1,dp)/real(nang,dp) | |
| 751 | 352 | zsample(1)=radius*cos(angle) | |
| 752 | 352 | zsample(2)=radius*sin(angle) | |
| 753 | call evaluate_axis_rhs_raw(tau, zsample, & | ||
| 754 | 396 | values(:, iangle, iring), radial_is_rho) | |
| 755 | enddo | ||
| 756 | enddo | ||
| 757 | |||
| 758 | 22 | mean=0.d0 | |
| 759 | 22 | m1c=0.d0 | |
| 760 | 22 | m1s=0.d0 | |
| 761 | 22 | m2c=0.d0 | |
| 762 | 22 | m2s=0.d0 | |
| 763 |
2/2✓ Branch 0 taken 44 times.
✓ Branch 1 taken 22 times.
|
66 | do iring=1,2 |
| 764 |
2/2✓ Branch 0 taken 352 times.
✓ Branch 1 taken 44 times.
|
418 | do iangle=1,nang |
| 765 | 352 | angle=2.d0*acos(-1.d0)*real(iangle-1,dp)/real(nang,dp) | |
| 766 |
2/2✓ Branch 0 taken 1760 times.
✓ Branch 1 taken 352 times.
|
2112 | mean(:,iring)=mean(:,iring)+values(:,iangle,iring) |
| 767 |
2/2✓ Branch 0 taken 1760 times.
✓ Branch 1 taken 352 times.
|
2112 | m1c(:,iring)=m1c(:,iring)+values(:,iangle,iring)*cos(angle) |
| 768 |
2/2✓ Branch 0 taken 1760 times.
✓ Branch 1 taken 352 times.
|
2112 | m1s(:,iring)=m1s(:,iring)+values(:,iangle,iring)*sin(angle) |
| 769 |
2/2✓ Branch 0 taken 1760 times.
✓ Branch 1 taken 352 times.
|
2112 | m2c(:,iring)=m2c(:,iring)+values(:,iangle,iring)*cos(2.d0*angle) |
| 770 |
2/2✓ Branch 0 taken 1760 times.
✓ Branch 1 taken 352 times.
|
2156 | m2s(:,iring)=m2s(:,iring)+values(:,iangle,iring)*sin(2.d0*angle) |
| 771 | enddo | ||
| 772 | enddo | ||
| 773 |
4/4✓ Branch 0 taken 44 times.
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 220 times.
✓ Branch 3 taken 44 times.
|
286 | mean=mean/real(nang,dp) |
| 774 |
4/4✓ Branch 0 taken 44 times.
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 220 times.
✓ Branch 3 taken 44 times.
|
286 | m1c=2.d0*m1c/real(nang,dp) |
| 775 |
4/4✓ Branch 0 taken 44 times.
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 220 times.
✓ Branch 3 taken 44 times.
|
286 | m1s=2.d0*m1s/real(nang,dp) |
| 776 |
4/4✓ Branch 0 taken 44 times.
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 220 times.
✓ Branch 3 taken 44 times.
|
286 | m2c=2.d0*m2c/real(nang,dp) |
| 777 |
4/4✓ Branch 0 taken 44 times.
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 220 times.
✓ Branch 3 taken 44 times.
|
286 | m2s=2.d0*m2s/real(nang,dp) |
| 778 | |||
| 779 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | reconstructed=(4.d0*mean(:,1)-mean(:,2))/3.d0 |
| 780 | coeff=(4.d0*m1c(:,1)/radius0 & | ||
| 781 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | -m1c(:,2)/(2.d0*radius0))/3.d0 |
| 782 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | reconstructed=reconstructed+coeff*z_axis(1) |
| 783 | coeff=(4.d0*m1s(:,1)/radius0 & | ||
| 784 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | -m1s(:,2)/(2.d0*radius0))/3.d0 |
| 785 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | reconstructed=reconstructed+coeff*z_axis(2) |
| 786 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | coeff=(mean(:,2)-mean(:,1))/(3.d0*sample_s) |
| 787 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | reconstructed=reconstructed+coeff*r2 |
| 788 | coeff=(4.d0*m2c(:,1)/sample_s & | ||
| 789 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | -m2c(:,2)/(4.d0*sample_s))/3.d0 |
| 790 | reconstructed=reconstructed+coeff* & | ||
| 791 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | (z_axis(1)*z_axis(1)-z_axis(2)*z_axis(2)) |
| 792 | coeff=(4.d0*m2s(:,1)/sample_s & | ||
| 793 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | -m2s(:,2)/(4.d0*sample_s))/3.d0 |
| 794 | reconstructed=reconstructed+coeff* & | ||
| 795 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | (2.d0*z_axis(1)*z_axis(2)) |
| 796 | |||
| 797 | ! Join the reconstruction to the exact pullback with a C1 | ||
| 798 | ! smoothstep in the outer half of the disk. | ||
| 799 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 20 times.
|
22 | if(r2.gt.0.25d0*sample_s) then |
| 800 | 2 | call evaluate_axis_rhs_raw(tau, z_axis, raw, radial_is_rho) | |
| 801 | 2 | blend_arg=2.d0*sqrt(r2/sample_s)-1.d0 | |
| 802 | 2 | blend=blend_arg*blend_arg*(3.d0-2.d0*blend_arg) | |
| 803 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 10 times.
|
12 | reconstructed=(1.d0-blend)*reconstructed+blend*raw |
| 804 | endif | ||
| 805 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 22 times.
|
132 | vz_axis=reconstructed |
| 806 | end subroutine evaluate_axis_rhs | ||
| 807 | |||
| 808 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 535076 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 535076 times.
|
535076 | subroutine evaluate_axis_rhs_raw(tau, z_axis, vz_axis, radial_is_rho) |
| 809 | real(dp), intent(in) :: tau, z_axis(:) | ||
| 810 | real(dp), intent(out) :: vz_axis(:) | ||
| 811 | logical, intent(in) :: radial_is_rho | ||
| 812 | real(dp) :: zpolar(5), vpolar(5) | ||
| 813 | |||
| 814 |
3/4✓ Branch 0 taken 1070152 times.
✓ Branch 1 taken 535076 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 535076 times.
|
1605228 | if (.not. (sum(z_axis(1:2)**2) > 0.0_dp)) error stop & |
| 815 | ✗ | 'raw axis-RHS evaluation requires nonzero radius' | |
| 816 |
2/8✗ Branch 0 not taken.
✓ Branch 1 taken 535076 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 535076 times.
|
535076 | call axis_cartesian_to_polar(z_axis, zpolar, radial_is_rho) |
| 817 | 535076 | call velo_can(tau, zpolar, vpolar) | |
| 818 |
9/16✗ Branch 0 not taken.
✓ Branch 1 taken 535076 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 534720 times.
✓ Branch 7 taken 356 times.
✓ Branch 8 taken 534720 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 535076 times.
✓ Branch 13 taken 534720 times.
✓ Branch 14 taken 356 times.
✓ Branch 15 taken 2673600 times.
✓ Branch 16 taken 534720 times.
|
3743396 | call axis_polar_velocity(z_axis, vpolar, vz_axis, radial_is_rho) |
| 819 | 535076 | end subroutine evaluate_axis_rhs_raw | |
| 820 | ! | ||
| 821 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 822 | ! | ||
| 823 | 3759304 | subroutine integrate_axis_segment(z,tau1,tau2,relerr,near_axis, & | |
| 824 | axis_floor,step_limit) | ||
| 825 | use odeint_allroutines_sub, only : odeint_allroutines | ||
| 826 | implicit none | ||
| 827 | logical, intent(in) :: near_axis | ||
| 828 | integer, intent(in), optional :: step_limit | ||
| 829 | real(dp), intent(in) :: tau1,tau2,relerr,axis_floor | ||
| 830 | real(dp), intent(inout) :: z(5) | ||
| 831 | real(dp) :: context | ||
| 832 | |||
| 833 | 3759304 | context=axis_floor | |
| 834 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 3759301 times.
|
3759304 | if(present(step_limit)) then |
| 835 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if(near_axis) then |
| 836 | call odeint_allroutines(z,5,context,tau1,tau2,relerr, & | ||
| 837 | ✗ | velo_axis_context,step_limit=step_limit) | |
| 838 | else | ||
| 839 | call odeint_allroutines(z,5,context,tau1,tau2,relerr, & | ||
| 840 | 3 | velo_can_context,step_limit=step_limit) | |
| 841 | endif | ||
| 842 |
2/2✓ Branch 0 taken 89120 times.
✓ Branch 1 taken 3670181 times.
|
3759301 | else if(near_axis) then |
| 843 | call odeint_allroutines(z,5,context,tau1,tau2,relerr, & | ||
| 844 | 89120 | velo_axis_context) | |
| 845 | else | ||
| 846 | 3670181 | call odeint_allroutines(z,5,tau1,tau2,relerr,velo_can) | |
| 847 | endif | ||
| 848 | 3759304 | end subroutine integrate_axis_segment | |
| 849 | |||
| 850 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
|
18 | subroutine velo_can_context(tau,z,vz,context) |
| 851 | implicit none | ||
| 852 | class(*), intent(in) :: context | ||
| 853 | real(dp), intent(in) :: tau,z(:) | ||
| 854 | real(dp), intent(out) :: vz(:) | ||
| 855 | |||
| 856 | select type(context) | ||
| 857 | type is(integer) | ||
| 858 | end select | ||
| 859 | 18 | call velo_can(tau,z,vz) | |
| 860 | 18 | end subroutine velo_can_context | |
| 861 | |||
| 862 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 534720 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 534720 times.
|
534720 | subroutine velo_axis_context(tau,z,vz,context) |
| 863 | implicit none | ||
| 864 | class(*), intent(in) :: context | ||
| 865 | real(dp), intent(in) :: tau,z(:) | ||
| 866 | real(dp), intent(out) :: vz(:) | ||
| 867 | |||
| 868 | select type(context) | ||
| 869 | type is(real(dp)) | ||
| 870 | 534720 | call velo_axis_regularized(tau,z,vz,context) | |
| 871 | class default | ||
| 872 | ✗ | error stop 'invalid axis ODE context' | |
| 873 | end select | ||
| 874 | 534720 | end subroutine velo_axis_context | |
| 875 | |||
| 876 | 3759308 | subroutine orbit_timestep_axis(z, dtau, dtaumin, relerr, ierr, step_limit, & | |
| 877 | axis_floor, near_axis_limit) | ||
| 878 | use magfie_sub, only: axis_radial_coordinate | ||
| 879 | use finite_check, only: is_finite | ||
| 880 | use odeint_allroutines_sub, only: odeint_has_failed | ||
| 881 | use chamb_sub, only: chamb_can | ||
| 882 | ! | ||
| 883 | implicit none | ||
| 884 | ! | ||
| 885 | integer, parameter :: ndim = 5, nstepmax = 1000000 | ||
| 886 | logical :: near_axis, radial_is_rho | ||
| 887 | integer :: ierr | ||
| 888 | integer, intent(in), optional :: step_limit | ||
| 889 | real(dp), intent(in), optional :: axis_floor, near_axis_limit | ||
| 890 | real(dp) :: dtau, dtaumin, phi, tau1, tau2, axis_floor_local | ||
| 891 | real(dp) :: relerr, snear_axis, near_radius_squared | ||
| 892 | real(dp) :: z_mapped(5) | ||
| 893 | ! | ||
| 894 | real(dp), dimension(2) :: y | ||
| 895 | real(dp), dimension(ndim) :: z, z_initial | ||
| 896 | ! | ||
| 897 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 3759304 times.
|
3759308 | if (abs(dtaumin*nstepmax) .le. abs(dtau)) then |
| 898 | 4 | ierr = 2 | |
| 899 | 4 | print *, 'orbit_timestep: number of steps exceeds nstepmax' | |
| 900 | 4 | return | |
| 901 | end if | ||
| 902 | ! | ||
| 903 | 3759304 | ierr = 0 | |
| 904 | 3759304 | axis_floor_local = 1.d-8 | |
| 905 | 3759304 | snear_axis = 0.01d0 | |
| 906 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3759303 times.
|
3759304 | if (present(axis_floor)) axis_floor_local = axis_floor |
| 907 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3759303 times.
|
3759304 | if (present(near_axis_limit)) snear_axis = near_axis_limit |
| 908 | if (.not. is_finite(axis_floor_local) .or. & | ||
| 909 | axis_floor_local .le. 0.d0 .or. & | ||
| 910 | .not. is_finite(snear_axis) .or. & | ||
| 911 |
5/10✓ Branch 0 taken 3759304 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3759304 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3759304 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3759304 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 3759304 times.
|
3759304 | snear_axis .le. axis_floor_local .or. snear_axis .gt. 1.d0) then |
| 912 | ✗ | ierr = 2 | |
| 913 | ✗ | return | |
| 914 | end if | ||
| 915 | 3759304 | call axis_radial_coordinate(radial_is_rho, ierr) | |
| 916 |
1/2✓ Branch 0 taken 3759304 times.
✗ Branch 1 not taken.
|
3759304 | if (ierr /= 0) return |
| 917 | 3759304 | near_radius_squared = axis_native_radius_squared(snear_axis, radial_is_rho) | |
| 918 | 3759304 | z_initial = z | |
| 919 | 3759304 | y(1) = z(1) | |
| 920 | 3759304 | y(2) = z(2) | |
| 921 | 3759304 | phi = z(3) | |
| 922 | ! | ||
| 923 | 3759304 | call chamb_can(y, phi, ierr) | |
| 924 | ! | ||
| 925 |
1/2✓ Branch 0 taken 3759304 times.
✗ Branch 1 not taken.
|
3759304 | if (ierr .eq. 1) return |
| 926 | 3759304 | tau1 = 0.d0 | |
| 927 | 3759304 | tau2 = dtaumin | |
| 928 | ! | ||
| 929 | 3759304 | near_axis = .false. | |
| 930 | ! | ||
| 931 |
3/8✓ Branch 0 taken 3759304 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3759304 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3759304 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
3759304 | do while ((dtau > 0 .and. (tau2 .lt. dtau)) .or. (dtau < 0 .and. (tau2 .gt. dtau))) |
| 932 | ✗ | if (near_axis) then | |
| 933 | ✗ | if (z(1)**2 + z(2)**2 .gt. near_radius_squared) then | |
| 934 | ✗ | near_axis = .false. | |
| 935 | ✗ | call axis_cartesian_to_polar(z, z_mapped, radial_is_rho) | |
| 936 | ✗ | z = z_mapped | |
| 937 | ! | ||
| 938 | call integrate_axis_segment(z, tau1, tau2, relerr, .false., & | ||
| 939 | ✗ | axis_floor_local, step_limit) | |
| 940 | ✗ | if (odeint_has_failed() .or. .not. all(is_finite(z))) then | |
| 941 | ✗ | z = z_initial | |
| 942 | ✗ | ierr = 2 | |
| 943 | ✗ | return | |
| 944 | end if | ||
| 945 | ! | ||
| 946 | ✗ | y(1) = z(1) | |
| 947 | ✗ | y(2) = z(2) | |
| 948 | ✗ | phi = z(3) | |
| 949 | ! | ||
| 950 | ✗ | call chamb_can(y, phi, ierr) | |
| 951 | ! | ||
| 952 | ✗ | if (ierr .eq. 1) return | |
| 953 | else | ||
| 954 | ! | ||
| 955 | call integrate_axis_segment(z, tau1, tau2, relerr, .true., & | ||
| 956 | ✗ | axis_floor_local, step_limit) | |
| 957 | ✗ | if (odeint_has_failed() .or. .not. all(is_finite(z))) then | |
| 958 | ✗ | z = z_initial | |
| 959 | ✗ | ierr = 2 | |
| 960 | ✗ | return | |
| 961 | end if | ||
| 962 | ! | ||
| 963 | end if | ||
| 964 | else | ||
| 965 | ✗ | if (z(1) .lt. snear_axis) then | |
| 966 | ✗ | near_axis = .true. | |
| 967 | ✗ | call axis_polar_to_cartesian(z, z_mapped, radial_is_rho) | |
| 968 | ✗ | z = z_mapped | |
| 969 | ! | ||
| 970 | call integrate_axis_segment(z, tau1, tau2, relerr, .true., & | ||
| 971 | ✗ | axis_floor_local, step_limit) | |
| 972 | ✗ | if (odeint_has_failed() .or. .not. all(is_finite(z))) then | |
| 973 | ✗ | z = z_initial | |
| 974 | ✗ | ierr = 2 | |
| 975 | ✗ | return | |
| 976 | end if | ||
| 977 | ! | ||
| 978 | else | ||
| 979 | ! | ||
| 980 | call integrate_axis_segment(z, tau1, tau2, relerr, .false., & | ||
| 981 | ✗ | axis_floor_local, step_limit) | |
| 982 | ✗ | if (odeint_has_failed() .or. .not. all(is_finite(z))) then | |
| 983 | ✗ | z = z_initial | |
| 984 | ✗ | ierr = 2 | |
| 985 | ✗ | return | |
| 986 | end if | ||
| 987 | ! | ||
| 988 | ✗ | y(1) = z(1) | |
| 989 | ✗ | y(2) = z(2) | |
| 990 | ✗ | phi = z(3) | |
| 991 | ! | ||
| 992 | ✗ | call chamb_can(y, phi, ierr) | |
| 993 | ! | ||
| 994 | ✗ | if (ierr .eq. 1) return | |
| 995 | ! | ||
| 996 | end if | ||
| 997 | end if | ||
| 998 | ✗ | tau1 = tau2 | |
| 999 | ✗ | tau2 = tau2 + dtaumin | |
| 1000 | end do | ||
| 1001 | ! | ||
| 1002 | 3759304 | tau2 = dtau | |
| 1003 | ! | ||
| 1004 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3759304 times.
|
3759304 | if (near_axis) then |
| 1005 | ✗ | if (z(1)**2 + z(2)**2 .gt. near_radius_squared) then | |
| 1006 | ✗ | near_axis = .false. | |
| 1007 | ✗ | call axis_cartesian_to_polar(z, z_mapped, radial_is_rho) | |
| 1008 | ✗ | z = z_mapped | |
| 1009 | ! | ||
| 1010 | call integrate_axis_segment(z, tau1, tau2, relerr, .false., & | ||
| 1011 | ✗ | axis_floor_local, step_limit) | |
| 1012 | ✗ | if (odeint_has_failed() .or. .not. all(is_finite(z))) then | |
| 1013 | ✗ | z = z_initial | |
| 1014 | ✗ | ierr = 2 | |
| 1015 | ✗ | return | |
| 1016 | end if | ||
| 1017 | ! | ||
| 1018 | ✗ | y(1) = z(1) | |
| 1019 | ✗ | y(2) = z(2) | |
| 1020 | ✗ | phi = z(3) | |
| 1021 | ! | ||
| 1022 | ✗ | call chamb_can(y, phi, ierr) | |
| 1023 | ! | ||
| 1024 | ✗ | if (ierr .eq. 1) return | |
| 1025 | else | ||
| 1026 | ! | ||
| 1027 | call integrate_axis_segment(z, tau1, tau2, relerr, .true., & | ||
| 1028 | ✗ | axis_floor_local, step_limit) | |
| 1029 | ✗ | if (odeint_has_failed() .or. .not. all(is_finite(z))) then | |
| 1030 | ✗ | z = z_initial | |
| 1031 | ✗ | ierr = 2 | |
| 1032 | ✗ | return | |
| 1033 | end if | ||
| 1034 | ! | ||
| 1035 | end if | ||
| 1036 | else | ||
| 1037 |
2/2✓ Branch 0 taken 89120 times.
✓ Branch 1 taken 3670184 times.
|
3759304 | if (z(1) .lt. snear_axis) then |
| 1038 | 89120 | near_axis = .true. | |
| 1039 | 89120 | call axis_polar_to_cartesian(z, z_mapped, radial_is_rho) | |
| 1040 | 89120 | z = z_mapped | |
| 1041 | ! | ||
| 1042 | call integrate_axis_segment(z, tau1, tau2, relerr, .true., & | ||
| 1043 | 89120 | axis_floor_local, step_limit) | |
| 1044 |
5/8✓ Branch 0 taken 445600 times.
✓ Branch 1 taken 89120 times.
✓ Branch 2 taken 445600 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 89120 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 89120 times.
|
534720 | if (odeint_has_failed() .or. .not. all(is_finite(z))) then |
| 1045 | ✗ | z = z_initial | |
| 1046 | ✗ | ierr = 2 | |
| 1047 | ✗ | return | |
| 1048 | end if | ||
| 1049 | ! | ||
| 1050 | else | ||
| 1051 | ! | ||
| 1052 | call integrate_axis_segment(z, tau1, tau2, relerr, .false., & | ||
| 1053 | 3670184 | axis_floor_local, step_limit) | |
| 1054 |
5/8✓ Branch 0 taken 18350920 times.
✓ Branch 1 taken 3670184 times.
✓ Branch 2 taken 18350920 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3670184 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 3670184 times.
|
22021104 | if (odeint_has_failed() .or. .not. all(is_finite(z))) then |
| 1055 | ✗ | z = z_initial | |
| 1056 | ✗ | ierr = 2 | |
| 1057 | ✗ | return | |
| 1058 | end if | ||
| 1059 | ! | ||
| 1060 | 3670184 | y(1) = z(1) | |
| 1061 | 3670184 | y(2) = z(2) | |
| 1062 | 3670184 | phi = z(3) | |
| 1063 | ! | ||
| 1064 | 3670184 | call chamb_can(y, phi, ierr) | |
| 1065 | ! | ||
| 1066 |
2/2✓ Branch 0 taken 3670180 times.
✓ Branch 1 taken 4 times.
|
3670184 | if (ierr .eq. 1) return |
| 1067 | ! | ||
| 1068 | end if | ||
| 1069 | end if | ||
| 1070 | ! | ||
| 1071 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3670180 times.
|
3670180 | if (near_axis) then |
| 1072 | 89120 | call axis_cartesian_to_polar(z, z_mapped, radial_is_rho) | |
| 1073 | 89120 | z = z_mapped | |
| 1074 | end if | ||
| 1075 | ! | ||
| 1076 | 3759308 | end subroutine orbit_timestep_axis | |
| 1077 | ! | ||
| 1078 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 1079 | ✗ | end module alpha_lifetime_sub | |
| 1080 |