| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module orbit_frequencies | ||
| 2 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 3 | use simple, only: tracer_t, init_sympl, tstep | ||
| 4 | use orbit_symplectic_base, only: SYMPLECTIC_STEP_OUTSIDE_DOMAIN, & | ||
| 5 | SYMPLECTIC_STEP_MAXITER, & | ||
| 6 | SYMPLECTIC_STEP_BOUNDARY | ||
| 7 | use util, only: twopi | ||
| 8 | |||
| 9 | implicit none | ||
| 10 | private | ||
| 11 | |||
| 12 | integer, parameter, public :: ORBIT_CLASS_UNKNOWN = 0 | ||
| 13 | integer, parameter, public :: ORBIT_CLASS_TRAPPED = 1 | ||
| 14 | integer, parameter, public :: ORBIT_CLASS_PASSING = 2 | ||
| 15 | |||
| 16 | integer, parameter, public :: FREQ_SUCCESS = 0 | ||
| 17 | integer, parameter, public :: FREQ_INVALID_INPUT = 1 | ||
| 18 | integer, parameter, public :: FREQ_ORBIT_LOST = 2 | ||
| 19 | integer, parameter, public :: FREQ_INTEGRATOR_ERROR = 3 | ||
| 20 | integer, parameter, public :: FREQ_MAX_STEPS = 4 | ||
| 21 | |||
| 22 | type, public :: frequency_options_t | ||
| 23 | integer :: n_periods = 1 | ||
| 24 | integer :: max_steps = 10000000 | ||
| 25 | end type frequency_options_t | ||
| 26 | |||
| 27 | type, public :: frequency_result_t | ||
| 28 | integer :: status = FREQ_INVALID_INPUT | ||
| 29 | integer :: orbit_class = ORBIT_CLASS_UNKNOWN | ||
| 30 | integer :: parallel_direction = 0 | ||
| 31 | integer :: n_periods = 0 | ||
| 32 | integer :: n_steps = 0 | ||
| 33 | real(dp) :: period = 0.0_dp | ||
| 34 | real(dp) :: period_std = 0.0_dp | ||
| 35 | real(dp) :: delta_phi = 0.0_dp | ||
| 36 | real(dp) :: delta_phi_std = 0.0_dp | ||
| 37 | real(dp) :: omega_b = 0.0_dp | ||
| 38 | real(dp) :: omega_phi = 0.0_dp | ||
| 39 | end type frequency_result_t | ||
| 40 | |||
| 41 | public :: compute_canonical_frequencies | ||
| 42 | |||
| 43 | contains | ||
| 44 | |||
| 45 | 6 | subroutine compute_canonical_frequencies(tracer, initial_state, options, result) | |
| 46 | type(tracer_t), intent(in) :: tracer | ||
| 47 | real(dp), intent(in) :: initial_state(5) | ||
| 48 | type(frequency_options_t), intent(in) :: options | ||
| 49 | type(frequency_result_t), intent(out) :: result | ||
| 50 | |||
| 51 | type(tracer_t) :: work | ||
| 52 | real(dp) :: z(5), z_previous(5) | ||
| 53 | real(dp) :: event_time, event_phi, previous_event_time, previous_event_phi | ||
| 54 | real(dp) :: target_theta, fraction, step_time | ||
| 55 | 3 | real(dp), allocatable :: periods(:), displacements(:) | |
| 56 | integer :: ierr, step, n_events, direction | ||
| 57 | logical :: have_previous_event, event_found | ||
| 58 | |||
| 59 | result = frequency_result_t() | ||
| 60 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
4 | if (options%n_periods < 1 .or. options%max_steps < 1) return |
| 61 |
2/4✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
2 | if (tracer%dtaumin <= 0.0_dp .or. tracer%v0 <= 0.0_dp) return |
| 62 | |||
| 63 |
9/18✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
|
8 | allocate (periods(options%n_periods), displacements(options%n_periods)) |
| 64 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
|
8 | periods = 0.0_dp |
| 65 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
|
8 | displacements = 0.0_dp |
| 66 | 2 | work = tracer | |
| 67 | 2 | z = initial_state | |
| 68 | 2 | step_time = work%dtaumin/work%v0 | |
| 69 | |||
| 70 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (work%integmode > 0) then |
| 71 | call init_sympl(work%si, work%f, z, work%dtaumin, work%dtaumin, & | ||
| 72 | 2 | work%relerr, work%integmode) | |
| 73 | else | ||
| 74 | ✗ | work%dtau = work%dtaumin | |
| 75 | end if | ||
| 76 | |||
| 77 | have_previous_event = .false. | ||
| 78 | n_events = 0 | ||
| 79 | direction = 0 | ||
| 80 | target_theta = 0.0_dp | ||
| 81 | |||
| 82 |
1/2✓ Branch 0 taken 68996 times.
✗ Branch 1 not taken.
|
68996 | do step = 1, options%max_steps |
| 83 | 68996 | z_previous = z | |
| 84 |
1/2✓ Branch 0 taken 68996 times.
✗ Branch 1 not taken.
|
68996 | if (work%integmode > 0) then |
| 85 | 68996 | call tstep(work%si, work%f, z, ierr) | |
| 86 | else | ||
| 87 | ✗ | call tstep(work, z, ierr) | |
| 88 | end if | ||
| 89 | 68996 | result%n_steps = step | |
| 90 | |||
| 91 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 68996 times.
|
68996 | if (ierr /= 0) then |
| 92 | ! A step that stops on the outer radial boundary is a physical | ||
| 93 | ! loss through the plasma edge, exactly like leaving the domain; | ||
| 94 | ! only genuinely numerical failures report an integrator error. | ||
| 95 | if (ierr == SYMPLECTIC_STEP_OUTSIDE_DOMAIN & | ||
| 96 | .or. ierr == SYMPLECTIC_STEP_MAXITER & | ||
| 97 | ✗ | .or. ierr == SYMPLECTIC_STEP_BOUNDARY) then | |
| 98 | ✗ | result%status = FREQ_ORBIT_LOST | |
| 99 | else | ||
| 100 | ✗ | result%status = FREQ_INTEGRATOR_ERROR | |
| 101 | end if | ||
| 102 | ✗ | return | |
| 103 | end if | ||
| 104 | |||
| 105 | 68996 | event_found = .false. | |
| 106 | 68996 | fraction = 0.0_dp | |
| 107 | |||
| 108 |
2/2✓ Branch 0 taken 12638 times.
✓ Branch 1 taken 56358 times.
|
68996 | if (result%orbit_class == ORBIT_CLASS_UNKNOWN) then |
| 109 |
4/4✓ Branch 0 taken 8018 times.
✓ Branch 1 taken 4620 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 8017 times.
|
12638 | if (z_previous(5) < 0.0_dp .and. z(5) >= 0.0_dp) then |
| 110 | 1 | result%orbit_class = ORBIT_CLASS_TRAPPED | |
| 111 | 1 | result%parallel_direction = 0 | |
| 112 | 1 | fraction = zero_crossing_fraction(z_previous(5), z(5)) | |
| 113 | 1 | event_found = .true. | |
| 114 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 12636 times.
|
12637 | else if (abs(z(2) - initial_state(2)) >= twopi) then |
| 115 | 1 | result%orbit_class = ORBIT_CLASS_PASSING | |
| 116 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | direction = merge(1, -1, z(2) > initial_state(2)) |
| 117 | 1 | result%parallel_direction = direction | |
| 118 | 1 | target_theta = initial_state(2) + real(direction, dp)*twopi | |
| 119 | fraction = zero_crossing_fraction(z_previous(2) - target_theta, & | ||
| 120 | 2 | z(2) - target_theta) | |
| 121 | 1 | event_found = .true. | |
| 122 | 1 | target_theta = target_theta + real(direction, dp)*twopi | |
| 123 | end if | ||
| 124 |
2/2✓ Branch 0 taken 47689 times.
✓ Branch 1 taken 8669 times.
|
56358 | else if (result%orbit_class == ORBIT_CLASS_TRAPPED) then |
| 125 |
4/4✓ Branch 0 taken 24022 times.
✓ Branch 1 taken 23667 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 24019 times.
|
47689 | if (z_previous(5) < 0.0_dp .and. z(5) >= 0.0_dp) then |
| 126 | 3 | fraction = zero_crossing_fraction(z_previous(5), z(5)) | |
| 127 | 3 | event_found = .true. | |
| 128 | end if | ||
| 129 | else | ||
| 130 |
1/2✓ Branch 0 taken 8669 times.
✗ Branch 1 not taken.
|
8669 | if (direction > 0) then |
| 131 |
3/4✓ Branch 0 taken 8669 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 8666 times.
|
8669 | if (z_previous(2) < target_theta .and. z(2) >= target_theta) then |
| 132 | fraction = zero_crossing_fraction( & | ||
| 133 | z_previous(2) - target_theta, z(2) - target_theta & | ||
| 134 | 6 | ) | |
| 135 | event_found = .true. | ||
| 136 | end if | ||
| 137 | else | ||
| 138 | ✗ | if (z_previous(2) > target_theta .and. z(2) <= target_theta) then | |
| 139 | fraction = zero_crossing_fraction( & | ||
| 140 | z_previous(2) - target_theta, z(2) - target_theta & | ||
| 141 | ✗ | ) | |
| 142 | event_found = .true. | ||
| 143 | end if | ||
| 144 | end if | ||
| 145 | 3 | if (event_found) target_theta = target_theta + real(direction, dp)*twopi | |
| 146 | end if | ||
| 147 | |||
| 148 |
2/2✓ Branch 0 taken 68988 times.
✓ Branch 1 taken 8 times.
|
68996 | if (.not. event_found) cycle |
| 149 | |||
| 150 | 8 | event_time = (real(step - 1, dp) + fraction)*step_time | |
| 151 | 8 | event_phi = z_previous(3) + fraction*(z(3) - z_previous(3)) | |
| 152 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
|
8 | if (have_previous_event) then |
| 153 | 6 | n_events = n_events + 1 | |
| 154 | 6 | periods(n_events) = event_time - previous_event_time | |
| 155 | 6 | displacements(n_events) = event_phi - previous_event_phi | |
| 156 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 4 times.
|
6 | if (n_events == options%n_periods) exit |
| 157 | end if | ||
| 158 | previous_event_time = event_time | ||
| 159 | previous_event_phi = event_phi | ||
| 160 | 68990 | have_previous_event = .true. | |
| 161 | end do | ||
| 162 | |||
| 163 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (n_events /= options%n_periods) then |
| 164 | ✗ | result%status = FREQ_MAX_STEPS | |
| 165 | ✗ | result%n_periods = n_events | |
| 166 | ✗ | return | |
| 167 | end if | ||
| 168 | |||
| 169 | 2 | result%n_periods = n_events | |
| 170 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
|
8 | result%period = sum(periods)/real(n_events, dp) |
| 171 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
|
8 | result%delta_phi = sum(displacements)/real(n_events, dp) |
| 172 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (n_events > 1) then |
| 173 | 2 | result%period_std = sample_std(periods, result%period) | |
| 174 | 2 | result%delta_phi_std = sample_std(displacements, result%delta_phi) | |
| 175 | end if | ||
| 176 | 2 | result%omega_b = twopi/result%period | |
| 177 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✓ Branch 3 taken 2 times.
|
16 | result%omega_phi = sum(displacements)/sum(periods) |
| 178 | 2 | result%status = FREQ_SUCCESS | |
| 179 |
18/22✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 18 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 192 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 9 times.
✓ Branch 7 taken 3 times.
✓ Branch 8 taken 9 times.
✓ Branch 9 taken 3 times.
✓ Branch 10 taken 9 times.
✓ Branch 11 taken 3 times.
✓ Branch 12 taken 2 times.
✓ Branch 13 taken 1 time.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 1 time.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 time.
|
252 | end subroutine compute_canonical_frequencies |
| 180 | |||
| 181 | 8 | pure real(dp) function zero_crossing_fraction(left, right) result(fraction) | |
| 182 | real(dp), intent(in) :: left, right | ||
| 183 | real(dp) :: denominator | ||
| 184 | |||
| 185 | 8 | denominator = right - left | |
| 186 |
4/10✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
8 | if (abs(denominator) <= tiny(denominator)) then |
| 187 | fraction = 0.5_dp | ||
| 188 | else | ||
| 189 | 8 | fraction = -left/denominator | |
| 190 | end if | ||
| 191 | 8 | fraction = max(0.0_dp, min(1.0_dp, fraction)) | |
| 192 | end function zero_crossing_fraction | ||
| 193 | |||
| 194 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | pure real(dp) function sample_std(values, mean_value) result(std_value) |
| 195 | real(dp), intent(in) :: values(:), mean_value | ||
| 196 | |||
| 197 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 4 times.
|
16 | std_value = sqrt(sum((values - mean_value)**2)/real(size(values) - 1, dp)) |
| 198 | 4 | end function sample_std | |
| 199 | |||
| 200 | ✗ | end module orbit_frequencies | |
| 201 |