| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module exchange_get_cancoord_mod | ||
| 2 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 3 | |||
| 4 | implicit none | ||
| 5 | private | ||
| 6 | |||
| 7 | logical, public :: onlytheta | ||
| 8 | real(dp), public :: vartheta_c, varphi_c, sqg, aiota, Bcovar_vartheta, & | ||
| 9 | Bcovar_varphi, A_theta, A_phi, theta, Bctrvr_vartheta, & | ||
| 10 | Bctrvr_varphi | ||
| 11 | |||
| 12 | !$omp threadprivate(onlytheta, vartheta_c, varphi_c, sqg, aiota) | ||
| 13 | !$omp threadprivate(Bcovar_vartheta, Bcovar_varphi, A_theta, A_phi) | ||
| 14 | !$omp threadprivate(theta, Bctrvr_vartheta, Bctrvr_varphi) | ||
| 15 | |||
| 16 | end module exchange_get_cancoord_mod | ||
| 17 | |||
| 18 | module get_can_sub | ||
| 19 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 20 | use spl_three_to_five_sub | ||
| 21 | use stencil_utils | ||
| 22 | use field, only: magnetic_field_t, vmec_field_t, create_vmec_field, field_clone | ||
| 23 | use field_newton, only: newton_theta_from_canonical | ||
| 24 | use interpolate, only: BatchSplineData1D, BatchSplineData3D, & | ||
| 25 | construct_batch_splines_1d, construct_batch_splines_3d, & | ||
| 26 | evaluate_batch_splines_1d_der2, & | ||
| 27 | evaluate_batch_splines_3d_der, & | ||
| 28 | evaluate_batch_splines_3d_der2, & | ||
| 29 | evaluate_batch_splines_3d_der2_rmix, & | ||
| 30 | destroy_batch_splines_1d, destroy_batch_splines_3d | ||
| 31 | |||
| 32 | implicit none | ||
| 33 | private | ||
| 34 | |||
| 35 | public :: get_canonical_coordinates, get_canonical_coordinates_with_field | ||
| 36 | public :: splint_can_coord | ||
| 37 | public :: can_to_vmec, vmec_to_can, vmec_to_cyl | ||
| 38 | public :: deallocate_can_coord | ||
| 39 | public :: reset_canflux_batch_splines | ||
| 40 | |||
| 41 | ! Constants | ||
| 42 | real(dp), parameter :: TWOPI = 2.0_dp*3.14159265358979_dp | ||
| 43 | |||
| 44 | ! Module variable to store the field for use in subroutines | ||
| 45 | class(magnetic_field_t), allocatable :: current_field | ||
| 46 | !$omp threadprivate(current_field) | ||
| 47 | |||
| 48 | ! Batch spline for A_phi (vector potential) | ||
| 49 | type(BatchSplineData1D), save :: aphi_batch_spline | ||
| 50 | logical, save :: aphi_batch_spline_ready = .false. | ||
| 51 | |||
| 52 | ! Batch spline for G_c (generating function) | ||
| 53 | type(BatchSplineData3D), save :: G_batch_spline | ||
| 54 | logical, save :: G_batch_spline_ready = .false. | ||
| 55 | |||
| 56 | ! Batch splines for sqg_c, B_vartheta_c, B_varphi_c (separate NQ=1 splines) | ||
| 57 | type(BatchSplineData3D), save :: sqg_batch_spline | ||
| 58 | type(BatchSplineData3D), save :: Bt_batch_spline | ||
| 59 | type(BatchSplineData3D), save :: Bp_batch_spline | ||
| 60 | logical, save :: sqg_batch_spline_ready = .false. | ||
| 61 | logical, save :: Bt_batch_spline_ready = .false. | ||
| 62 | logical, save :: Bp_batch_spline_ready = .false. | ||
| 63 | |||
| 64 | contains | ||
| 65 | |||
| 66 | 12 | subroutine get_canonical_coordinates_with_field(field) | |
| 67 | implicit none | ||
| 68 | |||
| 69 | class(magnetic_field_t), intent(in) :: field | ||
| 70 | |||
| 71 | ! Store field in module variable for use in nested subroutines | ||
| 72 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
6 | call field_clone(field, current_field) |
| 73 | |||
| 74 | 6 | call reset_canflux_batch_splines | |
| 75 | |||
| 76 | ! Call the actual implementation | ||
| 77 | 6 | call get_canonical_coordinates_impl | |
| 78 | |||
| 79 | 6 | end subroutine get_canonical_coordinates_with_field | |
| 80 | |||
| 81 | 1 | subroutine get_canonical_coordinates | |
| 82 | ! Backward compatibility wrapper - uses VMEC field by default | ||
| 83 | 1 | type(vmec_field_t) :: vmec_field | |
| 84 | |||
| 85 | 1 | call create_vmec_field(vmec_field) | |
| 86 | 1 | call get_canonical_coordinates_with_field(vmec_field) | |
| 87 |
3/6✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
✓ Branch 5 taken 1 time.
✗ Branch 6 not taken.
|
1 | end subroutine get_canonical_coordinates |
| 88 | |||
| 89 | 6 | subroutine reset_canflux_batch_splines | |
| 90 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (aphi_batch_spline_ready) then |
| 91 | ✗ | call destroy_batch_splines_1d(aphi_batch_spline) | |
| 92 | ✗ | aphi_batch_spline_ready = .false. | |
| 93 | end if | ||
| 94 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (G_batch_spline_ready) then |
| 95 | ✗ | call destroy_batch_splines_3d(G_batch_spline) | |
| 96 | ✗ | G_batch_spline_ready = .false. | |
| 97 | end if | ||
| 98 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (sqg_batch_spline_ready) then |
| 99 | ✗ | call destroy_batch_splines_3d(sqg_batch_spline) | |
| 100 | ✗ | sqg_batch_spline_ready = .false. | |
| 101 | end if | ||
| 102 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (Bt_batch_spline_ready) then |
| 103 | ✗ | call destroy_batch_splines_3d(Bt_batch_spline) | |
| 104 | ✗ | Bt_batch_spline_ready = .false. | |
| 105 | end if | ||
| 106 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (Bp_batch_spline_ready) then |
| 107 | ✗ | call destroy_batch_splines_3d(Bp_batch_spline) | |
| 108 | ✗ | Bp_batch_spline_ready = .false. | |
| 109 | end if | ||
| 110 | 6 | end subroutine reset_canflux_batch_splines | |
| 111 | |||
| 112 | 6 | subroutine get_canonical_coordinates_impl | |
| 113 | use canonical_coordinates_mod, only: ns_c, n_theta_c, n_phi_c, & | ||
| 114 | hs_c, h_theta_c, h_phi_c, & | ||
| 115 | ns_s_c, ns_tp_c, & | ||
| 116 | nh_stencil, G_c, sqg_c, & | ||
| 117 | B_vartheta_c, B_varphi_c | ||
| 118 | use vector_potentail_mod, only: ns, hs | ||
| 119 | use exchange_get_cancoord_mod, only: vartheta_c, varphi_c, sqg, aiota, & | ||
| 120 | Bcovar_vartheta, Bcovar_varphi, & | ||
| 121 | onlytheta | ||
| 122 | use new_vmec_stuff_mod, only: n_theta, n_phi, h_theta, h_phi, ns_s, ns_tp | ||
| 123 | use odeint_allroutines_sub, only: odeint_allroutines | ||
| 124 | |||
| 125 | implicit none | ||
| 126 | |||
| 127 | real(dp), parameter :: relerr = 1.0e-10_dp | ||
| 128 | integer :: i_theta, i_phi, i_sten, ndim, is_beg | ||
| 129 | 6 | integer, dimension(:), allocatable :: ipoi_t, ipoi_p | |
| 130 | 6 | real(dp), dimension(:), allocatable :: y, dy | |
| 131 | real(dp) :: dstencil_theta(-nh_stencil:nh_stencil), & | ||
| 132 | dstencil_phi(-nh_stencil:nh_stencil) | ||
| 133 | |||
| 134 | real(dp) :: r, r1, r2, G_beg, dG_c_dt, dG_c_dp | ||
| 135 | integer :: is | ||
| 136 | integer :: i_ctr | ||
| 137 | |||
| 138 | 6 | ns_c = ns | |
| 139 | 6 | n_theta_c = n_theta | |
| 140 | 6 | n_phi_c = n_phi | |
| 141 | 6 | h_theta_c = h_theta | |
| 142 | 6 | h_phi_c = h_phi | |
| 143 | 6 | hs_c = hs | |
| 144 | |||
| 145 | ! Initialize derivative stencils using stencil_utils module | ||
| 146 | 6 | call init_derivative_stencil(nh_stencil, h_theta_c, dstencil_theta) | |
| 147 | 6 | call init_derivative_stencil(nh_stencil, h_phi_c, dstencil_phi) | |
| 148 | |||
| 149 |
7/14✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 6 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 6 times.
|
18 | allocate (ipoi_t(1 - nh_stencil:n_theta_c + nh_stencil)) |
| 150 |
7/14✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 6 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 6 times.
|
18 | allocate (ipoi_p(1 - nh_stencil:n_phi_c + nh_stencil)) |
| 151 | |||
| 152 |
2/2✓ Branch 0 taken 160 times.
✓ Branch 1 taken 6 times.
|
166 | do i_theta = 1, n_theta_c |
| 153 | 166 | ipoi_t(i_theta) = i_theta | |
| 154 | end do | ||
| 155 | |||
| 156 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 358 times.
|
364 | do i_phi = 1, n_phi_c |
| 157 | 364 | ipoi_p(i_phi) = i_phi | |
| 158 | end do | ||
| 159 | |||
| 160 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 6 times.
|
24 | do i_sten = 1, nh_stencil |
| 161 | 18 | ipoi_t(1 - i_sten) = ipoi_t(n_theta - i_sten) | |
| 162 | 18 | ipoi_t(n_theta_c + i_sten) = ipoi_t(1 + i_sten) | |
| 163 | 18 | ipoi_p(1 - i_sten) = ipoi_p(n_phi_c - i_sten) | |
| 164 | 24 | ipoi_p(n_phi_c + i_sten) = ipoi_p(1 + i_sten) | |
| 165 | end do | ||
| 166 | |||
| 167 |
11/22✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 6 times.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 6 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 6 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 6 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 6 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 6 times.
|
30 | allocate (G_c(ns_c, n_theta_c, n_phi_c)) |
| 168 |
9/18✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 6 times.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 6 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 6 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 6 times.
|
24 | allocate (sqg_c(ns_c, n_theta_c, n_phi_c)) |
| 169 |
9/18✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 6 times.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 6 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 6 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 6 times.
|
24 | allocate (B_vartheta_c(ns_c, n_theta_c, n_phi_c)) |
| 170 |
9/18✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 6 times.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 6 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 6 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 6 times.
|
24 | allocate (B_varphi_c(ns_c, n_theta_c, n_phi_c)) |
| 171 | |||
| 172 | 6 | onlytheta = .false. | |
| 173 | 6 | ndim = 1 | |
| 174 | 6 | is_beg = 1 | |
| 175 | 6 | G_beg = 1.0e-8_dp | |
| 176 | |||
| 177 | 6 | i_ctr = 0 | |
| 178 | 6 | !$omp parallel private(y, dy, i_theta, i_phi, is, r1, r2, r, dG_c_dt, dG_c_dp) | |
| 179 | !$omp critical | ||
| 180 | allocate (y(ndim), dy(ndim)) | ||
| 181 | !$omp end critical | ||
| 182 | |||
| 183 | !$omp do | ||
| 184 | do i_theta = 1, n_theta_c | ||
| 185 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 186 | !$omp critical | ||
| 187 | i_ctr = i_ctr + 1 | ||
| 188 | call print_progress('integrate ODE: ', i_ctr, n_theta_c) | ||
| 189 | !$omp end critical | ||
| 190 | #endif | ||
| 191 | vartheta_c = h_theta_c*real(i_theta - 1, dp) | ||
| 192 | do i_phi = 1, n_phi_c | ||
| 193 | varphi_c = h_phi_c*real(i_phi - 1, dp) | ||
| 194 | |||
| 195 | G_c(is_beg, i_theta, i_phi) = G_beg | ||
| 196 | y(1) = G_beg | ||
| 197 | |||
| 198 | do is = is_beg - 1, 2, -1 | ||
| 199 | r1 = hs_c*real(is, dp) | ||
| 200 | r2 = hs_c*real(is - 1, dp) | ||
| 201 | |||
| 202 | call odeint_allroutines(y, ndim, r1, r2, relerr, rhs_cancoord) | ||
| 203 | |||
| 204 | G_c(is, i_theta, i_phi) = y(1) | ||
| 205 | end do | ||
| 206 | |||
| 207 | y(1) = G_beg | ||
| 208 | |||
| 209 | do is = is_beg + 1, ns_c | ||
| 210 | r1 = hs_c*real(is - 2, dp) | ||
| 211 | r2 = hs_c*real(is - 1, dp) | ||
| 212 | if (is == 2) r1 = 1.0e-8_dp | ||
| 213 | |||
| 214 | call odeint_allroutines(y, ndim, r1, r2, relerr, rhs_cancoord) | ||
| 215 | |||
| 216 | G_c(is, i_theta, i_phi) = y(1) | ||
| 217 | end do | ||
| 218 | end do | ||
| 219 | end do | ||
| 220 | !$omp end do | ||
| 221 | |||
| 222 | i_ctr = 0 | ||
| 223 | !$omp barrier | ||
| 224 | !$omp do | ||
| 225 | do i_theta = 1, n_theta_c | ||
| 226 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 227 | !$omp critical | ||
| 228 | i_ctr = i_ctr + 1 | ||
| 229 | call print_progress('compute components: ', i_ctr, n_theta_c) | ||
| 230 | !$omp end critical | ||
| 231 | #endif | ||
| 232 | vartheta_c = h_theta_c*real(i_theta - 1, dp) | ||
| 233 | do i_phi = 1, n_phi_c | ||
| 234 | varphi_c = h_phi_c*real(i_phi - 1, dp) | ||
| 235 | |||
| 236 | do is = 2, ns_c | ||
| 237 | r = hs_c*real(is - 1, dp) | ||
| 238 | y(1) = G_c(is, i_theta, i_phi) | ||
| 239 | |||
| 240 | call rhs_cancoord(r, y, dy) | ||
| 241 | |||
| 242 | dG_c_dt = sum(dstencil_theta*G_c(is, & | ||
| 243 | ipoi_t(i_theta - & | ||
| 244 | nh_stencil:i_theta + & | ||
| 245 | nh_stencil), i_phi)) | ||
| 246 | dG_c_dp = sum(dstencil_phi*G_c(is, i_theta, & | ||
| 247 | ipoi_p(i_phi - nh_stencil:i_phi + & | ||
| 248 | nh_stencil))) | ||
| 249 | sqg_c(is, i_theta, i_phi) = sqg*(1.0_dp + aiota*dG_c_dt + dG_c_dp) | ||
| 250 | B_vartheta_c(is, i_theta, i_phi) = Bcovar_vartheta + & | ||
| 251 | (aiota*Bcovar_vartheta + & | ||
| 252 | Bcovar_varphi)*dG_c_dt | ||
| 253 | B_varphi_c(is, i_theta, i_phi) = Bcovar_varphi + & | ||
| 254 | (aiota*Bcovar_vartheta + & | ||
| 255 | Bcovar_varphi)*dG_c_dp | ||
| 256 | end do | ||
| 257 | ! Extrapolate on-axis point (is=1) with parabola | ||
| 258 | sqg_c(1, i_theta, i_phi) = 3.0_dp*(sqg_c(2, i_theta, i_phi) & | ||
| 259 | - sqg_c(3, i_theta, i_phi)) + & | ||
| 260 | sqg_c(4, i_theta, i_phi) | ||
| 261 | B_vartheta_c(1, i_theta, i_phi) = 0.0_dp | ||
| 262 | B_varphi_c(1, i_theta, i_phi) = 3.0_dp*(B_varphi_c(2, i_theta, i_phi) & | ||
| 263 | - B_varphi_c(3, i_theta, & | ||
| 264 | i_phi)) + B_varphi_c(4, & | ||
| 265 | i_theta, i_phi) | ||
| 266 | end do | ||
| 267 | end do | ||
| 268 | !$omp end do | ||
| 269 | |||
| 270 | !$omp critical | ||
| 271 | deallocate (y, dy) | ||
| 272 | !$omp end critical | ||
| 273 | !$omp end parallel | ||
| 274 | |||
| 275 | 6 | ns_s_c = ns_s | |
| 276 | 6 | ns_tp_c = ns_tp | |
| 277 | |||
| 278 | 6 | onlytheta = .true. | |
| 279 | |||
| 280 | ! Build batch splines from computed grids | ||
| 281 | 6 | call build_canflux_aphi_batch_spline | |
| 282 | 6 | call build_canflux_G_batch_spline | |
| 283 | 6 | call build_canflux_sqg_Bt_Bp_batch_spline | |
| 284 | |||
| 285 |
6/12✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 6 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 6 times.
|
6 | deallocate (ipoi_t, ipoi_p, sqg_c, B_vartheta_c, B_varphi_c, G_c) |
| 286 | |||
| 287 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 6 times.
|
6 | end subroutine get_canonical_coordinates_impl |
| 288 | |||
| 289 | 6 | subroutine build_canflux_aphi_batch_spline | |
| 290 | use vector_potentail_mod, only: ns, hs, sA_phi | ||
| 291 | use new_vmec_stuff_mod, only: ns_A | ||
| 292 | |||
| 293 | integer :: order | ||
| 294 | |||
| 295 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (aphi_batch_spline_ready) then |
| 296 | ✗ | call destroy_batch_splines_1d(aphi_batch_spline) | |
| 297 | ✗ | aphi_batch_spline_ready = .false. | |
| 298 | end if | ||
| 299 | |||
| 300 | 6 | order = ns_A | |
| 301 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (order < 3 .or. order > 5) then |
| 302 | ✗ | error stop "build_canflux_aphi_batch_spline: spline order must be 3..5" | |
| 303 | end if | ||
| 304 | |||
| 305 | 6 | aphi_batch_spline%order = order | |
| 306 | 6 | aphi_batch_spline%num_points = ns | |
| 307 | 6 | aphi_batch_spline%periodic = .false. | |
| 308 | 6 | aphi_batch_spline%x_min = 0.0_dp | |
| 309 | 6 | aphi_batch_spline%h_step = hs | |
| 310 | 6 | aphi_batch_spline%num_quantities = 1 | |
| 311 | |||
| 312 |
9/18✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 6 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 6 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 6 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 6 times.
|
24 | allocate (aphi_batch_spline%coeff(1, 0:order, ns)) |
| 313 |
4/4✓ Branch 0 taken 300 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 1800 times.
✓ Branch 3 taken 300 times.
|
2106 | aphi_batch_spline%coeff(1, 0:order, :) = sA_phi(1:order + 1, :) |
| 314 | |||
| 315 | 6 | aphi_batch_spline_ready = .true. | |
| 316 | 6 | end subroutine build_canflux_aphi_batch_spline | |
| 317 | |||
| 318 | 6 | subroutine build_canflux_G_batch_spline | |
| 319 | use canonical_coordinates_mod, only: ns_c, n_theta_c, n_phi_c, & | ||
| 320 | hs_c, h_theta_c, h_phi_c, & | ||
| 321 | ns_s_c, ns_tp_c, G_c | ||
| 322 | |||
| 323 | integer :: order(3) | ||
| 324 | real(dp) :: x_min(3), x_max(3) | ||
| 325 | logical :: periodic(3) | ||
| 326 | 6 | real(dp), allocatable :: y_batch(:, :, :, :) | |
| 327 | |||
| 328 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (G_batch_spline_ready) then |
| 329 | ✗ | call destroy_batch_splines_3d(G_batch_spline) | |
| 330 | ✗ | G_batch_spline_ready = .false. | |
| 331 | end if | ||
| 332 | |||
| 333 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 6 times.
|
24 | order = [ns_s_c, ns_tp_c, ns_tp_c] |
| 334 |
5/6✓ Branch 0 taken 6 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 18 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
|
48 | if (minval(order) < 3 .or. maxval(order) > 5) then |
| 335 | ✗ | error stop "build_canflux_G_batch_spline: spline order must be 3..5" | |
| 336 | end if | ||
| 337 | |||
| 338 | 6 | x_min = [0.0_dp, 0.0_dp, 0.0_dp] | |
| 339 | 6 | x_max(1) = hs_c*real(ns_c - 1, dp) | |
| 340 | 6 | x_max(2) = h_theta_c*real(n_theta_c - 1, dp) | |
| 341 | 6 | x_max(3) = h_phi_c*real(n_phi_c - 1, dp) | |
| 342 | |||
| 343 | 6 | periodic = [.false., .true., .true.] | |
| 344 | |||
| 345 |
10/20✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 6 times.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 6 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 6 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 6 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 6 times.
|
30 | allocate (y_batch(ns_c, n_theta_c, n_phi_c, 1)) |
| 346 |
6/6✓ Branch 0 taken 358 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 10144 times.
✓ Branch 3 taken 358 times.
✓ Branch 4 taken 507200 times.
✓ Branch 5 taken 10144 times.
|
517708 | y_batch(:, :, :, 1) = G_c |
| 347 | |||
| 348 | call construct_batch_splines_3d(x_min, x_max, y_batch, order, periodic, & | ||
| 349 | 6 | G_batch_spline) | |
| 350 | 6 | G_batch_spline_ready = .true. | |
| 351 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | deallocate (y_batch) |
| 352 | 6 | end subroutine build_canflux_G_batch_spline | |
| 353 | |||
| 354 | 6 | subroutine build_canflux_sqg_Bt_Bp_batch_spline | |
| 355 | use canonical_coordinates_mod, only: ns_c, n_theta_c, n_phi_c, & | ||
| 356 | hs_c, h_theta_c, h_phi_c, & | ||
| 357 | ns_s_c, ns_tp_c, & | ||
| 358 | sqg_c, B_vartheta_c, B_varphi_c | ||
| 359 | |||
| 360 | integer :: order(3) | ||
| 361 | real(dp) :: x_min(3), x_max(3) | ||
| 362 | logical :: periodic(3) | ||
| 363 | 6 | real(dp), allocatable :: y_batch(:, :, :, :) | |
| 364 | |||
| 365 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (sqg_batch_spline_ready) then |
| 366 | ✗ | call destroy_batch_splines_3d(sqg_batch_spline) | |
| 367 | ✗ | sqg_batch_spline_ready = .false. | |
| 368 | end if | ||
| 369 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (Bt_batch_spline_ready) then |
| 370 | ✗ | call destroy_batch_splines_3d(Bt_batch_spline) | |
| 371 | ✗ | Bt_batch_spline_ready = .false. | |
| 372 | end if | ||
| 373 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (Bp_batch_spline_ready) then |
| 374 | ✗ | call destroy_batch_splines_3d(Bp_batch_spline) | |
| 375 | ✗ | Bp_batch_spline_ready = .false. | |
| 376 | end if | ||
| 377 | |||
| 378 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 6 times.
|
24 | order = [ns_s_c, ns_tp_c, ns_tp_c] |
| 379 |
5/6✓ Branch 0 taken 6 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 18 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
|
48 | if (minval(order) < 3 .or. maxval(order) > 5) then |
| 380 | ✗ | error stop "build_canflux_sqg_Bt_Bp_batch_spline: spline order must be 3..5" | |
| 381 | end if | ||
| 382 | |||
| 383 | 6 | x_min = [0.0_dp, 0.0_dp, 0.0_dp] | |
| 384 | 6 | x_max(1) = hs_c*real(ns_c - 1, dp) | |
| 385 | 6 | x_max(2) = h_theta_c*real(n_theta_c - 1, dp) | |
| 386 | 6 | x_max(3) = h_phi_c*real(n_phi_c - 1, dp) | |
| 387 | |||
| 388 | 6 | periodic = [.false., .true., .true.] | |
| 389 | |||
| 390 |
10/20✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 6 times.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 6 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 6 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 6 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 6 times.
|
30 | allocate (y_batch(ns_c, n_theta_c, n_phi_c, 1)) |
| 391 | |||
| 392 |
6/6✓ Branch 0 taken 358 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 10144 times.
✓ Branch 3 taken 358 times.
✓ Branch 4 taken 507200 times.
✓ Branch 5 taken 10144 times.
|
517708 | y_batch(:, :, :, 1) = sqg_c |
| 393 | call construct_batch_splines_3d(x_min, x_max, y_batch, order, periodic, & | ||
| 394 | 6 | sqg_batch_spline) | |
| 395 | 6 | sqg_batch_spline_ready = .true. | |
| 396 | |||
| 397 |
6/6✓ Branch 0 taken 358 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 10144 times.
✓ Branch 3 taken 358 times.
✓ Branch 4 taken 507200 times.
✓ Branch 5 taken 10144 times.
|
517708 | y_batch(:, :, :, 1) = B_vartheta_c |
| 398 | call construct_batch_splines_3d(x_min, x_max, y_batch, order, periodic, & | ||
| 399 | 6 | Bt_batch_spline) | |
| 400 | 6 | Bt_batch_spline_ready = .true. | |
| 401 | |||
| 402 |
6/6✓ Branch 0 taken 358 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 10144 times.
✓ Branch 3 taken 358 times.
✓ Branch 4 taken 507200 times.
✓ Branch 5 taken 10144 times.
|
517708 | y_batch(:, :, :, 1) = B_varphi_c |
| 403 | call construct_batch_splines_3d(x_min, x_max, y_batch, order, periodic, & | ||
| 404 | 6 | Bp_batch_spline) | |
| 405 | 6 | Bp_batch_spline_ready = .true. | |
| 406 | |||
| 407 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | deallocate (y_batch) |
| 408 | 6 | end subroutine build_canflux_sqg_Bt_Bp_batch_spline | |
| 409 | |||
| 410 | 18807970 | subroutine rhs_cancoord(r, y, dy) | |
| 411 | use exchange_get_cancoord_mod, only: vartheta_c, varphi_c, sqg, aiota, & | ||
| 412 | Bcovar_vartheta, Bcovar_varphi, & | ||
| 413 | theta, onlytheta | ||
| 414 | use spline_vmec_sub | ||
| 415 | use vmec_field_eval | ||
| 416 | |||
| 417 | implicit none | ||
| 418 | |||
| 419 | real(dp), intent(in) :: r | ||
| 420 | real(dp), intent(in) :: y(:) | ||
| 421 | real(dp), intent(out) :: dy(:) | ||
| 422 | |||
| 423 | real(dp), parameter :: epserr = 1.0e-14_dp | ||
| 424 | integer :: iter | ||
| 425 | real(dp) :: s, varphi, A_theta, A_phi, dA_theta_ds, dA_phi_ds, & | ||
| 426 | alam, dl_ds, dl_dt, dl_dp, Bctrvr_vartheta, Bctrvr_varphi, Bcovar_r | ||
| 427 | logical :: converged | ||
| 428 | |||
| 429 | real(dp) :: vartheta, daiota_ds, deltheta | ||
| 430 | |||
| 431 | 9403985 | s = r**2 | |
| 432 | |||
| 433 |
1/2✓ Branch 0 taken 9403985 times.
✗ Branch 1 not taken.
|
9403985 | if (allocated(current_field)) then |
| 434 | 9403985 | call vmec_iota_interpolate_with_field(current_field, s, aiota, daiota_ds) | |
| 435 | else | ||
| 436 | ✗ | call vmec_iota_interpolate(s, aiota, daiota_ds) | |
| 437 | end if | ||
| 438 | |||
| 439 | 9403985 | vartheta = vartheta_c + aiota*y(1) | |
| 440 | 9403985 | varphi = varphi_c + y(1) | |
| 441 | |||
| 442 | ! Newton iteration to find field-specific theta from canonical theta | ||
| 443 |
1/2✓ Branch 0 taken 9403985 times.
✗ Branch 1 not taken.
|
9403985 | if (allocated(current_field)) then |
| 444 | 9403985 | theta = vartheta | |
| 445 | call newton_theta_from_canonical(current_field, s, vartheta, varphi, & | ||
| 446 | 9403985 | theta, converged) | |
| 447 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9403985 times.
|
9403985 | if (.not. converged) then |
| 448 | ✗ | print *, 'WARNING: Newton iteration failed in rhs_cancoord' | |
| 449 | end if | ||
| 450 | else | ||
| 451 | ✗ | theta = vartheta | |
| 452 | ✗ | do iter = 1, 100 | |
| 453 | ✗ | call vmec_lambda_interpolate(s, theta, varphi, alam, dl_dt) | |
| 454 | ✗ | deltheta = (vartheta - theta - alam)/(1.0_dp + dl_dt) | |
| 455 | ✗ | theta = theta + deltheta | |
| 456 | ✗ | if (abs(deltheta) < epserr) exit | |
| 457 | end do | ||
| 458 | end if | ||
| 459 | |||
| 460 |
2/2✓ Branch 0 taken 33 times.
✓ Branch 1 taken 9403952 times.
|
9403985 | if (onlytheta) return |
| 461 | |||
| 462 |
1/2✓ Branch 0 taken 9403952 times.
✗ Branch 1 not taken.
|
9403952 | if (allocated(current_field)) then |
| 463 | call vmec_field_evaluate_with_field(current_field, s, theta, varphi, & | ||
| 464 | A_theta, A_phi, dA_theta_ds, & | ||
| 465 | dA_phi_ds, aiota, & | ||
| 466 | sqg, alam, dl_ds, dl_dt, dl_dp, & | ||
| 467 | Bctrvr_vartheta, Bctrvr_varphi, & | ||
| 468 | Bcovar_r, Bcovar_vartheta, & | ||
| 469 | 9403952 | Bcovar_varphi) | |
| 470 | else | ||
| 471 | call vmec_field_evaluate(s, theta, varphi, A_theta, A_phi, & | ||
| 472 | dA_theta_ds, dA_phi_ds, aiota, & | ||
| 473 | sqg, alam, dl_ds, dl_dt, dl_dp, Bctrvr_vartheta, & | ||
| 474 | Bctrvr_varphi, & | ||
| 475 | ✗ | Bcovar_r, Bcovar_vartheta, Bcovar_varphi) | |
| 476 | end if | ||
| 477 | |||
| 478 | dy(1) = -(Bcovar_r + daiota_ds*Bcovar_vartheta*y(1))/ & | ||
| 479 | 9403952 | (aiota*Bcovar_vartheta + Bcovar_varphi) | |
| 480 | 9403952 | dy(1) = 2.0_dp*r*dy(1) | |
| 481 | |||
| 482 | 9403985 | end subroutine rhs_cancoord | |
| 483 | |||
| 484 | subroutine print_progress(message, progress, total) | ||
| 485 | character(*), intent(in) :: message | ||
| 486 | integer, intent(in) :: progress, total | ||
| 487 | |||
| 488 | #ifndef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 489 | return | ||
| 490 | #endif | ||
| 491 | |||
| 492 | write (*, '(A, I4, A, I4)', advance='no') message, progress, ' of ', total | ||
| 493 | |||
| 494 | if (progress < total) then | ||
| 495 | write (*, '(A)', advance="no") char(13) | ||
| 496 | else | ||
| 497 | write (*, *) | ||
| 498 | end if | ||
| 499 | 9403985 | end subroutine print_progress | |
| 500 | |||
| 501 | 1407529 | subroutine splint_can_coord(fullset, mode_secders, r, vartheta_c, varphi_c, & | |
| 502 | A_theta, A_phi, dA_theta_dr, dA_phi_dr, & | ||
| 503 | d2A_phi_dr2, d3A_phi_dr3, & | ||
| 504 | sqg_c, dsqg_c_dr, dsqg_c_dt, dsqg_c_dp, & | ||
| 505 | B_vartheta_c, dB_vartheta_c_dr, & | ||
| 506 | dB_vartheta_c_dt, dB_vartheta_c_dp, & | ||
| 507 | B_varphi_c, dB_varphi_c_dr, & | ||
| 508 | dB_varphi_c_dt, dB_varphi_c_dp, & | ||
| 509 | d2sqg_rr, d2sqg_rt, d2sqg_rp, & | ||
| 510 | d2sqg_tt, d2sqg_tp, d2sqg_pp, & | ||
| 511 | d2bth_rr, d2bth_rt, d2bth_rp, & | ||
| 512 | d2bth_tt, d2bth_tp, d2bth_pp, & | ||
| 513 | d2bph_rr, d2bph_rt, d2bph_rp, & | ||
| 514 | d2bph_tt, d2bph_tp, d2bph_pp, G_c) | ||
| 515 | |||
| 516 | use vector_potentail_mod, only: torflux | ||
| 517 | use new_vmec_stuff_mod, only: nper | ||
| 518 | use chamb_mod, only: rnegflag | ||
| 519 | use diag_mod, only: dodiag, icounter | ||
| 520 | |||
| 521 | implicit none | ||
| 522 | |||
| 523 | logical, intent(in) :: fullset | ||
| 524 | integer, intent(in) :: mode_secders | ||
| 525 | real(dp), intent(in) :: r | ||
| 526 | real(dp), intent(in) :: vartheta_c, varphi_c | ||
| 527 | |||
| 528 | real(dp), intent(out) :: A_phi, A_theta, dA_phi_dr, dA_theta_dr | ||
| 529 | real(dp), intent(out) :: d2A_phi_dr2, d3A_phi_dr3 | ||
| 530 | real(dp), intent(out) :: sqg_c, dsqg_c_dr, dsqg_c_dt, dsqg_c_dp | ||
| 531 | real(dp), intent(out) :: B_vartheta_c, dB_vartheta_c_dr | ||
| 532 | real(dp), intent(out) :: dB_vartheta_c_dt, dB_vartheta_c_dp | ||
| 533 | real(dp), intent(out) :: B_varphi_c, dB_varphi_c_dr | ||
| 534 | real(dp), intent(out) :: dB_varphi_c_dt, dB_varphi_c_dp | ||
| 535 | real(dp), intent(out) :: d2sqg_rr, d2sqg_rt, d2sqg_rp | ||
| 536 | real(dp), intent(out) :: d2sqg_tt, d2sqg_tp, d2sqg_pp | ||
| 537 | real(dp), intent(out) :: d2bth_rr, d2bth_rt, d2bth_rp | ||
| 538 | real(dp), intent(out) :: d2bth_tt, d2bth_tp, d2bth_pp | ||
| 539 | real(dp), intent(out) :: d2bph_rr, d2bph_rt, d2bph_rp | ||
| 540 | real(dp), intent(out) :: d2bph_tt, d2bph_tp, d2bph_pp | ||
| 541 | real(dp), intent(out) :: G_c | ||
| 542 | |||
| 543 | real(dp) :: r_eval | ||
| 544 | real(dp) :: rho_tor, drhods, drhods2, d2rhods2m | ||
| 545 | real(dp) :: x_eval(3) | ||
| 546 | real(dp) :: yq(1), dyq(3, 1), d2yq(6, 1) | ||
| 547 | real(dp) :: d2yq_rmix(3, 1) | ||
| 548 | real(dp) :: y_G(1), dy_G(3, 1) | ||
| 549 | real(dp) :: y1d(1), dy1d(1), d2y1d(1) | ||
| 550 | real(dp) :: theta_wrapped, phi_wrapped | ||
| 551 | real(dp) :: qua, dqua_dr, dqua_dt, dqua_dp | ||
| 552 | real(dp) :: d2qua_dr2, d2qua_drdt, d2qua_drdp | ||
| 553 | real(dp) :: d2qua_dt2, d2qua_dtdp, d2qua_dp2 | ||
| 554 | |||
| 555 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1407529 times.
|
1407529 | if (dodiag) then |
| 556 | ✗ | !$omp atomic | |
| 557 | icounter = icounter + 1 | ||
| 558 | end if | ||
| 559 | 1407529 | r_eval = r | |
| 560 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1407529 times.
|
1407529 | if (r_eval <= 0.0_dp) then |
| 561 | ✗ | rnegflag = .true. | |
| 562 | ✗ | r_eval = abs(r_eval) | |
| 563 | end if | ||
| 564 | |||
| 565 | 1407529 | A_theta = torflux*r_eval | |
| 566 | 1407529 | dA_theta_dr = torflux | |
| 567 | |||
| 568 | ! Interpolate A_phi using batch spline (1D) | ||
| 569 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1407529 times.
|
1407529 | if (.not. aphi_batch_spline_ready) then |
| 570 | ✗ | error stop "splint_can_coord: Aphi batch spline not initialized" | |
| 571 | end if | ||
| 572 | |||
| 573 | call evaluate_batch_splines_1d_der2(aphi_batch_spline, r_eval, & | ||
| 574 | 1407529 | y1d, dy1d, d2y1d) | |
| 575 | 1407529 | d3A_phi_dr3 = 0.0_dp | |
| 576 | 1407529 | A_phi = y1d(1) | |
| 577 | 1407529 | dA_phi_dr = dy1d(1) | |
| 578 | 1407529 | d2A_phi_dr2 = d2y1d(1) | |
| 579 | |||
| 580 | ! Prepare coordinates for 3D interpolation | ||
| 581 | 1407529 | rho_tor = sqrt(r_eval) | |
| 582 |
3/4✓ Branch 0 taken 1407529 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 519485 times.
✓ Branch 3 taken 888044 times.
|
1407529 | theta_wrapped = modulo(vartheta_c, TWOPI) |
| 583 |
3/4✓ Branch 0 taken 1407529 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 628601 times.
✓ Branch 3 taken 778928 times.
|
1407529 | phi_wrapped = modulo(varphi_c, TWOPI/real(nper, dp)) |
| 584 | |||
| 585 | 1407529 | x_eval(1) = rho_tor | |
| 586 | 1407529 | x_eval(2) = theta_wrapped | |
| 587 | 1407529 | x_eval(3) = phi_wrapped | |
| 588 | |||
| 589 | ! Chain rule coefficients for rho -> s conversion | ||
| 590 | ! rho = sqrt(s), drho/ds = 0.5/rho, d2rho/ds2 = -0.25/rho^3 | ||
| 591 | 1407529 | drhods = 0.5_dp/rho_tor | |
| 592 | 1407529 | drhods2 = drhods**2 | |
| 593 | 1407529 | d2rhods2m = drhods2/rho_tor ! -d2rho/ds2 (negated for chain rule) | |
| 594 | |||
| 595 | ! Interpolate G if needed | ||
| 596 |
2/2✓ Branch 0 taken 33 times.
✓ Branch 1 taken 1407496 times.
|
1407529 | if (fullset) then |
| 597 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 33 times.
|
33 | if (.not. G_batch_spline_ready) then |
| 598 | ✗ | error stop "splint_can_coord: G batch spline not initialized" | |
| 599 | end if | ||
| 600 | 33 | call evaluate_batch_splines_3d_der(G_batch_spline, x_eval, y_G, dy_G) | |
| 601 | 33 | G_c = y_G(1) | |
| 602 | else | ||
| 603 | 1407496 | G_c = 0.0_dp | |
| 604 | end if | ||
| 605 | |||
| 606 | ! Interpolate sqg, B_vartheta, B_varphi (separate NQ=1 splines) | ||
| 607 |
3/6✓ Branch 0 taken 1407529 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1407529 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1407529 times.
|
1407529 | if (.not. (sqg_batch_spline_ready .and. Bt_batch_spline_ready .and. & |
| 608 | Bp_batch_spline_ready)) then | ||
| 609 | ✗ | error stop "splint_can_coord: sqg/Bt/Bp batch splines not initialized" | |
| 610 | end if | ||
| 611 | |||
| 612 |
2/2✓ Branch 0 taken 1089188 times.
✓ Branch 1 taken 318341 times.
|
1407529 | if (mode_secders == 2 .or. mode_secders == 3) then |
| 613 |
1/2✓ Branch 0 taken 1089188 times.
✗ Branch 1 not taken.
|
1089188 | if (mode_secders == 2) then |
| 614 | call evaluate_batch_splines_3d_der2(sqg_batch_spline, x_eval, yq, & | ||
| 615 | 1089188 | dyq, d2yq) | |
| 616 | else | ||
| 617 | call evaluate_batch_splines_3d_der2_rmix(sqg_batch_spline, x_eval, yq, & | ||
| 618 | ✗ | dyq, d2yq_rmix) | |
| 619 | ✗ | d2yq(1:3, 1) = d2yq_rmix(:, 1) | |
| 620 | ✗ | d2yq(4:6, 1) = 0.0_dp | |
| 621 | end if | ||
| 622 | |||
| 623 | 1089188 | qua = yq(1) | |
| 624 | 1089188 | dqua_dr = dyq(1, 1) | |
| 625 | 1089188 | dqua_dt = dyq(2, 1) | |
| 626 | 1089188 | dqua_dp = dyq(3, 1) | |
| 627 | 1089188 | d2qua_dr2 = d2yq(1, 1) | |
| 628 | 1089188 | d2qua_drdt = d2yq(2, 1) | |
| 629 | 1089188 | d2qua_drdp = d2yq(3, 1) | |
| 630 | 1089188 | d2qua_dt2 = d2yq(4, 1) | |
| 631 | 1089188 | d2qua_dtdp = d2yq(5, 1) | |
| 632 | 1089188 | d2qua_dp2 = d2yq(6, 1) | |
| 633 | |||
| 634 | 1089188 | d2qua_dr2 = d2qua_dr2*drhods2 - dqua_dr*d2rhods2m | |
| 635 | 1089188 | dqua_dr = dqua_dr*drhods | |
| 636 | 1089188 | d2qua_drdt = d2qua_drdt*drhods | |
| 637 | 1089188 | d2qua_drdp = d2qua_drdp*drhods | |
| 638 | |||
| 639 | 1089188 | sqg_c = qua | |
| 640 | 1089188 | dsqg_c_dr = dqua_dr | |
| 641 | 1089188 | dsqg_c_dt = dqua_dt | |
| 642 | 1089188 | dsqg_c_dp = dqua_dp | |
| 643 | 1089188 | d2sqg_rr = d2qua_dr2 | |
| 644 | 1089188 | d2sqg_rt = d2qua_drdt | |
| 645 | 1089188 | d2sqg_rp = d2qua_drdp | |
| 646 | 1089188 | d2sqg_tt = d2qua_dt2 | |
| 647 | 1089188 | d2sqg_tp = d2qua_dtdp | |
| 648 | 1089188 | d2sqg_pp = d2qua_dp2 | |
| 649 | |||
| 650 |
1/2✓ Branch 0 taken 1089188 times.
✗ Branch 1 not taken.
|
1089188 | if (mode_secders == 2) then |
| 651 | call evaluate_batch_splines_3d_der2(Bt_batch_spline, x_eval, yq, & | ||
| 652 | 1089188 | dyq, d2yq) | |
| 653 | else | ||
| 654 | call evaluate_batch_splines_3d_der2_rmix(Bt_batch_spline, x_eval, yq, & | ||
| 655 | ✗ | dyq, d2yq_rmix) | |
| 656 | ✗ | d2yq(1:3, 1) = d2yq_rmix(:, 1) | |
| 657 | ✗ | d2yq(4:6, 1) = 0.0_dp | |
| 658 | end if | ||
| 659 | |||
| 660 | 1089188 | qua = yq(1) | |
| 661 | 1089188 | dqua_dr = dyq(1, 1) | |
| 662 | 1089188 | dqua_dt = dyq(2, 1) | |
| 663 | 1089188 | dqua_dp = dyq(3, 1) | |
| 664 | 1089188 | d2qua_dr2 = d2yq(1, 1) | |
| 665 | 1089188 | d2qua_drdt = d2yq(2, 1) | |
| 666 | 1089188 | d2qua_drdp = d2yq(3, 1) | |
| 667 | 1089188 | d2qua_dt2 = d2yq(4, 1) | |
| 668 | 1089188 | d2qua_dtdp = d2yq(5, 1) | |
| 669 | 1089188 | d2qua_dp2 = d2yq(6, 1) | |
| 670 | |||
| 671 | 1089188 | d2qua_dr2 = d2qua_dr2*drhods2 - dqua_dr*d2rhods2m | |
| 672 | 1089188 | dqua_dr = dqua_dr*drhods | |
| 673 | 1089188 | d2qua_drdt = d2qua_drdt*drhods | |
| 674 | 1089188 | d2qua_drdp = d2qua_drdp*drhods | |
| 675 | |||
| 676 | 1089188 | B_vartheta_c = qua | |
| 677 | 1089188 | dB_vartheta_c_dr = dqua_dr | |
| 678 | 1089188 | dB_vartheta_c_dt = dqua_dt | |
| 679 | 1089188 | dB_vartheta_c_dp = dqua_dp | |
| 680 | 1089188 | d2bth_rr = d2qua_dr2 | |
| 681 | 1089188 | d2bth_rt = d2qua_drdt | |
| 682 | 1089188 | d2bth_rp = d2qua_drdp | |
| 683 | 1089188 | d2bth_tt = d2qua_dt2 | |
| 684 | 1089188 | d2bth_tp = d2qua_dtdp | |
| 685 | 1089188 | d2bth_pp = d2qua_dp2 | |
| 686 | |||
| 687 |
1/2✓ Branch 0 taken 1089188 times.
✗ Branch 1 not taken.
|
1089188 | if (mode_secders == 2) then |
| 688 | call evaluate_batch_splines_3d_der2(Bp_batch_spline, x_eval, yq, & | ||
| 689 | 1089188 | dyq, d2yq) | |
| 690 | else | ||
| 691 | call evaluate_batch_splines_3d_der2_rmix(Bp_batch_spline, x_eval, yq, & | ||
| 692 | ✗ | dyq, d2yq_rmix) | |
| 693 | ✗ | d2yq(1:3, 1) = d2yq_rmix(:, 1) | |
| 694 | ✗ | d2yq(4:6, 1) = 0.0_dp | |
| 695 | end if | ||
| 696 | |||
| 697 | 1089188 | qua = yq(1) | |
| 698 | 1089188 | dqua_dr = dyq(1, 1) | |
| 699 | 1089188 | dqua_dt = dyq(2, 1) | |
| 700 | 1089188 | dqua_dp = dyq(3, 1) | |
| 701 | 1089188 | d2qua_dr2 = d2yq(1, 1) | |
| 702 | 1089188 | d2qua_drdt = d2yq(2, 1) | |
| 703 | 1089188 | d2qua_drdp = d2yq(3, 1) | |
| 704 | 1089188 | d2qua_dt2 = d2yq(4, 1) | |
| 705 | 1089188 | d2qua_dtdp = d2yq(5, 1) | |
| 706 | 1089188 | d2qua_dp2 = d2yq(6, 1) | |
| 707 | |||
| 708 | 1089188 | d2qua_dr2 = d2qua_dr2*drhods2 - dqua_dr*d2rhods2m | |
| 709 | 1089188 | dqua_dr = dqua_dr*drhods | |
| 710 | 1089188 | d2qua_drdt = d2qua_drdt*drhods | |
| 711 | 1089188 | d2qua_drdp = d2qua_drdp*drhods | |
| 712 | |||
| 713 | 1089188 | B_varphi_c = qua | |
| 714 | 1089188 | dB_varphi_c_dr = dqua_dr | |
| 715 | 1089188 | dB_varphi_c_dt = dqua_dt | |
| 716 | 1089188 | dB_varphi_c_dp = dqua_dp | |
| 717 | 1089188 | d2bph_rr = d2qua_dr2 | |
| 718 | 1089188 | d2bph_rt = d2qua_drdt | |
| 719 | 1089188 | d2bph_rp = d2qua_drdp | |
| 720 | 1089188 | d2bph_tt = d2qua_dt2 | |
| 721 | 1089188 | d2bph_tp = d2qua_dtdp | |
| 722 | 1089188 | d2bph_pp = d2qua_dp2 | |
| 723 | |||
| 724 | else | ||
| 725 | 318341 | call evaluate_batch_splines_3d_der(sqg_batch_spline, x_eval, yq, dyq) | |
| 726 | 318341 | sqg_c = yq(1) | |
| 727 | 318341 | dsqg_c_dr = dyq(1, 1)*drhods | |
| 728 | 318341 | dsqg_c_dt = dyq(2, 1) | |
| 729 | 318341 | dsqg_c_dp = dyq(3, 1) | |
| 730 | |||
| 731 | 318341 | call evaluate_batch_splines_3d_der(Bt_batch_spline, x_eval, yq, dyq) | |
| 732 | 318341 | B_vartheta_c = yq(1) | |
| 733 | 318341 | dB_vartheta_c_dr = dyq(1, 1)*drhods | |
| 734 | 318341 | dB_vartheta_c_dt = dyq(2, 1) | |
| 735 | 318341 | dB_vartheta_c_dp = dyq(3, 1) | |
| 736 | |||
| 737 | 318341 | call evaluate_batch_splines_3d_der(Bp_batch_spline, x_eval, yq, dyq) | |
| 738 | 318341 | B_varphi_c = yq(1) | |
| 739 | 318341 | dB_varphi_c_dr = dyq(1, 1)*drhods | |
| 740 | 318341 | dB_varphi_c_dt = dyq(2, 1) | |
| 741 | 318341 | dB_varphi_c_dp = dyq(3, 1) | |
| 742 | |||
| 743 | 318341 | d2sqg_rr = 0.0_dp | |
| 744 | 318341 | d2sqg_rt = 0.0_dp | |
| 745 | 318341 | d2sqg_rp = 0.0_dp | |
| 746 | 318341 | d2sqg_tt = 0.0_dp | |
| 747 | 318341 | d2sqg_tp = 0.0_dp | |
| 748 | 318341 | d2sqg_pp = 0.0_dp | |
| 749 | 318341 | d2bth_rr = 0.0_dp | |
| 750 | 318341 | d2bth_rt = 0.0_dp | |
| 751 | 318341 | d2bth_rp = 0.0_dp | |
| 752 | 318341 | d2bth_tt = 0.0_dp | |
| 753 | 318341 | d2bth_tp = 0.0_dp | |
| 754 | 318341 | d2bth_pp = 0.0_dp | |
| 755 | 318341 | d2bph_rr = 0.0_dp | |
| 756 | 318341 | d2bph_rt = 0.0_dp | |
| 757 | 318341 | d2bph_rp = 0.0_dp | |
| 758 | 318341 | d2bph_tt = 0.0_dp | |
| 759 | 318341 | d2bph_tp = 0.0_dp | |
| 760 | 318341 | d2bph_pp = 0.0_dp | |
| 761 | |||
| 762 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 318341 times.
|
318341 | if (mode_secders == 1) then |
| 763 | call evaluate_batch_splines_3d_der2(sqg_batch_spline, x_eval, yq, dyq, & | ||
| 764 | ✗ | d2yq) | |
| 765 | ✗ | d2sqg_rr = d2yq(1, 1)*drhods2 - dyq(1, 1)*d2rhods2m | |
| 766 | |||
| 767 | call evaluate_batch_splines_3d_der2(Bt_batch_spline, x_eval, yq, dyq, & | ||
| 768 | ✗ | d2yq) | |
| 769 | ✗ | d2bth_rr = d2yq(1, 1)*drhods2 - dyq(1, 1)*d2rhods2m | |
| 770 | |||
| 771 | call evaluate_batch_splines_3d_der2(Bp_batch_spline, x_eval, yq, dyq, & | ||
| 772 | ✗ | d2yq) | |
| 773 | ✗ | d2bph_rr = d2yq(1, 1)*drhods2 - dyq(1, 1)*d2rhods2m | |
| 774 | end if | ||
| 775 | end if | ||
| 776 | |||
| 777 | 1407529 | end subroutine splint_can_coord | |
| 778 | |||
| 779 | 33 | subroutine can_to_vmec(r, vartheta_c_in, varphi_c_in, theta_vmec, varphi_vmec) | |
| 780 | use exchange_get_cancoord_mod, only: vartheta_c, varphi_c, theta | ||
| 781 | |||
| 782 | implicit none | ||
| 783 | |||
| 784 | real(dp), intent(in) :: r, vartheta_c_in, varphi_c_in | ||
| 785 | real(dp), intent(out) :: theta_vmec, varphi_vmec | ||
| 786 | |||
| 787 | logical :: fullset | ||
| 788 | integer :: mode_secders | ||
| 789 | real(dp) :: r_local | ||
| 790 | real(dp) :: A_phi, A_theta, dA_phi_dr, dA_theta_dr, d2A_phi_dr2, d3A_phi_dr3 | ||
| 791 | real(dp) :: sqg_c, dsqg_c_dr, dsqg_c_dt, dsqg_c_dp | ||
| 792 | real(dp) :: B_vartheta_c, dB_vartheta_c_dr, dB_vartheta_c_dt, dB_vartheta_c_dp | ||
| 793 | real(dp) :: B_varphi_c, dB_varphi_c_dr, dB_varphi_c_dt, dB_varphi_c_dp | ||
| 794 | real(dp) :: G_c | ||
| 795 | real(dp) :: d2sqg_rr, d2sqg_rt, d2sqg_rp, d2sqg_tt, d2sqg_tp, d2sqg_pp | ||
| 796 | real(dp) :: d2bth_rr, d2bth_rt, d2bth_rp, d2bth_tt, d2bth_tp, d2bth_pp | ||
| 797 | real(dp) :: d2bph_rr, d2bph_rt, d2bph_rp, d2bph_tt, d2bph_tp, d2bph_pp | ||
| 798 | real(dp), dimension(1) :: y, dy | ||
| 799 | |||
| 800 | 33 | fullset = .true. | |
| 801 | 33 | mode_secders = 0 | |
| 802 | 33 | r_local = r | |
| 803 | |||
| 804 | call splint_can_coord(fullset, mode_secders, r_local, vartheta_c_in, & | ||
| 805 | varphi_c_in, A_theta, A_phi, dA_theta_dr, dA_phi_dr, & | ||
| 806 | d2A_phi_dr2, d3A_phi_dr3, & | ||
| 807 | sqg_c, dsqg_c_dr, dsqg_c_dt, dsqg_c_dp, & | ||
| 808 | B_vartheta_c, dB_vartheta_c_dr, dB_vartheta_c_dt, & | ||
| 809 | dB_vartheta_c_dp, & | ||
| 810 | B_varphi_c, dB_varphi_c_dr, dB_varphi_c_dt, & | ||
| 811 | dB_varphi_c_dp, & | ||
| 812 | d2sqg_rr, d2sqg_rt, d2sqg_rp, d2sqg_tt, d2sqg_tp, & | ||
| 813 | d2sqg_pp, & | ||
| 814 | d2bth_rr, d2bth_rt, d2bth_rp, d2bth_tt, d2bth_tp, & | ||
| 815 | d2bth_pp, & | ||
| 816 | d2bph_rr, d2bph_rt, d2bph_rp, d2bph_tt, d2bph_tp, & | ||
| 817 | 33 | d2bph_pp, G_c) | |
| 818 | |||
| 819 | 33 | vartheta_c = vartheta_c_in | |
| 820 | 33 | varphi_c = varphi_c_in | |
| 821 | 33 | y(1) = G_c | |
| 822 | |||
| 823 | ! Transform from r (toroidal flux) to rho_tor for ODE integration | ||
| 824 | 33 | call rhs_cancoord(sqrt(r_local), y, dy) | |
| 825 | |||
| 826 | 33 | theta_vmec = theta | |
| 827 | 33 | varphi_vmec = varphi_c_in + G_c | |
| 828 | |||
| 829 | 33 | end subroutine can_to_vmec | |
| 830 | |||
| 831 | ✗ | subroutine deallocate_can_coord | |
| 832 | ✗ | call reset_canflux_batch_splines | |
| 833 | ✗ | end subroutine deallocate_can_coord | |
| 834 | |||
| 835 | 40 | subroutine vmec_to_can(r, theta, varphi, vartheta_c, varphi_c) | |
| 836 | ! Input : r,theta,varphi - VMEC coordinates | ||
| 837 | ! Output: vartheta_c,varphi_c - canonical coordinates | ||
| 838 | |||
| 839 | use spline_vmec_sub | ||
| 840 | use new_vmec_stuff_mod, only: nper | ||
| 841 | use vector_potentail_mod, only: torflux | ||
| 842 | use chamb_mod, only: rnegflag | ||
| 843 | use vmec_field_eval | ||
| 844 | |||
| 845 | implicit none | ||
| 846 | |||
| 847 | real(dp), parameter :: epserr = 1.0e-14_dp | ||
| 848 | integer, parameter :: niter = 100 | ||
| 849 | integer :: iter | ||
| 850 | real(dp), intent(in) :: r, theta, varphi | ||
| 851 | real(dp), intent(out) :: vartheta_c, varphi_c | ||
| 852 | real(dp) :: delthe, delphi, alam, dl_dt, vartheta | ||
| 853 | real(dp) :: rho_tor, x_eval(3), y_G(1), dy_G(3, 1) | ||
| 854 | real(dp) :: G_c, dG_c_dt, dG_c_dp, aiota | ||
| 855 | real(dp) :: ts, ps, dts_dtc, dts_dpc, dps_dtc, dps_dpc, det | ||
| 856 | real(dp) :: y1d(1), dy1d(1), d2y1d(1), dA_phi_dr, dA_theta_dr | ||
| 857 | real(dp) :: r_local | ||
| 858 | |||
| 859 | 40 | r_local = r | |
| 860 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
|
40 | if (r_local <= 0.0_dp) then |
| 861 | ✗ | rnegflag = .true. | |
| 862 | ✗ | r_local = abs(r_local) | |
| 863 | end if | ||
| 864 | |||
| 865 |
1/2✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
|
40 | if (allocated(current_field)) then |
| 866 | call vmec_lambda_interpolate_with_field(current_field, r_local, theta, & | ||
| 867 | 40 | varphi, alam, dl_dt) | |
| 868 | else | ||
| 869 | ✗ | call vmec_lambda_interpolate(r_local, theta, varphi, alam, dl_dt) | |
| 870 | end if | ||
| 871 | |||
| 872 | 40 | vartheta = theta + alam | |
| 873 | |||
| 874 | 40 | vartheta_c = vartheta | |
| 875 | 40 | varphi_c = varphi | |
| 876 | |||
| 877 | ! Get iota from A_phi interpolation | ||
| 878 | 40 | dA_theta_dr = torflux | |
| 879 | call evaluate_batch_splines_1d_der2(aphi_batch_spline, r_local, y1d, & | ||
| 880 | 40 | dy1d, d2y1d) | |
| 881 | 40 | dA_phi_dr = dy1d(1) | |
| 882 | 40 | aiota = -dA_phi_dr/dA_theta_dr | |
| 883 | |||
| 884 |
2/2✓ Branch 0 taken 259 times.
✓ Branch 1 taken 1 time.
|
260 | do iter = 1, niter |
| 885 | 259 | rho_tor = sqrt(r_local) | |
| 886 | 259 | x_eval(1) = rho_tor | |
| 887 |
2/4✓ Branch 0 taken 259 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 259 times.
|
259 | x_eval(2) = modulo(vartheta_c, TWOPI) |
| 888 |
2/4✓ Branch 0 taken 259 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 259 times.
|
259 | x_eval(3) = modulo(varphi_c, TWOPI/real(nper, dp)) |
| 889 | |||
| 890 | 259 | call evaluate_batch_splines_3d_der(G_batch_spline, x_eval, y_G, dy_G) | |
| 891 | 259 | G_c = y_G(1) | |
| 892 | 259 | dG_c_dt = dy_G(2, 1) | |
| 893 | 259 | dG_c_dp = dy_G(3, 1) | |
| 894 | |||
| 895 | 259 | ts = vartheta_c + aiota*G_c - vartheta | |
| 896 | 259 | ps = varphi_c + G_c - varphi | |
| 897 | 259 | dts_dtc = 1.0_dp + aiota*dG_c_dt | |
| 898 | 259 | dts_dpc = aiota*dG_c_dp | |
| 899 | 259 | dps_dtc = dG_c_dt | |
| 900 | 259 | dps_dpc = 1.0_dp + dG_c_dp | |
| 901 | 259 | det = 1.0_dp + aiota*dG_c_dt + dG_c_dp | |
| 902 | |||
| 903 | 259 | delthe = (ps*dts_dpc - ts*dps_dpc)/det | |
| 904 | 259 | delphi = (ts*dps_dtc - ps*dts_dtc)/det | |
| 905 | |||
| 906 | 259 | vartheta_c = vartheta_c + delthe | |
| 907 | 259 | varphi_c = varphi_c + delphi | |
| 908 |
2/2✓ Branch 0 taken 220 times.
✓ Branch 1 taken 39 times.
|
260 | if (abs(delthe) + abs(delphi) < epserr) exit |
| 909 | end do | ||
| 910 | |||
| 911 | 40 | end subroutine vmec_to_can | |
| 912 | |||
| 913 | ✗ | subroutine vmec_to_cyl(s, theta, varphi, Rcyl, Zcyl) | |
| 914 | 40 | use spline_vmec_sub | |
| 915 | use vmec_field_eval | ||
| 916 | |||
| 917 | real(dp), intent(in) :: s, theta, varphi | ||
| 918 | real(dp), intent(out) :: Rcyl, Zcyl | ||
| 919 | |||
| 920 | real(dp) :: A_phi, A_theta, dA_phi_ds, dA_theta_ds, aiota, & | ||
| 921 | R, Z, alam, dR_ds, dR_dt, dR_dp, dZ_ds, dZ_dt, dZ_dp, dl_ds, & | ||
| 922 | dl_dt, dl_dp | ||
| 923 | |||
| 924 | ✗ | if (allocated(current_field)) then | |
| 925 | call vmec_data_interpolate_with_field(current_field, s, theta, varphi, & | ||
| 926 | A_phi, A_theta, dA_phi_ds, & | ||
| 927 | dA_theta_ds, aiota, & | ||
| 928 | R, Z, alam, dR_ds, dR_dt, dR_dp, & | ||
| 929 | dZ_ds, dZ_dt, dZ_dp, & | ||
| 930 | ✗ | dl_ds, dl_dt, dl_dp) | |
| 931 | else | ||
| 932 | call vmec_data_interpolate(s, theta, varphi, A_phi, A_theta, & | ||
| 933 | dA_phi_ds, dA_theta_ds, aiota, & | ||
| 934 | R, Z, alam, dR_ds, dR_dt, dR_dp, dZ_ds, & | ||
| 935 | dZ_dt, dZ_dp, & | ||
| 936 | ✗ | dl_ds, dl_dt, dl_dp) | |
| 937 | end if | ||
| 938 | |||
| 939 | ✗ | Rcyl = R | |
| 940 | ✗ | Zcyl = Z | |
| 941 | ✗ | end subroutine vmec_to_cyl | |
| 942 | |||
| 943 | ✗ | end module get_can_sub | |
| 944 |