| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module field_can_albert | ||
| 2 | !> Albert canonical coordinates extending Meiss coordinates. | ||
| 3 | !> | ||
| 4 | !> This module builds on field_can_meiss to provide Albert-type canonical | ||
| 5 | !> coordinates (psi, theta_c, phi_c) where psi is proportional to the | ||
| 6 | !> poloidal magnetic flux A_theta, giving the simple form A = (0, psi, Aphi). | ||
| 7 | !> | ||
| 8 | !> The Albert coordinates differ from Meiss in the radial coordinate: | ||
| 9 | !> - Meiss: r = sqrt(s), fixed radial coordinate | ||
| 10 | !> - Albert: psi ~ A_theta, flux-based radial coordinate | ||
| 11 | !> | ||
| 12 | !> The transformation is built by: | ||
| 13 | !> 1. Computing Meiss coordinates via field_can_meiss | ||
| 14 | !> 2. Regridding onto a uniform psi grid via psi_transform | ||
| 15 | !> | ||
| 16 | !> Key subroutines: | ||
| 17 | !> - get_albert_coordinates: Build full transformation (calls Meiss first) | ||
| 18 | !> - evaluate_albert: Fast splined field evaluation in Albert coords | ||
| 19 | !> - integ_to_ref_albert: Convert integrator coords to reference (s,th,ph) | ||
| 20 | !> | ||
| 21 | !> The Albert form simplifies the symplectic integrator since dA_theta/dr = 0. | ||
| 22 | |||
| 23 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 24 | use interpolate, only: & | ||
| 25 | BatchSplineData3D, construct_batch_splines_3d, & | ||
| 26 | evaluate_batch_splines_3d, evaluate_batch_splines_3d_der, & | ||
| 27 | evaluate_batch_splines_3d_der2 | ||
| 28 | use field_can_base, only: count_field_evaluation, field_can_t | ||
| 29 | use field_can_meiss, only: xmin, xmax, n_r, n_th, n_phi, order, periodic, twopi, & | ||
| 30 | get_grid_point, & | ||
| 31 | init_albert => init_meiss, init_transformation, spline_transformation, & | ||
| 32 | init_canonical_field_components | ||
| 33 | use psi_transform, only: grid_r_to_psi | ||
| 34 | |||
| 35 | implicit none | ||
| 36 | |||
| 37 | ! For splining psi | ||
| 38 | real(dp) :: psi_inner, psi_outer | ||
| 39 | real(dp), dimension(:, :, :), allocatable :: psi_of_x | ||
| 40 | real(dp), dimension(:), allocatable :: psi_grid | ||
| 41 | logical :: dpsi_dr_positive | ||
| 42 | |||
| 43 | ! For splining field components over canonical coordinates | ||
| 44 | real(dp), dimension(:, :, :), allocatable :: r_of_xc, & | ||
| 45 | Aph_of_xc, hth_of_xc, hph_of_xc, Bmod_of_xc | ||
| 46 | |||
| 47 | ! Batch spline for r_of_xc transformation (1 component: r) | ||
| 48 | type(BatchSplineData3D) :: spl_r_batch | ||
| 49 | |||
| 50 | ! Batch spline for optimized field evaluation (4 components: Aphi, hth, hph, Bmod) | ||
| 51 | type(BatchSplineData3D) :: spl_albert_batch | ||
| 52 | |||
| 53 | real(dp) :: Ath_norm | ||
| 54 | |||
| 55 | contains | ||
| 56 | |||
| 57 | 4 | subroutine evaluate_albert(f, r, th_c, ph_c, mode_secders) | |
| 58 | type(field_can_t), intent(inout) :: f | ||
| 59 | real(dp), intent(in) :: r, th_c, ph_c | ||
| 60 | integer, intent(in) :: mode_secders | ||
| 61 | |||
| 62 | real(dp) :: x(3) | ||
| 63 | |||
| 64 | 4 | call count_field_evaluation(mode_secders) | |
| 65 | |||
| 66 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 4 times.
|
16 | x = [r, th_c, ph_c] |
| 67 | |||
| 68 | 4 | f%Ath = Ath_norm*x(1) | |
| 69 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 4 times.
|
16 | f%dAth = [Ath_norm, 0d0, 0d0] |
| 70 | |||
| 71 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (mode_secders > 0) then |
| 72 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 4 times.
|
28 | f%d2Ath = 0d0 |
| 73 | 4 | call evaluate_albert_batch_der2(f, x) | |
| 74 | 4 | return | |
| 75 | end if | ||
| 76 | |||
| 77 | ✗ | call evaluate_albert_batch_der(f, x) | |
| 78 | end subroutine evaluate_albert | ||
| 79 | |||
| 80 | 2 | subroutine integ_to_ref_albert(xinteg, xref) | |
| 81 | use field_can_meiss, only: integ_to_ref_meiss | ||
| 82 | |||
| 83 | real(dp), intent(in) :: xinteg(3) | ||
| 84 | real(dp), intent(out) :: xref(3) | ||
| 85 | real(dp) :: xmeiss(3), y_batch(1) | ||
| 86 | |||
| 87 | 2 | call evaluate_batch_splines_3d(spl_r_batch, xinteg, y_batch) | |
| 88 | 2 | xmeiss(1) = y_batch(1) ! r component | |
| 89 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
|
6 | xmeiss(2:3) = xinteg(2:3) |
| 90 | 2 | call integ_to_ref_meiss(xmeiss, xref) | |
| 91 | 2 | end subroutine integ_to_ref_albert | |
| 92 | |||
| 93 | !> Reference radius of integ_to_ref_albert without the angular splines. | ||
| 94 | 3 | function integ_to_ref_radius_albert(xinteg) result(r_ref) | |
| 95 | 2 | use field_can_meiss, only: integ_to_ref_radius_meiss | |
| 96 | |||
| 97 | real(dp), intent(in) :: xinteg(3) | ||
| 98 | real(dp) :: r_ref | ||
| 99 | real(dp) :: xmeiss(3), y_batch(1) | ||
| 100 | |||
| 101 | 3 | call evaluate_batch_splines_3d(spl_r_batch, xinteg, y_batch) | |
| 102 | 3 | xmeiss(1) = y_batch(1) | |
| 103 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 3 times.
|
9 | xmeiss(2:3) = xinteg(2:3) |
| 104 | 3 | r_ref = integ_to_ref_radius_meiss(xmeiss) | |
| 105 | 3 | end function integ_to_ref_radius_albert | |
| 106 | |||
| 107 | 7 | subroutine ref_to_integ_albert(xref, xinteg) | |
| 108 | 3 | use field_can_meiss, only: ref_to_integ_meiss, spl_field_batch | |
| 109 | |||
| 110 | real(dp), intent(in) :: xref(3) | ||
| 111 | real(dp), intent(out) :: xinteg(3) | ||
| 112 | |||
| 113 | real(dp) :: Ath, xmeiss(3), y_batch_local(5) | ||
| 114 | |||
| 115 | 7 | call ref_to_integ_meiss(xref, xmeiss) | |
| 116 | 7 | call evaluate_batch_splines_3d(spl_field_batch, xmeiss, y_batch_local) | |
| 117 | 7 | Ath = y_batch_local(1) ! Extract Ath component | |
| 118 | 7 | xinteg(1) = Ath/Ath_norm | |
| 119 |
2/2✓ Branch 0 taken 14 times.
✓ Branch 1 taken 7 times.
|
21 | xinteg(2:3) = xmeiss(2:3) |
| 120 | 7 | end subroutine ref_to_integ_albert | |
| 121 | |||
| 122 | 3 | subroutine get_albert_coordinates | |
| 123 | 3 | call cleanup_albert | |
| 124 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 125 | print *, 'field_can_meiss.init_transformation' | ||
| 126 | #endif | ||
| 127 | 3 | call init_transformation | |
| 128 | |||
| 129 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 130 | print *, 'field_can_meiss.spline_transformation' | ||
| 131 | #endif | ||
| 132 | 3 | call spline_transformation | |
| 133 | |||
| 134 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 135 | print *, 'field_can_meiss.init_canonical_field_components' | ||
| 136 | #endif | ||
| 137 | 3 | call init_canonical_field_components | |
| 138 | |||
| 139 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 140 | print *, 'field_can_albert.init_splines_with_psi' | ||
| 141 | #endif | ||
| 142 | 3 | call init_splines_with_psi | |
| 143 | 7 | end subroutine get_albert_coordinates | |
| 144 | |||
| 145 | 78 | subroutine cleanup_albert | |
| 146 | use interpolate, only: destroy_batch_splines_3d | ||
| 147 | |||
| 148 | 78 | call destroy_batch_splines_3d(spl_r_batch) | |
| 149 | 78 | call destroy_batch_splines_3d(spl_albert_batch) | |
| 150 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 77 times.
|
78 | if (allocated(psi_of_x)) deallocate(psi_of_x) |
| 151 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 77 times.
|
78 | if (allocated(psi_grid)) deallocate(psi_grid) |
| 152 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 77 times.
|
78 | if (allocated(r_of_xc)) deallocate(r_of_xc) |
| 153 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 77 times.
|
78 | if (allocated(Aph_of_xc)) deallocate(Aph_of_xc) |
| 154 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 77 times.
|
78 | if (allocated(hth_of_xc)) deallocate(hth_of_xc) |
| 155 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 77 times.
|
78 | if (allocated(hph_of_xc)) deallocate(hph_of_xc) |
| 156 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 77 times.
|
78 | if (allocated(Bmod_of_xc)) deallocate(Bmod_of_xc) |
| 157 | 78 | end subroutine cleanup_albert | |
| 158 | |||
| 159 | 3 | subroutine init_splines_with_psi | |
| 160 | use psi_transform, only: grid_r_to_psi | ||
| 161 | use field_can_meiss, only: spl_field_batch | ||
| 162 | 3 | real(dp), dimension(:, :, :, :), allocatable :: y_batch | |
| 163 | real(dp), dimension(:, :, :), allocatable :: Aphi_grid, hth_grid, hph_grid, Bmod_grid | ||
| 164 | real(dp) :: x_grid(3), y_batch_temp(5) | ||
| 165 | integer :: i_r, i_th, i_phi, dims(3) | ||
| 166 | |||
| 167 | allocate ( & | ||
| 168 | ✗ | r_of_xc(n_r, n_th, n_phi), & | |
| 169 | ✗ | Aph_of_xc(n_r, n_th, n_phi), & | |
| 170 | ✗ | hth_of_xc(n_r, n_th, n_phi), & | |
| 171 | ✗ | hph_of_xc(n_r, n_th, n_phi), & | |
| 172 | ✗ | Bmod_of_xc(n_r, n_th, n_phi) & | |
| 173 |
47/94✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 3 times.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 3 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 3 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 3 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 3 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 3 times.
✓ Branch 26 taken 3 times.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✓ Branch 29 taken 3 times.
✓ Branch 30 taken 3 times.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✓ Branch 33 taken 3 times.
✓ Branch 34 taken 3 times.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✓ Branch 37 taken 3 times.
✗ Branch 39 not taken.
✓ Branch 40 taken 3 times.
✗ Branch 41 not taken.
✓ Branch 42 taken 3 times.
✗ Branch 44 not taken.
✓ Branch 45 taken 3 times.
✓ Branch 46 taken 3 times.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✓ Branch 49 taken 3 times.
✓ Branch 50 taken 3 times.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✓ Branch 53 taken 3 times.
✓ Branch 54 taken 3 times.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✓ Branch 57 taken 3 times.
✗ Branch 59 not taken.
✓ Branch 60 taken 3 times.
✗ Branch 61 not taken.
✓ Branch 62 taken 3 times.
✗ Branch 64 not taken.
✓ Branch 65 taken 3 times.
✓ Branch 66 taken 3 times.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✓ Branch 69 taken 3 times.
✓ Branch 70 taken 3 times.
✗ Branch 71 not taken.
✗ Branch 72 not taken.
✓ Branch 73 taken 3 times.
✓ Branch 74 taken 3 times.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✓ Branch 77 taken 3 times.
✗ Branch 79 not taken.
✓ Branch 80 taken 3 times.
✗ Branch 81 not taken.
✓ Branch 82 taken 3 times.
✗ Branch 84 not taken.
✓ Branch 85 taken 3 times.
✓ Branch 86 taken 3 times.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✓ Branch 89 taken 3 times.
✓ Branch 90 taken 3 times.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✓ Branch 93 taken 3 times.
✓ Branch 94 taken 3 times.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✓ Branch 97 taken 3 times.
✗ Branch 99 not taken.
✓ Branch 100 taken 3 times.
✗ Branch 101 not taken.
✓ Branch 102 taken 3 times.
|
51 | ) |
| 174 | |||
| 175 | 3 | call init_psi_grid | |
| 176 | |||
| 177 | ! For Albert coordinates, we need to reconstruct the field components from Meiss | ||
| 178 | ! This requires grid evaluation - let's compute them explicitly | ||
| 179 |
34/68✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 3 times.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 3 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 3 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 3 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 3 times.
✓ Branch 24 taken 3 times.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 3 times.
✓ Branch 28 taken 3 times.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✓ Branch 31 taken 3 times.
✓ Branch 32 taken 3 times.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✓ Branch 35 taken 3 times.
✗ Branch 37 not taken.
✓ Branch 38 taken 3 times.
✗ Branch 40 not taken.
✓ Branch 41 taken 3 times.
✓ Branch 42 taken 3 times.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✓ Branch 45 taken 3 times.
✓ Branch 46 taken 3 times.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✓ Branch 49 taken 3 times.
✓ Branch 50 taken 3 times.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✓ Branch 53 taken 3 times.
✗ Branch 55 not taken.
✓ Branch 56 taken 3 times.
✗ Branch 58 not taken.
✓ Branch 59 taken 3 times.
✓ Branch 60 taken 3 times.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✓ Branch 63 taken 3 times.
✓ Branch 64 taken 3 times.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✓ Branch 67 taken 3 times.
✓ Branch 68 taken 3 times.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✓ Branch 71 taken 3 times.
✗ Branch 73 not taken.
✓ Branch 74 taken 3 times.
|
42 | allocate(Aphi_grid(n_r,n_th,n_phi), hth_grid(n_r,n_th,n_phi), hph_grid(n_r,n_th,n_phi), Bmod_grid(n_r,n_th,n_phi)) |
| 180 | |||
| 181 | ! Evaluate Meiss batch spline on grid to get field components | ||
| 182 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 37 times.
|
40 | do i_phi = 1, n_phi |
| 183 |
2/2✓ Branch 0 taken 468 times.
✓ Branch 1 taken 37 times.
|
508 | do i_th = 1, n_th |
| 184 |
2/2✓ Branch 0 taken 6120 times.
✓ Branch 1 taken 468 times.
|
6625 | do i_r = 1, n_r |
| 185 | x_grid = [xmin(1) + (i_r - 1)*(xmax(1) - xmin(1))/(n_r - 1), & | ||
| 186 | xmin(2) + (i_th - 1)*(xmax(2) - xmin(2))/(n_th - 1), & | ||
| 187 |
2/2✓ Branch 0 taken 18360 times.
✓ Branch 1 taken 6120 times.
|
24480 | xmin(3) + (i_phi - 1)*(xmax(3) - xmin(3))/(n_phi - 1)] |
| 188 | 6120 | call evaluate_batch_splines_3d(spl_field_batch, x_grid, y_batch_temp) | |
| 189 | 6120 | Aphi_grid(i_r, i_th, i_phi) = y_batch_temp(2) ! Aph component | |
| 190 | 6120 | hth_grid(i_r, i_th, i_phi) = y_batch_temp(3) ! hth component | |
| 191 | 6120 | hph_grid(i_r, i_th, i_phi) = y_batch_temp(4) ! hph component | |
| 192 | 6588 | Bmod_grid(i_r, i_th, i_phi) = y_batch_temp(5) ! Bmod component | |
| 193 | end do | ||
| 194 | end do | ||
| 195 | end do | ||
| 196 | |||
| 197 | ! Center Aphi around zero | ||
| 198 |
12/12✓ Branch 0 taken 37 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 468 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 6120 times.
✓ Branch 5 taken 468 times.
✓ Branch 6 taken 37 times.
✓ Branch 7 taken 3 times.
✓ Branch 8 taken 468 times.
✓ Branch 9 taken 37 times.
✓ Branch 10 taken 6120 times.
✓ Branch 11 taken 468 times.
|
13253 | Aphi_grid = Aphi_grid - 0.5d0*sum(Aphi_grid)/real(n_r*n_th*n_phi, dp) |
| 199 | |||
| 200 | call grid_r_to_psi(xmin(1), xmax(1), psi_inner, psi_outer, psi_of_x, & | ||
| 201 | Aphi_grid, hth_grid, hph_grid, Bmod_grid, r_of_xc, Aph_of_xc, & | ||
| 202 | 3 | hth_of_xc, hph_of_xc, Bmod_of_xc) | |
| 203 | |||
| 204 | ! Construct batch spline for r_of_xc (1 component: r) | ||
| 205 | 3 | block | |
| 206 | real(dp), dimension(:, :, :, :), allocatable :: y_r_batch | ||
| 207 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
|
12 | dims = shape(r_of_xc) |
| 208 |
10/20✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 3 times.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 3 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 3 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 3 times.
|
15 | allocate (y_r_batch(dims(1), dims(2), dims(3), 1)) |
| 209 |
6/6✓ Branch 0 taken 37 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 468 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 6120 times.
✓ Branch 5 taken 468 times.
|
6628 | y_r_batch(:, :, :, 1) = r_of_xc |
| 210 | call construct_batch_splines_3d([psi_inner, xmin(2), xmin(3)], & | ||
| 211 |
5/6✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3 times.
|
24 | [psi_outer, xmax(2), xmax(3)], y_r_batch, order, periodic, spl_r_batch) |
| 212 | end block | ||
| 213 | |||
| 214 | ! Construct batch spline for 4 Albert field components: [Aphi, hth, hph, Bmod] | ||
| 215 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
|
12 | dims = shape(Aph_of_xc) |
| 216 |
11/22✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 3 times.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 3 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 3 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 3 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 3 times.
|
18 | allocate (y_batch(dims(1), dims(2), dims(3), 4)) |
| 217 | |||
| 218 |
6/6✓ Branch 0 taken 37 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 468 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 6120 times.
✓ Branch 5 taken 468 times.
|
6628 | y_batch(:, :, :, 1) = Aph_of_xc |
| 219 |
6/6✓ Branch 0 taken 37 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 468 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 6120 times.
✓ Branch 5 taken 468 times.
|
6628 | y_batch(:, :, :, 2) = hth_of_xc |
| 220 |
6/6✓ Branch 0 taken 37 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 468 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 6120 times.
✓ Branch 5 taken 468 times.
|
6628 | y_batch(:, :, :, 3) = hph_of_xc |
| 221 |
6/6✓ Branch 0 taken 37 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 468 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 6120 times.
✓ Branch 5 taken 468 times.
|
6628 | y_batch(:, :, :, 4) = Bmod_of_xc |
| 222 | |||
| 223 | call construct_batch_splines_3d([psi_inner, xmin(2), xmin(3)], & | ||
| 224 |
4/4✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 3 times.
|
21 | [psi_outer, xmax(2), xmax(3)], y_batch, order, periodic, spl_albert_batch) |
| 225 |
5/10✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
|
3 | end subroutine init_splines_with_psi |
| 226 | |||
| 227 | 3 | subroutine init_psi_grid | |
| 228 | 3 | use field_can_meiss, only: spl_field_batch | |
| 229 | real(dp) :: x(3), y_batch_local(5) | ||
| 230 | integer :: i_r, i_th, i_phi | ||
| 231 | |||
| 232 |
18/36✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 3 times.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 3 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 3 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 3 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 3 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 3 times.
✓ Branch 27 taken 3 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 3 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 3 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✓ Branch 34 taken 3 times.
✗ Branch 36 not taken.
✓ Branch 37 taken 3 times.
✗ Branch 39 not taken.
✓ Branch 40 taken 3 times.
|
21 | allocate (psi_of_x(n_r, n_th, n_phi), psi_grid(n_r)) |
| 233 | |||
| 234 | ! Evaluate Meiss batch spline to get Ath (component 1) on grid | ||
| 235 |
2/2✓ Branch 0 taken 37 times.
✓ Branch 1 taken 3 times.
|
40 | do i_phi = 1, n_phi |
| 236 |
2/2✓ Branch 0 taken 468 times.
✓ Branch 1 taken 37 times.
|
508 | do i_th = 1, n_th |
| 237 |
2/2✓ Branch 0 taken 6120 times.
✓ Branch 1 taken 468 times.
|
6625 | do i_r = 1, n_r |
| 238 | x = [xmin(1) + (i_r - 1)*(xmax(1) - xmin(1))/(n_r - 1), & | ||
| 239 | xmin(2) + (i_th - 1)*(xmax(2) - xmin(2))/(n_th - 1), & | ||
| 240 |
2/2✓ Branch 0 taken 18360 times.
✓ Branch 1 taken 6120 times.
|
24480 | xmin(3) + (i_phi - 1)*(xmax(3) - xmin(3))/(n_phi - 1)] |
| 241 | 6120 | call evaluate_batch_splines_3d(spl_field_batch, x, y_batch_local) | |
| 242 | 6588 | psi_of_x(i_r, i_th, i_phi) = y_batch_local(1) ! Ath component | |
| 243 | end do | ||
| 244 | end do | ||
| 245 | end do | ||
| 246 | |||
| 247 |
12/16✓ Branch 0 taken 37 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 468 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 6120 times.
✓ Branch 5 taken 468 times.
✓ Branch 6 taken 6117 times.
✓ Branch 7 taken 3 times.
✓ Branch 8 taken 31 times.
✓ Branch 9 taken 6086 times.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
6628 | Ath_norm = sign(maxval(abs(psi_of_x)), psi_of_x(n_r, n_th/2, n_phi/2)) |
| 248 |
6/6✓ Branch 0 taken 37 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 468 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 6120 times.
✓ Branch 5 taken 468 times.
|
6628 | psi_of_x = psi_of_x/Ath_norm |
| 249 | |||
| 250 | ! Here we use the "safe side" approach (new grid is fully within the old grid). | ||
| 251 | ! For the risky approach (old grid within the new grid) exchange | ||
| 252 | ! "minval" and "maxval". | ||
| 253 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | if (psi_of_x(n_r, n_th/2, n_phi/2) > psi_of_x(1, n_th/2, n_phi/2)) then |
| 254 | 3 | dpsi_dr_positive = .true. | |
| 255 |
10/14✓ Branch 0 taken 37 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 468 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 465 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 9 times.
✓ Branch 7 taken 456 times.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 3 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
508 | psi_inner = maxval(psi_of_x(1, :, :)) |
| 256 |
10/14✓ Branch 0 taken 37 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 468 times.
✓ Branch 3 taken 37 times.
✓ Branch 4 taken 465 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 29 times.
✓ Branch 7 taken 436 times.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 3 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
508 | psi_outer = minval(psi_of_x(n_r, :, :)) |
| 257 | else | ||
| 258 | ✗ | dpsi_dr_positive = .false. | |
| 259 | ✗ | psi_inner = maxval(psi_of_x(n_r, :, :)) | |
| 260 | ✗ | psi_outer = minval(psi_of_x(1, :, :)) | |
| 261 | end if | ||
| 262 | |||
| 263 |
2/2✓ Branch 0 taken 31 times.
✓ Branch 1 taken 3 times.
|
34 | do i_r = 1, n_r |
| 264 | 34 | psi_grid(i_r) = psi_inner + (psi_outer - psi_inner)*(i_r - 1)/(n_r - 1) | |
| 265 | end do | ||
| 266 | 3 | end subroutine init_psi_grid | |
| 267 | |||
| 268 | ✗ | subroutine magfie_albert(x, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) | |
| 269 | ! Computes magnetic field and derivatives with bmod in units of the magnetic code | ||
| 270 | ! | ||
| 271 | ! Input parameters: | ||
| 272 | ! formal: x - array of canonicalized coordinates r, th, ph | ||
| 273 | ! Output parameters: | ||
| 274 | ! formal: bmod - magnetic field module | ||
| 275 | ! sqrtg - metric determinant | ||
| 276 | ! bder - covariant components of (grad B)/B | ||
| 277 | ! hcovar - covariant components of \bB/B | ||
| 278 | ! hctrvr - contravariant components of \bB/B | ||
| 279 | ! hcurl - contravariant components of curl (\bB/B) | ||
| 280 | real(dp), intent(in) :: x(3) | ||
| 281 | real(dp), intent(out) :: bmod, sqrtg | ||
| 282 | real(dp), dimension(3), intent(out) :: bder, hcovar, hctrvr, hcurl | ||
| 283 | |||
| 284 | type(field_can_t) :: f | ||
| 285 | real(dp) :: sqrtg_bmod | ||
| 286 | |||
| 287 | ✗ | call evaluate_albert(f, x(1), x(2), x(3), 0) | |
| 288 | |||
| 289 | ✗ | bmod = f%Bmod | |
| 290 | |||
| 291 | ✗ | sqrtg_bmod = f%hph*Ath_norm - f%hth*f%dAph(1) | |
| 292 | ✗ | sqrtg = sqrtg_bmod/bmod | |
| 293 | ✗ | bder = f%dBmod/bmod | |
| 294 | |||
| 295 | ✗ | hcovar(1) = 0.d0 | |
| 296 | ✗ | hcovar(2) = f%hth | |
| 297 | ✗ | hcovar(3) = f%hph | |
| 298 | |||
| 299 | ✗ | hctrvr(1) = f%dAph(2)/sqrtg_bmod | |
| 300 | ✗ | hctrvr(2) = -f%dAph(1)/sqrtg_bmod | |
| 301 | ✗ | hctrvr(3) = Ath_norm/sqrtg_bmod | |
| 302 | |||
| 303 | ✗ | hcurl(1) = (f%dhph(2) - f%dhth(3))/sqrtg | |
| 304 | ✗ | hcurl(2) = -f%dhph(1)/sqrtg | |
| 305 | ✗ | hcurl(3) = f%dhth(1)/sqrtg | |
| 306 |
0/4✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
3 | end subroutine magfie_albert |
| 307 | |||
| 308 | ! Batch evaluation helper for first derivatives | ||
| 309 | ✗ | subroutine evaluate_albert_batch_der(f, x) | |
| 310 | type(field_can_t), intent(inout) :: f | ||
| 311 | real(dp), intent(in) :: x(3) | ||
| 312 | |||
| 313 | real(dp) :: y_batch(4), dy_batch(3, 4) | ||
| 314 | |||
| 315 | ✗ | call evaluate_batch_splines_3d_der(spl_albert_batch, x, y_batch, dy_batch) | |
| 316 | |||
| 317 | ! Unpack results: order is [Aphi, hth, hph, Bmod] | ||
| 318 | ✗ | f%Aph = y_batch(1) | |
| 319 | ✗ | f%hth = y_batch(2) | |
| 320 | ✗ | f%hph = y_batch(3) | |
| 321 | ✗ | f%Bmod = y_batch(4) | |
| 322 | |||
| 323 | ✗ | f%dAph = dy_batch(:, 1) | |
| 324 | ✗ | f%dhth = dy_batch(:, 2) | |
| 325 | ✗ | f%dhph = dy_batch(:, 3) | |
| 326 | ✗ | f%dBmod = dy_batch(:, 4) | |
| 327 | ✗ | end subroutine evaluate_albert_batch_der | |
| 328 | |||
| 329 | ! Batch evaluation helper for second derivatives | ||
| 330 | 4 | subroutine evaluate_albert_batch_der2(f, x) | |
| 331 | type(field_can_t), intent(inout) :: f | ||
| 332 | real(dp), intent(in) :: x(3) | ||
| 333 | |||
| 334 | real(dp) :: y_batch(4), dy_batch(3, 4), d2y_batch(6, 4) | ||
| 335 | |||
| 336 | 4 | call evaluate_batch_splines_3d_der2(spl_albert_batch, x, y_batch, dy_batch, d2y_batch) | |
| 337 | |||
| 338 | ! Unpack results: order is [Aphi, hth, hph, Bmod] | ||
| 339 | 4 | f%Aph = y_batch(1) | |
| 340 | 4 | f%hth = y_batch(2) | |
| 341 | 4 | f%hph = y_batch(3) | |
| 342 | 4 | f%Bmod = y_batch(4) | |
| 343 | |||
| 344 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 12 times.
|
16 | f%dAph = dy_batch(:, 1) |
| 345 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 12 times.
|
16 | f%dhth = dy_batch(:, 2) |
| 346 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 12 times.
|
16 | f%dhph = dy_batch(:, 3) |
| 347 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 12 times.
|
16 | f%dBmod = dy_batch(:, 4) |
| 348 | |||
| 349 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 24 times.
|
28 | f%d2Aph = d2y_batch(:, 1) |
| 350 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 24 times.
|
28 | f%d2hth = d2y_batch(:, 2) |
| 351 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 24 times.
|
28 | f%d2hph = d2y_batch(:, 3) |
| 352 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 4 times.
|
28 | f%d2Bmod = d2y_batch(:, 4) |
| 353 | 4 | end subroutine evaluate_albert_batch_der2 | |
| 354 | |||
| 355 | end module field_can_albert | ||
| 356 |