| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module samplers | ||
| 2 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 3 | use util | ||
| 4 | |||
| 5 | implicit none | ||
| 6 | |||
| 7 | character(len=*), parameter :: START_FILE = 'start.dat' | ||
| 8 | character(len=*), parameter :: START_FILE_ANTS = 'start_ants.dat' | ||
| 9 | character(len=*), parameter :: START_FILE_BATCH = 'batch.dat' | ||
| 10 | |||
| 11 | ! Interface ################################ | ||
| 12 | INTERFACE sample | ||
| 13 | MODULE PROCEDURE sample_read | ||
| 14 | MODULE PROCEDURE sample_surface_fieldline | ||
| 15 | MODULE PROCEDURE sample_grid | ||
| 16 | MODULE PROCEDURE sample_volume_single | ||
| 17 | MODULE PROCEDURE sample_random_batch | ||
| 18 | MODULE PROCEDURE sample_points_ants | ||
| 19 | END INTERFACE sample | ||
| 20 | |||
| 21 | contains | ||
| 22 | ! Functions ################################# | ||
| 23 | 57 | subroutine init_starting_surf | |
| 24 | use alpha_lifetime_sub, only: integrate_mfl_can | ||
| 25 | use params, only: dphi, nper, npoiper, phibeg, thetabeg, volstart, & | ||
| 26 | xstart, sbeg, bmin, bmax, bmod00 | ||
| 27 | |||
| 28 | integer :: ierr = 0 | ||
| 29 | 114 | real(dp), dimension(npoiper*nper) :: bstart | |
| 30 | |||
| 31 |
4/4✓ Branch 0 taken 57 times.
✓ Branch 1 taken 2200280 times.
✓ Branch 2 taken 6600840 times.
✓ Branch 3 taken 2200280 times.
|
8801177 | xstart = 0.d0 |
| 32 |
2/2✓ Branch 0 taken 2200280 times.
✓ Branch 1 taken 57 times.
|
2200337 | bstart = 0.d0 |
| 33 |
2/2✓ Branch 0 taken 2200280 times.
✓ Branch 1 taken 57 times.
|
2200337 | volstart = 0.d0 |
| 34 | |||
| 35 | ! For VMEC-backed runs the driver calls this while VMEC magfie is active, | ||
| 36 | ! so xstart can be copied directly to reference-coordinate zstart. The | ||
| 37 | ! volstart integral gives volume-weighted sampling on this one surface. | ||
| 38 | call integrate_mfl_can( & | ||
| 39 | npoiper*nper, dphi, sbeg(1), phibeg, thetabeg, & | ||
| 40 | 57 | xstart, bstart, volstart, bmod00, ierr) | |
| 41 | |||
| 42 |
1/2✓ Branch 0 taken 57 times.
✗ Branch 1 not taken.
|
57 | if (ierr .ne. 0) then |
| 43 | ✗ | print *, 'starting field line has points outside the chamber' | |
| 44 | ✗ | stop | |
| 45 | end if | ||
| 46 | |||
| 47 | ! maximum value of B module: | ||
| 48 |
9/10✓ Branch 0 taken 63 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 55 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2200272 times.
✓ Branch 7 taken 57 times.
✓ Branch 8 taken 4284 times.
✓ Branch 9 taken 2195988 times.
|
2200394 | bmax = maxval(bstart) |
| 49 |
9/10✓ Branch 0 taken 63 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 55 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2200272 times.
✓ Branch 7 taken 57 times.
✓ Branch 8 taken 340 times.
✓ Branch 9 taken 2199932 times.
|
2200394 | bmin = minval(bstart) |
| 50 | |||
| 51 | 57 | print *, 'bmod00 = ', bmod00, 'bmin = ', bmin, 'bmax = ', bmax | |
| 52 | 57 | end subroutine init_starting_surf | |
| 53 | |||
| 54 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
|
15 | subroutine load_starting_points(zstart, filename) |
| 55 | real(dp), dimension(:, :), intent(inout) :: zstart | ||
| 56 | character(len=*), intent(in) :: filename | ||
| 57 | integer :: ipart | ||
| 58 | |||
| 59 | 15 | open (1, file=filename, recl=1024) | |
| 60 |
2/2✓ Branch 0 taken 115 times.
✓ Branch 1 taken 15 times.
|
130 | do ipart = 1, size(zstart, 2) |
| 61 | 130 | read (1, *) zstart(:, ipart) | |
| 62 | end do | ||
| 63 | 15 | close (1) | |
| 64 | 57 | end subroutine load_starting_points | |
| 65 | |||
| 66 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
|
52 | subroutine save_starting_points(zstart) |
| 67 | real(dp), dimension(:, :), intent(in) :: zstart | ||
| 68 | integer :: ipart | ||
| 69 | |||
| 70 | 52 | open (1, file=START_FILE, recl=1024) | |
| 71 |
2/2✓ Branch 0 taken 700697 times.
✓ Branch 1 taken 52 times.
|
700749 | do ipart = 1, size(zstart, 2) |
| 72 | 700749 | write (1, *) zstart(:, ipart) | |
| 73 | end do | ||
| 74 | 52 | close (1) | |
| 75 | 52 | end subroutine save_starting_points | |
| 76 | |||
| 77 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 14 times.
|
14 | subroutine sample_read(zstart, filename) |
| 78 | real(dp), dimension(:, :), intent(inout) :: zstart | ||
| 79 | character(len=*), intent(in) :: filename | ||
| 80 | |||
| 81 | 14 | call load_starting_points(zstart, filename) | |
| 82 | 14 | end subroutine | |
| 83 | |||
| 84 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | subroutine sample_read_integration(zstart, filename) |
| 85 | ! Read (s, theta_B, zeta_B, |v|/v0, v_parallel/|v|) and convert only | ||
| 86 | ! the position to SIMPLE's stored reference-coordinate convention. | ||
| 87 | ! This permits exact reuse of externally generated Boozer particles. | ||
| 88 | use field_can_mod, only: integ_to_ref | ||
| 89 | |||
| 90 | real(dp), dimension(:, :), intent(inout) :: zstart | ||
| 91 | character(len=*), intent(in) :: filename | ||
| 92 | real(dp) :: xinteg(3) | ||
| 93 | integer :: ipart | ||
| 94 | |||
| 95 | 1 | call load_starting_points(zstart, filename) | |
| 96 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
|
4 | do ipart = 1, size(zstart, 2) |
| 97 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
|
12 | xinteg = zstart(1:3, ipart) |
| 98 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
4 | call integ_to_ref(xinteg, zstart(1:3, ipart)) |
| 99 | end do | ||
| 100 | 1 | end subroutine sample_read_integration | |
| 101 | |||
| 102 | ! Samplers ################################ | ||
| 103 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | subroutine sample_volume_single(zstart, s_inner, s_outer) |
| 104 | 1 | use params, only: isw_field_type, num_surf | |
| 105 | use field_can_mod, only: integ_to_ref | ||
| 106 | |||
| 107 | real(dp), intent(in) :: s_inner | ||
| 108 | real(dp), intent(in) :: s_outer | ||
| 109 | real(dp), parameter :: s_min = 0.01d0 | ||
| 110 | real(dp) :: tmp_rand, s_lo, s_hi | ||
| 111 | real(dp) :: r, vartheta, varphi | ||
| 112 | real(dp), dimension(:, :), intent(inout) :: zstart | ||
| 113 | integer :: ipart | ||
| 114 | |||
| 115 | ! If user wants to do volume with 0 or 1 surfaces, | ||
| 116 | ! we "add" the constraints, therefore having 2 surfaces. | ||
| 117 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (2 /= num_surf) then |
| 118 | 1 | num_surf = 2 | |
| 119 | end if | ||
| 120 | |||
| 121 | ! Clamp lower bound to s_min to avoid axis singularity | ||
| 122 | 1 | s_lo = max(s_inner, s_min) | |
| 123 | 1 | s_hi = max(s_outer, s_min) | |
| 124 | |||
| 125 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | do ipart = 1, size(zstart, 2) |
| 126 | 2 | call random_number(tmp_rand) | |
| 127 | 2 | r = tmp_rand*(s_hi - s_lo) + s_lo | |
| 128 | |||
| 129 | 2 | call random_number(tmp_rand) | |
| 130 | 2 | vartheta = twopi*tmp_rand | |
| 131 | 2 | call random_number(tmp_rand) | |
| 132 | 2 | varphi = twopi*tmp_rand | |
| 133 | ! we store starting points in reference coordinates: | ||
| 134 |
4/8✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
8 | call integ_to_ref([r, vartheta, varphi], zstart(1:3, ipart)) |
| 135 | ! normalized velocity module z(4) = v / v_0: | ||
| 136 | 2 | zstart(4, ipart) = 1.d0 | |
| 137 | ! starting pitch z(5)=v_\parallel / v: | ||
| 138 | 2 | call random_number(tmp_rand) | |
| 139 | 3 | zstart(5, ipart) = 2.d0*(tmp_rand - 0.5d0) | |
| 140 | end do | ||
| 141 | |||
| 142 | 1 | call save_starting_points(zstart) | |
| 143 | |||
| 144 | 1 | end subroutine sample_volume_single | |
| 145 | |||
| 146 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
|
9 | subroutine sample_fusion_birth(zstart, s_inner, s_outer, integration_coords) |
| 147 | use params, only: fusion_grid, num_surf, facE_al, n_d, n_e, & | ||
| 148 | 1 | am1, am2, Z1, Z2, sbeg | |
| 149 | use fusion_source, only: build_birth_cdf, draw_birth, alpha_energy_joule | ||
| 150 | use field_can_mod, only: integ_to_ref | ||
| 151 | real(dp), intent(out) :: zstart(:, :) | ||
| 152 | real(dp), value, intent(in) :: s_inner, s_outer | ||
| 153 | logical, intent(in) :: integration_coords | ||
| 154 | 9 | real(dp), allocatable :: cdf(:) | |
| 155 | real(dp) :: power, x(3) | ||
| 156 | integer :: ipart, unit, n | ||
| 157 | |||
| 158 |
4/6✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 8 times.
|
9 | if (facE_al /= 1.0_dp .or. n_d /= 4.0_dp .or. n_e /= 2) & |
| 159 | 1 | error stop 'fusion births require 3.5 MeV alphas: facE_al=1, n_d=4, n_e=2' | |
| 160 | if (am1 /= 2.0_dp .or. am2 /= 3.0_dp .or. & | ||
| 161 |
4/8✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 8 times.
|
8 | Z1 /= 1.0_dp .or. Z2 /= 1.0_dp) & |
| 162 | ✗ | error stop 'fusion births require species 1=D and species 2=T' | |
| 163 | 8 | n = size(zstart, 2) | |
| 164 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | if (n < 1) error stop 'fusion births require at least one marker' |
| 165 | 8 | call build_birth_cdf(s_inner, s_outer, fusion_grid, cdf, power) | |
| 166 |
2/2✓ Branch 1 taken 700008 times.
✓ Branch 2 taken 8 times.
|
700016 | do ipart = 1, n |
| 167 |
2/8✗ Branch 0 not taken.
✓ Branch 1 taken 700008 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 700008 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
700008 | call draw_birth(s_inner, s_outer, fusion_grid, cdf, zstart(:, ipart)) |
| 168 |
2/2✓ Branch 0 taken 100000 times.
✓ Branch 1 taken 600008 times.
|
700016 | if (integration_coords) then |
| 169 |
2/2✓ Branch 0 taken 300000 times.
✓ Branch 1 taken 100000 times.
|
400000 | x = zstart(1:3, ipart) |
| 170 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 100000 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 100000 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
100000 | call integ_to_ref(x, zstart(1:3, ipart)) |
| 171 | end if | ||
| 172 | end do | ||
| 173 | 8 | num_surf = 2 | |
| 174 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 8 times.
|
24 | sbeg(1:2) = [s_inner, s_outer] |
| 175 | 8 | call save_starting_points(zstart) | |
| 176 | 8 | open (newunit=unit, file='alpha_source.dat', status='replace') | |
| 177 | 8 | write (unit, '(a)') '# Equal-weight D-T births; rows of start.dat are markers.' | |
| 178 | 8 | write (unit, '(a)') '# Midpoint source cells, uniform positions within cells.' | |
| 179 | 8 | write (unit, '(a)') '# alpha_power_W uses 3.5 MeV per alpha, not 3.52 MeV.' | |
| 180 | 8 | write (unit, '(a,es24.16)') 'alpha_energy_MeV ', & | |
| 181 | 16 | alpha_energy_joule/1.602176634e-13_dp | |
| 182 | 8 | write (unit, '(a,es24.16)') 'alpha_power_W ', power | |
| 183 | 8 | write (unit, '(a,es24.16)') 'alpha_birth_rate_per_s ', power/alpha_energy_joule | |
| 184 | 8 | write (unit, '(a,es24.16)') 'marker_rate_per_s ', power/alpha_energy_joule/n | |
| 185 | 8 | write (unit, '(a,es24.16)') 'marker_power_W ', power/n | |
| 186 | 8 | write (unit, '(a,i0)') 'marker_count ', n | |
| 187 | 8 | write (unit, '(a,2es24.16)') 's_bounds ', s_inner, s_outer | |
| 188 | 8 | write (unit, '(a,3i8)') 'fusion_grid ', fusion_grid | |
| 189 | 8 | close (unit) | |
| 190 | 8 | print *, 'Total alpha power [W] = ', power | |
| 191 | 8 | print *, 'Equal marker power [W] = ', power/n, ' (alpha_source.dat)' | |
| 192 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
8 | end subroutine sample_fusion_birth |
| 193 | |||
| 194 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
|
18 | subroutine sample_spectre_surface(zstart) |
| 195 | !> Place markers on the SPECTRE control surface rho_g = sbeg(1), drawn | ||
| 196 | !> uniformly over the real surface: (theta, zeta) grid nodes are picked | ||
| 197 | !> with probability proportional to the surface Jacobian |e_theta x e_zeta|. | ||
| 198 | !> The chart is the integration coordinate, so zstart is (rho_g, theta, | ||
| 199 | !> zeta) directly. Also sets bmin/bmax/bmod00 for pitch classification. | ||
| 200 | 9 | use params, only: sbeg, bmin, bmax, bmod00, & | |
| 201 | spectre_sbeg_is_toroidal_flux | ||
| 202 | use magfie_sub, only: magfie, spectre_field | ||
| 203 | use binsrc_sub, only: binsrc | ||
| 204 | |||
| 205 | real(dp), dimension(:, :), intent(inout) :: zstart | ||
| 206 | |||
| 207 | integer, parameter :: ntheta = 64, nzeta = 64 | ||
| 208 | integer :: ngrid, it, iz, k, ipart, ierr | ||
| 209 | real(dp) :: rho_g, th, ze, area, wsum, bsum, bmod, xi | ||
| 210 | real(dp) :: u(3), e_cov(3, 3), cross(3) | ||
| 211 | real(dp) :: sqrtg, bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 212 | 18 | real(dp), allocatable :: theta_g(:), zeta_g(:), wcum(:) | |
| 213 | |||
| 214 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
|
18 | if (.not. allocated(spectre_field)) & |
| 215 | ✗ | error stop 'sample_spectre_surface: spectre_field not set' | |
| 216 | |||
| 217 | 18 | ngrid = ntheta*nzeta | |
| 218 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 17 times.
|
18 | if (spectre_sbeg_is_toroidal_flux) then |
| 219 | 1 | call spectre_field%axis_rho_from_toroidal_flux(sbeg(1), rho_g, ierr) | |
| 220 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (ierr /= 0) & |
| 221 | ✗ | error stop 'sample_spectre_surface: target flux outside axis volume' | |
| 222 | else | ||
| 223 | 17 | rho_g = sbeg(1) | |
| 224 | end if | ||
| 225 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 18 times.
|
18 | allocate (theta_g(ngrid), zeta_g(ngrid), wcum(ngrid)) |
| 226 | |||
| 227 | 18 | wsum = 0.0d0 | |
| 228 | 18 | bsum = 0.0d0 | |
| 229 | 18 | bmin = huge(1.0d0) | |
| 230 | 18 | bmax = -huge(1.0d0) | |
| 231 | 18 | k = 0 | |
| 232 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 1152 times.
|
1170 | do it = 1, ntheta |
| 233 |
2/2✓ Branch 0 taken 73728 times.
✓ Branch 1 taken 1152 times.
|
74898 | do iz = 1, nzeta |
| 234 | 73728 | k = k + 1 | |
| 235 | 73728 | th = twopi*(real(it, dp) - 0.5d0)/real(ntheta, dp) | |
| 236 | 73728 | ze = twopi*(real(iz, dp) - 0.5d0)/real(nzeta, dp) | |
| 237 |
2/2✓ Branch 0 taken 221184 times.
✓ Branch 1 taken 73728 times.
|
294912 | u = [rho_g, th, ze] |
| 238 | 73728 | call spectre_field%coords%covariant_basis(u, e_cov) | |
| 239 | 73728 | cross(1) = e_cov(2, 2)*e_cov(3, 3) - e_cov(3, 2)*e_cov(2, 3) | |
| 240 | 73728 | cross(2) = e_cov(3, 2)*e_cov(1, 3) - e_cov(1, 2)*e_cov(3, 3) | |
| 241 | 73728 | cross(3) = e_cov(1, 2)*e_cov(2, 3) - e_cov(2, 2)*e_cov(1, 3) | |
| 242 |
2/2✓ Branch 0 taken 221184 times.
✓ Branch 1 taken 73728 times.
|
294912 | area = sqrt(sum(cross**2)) |
| 243 | 73728 | call magfie(u, bmod, sqrtg, bder, hcovar, hctrvr, hcurl) | |
| 244 | 73728 | bmin = min(bmin, bmod) | |
| 245 | 73728 | bmax = max(bmax, bmod) | |
| 246 | 73728 | bsum = bsum + bmod | |
| 247 | 73728 | wsum = wsum + area | |
| 248 | 73728 | theta_g(k) = th | |
| 249 | 73728 | zeta_g(k) = ze | |
| 250 | 74880 | wcum(k) = wsum | |
| 251 | end do | ||
| 252 | end do | ||
| 253 |
2/2✓ Branch 0 taken 73728 times.
✓ Branch 1 taken 18 times.
|
73746 | wcum = wcum/wsum |
| 254 | 18 | bmod00 = bsum/real(ngrid, dp) | |
| 255 | |||
| 256 |
2/2✓ Branch 0 taken 476 times.
✓ Branch 1 taken 18 times.
|
494 | do ipart = 1, size(zstart, 2) |
| 257 | 476 | call random_number(xi) | |
| 258 | 476 | call binsrc(wcum, 1, ngrid, xi, k) | |
| 259 | 476 | zstart(1, ipart) = rho_g | |
| 260 | 476 | zstart(2, ipart) = theta_g(k) | |
| 261 | 476 | zstart(3, ipart) = zeta_g(k) | |
| 262 | 476 | zstart(4, ipart) = 1.0d0 | |
| 263 | 476 | call random_number(xi) | |
| 264 | 494 | zstart(5, ipart) = 2.0d0*(xi - 0.5d0) | |
| 265 | end do | ||
| 266 | |||
| 267 | 18 | print *, 'sample_spectre_surface: rho_g = ', rho_g, ' target = ', sbeg(1), & | |
| 268 | 18 | ' flux_label = ', spectre_sbeg_is_toroidal_flux, ' bmod00 = ', bmod00, & | |
| 269 | 36 | ' bmin = ', bmin, ' bmax = ', bmax | |
| 270 | 18 | call save_starting_points(zstart) | |
| 271 |
3/6✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
|
18 | end subroutine sample_spectre_surface |
| 272 | |||
| 273 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | subroutine init_spectre_start_bounds(zstart) |
| 274 | 18 | use params, only: bmin, bmax, bmod00 | |
| 275 | use magfie_sub, only: magfie | ||
| 276 | |||
| 277 | real(dp), dimension(:, :), intent(in) :: zstart | ||
| 278 | |||
| 279 | integer, parameter :: nrho = 8, ntheta = 64, nzeta = 64 | ||
| 280 | integer :: ir, it, iz, nr | ||
| 281 | real(dp) :: rho_lo, rho_hi, rho, theta, zeta, bmod, bsum | ||
| 282 | real(dp) :: sqrtg, bder(3), hcovar(3), hctrvr(3), hcurl(3) | ||
| 283 | |||
| 284 |
6/10✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 3 times.
✓ Branch 7 taken 1 time.
✓ Branch 8 taken 1 time.
✓ Branch 9 taken 2 times.
|
5 | rho_lo = minval(zstart(1, :)) |
| 285 |
5/10✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 3 times.
✓ Branch 7 taken 1 time.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
|
5 | rho_hi = maxval(zstart(1, :)) |
| 286 | nr = merge(1, nrho, abs(rho_hi - rho_lo) <= & | ||
| 287 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | epsilon(1.0d0)*max(1.0d0, abs(rho_lo), abs(rho_hi))) |
| 288 | 1 | bmin = huge(1.0d0) | |
| 289 | 1 | bmax = -huge(1.0d0) | |
| 290 | 1 | bsum = 0.0d0 | |
| 291 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 time.
|
9 | do ir = 1, nr |
| 292 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
8 | if (nr == 1) then |
| 293 | rho = rho_lo | ||
| 294 | else | ||
| 295 | 8 | rho = rho_lo + (rho_hi - rho_lo)*real(ir - 1, dp)/real(nr - 1, dp) | |
| 296 | end if | ||
| 297 |
2/2✓ Branch 0 taken 512 times.
✓ Branch 1 taken 8 times.
|
521 | do it = 1, ntheta |
| 298 | 512 | theta = twopi*(real(it, dp) - 0.5d0)/real(ntheta, dp) | |
| 299 |
2/2✓ Branch 0 taken 32768 times.
✓ Branch 1 taken 512 times.
|
33288 | do iz = 1, nzeta |
| 300 | 32768 | zeta = twopi*(real(iz, dp) - 0.5d0)/real(nzeta, dp) | |
| 301 | call magfie([rho, theta, zeta], bmod, sqrtg, bder, hcovar, & | ||
| 302 |
2/2✓ Branch 0 taken 98304 times.
✓ Branch 1 taken 32768 times.
|
131072 | hctrvr, hcurl) |
| 303 | 32768 | bmin = min(bmin, bmod) | |
| 304 | 32768 | bmax = max(bmax, bmod) | |
| 305 | 33280 | bsum = bsum + bmod | |
| 306 | end do | ||
| 307 | end do | ||
| 308 | end do | ||
| 309 | 1 | bmod00 = bsum/real(nr*ntheta*nzeta, dp) | |
| 310 | |||
| 311 | 1 | print *, 'init_spectre_start_bounds: rho = ', rho_lo, rho_hi, & | |
| 312 | 2 | ' bmod00 = ', bmod00, ' bmin = ', bmin, ' bmax = ', bmax | |
| 313 | 1 | end subroutine init_spectre_start_bounds | |
| 314 | |||
| 315 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 20 times.
|
20 | subroutine sample_surface_fieldline(zstart) |
| 316 | real(dp), dimension(:, :), intent(inout) :: zstart | ||
| 317 | |||
| 318 | 20 | call sample_surface_fieldline_impl(zstart, .false.) | |
| 319 | 1 | end subroutine sample_surface_fieldline | |
| 320 | |||
| 321 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | subroutine sample_surface_fieldline_from_integ(zstart) |
| 322 | real(dp), dimension(:, :), intent(inout) :: zstart | ||
| 323 | |||
| 324 | 3 | call sample_surface_fieldline_impl(zstart, .true.) | |
| 325 | 3 | end subroutine sample_surface_fieldline_from_integ | |
| 326 | |||
| 327 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23 times.
|
23 | subroutine sample_surface_fieldline_impl(zstart, xstart_is_integ_coords) |
| 328 | use params, only: volstart, ibins, xstart, npoiper, nper | ||
| 329 | use binsrc_sub, only: binsrc | ||
| 330 | use field_can_mod, only: integ_to_ref | ||
| 331 | |||
| 332 | real(dp), dimension(:, :), intent(inout) :: zstart | ||
| 333 | logical, intent(in) :: xstart_is_integ_coords | ||
| 334 | |||
| 335 | real(dp) :: xi | ||
| 336 | integer :: ipart, i | ||
| 337 | |||
| 338 |
2/2✓ Branch 0 taken 199 times.
✓ Branch 1 taken 23 times.
|
222 | do ipart = 1, size(zstart, 2) |
| 339 | 199 | call random_number(xi) | |
| 340 | 199 | call binsrc(volstart, 1, npoiper*nper, xi, i) | |
| 341 | 199 | ibins = i | |
| 342 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 175 times.
|
199 | if (xstart_is_integ_coords) then |
| 343 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
24 | call integ_to_ref(xstart(:, i), zstart(1:3, ipart)) |
| 344 | else | ||
| 345 |
2/2✓ Branch 0 taken 525 times.
✓ Branch 1 taken 175 times.
|
700 | zstart(1:3, ipart) = xstart(:, i) |
| 346 | end if | ||
| 347 | 199 | zstart(4, ipart) = 1.d0 ! normalized velocity module z(4) = v / v_0 | |
| 348 | 199 | call random_number(xi) | |
| 349 | 222 | zstart(5, ipart) = 2.d0*(xi - 0.5d0) ! starting pitch z(5)=v_\parallel / v | |
| 350 | end do | ||
| 351 | |||
| 352 | 23 | call save_starting_points(zstart) | |
| 353 | |||
| 354 | 23 | end subroutine sample_surface_fieldline_impl | |
| 355 | |||
| 356 | 1 | subroutine sample_grid(zstart, grid_density, xstart_is_integ_coords) | |
| 357 | use params, only: ntestpart, zstart_dim1, zend, times_lost, & | ||
| 358 | orbit_exit_code, boundary_event_radial_residual, & | ||
| 359 | 23 | boundary_event_time_width, trap_par, perp_inv, iclass, sbeg | |
| 360 | use util, only: pi | ||
| 361 | use field_can_mod, only: integ_to_ref | ||
| 362 | |||
| 363 | real(dp), dimension(:, :), allocatable, intent(inout) :: zstart | ||
| 364 | real(dp), intent(in) :: grid_density | ||
| 365 | logical, intent(in), optional :: xstart_is_integ_coords | ||
| 366 | real(dp) :: xi, xsize_real | ||
| 367 | real(dp) :: xinteg(3) | ||
| 368 | integer :: ngrid, ipart, jpart, lidx | ||
| 369 | logical :: convert_surface_starts | ||
| 370 | |||
| 371 | 1 | convert_surface_starts = .false. | |
| 372 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (present(xstart_is_integ_coords)) then |
| 373 | 1 | convert_surface_starts = xstart_is_integ_coords | |
| 374 | end if | ||
| 375 | |||
| 376 | 1 | xsize_real = (2*pi)*grid_density !angle density | |
| 377 | 1 | ngrid = int((1/grid_density) - 1) | |
| 378 | 1 | ntestpart = ngrid**2 !number of total angle points | |
| 379 | |||
| 380 | ! Resize particle coord. arrays and result memory. | ||
| 381 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (allocated(zstart)) deallocate (zstart) |
| 382 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (allocated(zend)) deallocate (zend) |
| 383 |
16/32✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 time.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 time.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 time.
✗ Branch 18 not taken.
✓ Branch 19 taken 1 time.
✗ Branch 21 not taken.
✓ Branch 22 taken 1 time.
✓ Branch 23 taken 1 time.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 1 time.
✓ Branch 27 taken 1 time.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 30 taken 1 time.
✗ Branch 32 not taken.
✓ Branch 33 taken 1 time.
✗ Branch 35 not taken.
✓ Branch 36 taken 1 time.
|
6 | allocate (zstart(zstart_dim1, ntestpart), zend(zstart_dim1, ntestpart)) |
| 384 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (allocated(times_lost)) deallocate (times_lost) |
| 385 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (allocated(orbit_exit_code)) deallocate (orbit_exit_code) |
| 386 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (allocated(boundary_event_radial_residual)) & |
| 387 | 1 | deallocate (boundary_event_radial_residual) | |
| 388 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (allocated(boundary_event_time_width)) deallocate (boundary_event_time_width) |
| 389 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (allocated(trap_par)) deallocate (trap_par) |
| 390 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (allocated(perp_inv)) deallocate (perp_inv) |
| 391 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (allocated(iclass)) deallocate (iclass) |
| 392 | ✗ | allocate(times_lost(ntestpart), orbit_exit_code(ntestpart), & | |
| 393 | ✗ | boundary_event_radial_residual(ntestpart), & | |
| 394 | ✗ | boundary_event_time_width(ntestpart), & | |
| 395 |
41/82✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 time.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 time.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 time.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 time.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 time.
✓ Branch 18 taken 1 time.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 time.
✗ Branch 21 not taken.
✓ Branch 22 taken 1 time.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 1 time.
✗ Branch 27 not taken.
✓ Branch 28 taken 1 time.
✗ Branch 29 not taken.
✓ Branch 30 taken 1 time.
✗ Branch 32 not taken.
✓ Branch 33 taken 1 time.
✓ Branch 34 taken 1 time.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✓ Branch 37 taken 1 time.
✗ Branch 39 not taken.
✓ Branch 40 taken 1 time.
✗ Branch 41 not taken.
✓ Branch 42 taken 1 time.
✗ Branch 44 not taken.
✓ Branch 45 taken 1 time.
✓ Branch 46 taken 1 time.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✓ Branch 49 taken 1 time.
✗ Branch 51 not taken.
✓ Branch 52 taken 1 time.
✗ Branch 53 not taken.
✓ Branch 54 taken 1 time.
✗ Branch 56 not taken.
✓ Branch 57 taken 1 time.
✓ Branch 58 taken 1 time.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✓ Branch 61 taken 1 time.
✗ Branch 63 not taken.
✓ Branch 64 taken 1 time.
✗ Branch 66 not taken.
✓ Branch 67 taken 1 time.
✗ Branch 69 not taken.
✓ Branch 70 taken 1 time.
✓ Branch 71 taken 1 time.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✓ Branch 74 taken 1 time.
✗ Branch 76 not taken.
✓ Branch 77 taken 1 time.
✗ Branch 79 not taken.
✓ Branch 80 taken 1 time.
✗ Branch 82 not taken.
✓ Branch 83 taken 1 time.
✓ Branch 84 taken 1 time.
✗ Branch 85 not taken.
✓ Branch 86 taken 1 time.
✗ Branch 87 not taken.
✓ Branch 88 taken 1 time.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✓ Branch 91 taken 1 time.
✗ Branch 93 not taken.
✓ Branch 94 taken 1 time.
✗ Branch 96 not taken.
✓ Branch 97 taken 1 time.
|
11 | trap_par(ntestpart), perp_inv(ntestpart), iclass(3,ntestpart)) |
| 396 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9 times.
|
10 | orbit_exit_code = 0 |
| 397 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9 times.
|
10 | boundary_event_radial_residual = -1d0 |
| 398 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9 times.
|
10 | boundary_event_time_width = -1d0 |
| 399 | |||
| 400 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
|
4 | do ipart = 1, ngrid |
| 401 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
|
13 | do jpart = 1, ngrid |
| 402 | 9 | lidx = (jpart - 1)*ngrid + ipart | |
| 403 |
2/2✓ Branch 0 taken 27 times.
✓ Branch 1 taken 9 times.
|
36 | xinteg = [sbeg(1), xsize_real*ipart, xsize_real*jpart] |
| 404 |
1/2✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
|
9 | if (convert_surface_starts) then |
| 405 | 9 | call integ_to_ref(xinteg, zstart(1:3, lidx)) | |
| 406 | else | ||
| 407 | ✗ | zstart(1:3, lidx) = xinteg | |
| 408 | end if | ||
| 409 | 9 | zstart(4, lidx) = 1.d0 ! normalized velocity module z(4) = v / v_0 | |
| 410 | 9 | call random_number(xi) | |
| 411 | 12 | zstart(5, lidx) = 2.d0*(xi - 0.5d0) ! starting pitch z(5)=v_\parallel / v | |
| 412 | end do | ||
| 413 | end do | ||
| 414 | |||
| 415 | 1 | call save_starting_points(zstart) | |
| 416 | |||
| 417 | 1 | end subroutine sample_grid | |
| 418 | |||
| 419 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | subroutine sample_random_batch(zstart, reuse_existing) |
| 420 | ! init_batch owns the persisted source-row indices; batch.dat contains | ||
| 421 | ! indices, never particle coordinates. Read selected rows from start.dat. | ||
| 422 | 1 | use params, only: idx | |
| 423 | |||
| 424 | real(dp), intent(inout) :: zstart(:, :) | ||
| 425 | logical, intent(in) :: reuse_existing | ||
| 426 | 6 | real(dp) :: source_row(size(zstart, 1)) | |
| 427 | integer :: unit, source_index, selected, status | ||
| 428 | |||
| 429 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if (.not. allocated(idx)) then |
| 430 | ✗ | if (reuse_existing) error stop 'No selected indices for reused batch' | |
| 431 | ✗ | call load_starting_points(zstart, START_FILE) | |
| 432 | return | ||
| 433 | end if | ||
| 434 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if (size(idx) /= size(zstart, 2)) error stop 'Batch size mismatch' |
| 435 | 3 | open (newunit=unit, file=START_FILE, status='old', action='read') | |
| 436 | 3 | selected = 1 | |
| 437 |
1/2✓ Branch 0 taken 36 times.
✗ Branch 1 not taken.
|
36 | do source_index = 1, idx(size(idx)) |
| 438 | 36 | read (unit, *, iostat=status) source_row | |
| 439 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
|
36 | if (status /= 0) error stop 'Cannot read selected start.dat source row' |
| 440 |
2/2✓ Branch 0 taken 15 times.
✓ Branch 1 taken 21 times.
|
36 | if (source_index /= idx(selected)) cycle |
| 441 |
2/2✓ Branch 0 taken 75 times.
✓ Branch 1 taken 15 times.
|
90 | zstart(:, selected) = source_row |
| 442 | 15 | selected = selected + 1 | |
| 443 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 3 times.
|
15 | if (selected > size(idx)) exit |
| 444 | end do | ||
| 445 | 3 | close (unit) | |
| 446 | 3 | end subroutine sample_random_batch | |
| 447 | |||
| 448 | ✗ | subroutine sample_points_ants(use_special_ants_file) | |
| 449 | 3 | use parse_ants, only: process_line | |
| 450 | use get_can_sub, only: vmec_to_can | ||
| 451 | use params, only: ntestpart, zstart ! ANTS sampler uses global zstart | ||
| 452 | |||
| 453 | logical, intent(in) :: use_special_ants_file | ||
| 454 | |||
| 455 | integer, parameter :: maxlen = 4096 | ||
| 456 | character(len=maxlen) :: line | ||
| 457 | real(8) :: v_par, v_perp, u, v, s | ||
| 458 | real(8) :: th, ph | ||
| 459 | integer :: ipart | ||
| 460 | |||
| 461 | ✗ | do ipart = 1, ntestpart | |
| 462 | ✗ | if (use_special_ants_file) then | |
| 463 | ✗ | open (1, file=START_FILE_ANTS, recl=1024) | |
| 464 | ✗ | read (1, '(A)') line | |
| 465 | ✗ | close (1) | |
| 466 | else | ||
| 467 | ✗ | open (1, file=START_FILE, recl=1024) | |
| 468 | ✗ | read (1, '(A)') line | |
| 469 | ✗ | close (1) | |
| 470 | end if | ||
| 471 | |||
| 472 | ✗ | call process_line(line, v_par, v_perp, u, v, s) | |
| 473 | ! In the test case, u runs from 0 to 1 and v from 0 to 4 | ||
| 474 | ✗ | th = 2d0*pi*u | |
| 475 | ✗ | ph = 2d0*pi*v/4d0 | |
| 476 | ✗ | zstart(1, ipart) = s | |
| 477 | ✗ | zstart(2, ipart) = th | |
| 478 | ✗ | zstart(3, ipart) = ph | |
| 479 | ✗ | zstart(4, ipart) = 1.d0 | |
| 480 | ✗ | zstart(5, ipart) = v_par/sqrt(v_par**2 + v_perp**2) | |
| 481 | end do | ||
| 482 | ✗ | end subroutine sample_points_ants | |
| 483 | |||
| 484 | end module samplers | ||
| 485 |