| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module simple | ||
| 2 | use util, only: c, e_charge, p_mass, ev, twopi | ||
| 3 | use new_vmec_stuff_mod, only : netcdffile, multharm, ns_s, ns_tp, & | ||
| 4 | vmec_B_scale, vmec_RZ_scale | ||
| 5 | |||
| 6 | use parmot_mod, only : rmu, ro0 | ||
| 7 | use velo_mod, only : isw_field_type | ||
| 8 | use field_can_mod, only : field_can_t | ||
| 9 | use orbit_symplectic, only : symplectic_integrator_t, multistage_integrator_t, & | ||
| 10 | orbit_sympl_init, orbit_timestep_sympl | ||
| 11 | use field, only : vmec_field_t | ||
| 12 | use field_can_mod, only : eval_field => evaluate, init_field_can, field_can_t, get_val | ||
| 13 | use orbit_fo_boris, only : fo_state_t, fo_init, fo_step, fo_to_gc | ||
| 14 | use diag_mod, only : icounter | ||
| 15 | use chamb_sub, only : chamb_can | ||
| 16 | |||
| 17 | implicit none | ||
| 18 | |||
| 19 | ! Define real(dp) kind parameter | ||
| 20 | integer, parameter :: dp = kind(1.0d0) | ||
| 21 | save | ||
| 22 | |||
| 23 | public | ||
| 24 | |||
| 25 | integer, parameter :: ORBIT_FO_LOSS = 2 | ||
| 26 | integer, parameter :: ORBIT_FO_NUMERICAL = 3 | ||
| 27 | |||
| 28 | type :: tracer_t | ||
| 29 | real(dp) :: fper | ||
| 30 | real(dp) :: dtau, dtaumin, v0 | ||
| 31 | integer :: n_e, n_d | ||
| 32 | |||
| 33 | ! integmode: 0 = RK, 1 = Euler1, 2 = Euler2, 3 = Midpoint, 4-7 = Gauss1-4, | ||
| 34 | ! 15 = Lobatto3, 20 = Gauss-Radau, 21 = Bulirsch-Stoer, 22 = Cash-Karp, | ||
| 35 | ! 24 = TDRK (fixed step), 25 = TDRK (adaptive) | ||
| 36 | integer :: integmode = 0 | ||
| 37 | real(dp) :: relerr | ||
| 38 | |||
| 39 | type(field_can_t) :: f | ||
| 40 | type(symplectic_integrator_t) :: si | ||
| 41 | type(multistage_integrator_t) :: mi | ||
| 42 | type(fo_state_t) :: fo ! full-orbit Boris state (orbit_model=ORBIT_FULL_ORBIT) | ||
| 43 | end type tracer_t | ||
| 44 | |||
| 45 | interface tstep | ||
| 46 | module procedure timestep | ||
| 47 | module procedure timestep_z | ||
| 48 | module procedure timestep_sympl_z | ||
| 49 | end interface tstep | ||
| 50 | |||
| 51 | public :: encode_symplectic_state | ||
| 52 | |||
| 53 | contains | ||
| 54 | |||
| 55 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 641 times.
|
641 | subroutine encode_symplectic_state(f, z0, zsympl, h0, jperp, psi_star) |
| 56 | !> Convert the public GC state to the normalization used by the symplectic | ||
| 57 | !> Hamiltonian. This is the single source of truth used by initialization, | ||
| 58 | !> reseeding, and the invariant API. | ||
| 59 | type(field_can_t), intent(inout) :: f | ||
| 60 | real(dp), intent(in) :: z0(:) | ||
| 61 | real(dp), intent(out) :: zsympl(4) | ||
| 62 | real(dp), intent(out), optional :: h0, jperp, psi_star | ||
| 63 | |||
| 64 | real(dp) :: jperp_local, psi_star_local | ||
| 65 | |||
| 66 | 641 | call eval_field(f, z0(1), z0(2), z0(3), 0) | |
| 67 | 641 | jperp_local = z0(4)**2*(1.0_dp - z0(5)**2)/f%Bmod | |
| 68 | 641 | psi_star_local = f%Aph + ro0*z0(4)*z0(5)*f%hph | |
| 69 | |||
| 70 | 641 | f%mu = jperp_local | |
| 71 | 641 | f%ro0 = ro0/sqrt(2.0_dp) | |
| 72 | 641 | f%vpar = z0(4)*z0(5)*sqrt(2.0_dp) | |
| 73 |
2/2✓ Branch 0 taken 1923 times.
✓ Branch 1 taken 641 times.
|
2564 | zsympl(1:3) = z0(1:3) |
| 74 | 641 | zsympl(4) = f%vpar*f%hph + f%Aph/f%ro0 | |
| 75 | |||
| 76 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 638 times.
|
641 | if (present(h0)) h0 = z0(4)**2 |
| 77 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 638 times.
|
641 | if (present(jperp)) jperp = jperp_local |
| 78 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 638 times.
|
641 | if (present(psi_star)) psi_star = psi_star_local |
| 79 | 641 | end subroutine encode_symplectic_state | |
| 80 | |||
| 81 | |||
| 82 | 73 | subroutine init_vmec(vmec_file, ans_s, ans_tp, amultharm, fper) | |
| 83 | use spline_vmec_sub, only : spline_vmec_data, volume_and_B00 | ||
| 84 | use vmecin_sub, only : stevvo | ||
| 85 | |||
| 86 | character(*), intent(in) :: vmec_file | ||
| 87 | integer, intent(in) :: ans_s, ans_tp, amultharm | ||
| 88 | real(dp), intent(out) :: fper | ||
| 89 | |||
| 90 | integer :: L1i | ||
| 91 | real(dp) :: RT0, R0i, cbfi, bz0i, bf0, volume, B00 | ||
| 92 | |||
| 93 | ! TODO: Remove side effects | ||
| 94 | 73 | netcdffile = vmec_file | |
| 95 | 73 | ns_s = ans_s | |
| 96 | 73 | ns_tp = ans_tp | |
| 97 | 73 | multharm = amultharm | |
| 98 | |||
| 99 | 73 | call spline_vmec_data ! initialize splines for VMEC field | |
| 100 | 73 | call stevvo(RT0, R0i, L1i, cbfi, bz0i, bf0) ! initialize periods and major radius | |
| 101 | 73 | fper = twopi/dble(L1i) !<= field period | |
| 102 | 73 | print *, 'R0 = ', RT0, ' cm, fper = ', fper | |
| 103 | 73 | call volume_and_B00(volume,B00) | |
| 104 | 73 | print *,'volume = ',volume,' cm^3, B_00 = ',B00,' G' | |
| 105 |
1/2✓ Branch 0 taken 73 times.
✗ Branch 1 not taken.
|
73 | end subroutine init_vmec |
| 106 | |||
| 107 | |||
| 108 | 8 | subroutine init_params(self, Z_charge, m_mass, E_kin, npoints, store_step, relerr) | |
| 109 | ! Initializes normalization for velocity and Larmor radius based on kinetic energy | ||
| 110 | ! of plasma particles (= temperature for thermal particles). | ||
| 111 | use new_vmec_stuff_mod, only : rmajor | ||
| 112 | |||
| 113 | type(tracer_t) :: self | ||
| 114 | integer, intent(in), optional :: Z_charge, m_mass | ||
| 115 | real(dp), intent(in), optional :: E_kin | ||
| 116 | integer, intent(in), optional :: npoints ! Integrator resolution. Number of | ||
| 117 | ! integration points for strongly passing particles around the torus | ||
| 118 | |||
| 119 | integer, intent(in), optional :: store_step ! Store every X timesteps | ||
| 120 | real(dp), intent(in), optional :: relerr ! Relative error | ||
| 121 | |||
| 122 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
8 | if (present(Z_charge)) self%n_e = Z_charge |
| 123 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
8 | if (present(m_mass)) self%n_d = m_mass |
| 124 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
8 | if (present(relerr)) self%relerr = relerr |
| 125 | |||
| 126 | ! Neglect relativistic effects by large inverse relativistic temperature | ||
| 127 | 8 | rmu=1d8 | |
| 128 | |||
| 129 | ! Reference velocity and normalized Larmor radius | ||
| 130 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
8 | if (present(E_kin)) then |
| 131 | 8 | self%v0 = sqrt(2.d0*E_kin*ev/(self%n_d*p_mass)) | |
| 132 | else | ||
| 133 | ✗ | self%v0 = sqrt(2.d0*3.5d6*ev/(self%n_d*p_mass)) | |
| 134 | end if | ||
| 135 | |||
| 136 | ! Larmor radius times magnetic field: | ||
| 137 | 8 | ro0=self%v0*self%n_d*p_mass*c/(self%n_e*e_charge) | |
| 138 | |||
| 139 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
8 | if (present(npoints)) then |
| 140 | 8 | self%dtaumin=twopi*rmajor*1.0d2/npoints | |
| 141 | else | ||
| 142 | ✗ | self%dtaumin=twopi*rmajor*1.0d2/256.0d0 | |
| 143 | end if | ||
| 144 | |||
| 145 | ! timestep where to get results | ||
| 146 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
8 | if (present(store_step)) then |
| 147 | 8 | self%dtau = store_step*self%dtaumin | |
| 148 | else | ||
| 149 | ✗ | self%dtau = self%dtaumin | |
| 150 | end if | ||
| 151 | |||
| 152 | 8 | end subroutine init_params | |
| 153 | |||
| 154 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 636 times.
|
636 | subroutine init_sympl(si, f, z0, dtau, dtaumin, rtol_init, mode_init) |
| 155 | ! | ||
| 156 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 157 | type(field_can_t), intent(inout) :: f | ||
| 158 | real(dp), intent(in) :: z0(:) | ||
| 159 | real(dp), intent(in) :: dtau, dtaumin | ||
| 160 | real(dp), intent(in) :: rtol_init | ||
| 161 | integer, intent(in) :: mode_init ! 1 = expl.-impl. Euler, 2 = impl.-expl. Euler, 3 = Midpoint | ||
| 162 | |||
| 163 | real(dp) :: z(4) | ||
| 164 | |||
| 165 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 636 times.
|
636 | if(min(dabs(mod(dtau, dtaumin)), dabs(mod(dtau, dtaumin)-dtaumin)) > 1d-9*dtaumin) then |
| 166 | ✗ | print *, 'dtau = ', dtau, ' dtaumin = ', dtaumin | |
| 167 | ✗ | error stop 'orbit_sympl_init - error: dtau/dtaumin not integer' | |
| 168 | endif | ||
| 169 | |||
| 170 | ! Initialize symplectic integrator | ||
| 171 | 636 | si%pabs = z0(4) | |
| 172 | 636 | call encode_symplectic_state(f, z0, z) | |
| 173 | |||
| 174 | ! factor 1/sqrt(2) due to velocity normalisation different from other modules | ||
| 175 | call orbit_sympl_init(si, f, z, dtaumin/dsqrt(2d0), nint(dtau/dtaumin), & | ||
| 176 | 636 | rtol_init, mode_init) | |
| 177 | 636 | end subroutine init_sympl | |
| 178 | |||
| 179 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | subroutine reseed_sympl(si, f, z0, momentum_reference) |
| 180 | !> Re-seed one particle's symplectic state from standard GC coordinates. | ||
| 181 | !> | ||
| 182 | !> Unlike init_sympl/orbit_sympl_init this routine does not select the | ||
| 183 | !> module-global integration procedure. It is therefore safe to call from | ||
| 184 | !> an OpenMP worker after a local RK recovery step: every write is confined | ||
| 185 | !> to the caller's per-particle integrator and field objects. | ||
| 186 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 187 | type(field_can_t), intent(inout) :: f | ||
| 188 | real(dp), intent(in) :: z0(:) | ||
| 189 | real(dp), intent(in), optional :: momentum_reference | ||
| 190 | 1 | real(dp) :: z_seed(size(z0)) | |
| 191 | |||
| 192 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1 time.
|
6 | z_seed = z0 |
| 193 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (present(momentum_reference)) z_seed(4) = momentum_reference |
| 194 | 1 | si%pabs = z_seed(4) | |
| 195 | 1 | call encode_symplectic_state(f, z_seed, si%z) | |
| 196 | 1 | call get_val(f, si%z(4)) | |
| 197 | 1 | si%pthold = f%pth | |
| 198 | 1 | si%last_step_fraction = 1.d0 | |
| 199 | 1 | si%last_event_radial_residual = 0.d0 | |
| 200 | 1 | si%last_event_fraction_width = 0.d0 | |
| 201 | 1 | end subroutine reseed_sympl | |
| 202 | |||
| 203 | ! Seed the full-orbit Boris pusher from a GC start record, same normalization and | ||
| 204 | ! gyrophase reference as the symplectic GC seed so the two integrators start from | ||
| 205 | ! the identical particle. fo_init places the Larmor offset itself. | ||
| 206 |
7/8✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 3 times.
✓ Branch 7 taken 1 time.
|
10 | subroutine init_fo(fo, z0, dtaumin) |
| 207 | use orbit_fo_field, only: fo_eval_field | ||
| 208 | use params, only: orbit_coord | ||
| 209 | type(fo_state_t), intent(out) :: fo | ||
| 210 | real(dp), intent(in) :: z0(:) | ||
| 211 | real(dp), intent(in) :: dtaumin | ||
| 212 | real(dp) :: ro0_bar, mu, vpar_bar, vperp0 | ||
| 213 | real(dp) :: u_gc(3), Bmod, Acov(3), dA(3,3), dBmod(3), hcov(3) | ||
| 214 | |||
| 215 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (orbit_coord /= 1) error stop & |
| 216 | ✗ | 'full-orbit tracing supports only orbit_coord=1 (Boozer)' | |
| 217 |
2/4✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
|
1 | if (z0(1) <= 0d0 .or. z0(1) >= 1d0) error stop & |
| 218 | ✗ | 'full-orbit initialization requires 0 < s < 1' | |
| 219 | |||
| 220 | ! |B| at the guiding centre from the chartmap field (rho = sqrt(s)). | ||
| 221 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
|
4 | u_gc = [dsqrt(z0(1)), z0(2), z0(3)] |
| 222 | 1 | call fo_eval_field(u_gc, Acov, dA, Bmod, dBmod, hcov) | |
| 223 | 1 | mu = .5d0*z0(4)**2*(1.d0 - z0(5)**2)/Bmod*2d0 | |
| 224 | 1 | ro0_bar = ro0/dsqrt(2d0) | |
| 225 | 1 | vpar_bar = z0(4)*z0(5)*dsqrt(2d0) | |
| 226 | 1 | vperp0 = dsqrt(max(2d0*mu*Bmod, 0d0)) | |
| 227 | call fo_init(fo, z0(1:3), vpar_bar, vperp0, mu, 1d0, 1d0, & | ||
| 228 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1 | dtaumin/dsqrt(2d0), ro0_bar, z0(4)) |
| 229 | 1 | end subroutine init_fo | |
| 230 | |||
| 231 | ! Advance the full-orbit Boris pusher one normalized step and write back the | ||
| 232 | ! standard SIMPLE z(1:5): z(1)=guiding-centre s, z(2:3)=angles, z(4)=pabs, | ||
| 233 | ! z(5)=lambda. The step itself only locates the field (FO_OK / FO_LOCATE_FAIL); | ||
| 234 | ! the ONLY confinement loss is the guiding centre crossing s>=1, decided in | ||
| 235 | ! fo_to_gc (FO_LOSS -> ierr=2, counted fo_loss). A field-locate non-convergence | ||
| 236 | ! (FO_LOCATE_FAIL -> ierr=3, counted fo_fault) is a numerical fault, reported | ||
| 237 | ! but NEVER counted as a physical loss. | ||
| 238 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | subroutine orbit_timestep_fo(fo, z, ierr) |
| 239 | 1 | use diag_counters, only: count_event, EVT_FO_LOSS, EVT_FO_FAULT | |
| 240 | use orbit_fo_boris, only: FO_OK, FO_LOSS | ||
| 241 | type(fo_state_t), intent(inout) :: fo | ||
| 242 | real(dp), intent(inout) :: z(:) | ||
| 243 | integer, intent(out) :: ierr | ||
| 244 | type(fo_state_t) :: fo_start | ||
| 245 | 2 | real(dp) :: s, th, ph, vpar, z_start(size(z)) | |
| 246 | integer :: status | ||
| 247 | |||
| 248 | 1 | fo_start = fo | |
| 249 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1 time.
|
6 | z_start = z |
| 250 | 1 | call fo_step(fo, status) | |
| 251 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (status /= FO_OK) then |
| 252 | ! Inversion unresolved (near-axis below chartmap resolution, or a field-period | ||
| 253 | ! seam). NOT a physical loss; fo_step left the state at the last resolved | ||
| 254 | ! position. Roll back both public and hidden pusher state so warning mode | ||
| 255 | ! can consume a true held interval and retry from a meaningful position. | ||
| 256 | ✗ | fo = fo_start | |
| 257 | ✗ | z = z_start | |
| 258 | ✗ | ierr = ORBIT_FO_NUMERICAL | |
| 259 | ✗ | call count_event(EVT_FO_FAULT) | |
| 260 | ✗ | return | |
| 261 | end if | ||
| 262 | 1 | call fo_to_gc(fo, s, th, ph, vpar, status) | |
| 263 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (status == FO_LOSS) then |
| 264 | ✗ | z(1) = s | |
| 265 | ✗ | z(2) = th | |
| 266 | ✗ | z(3) = ph | |
| 267 | ✗ | z(4) = fo%pabs | |
| 268 | ✗ | z(5) = vpar/(z(4)*dsqrt(2d0)) | |
| 269 | ✗ | ierr = ORBIT_FO_LOSS | |
| 270 | ✗ | call count_event(EVT_FO_LOSS) | |
| 271 | ✗ | return | |
| 272 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | else if (status /= FO_OK) then |
| 273 | ✗ | fo = fo_start | |
| 274 | ✗ | z = z_start | |
| 275 | ✗ | ierr = ORBIT_FO_NUMERICAL | |
| 276 | ✗ | call count_event(EVT_FO_FAULT) | |
| 277 | ✗ | return | |
| 278 | end if | ||
| 279 | 1 | z(1) = s; z(2) = th; z(3) = ph | |
| 280 | 1 | z(4) = fo%pabs | |
| 281 | 1 | z(5) = vpar/(z(4)*dsqrt(2d0)) | |
| 282 | 1 | ierr = 0 | |
| 283 | 2 | end subroutine orbit_timestep_fo | |
| 284 | |||
| 285 | 1 | subroutine orbit_timestep_fo_bridge(fo, z, interval, ierr, continue_state) | |
| 286 | !> Advance one failed guiding-centre interval with the physical Cartesian | ||
| 287 | !> Lorentz pusher. This is a narrow warning-mode bridge for reference fields | ||
| 288 | !> that are already splined on the chartmap; normal GC trajectories never | ||
| 289 | !> enter it. The Boris substeps resolve at most 0.1 rad of gyro-rotation. | ||
| 290 | 1 | use orbit_fo_boris, only: fo_init_reference, fo_step, & | |
| 291 | fo_to_reference_gc, FO_OK, FO_LOSS | ||
| 292 | use orbit_fo_field, only: fo_eval_reference_field, & | ||
| 293 | fo_reference_field_available | ||
| 294 | use field_can_mod, only: integ_to_ref, ref_to_integ | ||
| 295 | type(fo_state_t), intent(inout) :: fo | ||
| 296 | real(dp), intent(inout) :: z(5) | ||
| 297 | real(dp), intent(in) :: interval | ||
| 298 | integer, intent(out) :: ierr | ||
| 299 | logical, intent(in), optional :: continue_state | ||
| 300 | real(dp), parameter :: max_rotation_parameter = 0.05_dp | ||
| 301 | integer, parameter :: max_substeps = 100000 | ||
| 302 | real(dp) :: z_start(5), u_gc(3), hcov(3), Bmod, mu | ||
| 303 | real(dp) :: ro0_bar, vpar_bar, vperp, remaining, h, vpar_out | ||
| 304 | integer :: status, nstep | ||
| 305 | logical :: continuing | ||
| 306 | |||
| 307 | 1 | z_start = z | |
| 308 | 1 | ierr = 2 | |
| 309 | 1 | continuing = .false. | |
| 310 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (present(continue_state)) continuing = continue_state |
| 311 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (.not. fo_reference_field_available()) return |
| 312 |
2/4✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
|
1 | if (z(4) <= 0.0_dp .or. abs(z(5)) > 1.0_dp) return |
| 313 | |||
| 314 | 1 | ro0_bar = ro0/sqrt(2.0_dp) | |
| 315 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if (.not. continuing) then |
| 316 | 1 | call integ_to_ref(z(1:3), u_gc) | |
| 317 |
2/4✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
|
1 | if (u_gc(1) <= 0.0_dp .or. u_gc(1) >= 1.0_dp) return |
| 318 | 1 | call fo_eval_reference_field(u_gc, hcov, Bmod, status) | |
| 319 |
2/4✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
|
1 | if (status /= 0 .or. Bmod <= 0.0_dp) return |
| 320 | |||
| 321 | 1 | mu = z(4)**2*(1.0_dp - z(5)**2)/Bmod | |
| 322 | 1 | vpar_bar = sqrt(2.0_dp)*z(4)*z(5) | |
| 323 | 1 | vperp = sqrt(max(2.0_dp*mu*Bmod, 0.0_dp)) | |
| 324 | call fo_init_reference(fo, u_gc, vpar_bar, vperp, mu, 1.0_dp, & | ||
| 325 | 1 | 1.0_dp, 0.0_dp, ro0_bar, z(4)) | |
| 326 | ✗ | else if (.not. fo%reference_field) then | |
| 327 | return | ||
| 328 | end if | ||
| 329 | |||
| 330 | 1 | remaining = abs(interval)/sqrt(2.0_dp) | |
| 331 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | do nstep = 1, max_substeps |
| 332 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | if (remaining <= epsilon(1.0_dp)*max(abs(interval), 1.0_dp)) exit |
| 333 | 1 | call fo_eval_reference_field(fo%u, hcov, Bmod, status) | |
| 334 |
2/4✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
|
1 | if (status /= 0 .or. Bmod <= 0.0_dp) then |
| 335 | ✗ | z = z_start | |
| 336 | ✗ | return | |
| 337 | end if | ||
| 338 | 1 | h = min(remaining, 2.0_dp*max_rotation_parameter*ro0_bar/Bmod) | |
| 339 | 1 | fo%dt = sign(h, interval) | |
| 340 | 1 | call fo_step(fo, status) | |
| 341 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (status /= FO_OK) then |
| 342 | ✗ | z = z_start | |
| 343 | ✗ | return | |
| 344 | end if | ||
| 345 | 3 | remaining = remaining - h | |
| 346 | end do | ||
| 347 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (remaining > epsilon(1.0_dp)*max(abs(interval), 1.0_dp)) then |
| 348 | ✗ | z = z_start | |
| 349 | ✗ | return | |
| 350 | end if | ||
| 351 | |||
| 352 | 1 | call fo_to_reference_gc(fo, u_gc, vpar_out, status) | |
| 353 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (status /= FO_OK .and. status /= FO_LOSS) then |
| 354 | ✗ | z = z_start | |
| 355 | ✗ | return | |
| 356 | end if | ||
| 357 | 1 | call ref_to_integ(u_gc, z(1:3)) | |
| 358 | 1 | z(4) = fo%pabs | |
| 359 | 1 | z(5) = vpar_out/(sqrt(2.0_dp)*z(4)) | |
| 360 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (status == FO_LOSS) then |
| 361 | ✗ | ierr = 1 | |
| 362 | else | ||
| 363 | 1 | ierr = 0 | |
| 364 | end if | ||
| 365 | 1 | end subroutine orbit_timestep_fo_bridge | |
| 366 | |||
| 367 | ✗ | subroutine timestep(self, s, th, ph, lam, ierr) | |
| 368 | type(tracer_t), intent(inout) :: self | ||
| 369 | real(dp), intent(inout) :: s, th, ph, lam | ||
| 370 | integer, intent(out) :: ierr | ||
| 371 | |||
| 372 | real(dp), dimension(5) :: z | ||
| 373 | |||
| 374 | ✗ | z(1) = s | |
| 375 | ✗ | z(2) = th | |
| 376 | ✗ | z(3) = ph | |
| 377 | ✗ | z(4) = 1d0 | |
| 378 | ✗ | z(5) = lam | |
| 379 | |||
| 380 | ✗ | call timestep_z(self, z, ierr) | |
| 381 | |||
| 382 | ✗ | s = z(1) | |
| 383 | ✗ | th = z(2) | |
| 384 | ✗ | ph = z(3) | |
| 385 | ✗ | lam = z(5) | |
| 386 | 1 | end subroutine timestep | |
| 387 | |||
| 388 | ✗ | subroutine timestep_z(self, z, ierr) | |
| 389 | use alpha_lifetime_sub, only : orbit_timestep_axis | ||
| 390 | |||
| 391 | type(tracer_t), intent(inout) :: self | ||
| 392 | real(dp), intent(inout) :: z(:) | ||
| 393 | integer, intent(out) :: ierr | ||
| 394 | |||
| 395 | ✗ | call orbit_timestep_axis(z, self%dtau, self%dtau, self%relerr, ierr) | |
| 396 | ✗ | end subroutine timestep_z | |
| 397 | |||
| 398 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 68996 times.
|
68996 | subroutine timestep_sympl_z(si, f, z, ierr) |
| 399 | type(symplectic_integrator_t), intent(inout) :: si | ||
| 400 | type(field_can_t), intent(inout) :: f | ||
| 401 | real(dp), intent(inout) :: z(:) | ||
| 402 | integer, intent(out) :: ierr | ||
| 403 | |||
| 404 |
2/4✓ Branch 0 taken 68996 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 68996 times.
|
68996 | if (z(1) < 0.0 .or. z(1) > 1.0) then |
| 405 | ✗ | ierr = 1 | |
| 406 | ✗ | return | |
| 407 | end if | ||
| 408 | |||
| 409 | 68996 | call orbit_timestep_sympl(si, f, ierr) | |
| 410 | |||
| 411 |
2/2✓ Branch 0 taken 206988 times.
✓ Branch 1 taken 68996 times.
|
275984 | z(1:3) = si%z(1:3) |
| 412 | 68996 | z(5) = f%vpar/(si%pabs*dsqrt(2d0)) ! alambda | |
| 413 | |||
| 414 | end subroutine timestep_sympl_z | ||
| 415 | |||
| 416 | ✗ | end module simple | |
| 417 |