| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module profile_radius | ||
| 2 | !> Map the tracing state to the radius of the plasma profiles. | ||
| 3 | !> | ||
| 4 | !> &profiles are functions of the normalised toroidal flux s. The radial | ||
| 5 | !> state component z(1) is s only for integ_coords TEST, CANFLUX, VMEC and | ||
| 6 | !> BOOZER (VMEC input or Boozer chartmap). MEISS and ALBERT integrate in | ||
| 7 | !> their own radius and map back to the reference coordinates of their | ||
| 8 | !> source field; REFCOORDS integrates in the reference coordinates; SPECTRE | ||
| 9 | !> integrates in the stacked volume chart rho_g in [0, Mvol]. | ||
| 10 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 11 | |||
| 12 | implicit none | ||
| 13 | |||
| 14 | private | ||
| 15 | public :: init_profile_radius, set_profile_radius, profile_s | ||
| 16 | public :: RADIUS_S, RADIUS_RHO_TOR, RADIUS_SPECTRE, RADIUS_NONE | ||
| 17 | |||
| 18 | !> Reference radius is s. | ||
| 19 | integer, parameter :: RADIUS_S = 1 | ||
| 20 | !> Reference radius is rho = rho_scale*sqrt(s). | ||
| 21 | integer, parameter :: RADIUS_RHO_TOR = 2 | ||
| 22 | !> SPECTRE stacked chart rho_g: enclosed toroidal flux / tflux(Mvol). | ||
| 23 | integer, parameter :: RADIUS_SPECTRE = 3 | ||
| 24 | !> No toroidal-flux meaning (geometric or poloidal radius). | ||
| 25 | integer, parameter :: RADIUS_NONE = 4 | ||
| 26 | |||
| 27 | abstract interface | ||
| 28 | function radius_i(z) result(r) | ||
| 29 | import :: dp | ||
| 30 | real(dp), intent(in) :: z(3) | ||
| 31 | real(dp) :: r | ||
| 32 | end function radius_i | ||
| 33 | end interface | ||
| 34 | |||
| 35 | !> Set by set_profile_radius; called from profile_s only when needed. | ||
| 36 | procedure(radius_i), pointer, save :: reference_radius => null() | ||
| 37 | procedure(radius_i), pointer, save :: spectre_s_of_rho => null() | ||
| 38 | |||
| 39 | integer, save :: radius_kind = RADIUS_S | ||
| 40 | logical, save :: map_to_reference = .false. | ||
| 41 | real(dp), save :: rho_scale = 1.0_dp | ||
| 42 | |||
| 43 | contains | ||
| 44 | |||
| 45 | 15 | subroutine set_profile_radius(kind, use_integ_to_ref, rho_lcfs) | |
| 46 | use field_can_meiss, only: integ_to_ref_radius_meiss | ||
| 47 | integer, intent(in) :: kind | ||
| 48 | logical, intent(in) :: use_integ_to_ref | ||
| 49 | real(dp), intent(in), optional :: rho_lcfs | ||
| 50 | |||
| 51 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
10 | if (kind < RADIUS_S .or. kind > RADIUS_NONE) & |
| 52 | ✗ | error stop 'set_profile_radius: unknown radius kind' | |
| 53 | 10 | radius_kind = kind | |
| 54 | 10 | map_to_reference = use_integ_to_ref | |
| 55 | 10 | reference_radius => null() | |
| 56 | 10 | spectre_s_of_rho => null() | |
| 57 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 7 times.
|
10 | if (use_integ_to_ref) reference_radius => integ_to_ref_radius_meiss |
| 58 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9 times.
|
10 | if (kind == RADIUS_SPECTRE) spectre_s_of_rho => spectre_s |
| 59 | 10 | rho_scale = 1.0_dp | |
| 60 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9 times.
|
10 | if (present(rho_lcfs)) then |
| 61 |
2/4✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
|
1 | if (.not. (rho_lcfs > 0.0_dp .and. rho_lcfs <= 1.0_dp)) & |
| 62 | ✗ | error stop 'set_profile_radius: rho_lcfs must be in (0, 1]' | |
| 63 | 1 | rho_scale = rho_lcfs | |
| 64 | end if | ||
| 65 | 10 | end subroutine set_profile_radius | |
| 66 | |||
| 67 | !> Select the mapping for the active integ_coords and field. Call after | ||
| 68 | !> the field and canonical coordinates are initialised. | ||
| 69 | 9 | subroutine init_profile_radius | |
| 70 | 10 | use velo_mod, only: isw_field_type | |
| 71 | use field_can_albert, only: integ_to_ref_radius_albert | ||
| 72 | use magfie_sub, only: TEST, CANFLUX, VMEC, BOOZER, MEISS, ALBERT, & | ||
| 73 | REFCOORDS, SPECTRE, spectre_field | ||
| 74 | use field_can_meiss, only: field_noncan | ||
| 75 | use field_vmec, only: vmec_field_t | ||
| 76 | use field_splined, only: splined_field_t | ||
| 77 | use reference_coordinates, only: ref_coords | ||
| 78 | use params, only: coord_input | ||
| 79 | |||
| 80 | 3 | select case (isw_field_type) | |
| 81 | case (TEST, CANFLUX, VMEC, BOOZER) | ||
| 82 | ! Boozer chartmaps integrate in s = rho^2 (ref_to_integ_boozer_chartmap). | ||
| 83 | 3 | call set_profile_radius(RADIUS_S, .false.) | |
| 84 | case (REFCOORDS) | ||
| 85 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (.not. allocated(ref_coords)) & |
| 86 | ✗ | error stop 'init_profile_radius: no reference coordinates' | |
| 87 | 2 | call from_coordinates(ref_coords, .false.) | |
| 88 | case (MEISS, ALBERT) | ||
| 89 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if (.not. allocated(field_noncan)) & |
| 90 | ✗ | error stop 'init_profile_radius: canonical source field not set' | |
| 91 | 1 | select type (f => field_noncan) | |
| 92 | type is (vmec_field_t) | ||
| 93 | 2 | call set_profile_radius(RADIUS_S, .true.) | |
| 94 | type is (splined_field_t) | ||
| 95 | 1 | call from_coordinates(f%coords, .true.) | |
| 96 | class default | ||
| 97 | ✗ | call set_profile_radius(RADIUS_NONE, .true.) | |
| 98 | end select | ||
| 99 | case (SPECTRE) | ||
| 100 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (.not. allocated(spectre_field)) & |
| 101 | ✗ | error stop 'init_profile_radius: SPECTRE field not set' | |
| 102 | 1 | call set_profile_radius(RADIUS_SPECTRE, .false.) | |
| 103 | case default | ||
| 104 |
4/5✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 1 time.
✗ Branch 4 not taken.
|
9 | error stop 'init_profile_radius: unknown integ_coords' |
| 105 | end select | ||
| 106 |
4/4✓ Branch 0 taken 3 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 2 times.
|
9 | if (map_to_reference .and. isw_field_type == ALBERT) & |
| 107 | 1 | reference_radius => integ_to_ref_radius_albert | |
| 108 | |||
| 109 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
9 | if (radius_kind == RADIUS_NONE .and. .not. profiles_uniform()) then |
| 110 | ✗ | print *, 'init_profile_radius: the radial coordinate of integ_coords =', & | |
| 111 | ✗ | isw_field_type, ' with coord_input = ', trim(coord_input), & | |
| 112 | ✗ | ' has no toroidal-flux meaning.' | |
| 113 | error stop 'radially varying &profiles need a toroidal-flux radius; '// & | ||
| 114 | ✗ | 'use a rho_tor/psi_tor_norm chartmap, VMEC or radially uniform profiles' | |
| 115 | end if | ||
| 116 | |||
| 117 | contains | ||
| 118 | |||
| 119 | 3 | subroutine from_coordinates(cs, use_integ_to_ref) | |
| 120 | 9 | use libneo_coordinates, only: coordinate_system_t, RHO_TOR, & | |
| 121 | PSI_TOR_NORM, chartmap_coordinate_system_t | ||
| 122 | use chartmap_metadata, only: chartmap_metadata_t, read_chartmap_metadata | ||
| 123 | class(coordinate_system_t), intent(in) :: cs | ||
| 124 | logical, intent(in) :: use_integ_to_ref | ||
| 125 | type(chartmap_metadata_t) :: meta | ||
| 126 | |||
| 127 | select type (cs) | ||
| 128 | class is (chartmap_coordinate_system_t) | ||
| 129 | 6 | select case (cs%rho_convention) | |
| 130 | case (RHO_TOR) | ||
| 131 | 3 | call read_chartmap_metadata(coord_input, meta) | |
| 132 | ! Extended VMEC chartmaps put the LCFS at rho_lcfs < 1 and | ||
| 133 | ! scale rho = rho_lcfs*sqrt(s) inside it (libneo chartmap.py). | ||
| 134 | 3 | select case (trim(meta%rho_convention)) | |
| 135 | case ('vmec_extended', 'vmec_to_wall') | ||
| 136 | call set_profile_radius(RADIUS_RHO_TOR, use_integ_to_ref, & | ||
| 137 | ✗ | meta%rho_lcfs) | |
| 138 | case default | ||
| 139 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | call set_profile_radius(RADIUS_RHO_TOR, use_integ_to_ref) |
| 140 | end select | ||
| 141 | case (PSI_TOR_NORM) | ||
| 142 | ✗ | call set_profile_radius(RADIUS_S, use_integ_to_ref) | |
| 143 | case default | ||
| 144 |
1/3✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
3 | call set_profile_radius(RADIUS_NONE, use_integ_to_ref) |
| 145 | end select | ||
| 146 | class default | ||
| 147 | ! VMEC reference coordinates: the radius is s. | ||
| 148 | ✗ | call set_profile_radius(RADIUS_S, use_integ_to_ref) | |
| 149 | end select | ||
| 150 | 3 | end subroutine from_coordinates | |
| 151 | |||
| 152 | end subroutine init_profile_radius | ||
| 153 | |||
| 154 | ✗ | logical function profiles_uniform() | |
| 155 | use simple_profiles, only: get_plasma_params | ||
| 156 | integer, parameter :: n = 101 | ||
| 157 | real(dp) :: p0(5), p(5) | ||
| 158 | integer :: i | ||
| 159 | |||
| 160 | ✗ | call get_plasma_params(0.0_dp, p0(1), p0(2), p0(3), p0(4), p0(5)) | |
| 161 | ✗ | profiles_uniform = .true. | |
| 162 | ✗ | do i = 2, n | |
| 163 | call get_plasma_params(real(i - 1, dp)/real(n - 1, dp), & | ||
| 164 | ✗ | p(1), p(2), p(3), p(4), p(5)) | |
| 165 | ✗ | if (any(p /= p0)) then | |
| 166 | ✗ | profiles_uniform = .false. | |
| 167 | return | ||
| 168 | end if | ||
| 169 | end do | ||
| 170 | end function profiles_uniform | ||
| 171 | |||
| 172 | !> Normalised toroidal flux s at the tracing state z(1:3). | ||
| 173 | !> | ||
| 174 | !> This runs at every collision. It must not use or call procedures of | ||
| 175 | !> modules that pull in IEEE arithmetic: gfortran then wraps the function in | ||
| 176 | !> _gfortran_ieee_procedure_entry/exit (about 0.15 us per call, +11% of a | ||
| 177 | !> collisional run). The expensive branches are therefore reached through | ||
| 178 | !> the procedure pointers set in set_profile_radius. | ||
| 179 | 33 | function profile_s(z) result(s) | |
| 180 | real(dp), intent(in) :: z(3) | ||
| 181 | real(dp) :: s | ||
| 182 | real(dp) :: r | ||
| 183 | |||
| 184 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 24 times.
|
33 | if (map_to_reference) then |
| 185 | 9 | r = reference_radius(z) | |
| 186 | else | ||
| 187 | 24 | r = z(1) | |
| 188 | end if | ||
| 189 | 9 | select case (radius_kind) | |
| 190 | case (RADIUS_RHO_TOR) | ||
| 191 | 9 | s = (r/rho_scale)**2 | |
| 192 | case (RADIUS_SPECTRE) | ||
| 193 | ! SPECTRE integrates in rho_g itself, so r = z(1) here. | ||
| 194 | 5 | s = spectre_s_of_rho(z) | |
| 195 | case default | ||
| 196 | ! RADIUS_NONE is accepted only for radially uniform profiles. | ||
| 197 |
3/3✓ Branch 0 taken 9 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 19 times.
|
33 | s = r |
| 198 | end select | ||
| 199 | 33 | end function profile_s | |
| 200 | |||
| 201 | 5 | function spectre_s(z) result(s) | |
| 202 | use magfie_sub, only: spectre_field | ||
| 203 | use spectre_basis, only: eval_spectre_toroidal_flux_fraction | ||
| 204 | real(dp), intent(in) :: z(3) | ||
| 205 | real(dp) :: s | ||
| 206 | real(dp) :: rho_g, rho, s_local, fraction, derivative, inner | ||
| 207 | integer :: lvol, mvol, ierr | ||
| 208 | |||
| 209 | 5 | rho_g = z(1) | |
| 210 | 5 | mvol = spectre_field%data%Mvol | |
| 211 | 5 | rho = min(max(rho_g, 0.0_dp), real(mvol, dp)) | |
| 212 | 5 | lvol = min(int(rho) + 1, mvol) | |
| 213 | 5 | s_local = 2.0_dp*(rho - real(lvol - 1, dp)) - 1.0_dp | |
| 214 | call eval_spectre_toroidal_flux_fraction(spectre_field%data, lvol, s_local, & | ||
| 215 | 5 | fraction, derivative, ierr) | |
| 216 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if (ierr /= 0) error stop 'profile_s: SPECTRE toroidal flux unavailable' |
| 217 | ! tflux holds the cumulative toroidal flux enclosed by each interface. | ||
| 218 | 5 | inner = 0.0_dp | |
| 219 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
|
5 | if (lvol > 1) inner = spectre_field%data%tflux(lvol - 1) |
| 220 | s = (inner + fraction*(spectre_field%data%tflux(lvol) - inner)) & | ||
| 221 | 5 | /spectre_field%data%tflux(mvol) | |
| 222 | 5 | end function spectre_s | |
| 223 | |||
| 224 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
|
6 | end module profile_radius |
| 225 |