| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module field_can_spectre | ||
| 2 | !> Per-volume Meiss canonical coordinates for SPECTRE equilibria (#439). | ||
| 3 | !> | ||
| 4 | !> SPECTRE's A_s = 0 gauge leaves the covariant h_s /= 0 (metric off-diagonals), | ||
| 5 | !> so the phase-space 1-form keeps a rho_par*B_s ds term. Each SPECTRE volume is | ||
| 6 | !> smooth in its own stacked-rho slab r in [lvol-1, lvol]; a global chart would | ||
| 7 | !> integrate the Meiss construction ODE through the interface current sheets. | ||
| 8 | !> This module runs the existing Meiss machinery once per volume on the identity | ||
| 9 | !> radial scaling (the SPECTRE chart is its own radial coordinate) and harvests | ||
| 10 | !> the resulting batch splines into a per-volume slot. | ||
| 11 | !> | ||
| 12 | !> Construction stays in the field's SI covariant units. The Gaussian-CGS | ||
| 13 | !> conversion SIMPLE integrates in is applied only here at evaluation, uniformly | ||
| 14 | !> to values and derivatives because r, theta, zeta are all dimensionless: | ||
| 15 | !> Bmod Tesla->Gauss, covariant h meter->cm, covariant A (flux-like, [B]*L^2) | ||
| 16 | !> Tesla*m^2 -> Gauss*cm^2. | ||
| 17 | |||
| 18 | use, intrinsic :: iso_fortran_env, only: dp => real64, int64 | ||
| 19 | use field_base, only: magnetic_field_t | ||
| 20 | use field_spectre, only: spectre_field_t | ||
| 21 | use field_can_base, only: count_field_evaluation, field_can_t, twopi | ||
| 22 | use finite_check, only: is_finite | ||
| 23 | use field_can_meiss, only: meiss_volume_t, init_meiss, get_meiss_coordinates, & | ||
| 24 | harvest_meiss_volume, cleanup_meiss | ||
| 25 | use coordinate_scaling, only: identity_scaling_t | ||
| 26 | use interpolate, only: evaluate_batch_splines_3d, & | ||
| 27 | evaluate_batch_splines_3d_der, & | ||
| 28 | evaluate_batch_splines_3d_der2 | ||
| 29 | use magfie_sub, only: TESLA_TO_GAUSS, M_TO_CM | ||
| 30 | use diag_counters, only: count_event, EVT_SPECTRE_REF_INVERSE_MAXIT, & | ||
| 31 | EVT_SPECTRE_INVALID_STATE | ||
| 32 | |||
| 33 | implicit none | ||
| 34 | private | ||
| 35 | |||
| 36 | public :: construct_spectre_coordinates, evaluate_spectre_can, & | ||
| 37 | integ_to_ref_spectre, ref_to_integ_spectre, cleanup_spectre, & | ||
| 38 | spectre_volumes, spectre_mvol, set_spectre_volume_lock, & | ||
| 39 | set_spectre_construction_grid | ||
| 40 | |||
| 41 | !> Per-volume Meiss construction grid. phi = -1 selects the automatic policy | ||
| 42 | !> below. set_spectre_construction_grid | ||
| 43 | !> threads the namelist knob (params: spectre_ncon_r/th/phi) in before the build, | ||
| 44 | !> kept as a setter because field_can_* cannot use params (dependency cycle). | ||
| 45 | integer :: ncon_r = 48, ncon_th = 48, ncon_phi = -1, ncon_order = 5 | ||
| 46 | integer :: ncon_ode_step_limit = 1000000 | ||
| 47 | real(dp) :: ncon_ode_relerr = 1.0e-2_dp | ||
| 48 | |||
| 49 | !> Automatic phi resolution. Axisymmetric SPECTRE | ||
| 50 | !> fields (all n = 0) are phi-invariant, so the construction only needs a minimal | ||
| 51 | !> periodic-quintic phi grid; AXISYM_NPHI is the dominant memory saver for tokamak | ||
| 52 | !> cases (arrays and splines scale linearly in n_phi). It stays above the quintic | ||
| 53 | !> order+1 = 6 so spl_five_per is well conditioned, and preserves the symplectic | ||
| 54 | !> energy invariants; a positive spectre_ncon_phi overrides it (raise to NPHI_FULL | ||
| 55 | !> for high-order convergence, lower for RK45-only or memory-bound runs). | ||
| 56 | integer, parameter :: AXISYM_NPHI = 8, NPHI_FULL = 32 | ||
| 57 | real(dp), parameter :: AXIS_OFFSET = 1.0e-3_dp | ||
| 58 | !> C1 radial extension used by canonical field evaluation at volume edges. | ||
| 59 | !> This boundary regularization is retained for single-volume tracing and | ||
| 60 | !> per-volume field diagnostics; it supplies no interface crossing model. | ||
| 61 | real(dp), parameter :: EDGE_BAND = 5.0e-2_dp | ||
| 62 | real(dp), parameter :: A_SCALE = TESLA_TO_GAUSS*M_TO_CM**2 | ||
| 63 | real(dp), parameter :: H_SCALE = M_TO_CM | ||
| 64 | real(dp), parameter :: B_SCALE = TESLA_TO_GAUSS | ||
| 65 | |||
| 66 | type(meiss_volume_t), allocatable :: spectre_volumes(:) | ||
| 67 | integer :: spectre_mvol = 0 | ||
| 68 | |||
| 69 | !> Pin Newton trials and coordinate transforms to the guiding centre owner. | ||
| 70 | integer :: lock_lvol = 0 | ||
| 71 | !$omp threadprivate(lock_lvol) | ||
| 72 | |||
| 73 | contains | ||
| 74 | |||
| 75 | |||
| 76 | 649 | subroutine set_spectre_volume_lock(lvol) | |
| 77 | !> lvol > 0 pins field evaluation to that volume; 0 restores dispatch on | ||
| 78 | !> int(r). | ||
| 79 | integer, intent(in) :: lvol | ||
| 80 | |||
| 81 | 649 | lock_lvol = lvol | |
| 82 | 649 | end subroutine set_spectre_volume_lock | |
| 83 | |||
| 84 | 6393862 | integer function active_volume(r) result(lvol) | |
| 85 | real(dp), intent(in) :: r | ||
| 86 | |||
| 87 | 6393862 | if (lock_lvol > 0) then | |
| 88 | 3990915 | lvol = min(lock_lvol, spectre_mvol) | |
| 89 | else | ||
| 90 | 2402947 | lvol = volume_index(r) | |
| 91 | end if | ||
| 92 | end function active_volume | ||
| 93 | |||
| 94 | 19 | subroutine set_spectre_construction_grid(n_r, n_th, n_phi, spline_order, & | |
| 95 | ode_step_limit, ode_relerr) | ||
| 96 | !> Set the per-volume Meiss construction grid before the build (n_phi <= 0 | ||
| 97 | !> selects the automatic phi policy). Called from init_spectre_field. | ||
| 98 | integer, intent(in) :: n_r, n_th, n_phi | ||
| 99 | integer, intent(in), optional :: spline_order | ||
| 100 | integer, intent(in), optional :: ode_step_limit | ||
| 101 | real(dp), intent(in), optional :: ode_relerr | ||
| 102 | |||
| 103 | 19 | ncon_r = n_r | |
| 104 | 19 | ncon_th = n_th | |
| 105 | 19 | ncon_phi = n_phi | |
| 106 | 19 | ncon_order = 5 | |
| 107 |
1/2✓ Branch 0 taken 19 times.
✗ Branch 1 not taken.
|
19 | if (present(spline_order)) ncon_order = spline_order |
| 108 | 19 | ncon_ode_step_limit = 1000000 | |
| 109 |
1/2✓ Branch 0 taken 19 times.
✗ Branch 1 not taken.
|
19 | if (present(ode_step_limit)) ncon_ode_step_limit = ode_step_limit |
| 110 | 19 | ncon_ode_relerr = 1.0e-2_dp | |
| 111 |
1/2✓ Branch 0 taken 19 times.
✗ Branch 1 not taken.
|
19 | if (present(ode_relerr)) ncon_ode_relerr = ode_relerr |
| 112 | 19 | end subroutine set_spectre_construction_grid | |
| 113 | |||
| 114 | 20 | subroutine construct_spectre_coordinates(field_noncan, ierr) | |
| 115 | !> Build one Meiss chart per SPECTRE volume on r in [lvol-1, lvol]. The axis | ||
| 116 | !> volume starts at AXIS_OFFSET so the construction never touches the r = 0 | ||
| 117 | !> coordinate singularity. | ||
| 118 | class(magnetic_field_t), intent(in) :: field_noncan | ||
| 119 | integer, intent(out), optional :: ierr | ||
| 120 | |||
| 121 | integer :: lvol, nphi_eff, meiss_ierr | ||
| 122 | real(dp) :: rmin, rmax | ||
| 123 | logical :: axisym | ||
| 124 | |||
| 125 | 20 | call cleanup_spectre | |
| 126 | |||
| 127 | select type (fld => field_noncan) | ||
| 128 | type is (spectre_field_t) | ||
| 129 | 20 | spectre_mvol = fld%data%Mvol | |
| 130 |
4/4✓ Branch 0 taken 79 times.
✓ Branch 1 taken 13 times.
✓ Branch 2 taken 72 times.
✓ Branch 3 taken 7 times.
|
92 | axisym = all(fld%data%in == 0) |
| 131 | class default | ||
| 132 | ✗ | error stop 'construct_spectre_coordinates: field is not spectre_field_t' | |
| 133 | end select | ||
| 134 | |||
| 135 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
20 | if (ncon_phi > 0) then |
| 136 | 10 | nphi_eff = ncon_phi | |
| 137 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 4 times.
|
10 | else if (axisym) then |
| 138 | 6 | nphi_eff = AXISYM_NPHI | |
| 139 | else | ||
| 140 | 4 | nphi_eff = NPHI_FULL | |
| 141 | end if | ||
| 142 | |||
| 143 | print '(A,I0,A,I0,A,L1,A,I0,A,L1,A,I0,A,I0,A,ES9.2)', & | ||
| 144 | 20 | ' spectre_construction_grid: n_r=', ncon_r, ' n_th=', ncon_th, & | |
| 145 | 20 | ' axisym=', axisym, ' n_phi=', nphi_eff, & | |
| 146 | 20 | ' phi_auto=', ncon_phi <= 0, ' order=', ncon_order, & | |
| 147 | 40 | ' ode_steps=', ncon_ode_step_limit, ' ode_relerr=', ncon_ode_relerr | |
| 148 | |||
| 149 |
37/44✗ Branch 0 not taken.
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 20 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 20 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 20 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 20 times.
✓ Branch 17 taken 60 times.
✓ Branch 18 taken 20 times.
✓ Branch 19 taken 60 times.
✓ Branch 20 taken 20 times.
✓ Branch 21 taken 60 times.
✓ Branch 22 taken 20 times.
✓ Branch 23 taken 60 times.
✓ Branch 24 taken 20 times.
✓ Branch 25 taken 60 times.
✓ Branch 26 taken 20 times.
✓ Branch 27 taken 60 times.
✓ Branch 28 taken 20 times.
✓ Branch 29 taken 60 times.
✓ Branch 30 taken 20 times.
✓ Branch 31 taken 60 times.
✓ Branch 32 taken 20 times.
✓ Branch 33 taken 60 times.
✓ Branch 34 taken 20 times.
✓ Branch 35 taken 60 times.
✓ Branch 36 taken 20 times.
✓ Branch 37 taken 60 times.
✓ Branch 38 taken 20 times.
✓ Branch 39 taken 60 times.
✓ Branch 40 taken 20 times.
✓ Branch 41 taken 60 times.
✓ Branch 42 taken 20 times.
✓ Branch 43 taken 60 times.
✓ Branch 44 taken 20 times.
✓ Branch 45 taken 34 times.
✓ Branch 46 taken 20 times.
|
934 | allocate (spectre_volumes(spectre_mvol)) |
| 150 | |||
| 151 |
2/2✓ Branch 0 taken 33 times.
✓ Branch 1 taken 19 times.
|
52 | do lvol = 1, spectre_mvol |
| 152 | 33 | rmin = real(lvol - 1, dp) | |
| 153 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 13 times.
|
33 | if (lvol == 1) rmin = rmin + AXIS_OFFSET |
| 154 | 33 | rmax = real(lvol, dp) | |
| 155 | |||
| 156 | call init_meiss(field_noncan, ncon_r, ncon_th, nphi_eff, rmin, rmax, & | ||
| 157 | 0.0_dp, twopi, identity_scaling_t(), & | ||
| 158 | transformation_relerr_=ncon_ode_relerr, & | ||
| 159 | spline_order_=ncon_order, & | ||
| 160 | 33 | ode_step_limit_=ncon_ode_step_limit) | |
| 161 | 33 | call get_meiss_coordinates(meiss_ierr) | |
| 162 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 32 times.
|
33 | if (meiss_ierr /= 0) then |
| 163 | 1 | call cleanup_spectre | |
| 164 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (present(ierr)) then |
| 165 | 1 | ierr = lvol | |
| 166 | return | ||
| 167 | end if | ||
| 168 | ✗ | error stop 'SPECTRE canonical construction ODE failed' | |
| 169 | end if | ||
| 170 | 84 | call harvest_meiss_volume(spectre_volumes(lvol)) | |
| 171 | end do | ||
| 172 | |||
| 173 | 19 | call cleanup_meiss | |
| 174 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 17 times.
|
19 | if (present(ierr)) ierr = 0 |
| 175 | end subroutine construct_spectre_coordinates | ||
| 176 | |||
| 177 | |||
| 178 | integer function volume_index(r) result(lvol) | ||
| 179 | real(dp), intent(in) :: r | ||
| 180 | |||
| 181 | 2402947 | lvol = min(int(r) + 1, spectre_mvol) | |
| 182 | 2402947 | lvol = max(1, lvol) | |
| 183 | end function volume_index | ||
| 184 | |||
| 185 | |||
| 186 | 824511 | pure function clamp_to_volume(r, vol) result(rc) | |
| 187 | real(dp), intent(in) :: r | ||
| 188 | type(meiss_volume_t), intent(in) :: vol | ||
| 189 | real(dp) :: rc | ||
| 190 | |||
| 191 | ! Angle-transform evaluation: the batch splines carry no valid | ||
| 192 | ! polynomial extension past the construction domain, so pin the radius | ||
| 193 | ! to the edge. Field evaluation instead goes through band_clamp and | ||
| 194 | ! extend_band (see EDGE_BAND). | ||
| 195 | 824511 | rc = max(vol%rmin, min(vol%rmax, r)) | |
| 196 | end function clamp_to_volume | ||
| 197 | |||
| 198 | |||
| 199 | 5569687 | pure function band_clamp(r, vol) result(rc) | |
| 200 | !> Evaluation radius for the canonical field: pinned EDGE_BAND inside the | ||
| 201 | !> volume; evaluate_spectre_can extends the values over the remaining | ||
| 202 | !> distance (see EDGE_BAND and extend_band). | ||
| 203 | real(dp), intent(in) :: r | ||
| 204 | type(meiss_volume_t), intent(in) :: vol | ||
| 205 | real(dp) :: rc | ||
| 206 | |||
| 207 | 5569687 | rc = max(vol%rmin + EDGE_BAND, min(vol%rmax - EDGE_BAND, r)) | |
| 208 | end function band_clamp | ||
| 209 | |||
| 210 | |||
| 211 | 8480 | pure subroutine extend_band(y, dy, d2y, dr_out, secders) | |
| 212 | !> Radial extension of the harvested field over the signed distance | ||
| 213 | !> dr_out from the band edge (see EDGE_BAND). Per quantity the value | ||
| 214 | !> follows y + y_r*G(dr_out): G = dr_out (linear, C-infinity) for Ath, | ||
| 215 | !> Aph, Bmod; for hth, hph the slope G' blends quadratically from 1 to 0 | ||
| 216 | !> over one band width and stays 0 beyond, so the h-components go over | ||
| 217 | !> smoothly (C1) into radially frozen values at the interface. All first | ||
| 218 | !> derivatives are updated consistently with the extended values; only | ||
| 219 | !> curvatures unavailable without third derivatives stay frozen, which | ||
| 220 | !> degrades the Newton Jacobian alone. | ||
| 221 | real(dp), intent(inout) :: y(5), dy(3, 5), d2y(6, 5) | ||
| 222 | real(dp), intent(in) :: dr_out | ||
| 223 | logical, intent(in) :: secders | ||
| 224 | |||
| 225 | integer :: k | ||
| 226 | real(dp) :: a, s, gval, gp, gpp, y1 | ||
| 227 | |||
| 228 | 8480 | a = abs(dr_out) | |
| 229 | 8480 | s = sign(1.0_dp, dr_out) | |
| 230 | |||
| 231 |
2/2✓ Branch 0 taken 42400 times.
✓ Branch 1 taken 8480 times.
|
50880 | do k = 1, 5 |
| 232 |
2/2✓ Branch 0 taken 16960 times.
✓ Branch 1 taken 25440 times.
|
42400 | if (k == 3 .or. k == 4) then |
| 233 |
2/2✓ Branch 0 taken 14532 times.
✓ Branch 1 taken 2428 times.
|
16960 | if (a < EDGE_BAND) then |
| 234 | 14532 | gval = dr_out - s*0.5_dp*a*a/EDGE_BAND | |
| 235 | 14532 | gp = 1.0_dp - a/EDGE_BAND | |
| 236 | 14532 | gpp = -s/EDGE_BAND | |
| 237 | else | ||
| 238 | 2428 | gval = s*0.5_dp*EDGE_BAND | |
| 239 | 2428 | gp = 0.0_dp | |
| 240 | 2428 | gpp = 0.0_dp | |
| 241 | end if | ||
| 242 | else | ||
| 243 | gval = dr_out | ||
| 244 | gp = 1.0_dp | ||
| 245 | gpp = 0.0_dp | ||
| 246 | end if | ||
| 247 | |||
| 248 | 42400 | y1 = dy(1, k) | |
| 249 | 42400 | y(k) = y(k) + y1*gval | |
| 250 |
2/2✓ Branch 0 taken 42240 times.
✓ Branch 1 taken 160 times.
|
42400 | if (secders) then |
| 251 | 42240 | dy(2, k) = dy(2, k) + d2y(2, k)*gval | |
| 252 | 42240 | dy(3, k) = dy(3, k) + d2y(3, k)*gval | |
| 253 | 42240 | d2y(1, k) = y1*gpp | |
| 254 | 42240 | d2y(2, k) = d2y(2, k)*gp | |
| 255 | 42240 | d2y(3, k) = d2y(3, k)*gp | |
| 256 | end if | ||
| 257 | 50880 | dy(1, k) = y1*gp | |
| 258 | end do | ||
| 259 | 8480 | end subroutine extend_band | |
| 260 | |||
| 261 | |||
| 262 | 5569687 | subroutine evaluate_spectre_can(f, r, th_c, ph_c, mode_secders) | |
| 263 | type(field_can_t), intent(inout) :: f | ||
| 264 | real(dp), intent(in) :: r, th_c, ph_c | ||
| 265 | integer, intent(in) :: mode_secders | ||
| 266 | |||
| 267 | integer :: lvol, k | ||
| 268 | real(dp) :: x(3), dr_out | ||
| 269 | real(dp) :: y(5), dy(3, 5), d2y(6, 5) | ||
| 270 | |||
| 271 |
2/2✓ Branch 1 taken 3168972 times.
✓ Branch 2 taken 2400715 times.
|
5569687 | call count_field_evaluation(mode_secders) |
| 272 | |||
| 273 | 5569687 | lvol = active_volume(r) | |
| 274 |
2/2✓ Branch 0 taken 16709061 times.
✓ Branch 1 taken 5569687 times.
|
22278748 | x = [band_clamp(r, spectre_volumes(lvol)), th_c, ph_c] |
| 275 | 5569687 | dr_out = r - x(1) | |
| 276 | |||
| 277 |
2/2✓ Branch 0 taken 5568970 times.
✓ Branch 1 taken 717 times.
|
5569687 | if (mode_secders > 0) then |
| 278 | call evaluate_batch_splines_3d_der2(spectre_volumes(lvol)%spl_field, & | ||
| 279 | 5568970 | x, y, dy, d2y) | |
| 280 |
2/2✓ Branch 0 taken 8448 times.
✓ Branch 1 taken 5560522 times.
|
5568970 | if (dr_out /= 0.0_dp) then |
| 281 | 8448 | call extend_band(y, dy, d2y, dr_out, .true.) | |
| 282 | end if | ||
| 283 |
2/2✓ Branch 0 taken 5568970 times.
✓ Branch 1 taken 33413820 times.
|
38982790 | f%d2Ath = d2y(:, 1)*A_SCALE |
| 284 |
2/2✓ Branch 0 taken 5568970 times.
✓ Branch 1 taken 33413820 times.
|
38982790 | f%d2Aph = d2y(:, 2)*A_SCALE |
| 285 |
2/2✓ Branch 0 taken 5568970 times.
✓ Branch 1 taken 33413820 times.
|
38982790 | f%d2hth = d2y(:, 3)*H_SCALE |
| 286 |
2/2✓ Branch 0 taken 5568970 times.
✓ Branch 1 taken 33413820 times.
|
38982790 | f%d2hph = d2y(:, 4)*H_SCALE |
| 287 |
2/2✓ Branch 0 taken 33413820 times.
✓ Branch 1 taken 5568970 times.
|
38982790 | f%d2Bmod = d2y(:, 5)*B_SCALE |
| 288 | else | ||
| 289 | call evaluate_batch_splines_3d_der(spectre_volumes(lvol)%spl_field, & | ||
| 290 | 717 | x, y, dy) | |
| 291 |
2/2✓ Branch 0 taken 32 times.
✓ Branch 1 taken 685 times.
|
717 | if (dr_out /= 0.0_dp) then |
| 292 | 32 | call extend_band(y, dy, d2y, dr_out, .false.) | |
| 293 | end if | ||
| 294 | end if | ||
| 295 | |||
| 296 | 5569687 | f%Ath = y(1)*A_SCALE | |
| 297 | 5569687 | f%Aph = y(2)*A_SCALE | |
| 298 | 5569687 | f%hth = y(3)*H_SCALE | |
| 299 | 5569687 | f%hph = y(4)*H_SCALE | |
| 300 | 5569687 | f%Bmod = y(5)*B_SCALE | |
| 301 | |||
| 302 |
2/2✓ Branch 0 taken 5569687 times.
✓ Branch 1 taken 16709061 times.
|
22278748 | f%dAth = dy(:, 1)*A_SCALE |
| 303 |
2/2✓ Branch 0 taken 5569687 times.
✓ Branch 1 taken 16709061 times.
|
22278748 | f%dAph = dy(:, 2)*A_SCALE |
| 304 |
2/2✓ Branch 0 taken 5569687 times.
✓ Branch 1 taken 16709061 times.
|
22278748 | f%dhth = dy(:, 3)*H_SCALE |
| 305 |
2/2✓ Branch 0 taken 5569687 times.
✓ Branch 1 taken 16709061 times.
|
22278748 | f%dhph = dy(:, 4)*H_SCALE |
| 306 |
2/2✓ Branch 0 taken 16709061 times.
✓ Branch 1 taken 5569687 times.
|
22278748 | f%dBmod = dy(:, 5)*B_SCALE |
| 307 | 5569687 | end subroutine evaluate_spectre_can | |
| 308 | |||
| 309 | |||
| 310 | 823854 | subroutine integ_to_ref_spectre(xinteg, xref) | |
| 311 | real(dp), intent(in) :: xinteg(3) | ||
| 312 | real(dp), intent(out) :: xref(3) | ||
| 313 | |||
| 314 | integer :: lvol | ||
| 315 | real(dp) :: x(3), y(2) | ||
| 316 | |||
| 317 |
2/2✓ Branch 0 taken 821942 times.
✓ Branch 1 taken 1912 times.
|
823854 | lvol = active_volume(xinteg(1)) |
| 318 |
2/2✓ Branch 0 taken 2471562 times.
✓ Branch 1 taken 823854 times.
|
3295416 | x = [clamp_to_volume(xinteg(1), spectre_volumes(lvol)), xinteg(2), xinteg(3)] |
| 319 | 823854 | call evaluate_batch_splines_3d(spectre_volumes(lvol)%spl_transform, x, y) | |
| 320 | |||
| 321 | 823854 | xref(1) = xinteg(1) | |
| 322 |
4/4✓ Branch 0 taken 823853 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 46101 times.
✓ Branch 3 taken 777752 times.
|
823854 | xref(2) = modulo(xinteg(2), twopi) |
| 323 |
3/4✓ Branch 0 taken 823854 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1596 times.
✓ Branch 3 taken 822258 times.
|
823854 | xref(3) = modulo(xinteg(3) + y(1), twopi) |
| 324 | 823854 | end subroutine integ_to_ref_spectre | |
| 325 | |||
| 326 | |||
| 327 | 321 | subroutine ref_to_integ_spectre(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 | integer :: lvol, i | ||
| 335 | real(dp) :: x(3), y(2), dy(3, 2), phi, phi_next, target | ||
| 336 | real(dp) :: residual, correction, best_phi, best_residual, denom | ||
| 337 | |||
| 338 |
4/6✓ Branch 0 taken 963 times.
✓ Branch 1 taken 321 times.
✓ Branch 2 taken 963 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 321 times.
|
1284 | if (.not. all(is_finite(xref))) then |
| 339 | ✗ | xinteg = xref | |
| 340 | ✗ | call count_event(EVT_SPECTRE_INVALID_STATE) | |
| 341 | 321 | return | |
| 342 | end if | ||
| 343 | |||
| 344 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 320 times.
|
321 | lvol = active_volume(xref(1)) |
| 345 | 321 | xinteg(1) = xref(1) | |
| 346 |
3/4✓ Branch 0 taken 320 times.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 320 times.
|
321 | xinteg(2) = modulo(xref(2), twopi) |
| 347 |
3/4✓ Branch 0 taken 319 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 319 times.
|
321 | target = modulo(xref(3), twopi) |
| 348 | 321 | phi = target | |
| 349 | 321 | best_phi = phi | |
| 350 | 321 | best_residual = huge(1.0_dp) | |
| 351 | |||
| 352 |
1/2✓ Branch 0 taken 657 times.
✗ Branch 1 not taken.
|
657 | do i = 1, MAX_ITER |
| 353 | x = [clamp_to_volume(xinteg(1), spectre_volumes(lvol)), xinteg(2), & | ||
| 354 |
6/6✓ Branch 0 taken 655 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 654 times.
✓ Branch 4 taken 1971 times.
✓ Branch 5 taken 657 times.
|
2628 | modulo(phi, twopi)] |
| 355 | call evaluate_batch_splines_3d_der(spectre_volumes(lvol)%spl_transform, & | ||
| 356 | 657 | x, y, dy) | |
| 357 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 657 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 657 times.
|
657 | if (.not. is_finite(y(1)) .or. & |
| 358 | .not. is_finite(dy(3, 1))) exit | ||
| 359 | residual = modulo(phi + y(1) - target + 0.5_dp*twopi, twopi) - & | ||
| 360 |
2/4✓ Branch 0 taken 657 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 657 times.
|
657 | 0.5_dp*twopi |
| 361 |
1/2✓ Branch 0 taken 657 times.
✗ Branch 1 not taken.
|
657 | if (abs(residual) < best_residual) then |
| 362 | 657 | best_residual = abs(residual) | |
| 363 | 657 | best_phi = phi | |
| 364 | end if | ||
| 365 | 657 | denom = 1.0_dp + dy(3, 1) | |
| 366 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 657 times.
|
657 | if (abs(denom) <= 100.0_dp*epsilon(1.0_dp)) exit |
| 367 | 657 | correction = max(-0.5_dp*twopi, min(0.5_dp*twopi, residual/denom)) | |
| 368 | 657 | phi_next = phi - correction | |
| 369 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 657 times.
|
657 | if (.not. is_finite(phi_next)) exit |
| 370 | 657 | phi = phi_next | |
| 371 |
2/2✓ Branch 0 taken 321 times.
✓ Branch 1 taken 336 times.
|
657 | if (abs(correction) < TOL) then |
| 372 | 321 | xinteg(3) = phi + twopi*anint((xref(3) - phi)/twopi) | |
| 373 | 321 | return | |
| 374 | end if | ||
| 375 | end do | ||
| 376 | ✗ | xinteg(3) = best_phi + twopi*anint((xref(3) - best_phi)/twopi) | |
| 377 | ✗ | call count_event(EVT_SPECTRE_REF_INVERSE_MAXIT) | |
| 378 | end subroutine ref_to_integ_spectre | ||
| 379 | |||
| 380 | |||
| 381 | 23 | subroutine cleanup_spectre | |
| 382 |
14/16✓ Branch 0 taken 3 times.
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 7 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 7 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 7 times.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 7 times.
✓ Branch 13 taken 7 times.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 3 times.
|
32 | if (allocated(spectre_volumes)) deallocate (spectre_volumes) |
| 383 | 23 | spectre_mvol = 0 | |
| 384 | 23 | end subroutine cleanup_spectre | |
| 385 | |||
| 386 |
1/2✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
|
20 | end module field_can_spectre |
| 387 |