| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module field_can_flux | ||
| 2 | |||
| 3 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 4 | use field_can_base, only: count_field_evaluation, field_can_t, twopi | ||
| 5 | |||
| 6 | implicit none | ||
| 7 | private | ||
| 8 | |||
| 9 | public :: evaluate_flux, integ_to_ref_flux, ref_to_integ_flux, eval_field_can | ||
| 10 | public :: twopi | ||
| 11 | |||
| 12 | contains | ||
| 13 | |||
| 14 | 2814928 | subroutine evaluate_flux(f, r, th_c, ph_c, mode_secders) | |
| 15 | type(field_can_t), intent(inout) :: f | ||
| 16 | real(dp), intent(in) :: r, th_c, ph_c | ||
| 17 | integer, intent(in) :: mode_secders | ||
| 18 | |||
| 19 | 1407464 | call eval_field_can(f, r, th_c, ph_c, mode_secders) | |
| 20 | |||
| 21 | 1407464 | call count_field_evaluation(mode_secders) | |
| 22 | 1407464 | end subroutine evaluate_flux | |
| 23 | |||
| 24 | 33 | subroutine integ_to_ref_flux(xinteg, xref) | |
| 25 | use get_can_sub, only: can_to_vmec | ||
| 26 | real(dp), intent(in) :: xinteg(3) | ||
| 27 | real(dp), intent(out) :: xref(3) | ||
| 28 | |||
| 29 | 33 | xref(1) = xinteg(1) | |
| 30 | 33 | call can_to_vmec(xinteg(1), xinteg(2), xinteg(3), xref(2), xref(3)) | |
| 31 | 33 | xref(2) = mod(xref(2), twopi) | |
| 32 | 33 | xref(3) = mod(xref(3), twopi) | |
| 33 | 33 | end subroutine integ_to_ref_flux | |
| 34 | |||
| 35 | 40 | subroutine ref_to_integ_flux(xref, xinteg) | |
| 36 | 33 | use get_can_sub, only: vmec_to_can | |
| 37 | real(dp), intent(in) :: xref(3) | ||
| 38 | real(dp), intent(out) :: xinteg(3) | ||
| 39 | |||
| 40 | 40 | xinteg(1) = xref(1) | |
| 41 | 40 | call vmec_to_can(xref(1), xref(2), xref(3), xinteg(2), xinteg(3)) | |
| 42 | 40 | xinteg(2) = mod(xinteg(2), twopi) | |
| 43 | 40 | xinteg(3) = mod(xinteg(3), twopi) | |
| 44 | 40 | end subroutine ref_to_integ_flux | |
| 45 | |||
| 46 | !> Evaluate magnetic field in canonical coordinates (r, th_c, ph_c) | ||
| 47 | !> Works for A_th linear in r (toroidal flux as radial variable) | ||
| 48 | !> mode_secders = 0: no second derivatives | ||
| 49 | !> mode_secders = 1: second derivatives only in d/dr^2 | ||
| 50 | !> mode_secders = 2: all second derivatives, including mixed | ||
| 51 | 1407464 | subroutine eval_field_can(f, r, th_c, ph_c, mode_secders) | |
| 52 | 40 | use get_can_sub, only: splint_can_coord | |
| 53 | |||
| 54 | type(field_can_t), intent(inout) :: f | ||
| 55 | real(dp), intent(in) :: r, th_c, ph_c | ||
| 56 | integer, intent(in) :: mode_secders | ||
| 57 | |||
| 58 | real(dp) :: Bctr_vartheta, Bctr_varphi, bmod2, sqg, dsqg(3), & | ||
| 59 | d2sqg(6), d3Aphdr3, dummy, Bth, Bph, dBth(3), dBph(3), d2Bth(6), & | ||
| 60 | d2Bph(6), twobmod, dbmod2(3) | ||
| 61 | |||
| 62 | ! Initialize to zero - no angular derivatives for straight field line Ath(r) Aph(r) | ||
| 63 |
2/2✓ Branch 0 taken 4222392 times.
✓ Branch 1 taken 1407464 times.
|
5629856 | f%dAth = 0.0_dp |
| 64 |
2/2✓ Branch 0 taken 4222392 times.
✓ Branch 1 taken 1407464 times.
|
5629856 | f%dAph = 0.0_dp |
| 65 | |||
| 66 | ! Initialize all 2nd derivatives to zero; mode_secders controls which are computed | ||
| 67 |
2/2✓ Branch 0 taken 8444784 times.
✓ Branch 1 taken 1407464 times.
|
9852248 | f%d2Ath = 0.0_dp |
| 68 |
2/2✓ Branch 0 taken 8444784 times.
✓ Branch 1 taken 1407464 times.
|
9852248 | f%d2Aph = 0.0_dp |
| 69 |
2/2✓ Branch 0 taken 8444784 times.
✓ Branch 1 taken 1407464 times.
|
9852248 | f%d2hth = 0.0_dp |
| 70 |
2/2✓ Branch 0 taken 8444784 times.
✓ Branch 1 taken 1407464 times.
|
9852248 | f%d2hph = 0.0_dp |
| 71 |
2/2✓ Branch 0 taken 8444784 times.
✓ Branch 1 taken 1407464 times.
|
9852248 | f%d2Bmod = 0.0_dp |
| 72 | |||
| 73 | call splint_can_coord(.false., mode_secders, r, th_c, ph_c, & | ||
| 74 | f%Ath, f%Aph, f%dAth(1), f%dAph(1), f%d2Aph(1), & | ||
| 75 | d3Aphdr3, & | ||
| 76 | sqg, dsqg(1), dsqg(2), dsqg(3), & | ||
| 77 | Bth, dBth(1), dBth(2), dBth(3), & | ||
| 78 | Bph, dBph(1), dBph(2), dBph(3), & | ||
| 79 | d2sqg(1), d2sqg(2), d2sqg(3), d2sqg(4), d2sqg(5), & | ||
| 80 | d2sqg(6), & | ||
| 81 | d2Bth(1), d2Bth(2), d2Bth(3), d2Bth(4), d2Bth(5), & | ||
| 82 | d2Bth(6), & | ||
| 83 | d2Bph(1), d2Bph(2), d2Bph(3), d2Bph(4), d2Bph(5), & | ||
| 84 | 1407464 | d2Bph(6), dummy) | |
| 85 | |||
| 86 | 1407464 | Bctr_vartheta = -f%dAph(1)/sqg | |
| 87 | 1407464 | Bctr_varphi = f%dAth(1)/sqg | |
| 88 | |||
| 89 | 1407464 | bmod2 = Bctr_vartheta*Bth + Bctr_varphi*Bph | |
| 90 | 1407464 | f%Bmod = sqrt(abs(bmod2)) | |
| 91 | 1407464 | twobmod = 2.0_dp*f%Bmod | |
| 92 | |||
| 93 | dbmod2(1) = (f%dAth(1)*dBph(1) - f%dAph(1)*dBth(1) & | ||
| 94 | 1407464 | - f%d2Aph(1)*Bth - bmod2*dsqg(1))/sqg | |
| 95 | 1407464 | dbmod2(2) = (f%dAth(1)*dBph(2) - f%dAph(1)*dBth(2) - bmod2*dsqg(2))/sqg | |
| 96 | 1407464 | dbmod2(3) = (f%dAth(1)*dBph(3) - f%dAph(1)*dBth(3) - bmod2*dsqg(3))/sqg | |
| 97 | |||
| 98 |
2/2✓ Branch 0 taken 4222392 times.
✓ Branch 1 taken 1407464 times.
|
5629856 | f%dBmod = dbmod2/twobmod |
| 99 | |||
| 100 | 1407464 | f%hth = Bth/f%Bmod | |
| 101 | 1407464 | f%hph = Bph/f%Bmod | |
| 102 |
2/2✓ Branch 0 taken 4222392 times.
✓ Branch 1 taken 1407464 times.
|
5629856 | f%dhth = dBth/f%Bmod - Bth*f%dBmod/bmod2 |
| 103 |
2/2✓ Branch 0 taken 4222392 times.
✓ Branch 1 taken 1407464 times.
|
5629856 | f%dhph = dBph/f%Bmod - Bph*f%dBmod/bmod2 |
| 104 | |||
| 105 |
2/2✓ Branch 0 taken 1089188 times.
✓ Branch 1 taken 318276 times.
|
1407464 | if (mode_secders > 0) then |
| 106 | ! d2dr2 | ||
| 107 | f%d2Bmod(1) = (d2Bph(1)*f%dAth(1) - d2Bth(1)*f%dAph(1) & | ||
| 108 | - 2.0_dp*dBth(1)*f%d2Aph(1) - f%Bmod*f%hth*d3Aphdr3 & | ||
| 109 | 1089188 | - 2.0_dp*dsqg(1)*dbmod2(1) - bmod2*d2sqg(1))/sqg | |
| 110 | |||
| 111 | 1089188 | f%d2Bmod(1) = f%d2Bmod(1)/twobmod - f%dBmod(1)**2/f%Bmod | |
| 112 | |||
| 113 | f%d2hth(1) = d2Bth(1)/f%Bmod - 2.0_dp*dBth(1)*f%dBmod(1)/bmod2 & | ||
| 114 | 1089188 | + Bth/bmod2*(2.0_dp*f%dBmod(1)**2/f%Bmod - f%d2Bmod(1)) | |
| 115 | f%d2hph(1) = d2Bph(1)/f%Bmod - 2.0_dp*dBph(1)*f%dBmod(1)/bmod2 & | ||
| 116 | 1089188 | + Bph/bmod2*(2.0_dp*f%dBmod(1)**2/f%Bmod - f%d2Bmod(1)) | |
| 117 | end if | ||
| 118 | |||
| 119 |
2/2✓ Branch 0 taken 1089188 times.
✓ Branch 1 taken 318276 times.
|
1407464 | if (mode_secders == 2) then |
| 120 | ! d2dth2, d2dph2 | ||
| 121 | f%d2Bmod(4) = (d2Bph(4)*f%dAth(1) - d2Bth(4)*f%dAph(1) & | ||
| 122 | 1089188 | - 2.0_dp*dsqg(2)*dbmod2(2) - bmod2*d2sqg(4))/sqg | |
| 123 | f%d2Bmod(6) = (d2Bph(6)*f%dAth(1) - d2Bth(6)*f%dAph(1) & | ||
| 124 | 1089188 | - 2.0_dp*dsqg(3)*dbmod2(3) - bmod2*d2sqg(6))/sqg | |
| 125 | |||
| 126 | 1089188 | f%d2Bmod(4) = f%d2Bmod(4)/twobmod - f%dBmod(2)**2/f%Bmod | |
| 127 | 1089188 | f%d2Bmod(6) = f%d2Bmod(6)/twobmod - f%dBmod(3)**2/f%Bmod | |
| 128 | |||
| 129 | f%d2hth((/4, 6/)) = d2Bth((/4, 6/))/f%Bmod & | ||
| 130 | - 2.0_dp*dBth((/2, 3/))*f%dBmod((/2, 3/))/bmod2 & | ||
| 131 | + Bth/bmod2*(2.0_dp*f%dBmod((/2, 3/))**2/f%Bmod - & | ||
| 132 |
2/2✓ Branch 0 taken 2178376 times.
✓ Branch 1 taken 1089188 times.
|
3267564 | f%d2Bmod((/4, 6/))) |
| 133 | f%d2hph((/4, 6/)) = d2Bph((/4, 6/))/f%Bmod & | ||
| 134 | - 2.0_dp*dBph((/2, 3/))*f%dBmod((/2, 3/))/bmod2 & | ||
| 135 | + Bph/bmod2*(2.0_dp*f%dBmod((/2, 3/))**2/f%Bmod - & | ||
| 136 |
2/2✓ Branch 0 taken 2178376 times.
✓ Branch 1 taken 1089188 times.
|
3267564 | f%d2Bmod((/4, 6/))) |
| 137 | end if | ||
| 138 | |||
| 139 |
2/2✓ Branch 0 taken 1089188 times.
✓ Branch 1 taken 318276 times.
|
1407464 | if (mode_secders == 2 .or. mode_secders == 3) then |
| 140 | ! d2drdth, d2drdph | ||
| 141 | f%d2Bmod(2) = (d2Bph(2)*f%dAth(1) - d2Bth(2)*f%dAph(1) & | ||
| 142 | - dBth(2)*f%d2Aph(1) - dsqg(1)*dbmod2(2) & | ||
| 143 | 1089188 | - dsqg(2)*dbmod2(1) - bmod2*d2sqg(2))/sqg | |
| 144 | f%d2Bmod(3) = (d2Bph(3)*f%dAth(1) - d2Bth(3)*f%dAph(1) & | ||
| 145 | - dBth(3)*f%d2Aph(1) - dsqg(1)*dbmod2(3) & | ||
| 146 | 1089188 | - dsqg(3)*dbmod2(1) - bmod2*d2sqg(3))/sqg | |
| 147 | |||
| 148 | 1089188 | f%d2Bmod(2) = f%d2Bmod(2)/twobmod - f%dBmod(1)*f%dBmod(2)/f%Bmod | |
| 149 | 1089188 | f%d2Bmod(3) = f%d2Bmod(3)/twobmod - f%dBmod(1)*f%dBmod(3)/f%Bmod | |
| 150 | |||
| 151 | f%d2hth(2) = d2Bth(2)/f%Bmod & | ||
| 152 | - (dBth(1)*f%dBmod(2) + dBth(2)*f%dBmod(1))/bmod2 & | ||
| 153 | 1089188 | + Bth/bmod2*(2.0_dp*f%dBmod(1)*f%dBmod(2)/f%Bmod - f%d2Bmod(2)) | |
| 154 | f%d2hph(2) = d2Bph(2)/f%Bmod & | ||
| 155 | - (dBph(1)*f%dBmod(2) + dBph(2)*f%dBmod(1))/bmod2 & | ||
| 156 | 1089188 | + Bph/bmod2*(2.0_dp*f%dBmod(1)*f%dBmod(2)/f%Bmod - f%d2Bmod(2)) | |
| 157 | |||
| 158 | f%d2hth(3) = d2Bth(3)/f%Bmod & | ||
| 159 | - (dBth(1)*f%dBmod(3) + dBth(3)*f%dBmod(1))/bmod2 & | ||
| 160 | 1089188 | + Bth/bmod2*(2.0_dp*f%dBmod(1)*f%dBmod(3)/f%Bmod - f%d2Bmod(3)) | |
| 161 | f%d2hph(3) = d2Bph(3)/f%Bmod & | ||
| 162 | - (dBph(1)*f%dBmod(3) + dBph(3)*f%dBmod(1))/bmod2 & | ||
| 163 | 1089188 | + Bph/bmod2*(2.0_dp*f%dBmod(1)*f%dBmod(3)/f%Bmod - f%d2Bmod(3)) | |
| 164 | |||
| 165 |
1/2✓ Branch 0 taken 1089188 times.
✗ Branch 1 not taken.
|
1089188 | if (mode_secders == 2) then |
| 166 | ! d2dthdph | ||
| 167 | f%d2Bmod(5) = (d2Bph(5)*f%dAth(1) - d2Bth(5)*f%dAph(1) & | ||
| 168 | - dsqg(2)*dbmod2(3) - dsqg(3)*dbmod2(2) - & | ||
| 169 | 1089188 | bmod2*d2sqg(5))/sqg | |
| 170 | 1089188 | f%d2Bmod(5) = f%d2Bmod(5)/twobmod - f%dBmod(2)*f%dBmod(3)/f%Bmod | |
| 171 | |||
| 172 | f%d2hth(5) = d2Bth(5)/f%Bmod & | ||
| 173 | - (dBth(2)*f%dBmod(3) + dBth(3)*f%dBmod(2))/bmod2 & | ||
| 174 | + Bth/bmod2*(2.0_dp*f%dBmod(2)*f%dBmod(3)/f%Bmod - & | ||
| 175 | 1089188 | f%d2Bmod(5)) | |
| 176 | f%d2hph(5) = d2Bph(5)/f%Bmod & | ||
| 177 | - (dBph(2)*f%dBmod(3) + dBph(3)*f%dBmod(2))/bmod2 & | ||
| 178 | + Bph/bmod2*(2.0_dp*f%dBmod(2)*f%dBmod(3)/f%Bmod - & | ||
| 179 | 1089188 | f%d2Bmod(5)) | |
| 180 | end if | ||
| 181 | end if | ||
| 182 | 1407464 | end subroutine eval_field_can | |
| 183 | |||
| 184 | end module field_can_flux | ||
| 185 |