| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module field_splined | ||
| 2 | !> Splined field decorator for magnetic fields. | ||
| 3 | !> | ||
| 4 | !> Wraps any magnetic_field_t and provides fast splined evaluation. | ||
| 5 | !> Splines are built in scaled coordinates (default: r = sqrt(s) for better | ||
| 6 | !> resolution near axis, matching libneo convention). | ||
| 7 | !> | ||
| 8 | !> The reference coordinate system (ref_coords) is used for: | ||
| 9 | !> - Mapping from reference coordinates to Cartesian (evaluate_cart) | ||
| 10 | !> - Computing covariant basis vectors (covariant_basis) | ||
| 11 | !> This allows proper support for VMEC, chartmap, or any coordinate system. | ||
| 12 | !> | ||
| 13 | !> Usage: | ||
| 14 | !> type(coils_field_t) :: raw_coils | ||
| 15 | !> type(splined_field_t) :: splined | ||
| 16 | !> call create_coils_field('coils.dat', raw_coils) | ||
| 17 | !> call create_splined_field(raw_coils, ref_coords, splined) | ||
| 18 | !> call splined%evaluate(x, Acov, hcov, Bmod) | ||
| 19 | !> | ||
| 20 | !> The scaling parameter controls how the first coordinate is transformed: | ||
| 21 | !> - sqrt_s_scaling_t (default): grid in r = sqrt(s), better axis resolution | ||
| 22 | !> | ||
| 23 | !> Limitations: | ||
| 24 | !> - Source field must evaluate in Cartesian coordinates | ||
| 25 | !> - sqgBctr output not supported (error stop if requested) | ||
| 26 | |||
| 27 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 28 | use interpolate, only: BatchSplineData3D, construct_batch_splines_3d, & | ||
| 29 | evaluate_batch_splines_3d, evaluate_batch_splines_3d_der, & | ||
| 30 | destroy_batch_splines_3d | ||
| 31 | use field_base, only: magnetic_field_t | ||
| 32 | use field_vmec, only: vmec_field_t | ||
| 33 | use field_boozer_chartmap, only: boozer_chartmap_field_t | ||
| 34 | use libneo_coordinates, only: coordinate_system_t, vmec_coordinate_system_t, & | ||
| 35 | chartmap_coordinate_system_t, RHO_TOR, RHO_POL, & | ||
| 36 | PSI_TOR_NORM, PSI_POL_NORM | ||
| 37 | use coordinate_scaling, only: coordinate_scaling_t, sqrt_s_scaling_t, & | ||
| 38 | identity_scaling_t | ||
| 39 | use cartesian_coordinates, only: cartesian_coordinate_system_t | ||
| 40 | |||
| 41 | implicit none | ||
| 42 | |||
| 43 | type, extends(magnetic_field_t) :: splined_field_t | ||
| 44 | type(BatchSplineData3D) :: splines | ||
| 45 | logical :: initialized = .false. | ||
| 46 | contains | ||
| 47 | procedure :: evaluate => splined_evaluate | ||
| 48 | procedure :: evaluate_with_der => splined_evaluate_with_der | ||
| 49 | final :: splined_field_cleanup | ||
| 50 | end type splined_field_t | ||
| 51 | |||
| 52 | contains | ||
| 53 | |||
| 54 | 1091921 | subroutine splined_evaluate(self, x, Acov, hcov, Bmod, sqgBctr) | |
| 55 | !> Evaluate splined field at x in scaled coordinates. | ||
| 56 | class(splined_field_t), intent(in) :: self | ||
| 57 | real(dp), intent(in) :: x(3) | ||
| 58 | real(dp), intent(out) :: Acov(3), hcov(3), Bmod | ||
| 59 | real(dp), intent(out), optional :: sqgBctr(3) | ||
| 60 | |||
| 61 | real(dp) :: y_batch(7) | ||
| 62 | |||
| 63 | 1091921 | call evaluate_batch_splines_3d(self%splines, x, y_batch) | |
| 64 | |||
| 65 | 1091921 | Acov(1) = y_batch(1) | |
| 66 | 1091921 | Acov(2) = y_batch(2) | |
| 67 | 1091921 | Acov(3) = y_batch(3) | |
| 68 | |||
| 69 | 1091921 | hcov(1) = y_batch(4) | |
| 70 | 1091921 | hcov(2) = y_batch(5) | |
| 71 | 1091921 | hcov(3) = y_batch(6) | |
| 72 | |||
| 73 | 1091921 | Bmod = y_batch(7) | |
| 74 | |||
| 75 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1091921 times.
|
1091921 | if (present(sqgBctr)) then |
| 76 | ✗ | error stop "sqgBctr not implemented for splined_field_t" | |
| 77 | end if | ||
| 78 | 1091921 | end subroutine splined_evaluate | |
| 79 | |||
| 80 | 9537 | subroutine splined_evaluate_with_der(self, x, Acov, hcov, Bmod, dAcov, dhcov, & | |
| 81 | dBmod, sqgBctr) | ||
| 82 | !> Evaluate splined field at x in scaled coordinates and return | ||
| 83 | !> first derivatives with respect to the same coordinates x. | ||
| 84 | class(splined_field_t), intent(in) :: self | ||
| 85 | real(dp), intent(in) :: x(3) | ||
| 86 | real(dp), intent(out) :: Acov(3), hcov(3), Bmod | ||
| 87 | real(dp), intent(out) :: dAcov(3, 3), dhcov(3, 3), dBmod(3) | ||
| 88 | real(dp), intent(out), optional :: sqgBctr(3) | ||
| 89 | |||
| 90 | real(dp) :: y_batch(7) | ||
| 91 | real(dp) :: dy_batch(3, 7) | ||
| 92 | |||
| 93 | 9537 | call evaluate_batch_splines_3d_der(self%splines, x, y_batch, dy_batch) | |
| 94 | |||
| 95 | 9537 | Acov(1) = y_batch(1) | |
| 96 | 9537 | Acov(2) = y_batch(2) | |
| 97 | 9537 | Acov(3) = y_batch(3) | |
| 98 | |||
| 99 | 9537 | hcov(1) = y_batch(4) | |
| 100 | 9537 | hcov(2) = y_batch(5) | |
| 101 | 9537 | hcov(3) = y_batch(6) | |
| 102 | |||
| 103 | 9537 | Bmod = y_batch(7) | |
| 104 | |||
| 105 |
2/2✓ Branch 0 taken 9537 times.
✓ Branch 1 taken 28611 times.
|
38148 | dAcov(:, 1) = dy_batch(:, 1) |
| 106 |
2/2✓ Branch 0 taken 9537 times.
✓ Branch 1 taken 28611 times.
|
38148 | dAcov(:, 2) = dy_batch(:, 2) |
| 107 |
2/2✓ Branch 0 taken 9537 times.
✓ Branch 1 taken 28611 times.
|
38148 | dAcov(:, 3) = dy_batch(:, 3) |
| 108 | |||
| 109 |
2/2✓ Branch 0 taken 9537 times.
✓ Branch 1 taken 28611 times.
|
38148 | dhcov(:, 1) = dy_batch(:, 4) |
| 110 |
2/2✓ Branch 0 taken 9537 times.
✓ Branch 1 taken 28611 times.
|
38148 | dhcov(:, 2) = dy_batch(:, 5) |
| 111 |
2/2✓ Branch 0 taken 9537 times.
✓ Branch 1 taken 28611 times.
|
38148 | dhcov(:, 3) = dy_batch(:, 6) |
| 112 | |||
| 113 |
2/2✓ Branch 0 taken 28611 times.
✓ Branch 1 taken 9537 times.
|
38148 | dBmod = dy_batch(:, 7) |
| 114 | |||
| 115 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9537 times.
|
9537 | if (present(sqgBctr)) then |
| 116 | ✗ | error stop "sqgBctr not implemented for splined_field_t" | |
| 117 | end if | ||
| 118 | 9537 | end subroutine splined_evaluate_with_der | |
| 119 | |||
| 120 | 11 | subroutine create_splined_field(source, ref_coords, field, scaling, & | |
| 121 | n_r, n_th, n_phi, xmin, xmax) | ||
| 122 | !> Create splined field by sampling source field onto ref_coords grid. | ||
| 123 | !> Grid is built in scaled coordinates (default: r = sqrt(s)). | ||
| 124 | class(magnetic_field_t), intent(in) :: source | ||
| 125 | class(coordinate_system_t), intent(in), target :: ref_coords | ||
| 126 | type(splined_field_t), intent(out) :: field | ||
| 127 | class(coordinate_scaling_t), intent(in), optional, target :: scaling | ||
| 128 | integer, intent(in), optional :: n_r, n_th, n_phi | ||
| 129 | real(dp), intent(in), optional :: xmin(3), xmax(3) | ||
| 130 | |||
| 131 | call build_splines(source, ref_coords, field%splines, scaling, & | ||
| 132 | 11 | n_r, n_th, n_phi, xmin, xmax) | |
| 133 | |||
| 134 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 11 times.
|
11 | allocate (field%coords, source=ref_coords) |
| 135 | 11 | field%initialized = .true. | |
| 136 | 11 | end subroutine create_splined_field | |
| 137 | |||
| 138 |
15/16✓ Branch 0 taken 33 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 33 times.
✓ Branch 3 taken 11 times.
✓ Branch 4 taken 33 times.
✓ Branch 5 taken 11 times.
✓ Branch 6 taken 33 times.
✓ Branch 7 taken 11 times.
✓ Branch 8 taken 33 times.
✓ Branch 9 taken 11 times.
✓ Branch 10 taken 33 times.
✓ Branch 11 taken 11 times.
✓ Branch 12 taken 33 times.
✓ Branch 13 taken 11 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 11 times.
|
242 | subroutine build_splines(source, ref_coords, spl, scaling, & |
| 139 | n_r_in, n_th_in, n_phi_in, xmin_in, xmax_in) | ||
| 140 | !> Build splines by sampling source field on uniform grid. | ||
| 141 | !> Grid is in scaled coordinates. Default scaling: r = sqrt(s). | ||
| 142 | class(magnetic_field_t), intent(in) :: source | ||
| 143 | class(coordinate_system_t), intent(in) :: ref_coords | ||
| 144 | type(BatchSplineData3D), intent(out) :: spl | ||
| 145 | class(coordinate_scaling_t), intent(in), optional, target :: scaling | ||
| 146 | integer, intent(in), optional :: n_r_in, n_th_in, n_phi_in | ||
| 147 | real(dp), intent(in), optional :: xmin_in(3), xmax_in(3) | ||
| 148 | |||
| 149 | class(coordinate_scaling_t), pointer :: scaling_ptr | ||
| 150 | type(sqrt_s_scaling_t), target :: sqrt_s_scaling | ||
| 151 | type(identity_scaling_t), target :: identity_scaling | ||
| 152 | |||
| 153 | integer :: n_r, n_th, n_phi | ||
| 154 | real(dp) :: xmin(3), xmax(3) | ||
| 155 | |||
| 156 | integer, parameter :: order(3) = [5, 5, 5] | ||
| 157 | logical, parameter :: periodic(3) = [.False., .True., .True.] | ||
| 158 | |||
| 159 | real(dp) :: h_r, h_th, h_phi | ||
| 160 | |||
| 161 | 11 | real(dp), dimension(:, :, :), allocatable :: Ar, Ath, Aphi, hr, hth, hphi | |
| 162 | 11 | real(dp), dimension(:, :, :), allocatable :: Bmod_arr | |
| 163 | 11 | real(dp), dimension(:, :, :, :), allocatable :: y_batch | |
| 164 | integer :: dims(3) | ||
| 165 | |||
| 166 | call configure_spline_grid(ref_coords, n_r_in, n_th_in, n_phi_in, xmin_in, & | ||
| 167 | 11 | xmax_in, n_r, n_th, n_phi, xmin, xmax) | |
| 168 | |||
| 169 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
11 | if (present(scaling)) then |
| 170 | ✗ | scaling_ptr => scaling | |
| 171 | else | ||
| 172 | scaling_ptr => get_default_scaling(ref_coords, sqrt_s_scaling, & | ||
| 173 | 11 | identity_scaling) | |
| 174 | end if | ||
| 175 | |||
| 176 | 11 | h_r = (xmax(1) - xmin(1))/(n_r - 1) | |
| 177 | 11 | h_th = (xmax(2) - xmin(2))/(n_th - 1) | |
| 178 | 11 | h_phi = (xmax(3) - xmin(3))/(n_phi - 1) | |
| 179 | |||
| 180 | call sample_spline_arrays(source, ref_coords, scaling_ptr, n_r, n_th, n_phi, & | ||
| 181 | xmin, h_r, h_th, h_phi, Ar, Ath, Aphi, hr, hth, & | ||
| 182 | hphi, Bmod_arr) | ||
| 183 | |||
| 184 |
2/2✓ Branch 0 taken 33 times.
✓ Branch 1 taken 11 times.
|
44 | dims = shape(Ar) |
| 185 |
11/22✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 11 times.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 11 times.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 11 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 11 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 11 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 11 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 11 times.
|
66 | allocate (y_batch(dims(1), dims(2), dims(3), 7)) |
| 186 | |||
| 187 |
6/6✓ Branch 0 taken 372 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 15220 times.
✓ Branch 3 taken 372 times.
✓ Branch 4 taken 645316 times.
✓ Branch 5 taken 15220 times.
|
660919 | y_batch(:, :, :, 1) = Ar |
| 188 |
6/6✓ Branch 0 taken 372 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 15220 times.
✓ Branch 3 taken 372 times.
✓ Branch 4 taken 645316 times.
✓ Branch 5 taken 15220 times.
|
660919 | y_batch(:, :, :, 2) = Ath |
| 189 |
6/6✓ Branch 0 taken 372 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 15220 times.
✓ Branch 3 taken 372 times.
✓ Branch 4 taken 645316 times.
✓ Branch 5 taken 15220 times.
|
660919 | y_batch(:, :, :, 3) = Aphi |
| 190 |
6/6✓ Branch 0 taken 372 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 15220 times.
✓ Branch 3 taken 372 times.
✓ Branch 4 taken 645316 times.
✓ Branch 5 taken 15220 times.
|
660919 | y_batch(:, :, :, 4) = hr |
| 191 |
6/6✓ Branch 0 taken 372 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 15220 times.
✓ Branch 3 taken 372 times.
✓ Branch 4 taken 645316 times.
✓ Branch 5 taken 15220 times.
|
660919 | y_batch(:, :, :, 5) = hth |
| 192 |
6/6✓ Branch 0 taken 372 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 15220 times.
✓ Branch 3 taken 372 times.
✓ Branch 4 taken 645316 times.
✓ Branch 5 taken 15220 times.
|
660919 | y_batch(:, :, :, 6) = hphi |
| 193 |
6/6✓ Branch 0 taken 372 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 15220 times.
✓ Branch 3 taken 372 times.
✓ Branch 4 taken 645316 times.
✓ Branch 5 taken 15220 times.
|
660919 | y_batch(:, :, :, 7) = Bmod_arr |
| 194 | |||
| 195 | 11 | call construct_batch_splines_3d(xmin, xmax, y_batch, order, periodic, spl) | |
| 196 | |||
| 197 |
8/16✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 11 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 11 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 11 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 11 times.
✗ Branch 17 not taken.
|
33 | end subroutine build_splines |
| 198 | |||
| 199 | 11 | subroutine configure_spline_grid(ref_coords, n_r_in, n_th_in, n_phi_in, xmin_in, & | |
| 200 | xmax_in, n_r, n_th, n_phi, xmin, xmax) | ||
| 201 | use new_vmec_stuff_mod, only: nper | ||
| 202 | use util, only: twopi | ||
| 203 | class(coordinate_system_t), intent(in) :: ref_coords | ||
| 204 | integer, intent(in), optional :: n_r_in, n_th_in, n_phi_in | ||
| 205 | real(dp), intent(in), optional :: xmin_in(3), xmax_in(3) | ||
| 206 | integer, intent(out) :: n_r, n_th, n_phi | ||
| 207 | real(dp), intent(out) :: xmin(3), xmax(3) | ||
| 208 | |||
| 209 | 11 | n_r = 62 | |
| 210 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 6 times.
|
11 | if (present(n_r_in)) n_r = n_r_in |
| 211 | 11 | n_th = 63 | |
| 212 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 6 times.
|
11 | if (present(n_th_in)) n_th = n_th_in |
| 213 | 11 | n_phi = 64 | |
| 214 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 6 times.
|
11 | if (present(n_phi_in)) n_phi = n_phi_in |
| 215 | |||
| 216 | select type (cs => ref_coords) | ||
| 217 | class is (chartmap_coordinate_system_t) | ||
| 218 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (.not. present(n_r_in)) n_r = cs%nrho |
| 219 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (.not. present(n_th_in)) n_th = cs%ntheta |
| 220 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (.not. present(n_phi_in)) n_phi = cs%nzeta |
| 221 | class default | ||
| 222 | continue | ||
| 223 | end select | ||
| 224 | |||
| 225 | 11 | xmin = [1d-12, 0d0, 0d0] | |
| 226 |
2/2✓ Branch 0 taken 33 times.
✓ Branch 1 taken 11 times.
|
44 | xmax = [1d0, twopi, twopi/nper] |
| 227 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
|
11 | if (present(xmin_in)) xmin = xmin_in |
| 228 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
|
11 | if (present(xmax_in)) xmax = xmax_in |
| 229 | |||
| 230 | select type (cs => ref_coords) | ||
| 231 | class is (chartmap_coordinate_system_t) | ||
| 232 | block | ||
| 233 | real(dp) :: rho_min, rho_max, rho_max_safe, h_rho | ||
| 234 | |||
| 235 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (cs%has_spl_rz) then |
| 236 | 4 | rho_min = cs%spl_rz%x_min(1) | |
| 237 | rho_max = cs%spl_rz%x_min(1) + cs%spl_rz%h_step(1) * & | ||
| 238 | 4 | real(cs%spl_rz%num_points(1) - 1, dp) | |
| 239 | 4 | h_rho = cs%spl_rz%h_step(1) | |
| 240 | else | ||
| 241 | ✗ | rho_min = cs%spl_cart%x_min(1) | |
| 242 | rho_max = cs%spl_cart%x_min(1) + cs%spl_cart%h_step(1) * & | ||
| 243 | ✗ | real(cs%spl_cart%num_points(1) - 1, dp) | |
| 244 | ✗ | h_rho = cs%spl_cart%h_step(1) | |
| 245 | end if | ||
| 246 | |||
| 247 | ! Sampling the VMEC field requires inverting VMEC coordinates. This | ||
| 248 | ! inversion can fail exactly at the outermost rho plane, so sample | ||
| 249 | ! slightly inside while keeping essentially the full range. | ||
| 250 | 4 | rho_max_safe = rho_max | |
| 251 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (h_rho > 0d0) rho_max_safe = rho_max - 0.1d0*h_rho |
| 252 | |||
| 253 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (.not. present(xmin_in)) xmin(1) = max(xmin(1), rho_min) |
| 254 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (.not. present(xmax_in)) xmax(1) = min(xmax(1), rho_max_safe) |
| 255 | end block | ||
| 256 | class default | ||
| 257 | continue | ||
| 258 | end select | ||
| 259 | 11 | end subroutine configure_spline_grid | |
| 260 | |||
| 261 | ✗ | subroutine print_progress(i, n) | |
| 262 | integer, intent(in) :: i, n | ||
| 263 | #ifndef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 264 | return | ||
| 265 | #endif | ||
| 266 | !$omp critical | ||
| 267 | if (mod(i, max(1, n/10)) == 0) then | ||
| 268 | write (*, '(a,f6.2,a)', advance='no') char(13), 100.0*i/n, ' %' | ||
| 269 | call flush (6) | ||
| 270 | end if | ||
| 271 | !$omp end critical | ||
| 272 | end subroutine print_progress | ||
| 273 | |||
| 274 | 11 | subroutine sample_spline_arrays(source, ref_coords, scaling_ptr, n_r, & | |
| 275 | n_th, n_phi, & | ||
| 276 | xmin, h_r, h_th, h_phi, Ar, Ath, & | ||
| 277 | Aphi, hr, hth, & | ||
| 278 | hphi, Bmod_arr) | ||
| 279 | class(magnetic_field_t), intent(in) :: source | ||
| 280 | class(coordinate_system_t), intent(in) :: ref_coords | ||
| 281 | class(coordinate_scaling_t), pointer, intent(in) :: scaling_ptr | ||
| 282 | integer, intent(in) :: n_r, n_th, n_phi | ||
| 283 | real(dp), intent(in) :: xmin(3), h_r, h_th, h_phi | ||
| 284 | real(dp), allocatable, intent(out) :: Ar(:, :, :), Ath(:, :, :), & | ||
| 285 | Aphi(:, :, :) | ||
| 286 | real(dp), allocatable, intent(out) :: hr(:, :, :), hth(:, :, :), & | ||
| 287 | hphi(:, :, :) | ||
| 288 | real(dp), allocatable, intent(out) :: Bmod_arr(:, :, :) | ||
| 289 | |||
| 290 | real(dp) :: x_scaled(3), x_ref(3), x_cart(3) | ||
| 291 | real(dp) :: e_cov(3, 3), scaling_jac(3) | ||
| 292 | real(dp) :: Bmod | ||
| 293 | real(dp) :: Acov(3), hcov(3) | ||
| 294 | logical :: needs_scaling | ||
| 295 | logical :: is_identity_scaling, is_sqrt_s_scaling | ||
| 296 | integer :: i_r, i_th, i_phi, i_ctr | ||
| 297 | |||
| 298 |
20/40✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 11 times.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 11 times.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 11 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 11 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 11 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 11 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 11 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 11 times.
✓ Branch 27 taken 11 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 30 taken 11 times.
✓ Branch 31 taken 11 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✓ Branch 34 taken 11 times.
✓ Branch 35 taken 11 times.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✓ Branch 38 taken 11 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 11 times.
✗ Branch 43 not taken.
✓ Branch 44 taken 11 times.
|
88 | allocate (Ar(n_r, n_th, n_phi), Ath(n_r, n_th, n_phi)) |
| 299 |
9/18✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 11 times.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 11 times.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 11 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 11 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 11 times.
|
44 | allocate (Aphi(n_r, n_th, n_phi)) |
| 300 |
18/36✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 11 times.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 11 times.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 11 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 11 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 11 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 11 times.
✓ Branch 23 taken 11 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 11 times.
✓ Branch 27 taken 11 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 30 taken 11 times.
✓ Branch 31 taken 11 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✓ Branch 34 taken 11 times.
✗ Branch 36 not taken.
✓ Branch 37 taken 11 times.
✗ Branch 39 not taken.
✓ Branch 40 taken 11 times.
|
77 | allocate (hr(n_r, n_th, n_phi), hth(n_r, n_th, n_phi)) |
| 301 |
9/18✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 11 times.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 11 times.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 11 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 11 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 11 times.
|
44 | allocate (hphi(n_r, n_th, n_phi)) |
| 302 |
9/18✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 11 times.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 11 times.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 11 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 11 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 11 times.
|
44 | allocate (Bmod_arr(n_r, n_th, n_phi)) |
| 303 | |||
| 304 | 11 | is_identity_scaling = .false. | |
| 305 | 11 | is_sqrt_s_scaling = .false. | |
| 306 | select type (scaling_ptr) | ||
| 307 | type is (identity_scaling_t) | ||
| 308 | 11 | is_identity_scaling = .true. | |
| 309 | type is (sqrt_s_scaling_t) | ||
| 310 | 11 | is_sqrt_s_scaling = .true. | |
| 311 | class default | ||
| 312 | ✗ | error stop 'sample_spline_arrays: Unsupported scaling type' | |
| 313 | end select | ||
| 314 | |||
| 315 | 11 | i_ctr = 0 | |
| 316 | |||
| 317 | #ifdef SIMPLE_NVHPC | ||
| 318 | do i_phi = 1, n_phi | ||
| 319 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 320 | i_ctr = i_ctr + 1 | ||
| 321 | call print_progress(i_ctr, n_phi) | ||
| 322 | #endif | ||
| 323 | do i_th = 1, n_th | ||
| 324 | do i_r = 1, n_r | ||
| 325 | x_scaled(1) = xmin(1) + h_r*dble(i_r - 1) | ||
| 326 | x_scaled(2) = xmin(2) + h_th*dble(i_th - 1) | ||
| 327 | x_scaled(3) = xmin(3) + h_phi*dble(i_phi - 1) | ||
| 328 | |||
| 329 | scaling_jac = 1d0 | ||
| 330 | if (is_identity_scaling) then | ||
| 331 | x_ref = x_scaled | ||
| 332 | else if (is_sqrt_s_scaling) then | ||
| 333 | x_ref(1) = x_scaled(1)**2 | ||
| 334 | x_ref(2) = x_scaled(2) | ||
| 335 | x_ref(3) = x_scaled(3) | ||
| 336 | scaling_jac(1) = 2d0*x_scaled(1) | ||
| 337 | end if | ||
| 338 | |||
| 339 | call evaluate_at_ref_coords(source, ref_coords, x_scaled, & | ||
| 340 | x_ref, x_cart, e_cov, Acov, & | ||
| 341 | hcov, Bmod, needs_scaling) | ||
| 342 | if (needs_scaling) then | ||
| 343 | Acov(1) = Acov(1)*scaling_jac(1) | ||
| 344 | hcov(1) = hcov(1)*scaling_jac(1) | ||
| 345 | end if | ||
| 346 | |||
| 347 | Ar(i_r, i_th, i_phi) = Acov(1) | ||
| 348 | Ath(i_r, i_th, i_phi) = Acov(2) | ||
| 349 | Aphi(i_r, i_th, i_phi) = Acov(3) | ||
| 350 | |||
| 351 | hr(i_r, i_th, i_phi) = hcov(1) | ||
| 352 | hth(i_r, i_th, i_phi) = hcov(2) | ||
| 353 | hphi(i_r, i_th, i_phi) = hcov(3) | ||
| 354 | |||
| 355 | Bmod_arr(i_r, i_th, i_phi) = Bmod | ||
| 356 | end do | ||
| 357 | end do | ||
| 358 | end do | ||
| 359 | #else | ||
| 360 | !$omp parallel private(i_r, i_th, i_phi, x_scaled, x_ref, x_cart, & | ||
| 361 | 11 | !$omp& e_cov, scaling_jac, Bmod, Acov, hcov, needs_scaling) | |
| 362 | !$omp do | ||
| 363 | do i_phi = 1, n_phi | ||
| 364 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 365 | !$omp atomic | ||
| 366 | i_ctr = i_ctr + 1 | ||
| 367 | call print_progress(i_ctr, n_phi) | ||
| 368 | #endif | ||
| 369 | do i_th = 1, n_th | ||
| 370 | do i_r = 1, n_r | ||
| 371 | x_scaled(1) = xmin(1) + h_r*dble(i_r - 1) | ||
| 372 | x_scaled(2) = xmin(2) + h_th*dble(i_th - 1) | ||
| 373 | x_scaled(3) = xmin(3) + h_phi*dble(i_phi - 1) | ||
| 374 | |||
| 375 | scaling_jac = 1d0 | ||
| 376 | if (is_identity_scaling) then | ||
| 377 | x_ref = x_scaled | ||
| 378 | else if (is_sqrt_s_scaling) then | ||
| 379 | x_ref(1) = x_scaled(1)**2 | ||
| 380 | x_ref(2) = x_scaled(2) | ||
| 381 | x_ref(3) = x_scaled(3) | ||
| 382 | scaling_jac(1) = 2d0*x_scaled(1) | ||
| 383 | end if | ||
| 384 | |||
| 385 | call evaluate_at_ref_coords(source, ref_coords, x_scaled, & | ||
| 386 | x_ref, x_cart, e_cov, Acov, & | ||
| 387 | hcov, Bmod, needs_scaling) | ||
| 388 | if (needs_scaling) then | ||
| 389 | Acov(1) = Acov(1)*scaling_jac(1) | ||
| 390 | hcov(1) = hcov(1)*scaling_jac(1) | ||
| 391 | end if | ||
| 392 | |||
| 393 | Ar(i_r, i_th, i_phi) = Acov(1) | ||
| 394 | Ath(i_r, i_th, i_phi) = Acov(2) | ||
| 395 | Aphi(i_r, i_th, i_phi) = Acov(3) | ||
| 396 | |||
| 397 | hr(i_r, i_th, i_phi) = hcov(1) | ||
| 398 | hth(i_r, i_th, i_phi) = hcov(2) | ||
| 399 | hphi(i_r, i_th, i_phi) = hcov(3) | ||
| 400 | |||
| 401 | Bmod_arr(i_r, i_th, i_phi) = Bmod | ||
| 402 | end do | ||
| 403 | end do | ||
| 404 | end do | ||
| 405 | !$omp end do | ||
| 406 | !$omp end parallel | ||
| 407 | #endif | ||
| 408 | |||
| 409 |
6/6✓ Branch 0 taken 372 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 15220 times.
✓ Branch 3 taken 372 times.
✓ Branch 4 taken 645316 times.
✓ Branch 5 taken 15220 times.
|
660919 | Ar = Ar - Ar(1, 1, 1) |
| 410 |
6/6✓ Branch 0 taken 372 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 15220 times.
✓ Branch 3 taken 372 times.
✓ Branch 4 taken 645316 times.
✓ Branch 5 taken 15220 times.
|
660919 | Ath = Ath - Ath(1, 1, 1) |
| 411 |
6/6✓ Branch 0 taken 372 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 15220 times.
✓ Branch 3 taken 372 times.
✓ Branch 4 taken 645316 times.
✓ Branch 5 taken 15220 times.
|
660919 | Aphi = Aphi - Aphi(1, 1, 1) |
| 412 | 11 | end subroutine sample_spline_arrays | |
| 413 | |||
| 414 | 645316 | subroutine evaluate_at_ref_coords(source, ref_coords, x_scaled, x_ref, x_cart, & | |
| 415 | e_cov, Acov, hcov, Bmod, needs_scaling) | ||
| 416 | !> Evaluate source field for spline sampling. | ||
| 417 | !> | ||
| 418 | !> Supported source modes: | ||
| 419 | !> 1. VMEC source: evaluate at x_ref = (s, theta, phi) and return covariant | ||
| 420 | !> components in reference coordinates. Caller must apply scaling. | ||
| 421 | !> 2. Cartesian source: evaluate at x_cart and transform to reference | ||
| 422 | !> coordinates (s, theta, phi). Caller must apply scaling. | ||
| 423 | !> | ||
| 424 | !> All sources now return components in reference coordinates (s, theta, phi). | ||
| 425 | !> Caller applies coordinate scaling to obtain components in scaled coords. | ||
| 426 | class(magnetic_field_t), intent(in) :: source | ||
| 427 | class(coordinate_system_t), intent(in) :: ref_coords | ||
| 428 | real(dp), intent(in) :: x_scaled(3) | ||
| 429 | real(dp), intent(in) :: x_ref(3) | ||
| 430 | real(dp), intent(out) :: x_cart(3), e_cov(3, 3) | ||
| 431 | real(dp), intent(out) :: Acov(3), hcov(3), Bmod | ||
| 432 | logical, intent(out) :: needs_scaling | ||
| 433 | |||
| 434 | real(dp) :: A_cart(3), h_cart(3) | ||
| 435 | real(dp) :: x_src(3) | ||
| 436 | |||
| 437 | select type (source) | ||
| 438 | type is (vmec_field_t) | ||
| 439 | select type (ref_coords) | ||
| 440 | type is (vmec_coordinate_system_t) | ||
| 441 | ! Same coordinate system: evaluate directly in VMEC coordinates. | ||
| 442 | 499968 | call source%evaluate(x_ref, Acov, hcov, Bmod) | |
| 443 | 499968 | x_cart = 0d0 | |
| 444 | 499968 | e_cov = 0d0 | |
| 445 | 499968 | needs_scaling = .true. | |
| 446 | 499968 | return | |
| 447 | class default | ||
| 448 | ! Different reference coordinates: evaluate VMEC field at the same | ||
| 449 | ! physical point and transform covariant components into ref_coords. | ||
| 450 | 74052 | call ref_coords%evaluate_cart(x_ref, x_cart) | |
| 451 | call find_vmec_coords_for_cart_point(source%coords, ref_coords, & | ||
| 452 | 74052 | x_ref, x_cart, x_src) | |
| 453 | |||
| 454 | 74052 | call source%evaluate(x_src, Acov, hcov, Bmod) | |
| 455 | |||
| 456 | 74052 | call source%coords%cov_to_cart(x_src, Acov, A_cart) | |
| 457 | 74052 | call source%coords%cov_to_cart(x_src, hcov, h_cart) | |
| 458 | |||
| 459 | 74052 | call ref_coords%covariant_basis(x_ref, e_cov) | |
| 460 |
4/4✓ Branch 0 taken 222156 times.
✓ Branch 1 taken 74052 times.
✓ Branch 2 taken 666468 times.
✓ Branch 3 taken 222156 times.
|
962676 | Acov = matmul(A_cart, e_cov) |
| 461 |
4/4✓ Branch 0 taken 222156 times.
✓ Branch 1 taken 74052 times.
✓ Branch 2 taken 666468 times.
✓ Branch 3 taken 222156 times.
|
962676 | hcov = matmul(h_cart, e_cov) |
| 462 | end select | ||
| 463 | 74052 | needs_scaling = .true. | |
| 464 | 74052 | return | |
| 465 | type is (boozer_chartmap_field_t) | ||
| 466 | ! Boozer chartmap evaluates directly in (rho, theta_B, phi_B), | ||
| 467 | ! which matches the chartmap reference coordinate system. | ||
| 468 | ✗ | call source%evaluate(x_ref, Acov, hcov, Bmod) | |
| 469 | ✗ | x_cart = 0d0 | |
| 470 | ✗ | e_cov = 0d0 | |
| 471 | ✗ | needs_scaling = .true. | |
| 472 | ✗ | return | |
| 473 | class default | ||
| 474 | continue | ||
| 475 | end select | ||
| 476 | |||
| 477 | select type (coords => source%coords) | ||
| 478 | type is (cartesian_coordinate_system_t) | ||
| 479 | 71296 | call ref_coords%evaluate_cart(x_ref, x_cart) | |
| 480 | 71296 | call ref_coords%covariant_basis(x_ref, e_cov) | |
| 481 | 71296 | call source%evaluate(x_cart, A_cart, h_cart, Bmod) | |
| 482 |
4/4✓ Branch 0 taken 213888 times.
✓ Branch 1 taken 71296 times.
✓ Branch 2 taken 641664 times.
✓ Branch 3 taken 213888 times.
|
926848 | Acov = matmul(A_cart, e_cov) |
| 483 |
4/4✓ Branch 0 taken 213888 times.
✓ Branch 1 taken 71296 times.
✓ Branch 2 taken 641664 times.
✓ Branch 3 taken 213888 times.
|
926848 | hcov = matmul(h_cart, e_cov) |
| 484 | 71296 | needs_scaling = .true. | |
| 485 | class default | ||
| 486 | error stop & | ||
| 487 | ✗ | "evaluate_at_ref_coords: unsupported source coordinate system" | |
| 488 | end select | ||
| 489 | end subroutine evaluate_at_ref_coords | ||
| 490 | |||
| 491 | 74052 | subroutine find_vmec_coords_for_cart_point(vmec_coords, ref_coords, u_ref, & | |
| 492 | x_target, u_vmec) | ||
| 493 | use util, only: twopi | ||
| 494 | !> Find VMEC coordinates u_vmec such that | ||
| 495 | !> vmec_coords%evaluate_cart(u_vmec) == x_target. | ||
| 496 | !> Prefer the native libneo inversion (from_cyl) to avoid fragile | ||
| 497 | !> finite-difference Newton iterations on periodic variables. | ||
| 498 | class(coordinate_system_t), intent(in) :: vmec_coords | ||
| 499 | class(coordinate_system_t), intent(in) :: ref_coords | ||
| 500 | real(dp), intent(in) :: u_ref(3) | ||
| 501 | real(dp), intent(in) :: x_target(3) | ||
| 502 | real(dp), intent(out) :: u_vmec(3) | ||
| 503 | |||
| 504 | real(dp) :: xcyl(3) | ||
| 505 | integer :: ierr | ||
| 506 | |||
| 507 | 74052 | associate (dummy => ref_coords) | |
| 508 | end associate | ||
| 509 | |||
| 510 | 74052 | xcyl(1) = sqrt(x_target(1)**2 + x_target(2)**2) | |
| 511 |
3/4✓ Branch 0 taken 71808 times.
✓ Branch 1 taken 2244 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 71808 times.
|
74052 | xcyl(2) = modulo(atan2(x_target(2), x_target(1)), twopi) |
| 512 | 74052 | xcyl(3) = x_target(3) | |
| 513 | |||
| 514 | select type (vcs => vmec_coords) | ||
| 515 | type is (vmec_coordinate_system_t) | ||
| 516 | 74052 | call vcs%from_cyl(xcyl, u_vmec, ierr) | |
| 517 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 74052 times.
|
74052 | if (ierr /= 0) then |
| 518 | ✗ | print *, "find_vmec_coords_for_cart_point: vmec from_cyl failed" | |
| 519 | ✗ | print *, " ierr = ", ierr | |
| 520 | ✗ | print *, " u_ref = ", u_ref | |
| 521 | ✗ | print *, " xcyl = ", xcyl | |
| 522 | ✗ | print *, " x_target = ", x_target | |
| 523 | ✗ | error stop | |
| 524 | end if | ||
| 525 | class default | ||
| 526 | ✗ | error stop "find_vmec_coords_for_cart_point: expected VMEC coordinates" | |
| 527 | end select | ||
| 528 | 74052 | end subroutine find_vmec_coords_for_cart_point | |
| 529 | |||
| 530 | 19 | subroutine splined_field_cleanup(self) | |
| 531 | type(splined_field_t), intent(inout) :: self | ||
| 532 | |||
| 533 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 12 times.
|
19 | if (self%initialized) then |
| 534 | 7 | call destroy_batch_splines_3d(self%splines) | |
| 535 | 7 | self%initialized = .false. | |
| 536 | end if | ||
| 537 | 19 | end subroutine splined_field_cleanup | |
| 538 | |||
| 539 | 11 | function get_default_scaling(ref_coords, sqrt_s_scaling, identity_scaling) & | |
| 540 | result(scaling_ptr) | ||
| 541 | !> Select default scaling based on coordinate system rho_convention. | ||
| 542 | !> RHO_TOR/RHO_POL: identity (coords already in rho = sqrt(s)) | ||
| 543 | !> PSI_TOR_NORM/PSI_POL_NORM: sqrt_s (coords in s, need to transform) | ||
| 544 | !> UNKNOWN/other: identity (rho is a geometric radius-like coordinate) | ||
| 545 | class(coordinate_system_t), intent(in) :: ref_coords | ||
| 546 | type(sqrt_s_scaling_t), intent(in), target :: sqrt_s_scaling | ||
| 547 | type(identity_scaling_t), intent(in), target :: identity_scaling | ||
| 548 | class(coordinate_scaling_t), pointer :: scaling_ptr | ||
| 549 | |||
| 550 | select type (cs => ref_coords) | ||
| 551 | class is (chartmap_coordinate_system_t) | ||
| 552 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if (cs%rho_convention == PSI_TOR_NORM .or. & |
| 553 | cs%rho_convention == PSI_POL_NORM) then | ||
| 554 | scaling_ptr => sqrt_s_scaling | ||
| 555 | else | ||
| 556 | 11 | scaling_ptr => identity_scaling | |
| 557 | end if | ||
| 558 | class default | ||
| 559 | scaling_ptr => sqrt_s_scaling | ||
| 560 | end select | ||
| 561 | 11 | end function get_default_scaling | |
| 562 | |||
| 563 |
16/24__field_splined_MOD___copy_field_splined_Splined_field_t:
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 3 times.
✓ Branch 5 taken 5 times.
✓ Branch 6 taken 3 times.
__field_splined_MOD___final_field_splined_Splined_field_t:
✗ Branch 0 not taken.
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 19 times.
✓ Branch 8 taken 19 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 19 times.
✓ Branch 11 taken 19 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 8 times.
✓ Branch 14 taken 11 times.
✓ Branch 15 taken 6 times.
✓ Branch 16 taken 2 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 19 times.
|
65 | end module field_splined |
| 564 |