| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module field_can_boozer | ||
| 2 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 3 | use field_can_base, only: apply_radial_electric_potential, & | ||
| 4 | count_field_evaluation, field_can_t, twopi | ||
| 5 | |||
| 6 | implicit none | ||
| 7 | |||
| 8 | public :: integ_to_ref_boozer_chartmap, ref_to_integ_boozer_chartmap | ||
| 9 | public :: evaluate_boozer_regular | ||
| 10 | |||
| 11 | contains | ||
| 12 | |||
| 13 | 122 | subroutine evaluate_boozer_regular(x, y, phi, mode, rf) | |
| 14 | !> Regular near-axis Boozer field at X = sqrt(s) cos(theta_B), | ||
| 15 | !> Y = sqrt(s) sin(theta_B), phi_B (libneo boozer_regular). mode 0, 1, 2 | ||
| 16 | !> returns values, first and second derivatives in (X, Y, phi); every | ||
| 17 | !> component is C^2 through the magnetic axis. Requires the field to be | ||
| 18 | !> initialised with field_can_mod regular_axis_field set. | ||
| 19 | use boozer_regular, only: regular_field_t, eval_boozer_regular | ||
| 20 | real(dp), intent(in) :: x, y, phi | ||
| 21 | integer, intent(in) :: mode | ||
| 22 | type(regular_field_t), intent(inout) :: rf | ||
| 23 | |||
| 24 | 122 | call eval_boozer_regular(x, y, phi, mode, rf) | |
| 25 | |||
| 26 | ! Only mode 2 carries second derivatives, the d2 class of the counter. | ||
| 27 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 120 times.
|
122 | if (mode >= 2) then |
| 28 | 2 | call count_field_evaluation(1) | |
| 29 | else | ||
| 30 | 120 | call count_field_evaluation(0) | |
| 31 | end if | ||
| 32 | 122 | end subroutine evaluate_boozer_regular | |
| 33 | |||
| 34 | 100174 | subroutine integ_to_ref_boozer_chartmap(xinteg, xref) | |
| 35 | !> Convert integrator coords (s, theta_B, phi_B) to | ||
| 36 | !> chartmap reference coords (rho, theta_B, phi_B). | ||
| 37 | real(dp), intent(in) :: xinteg(3) | ||
| 38 | real(dp), intent(out) :: xref(3) | ||
| 39 | |||
| 40 | 100174 | xref(1) = sqrt(max(xinteg(1), 0.0_dp)) | |
| 41 | 100174 | xref(2) = xinteg(2) | |
| 42 | 100174 | xref(3) = xinteg(3) | |
| 43 | 100174 | end subroutine integ_to_ref_boozer_chartmap | |
| 44 | |||
| 45 | 100026 | subroutine ref_to_integ_boozer_chartmap(xref, xinteg) | |
| 46 | !> Convert chartmap reference coords (rho, theta_B, phi_B) to | ||
| 47 | !> integrator coords (s, theta_B, phi_B). | ||
| 48 | real(dp), intent(in) :: xref(3) | ||
| 49 | real(dp), intent(out) :: xinteg(3) | ||
| 50 | |||
| 51 | 100026 | xinteg(1) = xref(1)**2 | |
| 52 | 100026 | xinteg(2) = xref(2) | |
| 53 | 100026 | xinteg(3) = xref(3) | |
| 54 | 100026 | end subroutine ref_to_integ_boozer_chartmap | |
| 55 | |||
| 56 | 20335296 | subroutine evaluate_boozer(f, r, th_c, ph_c, mode_secders) | |
| 57 | type(field_can_t), intent(inout) :: f | ||
| 58 | real(dp), intent(in) :: r, th_c, ph_c | ||
| 59 | integer, intent(in) :: mode_secders | ||
| 60 | |||
| 61 | 20335296 | call eval_field_booz(f, r, th_c, ph_c, mode_secders) | |
| 62 | 20335296 | call apply_radial_electric_potential(f, r) | |
| 63 | |||
| 64 | 20335296 | call count_field_evaluation(mode_secders) | |
| 65 | 20335296 | end subroutine evaluate_boozer | |
| 66 | |||
| 67 | 320193 | subroutine integ_to_ref_boozer(xinteg, xref) | |
| 68 | use boozer_sub, only: boozer_to_vmec | ||
| 69 | real(dp), intent(in) :: xinteg(3) | ||
| 70 | real(dp), intent(out) :: xref(3) | ||
| 71 | |||
| 72 | 320193 | xref(1) = xinteg(1) | |
| 73 | call boozer_to_vmec( & | ||
| 74 | 320193 | xinteg(1), mod(xinteg(2), twopi), mod(xinteg(3), twopi), xref(2), xref(3)) | |
| 75 | 320193 | xref(2) = mod(xref(2), twopi) | |
| 76 | 320193 | xref(3) = mod(xref(3), twopi) | |
| 77 | 320193 | end subroutine integ_to_ref_boozer | |
| 78 | |||
| 79 | 179 | subroutine ref_to_integ_boozer(xref, xinteg) | |
| 80 | use boozer_sub, only: vmec_to_boozer | ||
| 81 | real(dp), intent(in) :: xref(3) | ||
| 82 | real(dp), intent(out) :: xinteg(3) | ||
| 83 | |||
| 84 | 179 | xinteg(1) = xref(1) | |
| 85 | call vmec_to_boozer( & | ||
| 86 | 179 | xref(1), mod(xref(2), twopi), mod(xref(3), twopi), xinteg(2), xinteg(3)) | |
| 87 | 179 | xinteg(2) = mod(xinteg(2), twopi) | |
| 88 | 179 | xinteg(3) = mod(xinteg(3), twopi) | |
| 89 | 179 | end subroutine ref_to_integ_boozer | |
| 90 | |||
| 91 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 92 | ! | ||
| 93 | 20335312 | subroutine eval_field_booz(f, r, th_c, ph_c, mode_secders) | |
| 94 | use boozer_sub, only: splint_boozer_coord, boozer_state | ||
| 95 | ! | ||
| 96 | ! Evaluates magnetic field in Boozer canonical coordinates (r, th_c, ph_c) | ||
| 97 | ! and stores results in variable f | ||
| 98 | ! Works for A_th linear in r (toroidal flux as radial variable) | ||
| 99 | ! | ||
| 100 | ! Sign convention: A_theta = torflux*r is the toroidal flux psi_tor (on the | ||
| 101 | ! poloidal covariant component), with torflux = -phi_vmec/(2*pi) > 0 since | ||
| 102 | ! VMEC is left-handed (signgs = -1). A_phi = -chi is the toroidal covariant | ||
| 103 | ! component carrying the poloidal flux. See libneo spline_vmec_data / | ||
| 104 | ! vmecinm_m and gvec_simple/CONVENTIONS.md. | ||
| 105 | ! | ||
| 106 | ! mode_secders = 0: no second derivatives | ||
| 107 | ! mode_secders = 1: second derivatives only in d/dr^2 | ||
| 108 | ! mode_secders = 2: all second derivatives, including mixed | ||
| 109 | ! | ||
| 110 | |||
| 111 | type(field_can_t), intent(inout) :: f | ||
| 112 | double precision, intent(in) :: r, th_c, ph_c | ||
| 113 | integer, intent(in) :: mode_secders | ||
| 114 | |||
| 115 | double precision :: Bctr_vartheta, Bctr_varphi, bmod2, sqg, & | ||
| 116 | d3Aphdr3, dummy, dummy3(3), dummy6(6), & | ||
| 117 | Bth, Bph, dBth, dBph, d2Bth, d2Bph | ||
| 118 | |||
| 119 | ! Count evaluations for profiling | ||
| 120 | |||
| 121 | ! initialize to zero - no angular derivatives will be set due to straight field | ||
| 122 | ! line Ath(r) Aph(r) | ||
| 123 |
2/2✓ Branch 0 taken 20335312 times.
✓ Branch 1 taken 61005936 times.
|
81341248 | f%dAth = 0d0 |
| 124 |
2/2✓ Branch 0 taken 20335312 times.
✓ Branch 1 taken 61005936 times.
|
81341248 | f%dAph = 0d0 |
| 125 | |||
| 126 | ! initialize all 2nd derivatives to zero, as the mode decides, which one to use | ||
| 127 |
2/2✓ Branch 0 taken 20335312 times.
✓ Branch 1 taken 122011872 times.
|
142347184 | f%d2Ath = 0d0 |
| 128 |
2/2✓ Branch 0 taken 20335312 times.
✓ Branch 1 taken 122011872 times.
|
142347184 | f%d2Aph = 0d0 |
| 129 |
2/2✓ Branch 0 taken 20335312 times.
✓ Branch 1 taken 122011872 times.
|
142347184 | f%d2hth = 0d0 |
| 130 |
2/2✓ Branch 0 taken 20335312 times.
✓ Branch 1 taken 122011872 times.
|
142347184 | f%d2hph = 0d0 |
| 131 |
2/2✓ Branch 0 taken 122011872 times.
✓ Branch 1 taken 20335312 times.
|
142347184 | f%d2Bmod = 0d0 |
| 132 | |||
| 133 | call splint_boozer_coord(r, th_c, ph_c, mode_secders, & | ||
| 134 | f%Ath, f%Aph, f%dAth(1), f%dAph(1), f%d2Aph(1), & | ||
| 135 | d3Aphdr3, & | ||
| 136 | Bth, dBth, d2Bth, & | ||
| 137 | Bph, dBph, d2Bph, & | ||
| 138 | 20335312 | f%Bmod, f%dBmod, f%d2Bmod, dummy, dummy3, dummy6) | |
| 139 | |||
| 140 | 20335312 | bmod2 = f%Bmod**2 | |
| 141 | 20335312 | sqg = (-f%dAph(1)/f%dAth(1)*Bth + Bph)/bmod2*boozer_state%torflux | |
| 142 | 20335312 | Bctr_vartheta = -f%dAph(1)/sqg | |
| 143 | 20335312 | Bctr_varphi = f%dAth(1)/sqg | |
| 144 | |||
| 145 | 20335312 | f%hth = Bth/f%Bmod | |
| 146 | 20335312 | f%hph = Bph/f%Bmod | |
| 147 | 20335312 | f%dhth(1) = dBth/f%Bmod - Bth*f%dBmod(1)/bmod2 | |
| 148 | 20335312 | f%dhph(1) = dBph/f%Bmod - Bph*f%dBmod(1)/bmod2 | |
| 149 |
2/2✓ Branch 0 taken 20335312 times.
✓ Branch 1 taken 40670624 times.
|
61005936 | f%dhth(2:3) = -Bth*f%dBmod(2:3)/bmod2 |
| 150 |
2/2✓ Branch 0 taken 40670624 times.
✓ Branch 1 taken 20335312 times.
|
61005936 | f%dhph(2:3) = -Bph*f%dBmod(2:3)/bmod2 |
| 151 | |||
| 152 |
2/2✓ Branch 0 taken 13096299 times.
✓ Branch 1 taken 7239013 times.
|
20335312 | if (mode_secders > 0) then |
| 153 | f%d2hth(1) = d2Bth/f%Bmod - 2d0*dBth*f%dBmod(1)/bmod2 + & | ||
| 154 | 13096299 | Bth/bmod2*(2d0*f%dBmod(1)**2/f%Bmod - f%d2Bmod(1)) | |
| 155 | f%d2hph(1) = d2Bph/f%Bmod - 2d0*dBph*f%dBmod(1)/bmod2 + & | ||
| 156 | 13096299 | Bph/bmod2*(2d0*f%dBmod(1)**2/f%Bmod - f%d2Bmod(1)) | |
| 157 | end if | ||
| 158 | |||
| 159 |
2/2✓ Branch 0 taken 13096299 times.
✓ Branch 1 taken 7239013 times.
|
20335312 | if (mode_secders .eq. 2) then |
| 160 | f%d2hth((/4, 6/)) = Bth/bmod2*(2d0*f%dBmod((/2, 3/))**2/f%Bmod - & | ||
| 161 |
2/2✓ Branch 0 taken 26192598 times.
✓ Branch 1 taken 13096299 times.
|
39288897 | f%d2Bmod((/4, 6/))) |
| 162 | f%d2hph((/4, 6/)) = Bph/bmod2*(2d0*f%dBmod((/2, 3/))**2/f%Bmod - & | ||
| 163 |
2/2✓ Branch 0 taken 26192598 times.
✓ Branch 1 taken 13096299 times.
|
39288897 | f%d2Bmod((/4, 6/))) |
| 164 | |||
| 165 | f%d2hth(2) = -dBth*f%dBmod(2)/bmod2 & | ||
| 166 | 13096299 | + Bth/bmod2*(2d0*f%dBmod(1)*f%dBmod(2)/f%Bmod - f%d2Bmod(2)) | |
| 167 | f%d2hph(2) = -dBph*f%dBmod(2)/bmod2 & | ||
| 168 | 13096299 | + Bph/bmod2*(2d0*f%dBmod(1)*f%dBmod(2)/f%Bmod - f%d2Bmod(2)) | |
| 169 | |||
| 170 | f%d2hth(3) = -dBth*f%dBmod(3)/bmod2 & | ||
| 171 | 13096299 | + Bth/bmod2*(2d0*f%dBmod(1)*f%dBmod(3)/f%Bmod - f%d2Bmod(3)) | |
| 172 | f%d2hph(3) = -dBph*f%dBmod(3)/bmod2 & | ||
| 173 | 13096299 | + Bph/bmod2*(2d0*f%dBmod(1)*f%dBmod(3)/f%Bmod - f%d2Bmod(3)) | |
| 174 | |||
| 175 | 13096299 | f%d2hth(5) = Bth/bmod2*(2d0*f%dBmod(2)*f%dBmod(3)/f%Bmod - f%d2Bmod(5)) | |
| 176 | 13096299 | f%d2hph(5) = Bph/bmod2*(2d0*f%dBmod(2)*f%dBmod(3)/f%Bmod - f%d2Bmod(5)) | |
| 177 | end if | ||
| 178 | 20335312 | end subroutine eval_field_booz | |
| 179 | |||
| 180 | 314696 | subroutine eval_field_booz_device(f, r, th_c, ph_c, mode_secders) | |
| 181 | ! Device-only counterpart of eval_field_booz. Initialization checks and | ||
| 182 | ! host diagnostic counters stay in splint_boozer_coord's host wrapper. | ||
| 183 | ! BOOZER_SECDERS_RADIAL_MIXED computes only rr, r-theta, and r-phi. | ||
| 184 | use boozer_sub, only: splint_boozer_coord_device, boozer_state, & | ||
| 185 | BOOZER_SECDERS_RADIAL_MIXED | ||
| 186 | |||
| 187 | type(field_can_t), intent(inout) :: f | ||
| 188 | real(dp), intent(in) :: r, th_c, ph_c | ||
| 189 | integer, intent(in) :: mode_secders | ||
| 190 | !$acc routine seq | ||
| 191 | |||
| 192 | real(dp) :: bmod2, & | ||
| 193 | d3Aphdr3, dummy, dummy3(3), dummy6(6), & | ||
| 194 | Bth, Bph, dBth, dBph, d2Bth, d2Bph | ||
| 195 | |||
| 196 |
2/2✓ Branch 0 taken 314696 times.
✓ Branch 1 taken 944088 times.
|
1258784 | f%dAth = 0.0_dp |
| 197 |
2/2✓ Branch 0 taken 944088 times.
✓ Branch 1 taken 314696 times.
|
1258784 | f%dAph = 0.0_dp |
| 198 |
2/2✓ Branch 0 taken 231476 times.
✓ Branch 1 taken 83220 times.
|
314696 | if (mode_secders == BOOZER_SECDERS_RADIAL_MIXED) then |
| 199 |
2/2✓ Branch 0 taken 83220 times.
✓ Branch 1 taken 249660 times.
|
332880 | f%d2Ath(1:3) = 0.0_dp |
| 200 |
2/2✓ Branch 0 taken 83220 times.
✓ Branch 1 taken 249660 times.
|
332880 | f%d2Aph(1:3) = 0.0_dp |
| 201 |
2/2✓ Branch 0 taken 83220 times.
✓ Branch 1 taken 249660 times.
|
332880 | f%d2hth(1:3) = 0.0_dp |
| 202 |
2/2✓ Branch 0 taken 249660 times.
✓ Branch 1 taken 83220 times.
|
332880 | f%d2hph(1:3) = 0.0_dp |
| 203 | else | ||
| 204 |
2/2✓ Branch 0 taken 231476 times.
✓ Branch 1 taken 1388856 times.
|
1620332 | f%d2Ath = 0.0_dp |
| 205 |
2/2✓ Branch 0 taken 231476 times.
✓ Branch 1 taken 1388856 times.
|
1620332 | f%d2Aph = 0.0_dp |
| 206 |
2/2✓ Branch 0 taken 231476 times.
✓ Branch 1 taken 1388856 times.
|
1620332 | f%d2hth = 0.0_dp |
| 207 |
2/2✓ Branch 0 taken 1388856 times.
✓ Branch 1 taken 231476 times.
|
1620332 | f%d2hph = 0.0_dp |
| 208 | end if | ||
| 209 | |||
| 210 | call splint_boozer_coord_device(r, th_c, ph_c, mode_secders, & | ||
| 211 | f%Ath, f%Aph, f%dAth(1), f%dAph(1), & | ||
| 212 | f%d2Aph(1), d3Aphdr3, & | ||
| 213 | Bth, dBth, d2Bth, Bph, dBph, d2Bph, & | ||
| 214 | f%Bmod, f%dBmod, f%d2Bmod, & | ||
| 215 | 314696 | dummy, dummy3, dummy6) | |
| 216 | |||
| 217 | 314696 | bmod2 = f%Bmod**2 | |
| 218 | 314696 | f%hth = Bth/f%Bmod | |
| 219 | 314696 | f%hph = Bph/f%Bmod | |
| 220 | 314696 | f%dhth(1) = dBth/f%Bmod - Bth*f%dBmod(1)/bmod2 | |
| 221 | 314696 | f%dhph(1) = dBph/f%Bmod - Bph*f%dBmod(1)/bmod2 | |
| 222 |
2/2✓ Branch 0 taken 314696 times.
✓ Branch 1 taken 629392 times.
|
944088 | f%dhth(2:3) = -Bth*f%dBmod(2:3)/bmod2 |
| 223 |
2/2✓ Branch 0 taken 629392 times.
✓ Branch 1 taken 314696 times.
|
944088 | f%dhph(2:3) = -Bph*f%dBmod(2:3)/bmod2 |
| 224 | |||
| 225 |
2/2✓ Branch 0 taken 178158 times.
✓ Branch 1 taken 136538 times.
|
314696 | if (mode_secders > 0) then |
| 226 | f%d2hth(1) = d2Bth/f%Bmod - 2.0_dp*dBth*f%dBmod(1)/bmod2 + & | ||
| 227 | 178158 | Bth/bmod2*(2.0_dp*f%dBmod(1)**2/f%Bmod - f%d2Bmod(1)) | |
| 228 | f%d2hph(1) = d2Bph/f%Bmod - 2.0_dp*dBph*f%dBmod(1)/bmod2 + & | ||
| 229 | 178158 | Bph/bmod2*(2.0_dp*f%dBmod(1)**2/f%Bmod - f%d2Bmod(1)) | |
| 230 | end if | ||
| 231 | |||
| 232 |
2/2✓ Branch 0 taken 178158 times.
✓ Branch 1 taken 136538 times.
|
314696 | if (mode_secders == 2 .or. & |
| 233 | mode_secders == BOOZER_SECDERS_RADIAL_MIXED) then | ||
| 234 | f%d2hth(2) = -dBth*f%dBmod(2)/bmod2 + & | ||
| 235 | 178158 | Bth/bmod2*(2.0_dp*f%dBmod(1)*f%dBmod(2)/f%Bmod - f%d2Bmod(2)) | |
| 236 | f%d2hph(2) = -dBph*f%dBmod(2)/bmod2 + & | ||
| 237 | 178158 | Bph/bmod2*(2.0_dp*f%dBmod(1)*f%dBmod(2)/f%Bmod - f%d2Bmod(2)) | |
| 238 | |||
| 239 | f%d2hth(3) = -dBth*f%dBmod(3)/bmod2 + & | ||
| 240 | 178158 | Bth/bmod2*(2.0_dp*f%dBmod(1)*f%dBmod(3)/f%Bmod - f%d2Bmod(3)) | |
| 241 | f%d2hph(3) = -dBph*f%dBmod(3)/bmod2 + & | ||
| 242 | 178158 | Bph/bmod2*(2.0_dp*f%dBmod(1)*f%dBmod(3)/f%Bmod - f%d2Bmod(3)) | |
| 243 | end if | ||
| 244 | |||
| 245 |
2/2✓ Branch 0 taken 94938 times.
✓ Branch 1 taken 219758 times.
|
314696 | if (mode_secders == 2) then |
| 246 | f%d2hth((/4, 6/)) = Bth/bmod2*( & | ||
| 247 |
2/2✓ Branch 0 taken 189876 times.
✓ Branch 1 taken 94938 times.
|
284814 | 2.0_dp*f%dBmod((/2, 3/))**2/f%Bmod - f%d2Bmod((/4, 6/))) |
| 248 | f%d2hph((/4, 6/)) = Bph/bmod2*( & | ||
| 249 |
2/2✓ Branch 0 taken 189876 times.
✓ Branch 1 taken 94938 times.
|
284814 | 2.0_dp*f%dBmod((/2, 3/))**2/f%Bmod - f%d2Bmod((/4, 6/))) |
| 250 | |||
| 251 | f%d2hth(5) = Bth/bmod2*( & | ||
| 252 | 94938 | 2.0_dp*f%dBmod(2)*f%dBmod(3)/f%Bmod - f%d2Bmod(5)) | |
| 253 | f%d2hph(5) = Bph/bmod2*( & | ||
| 254 | 94938 | 2.0_dp*f%dBmod(2)*f%dBmod(3)/f%Bmod - f%d2Bmod(5)) | |
| 255 | end if | ||
| 256 | 314696 | end subroutine eval_field_booz_device | |
| 257 | |||
| 258 | end module field_can_boozer | ||
| 259 |