| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module spectre_fo_field | ||
| 2 | !> Per-volume splines for the SPECTRE full orbit in logical coordinates. | ||
| 3 | !> | ||
| 4 | !> Volume k owns the slab rho in [k-1, k] of the stacked SPECTRE chart. On | ||
| 5 | !> that slab one quintic batch spline holds the cylindrical embedding R, Z | ||
| 6 | !> (cm) and the covariant vector potential A_theta, A_zeta (G cm^2) in | ||
| 7 | !> SPECTRE's A_rho = 0 gauge of that volume; zeta is the cylindrical angle, | ||
| 8 | !> periodic over one field period. The basis, the Jacobian and B from dA | ||
| 9 | !> follow from the jets, so the splined geometry and field are consistent | ||
| 10 | !> with each other. A regular (xi, eta, zeta) axis chart carries R, Z and | ||
| 11 | !> all three components of A; its jets join the polar representation before | ||
| 12 | !> B is evaluated. Each slab is sampled up to its own edges (the upper edge | ||
| 13 | !> just inside it), so the sheet at rho = k has the two one-sided limits of | ||
| 14 | !> volumes k and k+1. Outside its slab a spline continues its edge cell's | ||
| 15 | !> polynomial smoothly; the full orbit uses this only for trial stages. | ||
| 16 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 17 | use interpolate, only: BatchSplineData3D, construct_batch_splines_3d, & | ||
| 18 | destroy_batch_splines_3d, evaluate_batch_splines_3d_der | ||
| 19 | use magfie_sub, only: spectre_field, TESLA_TO_GAUSS, M_TO_CM | ||
| 20 | implicit none | ||
| 21 | private | ||
| 22 | |||
| 23 | integer, parameter :: NQ = 4, AXIS_NQ = 5, SPLINE_ORDER = 5 | ||
| 24 | real(dp), parameter, public :: FO_AXIS_RADIUS = 0.05_dp | ||
| 25 | real(dp), parameter :: AXIS_OUTER = 2*FO_AXIS_RADIUS | ||
| 26 | real(dp), parameter :: A_SCALE = TESLA_TO_GAUSS*M_TO_CM**2 | ||
| 27 | !> Upper-edge sample offset: rho = k itself belongs to volume k+1. | ||
| 28 | real(dp), parameter :: EDGE_INSET = 1.0e-12_dp | ||
| 29 | real(dp), parameter :: TWOPI = 8.0_dp*atan(1.0_dp) | ||
| 30 | |||
| 31 | type(BatchSplineData3D), allocatable :: vol_spl(:) | ||
| 32 | type(BatchSplineData3D) :: axis_spl | ||
| 33 | logical :: axis_built = .false. | ||
| 34 | integer :: fo_field_mvol = 0 | ||
| 35 | logical, public, protected :: fo_field_axisymmetric = .false. | ||
| 36 | |||
| 37 | public :: spectre_fo_field_build, spectre_fo_field_clear, & | ||
| 38 | spectre_fo_field_local, spectre_fo_field_locate, & | ||
| 39 | spectre_fo_local_t | ||
| 40 | |||
| 41 | !> Local geometry and field at one logical point of one volume's spline. | ||
| 42 | type :: spectre_fo_local_t | ||
| 43 | real(dp) :: x(3) = 0.0_dp !> Cartesian position (cm) | ||
| 44 | real(dp) :: e(3, 3) = 0.0_dp !> covariant basis e_i = dX/du^i | ||
| 45 | real(dp) :: dual(3, 3) = 0.0_dp !> contravariant basis (1/cm) | ||
| 46 | real(dp) :: jac = 0.0_dp !> signed Jacobian e_1.(e_2 x e_3) | ||
| 47 | real(dp) :: acov(3) = 0.0_dp !> covariant A (G cm^2), A_rho = 0 | ||
| 48 | real(dp) :: b(3) = 0.0_dp !> unit vector of the interpolated B | ||
| 49 | real(dp) :: bmod = 0.0_dp !> |B| (G) | ||
| 50 | end type spectre_fo_local_t | ||
| 51 | |||
| 52 | contains | ||
| 53 | |||
| 54 | 37 | subroutine spectre_fo_field_build(n_rho, n_theta, n_zeta) | |
| 55 | !> Spline the field that magfie_sub currently holds, every volume. | ||
| 56 | integer, intent(in) :: n_rho, n_theta, n_zeta | ||
| 57 | 37 | real(dp), allocatable :: y(:, :, :, :) | |
| 58 | real(dp) :: u(3), xcyl(3), acov(3), hcov(3), bmod, period, h(3) | ||
| 59 | integer :: k, i, j, l, nfp | ||
| 60 | |||
| 61 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 37 times.
|
37 | if (.not. allocated(spectre_field)) & |
| 62 | ✗ | error stop 'spectre_fo_field_build: no SPECTRE field' | |
| 63 |
3/6✓ Branch 0 taken 37 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 37 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 37 times.
|
37 | if (n_rho < SPLINE_ORDER + 1 .or. n_theta < SPLINE_ORDER + 2 .or. & |
| 64 | n_zeta < SPLINE_ORDER + 2) & | ||
| 65 | ✗ | error stop 'spectre_fo_field_build: grid too coarse for quintic splines' | |
| 66 | 37 | call spectre_fo_field_clear | |
| 67 | 37 | fo_field_mvol = spectre_field%data%Mvol | |
| 68 |
4/4✓ Branch 0 taken 137 times.
✓ Branch 1 taken 21 times.
✓ Branch 2 taken 121 times.
✓ Branch 3 taken 16 times.
|
158 | fo_field_axisymmetric = all(spectre_field%data%in == 0) |
| 69 | 37 | nfp = max(spectre_field%data%Nfp, 1) | |
| 70 | 37 | period = TWOPI/real(nfp, dp) | |
| 71 |
23/30✗ Branch 0 not taken.
✓ Branch 1 taken 37 times.
✓ Branch 2 taken 37 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 37 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 37 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 37 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 37 times.
✓ Branch 17 taken 37 times.
✓ Branch 18 taken 111 times.
✓ Branch 19 taken 37 times.
✓ Branch 20 taken 111 times.
✓ Branch 21 taken 37 times.
✓ Branch 22 taken 111 times.
✓ Branch 23 taken 37 times.
✓ Branch 24 taken 111 times.
✓ Branch 25 taken 37 times.
✓ Branch 26 taken 111 times.
✓ Branch 27 taken 37 times.
✓ Branch 28 taken 111 times.
✓ Branch 29 taken 111 times.
✓ Branch 30 taken 37 times.
✓ Branch 31 taken 73 times.
✓ Branch 32 taken 37 times.
|
1183 | allocate (vol_spl(fo_field_mvol)) |
| 72 |
11/22✗ Branch 0 not taken.
✓ Branch 1 taken 37 times.
✓ Branch 2 taken 37 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 37 times.
✓ Branch 6 taken 37 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 37 times.
✓ Branch 10 taken 37 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 37 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 37 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 37 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 37 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 37 times.
|
222 | allocate (y(n_rho, n_theta, n_zeta, NQ)) |
| 73 | h = [1.0_dp/real(n_rho - 1, dp), TWOPI/real(n_theta - 1, dp), & | ||
| 74 |
2/2✓ Branch 0 taken 111 times.
✓ Branch 1 taken 37 times.
|
148 | period/real(n_zeta - 1, dp)] |
| 75 |
2/2✓ Branch 0 taken 73 times.
✓ Branch 1 taken 37 times.
|
110 | do k = 1, fo_field_mvol |
| 76 | 73 | !$omp parallel do collapse(2) private(i, j, l, u, xcyl, acov, hcov, bmod) | |
| 77 | do l = 1, n_zeta | ||
| 78 | do j = 1, n_theta | ||
| 79 | do i = 1, n_rho | ||
| 80 | u = [real(k - 1, dp) + real(i - 1, dp)*h(1), & | ||
| 81 | real(j - 1, dp)*h(2), real(l - 1, dp)*h(3)] | ||
| 82 | ! Shared sheet geometry is sampled at the sheet itself; | ||
| 83 | ! only A needs a one-sided sample in the owning gauge. | ||
| 84 | call spectre_field%coords%evaluate_cyl(u, xcyl) | ||
| 85 | if (i == n_rho) u(1) = real(k, dp) - EDGE_INSET | ||
| 86 | if (i == 1 .and. k == 1) u(1) = EDGE_INSET | ||
| 87 | call spectre_field%evaluate(u, acov, hcov, bmod) | ||
| 88 | y(i, j, l, :) = [xcyl(1)*M_TO_CM, xcyl(3)*M_TO_CM, & | ||
| 89 | acov(2)*A_SCALE, acov(3)*A_SCALE] | ||
| 90 | end do | ||
| 91 | end do | ||
| 92 | end do | ||
| 93 | !$omp end parallel do | ||
| 94 | call construct_batch_splines_3d([real(k - 1, dp), 0.0_dp, 0.0_dp], & | ||
| 95 | [real(k, dp), TWOPI, period], y, & | ||
| 96 | [SPLINE_ORDER, SPLINE_ORDER, SPLINE_ORDER], & | ||
| 97 |
4/4✓ Branch 0 taken 219 times.
✓ Branch 1 taken 73 times.
✓ Branch 2 taken 219 times.
✓ Branch 3 taken 73 times.
|
548 | [.false., .true., .true.], vol_spl(k)) |
| 98 | end do | ||
| 99 | 37 | call build_axis_spline(n_rho, n_zeta, period) | |
| 100 |
1/2✓ Branch 0 taken 37 times.
✗ Branch 1 not taken.
|
37 | end subroutine spectre_fo_field_build |
| 101 | |||
| 102 | 37 | subroutine build_axis_spline(n_rho, n_zeta, period) | |
| 103 | ! Regular first-volume chart xi=rho*cos(theta), eta=rho*sin(theta). | ||
| 104 | ! Tabulate all three components of the potential in that chart. | ||
| 105 | integer, intent(in) :: n_rho, n_zeta | ||
| 106 | real(dp), intent(in) :: period | ||
| 107 | 37 | real(dp), allocatable :: values(:, :, :, :) | |
| 108 | real(dp) :: xi, eta, rho, u(3), cyl(3), a(3), hc(3), bmod, h | ||
| 109 | integer :: n, i, j, k | ||
| 110 | |||
| 111 | 37 | n = 2*max(3, n_rho/4) + 1 | |
| 112 | 37 | h = 2*AXIS_OUTER/real(n - 1, dp) | |
| 113 |
11/22✗ Branch 0 not taken.
✓ Branch 1 taken 37 times.
✓ Branch 2 taken 37 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 37 times.
✓ Branch 6 taken 37 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 37 times.
✓ Branch 10 taken 37 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 37 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 37 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 37 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 37 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 37 times.
|
222 | allocate (values(n, n, n_zeta, AXIS_NQ)) |
| 114 | !$omp parallel do collapse(2) & | ||
| 115 | 37 | !$omp private(i, j, k, xi, eta, rho, u, cyl, a, hc, bmod) | |
| 116 | do k = 1, n_zeta | ||
| 117 | do j = 1, n | ||
| 118 | eta = -AXIS_OUTER + (j - 1)*h | ||
| 119 | do i = 1, n | ||
| 120 | xi = -AXIS_OUTER + (i - 1)*h | ||
| 121 | rho = sqrt(xi*xi + eta*eta) | ||
| 122 | u = [rho, atan2(eta, xi), (k - 1)*period/(n_zeta - 1)] | ||
| 123 | call spectre_field%coords%evaluate_cyl(u, cyl) | ||
| 124 | if (rho < EDGE_INSET) u(1:2) = [EDGE_INSET, 0.0_dp] | ||
| 125 | call spectre_field%evaluate(u, a, hc, bmod) | ||
| 126 | values(i, j, k, 1:2) = [cyl(1), cyl(3)]*M_TO_CM | ||
| 127 | if (rho < EDGE_INSET) then | ||
| 128 | values(i, j, k, 3:4) = 0.0_dp | ||
| 129 | else | ||
| 130 | values(i, j, k, 3:4) = [-eta, xi]*a(2)*A_SCALE/rho**2 | ||
| 131 | end if | ||
| 132 | values(i, j, k, 5) = a(3)*A_SCALE | ||
| 133 | end do | ||
| 134 | end do | ||
| 135 | end do | ||
| 136 | !$omp end parallel do | ||
| 137 | call construct_batch_splines_3d([-AXIS_OUTER, -AXIS_OUTER, 0.0_dp], & | ||
| 138 | [AXIS_OUTER, AXIS_OUTER, period], values, & | ||
| 139 | [SPLINE_ORDER, SPLINE_ORDER, SPLINE_ORDER], & | ||
| 140 |
2/2✓ Branch 0 taken 111 times.
✓ Branch 1 taken 37 times.
|
148 | [.false., .false., .true.], axis_spl) |
| 141 | 37 | axis_built = .true. | |
| 142 |
1/2✓ Branch 0 taken 37 times.
✗ Branch 1 not taken.
|
37 | end subroutine build_axis_spline |
| 143 | |||
| 144 | 129 | subroutine spectre_fo_field_clear | |
| 145 | integer :: k | ||
| 146 | |||
| 147 | 129 | fo_field_axisymmetric = .false. | |
| 148 | 129 | fo_field_mvol = 0 | |
| 149 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 121 times.
|
129 | if (axis_built) call destroy_batch_splines_3d(axis_spl) |
| 150 | 129 | axis_built = .false. | |
| 151 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 121 times.
|
129 | if (.not. allocated(vol_spl)) return |
| 152 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 8 times.
|
28 | do k = 1, size(vol_spl) |
| 153 | 28 | call destroy_batch_splines_3d(vol_spl(k)) | |
| 154 | end do | ||
| 155 |
5/8✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 20 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 8 times.
|
28 | deallocate (vol_spl) |
| 156 | end subroutine spectre_fo_field_clear | ||
| 157 | |||
| 158 |
14/14✓ Branch 0 taken 415965 times.
✓ Branch 1 taken 138655 times.
✓ Branch 2 taken 415965 times.
✓ Branch 3 taken 138655 times.
✓ Branch 4 taken 1247895 times.
✓ Branch 5 taken 415965 times.
✓ Branch 6 taken 415965 times.
✓ Branch 7 taken 138655 times.
✓ Branch 8 taken 1247895 times.
✓ Branch 9 taken 415965 times.
✓ Branch 10 taken 415965 times.
✓ Branch 11 taken 138655 times.
✓ Branch 12 taken 415965 times.
✓ Branch 13 taken 138655 times.
|
4714270 | subroutine spectre_fo_field_local(u, owner, loc) |
| 159 | real(dp), intent(in) :: u(3) | ||
| 160 | integer, intent(in) :: owner | ||
| 161 | type(spectre_fo_local_t), intent(out) :: loc | ||
| 162 | real(dp) :: y(NQ), dy(3, NQ) | ||
| 163 | |||
| 164 |
4/4✓ Branch 0 taken 104193 times.
✓ Branch 1 taken 34462 times.
✓ Branch 2 taken 66 times.
✓ Branch 3 taken 104127 times.
|
138655 | if (owner == 1 .and. abs(u(1)) < AXIS_OUTER) then |
| 165 | 66 | call axis_local(u, loc) | |
| 166 | 66 | return | |
| 167 | end if | ||
| 168 | 138589 | call evaluate_batch_splines_3d_der(vol_spl(owner), u, y, dy) | |
| 169 | call local_from_values(u(3), y(1), dy(:, 1), y(2), dy(:, 2), y(3), & | ||
| 170 | 138589 | dy(:, 3), y(4), dy(:, 4), loc) | |
| 171 | end subroutine spectre_fo_field_local | ||
| 172 | |||
| 173 |
14/14✓ Branch 0 taken 138589 times.
✓ Branch 1 taken 415767 times.
✓ Branch 2 taken 138589 times.
✓ Branch 3 taken 415767 times.
✓ Branch 4 taken 1247301 times.
✓ Branch 5 taken 415767 times.
✓ Branch 6 taken 138589 times.
✓ Branch 7 taken 415767 times.
✓ Branch 8 taken 1247301 times.
✓ Branch 9 taken 415767 times.
✓ Branch 10 taken 138589 times.
✓ Branch 11 taken 415767 times.
✓ Branch 12 taken 415767 times.
✓ Branch 13 taken 138589 times.
|
5266382 | pure subroutine local_from_values(zeta, r, dr, z, dz, at, dat, az, daz, loc) |
| 174 | real(dp), intent(in) :: zeta, r, dr(3), z, dz(3), at, dat(3), az, daz(3) | ||
| 175 | type(spectre_fo_local_t), intent(out) :: loc | ||
| 176 | real(dp) :: sqgb(3) | ||
| 177 | |||
| 178 | 138589 | call geometry_from_values(zeta, r, dr, z, dz, loc) | |
| 179 |
2/2✓ Branch 0 taken 415767 times.
✓ Branch 1 taken 138589 times.
|
554356 | loc%acov = [0.0_dp, at, az] |
| 180 | ! Match libneo: |J| B^i = eps^ijk d_j A_k with A_rho = 0. | ||
| 181 |
2/2✓ Branch 0 taken 415767 times.
✓ Branch 1 taken 138589 times.
|
554356 | sqgb = [daz(2) - dat(3), -daz(1), dat(1)] |
| 182 | 138589 | call set_magnetic_field(loc, sqgb) | |
| 183 | 138589 | end subroutine local_from_values | |
| 184 | |||
| 185 |
14/14✓ Branch 0 taken 415965 times.
✓ Branch 1 taken 138655 times.
✓ Branch 2 taken 415965 times.
✓ Branch 3 taken 138655 times.
✓ Branch 4 taken 1247895 times.
✓ Branch 5 taken 415965 times.
✓ Branch 6 taken 415965 times.
✓ Branch 7 taken 138655 times.
✓ Branch 8 taken 1247895 times.
✓ Branch 9 taken 415965 times.
✓ Branch 10 taken 415965 times.
✓ Branch 11 taken 138655 times.
✓ Branch 12 taken 415965 times.
✓ Branch 13 taken 138655 times.
|
4714270 | pure subroutine geometry_from_values(zeta, r, dr, z, dz, loc) |
| 186 | real(dp), intent(in) :: zeta, r, dr(3), z, dz(3) | ||
| 187 | type(spectre_fo_local_t), intent(out) :: loc | ||
| 188 | real(dp) :: c, s, cof(3, 3) | ||
| 189 | integer :: i | ||
| 190 | |||
| 191 | 138655 | c = cos(zeta) | |
| 192 | 138655 | s = sin(zeta) | |
| 193 |
2/2✓ Branch 0 taken 415965 times.
✓ Branch 1 taken 138655 times.
|
554620 | loc%x = [r*c, r*s, z] |
| 194 |
2/2✓ Branch 0 taken 415965 times.
✓ Branch 1 taken 138655 times.
|
554620 | do i = 1, 3 |
| 195 |
2/2✓ Branch 0 taken 1247895 times.
✓ Branch 1 taken 415965 times.
|
1802515 | loc%e(:, i) = [dr(i)*c, dr(i)*s, dz(i)] |
| 196 | end do | ||
| 197 | 138655 | loc%e(1, 3) = loc%e(1, 3) - r*s | |
| 198 | 138655 | loc%e(2, 3) = loc%e(2, 3) + r*c | |
| 199 | 138655 | cof(:, 1) = cross(loc%e(:, 2), loc%e(:, 3)) | |
| 200 | 138655 | cof(:, 2) = cross(loc%e(:, 3), loc%e(:, 1)) | |
| 201 | 138655 | cof(:, 3) = cross(loc%e(:, 1), loc%e(:, 2)) | |
| 202 |
2/2✓ Branch 0 taken 415965 times.
✓ Branch 1 taken 138655 times.
|
554620 | loc%jac = dot_product(loc%e(:, 1), cof(:, 1)) |
| 203 |
4/4✓ Branch 0 taken 415965 times.
✓ Branch 1 taken 138655 times.
✓ Branch 2 taken 1247895 times.
✓ Branch 3 taken 415965 times.
|
1802515 | loc%dual = cof/loc%jac |
| 204 | 138655 | end subroutine geometry_from_values | |
| 205 | |||
| 206 | 138655 | pure subroutine set_magnetic_field(loc, sqgb) | |
| 207 | type(spectre_fo_local_t), intent(inout) :: loc | ||
| 208 | real(dp), intent(in) :: sqgb(3) | ||
| 209 |
8/8✓ Branch 0 taken 415965 times.
✓ Branch 1 taken 138655 times.
✓ Branch 2 taken 415965 times.
✓ Branch 3 taken 138655 times.
✓ Branch 4 taken 1247895 times.
✓ Branch 5 taken 415965 times.
✓ Branch 6 taken 415965 times.
✓ Branch 7 taken 138655 times.
|
2634445 | loc%b = matmul(loc%e, sqgb)/abs(loc%jac) |
| 210 |
5/6✓ Branch 0 taken 415965 times.
✓ Branch 1 taken 138655 times.
✓ Branch 2 taken 415965 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 219965 times.
✓ Branch 5 taken 196000 times.
|
554620 | loc%bmod = norm2(loc%b) |
| 211 |
2/2✓ Branch 0 taken 415965 times.
✓ Branch 1 taken 138655 times.
|
554620 | loc%b = loc%b/max(loc%bmod, tiny(1.0_dp)) |
| 212 | 138655 | end subroutine set_magnetic_field | |
| 213 | |||
| 214 | 259 | subroutine axis_values(q, y, dy) | |
| 215 | ! Join representations of smooth geometry and the potential, including | ||
| 216 | ! the derivative of the weight. Taking curl afterwards preserves the | ||
| 217 | ! field consistency and avoids an artificial chart-boundary sheet. | ||
| 218 | real(dp), intent(in) :: q(3) | ||
| 219 | real(dp), intent(out) :: y(AXIS_NQ), dy(3, AXIS_NQ) | ||
| 220 | real(dp) :: ya(AXIS_NQ), da(3, AXIS_NQ), yb(AXIS_NQ), db(3, AXIS_NQ) | ||
| 221 | real(dp) :: p(NQ), dpdu(3, NQ), u(3), rho, c, s, f, df(3), t, w, dw(3) | ||
| 222 | integer :: j | ||
| 223 | |||
| 224 |
5/6✓ Branch 0 taken 518 times.
✓ Branch 1 taken 259 times.
✓ Branch 2 taken 515 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 515 times.
|
777 | rho = norm2(q(1:2)) |
| 225 |
2/2✓ Branch 0 taken 133 times.
✓ Branch 1 taken 126 times.
|
259 | if (rho <= FO_AXIS_RADIUS) then |
| 226 | 133 | call evaluate_batch_splines_3d_der(axis_spl, q, y, dy) | |
| 227 | 161 | return | |
| 228 | end if | ||
| 229 | 126 | c = q(1)/rho | |
| 230 | 126 | s = q(2)/rho | |
| 231 |
2/2✓ Branch 0 taken 378 times.
✓ Branch 1 taken 126 times.
|
504 | u = [rho, atan2(q(2), q(1)), q(3)] |
| 232 | 126 | call evaluate_batch_splines_3d_der(vol_spl(1), u, p, dpdu) | |
| 233 |
2/2✓ Branch 0 taken 252 times.
✓ Branch 1 taken 126 times.
|
378 | yb(1:2) = p(1:2) |
| 234 |
2/2✓ Branch 0 taken 252 times.
✓ Branch 1 taken 126 times.
|
378 | db(1, 1:2) = c*dpdu(1, 1:2) - s/rho*dpdu(2, 1:2) |
| 235 |
2/2✓ Branch 0 taken 252 times.
✓ Branch 1 taken 126 times.
|
378 | db(2, 1:2) = s*dpdu(1, 1:2) + c/rho*dpdu(2, 1:2) |
| 236 |
2/2✓ Branch 0 taken 252 times.
✓ Branch 1 taken 126 times.
|
378 | db(3, 1:2) = dpdu(3, 1:2) |
| 237 | 126 | f = p(3)/rho | |
| 238 | df = [c*(dpdu(1, 3) - f)/rho - s*dpdu(2, 3)/rho**2, & | ||
| 239 |
2/2✓ Branch 0 taken 378 times.
✓ Branch 1 taken 126 times.
|
504 | s*(dpdu(1, 3) - f)/rho + c*dpdu(2, 3)/rho**2, dpdu(3, 3)/rho] |
| 240 |
2/2✓ Branch 0 taken 378 times.
✓ Branch 1 taken 126 times.
|
504 | yb(3:5) = [-s*f, c*f, p(4)] |
| 241 |
2/2✓ Branch 0 taken 378 times.
✓ Branch 1 taken 126 times.
|
504 | db(:, 3) = [s*c*f/rho, -c*c*f/rho, 0.0_dp] - s*df |
| 242 |
2/2✓ Branch 0 taken 378 times.
✓ Branch 1 taken 126 times.
|
504 | db(:, 4) = [s*s*f/rho, -c*s*f/rho, 0.0_dp] + c*df |
| 243 | db(:, 5) = [c*dpdu(1, 4) - s/rho*dpdu(2, 4), & | ||
| 244 |
2/2✓ Branch 0 taken 378 times.
✓ Branch 1 taken 126 times.
|
504 | s*dpdu(1, 4) + c/rho*dpdu(2, 4), dpdu(3, 4)] |
| 245 |
2/2✓ Branch 0 taken 28 times.
✓ Branch 1 taken 98 times.
|
126 | if (rho >= AXIS_OUTER) then |
| 246 | 28 | y = yb | |
| 247 | 28 | dy = db | |
| 248 | 28 | return | |
| 249 | end if | ||
| 250 | 98 | call evaluate_batch_splines_3d_der(axis_spl, q, ya, da) | |
| 251 | 98 | t = (rho - FO_AXIS_RADIUS)/(AXIS_OUTER - FO_AXIS_RADIUS) | |
| 252 | ! C4 weight matches the continuity of the quintic potential splines. | ||
| 253 | 98 | w = 1 - t**5*(126 - 420*t + 540*t*t - 315*t**3 + 70*t**4) | |
| 254 | dw = -630*t**4*(1 - t)**4/(AXIS_OUTER - FO_AXIS_RADIUS) & | ||
| 255 |
2/2✓ Branch 0 taken 294 times.
✓ Branch 1 taken 98 times.
|
392 | *[c, s, 0.0_dp] |
| 256 |
2/2✓ Branch 0 taken 490 times.
✓ Branch 1 taken 98 times.
|
588 | y = w*ya + (1 - w)*yb |
| 257 |
2/2✓ Branch 0 taken 490 times.
✓ Branch 1 taken 98 times.
|
588 | do j = 1, AXIS_NQ |
| 258 |
2/2✓ Branch 0 taken 1470 times.
✓ Branch 1 taken 490 times.
|
2058 | dy(:, j) = w*da(:, j) + (1 - w)*db(:, j) + dw*(ya(j) - yb(j)) |
| 259 | end do | ||
| 260 | end subroutine axis_values | ||
| 261 | |||
| 262 |
14/14✓ Branch 0 taken 198 times.
✓ Branch 1 taken 66 times.
✓ Branch 2 taken 198 times.
✓ Branch 3 taken 66 times.
✓ Branch 4 taken 594 times.
✓ Branch 5 taken 198 times.
✓ Branch 6 taken 198 times.
✓ Branch 7 taken 66 times.
✓ Branch 8 taken 594 times.
✓ Branch 9 taken 198 times.
✓ Branch 10 taken 198 times.
✓ Branch 11 taken 66 times.
✓ Branch 12 taken 198 times.
✓ Branch 13 taken 66 times.
|
2244 | subroutine axis_local(u, loc) |
| 263 | real(dp), intent(in) :: u(3) | ||
| 264 | type(spectre_fo_local_t), intent(out) :: loc | ||
| 265 | real(dp) :: y(AXIS_NQ), dy(3, AXIS_NQ), q(3), c, s | ||
| 266 | real(dp) :: basis(3, 2), dual(3, 2), curl_a(3) | ||
| 267 | |||
| 268 | 66 | c = cos(u(2)) | |
| 269 | 66 | s = sin(u(2)) | |
| 270 |
2/2✓ Branch 0 taken 198 times.
✓ Branch 1 taken 66 times.
|
264 | q = [u(1)*c, u(1)*s, u(3)] |
| 271 | 66 | call axis_values(q, y, dy) | |
| 272 | 66 | call geometry_from_values(u(3), y(1), dy(:, 1), y(2), dy(:, 2), loc) | |
| 273 | curl_a = [dy(2, 5) - dy(3, 4), dy(3, 3) - dy(1, 5), & | ||
| 274 |
2/2✓ Branch 0 taken 198 times.
✓ Branch 1 taken 66 times.
|
264 | dy(1, 4) - dy(2, 3)] |
| 275 | 66 | call set_magnetic_field(loc, curl_a) | |
| 276 | ! Expose the existing polar basis to callers; the magnetic field was | ||
| 277 | ! evaluated with the nonsingular Cartesian-chart Jacobian above. | ||
| 278 |
4/4✓ Branch 0 taken 132 times.
✓ Branch 1 taken 66 times.
✓ Branch 2 taken 396 times.
✓ Branch 3 taken 132 times.
|
594 | basis = loc%e(:, 1:2) |
| 279 |
4/4✓ Branch 0 taken 132 times.
✓ Branch 1 taken 66 times.
✓ Branch 2 taken 396 times.
✓ Branch 3 taken 132 times.
|
594 | dual = loc%dual(:, 1:2) |
| 280 |
2/2✓ Branch 0 taken 198 times.
✓ Branch 1 taken 66 times.
|
264 | loc%e(:, 1) = c*basis(:, 1) + s*basis(:, 2) |
| 281 |
2/2✓ Branch 0 taken 198 times.
✓ Branch 1 taken 66 times.
|
264 | loc%e(:, 2) = u(1)*(-s*basis(:, 1) + c*basis(:, 2)) |
| 282 |
2/2✓ Branch 0 taken 198 times.
✓ Branch 1 taken 66 times.
|
264 | loc%dual(:, 1) = c*dual(:, 1) + s*dual(:, 2) |
| 283 |
2/2✓ Branch 0 taken 198 times.
✓ Branch 1 taken 66 times.
|
264 | loc%dual(:, 2) = 0.0_dp |
| 284 |
1/2✓ Branch 0 taken 66 times.
✗ Branch 1 not taken.
|
66 | if (abs(u(1)) > tiny(1.0_dp)) & |
| 285 |
2/2✓ Branch 0 taken 198 times.
✓ Branch 1 taken 66 times.
|
264 | loc%dual(:, 2) = (-s*dual(:, 1) + c*dual(:, 2))/u(1) |
| 286 | 66 | loc%jac = loc%jac*u(1) | |
| 287 |
2/2✓ Branch 0 taken 198 times.
✓ Branch 1 taken 66 times.
|
264 | loc%acov = [c*y(3) + s*y(4), u(1)*(-s*y(3) + c*y(4)), y(5)] |
| 288 | 66 | end subroutine axis_local | |
| 289 | |||
| 290 | 70 | subroutine locate_axis(rt, zt, phi, v, status) | |
| 291 | real(dp), intent(in) :: rt, zt, phi | ||
| 292 | real(dp), intent(inout) :: v(3) | ||
| 293 | integer, intent(out) :: status | ||
| 294 | real(dp) :: q(3), y(AXIS_NQ), dy(3, AXIS_NQ), fr, fz, det | ||
| 295 | real(dp) :: delta(2), step, theta | ||
| 296 | integer :: it | ||
| 297 | |||
| 298 |
2/2✓ Branch 0 taken 210 times.
✓ Branch 1 taken 70 times.
|
280 | q = [v(1)*cos(v(2)), v(1)*sin(v(2)), phi] |
| 299 |
1/2✓ Branch 0 taken 193 times.
✗ Branch 1 not taken.
|
193 | do it = 1, 60 |
| 300 | 193 | call axis_values(q, y, dy) | |
| 301 | 193 | fr = y(1) - rt | |
| 302 | 193 | fz = y(2) - zt | |
| 303 | 193 | det = dy(1, 1)*dy(2, 2) - dy(2, 1)*dy(1, 2) | |
| 304 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 193 times.
|
193 | if (abs(det) <= tiny(1.0_dp)) exit |
| 305 | delta = [(fr*dy(2, 2) - fz*dy(2, 1))/det, & | ||
| 306 |
2/2✓ Branch 0 taken 386 times.
✓ Branch 1 taken 193 times.
|
579 | (fz*dy(1, 1) - fr*dy(1, 2))/det] |
| 307 |
6/8✗ Branch 0 not taken.
✓ Branch 1 taken 193 times.
✓ Branch 2 taken 193 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 386 times.
✓ Branch 5 taken 193 times.
✓ Branch 6 taken 269 times.
✓ Branch 7 taken 117 times.
|
772 | step = maxval(abs(delta)) |
| 308 |
4/4✓ Branch 0 taken 9 times.
✓ Branch 1 taken 184 times.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 18 times.
|
211 | if (step > AXIS_OUTER) delta = delta*AXIS_OUTER/step |
| 309 |
2/2✓ Branch 0 taken 386 times.
✓ Branch 1 taken 193 times.
|
579 | q(1:2) = q(1:2) - delta |
| 310 |
9/12✗ Branch 0 not taken.
✓ Branch 1 taken 193 times.
✓ Branch 2 taken 193 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 386 times.
✓ Branch 5 taken 193 times.
✓ Branch 6 taken 269 times.
✓ Branch 7 taken 117 times.
✓ Branch 8 taken 70 times.
✓ Branch 9 taken 123 times.
✓ Branch 10 taken 70 times.
✗ Branch 11 not taken.
|
772 | if (maxval(abs(delta)) <= 1.0e-12_dp .and. & |
| 311 | ✗ | abs(fr) + abs(fz) <= 1.0e-12_dp*max(rt, 1.0_dp)) then | |
| 312 | 70 | theta = v(2) | |
| 313 |
5/8✓ Branch 0 taken 140 times.
✓ Branch 1 taken 70 times.
✓ Branch 2 taken 140 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 140 times.
✓ Branch 6 taken 70 times.
✗ Branch 7 not taken.
|
210 | if (norm2(q(1:2)) > tiny(1.0_dp)) theta = atan2(q(2), q(1)) |
| 314 | v = [norm2(q(1:2)), & | ||
| 315 |
6/8✓ Branch 0 taken 140 times.
✓ Branch 1 taken 70 times.
✓ Branch 2 taken 140 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 140 times.
✓ Branch 6 taken 210 times.
✓ Branch 7 taken 70 times.
|
420 | theta + TWOPI*nint((v(2) - theta)/TWOPI), phi] |
| 316 | 70 | status = 0 | |
| 317 | 70 | return | |
| 318 | end if | ||
| 319 | end do | ||
| 320 | ✗ | status = 2 | |
| 321 | end subroutine locate_axis | ||
| 322 | |||
| 323 | 12211 | subroutine spectre_fo_field_locate(x, u, owner, status) | |
| 324 | !> Logical point of the Cartesian x (cm): zeta is the cylindrical angle, | ||
| 325 | !> (rho, theta) a Newton root of (R, Z) at that zeta, starting from u in | ||
| 326 | !> volume owner. A root outside the owner's slab moves to the adjacent | ||
| 327 | !> volume and is refined there. status = 0 on success, 1 if the point is | ||
| 328 | !> beyond the outer boundary, 2 if | ||
| 329 | !> Newton fails. | ||
| 330 | real(dp), intent(in) :: x(3) | ||
| 331 | real(dp), intent(inout) :: u(3) | ||
| 332 | integer, intent(inout) :: owner | ||
| 333 | integer, intent(out) :: status | ||
| 334 | real(dp) :: rt, zt, phi, v(3), y(NQ), dy(3, NQ), fr, fz, det, dr, dt | ||
| 335 | real(dp) :: scale | ||
| 336 | integer :: it, own, switches | ||
| 337 | |||
| 338 | 12211 | rt = sqrt(x(1)**2 + x(2)**2) | |
| 339 | 12211 | zt = x(3) | |
| 340 | 12211 | phi = atan2(x(2), x(1)) | |
| 341 | 12211 | v = u | |
| 342 | 12211 | v(3) = phi + TWOPI*nint((u(3) - phi)/TWOPI) | |
| 343 | 12211 | phi = v(3) | |
| 344 | 12211 | own = owner | |
| 345 | 12211 | scale = max(rt, 1.0_dp) | |
| 346 | 12211 | switches = 0 | |
| 347 | 12211 | status = 2 | |
| 348 |
1/2✓ Branch 0 taken 47679 times.
✗ Branch 1 not taken.
|
47679 | do it = 1, 60 |
| 349 |
4/4✓ Branch 0 taken 30894 times.
✓ Branch 1 taken 16785 times.
✓ Branch 2 taken 70 times.
✓ Branch 3 taken 30824 times.
|
47679 | if (own == 1 .and. v(1) < AXIS_OUTER) then |
| 350 | 70 | call locate_axis(rt, zt, phi, v, status) | |
| 351 |
1/2✓ Branch 0 taken 70 times.
✗ Branch 1 not taken.
|
152 | if (status /= 0) return |
| 352 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 6 times.
|
70 | if (v(1) < AXIS_OUTER) then |
| 353 | 64 | u = v | |
| 354 | 64 | owner = 1 | |
| 355 | 64 | return | |
| 356 | end if | ||
| 357 | 6 | status = 2 | |
| 358 | end if | ||
| 359 | 47615 | call evaluate_batch_splines_3d_der(vol_spl(own), v, y, dy) | |
| 360 | 47615 | fr = y(1) - rt | |
| 361 | 47615 | fz = y(2) - zt | |
| 362 | 47615 | det = dy(1, 1)*dy(2, 2) - dy(2, 1)*dy(1, 2) | |
| 363 |
1/2✓ Branch 0 taken 47615 times.
✗ Branch 1 not taken.
|
47615 | if (abs(det) <= tiny(1.0_dp)) return |
| 364 | 47615 | dr = (fr*dy(2, 2) - fz*dy(2, 1))/det | |
| 365 | 47615 | dt = (fz*dy(1, 1) - fr*dy(1, 2))/det | |
| 366 | 47615 | dr = max(-0.25_dp, min(0.25_dp, dr)) | |
| 367 | 47615 | v(1) = v(1) - dr | |
| 368 | 47615 | v(2) = v(2) - dt | |
| 369 |
4/4✓ Branch 0 taken 30830 times.
✓ Branch 1 taken 16785 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 30825 times.
|
47615 | if (own == 1 .and. v(1) < 0.0_dp) then |
| 370 | ! Polar chart flip: the axis-volume embedding uses rho**m | ||
| 371 | ! (rho**2 for m=0), so this is the same physical point. | ||
| 372 | 5 | v(1) = -v(1) | |
| 373 | 5 | v(2) = v(2) + 0.5_dp*TWOPI | |
| 374 | end if | ||
| 375 |
5/6✓ Branch 0 taken 1089 times.
✓ Branch 1 taken 46526 times.
✓ Branch 2 taken 1021 times.
✓ Branch 3 taken 68 times.
✓ Branch 4 taken 1021 times.
✗ Branch 5 not taken.
|
47615 | if (v(1) > real(own, dp) .and. own < fo_field_mvol .and. & |
| 376 | switches < 8) then | ||
| 377 | 1021 | own = own + 1 | |
| 378 | 1021 | switches = switches + 1 | |
| 379 |
3/4✓ Branch 0 taken 458 times.
✓ Branch 1 taken 46136 times.
✓ Branch 2 taken 458 times.
✗ Branch 3 not taken.
|
46594 | else if (v(1) < real(own - 1, dp) .and. own > 1 .and. & |
| 380 | switches < 8) then | ||
| 381 | 458 | own = own - 1 | |
| 382 | 458 | switches = switches + 1 | |
| 383 | end if | ||
| 384 | ! Coordinate and physical residuals must both converge. Requiring | ||
| 385 | ! a sub-ulp physical displacement can otherwise oscillate forever | ||
| 386 | ! after the cylindrical embedding has reached roundoff. | ||
| 387 |
4/4✓ Branch 0 taken 12216 times.
✓ Branch 1 taken 35399 times.
✓ Branch 2 taken 12147 times.
✓ Branch 3 taken 69 times.
|
47615 | if (abs(dr) <= 1.0e-12_dp .and. abs(fr) + abs(fz) <= & |
| 388 | 12147 | 1.0e-12_dp*scale) exit | |
| 389 | end do | ||
| 390 |
1/2✓ Branch 0 taken 12147 times.
✗ Branch 1 not taken.
|
12147 | if (it > 60) return |
| 391 |
5/6✓ Branch 0 taken 12129 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 7806 times.
✓ Branch 3 taken 4323 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 7806 times.
|
12147 | if (v(1) > real(fo_field_mvol, dp) .or. & |
| 392 | (own == 1 .and. v(1) <= 0.0_dp)) then | ||
| 393 | ! Preserve the converged exterior root so an entry map can | ||
| 394 | ! distinguish an initial particle outside the wall from failure. | ||
| 395 | 18 | u = v | |
| 396 | 18 | owner = own | |
| 397 | 18 | status = 1 | |
| 398 | 18 | return | |
| 399 | end if | ||
| 400 | ! Roundoff at a sheet: the root lies on the owner's edge. | ||
| 401 |
2/4✓ Branch 0 taken 12129 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 12129 times.
✗ Branch 3 not taken.
|
12129 | if (v(1) < real(own - 1, dp) - 1.0e-12_dp .or. & |
| 402 | v(1) > real(own, dp) + 1.0e-12_dp) return | ||
| 403 | 12129 | v(1) = max(real(own - 1, dp), min(real(own, dp), v(1))) | |
| 404 | 12129 | u = v | |
| 405 | 12129 | owner = own | |
| 406 | 12129 | status = 0 | |
| 407 | end subroutine spectre_fo_field_locate | ||
| 408 | |||
| 409 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 415965 times.
|
415965 | pure function cross(a, b) result(c) |
| 410 | real(dp), intent(in) :: a(3), b(3) | ||
| 411 | real(dp) :: c(3) | ||
| 412 | c = [a(2)*b(3) - a(3)*b(2), a(3)*b(1) - a(1)*b(3), & | ||
| 413 |
2/2✓ Branch 0 taken 1247895 times.
✓ Branch 1 taken 415965 times.
|
1663860 | a(1)*b(2) - a(2)*b(1)] |
| 414 | 415965 | end function cross | |
| 415 | ✗ | end module spectre_fo_field | |
| 416 |