GCC Code Coverage Report

Directory: ./
File: src/orbit_fo_boris.f90
Date: 2026-10-22 05:02:02
Coverage Exec Excl Total
Lines: 87.0% 280 0 322
Functions: 96.4% 27 0 28
Branches: 74.4% 189 0 254
Line Branch Exec Source
1 module orbit_fo_boris
2 ! Full-orbit (FO) Boris pusher: a gyro-resolved classical charged particle in
3 ! Cartesian (x, v), the full-orbit counterpart to SIMPLE's guiding-center (GC)
4 ! model and the ASCOT-style reference for it. The particle advances by the
5 ! explicit Boris drift-rotate-drift: the magnetic rotation is exact for constant
6 ! B over a step, the kinetic metric is the identity, the geodesic terms vanish,
7 ! and the magnetic axis is an ordinary point.
8 !
9 ! Field and geometry come from the chartmap (the Cartesian-side representation):
10 ! at the Cartesian point we invert to the logical chart u=(rho, theta_B, phi_B)
11 ! with the chartmap forward map (rho=sqrt(s)), evaluate the production Boozer
12 ! flux potential there (fo_eval_field: A_theta(s), A_phi(s), |B|, d|B|/du), and
13 ! build the Cartesian field as B = curl A so B^s = 0 exactly (B tangent to the
14 ! flux surface) with grad|B|_cart = Jc^{-T} d|B|/du. The magnetic axis (rho->0)
15 ! is healed by a pseudo-Cartesian chart w=(X,Y,phi)=(rho cos th, rho sin th, phi):
16 ! the polar basis dx/dtheta ~ rho makes det(Jc)->0, but the (X,Y) chart Jacobian
17 ! Jw is regular through rho=0, so the inverse Newton and the field assembly both
18 ! stay well-conditioned there. The chartmap also owns the loss boundary: the
19 ! guiding-centre crossing rho>=1 (out of the s<1 plasma) is the ONLY confinement
20 ! loss. A field-locate non-convergence is a numerical fault, retried/reported,
21 ! never a loss.
22 use, intrinsic :: iso_fortran_env, only: dp => real64
23 use reference_coordinates, only: ref_coords
24 use orbit_fo_field, only: fo_eval_field, fo_eval_reference_field
25 use libneo_coordinates, only: chartmap_coordinate_system_t, chartmap_from_cyl_ok
26 implicit none
27 private
28
29 real(dp), parameter :: c = 1.0_dp
30 real(dp), parameter :: twopi = 8.0_dp*atan(1.0_dp)
31
32 ! Inverse-Newton classification: a converged locate has residual below
33 ! NEWTON_ACCEPT_TOL; a loosely converged point counts as the edge only when its
34 ! residual is below EDGE_FRAC of a radial cell (a genuine gyro-overshoot loss sits
35 ! within a Larmor radius of RHO_EDGE), otherwise it is a stalled interior fault.
36 real(dp), parameter :: NEWTON_ACCEPT_TOL = 1.0e-6_dp, RHO_EDGE = 1.0_dp, &
37 EDGE_FRAC = 0.05_dp
38
39 ! A guiding-centre loss (u_gc >= 1) is confirmed only when the robustly-located
40 ! particle radius u_p is within this gap of the edge: the GC and particle differ by
41 ! one Larmor radius (rho*/a ~ 0.005-0.01 here), so 0.05 is several Larmor radii of
42 ! margin and rejects field-period-seam reconstruction glitches that put a mid-radius
43 ! particle's reconstructed GC spuriously at rho >= 1.
44 real(dp), parameter :: GC_PARTICLE_GAP = 0.05_dp
45
46 ! cart_field / locate status: regular interior point, physical edge loss, or a
47 ! numerical locate fault (NOT a loss).
48 integer, parameter, public :: FO_OK = 0, FO_LOSS = 1, FO_LOCATE_FAIL = 2
49
50 public :: fo_state_t, fo_init, fo_init_reference, fo_step, fo_energy, &
51 fo_mu, fo_to_gc, fo_to_reference_gc, accept_or_fail, fo_step_rkng, &
52 fo_step_rkng_adaptive
53
54
55 type :: fo_state_t
56 real(dp) :: x(3) = 0.0_dp ! Cartesian position (scaled cm)
57 real(dp) :: v(3) = 0.0_dp ! Cartesian velocity (normalized)
58 real(dp) :: u(3) = 0.0_dp ! last logical (rho, theta_B, phi_B)
59 real(dp) :: mu = 0.0_dp ! guiding-centre magnetic moment (diagnostic)
60 real(dp) :: dt = 0.0_dp
61 real(dp) :: mass = 1.0_dp
62 real(dp) :: charge = 1.0_dp
63 real(dp) :: ro0 = 1.0_dp
64 real(dp) :: pabs = 0.0_dp ! normalized speed (carried for z(4) write-back)
65 logical :: reference_field = .false.
66 end type fo_state_t
67
68 ! State for the RKNG right-hand side. fortnum's rkng interface passes only
69 ! (t, y, yp, ypp) plus an unlimited-polymorphic ctx, while the force needs the
70 ! warm-start logical coordinate and the reference-field flag, and it updates
71 ! that warm start as it goes. Carrying them here mirrors how the guiding-centre
72 ! quasi module holds si and f, and keeps the interface fortnum expects.
73 type(fo_state_t) :: rkng_st
74 integer :: rkng_status = FO_OK
75 !$omp threadprivate(rkng_st, rkng_status)
76
77 contains
78
79 ! Cartesian (wedge) -> logical chart (rho, theta_B, phi_B). Warm damped Newton on
80 ! the chartmap forward map x(u)=evaluate_cart(u) from the carried guess (a Larmor
81 ! step away, 1-2 iters, thread-safe read-only spline eval); on stall -- e.g. the
82 ! guess went stale across a field-period seam -- fall back to the chartmap's robust
83 ! multi-seed from_cart, which seeds zeta across [0, 2pi/nfp). status: FO_OK (located,
84 ! rho reported through u for the caller's guiding-centre loss test) or FO_LOCATE_FAIL
85 ! (no converged root -> numerical fault, never itself a loss).
86 1421432 subroutine invert_cart_warm(x, u_guess, u, status)
87 real(dp), intent(in) :: x(3), u_guess(3)
88 real(dp), intent(out) :: u(3)
89 integer, intent(out) :: status
90 real(dp) :: xc(3), Jc(3,3), rn
91 integer :: ierr
92
93 1421432 call invert_warm_newton(x, u_guess, u, status)
94
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1421432 times.
 1421432 if (status /= FO_LOCATE_FAIL) return
95 select type (cs => ref_coords) ! robust multi-seed fallback
96 class is (chartmap_coordinate_system_t)
97 ✗ call cs%from_cart(x, u, ierr)
98 class default
99 return
100 end select
101 ✗ if (ierr /= chartmap_from_cyl_ok) return ! genuine no-root: keep LOCATE_FAIL
102 ! Re-verify the fallback root with the same residual-vs-radial-cell criterion as
103 ! the warm path: from_cart clamps rho to [0,1], so a seam point it cannot solve
104 ! comes back pinned at rho=1 with a large residual. Accepting that as the edge
105 ! fakes a loss from mid-radius, so classify by the actual residual, not by rho.
106 ✗ call ref_coords%evaluate_cart(u, xc)
107 ✗ call ref_coords%covariant_basis(u, Jc)
108 ✗ rn = sqrt((xc(1) - x(1))**2 + (xc(2) - x(2))**2 + (xc(3) - x(3))**2)
109 status = accept_or_fail(u(1), rn, radial_scale(Jc), NEWTON_ACCEPT_TOL, RHO_EDGE, &
110 ✗ u_guess(1))
111 end subroutine invert_cart_warm
112
113 ! Cartesian (wedge) -> logical Newton. Seed rho and theta from the carried guess
114 ! (they do not jump between substeps) and the toroidal coordinate from the in-wedge
115 ! geometric angle atan2(y,x). The carried Boozer phi lives on the global
116 ! multi-period sheet and goes a full period 2*pi/nfp stale across a field-period
117 ! seam; the geometric angle is always in-wedge and differs from logical phi only by
118 ! the Boozer shift O(0.1 rad). One robust seed, no seam special case. (A warm phi
119 ! seed is faster away from seams but cannot be trusted at them: evaluate_cart wraps
120 ! phi mod 2*pi/nfp, so a stale guess can converge to a clamped-edge root and fake a
121 ! loss.) On stall the caller (invert_cart_warm) runs the multi-seed from_cart.
122 1421432 subroutine invert_warm_newton(x, u_guess, u, status)
123 real(dp), intent(in) :: x(3), u_guess(3)
124 real(dp), intent(out) :: u(3)
125 integer, intent(out) :: status
126
127
2/2
✓ Branch 0 taken 4264296 times.
✓ Branch 1 taken 1421432 times.
 5685728 call newton_from(x, [u_guess(1), u_guess(2), atan2(x(2), x(1))], u, status)
128 1421432 end subroutine invert_warm_newton
129
130 ! Damped Newton on the chartmap forward map x(u)=evaluate_cart(u) from an explicit
131 ! seed. Iterate in the pseudo-Cartesian chart w=(X,Y,phi)=(rho cos th, rho sin th,
132 ! phi): the polar (rho,theta) Newton is singular at the axis (dx/dtheta ~ rho,
133 ! det(Jc)->0), whereas the (X,Y) step stays regular and crosses the axis without
134 ! the reflect hack. w_to_u recovers rho>=0, theta=atan2 automatically.
135 1421432 subroutine newton_from(x, u_seed, u, status)
136 real(dp), intent(in) :: x(3), u_seed(3)
137 real(dp), intent(out) :: u(3)
138 integer, intent(out) :: status
139 integer, parameter :: maxit = 30, maxls = 30
140 ! The forward map is a deterministic spline, so a damped Newton converges to
141 ! ~machine precision; tol targets that. accept_or_fail classifies a stall.
142 real(dp), parameter :: tol = 1.0e-9_dp
143 real(dp) :: xc(3), Jc(3,3), Jw(3,3), Jinv(3,3), res(3)
144 real(dp) :: w(3), wt(3), ut(3), dw(3), cth, sth, rho, rn, rnew, alpha
145 integer :: it, ls, i
146
147 1421432 u = u_seed
148 1421432 w(1) = u(1)*cos(u(2)); w(2) = u(1)*sin(u(2)); w(3) = u(3)
149 1421432 call ref_coords%evaluate_cart(u, xc)
150
2/2
✓ Branch 0 taken 4264296 times.
✓ Branch 1 taken 1421432 times.
 5685728 res = xc - x
151 1421432 rn = sqrt(res(1)**2 + res(2)**2 + res(3)**2)
152
1/2
✓ Branch 0 taken 4413788 times.
✗ Branch 1 not taken.
 4413788 do it = 1, maxit
153
2/2
✓ Branch 0 taken 1421432 times.
✓ Branch 1 taken 2992356 times.
 4413788 if (rn < tol) then
154 1421432 status = accept_or_fail(u(1), rn, 0.0_dp, NEWTON_ACCEPT_TOL, RHO_EDGE, u_seed(1))
155 1421432 return
156 end if
157 2992356 call ref_coords%covariant_basis(u, Jc)
158 2992356 call pseudocart_basis(u, Jc, Jw, cth, sth, rho)
159
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2992356 times.
 2992356 if (.not. jacobian_ok(Jw)) then ! genuinely degenerate (off the chart)
160 ✗ status = FO_LOCATE_FAIL; return
161 end if
162 2992356 call inv3(Jw, Jinv)
163
2/2
✓ Branch 1 taken 8977068 times.
✓ Branch 2 taken 2992356 times.
 11969424 do i = 1, 3
164 11969424 dw(i) = -(Jinv(i,1)*res(1) + Jinv(i,2)*res(2) + Jinv(i,3)*res(3))
165 end do
166 ! Backtracking line search: Newton is not monotonic for a finite offset.
167 ! A trial that overshoots past rho=1 must be rejected, not evaluated: the
168 ! forward map clamps rho to the grid edge, so a past-edge trial returns a
169 ! point ON the edge whose residual can be smaller than the current interior
170 ! point -- the line search would accept it and the next step stalls on the
171 ! flat clamped region (the failure seen for mid-radius orbits crossing a
172 ! field-period seam, where the seam-corrupted step points outward). w_to_u
173 ! gives the trial rho before the clamp, so reject rho > 1 and keep backtracking
174 ! toward the true interior root. Loss is decided only on the guiding centre.
175 alpha = 1.0_dp
176
1/2
✓ Branch 0 taken 2992356 times.
✗ Branch 1 not taken.
 2992356 do ls = 1, maxls
177
2/2
✓ Branch 0 taken 8977068 times.
✓ Branch 1 taken 2992356 times.
 11969424 wt = w + alpha*dw
178 2992356 call w_to_u(wt, ut)
179
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2992356 times.
 2992356 if (ut(1) > 1.0_dp) then
180 ✗ alpha = 0.5_dp*alpha
181 ✗ cycle
182 end if
183 2992356 call ref_coords%evaluate_cart(ut, xc)
184
2/2
✓ Branch 0 taken 8977068 times.
✓ Branch 1 taken 2992356 times.
 11969424 res = xc - x
185 2992356 rnew = sqrt(res(1)**2 + res(2)**2 + res(3)**2)
186
1/2
✓ Branch 0 taken 2992356 times.
✗ Branch 1 not taken.
 2992356 if (rnew < rn) exit
187 2992356 alpha = 0.5_dp*alpha
188 end do
189
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2992356 times.
 2992356 if (rnew >= rn) then ! line search could not improve -> stalled at the floor
190 status = accept_or_fail(u(1), rn, radial_scale(Jc), NEWTON_ACCEPT_TOL, RHO_EDGE, &
191 ✗ u_seed(1))
192 ✗ return
193 end if
194 2992356 w = wt
195 2992356 u = ut
196 5984712 rn = rnew
197 end do
198 status = accept_or_fail(u(1), rn, radial_scale(Jc), NEWTON_ACCEPT_TOL, RHO_EDGE, &
199 ✗ u_seed(1))
200 end subroutine newton_from
201
202 ! Length of one unit-rho radial step |dx/drho| = |Jc(:,1)|, the chart scale used to
203 ! judge a stalled Newton: a residual that is a small fraction of a radial cell means
204 ! the target is essentially at the edge, a large fraction means an interior stall.
205 ✗ pure real(dp) function radial_scale(Jc) result(s)
206 real(dp), intent(in) :: Jc(3,3)
207 ✗ s = sqrt(Jc(1,1)**2 + Jc(2,1)**2 + Jc(3,1)**2)
208 end function radial_scale
209
210 ! Classify a finished Newton by the Cartesian residual rn, judged against the local
211 ! radial cell |dx/drho| (scale), not an absolute length. A point is located when rn
212 ! is at machine tolerance OR a small fraction (EDGE_FRAC) of a radial cell. The
213 ! relative test is what makes this scale-correct: on a reactor-size chartmap
214 ! (positions ~1e3 cm) and near the magnetic axis, where the chart is barely resolved
215 ! (innermost rho grid point ~1e-3) and |dx/drho| is large, the residual floors well
216 ! above any fixed accept_tol while the point still sits a tiny fraction of a cell
217 ! from its target. A genuine stall -- a residual that is a sizable fraction of a
218 ! radial cell -- is a fault, EXCEPT when it clamps to the edge (rho ~ rho_edge) for a
219 ! marker whose warm guess rho_guess was already within GC_PARTICLE_GAP of the edge:
220 ! that is a marker leaving the plasma, whose true position is past rho=1 where the
221 ! forward map cannot represent it, so the inverse stalls on the clamped edge. Accept
222 ! it as located at the edge so the push continues on the clamped-edge field and
223 ! fo_to_gc decides the loss on the guiding centre -- never a blanket confined fault
224 ! that would silently keep an exiting marker. A mid-radius seam glitch that clamps to
225 ! rho=1 has rho_guess well inside, fails the guard, and stays a fault, so it is never
226 ! turned into a spurious loss.
227 1421436 pure integer function accept_or_fail(rho, rn, scale, accept_tol, rho_edge, rho_guess) &
228 result(status)
229 real(dp), intent(in) :: rho, rn, scale, accept_tol, rho_edge, rho_guess
230
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1421433 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
 1421436 if (rn < accept_tol .or. (scale > 0.0_dp .and. rn < EDGE_FRAC*scale)) then
231 status = FO_OK
232
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 1 time.
 3 else if (rho >= rho_edge - GC_PARTICLE_GAP .and. &
233 rho_guess >= rho_edge - GC_PARTICLE_GAP) then
234 status = FO_OK
235 else
236 1421436 status = FO_LOCATE_FAIL
237 end if
238 1421436 end function accept_or_fail
239
240 ! Geometric field period 2*pi/nfp; the device is exactly nfp-fold symmetric about
241 ! Z, so a rotation by this angle maps one field period onto the next.
242 1421432 real(dp) function field_period()
243 integer :: nfp
244 1421432 nfp = 1
245 select type (cs => ref_coords)
246 class is (chartmap_coordinate_system_t)
247 ✗ nfp = cs%num_field_periods
248 end select
249 1421432 field_period = twopi/real(max(nfp, 1), dp)
250 1421432 end function field_period
251
252 2842864 pure function rotz(v, ca, sa) result(w)
253 real(dp), intent(in) :: v(3), ca, sa
254 real(dp) :: w(3)
255 2842864 w(1) = ca*v(1) - sa*v(2)
256 2842864 w(2) = sa*v(1) + ca*v(2)
257 2842864 w(3) = v(3)
258 end function rotz
259
260 ! Map a global Cartesian point into the fundamental field-period wedge by an
261 ! integer rotation about Z. The chartmap stores geometry only on one period, so
262 ! the inversion and field evaluation run in the wedge; (ca, sa) rotate the wedge
263 ! field vectors back to the global frame. This is what lets the Cartesian inverse
264 ! converge to machine precision on a multi-period (nfp>1) device.
265 1421432 subroutine to_wedge(x, xw, ca, sa)
266 real(dp), intent(in) :: x(3)
267 real(dp), intent(out) :: xw(3), ca, sa
268 real(dp) :: phi, period, alpha
269 1421432 period = field_period()
270 1421432 phi = atan2(x(2), x(1))
271
2/2
✓ Branch 0 taken 83256 times.
✓ Branch 1 taken 1338176 times.
 1421432 alpha = period*floor(phi/period)
272 1421432 ca = cos(alpha); sa = sin(alpha)
273 1421432 xw = rotz(x, ca, -sa) ! rotate by -alpha into the wedge
274 1421432 end subroutine to_wedge
275
276 ! Cartesian B, |B|, grad|B| at logical u from the chartmap field (field_can) and
277 ! geometry (ref_coords): B^i = |B| g^{ij} h_j, B_cart = Jc B^i; grad|B| covariant
278 ! d|B|/du -> Cartesian by Jc^{-T}. Jc returned for downstream Larmor offsets.
279 2842853 subroutine field_at_logical(u, Bvec, Bmod, gradB, Jc, status, reference_field)
280 real(dp), intent(in) :: u(3)
281 real(dp), intent(out) :: Bvec(3), Bmod, gradB(3), Jc(3,3)
282 integer, intent(out) :: status
283 logical, intent(in) :: reference_field
284 real(dp) :: ue(3), Acov(3), dA(3,3), dBmod(3), hcov(3)
285 real(dp) :: Jw(3,3), Jwinv(3,3), cth, sth, rho, detJw, Bw(3), dBw(3), Bn
286 real(dp) :: hw(3)
287 integer :: i
288
289 ! A particle may gyro-excurse a Larmor radius past s=1; evaluate the field at
290 ! the clamped edge there (field_can is undefined past the last closed surface).
291 1421432 ue = u
292 1421432 ue(1) = min(ue(1), 1.0_dp - 1.0e-9_dp)
293
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 1421421 times.
 1421432 if (reference_field) then
294 11 call fo_eval_reference_field(ue, hcov, Bmod, status)
295
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if (status /= 0) then
296 ✗ status = FO_LOCATE_FAIL
297 11 return
298 end if
299 11 call ref_coords%covariant_basis(ue, Jc)
300 ! Assemble the covector in the regular rotating (X,Y,phi) frame. Directly
301 ! inverting the polar Jc becomes ill-conditioned as rho->0 even though the
302 ! Cartesian field is smooth. The same frame is used by the chart inverse;
303 ! to_wedge supplies the exact field-period rotation on return to global
304 ! Cartesian coordinates.
305 11 call pseudocart_basis(ue, Jc, Jw, cth, sth, rho)
306
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if (.not. jacobian_ok(Jw)) then
307 ✗ status = FO_LOCATE_FAIL
308 ✗ return
309 end if
310 11 hw(1) = cth*hcov(1) - (sth/rho)*hcov(2)
311 11 hw(2) = sth*hcov(1) + (cth/rho)*hcov(2)
312 11 hw(3) = hcov(3)
313 11 call inv3(Jw, Jwinv)
314 11 Bvec = matmul(transpose(Jwinv), hw)
315
4/6
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 33 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 33 times.
 44 Bn = norm2(Bvec)
316
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if (Bn <= 0.0_dp) then
317 ✗ status = FO_LOCATE_FAIL
318 ✗ return
319 end if
320
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 11 times.
 44 Bvec = Bvec*(Bmod/Bn)
321 11 gradB = 0.0_dp
322 11 status = FO_OK
323 11 return
324 end if
325 1421421 call fo_eval_field(ue, Acov, dA, Bmod, dBmod, hcov)
326 1421421 call ref_coords%covariant_basis(ue, Jc)
327 ! Heal the axis: work in the pseudo-Cartesian chart w=(X,Y,phi). The polar
328 ! basis column Jc(:,2)=dx/dtheta ~ rho makes det(Jc)->0 at the axis; the (X,Y)
329 ! chart Jacobian Jw is regular there (det(Jw)=det(Jc)/rho, same sign).
330 1421421 call pseudocart_basis(ue, Jc, Jw, cth, sth, rho)
331
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1421421 times.
 1421421 if (.not. jacobian_ok(Jw)) then ! genuinely degenerate (off the chart)
332 ✗ status = FO_LOCATE_FAIL; return
333 end if
334 detJw = Jw(1,1)*(Jw(2,2)*Jw(3,3) - Jw(2,3)*Jw(3,2)) &
335 - Jw(1,2)*(Jw(2,1)*Jw(3,3) - Jw(2,3)*Jw(3,1)) &
336 1421421 + Jw(1,3)*(Jw(2,1)*Jw(3,2) - Jw(2,2)*Jw(3,1))
337 ! B = curl A with the Boozer flux-function potential A = (0, A_theta(s),
338 ! A_phi(s)): B^s = d_theta A_phi - d_phi A_theta = 0 EXACTLY (B tangent to the
339 ! flux surface). In the regular (X,Y,phi) chart the contravariant components
340 ! stay finite through the axis (det(Jw) is bounded, and dA_theta/drho ~ rho
341 ! cancels the surviving 1/rho). The signed det(Jw) keeps the left-handed chart
342 ! orientation: the unsigned sqrt(|det|) flips B and sends trapped bananas
343 ! outward. B^X = sin th dA_phi/drho /detJw, B^Y = -cos th dA_phi/drho /detJw,
344 ! B^phi = (dA_theta/drho)/rho /detJw. Renormalize to |B|.
345 1421421 Bw(1) = sth*dA(3,1)/detJw
346 1421421 Bw(2) = -cth*dA(3,1)/detJw
347 1421421 Bw(3) = (dA(2,1)/rho)/detJw
348
2/2
✓ Branch 0 taken 4264263 times.
✓ Branch 1 taken 1421421 times.
 5685684 do i = 1, 3
349 5685684 Bvec(i) = Jw(i,1)*Bw(1) + Jw(i,2)*Bw(2) + Jw(i,3)*Bw(3)
350 end do
351 1421421 Bn = sqrt(Bvec(1)**2 + Bvec(2)**2 + Bvec(3)**2)
352
2/2
✓ Branch 0 taken 4264263 times.
✓ Branch 1 taken 1421421 times.
 5685684 Bvec = Bvec*(Bmod/max(Bn, 1.0e-30_dp))
353 ! grad|B| in Cartesian: Jw^{-T} d|B|/dw, with d|B|/dw mapped from d|B|/du by the
354 ! pseudo-Cartesian chain rule (d|B|/dtheta ~ rho cancels the 1/rho).
355 1421421 dBw(1) = cth*dBmod(1) - (sth/rho)*dBmod(2)
356 1421421 dBw(2) = sth*dBmod(1) + (cth/rho)*dBmod(2)
357 1421421 dBw(3) = dBmod(3)
358 1421421 call inv3(Jw, Jwinv)
359
2/2
✓ Branch 1 taken 4264263 times.
✓ Branch 2 taken 1421421 times.
 5685684 do i = 1, 3
360 5685684 gradB(i) = Jwinv(1,i)*dBw(1) + Jwinv(2,i)*dBw(2) + Jwinv(3,i)*dBw(3)
361 end do
362 1421421 status = FO_OK
363 end subroutine field_at_logical
364
365 ! Cartesian B vector, |B|, and grad|B| at Cartesian x, from the chartmap field.
366 ! status: FO_OK (located, u_out valid) or FO_LOCATE_FAIL (numerical inversion
367 ! fault). On fault Bvec etc. are undefined and the caller must not push. Loss is
368 ! not decided here -- the field is defined through the clamped edge, and only the
369 ! guiding-centre crossing rho>=1 in fo_to_gc is a confinement loss.
370 ! Lorentz acceleration in general Runge-Kutta-Nystrom form.
371 !
372 ! Full-orbit motion is xdd = (q/mc) xd x B(x), which is genuinely
373 ! y'' = f(t, y, y') -- general RKN form. Unlike the guiding-centre case this
374 ! needs no reformulation and no second derivatives of the field: it is the one
375 ! place a Nystrom-family method applies to SIMPLE directly, and so the most
376 ! faithful test of the original suggestion.
377 !
378 ! The warm-start logical coordinate is advanced on every force evaluation, as
379 ! in fo_step, so the inversion stays close to its previous solution. A failed
380 ! inversion is recorded in rkng_status and the acceleration is returned as
381 ! zero: fortnum has no way to abort a fixed-step integration mid-flight, so the
382 ! failure is reported after the step rather than propagated through it.
383
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1073274 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1073274 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1073274 times.
 1073274 subroutine fo_rkng_force(t, y, yp, ypp, ctx)
384 real(dp), intent(in) :: t
385 real(dp), intent(in) :: y(:), yp(:)
386 real(dp), intent(out) :: ypp(:)
387 class(*), intent(in), optional :: ctx
388 real(dp) :: Bvec(3), Bmod, gradB(3), u(3), qcm
389 integer :: status
390
391
2/2
✓ Branch 0 taken 3219822 times.
✓ Branch 1 taken 1073274 times.
 4293096 ypp = 0.0_dp
392
1/2
✓ Branch 0 taken 1073274 times.
✗ Branch 1 not taken.
 1073274 if (rkng_status /= FO_OK) return
393
394 call cart_field(y(1:3), rkng_st%u, Bvec, Bmod, gradB, u, status, &
395
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1073274 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1073274 rkng_st%reference_field)
396
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1073274 times.
 1073274 if (status /= FO_OK) then
397 ✗ rkng_status = status
398 ✗ return
399 end if
400
2/2
✓ Branch 0 taken 3219822 times.
✓ Branch 1 taken 1073274 times.
 4293096 rkng_st%u = u
401
402 1073274 qcm = rkng_st%charge/(c*rkng_st%ro0*rkng_st%mass)
403
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1073274 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 5 taken 3219822 times.
✓ Branch 6 taken 1073274 times.
 5366370 ypp(1:3) = qcm*cross(yp(1:3), Bvec)
404 end subroutine fo_rkng_force
405
406 ! One full-orbit step with a general Runge-Kutta-Nystrom method, as an
407 ! alternative to the Boris pusher over the same dt and the same field path.
408 !
409 ! Boris is second order and volume preserving, and conserves energy to
410 ! round-off in a static magnetic field because its rotation is exact. RKNG is
411 ! higher order but not structure preserving, so the comparison to make is
412 ! energy drift at matched cost, not energy drift alone.
413 12000 subroutine fo_step_rkng(st, status, nsub)
414 use fortnum_ode_tdrk, only: rkng_integrate_fixed
415 use fortnum_status, only: fortnum_status_t, FORTNUM_OK
416 type(fo_state_t), intent(inout) :: st
417 integer, intent(out) :: status
418 integer, intent(in), optional :: nsub
419 real(dp) :: yend(3), ypend(3)
420 integer :: nfev, nsteps
421 type(fortnum_status_t) :: fstat
422
423 12000 nsteps = 1
424 12000 if (present(nsub)) nsteps = max(1, nsub)
425
426 12000 rkng_st = st
427 12000 rkng_status = FO_OK
428
429 call rkng_integrate_fixed(fo_rkng_force, 0.0_dp, st%dt, st%x, st%v, &
430 12000 nsteps, yend, ypend, nfev, fstat)
431
432
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12000 times.
 12000 if (rkng_status /= FO_OK) then
433 ✗ status = rkng_status ! leave st at the last resolved state
434 ✗ return
435 end if
436
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12000 times.
 12000 if (fstat%code /= FORTNUM_OK) then
437 ✗ status = FO_LOCATE_FAIL
438 ✗ return
439 end if
440
441
2/2
✓ Branch 0 taken 36000 times.
✓ Branch 1 taken 12000 times.
 48000 st%x = yend
442
2/2
✓ Branch 0 taken 36000 times.
✓ Branch 1 taken 12000 times.
 48000 st%v = ypend
443
2/2
✓ Branch 0 taken 36000 times.
✓ Branch 1 taken 12000 times.
 48000 st%u = rkng_st%u
444 12000 status = FO_OK
445
1/2
✓ Branch 0 taken 12000 times.
✗ Branch 1 not taken.
 12000 end subroutine fo_step_rkng
446
447 ! The same RKNG method under tolerance control rather than a fixed substep
448 ! count. rtol/atol replace nsub as the accuracy knob, which is what lets the
449 ! full-orbit arena be swept the same way the guiding-centre one is.
450 !
451 ! h_carry passes the accepted step size from one macro-step to the next: a
452 ! gyro-orbit's timescale barely changes between steps, so restarting the
453 ! controller from scratch every dt would waste the rejections it already paid
454 ! for. Pass 0 on the first call and feed the returned value back afterwards.
455 12000 subroutine fo_step_rkng_adaptive(st, rtol, atol, h_carry, status, nfev)
456 use fortnum_ode_tdrk, only: rkng_integrate_adaptive
457 use fortnum_status, only: fortnum_status_t, FORTNUM_OK
458 type(fo_state_t), intent(inout) :: st
459 real(dp), intent(in) :: rtol, atol
460 real(dp), intent(inout) :: h_carry
461 integer, intent(out) :: status
462 integer, intent(out), optional :: nfev
463 real(dp) :: yend(3), ypend(3), hlast
464 integer :: nf, naccept, nreject
465 type(fortnum_status_t) :: fstat
466
467 12000 rkng_st = st
468 12000 rkng_status = FO_OK
469
470 call rkng_integrate_adaptive(fo_rkng_force, 0.0_dp, st%dt, st%x, st%v, &
471 12000 rtol, atol, h_carry, yend, ypend, hlast, nf, naccept, nreject, fstat)
472
473
1/2
✓ Branch 0 taken 12000 times.
✗ Branch 1 not taken.
 12000 if (present(nfev)) nfev = nf
474
475
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12000 times.
 12000 if (rkng_status /= FO_OK) then
476 ✗ status = rkng_status ! leave st at the last resolved state
477 ✗ return
478 end if
479
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12000 times.
 12000 if (fstat%code /= FORTNUM_OK) then
480 ✗ status = FO_LOCATE_FAIL
481 ✗ return
482 end if
483
484
2/2
✓ Branch 0 taken 36000 times.
✓ Branch 1 taken 12000 times.
 48000 st%x = yend
485
2/2
✓ Branch 0 taken 36000 times.
✓ Branch 1 taken 12000 times.
 48000 st%v = ypend
486
2/2
✓ Branch 0 taken 36000 times.
✓ Branch 1 taken 12000 times.
 48000 st%u = rkng_st%u
487 12000 h_carry = hlast
488 12000 status = FO_OK
489 12000 end subroutine fo_step_rkng_adaptive
490
491 1133277 subroutine cart_field(x, u_guess, Bvec, Bmod, gradB, u_out, status, &
492 reference_field)
493 real(dp), intent(in) :: x(3), u_guess(3)
494 real(dp), intent(out) :: Bvec(3), Bmod, gradB(3), u_out(3)
495 integer, intent(out) :: status
496 logical, intent(in) :: reference_field
497 real(dp) :: xw(3), ca, sa, u(3), Jc(3,3), Bw(3), gw(3)
498
499 1133277 call to_wedge(x, xw, ca, sa)
500 1133277 call invert_cart_warm(xw, u_guess, u, status)
501
1/2
✓ Branch 0 taken 1133277 times.
✗ Branch 1 not taken.
 1133277 if (status /= FO_OK) return
502 1133277 u_out = u
503 1133277 call field_at_logical(u, Bw, Bmod, gw, Jc, status, reference_field)
504
1/2
✓ Branch 0 taken 1133277 times.
✗ Branch 1 not taken.
 1133277 if (status /= FO_OK) return
505 1133277 Bvec = rotz(Bw, ca, sa) ! wedge field vector -> global frame
506 1133277 gradB = rotz(gw, ca, sa)
507 end subroutine cart_field
508
509 ! Logical chart of a Cartesian point and the local covariant frame, for seeding
510 ! and Larmor offsets. status as in cart_field. u_guess warm-starts the inversion.
511 288155 subroutine locate(x, u_guess, u_out, bhat, eperp, Bmod, status, reference_field)
512 real(dp), intent(in) :: x(3), u_guess(3)
513 real(dp), intent(out) :: u_out(3), bhat(3), eperp(3), Bmod
514 integer, intent(out) :: status
515 logical, intent(in) :: reference_field
516 real(dp) :: xw(3), ca, sa, u(3), Jc(3,3), Bw(3), gw(3), bw_hat(3), ew(3)
517
518 288155 call to_wedge(x, xw, ca, sa)
519 288155 call invert_cart_warm(xw, u_guess, u, status)
520
1/2
✓ Branch 0 taken 288155 times.
✗ Branch 1 not taken.
 288155 if (status /= FO_OK) return
521 288155 u_out = u
522 ! Same axis-healed field assembly as the push (field_at_logical), so the seed
523 ! frame and the orbit-step field are identical.
524 288155 call field_at_logical(u, Bw, Bmod, gw, Jc, status, reference_field)
525
1/2
✓ Branch 0 taken 288155 times.
✗ Branch 1 not taken.
 288155 if (status /= FO_OK) return
526
2/2
✓ Branch 0 taken 864465 times.
✓ Branch 1 taken 288155 times.
 1152620 bw_hat = Bw/max(sqrt(Bw(1)**2 + Bw(2)**2 + Bw(3)**2), 1.0e-30_dp)
527 288155 call perp_ref(bw_hat, ew) ! arbitrary unit vector perpendicular to b
528 288155 bhat = rotz(bw_hat, ca, sa) ! wedge -> global frame
529 288155 eperp = rotz(ew, ca, sa)
530 288155 status = FO_OK
531 end subroutine locate
532
533 ! Seed from a guiding-centre start record u0=(s, theta_B, phi_B) with parallel
534 ! speed vpar0 and perpendicular speed vperp0. Place the particle a Larmor vector
535 ! off the guiding centre in Cartesian (regular through the axis) and seed
536 ! v = vpar0 bhat + vperp0 e_perp with e_perp the same gyrophase reference the
537 ! position offset uses.
538
6/6
✓ Branch 0 taken 48 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 48 times.
✓ Branch 3 taken 16 times.
✓ Branch 4 taken 48 times.
✓ Branch 5 taken 16 times.
 160 subroutine fo_init(st, x0_boozer, vpar0, vperp0, mu_in, mass, &
539 charge, dt, ro0_in, pabs)
540 type(fo_state_t), intent(out) :: st
541 real(dp), intent(in) :: x0_boozer(3), vpar0, vperp0, mu_in, mass, charge, &
542 dt, ro0_in, pabs
543 real(dp) :: u_gc(3)
544
545
2/2
✓ Branch 0 taken 48 times.
✓ Branch 1 taken 16 times.
 64 u_gc = [sqrt(max(x0_boozer(1), 0.0_dp)), x0_boozer(2), x0_boozer(3)]
546 call fo_init_logical(st, u_gc, vpar0, vperp0, mu_in, mass, charge, dt, &
547 16 ro0_in, pabs, .false.)
548 16 end subroutine fo_init
549
550
6/6
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 1 time.
 12 subroutine fo_init_reference(st, u_gc, vpar0, vperp0, mu_in, mass, &
551 charge, dt, ro0_in, pabs)
552 type(fo_state_t), intent(out) :: st
553 real(dp), intent(in) :: u_gc(3), vpar0, vperp0, mu_in, mass, charge, dt, &
554 ro0_in, pabs
555
556 call fo_init_logical(st, u_gc, vpar0, vperp0, mu_in, mass, charge, dt, &
557 1 ro0_in, pabs, .true.)
558 1 end subroutine fo_init_reference
559
560
6/6
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 17 times.
✓ Branch 2 taken 51 times.
✓ Branch 3 taken 17 times.
✓ Branch 4 taken 51 times.
✓ Branch 5 taken 17 times.
 170 subroutine fo_init_logical(st, u_gc, vpar0, vperp0, mu_in, mass, charge, &
561 dt, ro0_in, pabs, reference_field)
562 type(fo_state_t), intent(out) :: st
563 real(dp), intent(in) :: u_gc(3), vpar0, vperp0, mu_in, mass, charge, dt, &
564 ro0_in, pabs
565 logical, intent(in) :: reference_field
566 real(dp) :: xyz_gc(3), u_p(3), x_p(3), qc
567 real(dp) :: bhat(3), eperp(3), Bmod
568 integer :: status
569
570 17 st%mass = mass; st%charge = charge; st%dt = dt; st%ro0 = ro0_in
571 17 st%mu = mu_in; st%pabs = pabs
572 17 st%reference_field = reference_field
573
574 17 call ref_coords%evaluate_cart(u_gc, xyz_gc)
575 17 qc = charge/ro0_in
576
577 17 x_p = xyz_gc
578 17 u_p = u_gc
579
1/2
✓ Branch 0 taken 17 times.
✗ Branch 1 not taken.
 17 if (vperp0 > 0.0_dp) then
580 ! Larmor offset off the guiding centre. If the offset point falls outside the
581 ! chart (a near-edge marker whose gyro-circle pokes past s=1) or fails to
582 ! locate, fall back to seeding at the guiding centre: the offset is O(rho_L),
583 ! and a genuine edge orbit is then lost during integration, not at init. Never
584 ! abort -- this runs per particle inside the OpenMP loop.
585 call gc_to_particle(xyz_gc, u_gc, vperp0, mass, qc, x_p, u_p, status, &
586 17 reference_field)
587
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17 times.
 17 if (status /= FO_OK) then
588 ✗ x_p = xyz_gc
589 ✗ u_p = u_gc
590 end if
591 end if
592
593
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 17 times.
 68 st%x = x_p
594
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 17 times.
 68 st%u = u_p
595 17 call locate(x_p, u_p, u_p, bhat, eperp, Bmod, status, reference_field)
596
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17 times.
 17 if (status /= FO_OK) then
597 ! Cannot even seed the frame at the guiding centre: leave v=0 so the first
598 ! orbit step reports a locate fault (counted confined), never a crash.
599 ✗ st%v = 0.0_dp
600 ✗ return
601 end if
602
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 17 times.
 68 st%u = u_p
603
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 17 times.
 68 st%v = vpar0*bhat
604
3/4
✓ Branch 0 taken 17 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 51 times.
✓ Branch 3 taken 17 times.
 68 if (vperp0 > 0.0_dp) st%v = st%v + vperp0*eperp
605 end subroutine fo_init_logical
606
607 ! Cartesian guiding centre x_gc -> particle position a Larmor vector off it,
608 ! solved by the fixed point x_p with cart(x_p) - rho(x_p) = x_gc, rho the Larmor
609 ! vector built from the perpendicular speed at x_p (same gyrophase reference as
610 ! the velocity seed), so the seed offset and the GC reconstruction are inverses.
611 17 subroutine gc_to_particle(xyz_gc, u_gc, vperp0, mass, qc, x_p, u_p, status, &
612 reference_field)
613 real(dp), intent(in) :: xyz_gc(3), u_gc(3), vperp0, mass, qc
614 real(dp), intent(out) :: x_p(3), u_p(3)
615 integer, intent(out) :: status
616 logical, intent(in) :: reference_field
617 integer, parameter :: maxfp = 50
618 real(dp), parameter :: tol = 1.0e-10_dp
619 real(dp) :: bhat(3), eperp(3), Bmod
620 real(dp) :: rho_l(3), xnew(3)
621 integer :: it
622
623 17 x_p = xyz_gc
624
1/2
✓ Branch 0 taken 109 times.
✗ Branch 1 not taken.
 109 do it = 1, maxfp
625 109 call locate(x_p, u_gc, u_p, bhat, eperp, Bmod, status, reference_field)
626
1/2
✓ Branch 0 taken 109 times.
✗ Branch 1 not taken.
 126 if (status /= FO_OK) return
627 ! rho = (m/(qc|B|)) bhat x v_perp, v_perp = vperp0 eperp (Cartesian).
628
4/4
✓ Branch 0 taken 327 times.
✓ Branch 1 taken 109 times.
✓ Branch 3 taken 327 times.
✓ Branch 4 taken 109 times.
 763 rho_l = (mass/(qc*Bmod))*cross(bhat, vperp0*eperp)
629
2/2
✓ Branch 0 taken 327 times.
✓ Branch 1 taken 109 times.
 436 xnew = xyz_gc + rho_l
630
8/10
✗ Branch 0 not taken.
✓ Branch 1 taken 109 times.
✓ Branch 2 taken 109 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 327 times.
✓ Branch 5 taken 109 times.
✓ Branch 6 taken 109 times.
✓ Branch 7 taken 218 times.
✓ Branch 8 taken 17 times.
✓ Branch 9 taken 92 times.
 545 if (maxval(abs(xnew - x_p)) < tol) then
631 17 x_p = xnew
632 17 call locate(x_p, u_gc, u_p, bhat, eperp, Bmod, status, reference_field)
633 17 return
634 end if
635 201 x_p = xnew
636 end do
637 ✗ status = FO_OK
638 end subroutine gc_to_particle
639
640 60003 subroutine fo_step(st, status)
641 type(fo_state_t), intent(inout) :: st
642 integer, intent(out) :: status
643 real(dp) :: x(3), v(3), Bvec(3), Bmod, gradB(3), u(3)
644 real(dp) :: tvec(3), svec(3), vp(3), tmag2, qcm
645
646
2/2
✓ Branch 0 taken 180009 times.
✓ Branch 1 taken 60003 times.
 240012 x = st%x
647
2/2
✓ Branch 0 taken 180009 times.
✓ Branch 1 taken 60003 times.
 240012 v = st%v
648 60003 qcm = st%charge/(c*st%ro0*st%mass) ! rotation: dv/dt = qcm v x B
649
650
2/2
✓ Branch 0 taken 180009 times.
✓ Branch 1 taken 60003 times.
 240012 x = x + 0.5_dp*st%dt*v
651 call cart_field(x, st%u, Bvec, Bmod, gradB, u, status, &
652 60003 st%reference_field)
653
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 60003 times.
 60003 if (status /= FO_OK) return ! unresolved: leave st at the last resolved state
654
2/2
✓ Branch 0 taken 180009 times.
✓ Branch 1 taken 60003 times.
 240012 st%u = u
655
656 ! exact magnetic rotation (constant B over the step).
657
2/2
✓ Branch 0 taken 180009 times.
✓ Branch 1 taken 60003 times.
 240012 tvec = qcm*Bvec*0.5_dp*st%dt
658 60003 tmag2 = tvec(1)**2 + tvec(2)**2 + tvec(3)**2
659
2/2
✓ Branch 0 taken 180009 times.
✓ Branch 1 taken 60003 times.
 240012 svec = 2.0_dp*tvec/(1.0_dp + tmag2)
660
2/2
✓ Branch 1 taken 180009 times.
✓ Branch 2 taken 60003 times.
 300015 vp = v + cross(v, tvec)
661
2/2
✓ Branch 1 taken 180009 times.
✓ Branch 2 taken 60003 times.
 240012 v = v + cross(vp, svec)
662
663
2/2
✓ Branch 0 taken 180009 times.
✓ Branch 1 taken 60003 times.
 240012 x = x + 0.5_dp*st%dt*v
664
665
2/2
✓ Branch 0 taken 180009 times.
✓ Branch 1 taken 60003 times.
 240012 st%x = x
666
2/2
✓ Branch 0 taken 180009 times.
✓ Branch 1 taken 60003 times.
 240012 st%v = v
667 60003 status = FO_OK
668 end subroutine fo_step
669
670 84015 function fo_energy(st) result(energy)
671 type(fo_state_t), intent(in) :: st
672 real(dp) :: energy
673 84015 energy = 0.5_dp*st%mass*(st%v(1)**2 + st%v(2)**2 + st%v(3)**2)
674 84015 end function fo_energy
675
676 ! Guiding-centre magnetic moment mu = m v_perp^2/(2|B_gc|): evaluate at the
677 ! Larmor-corrected guiding centre, not the raw particle point, so the gyro
678 ! ripple O(rho/L) is removed and mu is conserved to O((rho/L)^2). Diagnostic only.
679 60003 function fo_mu(st) result(mu)
680 type(fo_state_t), intent(in) :: st
681 real(dp) :: mu, s, th, ph, vpar, Bgc
682 integer :: status
683 60003 call fo_to_gc(st, s, th, ph, vpar, status, Bmod_gc=Bgc)
684
1/2
✓ Branch 0 taken 60003 times.
✗ Branch 1 not taken.
 60003 if (status /= FO_OK) then
685 60003 mu = 0.0_dp; return
686 end if
687 mu = 0.5_dp*st%mass*max(st%v(1)**2 + st%v(2)**2 + st%v(3)**2 - vpar**2, &
688 60003 0.0_dp)/max(Bgc, 1.0e-30_dp)
689 60003 end function fo_mu
690
691 ! Guiding-centre reduction for output: remove the Larmor vector in Cartesian and
692 ! report the centre in (s, theta_B, phi_B) with the parallel speed at the centre.
693 ! status: FO_OK / FO_LOSS / FO_LOCATE_FAIL.
694 144005 subroutine fo_to_gc(st, s, th, ph, vpar, status, Bmod_gc)
695 type(fo_state_t), intent(in) :: st
696 real(dp), intent(out) :: s, th, ph, vpar
697 integer, intent(out) :: status
698 real(dp), intent(out), optional :: Bmod_gc
699 real(dp) :: u_gc(3), Bmod_local
700
701 144005 call fo_to_reference_gc(st, u_gc, vpar, status, Bmod_local)
702
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 144005 times.
 144005 if (status /= FO_OK .and. status /= FO_LOSS) then
703 ✗ s = 0.0_dp
704 ✗ th = 0.0_dp
705 ✗ ph = 0.0_dp
706 ✗ if (present(Bmod_gc)) Bmod_gc = 0.0_dp
707 ✗ return
708 end if
709 144005 s = u_gc(1)**2
710 144005 th = u_gc(2)
711 144005 ph = u_gc(3)
712
2/2
✓ Branch 0 taken 60003 times.
✓ Branch 1 taken 84002 times.
 144005 if (present(Bmod_gc)) Bmod_gc = Bmod_local
713 end subroutine fo_to_gc
714
715 288012 subroutine fo_to_reference_gc(st, u_gc, vpar, status, Bmod_gc)
716 type(fo_state_t), intent(in) :: st
717 real(dp), intent(out) :: u_gc(3), vpar
718 integer, intent(out) :: status
719 real(dp), intent(out), optional :: Bmod_gc
720 real(dp) :: u_p(3), x_gc(3), qc
721 real(dp) :: bhat(3), eperp(3), Bmod
722 real(dp) :: vpar_p, vperp_cart(3), rho_l(3)
723
724 144006 u_gc = 0.0_dp
725 144006 vpar = 0.0_dp
726
2/2
✓ Branch 0 taken 144005 times.
✓ Branch 1 taken 1 time.
 144006 if (present(Bmod_gc)) Bmod_gc = 0.0_dp
727
728 call locate(st%x, st%u, u_p, bhat, eperp, Bmod, status, &
729 144006 st%reference_field)
730
1/2
✓ Branch 0 taken 144006 times.
✗ Branch 1 not taken.
 144006 if (status /= FO_OK) return
731 144006 qc = st%charge/st%ro0
732
733 ! Larmor vector from the particle's perpendicular velocity at x (Cartesian):
734 ! rho = (m/(qc|B|)) bhat x v_perp; x_gc = x - rho.
735 144006 vpar_p = st%v(1)*bhat(1) + st%v(2)*bhat(2) + st%v(3)*bhat(3)
736
2/2
✓ Branch 0 taken 432018 times.
✓ Branch 1 taken 144006 times.
 576024 vperp_cart = st%v - vpar_p*bhat
737
2/2
✓ Branch 1 taken 432018 times.
✓ Branch 2 taken 144006 times.
 576024 rho_l = (st%mass/(qc*Bmod))*cross(bhat, vperp_cart)
738
2/2
✓ Branch 0 taken 432018 times.
✓ Branch 1 taken 144006 times.
 576024 x_gc = st%x - rho_l
739
740 call locate(x_gc, u_p, u_gc, bhat, eperp, Bmod, status, &
741 144006 st%reference_field)
742
1/2
✓ Branch 0 taken 144006 times.
✗ Branch 1 not taken.
 144006 if (status /= FO_OK) return
743 144006 vpar = st%v(1)*bhat(1) + st%v(2)*bhat(2) + st%v(3)*bhat(3)
744
2/2
✓ Branch 0 taken 144005 times.
✓ Branch 1 taken 1 time.
 144006 if (present(Bmod_gc)) Bmod_gc = Bmod
745 ! Confinement loss: the Larmor-corrected guiding centre crosses the last closed
746 ! surface (u_gc >= 1). The GC must be locatable, so the loss keys on u_gc rather
747 ! than the particle (a particle gyro-excursed past s=1 is off-chart and cannot be
748 ! inverted). But u_gc is a second, cold-guess locate of x_gc and carries the
749 ! residual field-period-seam noise, which occasionally returns rho >= 1 for a
750 ! particle that is in fact at mid-radius. Reject that with the robust warm-started
751 ! particle locate u_p: a real loss has the particle within ~a Larmor radius of the
752 ! edge (|x_gc - x| = |rho_l| is a Larmor radius), so u_gc >= 1 while u_p is well
753 ! inside is a reconstruction glitch, not a loss. The field, integrator and energy
754 ! match the ASCOT full-orbit reference, so the loss detector must be this clean.
755
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 144006 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 144006 if (u_gc(1) >= 1.0_dp .and. u_p(1) >= 1.0_dp - GC_PARTICLE_GAP) status = FO_LOSS
756 end subroutine fo_to_reference_gc
757
758 ! Pseudo-Cartesian near-axis chart w=(X,Y,phi)=(rho cos th, rho sin th, phi).
759 ! The chartmap polar chart (rho,theta) is singular at the magnetic axis: the
760 ! covariant basis column Jc(:,2)=dx/dtheta ~ rho vanishes, so det(Jc)->0 and both
761 ! the inverse Newton (ill-conditioned in theta) and the field assembly degrade.
762 ! The (X,Y) basis stays regular through rho=0 (Pfefferle et al.,
763 ! arXiv:1412.5464; libneo flux_pseudocartesian). Returns the regular chart
764 ! Jacobian Jw(a,i)=dx_a/dw_i and the trig used to map field components.
765 4413788 subroutine pseudocart_basis(u, Jc, Jw, cth, sth, rho)
766 real(dp), intent(in) :: u(3), Jc(3,3)
767 real(dp), intent(out) :: Jw(3,3), cth, sth, rho
768 integer :: a
769 4413788 rho = max(u(1), 1.0e-30_dp)
770 4413788 cth = cos(u(2)); sth = sin(u(2))
771
2/2
✓ Branch 0 taken 13241364 times.
✓ Branch 1 taken 4413788 times.
 17655152 do a = 1, 3
772 13241364 Jw(a,1) = Jc(a,1)*cth - Jc(a,2)*(sth/rho) ! e_X = dx/dX
773 13241364 Jw(a,2) = Jc(a,1)*sth + Jc(a,2)*(cth/rho) ! e_Y = dx/dY
774 17655152 Jw(a,3) = Jc(a,3) ! e_phi
775 end do
776 4413788 end subroutine pseudocart_basis
777
778 ! Pseudo-Cartesian w=(X,Y,phi) -> logical u=(rho,theta,phi). rho>=0 and the
779 ! atan2 branch make the axis an ordinary point (no reflect hack on the inverse).
780 2992356 pure subroutine w_to_u(w, u)
781 real(dp), intent(in) :: w(3)
782 real(dp), intent(out) :: u(3)
783 2992356 u(1) = sqrt(w(1)**2 + w(2)**2)
784 2992356 u(2) = atan2(w(2), w(1))
785 2992356 u(3) = w(3)
786 2992356 end subroutine w_to_u
787
788 ! An arbitrary unit vector perpendicular to b, regular everywhere. The gyrophase
789 ! reference is gauge: only b and |v_perp| are physical, and the guiding-centre
790 ! reduction recovers v_perp from v directly, not from this choice. Gram-Schmidt
791 ! off the least-aligned axis so the subtraction never cancels.
792 288155 pure subroutine perp_ref(b, e)
793 real(dp), intent(in) :: b(3)
794 real(dp), intent(out) :: e(3)
795 real(dp) :: a(3), d, n
796
1/2
✓ Branch 0 taken 288155 times.
✗ Branch 1 not taken.
 288155 if (abs(b(3)) < 0.9_dp) then
797 288155 a = [0.0_dp, 0.0_dp, 1.0_dp]
798 else
799 ✗ a = [1.0_dp, 0.0_dp, 0.0_dp]
800 end if
801 288155 d = a(1)*b(1) + a(2)*b(2) + a(3)*b(3)
802
2/2
✓ Branch 0 taken 864465 times.
✓ Branch 1 taken 288155 times.
 1152620 e = a - d*b
803 288155 n = sqrt(e(1)**2 + e(2)**2 + e(3)**2)
804
2/2
✓ Branch 0 taken 864465 times.
✓ Branch 1 taken 288155 times.
 1152620 e = e/max(n, 1.0e-30_dp)
805 288155 end subroutine perp_ref
806
807 ! A chart Jacobian is usable when its determinant is well above the round-off
808 ! floor relative to its size (the chartmap is singular at the magnetic axis,
809 ! rho->0). Rejecting near-singular Jc keeps the field push and the inversion off
810 ! the axis singularity instead of producing Inf/NaN.
811 4413788 pure logical function jacobian_ok(Jc)
812 real(dp), intent(in) :: Jc(3,3)
813 real(dp) :: det, scale
814 det = Jc(1,1)*(Jc(2,2)*Jc(3,3) - Jc(2,3)*Jc(3,2)) &
815 - Jc(1,2)*(Jc(2,1)*Jc(3,3) - Jc(2,3)*Jc(3,1)) &
816 4413788 + Jc(1,3)*(Jc(2,1)*Jc(3,2) - Jc(2,2)*Jc(3,1))
817
4/4
✓ Branch 0 taken 13241364 times.
✓ Branch 1 taken 4413788 times.
✓ Branch 2 taken 39724092 times.
✓ Branch 3 taken 13241364 times.
 57379244 scale = sqrt(sum(Jc**2))
818
2/4
✓ Branch 0 taken 4413788 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4413788 times.
 4413788 jacobian_ok = (det == det) .and. abs(det) > 1.0e-8_dp*max(scale, 1.0e-30_dp)**3
819 4413788 end function jacobian_ok
820
821
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1337395 times.
 1337395 pure function cross(a, b) result(cr)
822 real(dp), intent(in) :: a(3), b(3)
823 real(dp) :: cr(3)
824 1337395 cr(1) = a(2)*b(3) - a(3)*b(2)
825 1337395 cr(2) = a(3)*b(1) - a(1)*b(3)
826 1337395 cr(3) = a(1)*b(2) - a(2)*b(1)
827 1337395 end function cross
828
829 4413788 pure subroutine inv3(A, Ainv)
830 real(dp), intent(in) :: A(3,3)
831 real(dp), intent(out) :: Ainv(3,3)
832 real(dp) :: det
833 det = A(1,1)*(A(2,2)*A(3,3) - A(2,3)*A(3,2)) &
834 - A(1,2)*(A(2,1)*A(3,3) - A(2,3)*A(3,1)) &
835 4413788 + A(1,3)*(A(2,1)*A(3,2) - A(2,2)*A(3,1))
836 4413788 Ainv(1,1) = (A(2,2)*A(3,3) - A(2,3)*A(3,2))/det
837 4413788 Ainv(1,2) = (A(1,3)*A(3,2) - A(1,2)*A(3,3))/det
838 4413788 Ainv(1,3) = (A(1,2)*A(2,3) - A(1,3)*A(2,2))/det
839 4413788 Ainv(2,1) = (A(2,3)*A(3,1) - A(2,1)*A(3,3))/det
840 4413788 Ainv(2,2) = (A(1,1)*A(3,3) - A(1,3)*A(3,1))/det
841 4413788 Ainv(2,3) = (A(1,3)*A(2,1) - A(1,1)*A(2,3))/det
842 4413788 Ainv(3,1) = (A(2,1)*A(3,2) - A(2,2)*A(3,1))/det
843 4413788 Ainv(3,2) = (A(1,2)*A(3,1) - A(1,1)*A(3,2))/det
844 4413788 Ainv(3,3) = (A(1,1)*A(2,2) - A(1,2)*A(2,1))/det
845 4413788 end subroutine inv3
846
847
1/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1421432 times.
 1421432 end module orbit_fo_boris
848