| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module magfie_sub | ||
| 2 | use spline_vmec_sub, only: vmec_field | ||
| 3 | use field_can_meiss, only: magfie_meiss, field_noncan | ||
| 4 | use field_can_albert, only: magfie_albert | ||
| 5 | use magfie_can_boozer_sub, only: magfie_can, magfie_boozer | ||
| 6 | use field_base, only: magnetic_field_t | ||
| 7 | use field_splined, only: splined_field_t | ||
| 8 | use field_spectre, only: spectre_field_t | ||
| 9 | #ifdef SIMPLE_ENABLE_VMECPP | ||
| 10 | use field_vmecpp, only: magfie_vmecpp | ||
| 11 | #endif | ||
| 12 | use libneo_coordinates, only: coordinate_system_t, chartmap_coordinate_system_t, & | ||
| 13 | RHO_TOR, RHO_POL, PSI_TOR_NORM, PSI_POL_NORM | ||
| 14 | |||
| 15 | implicit none | ||
| 16 | |||
| 17 | ! Define real(dp) kind parameter | ||
| 18 | integer, parameter :: dp = kind(1.0d0) | ||
| 19 | |||
| 20 | ! SIMPLE runs in Gaussian-CGS; the SPECTRE field returns SI (Tesla, meter). | ||
| 21 | real(dp), parameter :: TESLA_TO_GAUSS = 1.0d4 | ||
| 22 | real(dp), parameter :: M_TO_CM = 1.0d2 | ||
| 23 | real(dp), parameter :: M3_TO_CM3 = 1.0d6 | ||
| 24 | |||
| 25 | abstract interface | ||
| 26 | subroutine magfie_base(x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl) | ||
| 27 | import :: dp | ||
| 28 | ! x(i) - set of 3 curvilinear space coordinates (input) | ||
| 29 | ! bmod - dimensionless magnetic field module: bmod=B/B_ref | ||
| 30 | ! sqrtg - Jacobian of space coordinates (square root of | ||
| 31 | ! metric tensor | ||
| 32 | ! bder - derivatives of logarithm of bmod over space coords | ||
| 33 | ! (covariant vector) | ||
| 34 | ! hcovar - covariant components of the unit vector along | ||
| 35 | ! the magnetic field | ||
| 36 | ! hctrvr - contravariant components of the unit vector along | ||
| 37 | ! the magnetic field | ||
| 38 | ! hcurl - contravariant components of the curl of this vector | ||
| 39 | real(dp), intent(in) :: x(3) | ||
| 40 | real(dp), intent(out) :: bmod,sqrtg | ||
| 41 | real(dp), intent(out) :: bder(3),hcovar(3),hctrvr(3),hcurl(3) | ||
| 42 | end subroutine magfie_base | ||
| 43 | end interface | ||
| 44 | |||
| 45 | procedure(magfie_base), pointer :: magfie => null() | ||
| 46 | |||
| 47 | integer, parameter :: TEST=-1, CANFLUX=0, VMEC=1, BOOZER=2, MEISS=3, ALBERT=4, & | ||
| 48 | REFCOORDS = 5, SPECTRE = 6, VMECPP = 7 | ||
| 49 | |||
| 50 | class(magnetic_field_t), allocatable :: refcoords_field | ||
| 51 | type(spectre_field_t), allocatable :: spectre_field | ||
| 52 | |||
| 53 | contains | ||
| 54 | |||
| 55 | 5 | logical function has_magfie_refcoords_field() | |
| 56 | 5 | has_magfie_refcoords_field = allocated(refcoords_field) | |
| 57 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 time.
|
5 | if (has_magfie_refcoords_field) return |
| 58 | 4 | has_magfie_refcoords_field = allocated(field_noncan) | |
| 59 | 4 | end function has_magfie_refcoords_field | |
| 60 | |||
| 61 | 104 | subroutine set_magfie_refcoords_field(field) | |
| 62 | class(magnetic_field_t), intent(in) :: field | ||
| 63 | |||
| 64 |
3/4✓ Branch 0 taken 1 time.
✓ Branch 1 taken 51 times.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
|
53 | if (allocated(refcoords_field)) deallocate (refcoords_field) |
| 65 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 52 times.
|
52 | allocate (refcoords_field, source=field) |
| 66 | 52 | end subroutine set_magfie_refcoords_field | |
| 67 | |||
| 68 | 49 | subroutine set_magfie_spectre_field(field) | |
| 69 | type(spectre_field_t), intent(in) :: field | ||
| 70 | |||
| 71 |
23/42✓ Branch 0 taken 10 times.
✓ Branch 1 taken 39 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 10 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 10 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 10 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 10 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 10 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 10 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 10 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 10 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 10 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 10 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 10 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 10 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 22 times.
✓ Branch 34 taken 10 times.
✓ Branch 35 taken 22 times.
✗ Branch 36 not taken.
✓ Branch 37 taken 22 times.
✗ Branch 38 not taken.
✓ Branch 39 taken 22 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 22 times.
✗ Branch 42 not taken.
|
71 | if (allocated(spectre_field)) deallocate (spectre_field) |
| 72 |
22/42✗ Branch 0 not taken.
✓ Branch 1 taken 49 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 49 times.
✓ Branch 6 taken 49 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 49 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 49 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 49 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 49 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 49 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 49 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 49 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 49 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 49 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 49 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 49 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 49 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 49 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 104 times.
✓ Branch 36 taken 49 times.
✓ Branch 37 taken 104 times.
✗ Branch 38 not taken.
✓ Branch 39 taken 104 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 104 times.
✗ Branch 42 not taken.
✓ Branch 43 taken 104 times.
✗ Branch 44 not taken.
|
153 | allocate (spectre_field, source=field) |
| 73 | 49 | end subroutine set_magfie_spectre_field | |
| 74 | |||
| 75 | 4294051 | subroutine axis_radial_coordinate(radial_is_rho, ierr) | |
| 76 | ! The RK axis map must follow the active evaluator's integration chart. | ||
| 77 | ! Non-reference evaluators retain their existing flux-coordinate recovery. | ||
| 78 | logical, intent(out) :: radial_is_rho | ||
| 79 | integer, intent(out) :: ierr | ||
| 80 | |||
| 81 | 4294051 | radial_is_rho = .false. | |
| 82 | 4294051 | ierr = 0 | |
| 83 |
2/2✓ Branch 0 taken 1120 times.
✓ Branch 1 taken 4292931 times.
|
4294051 | if (.not. associated(magfie, magfie_refcoords)) return |
| 84 | 1120 | ierr = 2 | |
| 85 |
1/2✓ Branch 0 taken 1120 times.
✗ Branch 1 not taken.
|
1120 | if (allocated(refcoords_field)) then |
| 86 | 1120 | call reference_axis_radial_coordinate(refcoords_field, radial_is_rho, ierr) | |
| 87 | ✗ | else if (allocated(field_noncan)) then | |
| 88 | ✗ | call reference_axis_radial_coordinate(field_noncan, radial_is_rho, ierr) | |
| 89 | end if | ||
| 90 | end subroutine axis_radial_coordinate | ||
| 91 | |||
| 92 | 1120 | subroutine reference_axis_radial_coordinate(field, radial_is_rho, ierr) | |
| 93 | use libneo_coordinates, only: vmec_coordinate_system_t | ||
| 94 | class(magnetic_field_t), intent(in) :: field | ||
| 95 | logical, intent(out) :: radial_is_rho | ||
| 96 | integer, intent(out) :: ierr | ||
| 97 | |||
| 98 | 1120 | radial_is_rho = .false. | |
| 99 | 1120 | ierr = 2 | |
| 100 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1120 times.
|
1120 | if (.not. allocated(field%coords)) return |
| 101 | select type (coords => field%coords) | ||
| 102 | class is (chartmap_coordinate_system_t) | ||
| 103 | 2000 | select case (coords%rho_convention) | |
| 104 | case (RHO_TOR, RHO_POL) | ||
| 105 | 1000 | radial_is_rho = .true. | |
| 106 |
1/2✓ Branch 0 taken 1000 times.
✗ Branch 1 not taken.
|
1000 | ierr = 0 |
| 107 | end select | ||
| 108 | class is (vmec_coordinate_system_t) | ||
| 109 | ! splined_field_t evaluates on its sqrt(s) grid; a raw VMEC field | ||
| 110 | ! evaluates in normalized toroidal flux itself. | ||
| 111 | select type (field) | ||
| 112 | type is (splined_field_t) | ||
| 113 | 119 | radial_is_rho = .true. | |
| 114 | end select | ||
| 115 | 120 | ierr = 0 | |
| 116 | end select | ||
| 117 | 1120 | end subroutine reference_axis_radial_coordinate | |
| 118 | |||
| 119 | 129 | subroutine init_magfie(id) | |
| 120 | integer, intent(in) :: id | ||
| 121 | |||
| 122 | 132 | select case(id) | |
| 123 | case(TEST) | ||
| 124 | 3 | magfie => magfie_test | |
| 125 | case(CANFLUX) | ||
| 126 | 1 | magfie => magfie_can | |
| 127 | case(VMEC) | ||
| 128 | 55 | magfie => magfie_vmec | |
| 129 | case(BOOZER) | ||
| 130 | 38 | magfie => magfie_boozer | |
| 131 | case(MEISS) | ||
| 132 | 2 | magfie => magfie_meiss | |
| 133 | case(ALBERT) | ||
| 134 | ✗ | magfie => magfie_albert | |
| 135 | case(REFCOORDS) | ||
| 136 | 4 | magfie => magfie_refcoords | |
| 137 | case(SPECTRE) | ||
| 138 | 26 | magfie => magfie_spectre | |
| 139 | #ifdef SIMPLE_ENABLE_VMECPP | ||
| 140 | case (VMECPP) | ||
| 141 | magfie => magfie_vmecpp | ||
| 142 | #endif | ||
| 143 | case default | ||
| 144 | ✗ | print *,'init_magfie: unknown id ', id | |
| 145 |
7/9✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 55 times.
✓ Branch 3 taken 38 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 26 times.
✗ Branch 8 not taken.
|
129 | error stop |
| 146 | end select | ||
| 147 | 1120 | end subroutine init_magfie | |
| 148 | |||
| 149 | |||
| 150 | 10347 | subroutine magfie_test(x, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) | |
| 151 | !> Magnetic field for analytic circular tokamak (TEST field). | ||
| 152 | !> Coordinates: x(1)=r (minor radius), x(2)=theta (poloidal), x(3)=phi (toroidal) | ||
| 153 | !> Uses same geometry as field_can_test: B0=1, R0=1, a=0.5, iota=1 | ||
| 154 | !> | ||
| 155 | !> WARNING: hcurl is set to zero (curvature drift not computed). | ||
| 156 | !> This is acceptable for symplectic integration (integmode > 0) which uses | ||
| 157 | !> field_can_test instead. For RK45 integration (integmode=0), curvature | ||
| 158 | !> drift would be missing - use symplectic integration with TEST field. | ||
| 159 | implicit none | ||
| 160 | |||
| 161 | real(dp), intent(in) :: x(3) | ||
| 162 | real(dp), intent(out) :: bmod, sqrtg, bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 163 | |||
| 164 | real(dp), parameter :: B0 = 1.0_dp, R0 = 1.0_dp, a = 0.5_dp, iota0 = 1.0_dp | ||
| 165 | real(dp) :: r, theta, cth, sth, R_cyl, dBmod_dr, dBmod_dth | ||
| 166 | |||
| 167 | 10347 | r = x(1) | |
| 168 | 10347 | theta = x(2) | |
| 169 | 10347 | cth = cos(theta) | |
| 170 | 10347 | sth = sin(theta) | |
| 171 | |||
| 172 | ! Major radius at this point | ||
| 173 | 10347 | R_cyl = R0 + r * cth | |
| 174 | |||
| 175 | ! Magnetic field magnitude: B = B0 * (1 - r/R0 * cos(theta)) | ||
| 176 | 10347 | bmod = B0 * (1.0_dp - r / R0 * cth) | |
| 177 | |||
| 178 | ! Jacobian sqrt(g) = r * R for circular tokamak | ||
| 179 | 10347 | sqrtg = r * R_cyl | |
| 180 | |||
| 181 | ! Derivatives of log(B) | ||
| 182 | 10347 | dBmod_dr = -B0 / R0 * cth | |
| 183 | 10347 | dBmod_dth = B0 * r / R0 * sth | |
| 184 | 10347 | bder(1) = dBmod_dr / bmod | |
| 185 | 10347 | bder(2) = dBmod_dth / bmod | |
| 186 | 10347 | bder(3) = 0.0_dp | |
| 187 | |||
| 188 | ! Covariant components of unit vector h = B/|B| | ||
| 189 | ! In (r, theta, phi) coordinates for circular tokamak with iota=1 | ||
| 190 | 10347 | hcovar(1) = 0.0_dp | |
| 191 | 10347 | hcovar(2) = iota0 * (1.0_dp - r**2 / a**2) * r**2 / R0 / bmod | |
| 192 | 10347 | hcovar(3) = R_cyl / bmod | |
| 193 | |||
| 194 | ! Contravariant components | ||
| 195 | 10347 | hctrvr(1) = 0.0_dp | |
| 196 | 10347 | hctrvr(2) = B0 * iota0 / (r * R_cyl * bmod) | |
| 197 | 10347 | hctrvr(3) = B0 / (r * R_cyl * bmod) | |
| 198 | |||
| 199 | ! Curl of h (simplified - not fully computed for TEST field) | ||
| 200 | 10347 | hcurl(1) = 0.0_dp | |
| 201 | 10347 | hcurl(2) = 0.0_dp | |
| 202 | 10347 | hcurl(3) = 0.0_dp | |
| 203 | |||
| 204 | 10347 | end subroutine magfie_test | |
| 205 | |||
| 206 | |||
| 207 | 9062 | subroutine magfie_refcoords(x, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) | |
| 208 | !> magfie in reference coordinates using analytic spline derivatives. | ||
| 209 | !> | ||
| 210 | !> Input x is in the splined_field_t coordinate system (typically | ||
| 211 | !> (rho, theta, phi) where rho = sqrt(s) for VMEC coordinate systems). | ||
| 212 | ! | ||
| 213 | implicit none | ||
| 214 | |||
| 215 | real(dp), intent(in) :: x(3) | ||
| 216 | real(dp), intent(out) :: bmod, sqrtg | ||
| 217 | real(dp), intent(out) :: bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 218 | |||
| 219 |
1/2✓ Branch 0 taken 9062 times.
✗ Branch 1 not taken.
|
9062 | if (allocated(refcoords_field)) then |
| 220 | select type (refcoords_field) | ||
| 221 | type is (splined_field_t) | ||
| 222 | call magfie_splined_refcoords(refcoords_field, x, bmod, sqrtg, bder, & | ||
| 223 | 9027 | hcovar, hctrvr, hcurl) | |
| 224 | 9027 | return | |
| 225 | class default | ||
| 226 | call magfie_generic_refcoords(refcoords_field, x, bmod, sqrtg, bder, & | ||
| 227 | 35 | hcovar, hctrvr, hcurl) | |
| 228 | 35 | return | |
| 229 | end select | ||
| 230 | end if | ||
| 231 | ✗ | if (allocated(field_noncan)) then | |
| 232 | select type (field_noncan) | ||
| 233 | type is (splined_field_t) | ||
| 234 | call magfie_splined_refcoords(field_noncan, x, bmod, sqrtg, bder, & | ||
| 235 | ✗ | hcovar, hctrvr, hcurl) | |
| 236 | ✗ | return | |
| 237 | class default | ||
| 238 | call magfie_generic_refcoords(field_noncan, x, bmod, sqrtg, bder, & | ||
| 239 | ✗ | hcovar, hctrvr, hcurl) | |
| 240 | ✗ | return | |
| 241 | end select | ||
| 242 | end if | ||
| 243 | ✗ | print *, 'magfie_refcoords: no splined reference field available' | |
| 244 | ✗ | error stop | |
| 245 | end subroutine magfie_refcoords | ||
| 246 | |||
| 247 | 9027 | subroutine magfie_splined_refcoords(field, x, bmod, sqrtg, bder, hcovar, & | |
| 248 | hctrvr, hcurl) | ||
| 249 | type(splined_field_t), intent(in) :: field | ||
| 250 | real(dp), intent(in) :: x(3) | ||
| 251 | real(dp), intent(out) :: bmod, sqrtg | ||
| 252 | real(dp), intent(out) :: bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 253 | real(dp) :: Acov(3), Bmod_local | ||
| 254 | real(dp) :: dAcov(3, 3), dhcov(3, 3), dBmod(3) | ||
| 255 | real(dp) :: g(3, 3), ginv_u(3, 3), ginv_x(3, 3) | ||
| 256 | real(dp) :: u_ref(3), J | ||
| 257 | real(dp) :: e_cov(3, 3), jac_signed | ||
| 258 | real(dp) :: sqrtg_abs | ||
| 259 | |||
| 260 | call field%evaluate_with_der(x, Acov, hcovar, Bmod_local, dAcov, dhcov, & | ||
| 261 | 9027 | dBmod) | |
| 262 | 9027 | bmod = Bmod_local | |
| 263 |
2/2✓ Branch 0 taken 27081 times.
✓ Branch 1 taken 9027 times.
|
36108 | bder = dBmod/max(bmod, 1.0d-30) |
| 264 | |||
| 265 | 9027 | call scaled_to_ref_coords(field%coords, x, u_ref, J) | |
| 266 | 9027 | call field%coords%metric_tensor(u_ref, g, ginv_u, sqrtg_abs) | |
| 267 | 9027 | call field%coords%covariant_basis(u_ref, e_cov) | |
| 268 | |||
| 269 | 9027 | jac_signed = signed_jacobian(e_cov) | |
| 270 | 9027 | sqrtg = jac_signed*J | |
| 271 | |||
| 272 | 9027 | call inverse_metric_scaled(J, ginv_u, ginv_x) | |
| 273 |
4/4✓ Branch 0 taken 27081 times.
✓ Branch 1 taken 9027 times.
✓ Branch 2 taken 81243 times.
✓ Branch 3 taken 27081 times.
|
117351 | hctrvr = matmul(ginv_x, hcovar) |
| 274 | |||
| 275 | 9027 | call compute_hcurl(sqrtg, dhcov, hcurl) | |
| 276 | 9027 | end subroutine magfie_splined_refcoords | |
| 277 | |||
| 278 | 35 | subroutine magfie_generic_refcoords(field, x, bmod, sqrtg, bder, hcovar, & | |
| 279 | hctrvr, hcurl) | ||
| 280 | class(magnetic_field_t), intent(in) :: field | ||
| 281 | real(dp), intent(in) :: x(3) | ||
| 282 | real(dp), intent(out) :: bmod, sqrtg | ||
| 283 | real(dp), intent(out) :: bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 284 | real(dp), parameter :: derivative_step = 1.0e-5_dp | ||
| 285 | real(dp) :: Acov(3), dBmod(3), dhcov(3, 3) | ||
| 286 | real(dp) :: ap(3), hp(3), bp, am(3), hm(3), bm | ||
| 287 | real(dp) :: xp(3), xm(3), denominator | ||
| 288 | real(dp) :: g(3, 3), ginv(3, 3), e_cov(3, 3) | ||
| 289 | real(dp) :: sqrtg_abs | ||
| 290 | integer :: i | ||
| 291 | |||
| 292 | 35 | call field%evaluate(x, Acov, hcovar, bmod) | |
| 293 |
2/2✓ Branch 1 taken 105 times.
✓ Branch 2 taken 35 times.
|
140 | do i = 1, 3 |
| 294 | 105 | xp = x | |
| 295 | 105 | xm = x | |
| 296 | 105 | xp(i) = x(i) + derivative_step | |
| 297 | 105 | xm(i) = x(i) - derivative_step | |
| 298 |
3/4✓ Branch 0 taken 35 times.
✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 35 times.
|
105 | if (i == 1 .and. xm(i) < 0.0_dp) xm(i) = x(i) |
| 299 | 105 | call field%evaluate(xp, ap, hp, bp) | |
| 300 | 105 | call field%evaluate(xm, am, hm, bm) | |
| 301 | 105 | denominator = xp(i) - xm(i) | |
| 302 | 105 | dBmod(i) = (bp - bm)/denominator | |
| 303 |
2/2✓ Branch 2 taken 315 times.
✓ Branch 3 taken 105 times.
|
455 | dhcov(i, :) = (hp - hm)/denominator |
| 304 | end do | ||
| 305 |
2/2✓ Branch 0 taken 105 times.
✓ Branch 1 taken 35 times.
|
140 | bder = dBmod/max(bmod, 1.0d-30) |
| 306 | 35 | call field%coords%metric_tensor(x, g, ginv, sqrtg_abs) | |
| 307 | 35 | call field%coords%covariant_basis(x, e_cov) | |
| 308 | 35 | sqrtg = sign(sqrtg_abs, signed_jacobian(e_cov)) | |
| 309 |
4/4✓ Branch 0 taken 105 times.
✓ Branch 1 taken 35 times.
✓ Branch 2 taken 315 times.
✓ Branch 3 taken 105 times.
|
455 | hctrvr = matmul(ginv, hcovar) |
| 310 | 35 | call compute_hcurl(sqrtg, dhcov, hcurl) | |
| 311 | 35 | end subroutine magfie_generic_refcoords | |
| 312 | |||
| 313 | 13 | subroutine evaluate_magfie_refcoords_field(x, Acov, hcov, Bmod, status) | |
| 314 | real(dp), intent(in) :: x(3) | ||
| 315 | real(dp), intent(out) :: Acov(3), hcov(3), Bmod | ||
| 316 | integer, intent(out) :: status | ||
| 317 | |||
| 318 | 13 | status = 0 | |
| 319 |
1/2✓ Branch 0 taken 13 times.
✗ Branch 1 not taken.
|
13 | if (allocated(refcoords_field)) then |
| 320 | 13 | call refcoords_field%evaluate(x, Acov, hcov, Bmod) | |
| 321 | 13 | return | |
| 322 | end if | ||
| 323 | ✗ | if (allocated(field_noncan)) then | |
| 324 | ✗ | call field_noncan%evaluate(x, Acov, hcov, Bmod) | |
| 325 | ✗ | return | |
| 326 | end if | ||
| 327 | ✗ | Acov = 0.0_dp | |
| 328 | ✗ | hcov = 0.0_dp | |
| 329 | ✗ | Bmod = 0.0_dp | |
| 330 | ✗ | status = 1 | |
| 331 | end subroutine evaluate_magfie_refcoords_field | ||
| 332 | |||
| 333 | |||
| 334 | 2519764 | subroutine magfie_spectre(x, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) | |
| 335 | !> magfie in the SPECTRE stacked-rho chart x = (rho_g, theta, zeta). | ||
| 336 | !> | ||
| 337 | !> The field returns |B|, covariant h = B/|B|, the exact sqrt(g) B^i (the | ||
| 338 | !> metric-free 2-form) and the analytic x-derivatives of |B| and h in SI units; | ||
| 339 | !> the coordinate system provides the analytic metric and its derivative. grad | ||
| 340 | !> log|B| (bder) and curl(h) (hcurl), the small drift corrections, follow | ||
| 341 | !> analytically with no finite differences. | ||
| 342 | implicit none | ||
| 343 | |||
| 344 | real(dp), intent(in) :: x(3) | ||
| 345 | real(dp), intent(out) :: bmod, sqrtg | ||
| 346 | real(dp), intent(out) :: bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 347 | |||
| 348 | real(dp) :: hcov(3), Bmod_si, sqgBctr(3), sqrtg_si | ||
| 349 | real(dp) :: dhcov_si(3, 3), dBmod_si(3), dhcov(3, 3) | ||
| 350 | |||
| 351 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2519764 times.
|
2519764 | if (.not. allocated(spectre_field)) then |
| 352 | ✗ | print *, 'magfie_spectre: spectre_field not set' | |
| 353 | ✗ | print *, 'Call set_magfie_spectre_field before init_magfie(SPECTRE)' | |
| 354 | ✗ | error stop | |
| 355 | end if | ||
| 356 | |||
| 357 | call spectre_field%evaluate_der(x, hcov, Bmod_si, sqgBctr, sqrtg_si, & | ||
| 358 | 2519764 | dhcov_si, dBmod_si) | |
| 359 | |||
| 360 | ! Gaussian-CGS boundary conversion (single place). The SPECTRE field is SI | ||
| 361 | ! (Tesla, meter); SIMPLE integrates in Gaussian-CGS. |B|: Tesla -> Gauss. | ||
| 362 | ! Lengths meter -> cm: the Jacobian sqrt(g) carries L**3, covariant h_i carry | ||
| 363 | ! L, contravariant h^i = B^i/|B| carries 1/L. bder = d log|B|/dx is scale-free | ||
| 364 | ! because x = (rho_g, theta, zeta) is dimensionless, so it needs no factor. | ||
| 365 | 2519764 | bmod = Bmod_si*TESLA_TO_GAUSS | |
| 366 | 2519764 | sqrtg = sqrtg_si*M3_TO_CM3 | |
| 367 |
2/2✓ Branch 0 taken 7559292 times.
✓ Branch 1 taken 2519764 times.
|
10079056 | hcovar = hcov*M_TO_CM |
| 368 |
2/2✓ Branch 0 taken 2519764 times.
✓ Branch 1 taken 7559292 times.
|
10079056 | hctrvr = sqgBctr/(sqrtg_si*Bmod_si)/M_TO_CM |
| 369 |
2/2✓ Branch 0 taken 2519764 times.
✓ Branch 1 taken 7559292 times.
|
10079056 | bder = dBmod_si/Bmod_si |
| 370 | |||
| 371 | ! hcurl^i = eps^ijk d_j h_k / sqrt(g). The covariant derivative d h_i/dx carries | ||
| 372 | ! the same M_TO_CM length factor as hcovar; compute_hcurl divides by the cm | ||
| 373 | ! sqrt(g), so hcurl comes out contravariant in Gaussian-CGS like hctrvr. | ||
| 374 | ! libneo returns dhcov_si(i, j) = d h_i/dx_j; compute_hcurl takes the | ||
| 375 | ! derivative index first. | ||
| 376 |
4/4✓ Branch 0 taken 7559292 times.
✓ Branch 1 taken 2519764 times.
✓ Branch 2 taken 22677876 times.
✓ Branch 3 taken 7559292 times.
|
32756932 | dhcov = transpose(dhcov_si)*M_TO_CM |
| 377 | 2519764 | call compute_hcurl(sqrtg, dhcov, hcurl) | |
| 378 | 2519764 | end subroutine magfie_spectre | |
| 379 | |||
| 380 | |||
| 381 | 9027 | subroutine scaled_to_ref_coords(coords, x_scaled, u_ref, J) | |
| 382 | class(coordinate_system_t), intent(in) :: coords | ||
| 383 | real(dp), intent(in) :: x_scaled(3) | ||
| 384 | real(dp), intent(out) :: u_ref(3) | ||
| 385 | real(dp), intent(out) :: J | ||
| 386 | |||
| 387 | select type (coords) | ||
| 388 | class is (chartmap_coordinate_system_t) | ||
| 389 | ! Chartmap methods (metric_tensor, covariant_basis) expect native rho coords. | ||
| 390 | ! The coordinate splines use rho, while A_phi is a separate s profile. | ||
| 391 | 7915 | u_ref = x_scaled | |
| 392 | 7915 | J = 1.0d0 | |
| 393 | class default | ||
| 394 | ! VMEC coordinate system: integration uses sqrt(s), VMEC methods expect s | ||
| 395 |
2/2✓ Branch 0 taken 3336 times.
✓ Branch 1 taken 1112 times.
|
4448 | u_ref = [x_scaled(1)**2, x_scaled(2), x_scaled(3)] |
| 396 | 1112 | J = 2.0d0*x_scaled(1) | |
| 397 | end select | ||
| 398 | 9027 | end subroutine scaled_to_ref_coords | |
| 399 | |||
| 400 | |||
| 401 | 9027 | subroutine inverse_metric_scaled(J, ginv_u, ginv_x) | |
| 402 | real(dp), intent(in) :: J | ||
| 403 | real(dp), intent(in) :: ginv_u(3, 3) | ||
| 404 | real(dp), intent(out) :: ginv_x(3, 3) | ||
| 405 | |||
| 406 | 9027 | ginv_x = ginv_u | |
| 407 | 9027 | ginv_x(1, 1) = ginv_u(1, 1)/(J*J) | |
| 408 | 9027 | ginv_x(1, 2) = ginv_u(1, 2)/J | |
| 409 | 9027 | ginv_x(1, 3) = ginv_u(1, 3)/J | |
| 410 | 9027 | ginv_x(2, 1) = ginv_u(2, 1)/J | |
| 411 | 9027 | ginv_x(3, 1) = ginv_u(3, 1)/J | |
| 412 | 9027 | end subroutine inverse_metric_scaled | |
| 413 | |||
| 414 | |||
| 415 | 2528826 | subroutine compute_hcurl(sqrtg, dh, hcurl) | |
| 416 | !> Contravariant curl h^i = eps^ijk d_j h_k / sqrt(g) with the derivative | ||
| 417 | !> index first: dh(j, k) = d h_k / d x_j. | ||
| 418 | real(dp), intent(in) :: sqrtg | ||
| 419 | real(dp), intent(in) :: dh(3, 3) | ||
| 420 | real(dp), intent(out) :: hcurl(3) | ||
| 421 | |||
| 422 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2528826 times.
|
2528826 | if (abs(sqrtg) <= 0.0d0) error stop 'compute_hcurl: sqrtg must be nonzero' |
| 423 | |||
| 424 | 2528826 | hcurl(1) = (dh(2, 3) - dh(3, 2))/sqrtg | |
| 425 | 2528826 | hcurl(2) = (dh(3, 1) - dh(1, 3))/sqrtg | |
| 426 | 2528826 | hcurl(3) = (dh(1, 2) - dh(2, 1))/sqrtg | |
| 427 | 2528826 | end subroutine compute_hcurl | |
| 428 | |||
| 429 | |||
| 430 | 9062 | pure function signed_jacobian(e_cov) result(jac) | |
| 431 | real(dp), intent(in) :: e_cov(3, 3) | ||
| 432 | real(dp) :: jac | ||
| 433 | real(dp) :: c(3) | ||
| 434 | |||
| 435 | 9062 | c = cross_product(e_cov(:, 2), e_cov(:, 3)) | |
| 436 |
2/2✓ Branch 0 taken 27186 times.
✓ Branch 1 taken 9062 times.
|
36248 | jac = dot_product(e_cov(:, 1), c) |
| 437 | 9062 | end function signed_jacobian | |
| 438 | |||
| 439 | |||
| 440 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9062 times.
|
9062 | pure function cross_product(a, b) result(c) |
| 441 | real(dp), intent(in) :: a(3), b(3) | ||
| 442 | real(dp) :: c(3) | ||
| 443 | |||
| 444 | 9062 | c(1) = a(2)*b(3) - a(3)*b(2) | |
| 445 | 9062 | c(2) = a(3)*b(1) - a(1)*b(3) | |
| 446 | 9062 | c(3) = a(1)*b(2) - a(2)*b(1) | |
| 447 | 9062 | end function cross_product | |
| 448 | |||
| 449 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 450 | ! | ||
| 451 | 141859767 | subroutine magfie_vmec(x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl) | |
| 452 | ! | ||
| 453 | ! Computes magnetic field module in units of the magnetic code - bmod, | ||
| 454 | ! square root of determinant of the metric tensor - sqrtg, | ||
| 455 | ! derivatives of the logarythm of the magnetic field module | ||
| 456 | ! over coordinates - bder, | ||
| 457 | ! covariant componets of the unit vector of the magnetic | ||
| 458 | ! field direction - hcovar, | ||
| 459 | ! contravariant components of this vector - hctrvr, | ||
| 460 | ! contravariant component of the curl of this vector - hcurl | ||
| 461 | ! Order of coordinates is the following: x(1)=s (normalized toroidal flux), | ||
| 462 | ! x(2)=theta (VMEC poloidal angle), x(3)=varphi (geometrical toroidal angle). | ||
| 463 | ! | ||
| 464 | ! Input parameters: | ||
| 465 | ! formal: x(3) - array of VMEC coordinates | ||
| 466 | ! Output parameters: | ||
| 467 | ! formal: bmod | ||
| 468 | ! sqrtg | ||
| 469 | ! bder(3) - derivatives of $\log(B)$ | ||
| 470 | ! hcovar(3) - covariant components of unit vector $\bh$ along $\bB$ | ||
| 471 | ! hctrvr(3) - contra-variant components of unit vector $\bh$ along $\bB$ | ||
| 472 | ! hcurl(3) - contra-variant components of curl of $\bh$ | ||
| 473 | ! | ||
| 474 | ! Called routines: vmec_field | ||
| 475 | ! | ||
| 476 | implicit none | ||
| 477 | ! | ||
| 478 | real(dp), parameter :: twopi=2.d0*3.14159265358979d0, hs=1.d-3, ht=hs*twopi, hp=ht/5.d0 | ||
| 479 | ! | ||
| 480 | real(dp), intent(out) :: bmod,sqrtg | ||
| 481 | real(dp) :: s,theta,varphi,A_theta,A_phi,dA_theta_ds,dA_phi_ds,aiota, & | ||
| 482 | sqg,alam,dl_ds,dl_dt,dl_dp,Bctrvr_vartheta,Bctrvr_varphi, & | ||
| 483 | Bcovar_r,Bcovar_vartheta,Bcovar_varphi | ||
| 484 | real(dp) :: cjac,bcov_s_vmec,bcov_t_vmec,bcov_p_vmec | ||
| 485 | real(dp) :: dhs_dt,dhs_dp,dht_ds,dht_dp,dhp_ds,dhp_dt | ||
| 486 | real(dp), dimension(3), intent(in) :: x | ||
| 487 | real(dp), dimension(3), intent(out) :: bder,hcovar,hctrvr,hcurl | ||
| 488 | ! | ||
| 489 | ! Begin derivatives over s | ||
| 490 | ! | ||
| 491 | 20265681 | theta=x(2) | |
| 492 | 20265681 | varphi=x(3) | |
| 493 | 20265681 | s=x(1)+hs | |
| 494 | ! | ||
| 495 | call vmec_field(s,theta,varphi,A_theta,A_phi,dA_theta_ds,dA_phi_ds,aiota, & | ||
| 496 | sqg,alam,dl_ds,dl_dt,dl_dp,Bctrvr_vartheta,Bctrvr_varphi, & | ||
| 497 | 20265681 | Bcovar_r,Bcovar_vartheta,Bcovar_varphi) | |
| 498 | ! | ||
| 499 | 20265681 | bmod=sqrt(Bctrvr_vartheta*Bcovar_vartheta+Bctrvr_varphi*Bcovar_varphi) | |
| 500 | 20265681 | bcov_s_vmec=Bcovar_r+dl_ds*Bcovar_vartheta | |
| 501 | 20265681 | bcov_t_vmec=(1.d0+dl_dt)*Bcovar_vartheta | |
| 502 | 20265681 | bcov_p_vmec=Bcovar_varphi+dl_dp*Bcovar_vartheta | |
| 503 | 20265681 | bder(1)=bmod | |
| 504 | 20265681 | dht_ds=bcov_t_vmec/bmod | |
| 505 | 20265681 | dhp_ds=bcov_p_vmec/bmod | |
| 506 | ! | ||
| 507 | 20265681 | s=x(1)-hs | |
| 508 | ! | ||
| 509 | call vmec_field(s,theta,varphi,A_theta,A_phi,dA_theta_ds,dA_phi_ds,aiota, & | ||
| 510 | sqg,alam,dl_ds,dl_dt,dl_dp,Bctrvr_vartheta,Bctrvr_varphi, & | ||
| 511 | 20265681 | Bcovar_r,Bcovar_vartheta,Bcovar_varphi) | |
| 512 | ! | ||
| 513 | 20265681 | bmod=sqrt(Bctrvr_vartheta*Bcovar_vartheta+Bctrvr_varphi*Bcovar_varphi) | |
| 514 | 20265681 | bcov_s_vmec=Bcovar_r+dl_ds*Bcovar_vartheta | |
| 515 | 20265681 | bcov_t_vmec=(1.d0+dl_dt)*Bcovar_vartheta | |
| 516 | 20265681 | bcov_p_vmec=Bcovar_varphi+dl_dp*Bcovar_vartheta | |
| 517 | 20265681 | bder(1)=(bder(1)-bmod)/(2.d0*hs) | |
| 518 | 20265681 | dht_ds=(dht_ds-bcov_t_vmec/bmod)/(2.d0*hs) | |
| 519 | 20265681 | dhp_ds=(dhp_ds-bcov_p_vmec/bmod)/(2.d0*hs) | |
| 520 | ! | ||
| 521 | ! End derivatives over s | ||
| 522 | ! | ||
| 523 | !------------------------- | ||
| 524 | ! | ||
| 525 | ! Begin derivatives over theta | ||
| 526 | ! | ||
| 527 | 20265681 | s=x(1) | |
| 528 | 20265681 | theta=x(2)+ht | |
| 529 | ! | ||
| 530 | call vmec_field(s,theta,varphi,A_theta,A_phi,dA_theta_ds,dA_phi_ds,aiota, & | ||
| 531 | sqg,alam,dl_ds,dl_dt,dl_dp,Bctrvr_vartheta,Bctrvr_varphi, & | ||
| 532 | 20265681 | Bcovar_r,Bcovar_vartheta,Bcovar_varphi) | |
| 533 | ! | ||
| 534 | 20265681 | bmod=sqrt(Bctrvr_vartheta*Bcovar_vartheta+Bctrvr_varphi*Bcovar_varphi) | |
| 535 | 20265681 | bcov_s_vmec=Bcovar_r+dl_ds*Bcovar_vartheta | |
| 536 | 20265681 | bcov_t_vmec=(1.d0+dl_dt)*Bcovar_vartheta | |
| 537 | 20265681 | bcov_p_vmec=Bcovar_varphi+dl_dp*Bcovar_vartheta | |
| 538 | 20265681 | bder(2)=bmod | |
| 539 | 20265681 | dhs_dt=bcov_s_vmec/bmod | |
| 540 | 20265681 | dhp_dt=bcov_p_vmec/bmod | |
| 541 | ! | ||
| 542 | 20265681 | theta=x(2)-ht | |
| 543 | ! | ||
| 544 | call vmec_field(s,theta,varphi,A_theta,A_phi,dA_theta_ds,dA_phi_ds,aiota, & | ||
| 545 | sqg,alam,dl_ds,dl_dt,dl_dp,Bctrvr_vartheta,Bctrvr_varphi, & | ||
| 546 | 20265681 | Bcovar_r,Bcovar_vartheta,Bcovar_varphi) | |
| 547 | ! | ||
| 548 | 20265681 | bmod=sqrt(Bctrvr_vartheta*Bcovar_vartheta+Bctrvr_varphi*Bcovar_varphi) | |
| 549 | 20265681 | bcov_s_vmec=Bcovar_r+dl_ds*Bcovar_vartheta | |
| 550 | 20265681 | bcov_t_vmec=(1.d0+dl_dt)*Bcovar_vartheta | |
| 551 | 20265681 | bcov_p_vmec=Bcovar_varphi+dl_dp*Bcovar_vartheta | |
| 552 | 20265681 | bder(2)=(bder(2)-bmod)/(2.d0*ht) | |
| 553 | 20265681 | dhs_dt=(dhs_dt-bcov_s_vmec/bmod)/(2.d0*ht) | |
| 554 | 20265681 | dhp_dt=(dhp_dt-bcov_p_vmec/bmod)/(2.d0*ht) | |
| 555 | ! | ||
| 556 | ! End derivatives over theta | ||
| 557 | ! | ||
| 558 | !------------------------- | ||
| 559 | ! | ||
| 560 | ! Begin derivatives over varphi | ||
| 561 | ! | ||
| 562 | 20265681 | theta=x(2) | |
| 563 | 20265681 | varphi=x(3)+hp | |
| 564 | ! | ||
| 565 | call vmec_field(s,theta,varphi,A_theta,A_phi,dA_theta_ds,dA_phi_ds,aiota, & | ||
| 566 | sqg,alam,dl_ds,dl_dt,dl_dp,Bctrvr_vartheta,Bctrvr_varphi, & | ||
| 567 | 20265681 | Bcovar_r,Bcovar_vartheta,Bcovar_varphi) | |
| 568 | ! | ||
| 569 | 20265681 | bmod=sqrt(Bctrvr_vartheta*Bcovar_vartheta+Bctrvr_varphi*Bcovar_varphi) | |
| 570 | 20265681 | bcov_s_vmec=Bcovar_r+dl_ds*Bcovar_vartheta | |
| 571 | 20265681 | bcov_t_vmec=(1.d0+dl_dt)*Bcovar_vartheta | |
| 572 | 20265681 | bcov_p_vmec=Bcovar_varphi+dl_dp*Bcovar_vartheta | |
| 573 | 20265681 | bder(3)=bmod | |
| 574 | 20265681 | dhs_dp=bcov_s_vmec/bmod | |
| 575 | 20265681 | dht_dp=bcov_t_vmec/bmod | |
| 576 | ! | ||
| 577 | 20265681 | varphi=x(3)-hp | |
| 578 | ! | ||
| 579 | call vmec_field(s,theta,varphi,A_theta,A_phi,dA_theta_ds,dA_phi_ds,aiota, & | ||
| 580 | sqg,alam,dl_ds,dl_dt,dl_dp,Bctrvr_vartheta,Bctrvr_varphi, & | ||
| 581 | 20265681 | Bcovar_r,Bcovar_vartheta,Bcovar_varphi) | |
| 582 | ! | ||
| 583 | 20265681 | bmod=sqrt(Bctrvr_vartheta*Bcovar_vartheta+Bctrvr_varphi*Bcovar_varphi) | |
| 584 | 20265681 | bcov_s_vmec=Bcovar_r+dl_ds*Bcovar_vartheta | |
| 585 | 20265681 | bcov_t_vmec=(1.d0+dl_dt)*Bcovar_vartheta | |
| 586 | 20265681 | bcov_p_vmec=Bcovar_varphi+dl_dp*Bcovar_vartheta | |
| 587 | 20265681 | bder(3)=(bder(3)-bmod)/(2.d0*hp) | |
| 588 | 20265681 | dhs_dp=(dhs_dp-bcov_s_vmec/bmod)/(2.d0*hp) | |
| 589 | 20265681 | dht_dp=(dht_dp-bcov_t_vmec/bmod)/(2.d0*hp) | |
| 590 | ! | ||
| 591 | ! End derivatives over varphi | ||
| 592 | ! | ||
| 593 | !------------------------- | ||
| 594 | ! | ||
| 595 | 20265681 | varphi=x(3) | |
| 596 | ! | ||
| 597 | call vmec_field(s,theta,varphi,A_theta,A_phi,dA_theta_ds,dA_phi_ds,aiota, & | ||
| 598 | sqg,alam,dl_ds,dl_dt,dl_dp,Bctrvr_vartheta,Bctrvr_varphi, & | ||
| 599 | 20265681 | Bcovar_r,Bcovar_vartheta,Bcovar_varphi) | |
| 600 | ! | ||
| 601 | 20265681 | bmod=sqrt(Bctrvr_vartheta*Bcovar_vartheta+Bctrvr_varphi*Bcovar_varphi) | |
| 602 | 20265681 | cjac=1.d0+dl_dt | |
| 603 | 20265681 | sqrtg=sqg*cjac | |
| 604 |
2/2✓ Branch 0 taken 60797043 times.
✓ Branch 1 taken 20265681 times.
|
81062724 | bder=bder/bmod |
| 605 | 20265681 | bcov_s_vmec=Bcovar_r+dl_ds*Bcovar_vartheta | |
| 606 | 20265681 | bcov_t_vmec=(1.d0+dl_dt)*Bcovar_vartheta | |
| 607 | 20265681 | bcov_p_vmec=Bcovar_varphi+dl_dp*Bcovar_vartheta | |
| 608 | 20265681 | hcovar(1)=bcov_s_vmec/bmod | |
| 609 | 20265681 | hcovar(2)=bcov_t_vmec/bmod | |
| 610 | 20265681 | hcovar(3)=bcov_p_vmec/bmod | |
| 611 | 20265681 | hctrvr(1)=0.d0 | |
| 612 | 20265681 | hctrvr(2)=(Bctrvr_vartheta-dl_dp*Bctrvr_varphi)/(cjac*bmod) | |
| 613 | 20265681 | hctrvr(3)=Bctrvr_varphi/bmod | |
| 614 | 20265681 | hcurl(1)=(dhp_dt-dht_dp)/sqrtg | |
| 615 | 20265681 | hcurl(2)=(dhs_dp-dhp_ds)/sqrtg | |
| 616 | 20265681 | hcurl(3)=(dht_ds-dhs_dt)/sqrtg | |
| 617 | ! | ||
| 618 | 20265681 | end subroutine magfie_vmec | |
| 619 | ! | ||
| 620 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 621 | ! | ||
| 622 | |||
| 623 |
9/12✓ Branch 1 taken 1000 times.
✓ Branch 2 taken 120 times.
✓ Branch 4 taken 120 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 119 times.
✓ Branch 7 taken 1 time.
✓ Branch 8 taken 9027 times.
✓ Branch 9 taken 35 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 13 taken 7915 times.
✓ Branch 14 taken 1112 times.
|
19209 | end module magfie_sub |
| 624 |