| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module field_can_meiss | ||
| 2 | !> Meiss canonical coordinates for magnetic field evaluation. | ||
| 3 | !> | ||
| 4 | !> This module implements Meiss-type canonical coordinates (r, theta_c, phi_c) | ||
| 5 | !> where the magnetic field has the simple covariant form h = (0, hth, hph). | ||
| 6 | !> The transformation to canonical coordinates is computed via ODE integration. | ||
| 7 | !> | ||
| 8 | !> Coordinate conventions: | ||
| 9 | !> - r: Integrator radial coordinate (= sqrt(s) by default) | ||
| 10 | !> - theta_c, phi_c: Canonical poloidal and toroidal angles | ||
| 11 | !> - s: Reference flux coordinate (normalized toroidal flux) | ||
| 12 | !> | ||
| 13 | !> Field source handling: | ||
| 14 | !> - vmec_field_t: Expects s-coordinates, returns components in s-coords. | ||
| 15 | !> This module converts r->s and applies Jacobian ds/dr = 2r. | ||
| 16 | !> - splined_field_t and others: Expected to operate in r-coordinates | ||
| 17 | !> directly, no coordinate transformation needed. | ||
| 18 | !> | ||
| 19 | !> The coordinate_scaling abstraction (default: sqrt_s_scaling_t) handles | ||
| 20 | !> the s <-> r transformation systematically. | ||
| 21 | !> | ||
| 22 | !> Key subroutines: | ||
| 23 | !> - init_meiss: Initialize with field source and grid parameters | ||
| 24 | !> - get_meiss_coordinates: Build canonical coordinate transformation | ||
| 25 | !> - evaluate_meiss: Fast splined field evaluation in canonical coords | ||
| 26 | !> - cleanup_meiss: Release resources | ||
| 27 | !> | ||
| 28 | !> Related modules: | ||
| 29 | !> - field_can_albert: Extends Meiss with psi-based radial coordinate | ||
| 30 | !> - coordinate_scaling: Abstraction for s <-> r transforms | ||
| 31 | |||
| 32 | use, intrinsic :: iso_fortran_env, only : dp => real64 | ||
| 33 | use interpolate, only : & | ||
| 34 | BatchSplineData3D, construct_batch_splines_3d, & | ||
| 35 | evaluate_batch_splines_3d, evaluate_batch_splines_3d_der, & | ||
| 36 | evaluate_batch_splines_3d_der2 | ||
| 37 | use util, only : twopi | ||
| 38 | use field_can_base, only: count_field_evaluation, field_can_t | ||
| 39 | use field, only : magnetic_field_t, vmec_field_t, field_clone | ||
| 40 | use field_splined, only: splined_field_t | ||
| 41 | use field_spectre, only: spectre_field_t | ||
| 42 | use coordinate_scaling, only : coordinate_scaling_t, sqrt_s_scaling_t, & | ||
| 43 | identity_scaling_t, coordinate_scaling_clone | ||
| 44 | use libneo_coordinates, only: chartmap_coordinate_system_t, UNKNOWN, & | ||
| 45 | RHO_TOR, RHO_POL, PSI_TOR_NORM, PSI_POL_NORM | ||
| 46 | |||
| 47 | implicit none | ||
| 48 | |||
| 49 | type :: grid_indices_t | ||
| 50 | integer :: i_th, i_phi | ||
| 51 | end type grid_indices_t | ||
| 52 | |||
| 53 | class(magnetic_field_t), allocatable :: field_noncan | ||
| 54 | integer, parameter :: default_n_r = 62, default_n_th = 63, default_n_phi = 64 | ||
| 55 | real(dp), parameter :: default_transformation_relerr = 1d-11 | ||
| 56 | real(dp), parameter :: default_xmin(3) = [1d-6, 0d0, 0d0] | ||
| 57 | integer :: n_r=default_n_r, n_th=default_n_th, n_phi=default_n_phi | ||
| 58 | real(dp) :: transformation_relerr = default_transformation_relerr | ||
| 59 | real(dp) :: xmin(3) = default_xmin | ||
| 60 | real(dp) :: xmax(3) = [1d0, twopi, twopi] | ||
| 61 | |||
| 62 | real(dp) :: h_r, h_th, h_phi | ||
| 63 | |||
| 64 | ! Batch spline for optimized field evaluation (5 components: Ath, Aph, hth, hph, Bmod) | ||
| 65 | type(BatchSplineData3D) :: spl_field_batch | ||
| 66 | logical :: batch_splines_initialized = .false. | ||
| 67 | |||
| 68 | ! For splining lambda (difference between canonical and toroidal cylinder angle) | ||
| 69 | ! and chi (gauge transformation) | ||
| 70 | real(dp), dimension(:,:,:), allocatable :: lam_phi, chi_gauge | ||
| 71 | ! Batch spline for transformation components (2 components: lam_phi, chi_gauge) | ||
| 72 | type(BatchSplineData3D) :: spl_transform_batch | ||
| 73 | logical :: transform_splines_initialized = .false. | ||
| 74 | |||
| 75 | integer, parameter :: default_spline_order = 5 | ||
| 76 | integer, parameter :: default_ode_step_limit = 1000000 | ||
| 77 | integer :: order(3) = [default_spline_order, default_spline_order, & | ||
| 78 | default_spline_order] | ||
| 79 | integer :: ode_step_limit = default_ode_step_limit | ||
| 80 | logical, parameter :: periodic(3) = [.False., .True., .True.] | ||
| 81 | |||
| 82 | ! Coordinate scaling for s <-> r transform (default: r = sqrt(s)) | ||
| 83 | class(coordinate_scaling_t), allocatable :: coord_scaling | ||
| 84 | |||
| 85 | ! Caller-owned harvest of one Meiss construction. field_can_spectre keeps one per | ||
| 86 | ! SPECTRE volume; harvest_meiss_volume moves the module-level splines into a slot. | ||
| 87 | type :: meiss_volume_t | ||
| 88 | type(BatchSplineData3D) :: spl_field | ||
| 89 | type(BatchSplineData3D) :: spl_transform | ||
| 90 | class(coordinate_scaling_t), allocatable :: scaling | ||
| 91 | real(dp) :: rmin | ||
| 92 | real(dp) :: rmax | ||
| 93 | end type meiss_volume_t | ||
| 94 | |||
| 95 | contains | ||
| 96 | |||
| 97 |
32/38✓ Branch 0 taken 32 times.
✓ Branch 1 taken 96 times.
✓ Branch 2 taken 32 times.
✓ Branch 3 taken 96 times.
✓ Branch 4 taken 32 times.
✓ Branch 5 taken 96 times.
✓ Branch 6 taken 32 times.
✓ Branch 7 taken 96 times.
✓ Branch 8 taken 32 times.
✓ Branch 9 taken 96 times.
✓ Branch 10 taken 32 times.
✓ Branch 11 taken 96 times.
✓ Branch 12 taken 32 times.
✓ Branch 13 taken 96 times.
✓ Branch 14 taken 32 times.
✓ Branch 15 taken 96 times.
✓ Branch 16 taken 32 times.
✓ Branch 17 taken 96 times.
✓ Branch 18 taken 32 times.
✓ Branch 19 taken 96 times.
✓ Branch 20 taken 32 times.
✓ Branch 21 taken 96 times.
✓ Branch 22 taken 32 times.
✓ Branch 23 taken 96 times.
✓ Branch 24 taken 32 times.
✓ Branch 25 taken 96 times.
✓ Branch 26 taken 96 times.
✓ Branch 27 taken 32 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 32 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 32 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 32 times.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✓ Branch 38 taken 32 times.
|
1792 | subroutine harvest_meiss_volume(vol) |
| 98 | !> Move the current module-level construction (field and transform batch | ||
| 99 | !> splines, radial domain, scaling) into a caller-owned slot. The batch spline | ||
| 100 | !> coefficient arrays are moved, not copied; the module init flags are cleared | ||
| 101 | !> so a following cleanup_meiss or construction re-uses the freed slots. | ||
| 102 | type(meiss_volume_t), intent(out) :: vol | ||
| 103 | |||
| 104 | 32 | call move_batch_spline_3d(spl_field_batch, vol%spl_field) | |
| 105 | 32 | call move_batch_spline_3d(spl_transform_batch, vol%spl_transform) | |
| 106 | 32 | batch_splines_initialized = .false. | |
| 107 | 32 | transform_splines_initialized = .false. | |
| 108 | |||
| 109 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
32 | call coordinate_scaling_clone(coord_scaling, vol%scaling) |
| 110 | 32 | vol%rmin = xmin(1) | |
| 111 | 32 | vol%rmax = xmax(1) | |
| 112 | 32 | end subroutine harvest_meiss_volume | |
| 113 | |||
| 114 | |||
| 115 |
15/16✓ Branch 0 taken 64 times.
✓ Branch 1 taken 192 times.
✓ Branch 2 taken 64 times.
✓ Branch 3 taken 192 times.
✓ Branch 4 taken 64 times.
✓ Branch 5 taken 192 times.
✓ Branch 6 taken 64 times.
✓ Branch 7 taken 192 times.
✓ Branch 8 taken 64 times.
✓ Branch 9 taken 192 times.
✓ Branch 10 taken 64 times.
✓ Branch 11 taken 192 times.
✓ Branch 12 taken 192 times.
✓ Branch 13 taken 64 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 64 times.
|
1792 | subroutine move_batch_spline_3d(src, dst) |
| 116 | type(BatchSplineData3D), intent(inout) :: src | ||
| 117 | type(BatchSplineData3D), intent(out) :: dst | ||
| 118 | |||
| 119 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 192 times.
|
256 | dst%order = src%order |
| 120 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 192 times.
|
256 | dst%num_points = src%num_points |
| 121 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 192 times.
|
256 | dst%periodic = src%periodic |
| 122 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 192 times.
|
256 | dst%h_step = src%h_step |
| 123 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 192 times.
|
256 | dst%x_min = src%x_min |
| 124 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 192 times.
|
256 | dst%inv_h_step = src%inv_h_step |
| 125 |
2/2✓ Branch 0 taken 192 times.
✓ Branch 1 taken 64 times.
|
256 | dst%period = src%period |
| 126 | 64 | dst%num_quantities = src%num_quantities | |
| 127 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 64 times.
|
64 | call move_alloc(src%coeff, dst%coeff) |
| 128 | 64 | end subroutine move_batch_spline_3d | |
| 129 | |||
| 130 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 32697907 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 32697907 times.
|
32697907 | subroutine rh_can_wrapper(r_c, z, dz, context) |
| 131 | real(dp), intent(in) :: r_c | ||
| 132 | real(dp), dimension(:), intent(in) :: z | ||
| 133 | real(dp), dimension(:), intent(out) :: dz | ||
| 134 | class(*), intent(in) :: context | ||
| 135 | |||
| 136 | integer :: i_th, i_phi | ||
| 137 | |||
| 138 | select type(context) | ||
| 139 | type is (grid_indices_t) | ||
| 140 | 32697907 | i_th = context%i_th | |
| 141 | 32697907 | i_phi = context%i_phi | |
| 142 | end select | ||
| 143 | |||
| 144 |
7/16✗ Branch 0 not taken.
✓ Branch 1 taken 32697907 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 32697907 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 32697907 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 32697907 times.
✓ Branch 13 taken 32697907 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 65395814 times.
✓ Branch 16 taken 32697907 times.
|
130791628 | call rh_can(r_c, z, dz, i_th, i_phi) |
| 145 | 32697907 | end subroutine rh_can_wrapper | |
| 146 | |||
| 147 | 92 | subroutine init_meiss(field_noncan_, n_r_, n_th_, n_phi_, rmin, rmax, thmin, thmax, & | |
| 148 | scaling, transformation_relerr_, spline_order_, ode_step_limit_) | ||
| 149 | use finite_check, only: is_finite | ||
| 150 | use, intrinsic :: iso_fortran_env, only: int64 | ||
| 151 | use new_vmec_stuff_mod, only : nper | ||
| 152 | |||
| 153 | class(magnetic_field_t), intent(in) :: field_noncan_ | ||
| 154 | integer, intent(in), optional :: n_r_, n_th_, n_phi_ | ||
| 155 | real(dp), intent(in), optional :: rmin, rmax, thmin, thmax | ||
| 156 | class(coordinate_scaling_t), intent(in), optional :: scaling | ||
| 157 | real(dp), intent(in), optional :: transformation_relerr_ | ||
| 158 | integer, intent(in), optional :: spline_order_ | ||
| 159 | integer, intent(in), optional :: ode_step_limit_ | ||
| 160 | |||
| 161 | 46 | call cleanup_meiss | |
| 162 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
46 | call field_clone(field_noncan_, field_noncan) |
| 163 | |||
| 164 | ! Initialize coordinate scaling based on field type | ||
| 165 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
46 | if (allocated(coord_scaling)) deallocate(coord_scaling) |
| 166 |
3/4✓ Branch 0 taken 33 times.
✓ Branch 1 taken 13 times.
✓ Branch 2 taken 33 times.
✗ Branch 3 not taken.
|
46 | if (present(scaling)) then |
| 167 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 33 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
33 | call coordinate_scaling_clone(scaling, coord_scaling) |
| 168 | else | ||
| 169 | ! Choose scaling based on field's coordinate system: | ||
| 170 | ! - VMEC field: uses s-coordinates, needs sqrt_s_scaling (r = sqrt(s)) | ||
| 171 | ! - Splined field with RHO_TOR/RHO_POL: uses rho = sqrt(s), needs sqrt_s_scaling | ||
| 172 | ! - Splined field with UNKNOWN/geometric: uses geometric rho, needs identity_scaling | ||
| 173 | ! - Splined field with PSI_TOR_NORM: uses s-coordinates, needs sqrt_s_scaling | ||
| 174 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
13 | call choose_default_scaling(field_noncan_, coord_scaling) |
| 175 | end if | ||
| 176 | |||
| 177 | 46 | n_r = default_n_r | |
| 178 | 46 | n_th = default_n_th | |
| 179 | 46 | n_phi = default_n_phi | |
| 180 | 46 | transformation_relerr = default_transformation_relerr | |
| 181 |
2/2✓ Branch 0 taken 138 times.
✓ Branch 1 taken 46 times.
|
184 | order = default_spline_order |
| 182 | 46 | ode_step_limit = default_ode_step_limit | |
| 183 | 46 | xmin = default_xmin | |
| 184 |
2/2✓ Branch 0 taken 138 times.
✓ Branch 1 taken 46 times.
|
184 | xmax = [1d0, twopi, twopi/nper] |
| 185 |
2/2✓ Branch 0 taken 43 times.
✓ Branch 1 taken 3 times.
|
46 | if (present(n_r_)) n_r = n_r_ |
| 186 |
2/2✓ Branch 0 taken 43 times.
✓ Branch 1 taken 3 times.
|
46 | if (present(n_th_)) n_th = n_th_ |
| 187 |
2/2✓ Branch 0 taken 43 times.
✓ Branch 1 taken 3 times.
|
46 | if (present(n_phi_)) n_phi = n_phi_ |
| 188 |
4/4✓ Branch 0 taken 33 times.
✓ Branch 1 taken 13 times.
✓ Branch 2 taken 99 times.
✓ Branch 3 taken 33 times.
|
145 | if (present(spline_order_)) order = spline_order_ |
| 189 |
2/2✓ Branch 0 taken 33 times.
✓ Branch 1 taken 13 times.
|
46 | if (present(ode_step_limit_)) ode_step_limit = ode_step_limit_ |
| 190 |
7/10✓ Branch 0 taken 138 times.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 138 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 138 times.
✓ Branch 5 taken 46 times.
✓ Branch 6 taken 138 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 46 times.
|
368 | if (any(order < 3) .or. any(order > 5)) then |
| 191 | ✗ | error stop 'canonical map spline order must be between 3 and 5' | |
| 192 | end if | ||
| 193 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
|
46 | if (ode_step_limit < 1) error stop 'canonical map ODE step limit must be positive' |
| 194 |
4/6✓ Branch 0 taken 138 times.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 138 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 46 times.
|
184 | if (any([n_r, n_th, n_phi] < order + 1)) then |
| 195 | ✗ | error stop 'canonical map grid dimensions must exceed the spline order' | |
| 196 | end if | ||
| 197 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
|
46 | if (max(n_r, n_th, n_phi) > 512) then |
| 198 | ✗ | error stop 'canonical map grid dimensions must not exceed 512' | |
| 199 | end if | ||
| 200 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
|
46 | if (int(n_r, int64)*int(n_th, int64)*int(n_phi, int64) > 2097152_int64) then |
| 201 | ✗ | error stop 'canonical map grid exceeds the 2097152-point limit' | |
| 202 | end if | ||
| 203 |
2/2✓ Branch 0 taken 41 times.
✓ Branch 1 taken 5 times.
|
46 | if (present(transformation_relerr_)) then |
| 204 |
2/4✓ Branch 0 taken 41 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 41 times.
|
41 | if (.not. is_finite(transformation_relerr_) .or. & |
| 205 | transformation_relerr_ <= 0d0) then | ||
| 206 | ✗ | error stop 'canonical map ODE tolerance must be finite and positive' | |
| 207 | end if | ||
| 208 | 41 | transformation_relerr = transformation_relerr_ | |
| 209 | end if | ||
| 210 | |||
| 211 |
2/2✓ Branch 0 taken 36 times.
✓ Branch 1 taken 10 times.
|
46 | if (present(rmin)) xmin(1) = rmin |
| 212 |
2/2✓ Branch 0 taken 36 times.
✓ Branch 1 taken 10 times.
|
46 | if (present(rmax)) xmax(1) = rmax |
| 213 |
2/2✓ Branch 0 taken 34 times.
✓ Branch 1 taken 12 times.
|
46 | if (present(thmin)) xmin(2) = thmin |
| 214 |
2/2✓ Branch 0 taken 34 times.
✓ Branch 1 taken 12 times.
|
46 | if (present(thmax)) xmax(2) = thmax |
| 215 | 46 | h_r = (xmax(1)-xmin(1))/(n_r-1) | |
| 216 | 46 | h_th = (xmax(2)-xmin(2))/(n_th-1) | |
| 217 | 46 | h_phi = (xmax(3)-xmin(3))/(n_phi-1) | |
| 218 | 46 | end subroutine init_meiss | |
| 219 | |||
| 220 | |||
| 221 | 14 | subroutine choose_default_scaling(field, scaling) | |
| 222 | !> Choose appropriate coordinate scaling based on field type. | ||
| 223 | !> | ||
| 224 | !> Scaling determines how the Meiss integrator coordinate r relates to flux: | ||
| 225 | !> - sqrt_s_scaling: r = sqrt(s), use when field is in s-coordinates (VMEC) | ||
| 226 | !> - identity_scaling: r = rho directly, use when field is already in rho coords | ||
| 227 | !> | ||
| 228 | !> For VMEC source field: uses sqrt_s_scaling (convert s -> r = sqrt(s)) | ||
| 229 | !> For RHO_TOR/RHO_POL chartmap: uses identity (rho IS sqrt(s)) | ||
| 230 | !> For PSI_TOR_NORM/PSI_POL_NORM: uses sqrt_s (convert s -> r = sqrt(s)) | ||
| 231 | !> For UNKNOWN (geometric): uses identity (rho is geometric) | ||
| 232 | class(magnetic_field_t), intent(in) :: field | ||
| 233 | class(coordinate_scaling_t), allocatable, intent(out) :: scaling | ||
| 234 | |||
| 235 | logical :: use_identity | ||
| 236 | |||
| 237 | 14 | use_identity = .false. | |
| 238 | |||
| 239 | select type (fld => field) | ||
| 240 | type is (splined_field_t) | ||
| 241 | 4 | select type (cs => fld%coords) | |
| 242 | class is (chartmap_coordinate_system_t) | ||
| 243 | ! RHO_TOR/RHO_POL: rho = sqrt(s), already in Meiss integrator coords | ||
| 244 | ! UNKNOWN: geometric rho, no flux meaning, use directly | ||
| 245 | ! PSI_TOR_NORM/PSI_POL_NORM: rho = s, needs sqrt transform | ||
| 246 | if (cs%rho_convention == RHO_TOR .or. & | ||
| 247 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | cs%rho_convention == RHO_POL .or. & |
| 248 | cs%rho_convention == UNKNOWN) then | ||
| 249 | use_identity = .true. | ||
| 250 | end if | ||
| 251 | end select | ||
| 252 | end select | ||
| 253 | |||
| 254 | if (use_identity) then | ||
| 255 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
|
4 | allocate(identity_scaling_t :: scaling) |
| 256 | else | ||
| 257 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
|
10 | allocate(sqrt_s_scaling_t :: scaling) |
| 258 | end if | ||
| 259 | 14 | end subroutine choose_default_scaling | |
| 260 | |||
| 261 | |||
| 262 | 141 | subroutine cleanup_meiss() | |
| 263 | ! Clean up batch splines to prevent memory leaks | ||
| 264 | use interpolate, only: destroy_batch_splines_3d | ||
| 265 | |||
| 266 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 133 times.
|
141 | if (batch_splines_initialized) then |
| 267 | 8 | call destroy_batch_splines_3d(spl_field_batch) | |
| 268 | 8 | batch_splines_initialized = .false. | |
| 269 | end if | ||
| 270 | |||
| 271 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 133 times.
|
141 | if (transform_splines_initialized) then |
| 272 | 8 | call destroy_batch_splines_3d(spl_transform_batch) | |
| 273 | 8 | transform_splines_initialized = .false. | |
| 274 | end if | ||
| 275 | |||
| 276 | ! Clean up allocated arrays | ||
| 277 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 140 times.
|
141 | if (allocated(lam_phi)) deallocate(lam_phi) |
| 278 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 140 times.
|
141 | if (allocated(chi_gauge)) deallocate(chi_gauge) |
| 279 |
3/4✓ Branch 0 taken 41 times.
✓ Branch 1 taken 100 times.
✓ Branch 2 taken 41 times.
✗ Branch 3 not taken.
|
182 | if (allocated(field_noncan)) deallocate(field_noncan) |
| 280 |
3/4✓ Branch 0 taken 41 times.
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 41 times.
|
141 | if (allocated(coord_scaling)) deallocate(coord_scaling) |
| 281 | 141 | end subroutine cleanup_meiss | |
| 282 | |||
| 283 | |||
| 284 | 12725 | subroutine evaluate_meiss(f, r, th_c, ph_c, mode_secders) | |
| 285 | type(field_can_t), intent(inout) :: f | ||
| 286 | real(dp), intent(in) :: r, th_c, ph_c | ||
| 287 | integer, intent(in) :: mode_secders | ||
| 288 | |||
| 289 | real(dp) :: x(3) | ||
| 290 | |||
| 291 | 12725 | call count_field_evaluation(mode_secders) | |
| 292 | |||
| 293 |
2/2✓ Branch 0 taken 38175 times.
✓ Branch 1 taken 12725 times.
|
50900 | x = [r, th_c, ph_c] |
| 294 | |||
| 295 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 12720 times.
|
12725 | if (mode_secders > 0) then |
| 296 | 5 | call evaluate_meiss_batch_der2(f, x) | |
| 297 | 5 | return | |
| 298 | end if | ||
| 299 | |||
| 300 | 12720 | call evaluate_meiss_batch_der(f, x) | |
| 301 | end subroutine evaluate_meiss | ||
| 302 | |||
| 303 | |||
| 304 | 4 | subroutine integ_to_ref_meiss(xinteg, xref) | |
| 305 | real(dp), intent(in) :: xinteg(3) | ||
| 306 | real(dp), intent(out) :: xref(3) | ||
| 307 | real(dp) :: y_batch(2) ! lam_phi, chi_gauge | ||
| 308 | |||
| 309 | 4 | call evaluate_batch_splines_3d(spl_transform_batch, xinteg, y_batch) | |
| 310 | 4 | call coord_scaling%inverse(xinteg, xref) ! r -> s (integrator -> reference) | |
| 311 |
2/4✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
|
4 | xref(2) = modulo(xinteg(2), twopi) |
| 312 |
2/4✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
|
4 | xref(3) = modulo(xinteg(3) + y_batch(1), twopi) ! y_batch(1) is lam_phi |
| 313 | 4 | end subroutine integ_to_ref_meiss | |
| 314 | |||
| 315 | |||
| 316 | !> Reference radius of integ_to_ref_meiss without the angular splines. | ||
| 317 | 9 | function integ_to_ref_radius_meiss(xinteg) result(r_ref) | |
| 318 | real(dp), intent(in) :: xinteg(3) | ||
| 319 | real(dp) :: r_ref | ||
| 320 | real(dp) :: xref(3) | ||
| 321 | |||
| 322 | 9 | call coord_scaling%inverse(xinteg, xref) | |
| 323 | 9 | r_ref = xref(1) | |
| 324 | 9 | end function integ_to_ref_radius_meiss | |
| 325 | |||
| 326 | |||
| 327 | 18 | subroutine ref_to_integ_meiss(xref, xinteg) | |
| 328 | real(dp), intent(in) :: xref(3) | ||
| 329 | real(dp), intent(out) :: xinteg(3) | ||
| 330 | |||
| 331 | real(dp), parameter :: TOL = 1d-12 | ||
| 332 | integer, parameter :: MAX_ITER = 16 | ||
| 333 | |||
| 334 | real(dp) :: y_batch(2), dy_batch(3, 2), phi_integ_prev | ||
| 335 | integer :: i | ||
| 336 | |||
| 337 | 18 | call coord_scaling%transform(xref, xinteg) ! s -> r (reference -> integrator) | |
| 338 |
2/4✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
|
18 | xinteg(2) = modulo(xref(2), twopi) |
| 339 |
2/4✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
|
18 | xinteg(3) = modulo(xref(3), twopi) |
| 340 | |||
| 341 |
1/2✓ Branch 0 taken 63 times.
✗ Branch 1 not taken.
|
63 | do i=1, MAX_ITER |
| 342 | 63 | call evaluate_batch_splines_3d_der(spl_transform_batch, xinteg, y_batch, dy_batch) | |
| 343 | 63 | phi_integ_prev = xinteg(3) | |
| 344 | ! y_batch(1) is lam_phi, dy_batch(3,1) is d(lam_phi)/d(phi) | ||
| 345 | 63 | xinteg(3) = phi_integ_prev - (phi_integ_prev + y_batch(1) - xref(3))/(1d0 + dy_batch(3,1)) | |
| 346 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 45 times.
|
63 | if (abs(xinteg(3) - phi_integ_prev) < TOL) return |
| 347 | enddo | ||
| 348 | ✗ | print *, 'WARNING: ref_to_integ_meiss did not converge after', MAX_ITER, 'iterations' | |
| 349 | end subroutine ref_to_integ_meiss | ||
| 350 | |||
| 351 | |||
| 352 | 42 | subroutine get_meiss_coordinates(ierr) | |
| 353 | integer, intent(out), optional :: ierr | ||
| 354 | integer :: status | ||
| 355 | |||
| 356 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 357 | print *, 'field_can_meiss.init_transformation' | ||
| 358 | #endif | ||
| 359 | 42 | call init_transformation(status) | |
| 360 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 41 times.
|
42 | if (status /= 0) then |
| 361 | 1 | call cleanup_meiss | |
| 362 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (present(ierr)) then |
| 363 | 1 | ierr = status | |
| 364 | 1 | return | |
| 365 | end if | ||
| 366 | ✗ | error stop 'Meiss transformation ODE failed' | |
| 367 | end if | ||
| 368 | |||
| 369 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 370 | print *, 'field_can_meiss.spline_transformation' | ||
| 371 | #endif | ||
| 372 | 41 | call spline_transformation | |
| 373 | |||
| 374 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 375 | print *, 'field_can_meiss.init_canonical_field_components' | ||
| 376 | #endif | ||
| 377 | 41 | call init_canonical_field_components | |
| 378 |
2/2✓ Branch 0 taken 32 times.
✓ Branch 1 taken 9 times.
|
41 | if (present(ierr)) ierr = 0 |
| 379 | end subroutine get_meiss_coordinates | ||
| 380 | |||
| 381 | |||
| 382 | 90 | subroutine init_transformation(ierr) | |
| 383 | integer, intent(out), optional :: ierr | ||
| 384 | integer :: status | ||
| 385 | |||
| 386 | 45 | call init_transformation_arrays() | |
| 387 | 45 | call compute_transformation(status) | |
| 388 |
2/2✓ Branch 0 taken 42 times.
✓ Branch 1 taken 3 times.
|
45 | if (present(ierr)) then |
| 389 | 42 | ierr = status | |
| 390 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | else if (status /= 0) then |
| 391 | ✗ | error stop 'Meiss transformation ODE failed' | |
| 392 | end if | ||
| 393 | 45 | end subroutine init_transformation | |
| 394 | |||
| 395 | 45 | subroutine init_transformation_arrays() | |
| 396 | !> Initialize transformation arrays - can be called separately for diagnostics | ||
| 397 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 45 times.
|
45 | if (allocated(lam_phi)) deallocate(lam_phi) |
| 398 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 45 times.
|
45 | if (allocated(chi_gauge)) deallocate(chi_gauge) |
| 399 |
20/40✗ Branch 0 not taken.
✓ Branch 1 taken 45 times.
✓ Branch 2 taken 45 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 45 times.
✓ Branch 6 taken 45 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 45 times.
✓ Branch 10 taken 45 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 45 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 45 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 45 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 45 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 45 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 45 times.
✓ Branch 27 taken 45 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 30 taken 45 times.
✓ Branch 31 taken 45 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✓ Branch 34 taken 45 times.
✓ Branch 35 taken 45 times.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✓ Branch 38 taken 45 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 45 times.
✗ Branch 43 not taken.
✓ Branch 44 taken 45 times.
|
360 | allocate(lam_phi(n_r, n_th, n_phi), chi_gauge(n_r, n_th, n_phi)) |
| 400 | |||
| 401 | ! Initialize to zero (boundary conditions) | ||
| 402 |
6/6✓ Branch 0 taken 45 times.
✓ Branch 1 taken 942 times.
✓ Branch 2 taken 39028 times.
✓ Branch 3 taken 942 times.
✓ Branch 4 taken 1884856 times.
✓ Branch 5 taken 39028 times.
|
1924871 | lam_phi = 0.0_dp |
| 403 |
6/6✓ Branch 0 taken 942 times.
✓ Branch 1 taken 45 times.
✓ Branch 2 taken 39028 times.
✓ Branch 3 taken 942 times.
✓ Branch 4 taken 1884856 times.
✓ Branch 5 taken 39028 times.
|
1924871 | chi_gauge = 0.0_dp |
| 404 | 45 | end subroutine init_transformation_arrays | |
| 405 | |||
| 406 | 45 | subroutine compute_transformation(ierr) | |
| 407 | !> Compute transformation data via integration (expensive operation) | ||
| 408 | integer, intent(out) :: ierr | ||
| 409 | integer :: i_ctr | ||
| 410 | 45 | integer, allocatable :: slice_status(:) | |
| 411 | logical :: parallel_slices | ||
| 412 | |||
| 413 | i_ctr = 0 | ||
| 414 |
8/14✗ Branch 0 not taken.
✓ Branch 1 taken 45 times.
✓ Branch 2 taken 45 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 45 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 45 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 45 times.
✓ Branch 14 taken 942 times.
✓ Branch 15 taken 45 times.
|
1077 | allocate(slice_status(n_phi), source=0) |
| 415 | 45 | parallel_slices = .true. | |
| 416 | select type (field_noncan) | ||
| 417 | type is (spectre_field_t) | ||
| 418 | parallel_slices = .false. | ||
| 419 | end select | ||
| 420 | |||
| 421 | if (parallel_slices) then | ||
| 422 | 12 | !$omp parallel do private(i_ctr) | |
| 423 | do i_ctr = 1, n_phi | ||
| 424 | call compute_phi_slice(i_ctr, slice_status(i_ctr)) | ||
| 425 | end do | ||
| 426 | !$omp end parallel do | ||
| 427 | else | ||
| 428 |
2/2✓ Branch 0 taken 672 times.
✓ Branch 1 taken 33 times.
|
705 | do i_ctr = 1, n_phi |
| 429 | 705 | call compute_phi_slice(i_ctr, slice_status(i_ctr)) | |
| 430 | end do | ||
| 431 | end if | ||
| 432 |
2/2✓ Branch 0 taken 942 times.
✓ Branch 1 taken 45 times.
|
987 | ierr = maxval(slice_status) |
| 433 |
1/2✓ Branch 0 taken 45 times.
✗ Branch 1 not taken.
|
45 | end subroutine compute_transformation |
| 434 | |||
| 435 | 942 | subroutine compute_phi_slice(i_phi, ierr) | |
| 436 | !> Compute transformation for single phi slice | ||
| 437 | integer, intent(in) :: i_phi | ||
| 438 | integer, intent(out) :: ierr | ||
| 439 | real(dp) :: y(2) | ||
| 440 | integer :: i_th | ||
| 441 | |||
| 442 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 443 | call print_progress(i_phi) | ||
| 444 | #endif | ||
| 445 | |||
| 446 | 942 | ierr = 0 | |
| 447 |
2/2✓ Branch 0 taken 38974 times.
✓ Branch 1 taken 934 times.
|
39908 | do i_th = 1, n_th |
| 448 | 38974 | lam_phi(1, i_th, i_phi) = 0d0 | |
| 449 | 38974 | chi_gauge(1, i_th, i_phi) = 0d0 | |
| 450 | 38974 | y = 0d0 | |
| 451 | |||
| 452 |
2/2✓ Branch 1 taken 832 times.
✓ Branch 2 taken 38142 times.
|
39908 | if (is_hr_zero_on_slice(i_th, i_phi)) then |
| 453 | 832 | call set_identity_slice(i_th, i_phi) | |
| 454 | else | ||
| 455 | 38142 | call integrate_radial_slice(i_th, i_phi, y, ierr) | |
| 456 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 38134 times.
|
38142 | if (ierr /= 0) return |
| 457 | endif | ||
| 458 | enddo | ||
| 459 | end subroutine compute_phi_slice | ||
| 460 | |||
| 461 | 832 | subroutine set_identity_slice(i_th, i_phi) | |
| 462 | !> Set identity transformation for radial slice where hr ≈ 0 | ||
| 463 | integer, intent(in) :: i_th, i_phi | ||
| 464 | integer :: i_r | ||
| 465 | |||
| 466 |
2/2✓ Branch 0 taken 34604 times.
✓ Branch 1 taken 832 times.
|
35436 | do i_r = 2, n_r |
| 467 | 34604 | lam_phi(i_r, i_th, i_phi) = 0d0 | |
| 468 | 35436 | chi_gauge(i_r, i_th, i_phi) = 0d0 | |
| 469 | enddo | ||
| 470 | 832 | end subroutine set_identity_slice | |
| 471 | |||
| 472 | 38142 | subroutine integrate_radial_slice(i_th, i_phi, y, ierr) | |
| 473 | !> Integrate along radial direction for given (i_th, i_phi) | ||
| 474 | integer, intent(in) :: i_th, i_phi | ||
| 475 | real(dp), intent(inout) :: y(2) | ||
| 476 | integer, intent(out) :: ierr | ||
| 477 | integer :: i_r | ||
| 478 | |||
| 479 | 38142 | ierr = 0 | |
| 480 |
2/2✓ Branch 0 taken 1810798 times.
✓ Branch 1 taken 38134 times.
|
1848932 | do i_r = 2, n_r |
| 481 | 1810798 | call integrate(i_r, i_th, i_phi, y, ierr) | |
| 482 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1810790 times.
|
1848932 | if (ierr /= 0) then |
| 483 | 8 | write (*, '(A,I0,A,I0,A,I0)') 'Meiss ODE failed at r index ', i_r, & | |
| 484 | 16 | ', theta index ', i_th, ', phi index ', i_phi | |
| 485 | 8 | return | |
| 486 | end if | ||
| 487 | enddo | ||
| 488 | end subroutine integrate_radial_slice | ||
| 489 | |||
| 490 | ✗ | subroutine print_progress(i_phi) | |
| 491 | !> Print integration progress | ||
| 492 | integer, intent(in) :: i_phi | ||
| 493 | |||
| 494 | #ifndef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 495 | return | ||
| 496 | #endif | ||
| 497 | |||
| 498 | !$omp critical | ||
| 499 | write(*,'(A, I4, A, I4)',advance='no') 'integrate ODE: ', i_phi, ' of ', n_phi | ||
| 500 | if (i_phi < n_phi) then | ||
| 501 | write(*, '(A)', advance="no") char(13) | ||
| 502 | else | ||
| 503 | write(*, *) | ||
| 504 | end if | ||
| 505 | !$omp end critical | ||
| 506 | end subroutine print_progress | ||
| 507 | |||
| 508 | 1810798 | subroutine integrate(i_r, i_th, i_phi, y, ierr) | |
| 509 | use odeint_allroutines_sub, only: odeint_allroutines | ||
| 510 | |||
| 511 | integer, intent(in) :: i_r, i_th, i_phi | ||
| 512 | real(dp), dimension(2), intent(inout) :: y | ||
| 513 | integer, intent(out) :: ierr | ||
| 514 | |||
| 515 | integer, parameter :: ndim = 2 | ||
| 516 | real(dp) :: r1, r2 | ||
| 517 | real(dp) :: relaxed_relerr | ||
| 518 | type(grid_indices_t) :: context | ||
| 519 | #ifdef SIMPLE_NVHPC | ||
| 520 | real(dp), parameter :: hp_threshold = 1d-12 | ||
| 521 | real(dp) :: phi_c, Ar_test, Ap_test, hr_test, hp_test | ||
| 522 | #endif | ||
| 523 | |||
| 524 | 1810798 | r1 = xmin(1) + h_r*(i_r-2) | |
| 525 | 1810798 | r2 = xmin(1) + h_r*(i_r-1) | |
| 526 | |||
| 527 | 1810798 | relaxed_relerr = transformation_relerr | |
| 528 | |||
| 529 | #ifdef SIMPLE_NVHPC | ||
| 530 | ! NVHPC/nvfortran: avoid ODE step-size underflow near singular points by | ||
| 531 | ! relaxing tolerance when hp is very small at the start of the step. | ||
| 532 | phi_c = xmin(3) + h_phi*(i_phi-1) | ||
| 533 | call ah_cov_on_slice(r1, modulo(phi_c + y(1), twopi), i_th, Ar_test, Ap_test, hr_test, hp_test) | ||
| 534 | if (abs(hp_test) < hp_threshold) then | ||
| 535 | relaxed_relerr = max(transformation_relerr, 1d-8) | ||
| 536 | end if | ||
| 537 | #endif | ||
| 538 | |||
| 539 | 1810798 | context = grid_indices_t(i_th, i_phi) | |
| 540 | call odeint_allroutines(y, ndim, context, r1, r2, relaxed_relerr, & | ||
| 541 | 1810798 | rh_can_wrapper, step_limit=ode_step_limit, ierr=ierr) | |
| 542 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1810790 times.
|
1810798 | if (ierr /= 0) return |
| 543 | |||
| 544 | 1810790 | lam_phi(i_r, i_th, i_phi) = y(1) | |
| 545 | 1810790 | chi_gauge(i_r, i_th, i_phi) = y(2) | |
| 546 | |||
| 547 | end subroutine integrate | ||
| 548 | |||
| 549 | 32697907 | subroutine rh_can(r_c, z, dz, i_th, i_phi) | |
| 550 | real(dp), intent(in) :: r_c | ||
| 551 | real(dp), dimension(2), intent(in) :: z ! lam_phi, chi_gauge | ||
| 552 | real(dp), dimension(2), intent(out) :: dz | ||
| 553 | integer, intent(in) :: i_th, i_phi | ||
| 554 | real(dp) :: phi_c | ||
| 555 | real(dp) :: Ar, Ap, hr, hp | ||
| 556 | real(dp), parameter :: hp_min = 1d-14 | ||
| 557 | real(dp) :: hp_safe | ||
| 558 | |||
| 559 | 32697907 | phi_c = xmin(3) + h_phi*(i_phi-1) | |
| 560 |
4/4✓ Branch 0 taken 32695782 times.
✓ Branch 1 taken 2125 times.
✓ Branch 2 taken 388690 times.
✓ Branch 3 taken 32307092 times.
|
32697907 | call ah_cov_on_slice(r_c, modulo(phi_c + z(1), twopi), i_th, Ar, Ap, hr, hp) |
| 561 | |||
| 562 | 32697907 | hp_safe = hp | |
| 563 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 32697907 times.
|
32697907 | if (abs(hp_safe) < hp_min) then |
| 564 | ✗ | hp_safe = sign(hp_min, hp_safe + 0.5d0*hp_min) | |
| 565 | end if | ||
| 566 | |||
| 567 | 32697907 | dz(1) = -hr / hp_safe | |
| 568 | 32697907 | dz(2) = Ar + Ap * dz(1) | |
| 569 | 32697907 | end subroutine rh_can | |
| 570 | |||
| 571 | |||
| 572 | 44 | subroutine spline_transformation | |
| 573 | real(dp), dimension(:,:,:,:), allocatable :: y_batch | ||
| 574 | integer :: dims(3) | ||
| 575 | |||
| 576 | ! Construct batch spline for transformation components (2: lam_phi, chi_gauge) | ||
| 577 |
2/2✓ Branch 0 taken 132 times.
✓ Branch 1 taken 44 times.
|
176 | dims = shape(lam_phi) |
| 578 |
11/22✗ Branch 0 not taken.
✓ Branch 1 taken 44 times.
✓ Branch 2 taken 44 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 44 times.
✓ Branch 6 taken 44 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 44 times.
✓ Branch 10 taken 44 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 44 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 44 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 44 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 44 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 44 times.
|
264 | allocate(y_batch(dims(1), dims(2), dims(3), 2)) |
| 579 | |||
| 580 |
6/6✓ Branch 0 taken 934 times.
✓ Branch 1 taken 44 times.
✓ Branch 2 taken 38964 times.
✓ Branch 3 taken 934 times.
✓ Branch 4 taken 1884344 times.
✓ Branch 5 taken 38964 times.
|
1924286 | y_batch(:,:,:,1) = lam_phi |
| 581 |
6/6✓ Branch 0 taken 934 times.
✓ Branch 1 taken 44 times.
✓ Branch 2 taken 38964 times.
✓ Branch 3 taken 934 times.
✓ Branch 4 taken 1884344 times.
✓ Branch 5 taken 38964 times.
|
1924286 | y_batch(:,:,:,2) = chi_gauge |
| 582 | |||
| 583 | 44 | call construct_batch_splines_3d(xmin, xmax, y_batch, order, periodic, spl_transform_batch) | |
| 584 | 44 | transform_splines_initialized = .true. | |
| 585 |
1/2✓ Branch 0 taken 44 times.
✗ Branch 1 not taken.
|
44 | end subroutine spline_transformation |
| 586 | |||
| 587 | |||
| 588 | 32763393 | subroutine ah_cov_on_slice(r, phi, i_th, Ar, Ap, hr, hp) | |
| 589 | real(dp), intent(in) :: r, phi | ||
| 590 | integer, intent(in) :: i_th | ||
| 591 | real(dp), intent(out) :: Ar, Ap, hr, hp | ||
| 592 | |||
| 593 | real(dp), dimension(3) :: Acov, hcov, x_integ, x_ref, jac | ||
| 594 | real(dp) :: Bmod | ||
| 595 | real(dp) :: th | ||
| 596 | |||
| 597 | 32763393 | th = xmin(2) + h_th*(i_th-1) | |
| 598 |
2/2✓ Branch 0 taken 98290179 times.
✓ Branch 1 taken 32763393 times.
|
131053572 | x_integ = [r, th, phi] |
| 599 | |||
| 600 | select type (fld => field_noncan) | ||
| 601 | type is (vmec_field_t) | ||
| 602 | ! VMEC field expects s-coordinates, so convert r->s and apply Jacobian | ||
| 603 | 22576153 | call coord_scaling%inverse(x_integ, x_ref, jac) | |
| 604 | 22576153 | call fld%evaluate(x_ref, Acov, hcov, Bmod) | |
| 605 | ! Convert output from s-coords to r-coords: A_r = A_s * ds/dr | ||
| 606 | 22576153 | Ar = Acov(1) * jac(1) | |
| 607 | 22576153 | Ap = Acov(3) | |
| 608 | 22576153 | hr = hcov(1) * jac(1) | |
| 609 | 45152306 | hp = hcov(3) | |
| 610 | type is (splined_field_t) | ||
| 611 | ! Splined field already stores data in scaled coordinates (r = sqrt(s)) | ||
| 612 | ! matching the Meiss integrator coordinates - no conversion needed | ||
| 613 | 1067988 | call fld%evaluate(x_integ, Acov, hcov, Bmod) | |
| 614 | 1067988 | Ar = Acov(1) | |
| 615 | 1067988 | Ap = Acov(3) | |
| 616 | 1067988 | hr = hcov(1) | |
| 617 | 1067988 | hp = hcov(3) | |
| 618 | type is (spectre_field_t) | ||
| 619 | ! SPECTRE chart (rho_g, theta, zeta) with identity scaling; covariant | ||
| 620 | ! components are used raw in SI, the Gaussian-CGS conversion is applied at | ||
| 621 | ! canonical evaluation. Here h_s = hcov(1) is nonzero (the reason the Meiss | ||
| 622 | ! angle transform is required), while A_s = Acov(1) = 0 in the SPECTRE gauge. | ||
| 623 | 9119252 | call bias_spectre_upper_endpoint(fld, x_integ) | |
| 624 | 9119252 | call fld%evaluate(x_integ, Acov, hcov, Bmod) | |
| 625 | 9119252 | Ar = Acov(1) | |
| 626 | 9119252 | Ap = Acov(3) | |
| 627 | 9119252 | hr = hcov(1) | |
| 628 | 18238504 | hp = hcov(3) | |
| 629 | class default | ||
| 630 | ✗ | error stop "ah_cov_on_slice: unsupported field type" | |
| 631 | end select | ||
| 632 | 32763393 | end subroutine ah_cov_on_slice | |
| 633 | |||
| 634 | 38974 | logical function is_hr_zero_on_slice(i_th, i_phi) result(is_zero) | |
| 635 | !> Check if hr ≈ 0 along the entire radial direction for given (i_th, i_phi) | ||
| 636 | !> Sample hr at several radial points to make this determination | ||
| 637 | integer, intent(in) :: i_th, i_phi | ||
| 638 | |||
| 639 | real(dp), parameter :: hr_slice_threshold = 1d-9 | ||
| 640 | integer, parameter :: n_sample_points = 32 | ||
| 641 | real(dp) :: r_sample, phi_c, Ar, Ap, hr, hp | ||
| 642 | integer :: i_sample | ||
| 643 | logical :: all_zero | ||
| 644 | |||
| 645 | 38974 | phi_c = xmin(3) + h_phi*(i_phi-1) | |
| 646 | 38974 | all_zero = .true. | |
| 647 | |||
| 648 | ! Sample hr at several radial points | ||
| 649 |
2/2✓ Branch 0 taken 64766 times.
✓ Branch 1 taken 832 times.
|
65598 | do i_sample = 1, n_sample_points |
| 650 | ! Distribute sample points across radial range | ||
| 651 | 64766 | r_sample = xmin(1) + (xmax(1) - xmin(1)) * real(i_sample-1, dp) / real(n_sample_points-1, dp) | |
| 652 | |||
| 653 | ! Evaluate field at this point (using lam_phi=0 for sampling) | ||
| 654 | 64766 | call ah_cov_on_slice(r_sample, phi_c, i_th, Ar, Ap, hr, hp) | |
| 655 | |||
| 656 |
2/2✓ Branch 1 taken 26624 times.
✓ Branch 2 taken 38142 times.
|
65598 | if (abs(hr) >= hr_slice_threshold) then |
| 657 | all_zero = .false. | ||
| 658 | exit ! Found non-zero hr, no need to check further | ||
| 659 | endif | ||
| 660 | enddo | ||
| 661 | |||
| 662 | 38974 | is_zero = all_zero | |
| 663 | 38974 | end function is_hr_zero_on_slice | |
| 664 | |||
| 665 | |||
| 666 | 44 | subroutine init_canonical_field_components | |
| 667 | real(dp) :: xcan(3) | ||
| 668 | 44 | real(dp), dimension(:,:,:,:), allocatable :: y_batch | |
| 669 | real(dp) :: xref(3), Acov(3), hcov(3), Bmod, lam, dlam(3), chi, dchi(3) | ||
| 670 | integer :: i_r, i_th, i_phi | ||
| 671 | |||
| 672 |
1/2✓ Branch 0 taken 44 times.
✗ Branch 1 not taken.
|
44 | if (allocated(lam_phi)) deallocate(lam_phi) |
| 673 |
1/2✓ Branch 0 taken 44 times.
✗ Branch 1 not taken.
|
44 | if (allocated(chi_gauge)) deallocate(chi_gauge) |
| 674 |
11/22✗ Branch 0 not taken.
✓ Branch 1 taken 44 times.
✓ Branch 2 taken 44 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 44 times.
✓ Branch 6 taken 44 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 44 times.
✓ Branch 10 taken 44 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 44 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 44 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 44 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 44 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 44 times.
|
264 | allocate(y_batch(n_r, n_th, n_phi, 5)) |
| 675 | |||
| 676 |
2/2✓ Branch 0 taken 934 times.
✓ Branch 1 taken 44 times.
|
978 | do i_phi=1,n_phi |
| 677 |
2/2✓ Branch 0 taken 38964 times.
✓ Branch 1 taken 934 times.
|
39942 | do i_th=1,n_th |
| 678 |
2/2✓ Branch 0 taken 1884344 times.
✓ Branch 1 taken 38964 times.
|
1924242 | do i_r=1,n_r |
| 679 |
2/2✓ Branch 1 taken 5653032 times.
✓ Branch 2 taken 1884344 times.
|
7537376 | xcan = get_grid_point(i_r, i_th, i_phi) |
| 680 | |||
| 681 | ! Use batch evaluation for both transformation components | ||
| 682 | block | ||
| 683 | real(dp) :: y_trans(2), dy_trans(3,2) | ||
| 684 | 1884344 | call evaluate_batch_splines_3d_der(spl_transform_batch, xcan, y_trans, dy_trans) | |
| 685 | 1884344 | lam = y_trans(1) ! lam_phi | |
| 686 | 1884344 | chi = y_trans(2) ! chi_gauge | |
| 687 |
2/2✓ Branch 0 taken 5653032 times.
✓ Branch 1 taken 1884344 times.
|
7537376 | dlam = dy_trans(:,1) ! derivatives of lam_phi |
| 688 |
2/2✓ Branch 0 taken 5653032 times.
✓ Branch 1 taken 1884344 times.
|
7537376 | dchi = dy_trans(:,2) ! derivatives of chi_gauge |
| 689 | end block | ||
| 690 | |||
| 691 | ! xcan is in integrator coords (r, th, phi) | ||
| 692 | ! xref here is the non-canonical angle coords (th, phi + lam) | ||
| 693 |
9/10✓ Branch 0 taken 1836156 times.
✓ Branch 1 taken 48188 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1836156 times.
✓ Branch 4 taken 1877733 times.
✓ Branch 5 taken 6611 times.
✓ Branch 6 taken 29471 times.
✓ Branch 7 taken 1848262 times.
✓ Branch 8 taken 5653032 times.
✓ Branch 9 taken 1884344 times.
|
7537376 | xref = [xcan(1), modulo(xcan(2), twopi), modulo(xcan(3) + lam, twopi)] |
| 694 | |||
| 695 | select type (fld => field_noncan) | ||
| 696 | type is (vmec_field_t) | ||
| 697 | ! VMEC field expects s-coordinates | ||
| 698 | 514392 | block | |
| 699 | real(dp) :: x_integ_local(3), x_ref_local(3), jac_local(3) | ||
| 700 | 514392 | x_integ_local = xref | |
| 701 | 514392 | call coord_scaling%inverse(x_integ_local, x_ref_local, jac_local) | |
| 702 | 514392 | call fld%evaluate(x_ref_local, Acov, hcov, Bmod) | |
| 703 | end block | ||
| 704 | type is (splined_field_t) | ||
| 705 | ! Splined field already in scaled coordinates | ||
| 706 | 14688 | call fld%evaluate(xref, Acov, hcov, Bmod) | |
| 707 | type is (spectre_field_t) | ||
| 708 | ! SPECTRE chart is its own radial coordinate (identity scaling) | ||
| 709 | 1355264 | call bias_spectre_upper_endpoint(fld, xref) | |
| 710 | 1355264 | call fld%evaluate(xref, Acov, hcov, Bmod) | |
| 711 | class default | ||
| 712 | ✗ | error stop "init_canonical_field_components: unsupported field type" | |
| 713 | end select | ||
| 714 | |||
| 715 | ! Note: We only use theta/phi components of Acov/hcov below | ||
| 716 | |||
| 717 | y_batch(i_r, i_th, i_phi, 1) = & | ||
| 718 | 1884344 | Acov(2) + Acov(3)*dlam(2) - dchi(2) | |
| 719 | y_batch(i_r, i_th, i_phi, 2) = & | ||
| 720 | 1884344 | Acov(3)*(1.0d0 + dlam(3)) - dchi(3) | |
| 721 | 1884344 | y_batch(i_r, i_th, i_phi, 3) = hcov(2) + hcov(3)*dlam(2) | |
| 722 | 1884344 | y_batch(i_r, i_th, i_phi, 4) = hcov(3)*(1.0d0 + dlam(3)) | |
| 723 | 1923308 | y_batch(i_r, i_th, i_phi, 5) = Bmod | |
| 724 | end do | ||
| 725 | end do | ||
| 726 | end do | ||
| 727 | |||
| 728 | 44 | call construct_batch_splines_3d(xmin, xmax, y_batch, order, periodic, spl_field_batch) | |
| 729 | 44 | batch_splines_initialized = .true. | |
| 730 |
1/2✓ Branch 0 taken 44 times.
✗ Branch 1 not taken.
|
44 | end subroutine init_canonical_field_components |
| 731 | |||
| 732 | |||
| 733 | 10474516 | subroutine bias_spectre_upper_endpoint(field, x) | |
| 734 | type(spectre_field_t), intent(in) :: field | ||
| 735 | real(dp), intent(inout) :: x(3) | ||
| 736 | |||
| 737 | ! Exact stacked interfaces belong to the outer volume. A lower-volume | ||
| 738 | ! construction must sample its upper endpoint as the one-sided inner limit. | ||
| 739 |
4/4✓ Branch 0 taken 57730 times.
✓ Branch 1 taken 10416786 times.
✓ Branch 2 taken 20610 times.
✓ Branch 3 taken 37120 times.
|
10474516 | if (x(1) >= xmax(1) .and. xmax(1) < real(field%data%Mvol, dp)) then |
| 740 | 20610 | x(1) = nearest(xmax(1), -1.0_dp) | |
| 741 | end if | ||
| 742 | 10474516 | end subroutine bias_spectre_upper_endpoint | |
| 743 | |||
| 744 | |||
| 745 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1884344 times.
|
1884344 | pure function get_grid_point(i_r, i_th, i_phi) |
| 746 | integer, intent(in) :: i_r, i_th, i_phi | ||
| 747 | real(dp) :: get_grid_point(3) | ||
| 748 | |||
| 749 | get_grid_point = [ & | ||
| 750 | xmin(1) + h_r*dble(i_r-1), & | ||
| 751 | xmin(2) + h_th*dble(i_th-1), & | ||
| 752 | xmin(3) + h_phi*dble(i_phi-1) & | ||
| 753 |
2/2✓ Branch 0 taken 5653032 times.
✓ Branch 1 taken 1884344 times.
|
7537376 | ] |
| 754 | end function get_grid_point | ||
| 755 | |||
| 756 | |||
| 757 | 12000 | subroutine magfie_meiss(x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl) | |
| 758 | ! Computes magnetic field and derivatives with bmod in units of the magnetic code | ||
| 759 | ! | ||
| 760 | ! Input parameters: | ||
| 761 | ! formal: x - array of canonicalized coordinates r, th, ph | ||
| 762 | ! Output parameters: | ||
| 763 | ! formal: bmod - magnetic field module | ||
| 764 | ! sqrtg - metric determinant | ||
| 765 | ! bder - covariant components of (grad B)/B | ||
| 766 | ! hcovar - covariant components of \bB/B | ||
| 767 | ! hctrvr - contravariant components of \bB/B | ||
| 768 | ! hcurl - contravariant components of curl (\bB/B) | ||
| 769 | real(dp), intent(in) :: x(3) | ||
| 770 | real(dp), intent(out) :: bmod, sqrtg | ||
| 771 | real(dp), dimension(3), intent(out) :: bder, hcovar, hctrvr, hcurl | ||
| 772 | |||
| 773 | type(field_can_t) :: f | ||
| 774 | real(dp) :: sqrtg_bmod | ||
| 775 | |||
| 776 | 12000 | call evaluate_meiss(f, x(1), x(2), x(3), 0) | |
| 777 | |||
| 778 | 12000 | bmod = f%Bmod | |
| 779 | |||
| 780 | 12000 | sqrtg_bmod = f%hph*f%dAth(1) - f%hth*f%dAph(1) | |
| 781 | 12000 | sqrtg = sqrtg_bmod/bmod | |
| 782 |
2/2✓ Branch 0 taken 36000 times.
✓ Branch 1 taken 12000 times.
|
48000 | bder = f%dBmod/bmod |
| 783 | |||
| 784 | 12000 | hcovar(1) = 0d0 | |
| 785 | 12000 | hcovar(2) = f%hth | |
| 786 | 12000 | hcovar(3) = f%hph | |
| 787 | |||
| 788 | 12000 | hctrvr(1) = (f%dAph(2) - f%dAth(3))/sqrtg_bmod | |
| 789 | 12000 | hctrvr(2) = -f%dAph(1)/sqrtg_bmod | |
| 790 | 12000 | hctrvr(3) = f%dAth(1)/sqrtg_bmod | |
| 791 | |||
| 792 | 12000 | hcurl(1) = (f%dhph(2) - f%dhth(3))/sqrtg | |
| 793 | 12000 | hcurl(2) = -f%dhph(1)/sqrtg | |
| 794 | 12000 | hcurl(3) = f%dhth(1)/sqrtg | |
| 795 |
4/4✓ Branch 0 taken 12000 times.
✓ Branch 1 taken 36000 times.
✓ Branch 2 taken 72000 times.
✓ Branch 3 taken 12000 times.
|
132000 | end subroutine magfie_meiss |
| 796 | |||
| 797 | |||
| 798 | ! Batch evaluation helper for first derivatives | ||
| 799 | 12720 | subroutine evaluate_meiss_batch_der(f, x) | |
| 800 | type(field_can_t), intent(inout) :: f | ||
| 801 | real(dp), intent(in) :: x(3) | ||
| 802 | |||
| 803 | real(dp) :: y_batch(5), dy_batch(3, 5) | ||
| 804 | |||
| 805 | 12720 | call evaluate_batch_splines_3d_der(spl_field_batch, x, y_batch, dy_batch) | |
| 806 | |||
| 807 | ! Unpack results: order is [Ath, Aph, hth, hph, Bmod] | ||
| 808 | 12720 | f%Ath = y_batch(1) | |
| 809 | 12720 | f%Aph = y_batch(2) | |
| 810 | 12720 | f%hth = y_batch(3) | |
| 811 | 12720 | f%hph = y_batch(4) | |
| 812 | 12720 | f%Bmod = y_batch(5) | |
| 813 | |||
| 814 |
2/2✓ Branch 0 taken 12720 times.
✓ Branch 1 taken 38160 times.
|
50880 | f%dAth = dy_batch(:, 1) |
| 815 |
2/2✓ Branch 0 taken 12720 times.
✓ Branch 1 taken 38160 times.
|
50880 | f%dAph = dy_batch(:, 2) |
| 816 |
2/2✓ Branch 0 taken 12720 times.
✓ Branch 1 taken 38160 times.
|
50880 | f%dhth = dy_batch(:, 3) |
| 817 |
2/2✓ Branch 0 taken 12720 times.
✓ Branch 1 taken 38160 times.
|
50880 | f%dhph = dy_batch(:, 4) |
| 818 |
2/2✓ Branch 0 taken 38160 times.
✓ Branch 1 taken 12720 times.
|
50880 | f%dBmod = dy_batch(:, 5) |
| 819 | 12720 | end subroutine evaluate_meiss_batch_der | |
| 820 | |||
| 821 | |||
| 822 | ! Batch evaluation helper for second derivatives | ||
| 823 | 5 | subroutine evaluate_meiss_batch_der2(f, x) | |
| 824 | type(field_can_t), intent(inout) :: f | ||
| 825 | real(dp), intent(in) :: x(3) | ||
| 826 | |||
| 827 | real(dp) :: y_batch(5), dy_batch(3, 5), d2y_batch(6, 5) | ||
| 828 | |||
| 829 | 5 | call evaluate_batch_splines_3d_der2(spl_field_batch, x, y_batch, dy_batch, d2y_batch) | |
| 830 | |||
| 831 | ! Unpack results: order is [Ath, Aph, hth, hph, Bmod] | ||
| 832 | 5 | f%Ath = y_batch(1) | |
| 833 | 5 | f%Aph = y_batch(2) | |
| 834 | 5 | f%hth = y_batch(3) | |
| 835 | 5 | f%hph = y_batch(4) | |
| 836 | 5 | f%Bmod = y_batch(5) | |
| 837 | |||
| 838 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 15 times.
|
20 | f%dAth = dy_batch(:, 1) |
| 839 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 15 times.
|
20 | f%dAph = dy_batch(:, 2) |
| 840 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 15 times.
|
20 | f%dhth = dy_batch(:, 3) |
| 841 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 15 times.
|
20 | f%dhph = dy_batch(:, 4) |
| 842 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 15 times.
|
20 | f%dBmod = dy_batch(:, 5) |
| 843 | |||
| 844 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 30 times.
|
35 | f%d2Ath = d2y_batch(:, 1) |
| 845 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 30 times.
|
35 | f%d2Aph = d2y_batch(:, 2) |
| 846 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 30 times.
|
35 | f%d2hth = d2y_batch(:, 3) |
| 847 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 30 times.
|
35 | f%d2hph = d2y_batch(:, 4) |
| 848 |
2/2✓ Branch 0 taken 30 times.
✓ Branch 1 taken 5 times.
|
35 | f%d2Bmod = d2y_batch(:, 5) |
| 849 | 5 | end subroutine evaluate_meiss_batch_der2 | |
| 850 | |||
| 851 | ✗ | end module field_can_meiss | |
| 852 |