| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module spectre_orbit | ||
| 2 | !> Single-volume SPECTRE RK45 stepping, with an absorbing outer boundary | ||
| 3 | !> and reflection before the singular coordinate axis. | ||
| 4 | |||
| 5 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 6 | use finite_check, only: is_finite | ||
| 7 | use odeint_allroutines_sub, only: odeint_allroutines, odeint_has_failed | ||
| 8 | use alpha_lifetime_sub, only: velo_can | ||
| 9 | |||
| 10 | implicit none | ||
| 11 | private | ||
| 12 | |||
| 13 | public :: spectre_orbit_state_t, spectre_event_t, spectre_state_reset, & | ||
| 14 | orbit_timestep_spectre, SPECTRE_OK, SPECTRE_BOUNDARY, SPECTRE_FAULT | ||
| 15 | |||
| 16 | integer, parameter :: SPECTRE_OK = 0 | ||
| 17 | integer, parameter :: SPECTRE_BOUNDARY = 88 | ||
| 18 | integer, parameter :: SPECTRE_FAULT = 2 | ||
| 19 | |||
| 20 | real(dp), parameter :: axis_offset = 1.0d-3 | ||
| 21 | real(dp), parameter :: rho_tol = 1.0d-10 | ||
| 22 | integer, parameter :: max_bisect = 80 | ||
| 23 | integer, parameter :: max_retry_depth = 8 | ||
| 24 | |||
| 25 | !> Boundary limits of the single-volume clamped RHS, private to each thread. | ||
| 26 | real(dp) :: ev_lo = 0.0_dp, ev_hi = 0.0_dp | ||
| 27 | !$omp threadprivate(ev_lo, ev_hi) | ||
| 28 | |||
| 29 | !> Keep the first accepted state strictly inside the outer boundary. | ||
| 30 | real(dp), parameter :: clamp_margin = 1.0d-9 | ||
| 31 | !> Preserve the existing single-volume boundary regularization. | ||
| 32 | real(dp), parameter :: drift_band = 2.0d-2 | ||
| 33 | !> Backstop for the same Bstar_par -> 0 breakdown when it falls inside the | ||
| 34 | !> band (it is a pitch-dependent phase-space surface, not a fixed radius): | ||
| 35 | !> cap the RHS magnitude so odeint keeps a finite step through the unphysical | ||
| 36 | !> layer instead of collapsing to zero step and exhausting its step budget. | ||
| 37 | real(dp), parameter :: drift_cap = 1.0d0 | ||
| 38 | |||
| 39 | type :: spectre_orbit_state_t | ||
| 40 | logical :: initialized = .false. | ||
| 41 | end type spectre_orbit_state_t | ||
| 42 | |||
| 43 | type :: spectre_event_t | ||
| 44 | real(dp) :: t_frac = 0.0_dp | ||
| 45 | end type spectre_event_t | ||
| 46 | |||
| 47 | contains | ||
| 48 | |||
| 49 | 68 | subroutine spectre_state_reset(state, mvol) | |
| 50 | type(spectre_orbit_state_t), intent(out) :: state | ||
| 51 | integer, intent(in) :: mvol | ||
| 52 | |||
| 53 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
|
68 | if (mvol /= 1) error stop 'SPECTRE orbit stepping requires one volume' |
| 54 | state%initialized = .false. | ||
| 55 | 68 | end subroutine spectre_state_reset | |
| 56 | |||
| 57 | 397982 | subroutine orbit_timestep_spectre(state, z, dtaumin, relerr, ierr, event) | |
| 58 | !> Advance z = (rho_g, theta, zeta, p, lambda) by one microstep. | ||
| 59 | type(spectre_orbit_state_t), intent(inout) :: state | ||
| 60 | real(dp), intent(inout) :: z(5) | ||
| 61 | real(dp), intent(in) :: dtaumin, relerr | ||
| 62 | integer, intent(out) :: ierr | ||
| 63 | type(spectre_event_t), intent(out) :: event | ||
| 64 | |||
| 65 | real(dp) :: z_start(5), z_end(5), z_hit(5), boundary | ||
| 66 | integer :: direction | ||
| 67 | |||
| 68 | 397982 | ierr = SPECTRE_OK | |
| 69 | 397982 | z_start = z | |
| 70 | |||
| 71 |
2/2✓ Branch 0 taken 68 times.
✓ Branch 1 taken 397914 times.
|
397982 | if (.not. state%initialized) then |
| 72 | 68 | ev_lo = axis_offset | |
| 73 | 68 | ev_hi = 1.0_dp | |
| 74 | 68 | call integrate_clamped(z_start, dtaumin, relerr, z_end, ierr) | |
| 75 |
1/2✓ Branch 0 taken 68 times.
✗ Branch 1 not taken.
|
397930 | if (ierr /= SPECTRE_OK) return |
| 76 | 68 | state%initialized = .true. | |
| 77 | 68 | z = z_end | |
| 78 | 68 | z(1) = max(0.0_dp, min(z(1), 1.0_dp - clamp_margin)) | |
| 79 | 68 | return | |
| 80 | end if | ||
| 81 | |||
| 82 | 397914 | ev_lo = 0.0_dp | |
| 83 | 397914 | ev_hi = 1.0_dp | |
| 84 | 397914 | call integrate_clamped(z_start, dtaumin, relerr, z_end, ierr) | |
| 85 |
1/2✓ Branch 0 taken 397914 times.
✗ Branch 1 not taken.
|
397914 | if (ierr /= SPECTRE_OK) return |
| 86 | |||
| 87 |
3/6✓ Branch 0 taken 397850 times.
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 397850 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
397914 | if (z_end(1) > 1.0_dp .or. & |
| 88 | (z_end(1) == 1.0_dp .and. z_start(1) < 1.0_dp)) then | ||
| 89 | 64 | boundary = 1.0_dp | |
| 90 | 64 | direction = 1 | |
| 91 |
3/6✓ Branch 0 taken 397794 times.
✓ Branch 1 taken 56 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 397794 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
397850 | else if (z_end(1) < 0.0_dp .or. & |
| 92 | (z_end(1) == 0.0_dp .and. z_start(1) > 0.0_dp)) then | ||
| 93 | 56 | boundary = 0.0_dp | |
| 94 | 56 | direction = -1 | |
| 95 | else | ||
| 96 | 397794 | z = z_end | |
| 97 | 397794 | return | |
| 98 | end if | ||
| 99 | |||
| 100 | call locate_boundary(z_start, dtaumin, relerr, boundary, z_hit, & | ||
| 101 | 120 | event%t_frac, ierr) | |
| 102 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 120 times.
|
120 | if (ierr /= SPECTRE_OK) then |
| 103 | ✗ | z = z_start | |
| 104 | ✗ | return | |
| 105 | end if | ||
| 106 | 120 | z = z_hit | |
| 107 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 56 times.
|
120 | if (direction > 0) then |
| 108 | 64 | z(1) = 1.0_dp | |
| 109 | 64 | ierr = SPECTRE_BOUNDARY | |
| 110 | else | ||
| 111 | 56 | z(1) = axis_offset | |
| 112 | 56 | z(5) = -z_hit(5) | |
| 113 | end if | ||
| 114 | end subroutine orbit_timestep_spectre | ||
| 115 | |||
| 116 | 120 | subroutine locate_boundary(z_start, dtaumin, relerr, boundary, z_hit, t_frac, & | |
| 117 | ierr) | ||
| 118 | !> Find the microstep time where the clamped-RHS trajectory reaches rho_g = | ||
| 119 | !> boundary to |rho_g - boundary| < rho_tol, by the Illinois variant of | ||
| 120 | !> false position. z_start is strictly inside the home volume and the full | ||
| 121 | !> step already overshot the boundary, so [0, dtaumin] brackets the root. | ||
| 122 | real(dp), intent(in) :: z_start(5), dtaumin, relerr, boundary | ||
| 123 | real(dp), intent(out) :: z_hit(5) | ||
| 124 | real(dp), intent(out) :: t_frac | ||
| 125 | integer, intent(out) :: ierr | ||
| 126 | |||
| 127 | real(dp) :: t_lo, t_hi, f_lo, f_hi, t_mid, f_mid, z_mid(5) | ||
| 128 | integer :: it, last_side | ||
| 129 | logical :: converged | ||
| 130 | |||
| 131 | 120 | t_lo = 0.0_dp | |
| 132 | 120 | ierr = SPECTRE_OK | |
| 133 | 120 | f_lo = z_start(1) - boundary | |
| 134 | 120 | t_hi = dtaumin | |
| 135 | 120 | call integrate_clamped(z_start, t_hi, relerr, z_mid, ierr) | |
| 136 |
1/2✓ Branch 0 taken 120 times.
✗ Branch 1 not taken.
|
120 | if (ierr /= SPECTRE_OK) return |
| 137 | 120 | f_hi = z_mid(1) - boundary | |
| 138 | 120 | z_hit = z_mid | |
| 139 | 120 | t_mid = t_hi | |
| 140 | 120 | last_side = 0 | |
| 141 | 120 | converged = abs(f_hi) < rho_tol | |
| 142 | |||
| 143 |
1/2✓ Branch 0 taken 506 times.
✗ Branch 1 not taken.
|
506 | do it = 1, max_bisect |
| 144 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 506 times.
|
506 | if (converged) exit |
| 145 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 506 times.
|
506 | if (abs(f_hi - f_lo) <= tiny(1.0_dp)) exit |
| 146 | 506 | t_mid = t_hi - f_hi*(t_hi - t_lo)/(f_hi - f_lo) | |
| 147 | 506 | call integrate_clamped(z_start, t_mid, relerr, z_mid, ierr) | |
| 148 |
1/2✓ Branch 0 taken 506 times.
✗ Branch 1 not taken.
|
506 | if (ierr /= SPECTRE_OK) return |
| 149 | 506 | f_mid = z_mid(1) - boundary | |
| 150 | 506 | z_hit = z_mid | |
| 151 |
2/2✓ Branch 0 taken 386 times.
✓ Branch 1 taken 120 times.
|
506 | if (abs(f_mid) < rho_tol) then |
| 152 | converged = .true. | ||
| 153 | exit | ||
| 154 | end if | ||
| 155 |
2/2✓ Branch 0 taken 168 times.
✓ Branch 1 taken 218 times.
|
386 | if (f_mid*f_lo > 0.0_dp) then |
| 156 | 168 | t_lo = t_mid | |
| 157 | 168 | f_lo = f_mid | |
| 158 |
2/2✓ Branch 0 taken 44 times.
✓ Branch 1 taken 124 times.
|
168 | if (last_side == 1) f_hi = 0.5_dp*f_hi |
| 159 | last_side = 1 | ||
| 160 | else | ||
| 161 | 218 | t_hi = t_mid | |
| 162 | 218 | f_hi = f_mid | |
| 163 |
2/2✓ Branch 0 taken 78 times.
✓ Branch 1 taken 140 times.
|
218 | if (last_side == -1) f_lo = 0.5_dp*f_lo |
| 164 | last_side = -1 | ||
| 165 | end if | ||
| 166 | end do | ||
| 167 | |||
| 168 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 120 times.
|
120 | if (.not. converged) then |
| 169 | ✗ | z_hit = z_start | |
| 170 | ✗ | t_frac = 0.0_dp | |
| 171 | ✗ | ierr = SPECTRE_FAULT | |
| 172 | ✗ | return | |
| 173 | end if | ||
| 174 | 120 | t_frac = t_mid/dtaumin | |
| 175 | end subroutine locate_boundary | ||
| 176 | |||
| 177 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2412638 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2412638 times.
|
2412638 | subroutine velo_can_clamped(tau, z, vz) |
| 178 | !> velo_can with rho_g clamped to the home volume [ev_lo, ev_hi). Inside | ||
| 179 | !> the volume the drift is exact; past a boundary the field is held at | ||
| 180 | !> its boundary value. The final state can move beyond the boundary, | ||
| 181 | !> allowing the caller to locate the loss or axis-reflection event. | ||
| 182 | real(dp), intent(in) :: tau | ||
| 183 | real(dp), intent(in) :: z(:) | ||
| 184 | real(dp), intent(out) :: vz(:) | ||
| 185 | |||
| 186 | 2412638 | real(dp) :: zc(size(z)), vmag | |
| 187 | |||
| 188 |
2/2✓ Branch 0 taken 12063190 times.
✓ Branch 1 taken 2412638 times.
|
14475828 | zc = z |
| 189 | zc(1) = max(max(ev_lo + drift_band, axis_offset), & | ||
| 190 | 2412638 | min(z(1), ev_hi - drift_band)) | |
| 191 | 2412638 | call velo_can(tau, zc, vz) | |
| 192 |
2/2✓ Branch 0 taken 12063190 times.
✓ Branch 1 taken 2412638 times.
|
14475828 | vmag = sqrt(sum(vz**2)) |
| 193 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2412638 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
2412638 | if (vmag > drift_cap) vz = vz*(drift_cap/vmag) |
| 194 | 2412638 | end subroutine velo_can_clamped | |
| 195 | |||
| 196 | 398608 | subroutine integrate_clamped(z_in, tau, relerr, z_out, ierr) | |
| 197 | real(dp), intent(in) :: z_in(5), tau, relerr | ||
| 198 | real(dp), intent(out) :: z_out(5) | ||
| 199 | integer, intent(out) :: ierr | ||
| 200 | |||
| 201 | 398608 | call integrate_clamped_retry(z_in, tau, relerr, 0, z_out, ierr) | |
| 202 | end subroutine integrate_clamped | ||
| 203 | |||
| 204 | 398608 | recursive subroutine integrate_clamped_retry(z_in, tau, relerr, depth, & | |
| 205 | z_out, ierr) | ||
| 206 | real(dp), intent(in) :: z_in(5), tau, relerr | ||
| 207 | integer, intent(in) :: depth | ||
| 208 | real(dp), intent(out) :: z_out(5) | ||
| 209 | integer, intent(out) :: ierr | ||
| 210 | |||
| 211 | real(dp) :: z_half(5) | ||
| 212 | integer, parameter :: ndim = 5 | ||
| 213 | |||
| 214 | 398608 | ierr = SPECTRE_OK | |
| 215 | 398608 | z_out = z_in | |
| 216 |
1/2✓ Branch 0 taken 398608 times.
✗ Branch 1 not taken.
|
797216 | if (tau <= 0.0_dp) return |
| 217 | 398608 | call odeint_allroutines(z_out, ndim, 0.0_dp, tau, relerr, velo_can_clamped) | |
| 218 |
5/8✓ Branch 0 taken 1993040 times.
✓ Branch 1 taken 398608 times.
✓ Branch 2 taken 1993040 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 398608 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 398608 times.
|
2391648 | if (.not. odeint_has_failed() .and. all(is_finite(z_out))) return |
| 219 | |||
| 220 | ✗ | z_out = z_in | |
| 221 | ✗ | ierr = SPECTRE_FAULT | |
| 222 | ✗ | if (depth >= max_retry_depth) return | |
| 223 | call integrate_clamped_retry(z_in, 0.5_dp*tau, relerr, depth + 1, & | ||
| 224 | ✗ | z_half, ierr) | |
| 225 | ✗ | if (ierr /= SPECTRE_OK) return | |
| 226 | call integrate_clamped_retry(z_half, 0.5_dp*tau, relerr, depth + 1, & | ||
| 227 | ✗ | z_out, ierr) | |
| 228 | ✗ | if (ierr /= SPECTRE_OK) z_out = z_in | |
| 229 | end subroutine integrate_clamped_retry | ||
| 230 | |||
| 231 | ✗ | end module spectre_orbit | |
| 232 |