| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | ! | ||
| 2 | module detect_oneline_mod | ||
| 3 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 4 | |||
| 5 | logical :: prop=.true. | ||
| 6 | integer, parameter :: nfp_dim=3, ipermin=10 | ||
| 7 | integer :: nfp_max,iper,igroup | ||
| 8 | double precision :: pi,twopi | ||
| 9 | integer, dimension(:), allocatable :: iret | ||
| 10 | double precision, dimension(:,:), allocatable :: fprs | ||
| 11 | !> Continuous scores from every banana tip collected so far. J_parallel | ||
| 12 | !> variation is the RMS relative drift rate per poloidal precession turn. | ||
| 13 | !> Rotation drift compares the mean tip-map rotation number in the first | ||
| 14 | !> and second halves of the trace. Sample counts distinguish a physical | ||
| 15 | !> zero from an unresolved score. | ||
| 16 | real(dp) :: jpar_variation_rate, rotation_number_drift, precession_turns | ||
| 17 | integer :: jpar_sample_count, rotation_half_count | ||
| 18 | !> Legacy margins retained to cross-check the published classifiers. | ||
| 19 | !> Production metrics consume the continuous scores above, not thresholds. | ||
| 20 | !> score_status: 0 unresolved; 1 resolved with a valid topology margin; | ||
| 21 | !> 2 early stochastic exit; 3 resolved but the ideal-orbit class came from | ||
| 22 | !> the recurrence test, which forms no margin. jpar_spread is valid for | ||
| 23 | !> status 1 and 3, topology_margin only for status 1. | ||
| 24 | double precision :: jpar_spread, jpar_ref, topology_margin | ||
| 25 | integer :: tip_count | ||
| 26 | integer :: score_status | ||
| 27 | !$omp threadprivate(prop, nfp_max, iper, igroup, iret, fprs, & | ||
| 28 | !$omp jpar_variation_rate, rotation_number_drift, & | ||
| 29 | !$omp precession_turns, jpar_sample_count, rotation_half_count, & | ||
| 30 | !$omp jpar_spread, jpar_ref, topology_margin, score_status, & | ||
| 31 | !$omp tip_count) | ||
| 32 | end module detect_oneline_mod | ||
| 33 | ! | ||
| 34 | module check_orbit_type_sub | ||
| 35 | |||
| 36 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 37 | use finite_check, only: is_finite | ||
| 38 | |||
| 39 | implicit none | ||
| 40 | |||
| 41 | contains | ||
| 42 | |||
| 43 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 44 | ! | ||
| 45 | ✗ | subroutine check_orbit_type(nturns,nfp,fpr_in,ideal,ijpar,ierr) | |
| 46 | ! | ||
| 47 | use detect_oneline_mod, only : prop,nfp_dim,nfp_max,iper, & | ||
| 48 | pi,twopi, & | ||
| 49 | fprs,igroup, & | ||
| 50 | ipermin, iret, & | ||
| 51 | jpar_variation_rate, rotation_number_drift, & | ||
| 52 | precession_turns, jpar_sample_count, & | ||
| 53 | rotation_half_count, & | ||
| 54 | jpar_spread, jpar_ref, topology_margin, & | ||
| 55 | score_status, tip_count | ||
| 56 | use sorting_mod, only : sortin | ||
| 57 | ! | ||
| 58 | implicit none | ||
| 59 | ! | ||
| 60 | logical :: did_recur | ||
| 61 | ! | ||
| 62 | integer :: nturns,ideal,ijpar | ||
| 63 | ! | ||
| 64 | integer :: nfp,iunit,i,j,k,ngroup,info,ierr,noline | ||
| 65 | integer :: ibeg,iend,npl,nturnm1 | ||
| 66 | ! | ||
| 67 | double precision :: dummy,dt1,dt2,derfix,per,ds1,ds2 | ||
| 68 | double precision :: theta_1,theta_2,theta_test,errint | ||
| 69 | double precision :: perpinv_beg, tol_perpinv, drift, margin | ||
| 70 | double precision, dimension(nfp_dim) :: fpr_in | ||
| 71 | ✗ | integer, dimension(:), allocatable :: ipoi | |
| 72 | ! | ||
| 73 | ✗ | tol_perpinv=15.d0 | |
| 74 | ! | ||
| 75 | ✗ | if(prop) then | |
| 76 | ✗ | prop=.false. | |
| 77 | ✗ | pi=atan(1.d0)*4.d0 | |
| 78 | ✗ | twopi=2.d0*pi | |
| 79 | ! | ||
| 80 | ! | ||
| 81 | ✗ | nfp_max=30000 | |
| 82 | ✗ | allocate(fprs(nfp_dim,nfp_max)) | |
| 83 | ✗ | allocate(iret(nfp_max)) | |
| 84 | ! | ||
| 85 | endif | ||
| 86 | ! | ||
| 87 | ✗ | ierr=0 | |
| 88 | ! | ||
| 89 | ✗ | if(nfp.eq.0) then | |
| 90 | ✗ | iper=0 | |
| 91 | ✗ | igroup=0 | |
| 92 | ✗ | ideal=0 | |
| 93 | ✗ | ijpar=0 | |
| 94 | ✗ | jpar_spread = 0.d0 | |
| 95 | ✗ | jpar_ref = 0.d0 | |
| 96 | ✗ | topology_margin = 0.d0 | |
| 97 | ✗ | jpar_variation_rate = 0.d0 | |
| 98 | ✗ | rotation_number_drift = 0.d0 | |
| 99 | ✗ | precession_turns = 0.d0 | |
| 100 | ✗ | jpar_sample_count = 0 | |
| 101 | ✗ | rotation_half_count = 0 | |
| 102 | ✗ | tip_count = 0 | |
| 103 | ✗ | score_status = 0 | |
| 104 | endif | ||
| 105 | ! | ||
| 106 | ✗ | if(igroup.ge.nturns) then | |
| 107 | ✗ | ierr=2 | |
| 108 | ✗ | return | |
| 109 | endif | ||
| 110 | ! | ||
| 111 | ✗ | nfp=nfp+1 | |
| 112 | ✗ | if(nfp.gt.nfp_max) then | |
| 113 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 114 | print *,'detect_oneline: nfp > nfp_max, nturns not reached' | ||
| 115 | #endif | ||
| 116 | ✗ | ierr=1 | |
| 117 | ✗ | return | |
| 118 | endif | ||
| 119 | ✗ | fprs(:,nfp)=fpr_in | |
| 120 | ! | ||
| 121 | ✗ | tip_count = nfp | |
| 122 | call compute_continuous_tip_scores(fprs, nfp, jpar_variation_rate, & | ||
| 123 | rotation_number_drift, precession_turns, jpar_sample_count, & | ||
| 124 | ✗ | rotation_half_count) | |
| 125 | ✗ | if (nfp .ge. 2) then | |
| 126 | ✗ | jpar_ref = abs(fprs(3, 2)) | |
| 127 | ✗ | jpar_spread = maxval(abs(fprs(3, 2:nfp) - fprs(3, 2))) | |
| 128 | end if | ||
| 129 | ! | ||
| 130 | ✗ | if(nfp.lt.3) return | |
| 131 | ! | ||
| 132 | ✗ | theta_1=fprs(2,1) | |
| 133 | ✗ | theta_2=modulo(fprs(2,2)-theta_1+pi,twopi)+theta_1-pi | |
| 134 | ✗ | theta_test=modulo(fprs(2,nfp)-theta_1+pi,twopi)+theta_1-pi | |
| 135 | ✗ | dt1=theta_test-theta_1 | |
| 136 | ✗ | dt2=theta_test-theta_2 | |
| 137 | ✗ | did_recur=dt1*dt2.lt.0.d0 | |
| 138 | ✗ | if(did_recur) then | |
| 139 | ✗ | if(iper.eq.0) then | |
| 140 | ✗ | iper=nfp | |
| 141 | ✗ | if(iper.le.ipermin) then | |
| 142 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 143 | print *,'detect_oneline: stochastic, iper < ipermin' | ||
| 144 | #endif | ||
| 145 | ✗ | ideal=2 | |
| 146 | ✗ | ijpar=2 | |
| 147 | ✗ | igroup=nturns | |
| 148 | ✗ | score_status = 2 | |
| 149 | ✗ | return | |
| 150 | endif | ||
| 151 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 152 | print *,'detect_oneline: return period',iper-1 | ||
| 153 | #endif | ||
| 154 | ✗ | igroup=1 | |
| 155 | else | ||
| 156 | ✗ | igroup=igroup+1 | |
| 157 | endif | ||
| 158 | ✗ | iret(igroup)=nfp-1 | |
| 159 | endif | ||
| 160 | ! | ||
| 161 | ✗ | if(igroup.eq.nturns) then | |
| 162 | !check all recurrencties: | ||
| 163 | ! | ||
| 164 | ✗ | call check_recurs(igroup,iret(1:igroup),noline) | |
| 165 | ! | ||
| 166 | ✗ | ideal=noline+1 | |
| 167 | ! | ||
| 168 | ! check sequencies in all intervals: | ||
| 169 | ✗ | margin = huge(1.d0) | |
| 170 | ✗ | if(ideal.eq.1) then | |
| 171 | ✗ | nturnm1=nturns-1 | |
| 172 | ✗ | allocate(ipoi(nturnm1)) | |
| 173 | ! | ||
| 174 | ✗ | call sortin(fprs(2,iret(1:nturnm1)+1),ipoi,nturnm1) | |
| 175 | ! | ||
| 176 | ✗ | do k=1,iret(nturns)-iret(nturnm1) | |
| 177 | ✗ | do j=2,nturnm1 | |
| 178 | ✗ | theta_1=fprs(2,iret(ipoi(j-1))+k+1) | |
| 179 | ✗ | theta_2=modulo(fprs(2,iret(ipoi(j))+k+1)-theta_1+pi,twopi)+theta_1-pi | |
| 180 | ✗ | margin = min(margin, theta_2 - theta_1) | |
| 181 | ✗ | if (theta_2 .lt. theta_1) ideal = 2 | |
| 182 | enddo | ||
| 183 | enddo | ||
| 184 | ✗ | deallocate(ipoi) | |
| 185 | endif | ||
| 186 | ! | ||
| 187 | ! Status 1 means topology_margin is a real margin. The recurrence test | ||
| 188 | ! decides without forming one, and the monotonicity loop is empty when | ||
| 189 | ! the last two return intervals coincide. | ||
| 190 | ✗ | if (margin .lt. huge(1.d0)) then | |
| 191 | ✗ | topology_margin = margin | |
| 192 | ✗ | score_status = 1 | |
| 193 | else | ||
| 194 | ✗ | score_status = 3 | |
| 195 | end if | ||
| 196 | ! | ||
| 197 | ! check J_perp change | ||
| 198 | ✗ | ijpar=1 | |
| 199 | ✗ | perpinv_beg=fprs(3,2) | |
| 200 | ✗ | jpar_ref = abs(perpinv_beg) | |
| 201 | ✗ | jpar_spread = 0.d0 | |
| 202 | ✗ | do k=2,iret(nturns) | |
| 203 | ✗ | drift = abs(fprs(3, k) - perpinv_beg) | |
| 204 | ✗ | jpar_spread = max(jpar_spread, drift) | |
| 205 | ✗ | if (drift .gt. tol_perpinv) ijpar = 2 | |
| 206 | enddo | ||
| 207 | endif | ||
| 208 | ! | ||
| 209 | end subroutine check_orbit_type | ||
| 210 | ! | ||
| 211 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 212 | ! | ||
| 213 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | pure subroutine compute_continuous_tip_scores(tips, n, jpar_rate, rotation_drift, & |
| 214 | achieved_turns, jpar_samples, rotation_samples) | ||
| 215 | real(dp), intent(in) :: tips(:, :) | ||
| 216 | integer, intent(in) :: n | ||
| 217 | real(dp), intent(out) :: jpar_rate, rotation_drift, achieved_turns | ||
| 218 | integer, intent(out) :: jpar_samples, rotation_samples | ||
| 219 | |||
| 220 | integer :: half, k, nincrements | ||
| 221 | real(dp) :: cumulative_turns, delta_theta, jpar_reference, pi_local | ||
| 222 | real(dp) :: rate, rotation_first, rotation_second, sum_rate_squared, twopi_local | ||
| 223 | |||
| 224 | 5 | pi_local = acos(-1.0_dp) | |
| 225 | 5 | twopi_local = 2.0_dp*pi_local | |
| 226 | 5 | jpar_rate = 0.0_dp | |
| 227 | 5 | rotation_drift = 0.0_dp | |
| 228 | 5 | achieved_turns = 0.0_dp | |
| 229 | 5 | jpar_samples = 0 | |
| 230 | 5 | rotation_samples = 0 | |
| 231 | |||
| 232 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 time.
|
5 | if (n .lt. 3) return |
| 233 |
8/8✓ Branch 0 taken 24 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 72 times.
✓ Branch 3 taken 23 times.
✓ Branch 4 taken 71 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 3 times.
✓ Branch 7 taken 1 time.
|
98 | if (.not. all(is_finite(tips(:, 1:n)))) return |
| 234 | |||
| 235 | 3 | jpar_reference = abs(tips(3, 2)) | |
| 236 | 3 | cumulative_turns = 0.0_dp | |
| 237 | 3 | sum_rate_squared = 0.0_dp | |
| 238 |
2/2✓ Branch 0 taken 14 times.
✓ Branch 1 taken 3 times.
|
17 | do k = 3, n |
| 239 |
3/4✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 12 times.
|
14 | delta_theta = modulo(tips(2, k) - tips(2, k - 1) + pi_local, twopi_local) - pi_local |
| 240 | 14 | cumulative_turns = cumulative_turns + abs(delta_theta)/twopi_local | |
| 241 |
2/2✓ Branch 0 taken 11 times.
✓ Branch 1 taken 3 times.
|
17 | if (jpar_reference .gt. tiny(1.0_dp)) then |
| 242 |
1/2✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
|
11 | if (cumulative_turns .gt. tiny(1.0_dp)) then |
| 243 | 11 | rate = (tips(3, k) - tips(3, 2))/jpar_reference/cumulative_turns | |
| 244 | 11 | sum_rate_squared = sum_rate_squared + rate*rate | |
| 245 | 11 | jpar_samples = jpar_samples + 1 | |
| 246 | end if | ||
| 247 | end if | ||
| 248 | end do | ||
| 249 | 3 | achieved_turns = cumulative_turns | |
| 250 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | if (jpar_samples .gt. 0) then |
| 251 | 2 | jpar_rate = sqrt(sum_rate_squared/real(jpar_samples, dp)) | |
| 252 | end if | ||
| 253 | |||
| 254 | 3 | nincrements = n - 1 | |
| 255 | 3 | half = nincrements/2 | |
| 256 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | if (half .lt. 2) return |
| 257 | |||
| 258 | 3 | rotation_first = 0.0_dp | |
| 259 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 3 times.
|
11 | do k = 2, half + 1 |
| 260 |
2/4✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
|
8 | delta_theta = modulo(tips(2, k) - tips(2, k - 1) + pi_local, twopi_local) - pi_local |
| 261 | 11 | rotation_first = rotation_first + delta_theta | |
| 262 | end do | ||
| 263 | 3 | rotation_second = 0.0_dp | |
| 264 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 3 times.
|
11 | do k = n - half + 1, n |
| 265 |
3/4✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
|
8 | delta_theta = modulo(tips(2, k) - tips(2, k - 1) + pi_local, twopi_local) - pi_local |
| 266 | 11 | rotation_second = rotation_second + delta_theta | |
| 267 | end do | ||
| 268 | 3 | rotation_samples = half | |
| 269 | 3 | rotation_first = rotation_first/(twopi_local*real(half, dp)) | |
| 270 | 3 | rotation_second = rotation_second/(twopi_local*real(half, dp)) | |
| 271 | 3 | rotation_drift = abs(rotation_first - rotation_second) | |
| 272 | end subroutine compute_continuous_tip_scores | ||
| 273 | ! | ||
| 274 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 275 | ! | ||
| 276 | ✗ | subroutine check_recurs(n,nrec,noline) | |
| 277 | ! | ||
| 278 | implicit none | ||
| 279 | ! | ||
| 280 | integer :: n,noline,k,kmax,mult,nfirst,nnext,maxmult | ||
| 281 | integer, dimension(n) :: nrec | ||
| 282 | ! | ||
| 283 | ✗ | noline=0 | |
| 284 | ! | ||
| 285 | ✗ | maxmult=n/2 | |
| 286 | ! | ||
| 287 | ✗ | do mult=1,maxmult | |
| 288 | ✗ | kmax=n/mult | |
| 289 | ✗ | if(kmax.lt.2) exit | |
| 290 | ✗ | nfirst=nrec(mult) | |
| 291 | ✗ | do k=2,kmax | |
| 292 | ✗ | nnext=nrec(mult*k)-nrec(mult*(k-1)) | |
| 293 | ✗ | if(nnext.gt.nfirst.or.nnext.lt.nfirst-1) then | |
| 294 | ✗ | noline=1 | |
| 295 | ✗ | return | |
| 296 | endif | ||
| 297 | enddo | ||
| 298 | enddo | ||
| 299 | ! | ||
| 300 | end subroutine check_recurs | ||
| 301 | ! | ||
| 302 | end module check_orbit_type_sub | ||
| 303 |