text:
\section{Introduction}

Global kinetic computations of quasi-steady plasma parameters in 3D toroidal fusion devices utilize the evaluation of the distribution function and/or its moments by direct modeling of particle orbits. This includes Monte Carlo transport simulations in given external fields~\citep{boozer_monte_1981,lotz_monte_1988, heikkinen_particle_2001, wakasa_study_2008, tribaldos_monte_2001, isaev_venusf_2006,allmaier_variance_2008, satake_neoclassical_2011, pfefferle_venus-levis_2014} as well as self-consistent turbulence models with particle codes.~\citep{jolliet_global_2007,mishchenko_global_2008,hager_gyrokinetic_2016,ku_fast_2018} Kinetic modeling of 3D plasma equilibria~\cite{albert_kinetic_2016} or edge plasmas~\cite{hager_gyrokinetic_2019} puts specific requirements on solving the guiding-center equations.~\citep{morozov_structure_1966,boozer_guiding_1980,littlejohn_variational_1983,cary_hamiltonian_2009} Namely, particle orbit integration should be computationally efficient, tolerant to statistical noise in the electromagnetic field and efficient in scoring statistical data from the orbits.
Geometric integrators~\citep{hairer_geometric_2006,morrison_structure_2017} address these targets by relaxing the requirement to the accuracy of guiding-center orbits while preserving physically \bluex{correct}\redx{consistent} long time orbit dynamics. In this context the word geometric refers to the preservation of the geometric structure of phase-space being a symplectic manifold.
The \blue{most well-known}\red{best-known} class of geometric integrators are symplectic integrators that rely on canonical coordinates in phase-space. Symplectic integrators are not directly applicable to guiding-center dynamics that are formulated in non-canonical coordinates. One way to circumvent this problem is the use of a (usually implicit) transformation to canonical coordinates.~\citep{zhang_canonicalization_2014,zhu_symmetric_2016,albert_symplectic_2020} Up to now such approaches rely either on magnetic flux coordinates or require a more expensive transformation of phase-space coordinates in the general case. The other
\red{well-}known alternative are variational integrators.~\citep{qin_variational_2008,qin_variational_2009} Such integrators do not assume canonical coordinates and include symplectic integrators as a special case. Stability problems of variational integrators arise for guiding-center orbit computations due to the degeneracy of their phase-space Lagrangian. This issue has only recently been resolved.~\citep{burby_toroidal_2017,kraus_projected_2017,ellison_degenerate_2018} The resulting integrators are fully implicit and/or require an augmented set of dynamical variables, and their competitiveness for practical
applications is still a topic open to investigation. \red{Yet another alternative
could arise from very recent results on slow-manifold methods~\cite{xiao_slow_2020}
that introduce a different construction of the guiding-center equations well suited for geometric integration.}

The method presented here is geometric in both\blue{,} a structure-preserving and a more literal sense, as it considers
intersections of orbits with a spatial mesh. \red{The underlying formulation and discretization
of fields exactly preserves the non-canonical Hamiltonian structure of the equations.
\redx{In contrast to usual geometric integrators, this is achieved by solving
the exact motion in simplified fields, i.e. the Hamiltonian flow in the
original fields is approximated by the exact Hamiltonian flow in piecewise
linear fields.}
The final algorithm is only quasi-geometric, as it relies on a series expansion in the orbit parameter.
\redx{Here the term quasi-geometric means that the error can be brought below
any given (computer) accuracy by a sufficiently high order.}
\bluex{but}\redx{Expansions to order 3 and 4 are} \bluex{is} shown to conserve invariants extremely well for at least $10^6$ toroidal turns in numerical experiments for typical fusion devices.}
\blue{It}\red{The approach} has been introduced in Ref.~\onlinecite{eder_three-dimensional_2019}
as a generalization of the 2D geometric integrator of Ref.~\onlinecite{kasilov_geometric_2016} for general 3D toroidal fields.
This approach has two features useful for application in global equilibrium and transport simulations: straightforward computation of
spatial distributions of macroscopic parameters and robustness
in the presence of noise in field data.  Moreover, the new method preserves total energy, magnetic moment and phase space volume. \red{In its present version the method is designed for (quasi-)static electromagnetic fields.}

The integration procedure is based on a special 3D discretization of space resulting in \blue{piecewise}\red{locally} linear guiding-center
equations while retaining the symplectic property of the original set.
\red{
Formally such an interpolation could be denoted by Whitney forms or finite element exterior calculus in a way similar
to existing work on charged particle orbits~\cite{squire_geometric_2012,xiao_explicit_2015,kraus_gempic_2017} and results in a divergence-free magnetic field.
}
Within this discretization, vectors and
scalars characterizing the electromagnetic field are represented by continuous piecewise linear functions
which reduces the cost of spatial interpolation as compared to high-order interpolation with continuous derivatives
(e.g. with help of 3D-splines) required for usual integration. In return, the integration procedure requires
accurate tracing of intersections of the orbit with spatial cell boundaries where the coefficients of the linear
guiding-center equation set
are discontinuous. However, tracing of the boundaries is also required in Monte Carlo procedures for
the computation of the spatial distribution of various velocity space moments of the distribution function which are
computed as path integrals over the test particle dwell time within spatial cells.

In the original implementation of Ref.~\onlinecite{eder_three-dimensional_2019}, the linear guiding-center equation
set was solved numerically by using a usual Runge-Kutta (RK) integrator.
In this case, tracing of orbit intersections with spatial cell boundaries requires a few Newton iterations for the computation
of the dwell time within the cell. Since the integration error of this set scales with the third power of
the Larmor radius, accurate results can still be obtained by a single RK integration step over the dwell time.
For the same reason, sufficient accuracy can be achieved also by using the polynomial expansion of the solution over \red{an orbit parameter which allows to compute the dwell time} and integrals of velocity powers analytically.

In \red{the }case of magnetic fields with spatial symmetry (e.g., tokamaks with toroidal symmetry) the \red{quasi-}geometric integration
accurately preserves the respective (toroidal) canonical momentum. Thus, the property of such systems to ideally confine
the orbits is retained. In \red{the }absence of spatial symmetry, the parallel adiabatic invariant, which is an approximate conserved
quantity in stellarators, is well preserved, since the \red{quasi-}geometric integration does not lead to a significant error accumulation.

It should be noted that the piecewise linear approximation of the guiding-center equations represents field lines as
polygonal chains in coordinate space. Such a representation may introduce artificial chaos in case of 3D fields
when using non-aligned coordinate systems.
In our earlier publication~\citep{eder_three-dimensional_2019}, some field line diffusion has been observed
in a perturbed tokamak field which seemed to be in agreement with the quasilinear estimate assuming a continuous safety
factor profile of the unperturbed field (which is, actually, not the case), and, therefore, this diffusion has
been attributed to the linearization.
However, further detailed studies and resulting improvements of the algorithm revealed that the level of
artificial diffusion observed earlier is not connected with the intrinsic deficiency of the method. In the improved
algorithm this diffusion is actually much smaller so that it provides a negligible correction to the existing neoclassical
transport even for \blue{rather}\red{the relatively} strongly perturbed non-aligned 3D fields and a coarse mesh.

The above mentioned characteristics of the method allow, in particular, its effective application to the computation
of neoclassical transport coefficients using the Monte Carlo method. For demonstration, the mono-energetic neoclasscial
diffusion coefficient is evaluated here for a quasi-isodynamic reactor-scale stellarator
field~\cite{drevlak_quasi-isodynamic_2014}. The results and performance of the new method are compared to usual (RK)
orbit integration methods.

The remainder of this article is organized as follows. In section~\ref{sec:derivation_geometric_integrator} the spatial discretization procedure
is introduced and the numerical solution of the resulting piecewise linear guiding-center equations is described.
In section~\ref{sec:orbits} single particle orbits obtained with the \red{quasi-}geometric integration method and respective invariants of motion for axisymmetric and non-axisymmetric geometries are analyzed in detail. Furthermore, the introduced artificial chaos is studied.
The \blue{method's }application \red{of the method }to the evaluation of transport coefficients is presented in section~\ref{sec:transport}, where the computational orbit-integration performance is benchmarked as well.
Finally, the conclusion from the present study and further outlook is given in section~\ref{sec:summary}.
\red{Details on the Hamiltonian structure of the underlying locally linear equations are discussed in Appendix~\ref{sec:appendix_hamiltonian} and a useful feature for scoring orbit statistics is presented in Appendix~\ref{sec:appendix_integrals}.}

\section{Derivation of the orbit integration method and numerical solution}
\label{sec:derivation_geometric_integrator}

\subsection{Locally linear Hamiltonian guiding-center equations}
\label{sec:piecewise}

\blue{As a starting point, the equations of guiding-center motion\citep{boozer_guiding_1980} in general curvilinear coordinates $x^i$
are}
\red{As a starting point, the equations of guiding-center motion, Eq.~(11) of Ref.~\onlinecite{littlejohn_variational_1983} equivalent to Eq.~(A7) of Ref.~\onlinecite{boozer_guiding_1980}, are expressed in general curvilinear coordinates $x^i$,}

\be{eqm_curv}
\dot x^i = \frac{v_\parallel \varepsilon^{ijk}}{\sqrt{g}\, B_\parallel^\ast}\difp{A^\ast_k}{x^j},
\qquad A^\ast_k = A_k + \frac{v_\parallel}{\omega_c} B_k.
\end{equation}

Here, $A_k$, $B_k$, $\omega_c$, $\Phi$ and $\sqrt{g}$ are the covariant components of the vector potential and the magnetic field, cyclotron frequency, electrostatic potential and the metric
determinant, respectively, and $\sqrt{g}\, B_\parallel^\ast = \varepsilon^{ijk} (B_i/B) \partial A^\ast_k/\partial x^j$.
Charge $e_\alpha$ and mass $m_\alpha$ of the considered species $\alpha$ enter
$\omega_c = e_\alpha B / (m_\alpha c)$ together with the magnetic field modulus
$B = \sqrt{B_k B^k}$ and the speed of light $c$.
The equations of motion are considered with the invariants $w=m_\alpha v^2/2+e_\alpha\Phi$ and $J_\perp=m_\alpha v_\perp^2/(2\omega_c)$ being total energy and perpendicular adiabatic invariant, respectively, and used as independent phase space variables. \red{The latter variable is related to the magnetic moment $\mu$ by a constant factor, $J_\perp=\mu m_\alpha c / e_\alpha$.}
The parallel velocity $v_\parallel$ in~\eq{eqm_curv} is not an independent variable but a known function of coordinates,
\be{defU}
v_\parallel^2=2U, \qquad
U=U({\bf x})=\frac{1}{m}\left(w-J_\perp\omega_c({\bf x})-e_\alpha\Phi({\bf x})\right).
\end{equation}
\redx{Due to the fact that the actual guiding-center orbits do not depend on the choice of their representing variables in phase space,
the parallel velocity $v_\parallel$ can also be treated as an independent variable. (Choosing $v_\parallel = v_\parallel \left({\bf x}, w, J_\perp \right)$ or $w = w \left({\bf x}, v_\parallel, J_\perp \right)$ is equivalent as long as Eq.~\eq{defU} is kept exact.)}
\bluex{Treating $v_\parallel$ as an independent variable now, i.e. replacing the first expression of~\eq{defU} by the differential equation}
\redx{By replacing the first expression of~\eq{defU} with the differential equation}
\be{}
\dot v_\parallel = \frac{1}{v_\parallel}\dot x^i \frac{\partial U} {\partial x^i},
\end{equation}
the set (\ref{eqm_curv}) turns into
\bea{setofgceqs}
B_\parallel^\ast\sqrt{g}\,\dot x^i = \frac{{\rm d} x^i}{{\rm d} \tau}
&=&
\varepsilon^{ijk}\left(
v_\parallel \difp{A_k}{x^j}+2U\difp{}{x^j}\frac{B_k}{\omega_c}+\frac{B_k}{\omega_c}\difp{U}{x^j}
\right),
\nonumber\\
B_\parallel^\ast\sqrt{g}\,\dot v_\parallel = \frac{{\rm d} v_\parallel}{{\rm d} \tau}
&=&
\varepsilon^{ijk}\difp{U}{x^i}
\left(
\difp{A_k}{x^j}+v_\parallel \difp{}{x^j}\frac{B_k}{\omega_c}
\right).
\end{eqnarray}
\redx{Here, $v_\parallel$ should be treated as an independent variable only at its explicit occurrences,
while the quantity $U({\bf x})$ should only be treated as a function of coordinates as defined by
the second expression of~\eq{defU}.}
Note that the invariants of motion remain in~\eq{setofgceqs} as parameters entering
the function of coordinates $U$\redx{.} \bluex{defined by the second expression of~\eq{defU}.}
\redx{As long as the contra-variant components of the phase space velocity defined by right hand sides
of set~\eq{setofgceqs} are used for obtaining the exact solution of this set, both
representations~\eq{defU} lead to the same result. However, if the derivatives of these components
are computed for the Jacobian as, e.g., in Appendix~\ref{sec:appendix_hamiltonian_Liouv}
the quantity $U$ should be treated as a function of parallel velocity defined by the first of~\eq{defU}.}
In~\eq{setofgceqs} the time variable is replaced
by an orbit parameter $\tau$ related to time by ${\rm d} t = B_\parallel^\ast\sqrt{g}{\rm d} \tau$,
and the time evolution is obtained implicitly from the integral $t(\tau)$.

The special form~\eq{setofgceqs} allows to reduce computational effort and noise sensitivity by independently approximating
the field quantities $A_k$, $B_k/\omega_c$, $\omega_c$ and $\Phi$ by continuous piecewise linear functions.
Thus, curvilinear coordinate
space is split into tetrahedral cells with exact field values on the cell's vertices.
Fig. \ref{fig:curvilinear_grid} depicts such a \redx{real space illustration of a }curvilinear field-aligned grid for the plasma core of a tokamak. \bluex{, where splitting along the symmetry direction is achieved with stackable hexahedra, each consisting of six tetrahedral cells.}
\redx{These tetrahedral cells must be specially oriented (explained below in section~\ref{sec:orbits_2d}) in order to preserve an invariant of motion in the case of axisymmetry.}

As a result \redx{of this piecewise field linearization}, in each cell, the equations of motion (\ref{setofgceqs})
turn into a set of four linear ODEs with constant coefficients
\be{standeqset}
\frac{{\rm d} z^i}{{\rm d} \tau} = a^i_l z^l + b^i,
\end{equation}
in phase-space variables $z^i=x^i$ for $i=1,2,3$ and $z^4=v_\parallel$. The matrix elements are
\bea{amatdef}
a^i_l &=& \varepsilon^{ijk}\left(
2\difp{U}{x^l}\difp{}{x^j}\frac{B_k}{\omega_c}+\difp{U}{x^j}\difp{}{x^l}\frac{B_k}{\omega_c}
\right)
\qquad\mbox{for}\qquad 1\le i,l \le 3,
\nonumber \\
a^i_4 &=& \varepsilon^{ijk} \difp{A_k}{x^j}
\qquad\mbox{for}\qquad 1\le i \le 3,
\nonumber \\
a^4_l &=& 0
\qquad\mbox{for}\qquad 1\le l \le 3,
\nonumber \\
a^4_4 &=& \varepsilon^{ijk}\difp{U}{x^i}
\difp{}{x^j}\frac{B_k}{\omega_c},
\end{eqnarray}
and components of vector $b^i$ are
\bea{bvecdef}
b^i &=& \varepsilon^{ijk}\left(
2U_0\difp{}{x^j}\frac{B_k}{\omega_c}+\left(\frac{B_k}{\omega_c}\right)_0\difp{U}{x^j}
\right)
\qquad\mbox{for}\qquad 1\le i \le 3,
\nonumber \\
b^4 &=& \varepsilon^{ijk}\difp{U}{x^i}
\difp{A_k}{x^j},
\end{eqnarray}
where quantities with zero mean the value at the origin of the coordinates \red{shifted in each
tetrahedral cell to one of the cell's vertices,}

\be{valorig}
U=U_0+x^i\difp{U}{x^i}, \qquad \frac{B_k}{\omega_c} = \left(\frac{B_k}{\omega_c}\right)_0
+x^i \difp{}{x^i}\frac{B_k}{\omega_c}.
\end{equation}

Since the piecewise constant coefficients of set~\eq{standeqset} are discontinuous at
the cell boundaries, orbit intersections with tetrahedra faces
must be computed exactly when integrating particle trajectories.

In fact, a linear approximation of field quantities which locally breaks the physical connection between them
does not destroy the Hamiltonian nature of the original set~\eq{eqm_curv}. Indeed, despite the
approximation made, equation set~\eq{setofgceqs} can still be cast to the non-canonical Hamiltonian
form
\be{hamform}
\frac{{\rm d} z^i}{{\rm d} \tau}=\Lambda^{ij}\difp{H}{z^j}, \qquad \Lambda^{ij}({\bf z})=\left\{z^i,z^j\right\}_\tau,
\end{equation}
where the Hamiltonian function is $H({\bf z})=v_\parallel^2/2-U({\bf x})$ and $\Lambda^{ij}({\bf z})$ is an antisymmetric Poisson matrix.
The latter is linked to Poisson brackets that are slightly re-defined from those in
Ref.~\onlinecite{grebogi_relativistic_1984} due to a different orbit parameter,
\be{poisson}
\left\{f,g\right\}_\tau = b_\ast^i\left(\difp{f}{x^i}\difp{g}{v_\parallel}-\difp{g}{x^i}\difp{f}{v_\parallel}\right)
+\varepsilon^{ijk}\difp{g}{x^i}\difp{f}{x^j}\frac{B_k}{\omega_c},
\qquad
b_\ast^i = \varepsilon^{ijk}\left(\difp{A_k}{x^j}+v_\parallel \difp{}{x^j}\frac{B_k}{\omega_c}\right).
\end{equation}
In the derivation above, one occurrence of $v_\parallel^2$ has been replaced by $2U({\bf x})$ to obtain equation set~\eq{setofgceqs}.
As long as this equality in the first of~\eq{defU} is kept exact by a numerical scheme,
\redx{this formal violation of the Poisson structure doesn't affect the final result.
The present method integrates equations of motion to computer accuracy and thus
exactly conserves invariants with respect to piecewise linear fields in order
to fulfil this requirement. Invariants with respect to original smooth fields
oscillate within fixed bounds in a similar way as for conventional symplectic
integrators~\cite{hairer_geometric_2006}.}
\bluex{a geometric Poisson integrator follows,
being a generalization of symplectic integrators in non-canonical coordinates. Similar to symplectic integrators, Poisson integrators
preserve the symplectic structure and invariants of motion within fixed bounds.~\cite{hairer_geometric_2006}}
\red{
Conservation properties and symplectic features of the locally linear
equation set \eq{standeqset} are discussed in more detail in Appendix~\ref{sec:appendix_hamiltonian}.
}

\begin{figure}[h]
	\centerline{\includegraphics[keepaspectratio,width=1.2\linewidth, trim=50 100 50 0, clip]{figure_01.eps}}
	\captionsetup{justification=raggedright,singlelinecheck=false,textfont=footnotesize,labelfont=footnotesize}
	\caption{Real space illustration of field-aligned grid for the plasma core of a toroidal fusion device. Magnetic field lines are traced in general curvilinear coordinates $(x^1$,$x^2$,$x^3)$ on $N_s$ flux surfaces. At equidistant spacing in toroidal and poloidal direction a 3D grid consisting of $N_s \times N_\vartheta \times N_\varphi$ hexahedra is created.\\
	(a) 2D poloidal projection of 3D grid with a marked hexahedron cell. (b) Magnification of marked hexahedron in (a) with indication how each hexahedron is split into six tetrahedral cells. \redx{The symmetry direction is along the coordinate $x^3$.}\blue{Resulting}\bluex{The resulting tetrahedra are }
\blue{oriented in symmetry direction $x^3$ such that linear interpolation of electromagnetic field enables exact field representation on the cell's vertices.}
(c) Six individual tetrahedral cells compose stackable hexahedron such that adjacent tetrahedra faces are congruent.}
	\label{fig:curvilinear_grid}
\end{figure}
\subsection{Numerical solution}
\label{sec:numerical_solution}

An approximate formal solution of set~\eq{standeqset} in a single cell is given as a polynomial series of the orbit parameter,
\be{explser}
{\bf z}(\tau)={\bf z}_0 + \sum\limits_{k=1}^K \frac{\tau^k}{k!} \left(\hat {\bf a}^{k-1}\cdot{\bf b} + \hat {\bf a}^k\cdot {\bf z}_0\right),
\end{equation}
where $\hat {\bf a}$ and ${\bf b}$ stand for matrix $a^i_l$ and vector $b^i$, respectively,
${\bf z}(0)={\bf z}_0$ is a starting point,
and the exact solution is obtained in the limit $K \rightarrow \infty$. \red{In the case $k=1$ the matrix $\hat {\bf a}^{k-1} = \hat {\bf a}^{0}$ is the identity matrix.}
\redx{It should be noted that for mild electric fields (validity domain of Eq.~\eq{eqm_curv})
where the potential energy $e_\alpha\Phi$ is of the order of the kinetic energy (sub-sonic rotations)
all elements of matrix $\hat{\bf a}$ except for $a^i_4$ scale linearly with the Larmor radius, and
in the zero Larmor radius limit the series expansion of the order $K =  2$ does already provide an exact solution.
Therefore, the series expansion~\eq{explser} converges rapidly in the case that the Larmor radius is small in comparison with
 the spatial scale of the electromagnetic field (see Eq.~\eq{RK4estimation} below).}

Since intersections with cell boundaries (tetrahedra faces) must be computed exactly,
the particle is pushed from cell boundary to cell boundary, if no intermediate position inside a cell,
e.g. after a pre-defined time step, is required deliberately.
\red{In the latter case, the positions before and after the time step are located inside the cells, while the time
step itself includes tracing of all orbit intersections with cell boundaries between these positions.}

An orbit intersection with a tetrahedron is found on exit as an intersection with one of
four planes
\be{planes}
F^\alpha({\bf z})\equiv n^{(\alpha)}_i \left(x^i-x^i_{(\alpha)}\right)=0, \qquad \alpha=1,\dots,4
\end{equation}
reached in the smallest positive time,\blue{after entry (``dwell time'')} \red{``exit time'', from the starting position $x^i(\tau_0)$ located either on the cell boundary or inside the cell }, see Fig.~\ref{fig:triangle01}.
Here, $n^{(\alpha)}_i$ and $x^i_{(\alpha)}$ are the (constant)
normal to the plane containing tetrahedron face $\alpha$ and coordinates of some vertex
on that face, respectively.
\red{In the case that the starting position is located on the cell boundary (due to boundary-boundary-pushings as in Fig.~\ref{fig:triangle01}), the exit time coincides with the ``dwell time'' of the particle inside the cell. If an intermediate stop (inside the cell) is deliberately required, the dwell time is the sum of the time to reach this stop (``entry time'') and the exit time.}

With substitution of the orbit, ${\bf z}={\bf z}(\tau)$, Eqs.~\eq{planes}
are nonlinear equations with respect to the orbit parameter, $F^\alpha({\bf z}(\tau))=0$,
which should be solved numerically. They become algebraic and can be solved analytically
if an approximate solution~\eq{explser} is used with $K\le 4$\red{, as explained below in detail.}
\begin{SCfigure}
	\centering{	\includegraphics[width=0.5\textwidth,keepaspectratio,trim=200 100 200 100, clip]{figure_02.eps}}
	\captionsetup{textfont=footnotesize,labelfont=footnotesize,justification = raggedright,singlelinecheck = false}
	 \caption{Illustration of intersections of the orbit $x^i(\tau)$ with planes confining the cell. \blue{For simplification, the }\red{The} three-dimensional tetrahedral cell is depicted as a two-dimensional triangle\blue{.}\red{ in the interest of simplification.} The particle enters the cell at $x^i(\tau_0)$ and leaves the cell at \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {1}}}. Other intersections at the points \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {2}}} and \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {3}}} are not realized. The correct orbit-tetrahedron intersection~\raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {1}}} is reached in the smallest positive time among all intersections.}
	\label{fig:triangle01}
\end{SCfigure}

Both cases have been implemented in a Fortran program with the name \blue{Geometric}\red{\textbf{G}uiding-center}
\textbf{OR}bit \textbf{I}ntegration with \textbf{L}ocal \textbf{L}inearization \textbf{A}pproach
(\bluex{\textit{Gorilla}}\redx{\textbf{GORILLA}}).

\redx{As a matter of completeness, it is worth mentioning that a third similar to Ref.~\onlinecite{kasilov_geometric_2016} variant exists
to solve the linear equation set~\eq{standeqset} exactly in terms of exponential functions
of eigenvalues in the eigenvector basis (as any other linear equation set). Formally, this solution
corresponds to the limit $K \rightarrow \infty$ in Eq.~\eq{explser} where infinite sums
can be expressed in terms of the exponential function $\exp(\tau \hat a)$. This method, however,
has two drawbacks. First, Eqs.~\eq{planes} turn into nonlinear transcendental equations
which have to be solved numerically up to computer accuracy. Second, the analytical solution results
in strong cancellation errors in the case of small Larmor radii which is exactly the case where the finite
series expansion results in negligible errors. Therefore, we refrain from such a solution.}

In its first variant (\bluex{\textit{Gorilla RK4}}\redx{GORILLA RK4}), the orbit intersections with tetrahedra faces are
computed numerically by solving equation set (\ref{standeqset}) in each cell with a single step of
the Runge-Kutta 4 method embedded into an iterative scheme
using Newton's
method to obtain the integration time $\Delta\tau$ required to reach the cell boundary
(\red{exit time $\tau_e=\Delta \tau$}\blue{dwell time $\tau_d=\Delta \tau$}).
\blue{Predominantly, this}\red{This} iterative scheme \red{predominantly }converges after roughly two Newton iterations, due to an
analytic estimation for the necessary initial
step length using an approximate parabolic solution of ODE set (\ref{standeqset}) taken in zeroth order
in Larmor radius, $a^i_l = 0$ for $1 \le i,j \le 3$ and $a^4_4 v_\parallel(\tau) = a^4_4 v_{\parallel,0}$.

Nevertheless, in \red{the }case of numerically challenging orbits (tangential to a tetrahedron face or
almost intersecting with a tetrahedron's edges or vertices),
several special cases appear \blue{where}\red{in which} the iterative scheme \blue{doesn't}\red{does not} converge, and those cases must be treated
separately in a computationally more expensive manner.
Since such cases appear only rarely, the \blue{additional}\red{additionally} required computational effort is negligible in
comparison to the standard \blue{one}\red{procedure}.

A single RK4 integration step per iteration is sufficient because the magnetic field is uniform
within a cell. Respectively, the error of the RK4 method strongly
scales with the Larmor radius $\rho$ and can be brought below computer accuracy by a moderate
grid refinement.
Namely, in a tokamak geometry the error of a single step traversing the cell can be estimated as
\be{RK4estimation}
\frac{\delta R(\Delta\tau)}{a} \sim \frac{\rho^3}{q^4 R^3} \Delta \varphi^5,
\end{equation}
with $R$, $a$, $q$ and $\Delta\varphi$ denoting major radius, plasma radius, safety
factor and toroidal cell length, respectively.

Due to the scaling of the error with $\rho$, for particles with mild energies (thermal electrons
and ions) a less elaborate algorithm (\bluex{\textit{Gorilla Poly}}\redx{GORILLA Poly}) can be realized.
By truncating the summation of Eq.~\eq{explser} at $K$~=~2, 3~or~4 one obtains
approximate solutions of various orders in Larmor radius.
With these solutions, equations for the \blue{dwell time $\tau_d$}\red{exit time $\tau_e$}
are algebraic equations which are solved analytically up to order $K=4$, by finding the
smallest positive root of
\be{root_explser}
{\bf n}^\alpha \cdot \left( {\bf z}_0 + \sum\limits_{k=1}^K
\frac{\tau_{\red{e}\blue{d}}^k}{k!} \left(\hat {\bf a}^{k-1}\cdot{\bf b} + \hat {\bf a}^k\cdot {\bf z}_0\right)
-{\bf z}^\alpha
\right) = 0,
\end{equation}
where ${\bf n}^\alpha = \left(n_1^{(\alpha)},n_2^{(\alpha)},n_3^{(\alpha)},0\right)$
are face normals (see above) and
${\bf z}^\alpha=\left(x^1_{(\alpha)},x^2_{(\alpha)},x^3_{(\alpha)},0\right)$.
Furthermore, by using the quadratic polynomial solution $K$~=~2 of Eq.~\eq{root_explser},
the appropriate orbit intersection plane $\alpha$ can be predicted for higher orders
$K$ = 3 or 4, which predominantly reduces the number of higher order root finding operations
from four to one.
In the numerical implementation,
series~\eq{explser} which contains matrix products is not directly evaluated as written above. Instead, various sub-products
are preliminarily evaluated and stored for grid cells\red{,} \blue{which minimizes}\red{minimizing} the number of matrix products
\blue{needed}\red{that need} to be evaluated directly.
The finite series solution~\eq{explser}
allows to analytically evaluate also various integrals
over the dwell time needed for scoring of macroscopic plasma plasma parameters in Monte Carlo
procedures.
In Appendix~\ref{sec:appendix_integrals}, such integrals of $v_\parallel$ and
of $v_\perp^2$ and $v_\parallel^2$ are given,
which respectively determine parallel plasma flow and
components of the pressure tensor in Chew-Goldberger-Low form
essential for computation of equilibrium plasma currents.

\section{Collisionless guiding-center orbits}
\label{sec:orbits}

\subsection{Guiding-center orbits in an axisymmetric tokamak field}
\label{sec:orbits_2d}

In this section the results of \red{quasi-}geometric orbit integration computed with \bluex{\textit{Gorilla}}\redx{GORILLA} and the comparison to an exact orbit computed with a usual adaptive RK4/5 integrator are presented for an axisymmetric tokamak field configuration of ASDEX Upgrade (shot 26884 at 4300 ms) described in Ref.~\onlinecite{heyn_quasilinear_2014}. The adaptive RK4/5 integrator requires high-order interpolation of electromagnetic fields with continuous derivatives, e.g. with help of 3D-splines, instead of continuous piecewise linear functions \red{as }in \red{the }case of the \red{quasi-}geometric integration method.

In axisymmetric configurations, the shape of the orbit is fully determined by three conservation laws,
$p_\varphi = \mathrm{const.}$, $J_\perp = \mathrm{const.}$ and $w =  \mathrm{const.}$, where
\be{ptor}
p_\varphi\redx{ = \frac{e}{c} A_\varphi^\ast} = m v_\parallel \frac{B_\varphi}{B} + \frac{e}{c}A_\varphi
\end{equation}
is \redx{the} canonical toroidal angular momentum. \redx{The conservation of $p_\varphi$ is obvious from Eq.~\eq{eqm_curv}
since after its substitution in $\dot p_\varphi=\dot x^i \partial p_\varphi / \partial x^i$ all the
terms in the resulting expression are proportional to partial derivatives over $\varphi$ of various $A^\ast_k$ components.}

In geometric/symplectic numerical integration schemes \bluex{this property is}\redx{these conservation properties are} retained [\onlinecite{hairer_geometric_2006}], which means that orbits must remain closed in the poloidal projection.

\bluex{Due to the special formulation of Eq.~\eq{defU}, where $v_\parallel$ is purely a function of position, the perpendicular adiabatic invariant and the total energy are conserved naturally.} \bluex{By properly orienting the tetrahedra with respect to the symmetry direction, axisymmetry in \red{the }case of 2D fields
is exactly preserved upon linearization.}
\redx{
In order to preserve $p_\varphi$ upon linearization, the tetrahedral cells of the method's underlying
grid must be specially oriented with respect to the symmetry direction. Namely, this is achieved with
stackable hexahedra each consisting of two triangular prisms which are both subsequently split into
3 tetrahedral cells.  The prisms are oriented such that all triangular prism faces lie on
$x^3= \mathrm{const.}$ planes where $x^3$ is the symmetry direction.
Thus, each tetrahedron face which lies on a $x^3= \mathrm{const.}$ plane is congruent with
all other tetrahedron faces opposing it in the symmetry direction.
This specific splitting realization can be seen in Fig.~\ref{fig:curvilinear_grid}.
}
Consequently, the canonical toroidal angular momentum remains invariant \redx{in the presented method}\redx{, since the linearization of the electromagnetic field within a tetrahedral cell does not introduce an interpolation error in the symmetry direction. Namely, partial derivatives of the field quantities with respect to the symmetry direction remain zero. In the following, cylindrical ($R$, $\varphi$, $Z$)
and symmetry flux coordinates ($s$, $\vartheta$, $\varphi$) of Ref.~\onlinecite{dhaeseleer_flux_2012} are used, where $\varphi = x^3$ in both coordinate systems.}

Fig. \ref{fig:orbit_comparison} depicts Poincaré plots ($\varphi = 0$) of trapped thermal ion orbits making $10^7$ toroidal turns which are integrated
by different methods from the same starting conditions. \blue{G}\red{Quasi-g}eometric integration using an iterative scheme with RK4 integration and Newton steps has been performed in cylindrical\bluex{($R$, $\varphi$, $Z$)} and symmetry flux coordinates\bluex{ ($s$, $\vartheta$, $\varphi$) of Ref.~\onlinecite{dhaeseleer_flux_2012}} and is compared to the exact orbit.
It can be seen, that the coarseness of the grid leads to slightly differently shaped orbits obtained in different coordinate systems, whereas the
effect of the integration error (\ref{RK4estimation}) is negligible.
\begin{figure}[h]
	\centerline{\includegraphics[keepaspectratio,width=1.2\linewidth]{figure_03.eps}}
	\captionsetup{justification=raggedright,singlelinecheck=false,textfont=footnotesize,labelfont=footnotesize}
	\caption{(a) Poincaré plot ($\varphi=0$, $10^7$ toroidal mappings) of a trapped 1.5~keV D ion in axisymmetric ASDEX Upgrade configuration with a tetrahedral grid size of 20x20x20. Two-dimensional Poincaré sections of orbits obtained with different integration methods are indicated with markers: Exact \blue{orbit:}\red{orbit (adaptive RK4/5):} $\blacktriangle$, \blue{3D Geometric Integrator }\bluex{\textit{Gorilla}}\redx{GORILLA} with cylindrical coordinates: $\blacksquare$, \bluex{\textit{Gorilla}}\redx{GORILLA} with symmetry flux coordinates: $\blacklozenge$. (b) and (c) are magnifications of the pertinent zones in (a).}
	\label{fig:orbit_comparison}
\end{figure}

In Fig.~\ref{fig:tokamak_passing_toroidal_angular} results of the \red{quasi-}geometric integration of a passing high energy ion (300~keV) using the polynomial series solution~\eq{explser} are shown for the same axisymmetric geometry and compared to the exact orbit.
Symmetry flux coordinates are used for the \red{quasi-}geometric integration, the results are then converted to cylindrical coordinates
in the first plot (a). In general, the quartic polynomial solution ($K$~=~4) of ODE set~\eq{standeqset} is equivalent to the
numerical solution using Runge-Kutta 4, thus, the result of the latter is omitted in the figure.
Moreover, it can be seen that even for high
energy ions the series expansions of third and fourth order are already accurate enough in order to fulfill the condition
$p_\varphi=\text{const.}$ over 10$^6$~toroidal mappings. However, the second order series expansion shows a convective behavior, due to a systematic error from solving the ODE set~\eq{standeqset}.
\blue{Nevertheless, for the electrons which have much smaller Larmor radii, the second order series expansion is sufficient.}\red{Nevertheless, the second order series expansion is sufficient for electrons which have much smaller Larmor radii.}
\begin{figure}[h]
	\centerline{\includegraphics[keepaspectratio,width=1.1\linewidth,trim=60 0 60 0, clip]{figure_04.eps}}
	\captionsetup{justification=raggedright,singlelinecheck=false,textfont=footnotesize,labelfont=footnotesize}
	\caption{(a \& b) Poincaré plot ($\varphi=0$, $10^6$ toroidal mappings) of a passing 300 keV D ion in axisymmetric ASDEX Upgrade configuration evaluated by \blue{the 3D Geometric Integrator }\bluex{\textit{Gorilla}}\redx{GORILLA} with the analytical solution in form of a polynomial series truncated at $K=2$ ($\Diamond$), $K=3$ ($\square$) and $K=4$ ($\triangledown$) compared to the exact orbit ($\bullet$). Cylindrical coordinates are used for the Poincaré plot in (a), whereas symmetry flux coordinates are used in (b).\\
		(c) Canonical toroidal angular momentum $p_\varphi$ normalized to the value at $t=0$ is evaluated at the Poincaré sections in (a): Exact result (solid) is compared to polynomial series truncated at $K=2$ (dash-dotted), $K=3$ (dotted) and $K=4$ (dashed).}
	\label{fig:tokamak_passing_toroidal_angular}
\end{figure}

To obtain the orbits of thermal ions shown in Fig.~\ref{fig:axisymmetric_noise_plot}, uniformly distributed axisymmetric random noise ($\xi = 0 \dots 1$) is added to the electrostatic potential $\Phi^{\text{noisy}} = \Phi(1+\epsilon_\Phi \xi)$, to the vector potential $A_k^{\text{noisy}} = A_k(1+\epsilon_A \xi)$ and simultaneously to both quantities, respectively. Here, $\epsilon$ is the relative magnitude of added noise. \red{The guiding-center orbits are evaluated in symmetry flux coordinates by \bluex{\textit{Gorilla}}\redx{GORILLA} with the analytical solution in form of a polynomial series truncated at $K=4$, the results are then converted to cylindrical coordinates.}
Even though relatively high noise (up to 30~\%) is added, the orbits keep a similar shape in comparison with the unperturbed ones and remain closed in the poloidal projection, meaning the condition $p_\varphi = \text{const.}$ is still fulfilled.
\begin{figure}[t]
	\centerline{\includegraphics[keepaspectratio,width=1.2\linewidth]{figure_05.eps}}
	\captionsetup{justification=raggedright,singlelinecheck=false,textfont=footnotesize,labelfont=footnotesize}
	\caption{Poincaré plot ($\varphi=0$) of trapped and passing 1.5 keV D ions in axisymmetric ASDEX Upgrade configuration with axisymmetric perturbation of electrostatic potential $\Phi$ and vector potential $A_k$. \red{Orbits are evaluated by \bluex{\textit{Gorilla}}\redx{GORILLA} with the analytical solution in form of a polynomial series truncated at $K=4$.}\\
The Poincaré plot for unperturbed electromagnetic fields is depicted in black as a reference. Uniformly distributed axisymmetric random noise ($\xi = 0 \dots 1$) is added in (a) to the electrostatic potential $\Phi^{\text{noisy}} = \Phi(1+\epsilon_\Phi \xi)$, in (b) to the vector potential	$A_k^{\text{noisy}} = A_k(1+\epsilon_A \xi)$ and in (c) to both quantities.
}
	\label{fig:axisymmetric_noise_plot}
\end{figure}

Therefore, the \red{quasi-}geometric integration method is suitable for self-consistent Monte Carlo modeling of 2D equilibrium plasma
parameters and electromagnetic fields computed from those parameters.
Further, it should be noted that the computational efficiency of the \red{quasi-}geometric integrator is not affected
by the presence of small-scale noise. However, this would be the case for an adaptive RK4/5 integration in combination
with high order smooth interpolation where the noise leads to small-scale oscillations.

The Poincaré plots ($\varphi = 0$) of Fig.~\ref{fig:tokamak_trapped_parallel_invariant} (a) \& (b) correspond to a trapped high energy ion (300~keV), where the orbits are integrated using the polynomial series expansion in several orders and compared to the exact orbit. Again, symmetry flux coordinates are used for the \red{quasi-}geometric integration, the results are then converted to cylindrical coordinates in the first plot (a). Despite the same starting conditions \blue{of}\red{for} all orbits, the magnification in (b) clearly depicts differences in the shape of the orbits caused by a finite grid size ($100 \times 100 \times 100$) and different orders of the series expansion.

Fig.~\ref{fig:tokamak_trapped_parallel_invariant} (c) shows the corresponding time evolution of the parallel adiabatic invariant~$J_\parallel$, which is defined by an integral over the distance $l$ passed along the field line during a single bounce period $t_b$ by a
trapped particle as follows,
\be{parinv}
J_\parallel = m \oint v_\parallel \mathrm{d}l = m \int_{0}^{t_b} v_\parallel^2 \left(t \right) \mathrm{d}t.
\end{equation}
For an exact orbit, $J_\parallel$ is a conserved quantity. Symplectic orbit integration does not lead to an error accumulation in the invariants of motion [\onlinecite{hairer_geometric_2006}], while systematic changes can only arise from numerical errors in solving the ODE. Even for high energy ions the series expansions of third and fourth order are already accurate enough in order to fulfill the condition $J_\parallel = \mathrm{const.}$ for $10^6$ bounce periods. Hence, in an axisymmetric configuration excellent long-term orbit dynamic is demonstrated by the \red{quasi-}geometric orbit integration method, as long as the ODE set~\eq{standeqset} is solved accurately.

\begin{figure}[h]
	\centerline{\includegraphics[keepaspectratio,width=1.2\linewidth]{figure_06.eps}}
	\captionsetup{justification=raggedright,singlelinecheck=false,textfont=footnotesize,labelfont=footnotesize}
	\caption{(a \& b) Poincaré plot ($\varphi=0$, $10^6$ bounce periods) of a trapped 300 keV D ion in axisymmetric ASDEX Upgrade configuration evaluated by \blue{the 3D Geometric Integrator }\bluex{\textit{Gorilla}}\redx{GORILLA} with the analytical solution in form of a polynomial series truncated at $K=2$ ($\Diamond$), $K=3$ ($\square$) and $K=4$ ($\triangledown$) compared to the exact orbit ($\bullet$). Cylindrical coordinates are used for the Poincaré plot in (a), whereas symmetry flux coordinates are used in (b). The magnification in (b) shows the banana tip and clearly depicts the difference of polynomial orders. \\
	(c) The parallel adiabatic invariant $J_\parallel$ normalized to the value at $t=0$ is depicted for $10^6$ bounce periods: Exact result (solid) is compared to polynomial series truncated at $K=2$ (dash-dotted), $K=3$ (dotted) and $K=4$ (dashed).}
	\label{fig:tokamak_trapped_parallel_invariant}
\end{figure}
\subsection{Guiding-center orbits in three-dimensional fields}
\label{sec:orbits_3d}

It can be seen from Eqs.~\eq{standeqset} and~\eq{amatdef} that in the field line limit $\omega_c\rightarrow\infty$
guiding-center orbits are straight within spatial cells and magnetic field lines are represented by polygonal chains, respectively.
In case of 3D magnetic fields described in non-aligned spatial coordinates, the existence of embedded KAM surfaces
is not obvious for such an approximate representation of the field even in \red{the }case \red{where }these surfaces exist in the exact
system. In order to study artificial chaos induced by the linearization,
\red{quasi-}geometric orbit integration has been performed for low energy particles with negligible FLR effects in symmetry flux coordinates
$(s,\vartheta,\varphi)$
associated with the axisymmetric tokamak field of the previous section with a harmonic perturbation added
to the toroidal co-variant component of the axisymmetric vector potential,
\be{example_onemode}
A_\varphi=\psi_{\rm pol}(s)(1+\varepsilon_M\cos(m_0\vartheta+n_0\varphi)).
\end{equation}
The harmonic indices $m_0=n_0=2$ used in the testing correspond to a non-resonant perturbation
which leads only to a corrugation of the magnetic surfaces\red{,} \blue{such}\red{with the effect} that they are \blue{not aligned}\red{no longer aligned} with the coordinate
surfaces $s=\mathrm{const.}$\blue{ anymore}.

\noindent
In Fig.~\ref{fig:tokamak_helical_pert_poincare}, Poincar\'e plots of magnetic field lines obtained by the \red{quasi-}geometric
integration method for this perturbed configuration are shown at the cross section $\varphi=0$ together with a cross section
of one exact corrugated flux surface. It can be seen that the orbits from the \red{quasi-}geometric integrator become more chaotic with
increasing perturbation amplitude $\varepsilon_M$.

\begin{figure}[h]
	\centerline{\includegraphics[keepaspectratio,width=1.2\linewidth]{figure_07.eps}}
	\captionsetup{justification=raggedright,singlelinecheck=false,textfont=footnotesize,labelfont=footnotesize}
	\caption{Poincaré plots of the orbits in zero Larmor radius limit (field lines)
	for $10^4$ (upper row) and $10^5$ (lower row) toroidal mappings and various perturbation amplitudes indicated
	in the titles.
	Orbits start at 34 equidistant flux surfaces between $s = 0.5$ and  $s=0.7$ and are evaluated by \bluex{\textit{Gorilla}}\redx{GORILLA}
	with polynomial series truncated at $K=2$ and angular grid size of $N_\vartheta = N_\varphi = 30$.
	Similar results are achieved for angular grid size with incommensurable dimensions, e.g. $N_\vartheta = 29, N_\varphi = 31$.
	Solid line shows a cross-section of one exact corrugated flux surface.}
	\label{fig:tokamak_helical_pert_poincare}
\end{figure}
Diffusive behavior of orbits can be characterized by the variance of
the normalized toroidal flux $s$, accumulated over the time for the ensemble of test particles
starting from the same perturbed flux surface. This variance is described by the magnetic field line diffusion
coefficient $D_M^{ss}$ as $\left\langle \delta s^2\right\rangle=2 D_M^{ss} N$ where $N$ is the number of
toroidal orbit turns.
The effective diffusion coefficient $D_M^{ss}$ computed from the
orbits has a strong inverse scaling with poloidal $N_\vartheta$ and toroidal $N_\varphi$ grid
sizes and, furthermore, shows in general a small magnitude of diffusion even at coarse angular grid resolution.
E.g., at an angular grid size of $N_\vartheta = N_\varphi = 30$ the effective diffusion coefficient $D_M^{ss}$ is in the order of $10^{-13}$, $10^{-11}$ and $10^{-9}$
for relative perturbation amplitudes of $\varepsilon_M = 1~\%$, $3~\%$ and $10~\%$, respectively.
This level of field line diffusion is roughly five orders of magnitude smaller than observed for the initial version of the
code~\citep{eder_three-dimensional_2019}, and, in the worst case of $\varepsilon_M = 0.1$, results in stochastic diffusion
of electrons with the coefficient $D_\perp \sim D_M^{ss} v_\parallel r^2 / R \sim 100$ cm$^2$s$^{-1}$. For smaller $\epsilon_M$
values of the typical order for external magnetic perturbations in tokamaks this numerical diffusion is below the level of
classical electron diffusion and, therefore, can be safely ignored.

In \red{the }case of field aligned coordinates, chaotization of passing orbits (lines of force of the effective field ${\bf B}^\ast$)
can only be caused by the cross-field drift.
Such a case is tested below for a strong violation of axial symmetry using as an example the stellarator field configuration
described in Ref.~\onlinecite{drevlak_quasi-isodynamic_2014}, namely,
a quasi-isodynamic reactor-scale device with five toroidal field periods and \red{a }major radius of 25 m.
Here the magnetic field has been normalized so that its modulus averaged over Boozer coordinate angles on the starting
surface is $B_{00}$~=~5~T.
Guiding-center orbits were computed with the \red{quasi-}geometric integration method in symmetry flux coordinates
for strongly passing electrons and ions with $v_\parallel / v = 0.9$  at the starting point on the flux surface $s = 0.6$ with an energy of 3~keV.
The numerical diffusion observed for a rather coarse angular grid with the size $N_\vartheta = N_\varphi = 30$
is roughly seven orders of
magnitude smaller than the minimum level of the neoclassical mono-energetic diffusion coefficient $D_{11}$
evaluated in section~\ref{sec:transport} for the same device with particles of the same energy.
\begin{figure}[h]
	\centerline{\includegraphics[keepaspectratio,width=1.2\linewidth]{figure_08.eps}}
	\captionsetup{justification=raggedright,singlelinecheck=false,textfont=footnotesize,labelfont=footnotesize}
	\caption{(a \& b) Poloidal projection of Poincaré sections at $v_\parallel=0$ switching sign from $-$ to $+$ of a trapped 3~keV~D ion in 3D stellarator field configuration. Orbits evaluated by \bluex{\textit{Gorilla}}\redx{GORILLA} with the analytical solution in form of a polynomial series truncated at $K=2$ ($\Diamond$), $K=3$ ($\square$) and $K=4$ ($\triangledown$) are compared to the exact orbit ($\bullet$). Poincaré plot is depicted in cylindrical coordinates (a) and symmetry flux coordinates (b).\\
				\blue{(c) Evolution of the parallel adiabatic invariant $J_\parallel$ normalized to the value at $t=0$ in (a) and (b) over $10^5$~bounce times. Lines in (c) represent linear fits to the $J_\parallel$ data points: Exact result (solid) is compared to polynomial series truncated at $K=2$ (dash-dotted), $K=3$ (dotted) and $K=4$ (dashed).}}
	\label{fig:hydra_trapped_poincare}
\end{figure}

For the visualization of a trapped particle orbit, we use orbit footprints on Poincaré sections defined by the condition
$v_\parallel(\tau)~=~0$, i.e. phase space hypersurfaces containing orbit turning points.
\blue{Out of two kinds}\red{From the two types} of these surfaces, those are chosen \blue{where}\red{in which} the sign of $v_\parallel$ changes from negative to positive.
Fig.~\ref{fig:hydra_trapped_poincare} depicts a poloidal projection of orbit footprints corresponding
to a trapped ion with an energy of 3 keV. The orbits have been integrated in symmetry flux coordinates using the polynomial
series expansion in several orders and compared to the exact orbit computed with an adaptive RK4/5 integrator.
Further, the poloidal coordinates of the footprints have also been converted to cylindrical coordinates and visualized in
both coordinate systems, respectively.
Despite the same starting conditions \blue{of}\red{for} all orbits, the magnification in (b) clearly depicts \red{slight }differences in the shape of the orbits caused by a finite grid size ($100 \times 100 \times 100$) and different orders of the series expansion.

\begin{figure}[h]
	\centerline{\includegraphics[keepaspectratio,width=1.2\linewidth]{figure_09.eps}}
	\captionsetup{justification=raggedright,singlelinecheck=false,textfont=footnotesize,labelfont=footnotesize}
	\caption{\red{
				(a) Evolution of the parallel adiabatic invariant $J_\parallel$ normalized to the value at $t=0$ of the trapped particle orbit in Fig.~\ref{fig:hydra_trapped_poincare} over $4 \cdot 10^5$~bounce times. Lines in (a) represent window filtered trend of the $J_\parallel$ data points: Exact result ($\triangle$) is compared to polynomial series truncated at $K=2$ ($\Diamond$), $K=3$ ($\square$) and $K=4$ ($\triangledown$).\\
				(b) Modulus of the relative error of the parallel adiabatic invariant $J_\parallel$ averaged over $10^5$~bounce times as a function of the angular grid size $N_\varphi = N_\vartheta$. Orbits are evaluated by \bluex{\textit{Gorilla}}\redx{GORILLA} with the analytical solution in form of a polynomial series truncated at $K=2$ ($\Diamond$), $K=3$ ($\square$) and $K=4$ ($\triangledown$). The fits of the results are depicted with lines in accordance with the legend.}
				}
	\label{fig:hydra_trapped_parallel_invariant}
\end{figure}

In contrast to the axisymmetric tokamak field of the previous section, the parallel adiabatic invariant $J_\parallel$ is not
an exact invariant in a stellarator. Nevertheless, it should be well preserved as long as the trapped orbit stays within the same
class (traverses the same number of field minima over its bounce period) which is the case here.
Fig.~\ref{fig:hydra_trapped_parallel_invariant} \blue{(c)}\red{(a)} shows the time evolution of $J_\parallel$ for \blue{$10^5$}\red{$4 \cdot 10^5$} bounce periods of the corresponding orbits \red{of Fig.~\ref{fig:hydra_trapped_poincare}}. It can be seen that the results for the exact orbit shows no visible deviation for this configuration.
\red{The truncation of the polynomial series at $K=2$ causes a violation of the system's Hamiltonian structure.
This is manifested by an attractor in phase space, what is clearly visible in Fig.~\ref{fig:hydra_trapped_parallel_invariant} (a). A detailed analysis of the corresponding Poincaré sections when the attractor is already fully established (last $10000$ evaluations) reveals that the orbit strictly follows a continuous curve, staying in the same class.

However, truncation at higher orders ($K=3$ \& $K=4$) violates the Hamiltonian structure only negligibly and thus does not lead to non-Hamiltonian features. Nevertheless, in contrast to the exact orbit, $J_\parallel$ evaluated by \bluex{\textit{Gorilla}}\redx{GORILLA} with these polynomial orders is not accurately conserved. In particular, the value of $J_\parallel$ meanders randomly around the exact value which is caused by the diffusive behavior of the orbit induced by the piecewise linearization of the electromagnetic field.

Fig.~\ref{fig:hydra_trapped_parallel_invariant} (b) shows the modulus of the relative error of the parallel adiabatic invariant $J_\parallel$ averaged over $10^5$~bounce times as a function of the angular grid size $N_\varphi = N_\vartheta$ varied from $10$ to $200$.
At a moderate angular grid size of $N_\vartheta = N_\varphi = 28$ (in accordance with the Nyquist-Shannon sampling theorem~\cite{nyquist_certain_1928,shannon_mathematical_1948} explained below) this mean relative error stays already below $10^{-3}$ for the orders $K=3$ \& $K=4$.
Thus, the finite grid size used in the \red{quasi-}geometric orbit computation does not lead to a significant error accumulation.
With regard to the comparatively large number of bounce times, even the order $K = 2$ which is naturally the fastest with respect
to CPU time yields quite accurate results at sufficient grid resolution.
}

\blue{Furthermore, the finite grid size used in the geometric orbit computation does not lead to a significant error accumulation even for the order $K = 2$ which is naturally the fastest with respect
to CPU time.}

\section{Monte Carlo evaluation of neoclassical transport coefficients, performance benchmark}
\label{sec:transport}

Evaluation of neoclassical transport coefficients using the Monte Carlo
method~\cite{boozer_monte_1981,lotz_monte_1988} is widely
used for stellarators and tokamaks with 3D perturbations of the magnetic
field~\cite{wakasa_study_2008,tribaldos_monte_2001,isaev_venusf_2006,allmaier_variance_2008,velasco_calculation_2011,satake_neoclassical_2011,pfefferle_venus-levis_2014}.
An advantage of this method in its original, full-$f$ form is the use of test particle guiding-center
orbits without requiring \red{the }model simplifications needed in (more efficient) local approaches. \blue{Therefore}\red{The} Monte Carlo methods \red{thus }provide
an unbiased reference point in cases where those simplifications affect the transport \red{in a manner }such as \red{that }for regimes with
significant role of the tangential magnetic \red{drift~\cite{matsuoka_effects_2015,huang_benchmark_2017,velasco_knosos_2020}.}
An obvious disadvantage is that for realistic magnetic configurations Monte Carlo methods are CPU-intensive
with most of the CPU time spent for the integration of the guiding-center motion. The application of the proposed \red{quasi-}geometric integration method for this purpose instead of the usual Runge-Kutta method results in a visible speed-up
of the computations without significantly biasing the results.
Here, this application is made for benchmarking purposes \blue{assuming}\red{on the assumption} that the inaccuracies in orbit integration
which are tolerable in computations of transport coefficients are also tolerable
in global modelling of macroscopic plasma parameters.

The proposed orbit integration method is applied
within a standard Monte Carlo algorithm~\cite{boozer_monte_1981}
using the Lorentz collision model for the evaluation of the mono-energetic radial diffusion coefficient $D_{11}$.
The latter is determined via the average square deviation
of the normalized toroidal flux $s$ from its starting value $s_0$ as follows,
\be{eq:def_d11}
D_{11} = \frac{1}{2 t} \langle\left(s(t) - s_0 \right)^2\rangle.
\end{equation}
Here, angle brackets $\langle\dots\rangle$ denote an ensemble average, $s(0)=s_0$, and the test particle tracing time $t$ is chosen to
be \blue{larger}\red{longer} than the local distribution function relaxation time $\tau_{\rm rel}$
\blue{and smaller}\red{but shorter} than the radial transport time, $t=10 \tau_{\rm rel}$.
A Monte Carlo collision operator identical to that of Ref.~\onlinecite{boozer_monte_1981} is applied here in-between
constant collisionless orbit integration steps $\Delta t$. These steps are small enough compared to the typical bounce time
$\tau_\mathrm{b}$ and collision time $\tau_\mathrm{c}$,
\be{eq:def_time_step_MC}
\Delta t = \min \left(\frac{\tau_\mathrm{b}}{20},\frac{\tau_{\mathrm{c}}}{20} \right).
\end{equation}
Here, $\tau_{\mathrm{c}} = 1/\nu$ and $\tau_\mathrm{b} = 2 \pi R_0 / (v N_\mathrm{tor.})$ with $\nu$,  $R_0$, $v$ and $N_\mathrm{tor}$ denoting collisional deflection frequency, major radius, particle velocity and
number of toroidal field periods, respectively. The relaxation time $\tau_{\rm rel}$ is determined as the largest of
$\tau_{\mathrm{c}}$ and $\tau_\mathrm{b}^2/\tau_{\mathrm{c}}$.

\begin{figure}[t]
	\centerline{\includegraphics[keepaspectratio,width=1.1\linewidth]{figure_10.eps}}
	\captionsetup{justification=raggedright,singlelinecheck=false,textfont=footnotesize,labelfont=footnotesize}
	\caption{Mono-energetic radial diffusion coefficients $D_{11}$ for electrons (top) and deuterium
	ions (bottom)
	as functions of (a) normalized collisionality $\nu^\ast$ and (b)  Mach number $v_E^\ast$.
	Lines of various styles (see the legends) - reference computation, markers - results of \red{quasi-}geometric
	integration with polynomial solution of the order $K$ for $K=2$ ($\Diamond$), $K=3$ ($\square$)
	and $K=4$ ($\triangledown$).
	Error bars indicate $95$~\% confidence interval.}
	\label{fig:transport_coefficient_electrons_dions}
\end{figure}
In the present example, the mono-energetic radial diffusion coefficient has been evaluated
for the quasi-isodynamic stellarator configuration~\cite{drevlak_quasi-isodynamic_2014}
used also for collisionless orbits in section~\ref{sec:orbits_3d}.
Guiding-center orbits were computed with the \red{quasi-}geometric integration method in symmetry flux coordinates
using polynomial series solutions of various orders $K$.
The grid size $N_s \times N_\vartheta \times N_\varphi = 100 \times 60 \times 60$ was selected to be appropriate \blue{to
minimize}\red{for minimizing} the numerical diffusion (see the previous section.) In a reference computation, guiding-center
equations~\eq{eqm_curv} in symmetry flux variables with electromagnetic field interpolated
by 3D qubic splines were integrated by an adaptive RK4/5 integrator.
\blue{In order to minimize statistical errors, computations have been performed for a large ensemble size
of $10000$ particles.}\red{Computations were performed for a large ensemble size of $10000$ particles in order to minimize statistical errors.}

The results for $D_{11}$ computed for 3 keV electrons and ions at $s_0=0.6$ are presented in
Fig.~\ref{fig:transport_coefficient_electrons_dions}.
Values of radial electric field $E_r$ and deflection frequency $\nu$, which determine
transport regimes, are respectively characterized here by two dimensionless
parameters~\cite{beidler_benchmarking_2011}, Mach number $v_E^\ast = c E_r /(vB_0 )$ and
collisionality $\nu^\ast = (R_0 \nu)/(\iota v)$, where $\iota$ is the rotational transform.
For the ions, in addition to the ${\bf E}\times{\bf B}$ rotation, also the tangential magnetic drift
plays a significant role which can be seen from the shift of the $D_{11}$ maximum on $v_E^\ast$
dependence. The results of \red{quasi-}geometric integration stay in agreement
with the reference computation within the 95 \% confidence interval in all cases even for the lowest
order polynomial solution $K=2$. \blue{Therefore, as shown below, a significant gain in the computation
time can be obtained in this kind of calculations.}\red{Therefore, in calculations of this kind a significant gain in the computation
time can be obtained as shown below.}

\begin{figure}[t]
	\centerline{\includegraphics[keepaspectratio,width=1.1\linewidth]{figure_11.eps}}
	\captionsetup{justification=raggedright,singlelinecheck=false,textfont=footnotesize,labelfont=footnotesize}
	\caption{Relative error of mono-energetic radial transport coefficient
	$D_{11}$ of electrons (top) and D~ions (bottom) vs. relative CPU time.
	\blue{Compared}\red{The compared} orbit integration methods are:
	Runge-Kutta 4 ($\star$),
	Adaptive RK4/5 with various relative errors indicated in the plot ($\times$),
	\red{quasi-}geometric integration with polynomial solution
	(\bluex{\textit{Gorilla Poly}}\redx{GORILLA Poly}) of the order $K=2$ ($\Diamond$), $K=3$ ($\square$) and $K=4$
	($\triangledown$), and with RK4 solution (\bluex{\textit{Gorilla RK4}}\redx{GORILLA RK4}, $\triangle$).
	\blue{Fits}\red{The fits} of \red{the }results are depicted with lines \blue{according to}\red{in accordance with} the legend.
	\blue{Random}\red{The random} error of the reference result, $D_{11,\rm{ref}}$,
	is depicted as a horizontal line limiting its 95~\% confidence interval.
	}
	\label{fig:cpu_speed_electrons_dions}
\end{figure}
Moreover, we compare the performance and scaling for parallel computation of guiding-center orbits
using the \red{quasi-}geometric orbit integration method with
computations using standard reference integrators (RK4 and adaptive RK4/5).
For this, different integrators have been used within $D_{11}$ computation described above
for a particular choice of dimensionless parameters, $v_E^\ast = 10^{-3}$ and $\nu^\ast =  10^{-3}$,
and an increased ensemble size of $30000$ test particles.
The numerical experiment has been performed on a single node of the COBRA cluster
of MPCDF with 40 CPU cores (Intel Xeon Gold 6126) running 80 concurrent threads with hyperthreading.

The reference value for the transport coefficient, $D_{11,\rm{ref}}$,
and the reference CPU time are obtained by orbit integration
with an adaptive RK4/5 integrator with a relative tolerance of $10^{-9}$.
The accuracy of the $D_{11}$ evaluation using different computation parameter settings
is represented by the relative error $\delta D_{11} / D_{11,\rm{ref}}$
where $\delta D_{11} = |D_{11} - D_{11,\rm{ref}}|$.
The CPU time purely used for orbit integration serves as
a measure for the computational effort of the methods.
This given CPU time does not contain any overhead operations, e.g. the construction of the grid, generation of random numbers for
pitch-angle scattering and computation of $D_{11}$ by evaluating
Eq.~\ref{eq:def_d11} with the help of a least-squares regression.

Fig.~\ref{fig:cpu_speed_electrons_dions} shows the relative error of the mono-energetic
radial transport coefficient versus the relative CPU time of computations using the
\red{quasi-}geometric orbit integration
method with the polynomial series solution of various orders, \bluex{\textit{Gorilla Poly}}\redx{GORILLA Poly},
and the iterative scheme with RK4 integration and Newton steps, \bluex{\textit{Gorilla RK4}}\redx{GORILLA RK4}.
Accuracy and CPU time of \red{quasi-}geometric orbit integrations have been varied by
mutually changing the angular grid size $N_\vartheta \times N_\varphi$
from $8 \times 8$ to $60 \times 60$ while keeping the radial grid size
constant at $N_s = 100$. In the stellarator configuration of
Ref.~\onlinecite{drevlak_quasi-isodynamic_2014} used here, the number of toroidal harmonic modes
per field period is $14$, leading to a minimum toroidal grid size $N_\varphi = 28$ in
order to satisfy the Nyquist-Shannon sampling theorem~\cite{nyquist_certain_1928,shannon_mathematical_1948}.
Therefore, regression lines are drawn for the range of data points with grid sizes
from $8 \times 8$ until $28 \times 28$, clearly showing a convergent behavior
of $D_{11}$ with increasing grid refinement.
Furthermore, the adaptive RK4/5 integration is additionally performed with
relative tolerances of $10^{-3}$, $10^{-6}$, $10^{-7}$ and $10^{-8}$, respectively.
Note that the computational
speed of the adaptive RK4/5 integration with a relative tolerance of $10^{-6}$ cannot be
increased by higher relative tolerances, e.g. $10^{-3}$, since the macroscopic Monte Carlo
time step, $\Delta t$ Eq.~\eq{eq:def_time_step_MC}, is already elapsed within a single
RK4/5 step with sufficient accuracy. Hence, also the non-adaptive Runge-Kutta~4 method
is tested, which naturally needs one field evaluation less per time step than RK4/5.
In all cases, the relative error of \red{the }RK4/5 and Runge-Kutta~4 results is determined here
mainly by statistical deviations, with a random error dominating the bias.

\blue{Besides}\red{Apart from} statistical errors due to Monte Carlo sampling, a limit for capturing
all toroidal and poloidal field harmonics is given by a minimum grid size of
two points per period due to the Nyquist–Shannon sampling theorem.
Fig.~\ref{fig:cpu_speed_electrons_dions}
visibly shows that statistical fluctuations already dominate the bias of all variants of the \red{quasi-}geometric
integration method above this sampling threshold, despite the large ensemble size of
30000 particles.
To avoid possible sampling artifacts at \red{an }even \blue{higher }higher particle count, we consider the \red{quasi-}geometric orbit integration method at the toroidal grid size $N_\varphi$ of
at minimum twice the number of toroidal modes in the magnetic field configurations.
The variant with the polynomial series solution truncated at $K=2$ (\bluex{\textit{Gorilla Poly 2}}\redx{GORILLA  Poly 2}) at
this grid resolution can be considered the fastest sufficiently accurate tested method to compute
$D_{11}$ for thermal ions and electrons. In \red{the }case of D ions with an energy of 3~keV this method
is one order of magnitude faster than the Runge-Kutta~4 integrator which is the fastest
reference method.

\section{Summary and outlook}
\label{sec:summary}

A \red{quasi-}geometric integration method for guiding-center orbits in general three-dimensional toroidal fields has been developed, implemented and presented here. This orbit integration procedure is based on a representation of the electromagnetic field by continuous piecewise linear functions using a spatial mesh.
Collisionless particle orbits in real space and magnetic coordinates and their respective invariants of motion have been studied in detail for axisymmetric and non-axisymmetric geometries. Due to the special formulation of the guiding-center equations, the magnetic moment and the total energy are conserved naturally. In \red{the }case of toroidal axisymmetry the canonical toroidal angular momentum is accurately preserved by the \red{quasi-}geometric method \redx{in third and fourth order series expansion
over the orbit parameter}, as well as the parallel \redx{adiabatic} invariant.
Thus, the property of such systems to ideally confine the orbits is retained.
\redx{In order to evaluate the limitations of these confinement properties, the kinetic energy of ions was chosen to be 300~keV. Otherwise systematic ODE integration errors originating from truncating the series expansion already at the second order would not be visible, since this error is proportional to the particle energy.}

For passing orbits in 3D fields, however, the piecewise linearization of the electromagnetic field
introduces some artificial chaotic diffusion which, nevertheless, could be made negligibly small by spatial grid refinement.
For trapped orbits in 3D fields, the approximate conservation of the parallel adiabatic invariant is not violated by significant
error accumulation.

To assess the method's performance, the mono-energetic radial transport coefficient, $D_{11}$, which gives a main contribution to neoclassical transport, has been evaluated for a quasi-isodynamic reactor-scale stellarator field~\cite{drevlak_quasi-isodynamic_2014} using the Monte Carlo method. For both, electrons and ions, the results obtained by \red{quasi-}geometric orbit integration are in good agreement with the results of adaptive RK4/5 integration with usual spline interpolation of electromagnetic fields.
In the performance benchmark, we observe that the \red{quasi-}geometric orbit integration method with the polynomial series solution truncated
at $K=2$ is the fastest sufficiently accurate tested method to compute $D_{11}$. For the case of D ions with an energy of
3~keV the guiding-center orbit integration is one order of magnitude faster than 4th order Runge-Kutta integration in splined fields.
\redx{Here, truncating the series expansion at the second order is not necessarily contradicting the result of collisionless guiding-center orbits, where a second order truncation leads to visible systematic errors at high kinetic energies. An appropriate choice of grid size and series expansion order strongly depends on the physical application, the kinetic particle energy and the complexity of the magnetic field (e.g. number of harmonic modes).}

For the application in global kinetic computations no extra effort is needed to obtain dwell times within spatial grid cells as
these are computed automatically in the present approach. Additionally, integrals of velocity powers over these dwell times
are available as analytical expressions.
The latter quantities are required for statistical scoring of orbits in Monte
Carlo computations of macroscopic parameters,
such as plasma response currents and charges caused by external non-axisymmetric
perturbations in tokamaks or parameters of the edge plasma in devices with 3D field geometry.
Moreover, similarly to the geometric integrator for axisymmetric two-dimensional fields described in Ref.~\onlinecite{kasilov_geometric_2016}, the presented method is less sensitive to noise
in the electromagnetic field than procedures relying upon high order polynomial interpolation.
These characteristics suggest additional overall performance enhancements in both numerical stability and computational efficiency, when the \red{quasi-}geometric orbit integration is applied to kinetic modeling.

\red{The applicability of the inherently low-order method to particle-in-cell
turbulence computations, where higher order schemes produce smoother solutions,
is still an open question.}
Finally, it should be mentioned that the method facilitates the coupling with kinetic neutral particle codes \blue{like}\red{such as}
EIRENE~\citep{reiter_eirene_2005}, where one needs to model particle conversion into neutrals and back with plasma
and neutral particles described in different coordinate systems. \blue{Necessary}\red{The necessary}
transformation of coordinates does not require solving any implicit dependencies (nonlinear equations),
since that is a linear operation in this approach and \red{is }therefore intrinsically fast.

\section*{Data Availability Statement}
The data that support the findings of this study are available from the corresponding author upon reasonable request.

\section*{Acknowledgements}

The authors would like to thank Michael Drevlak, Markus Meisterhofer, Artem Savchenko and Harry Mynick for respectively providing stellarator field configurations, Fortran advice, a tool for non-linear coordinate transformations and a tutorial for neoclassical transport in 3D systems, as well as Martin Heyn, Philipp Ulbl, Rico Buchholz, Patrick Lainer and Markus Richter for useful discussions. \red{Further thanks to three anonymous referees who have read the manuscript and suggested valuable improvements. }This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programmes 2014–2018 and 2019–2020 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The study was supported by the Reduced Complexity Models grant number ZT-I-0010 funded by the Helmholtz Association of German Research Centers. Support from NAWI Graz, and from the OeAD under the WTZ grant agreement with Ukraine No. UA 04/2017 is gratefully acknowledged.

\newpage

redpajama_set_name: RedPajamaArXiv


text:
Forget to take your medication? A new digital pill will alert you — and your doctor.
(Courtesy of University of Minnesota Health Oncology Service Line)
Peter Holley
Technology reporter
January 17, 2019 at 6:00 AM EST
When his chemotherapy patients leave the hospital to continue treatment at home, Edward Greeno faces a new challenge.
He can no longer ensure they're taking their medicine.
Greeno, the medical director of the Masonic Cancer Clinic at the University of Minnesota, has come to realize that some patients, like children hiding naughty behavior from a parent, will fudge the truth to avoid his disapproval, even when their health is at risk.
To combat patients' fibbing and forgetfulness, Greeno has begun deploying a new tool in recent months: a pill embedded with a tiny, ingestible sensor. The sensor transmits data from inside the patient's body to a wearable patch placed on their abdomen, which then connects to a mobile app that patients and doctors can access.
That data offers a new window into patients' health and behavior, Greeno said, allowing doctors to remotely monitor someone's heart rate, activity level and sleep cycle. The sensor, which is about the size of a grain of sand and dissolves in the gastrointestinal tract, also tells doctors when a patient has ingested their medication. The information is compiled in a database that doctors can access from their devices.
At the Masonic Cancer Clinic, as part of a pilot program, doctors have begun pairing the sensor — made by a California company called Proteus Digital Health — with a common chemotherapy drug used to treat stage 3 and 4 colorectal cancer patients.
"When we give people chemotherapy in the clinic with an intravenous drug, we're able to assess the dose and timing and make sure they're well enough to continue getting the treatment," Greeno said. "But when you send them home with a bottle of pills, you don't know when they're taking them or if they're well enough to take them."
"You'd think for cancer treatment that patients would be pretty diligent, and that's not always the case, for a variety of reasons," he added.
Giant Food Stores will place robotic assistants at 172 locations, company says
Ingestible technology has raised a variety of concerns among experts and ethicists since the Food and Drug Administration approved digital pills in 2017, from how physicians will use the wave of new data created by the technology to whether that data may be vulnerable to hackers.
"The health industry is behind the curve when it comes to cybersecurity," Jason Christopher, chief technology officer at the cybersecurity company Axio, told Forbes last year. "Forget health record databases — how do you patch a digital pill?"
Despite the concerns, many experts have pointed out that the technology will help doctors ensure patients are taking their prescribed medication. Not taking medicine — or "medication nonadherence" in the health-care world — is a "common and costly problem," according to a study cited by the National Center for Biotechnology Information.
"Approximately 30% to 50% of US adults are not adherent to long-term medications leading to an estimated $100 billion in preventable costs annually," according to the 2013 study.
Comparing medication nonadherence to weight loss, NCBI researchers say the problem is more complex than it might seem and often has multiple complicating factors. Greeno agreed, pointing out that his patients aren't always lying about their behavior. Sometimes older patients have trouble opening pill bottles. Other times, he said, they're so fatigued and foggy that they struggle to adhere to a rigorous schedule that requires them to ingest multiple pills each day. And sometimes, he said, they struggle to deal with a medication's side effects.
The problem, Greeno said, is that chemotherapy drugs can have a very narrow therapeutic window and often cost between $10,000 and $20,000 per month.
"I had one patient whose hands hurt and she couldn't open her pill bottle," Greeno said, noting that when the patient's daughter was around she would take her pills, but when her daughter was gone, she wouldn't. "Normally, I might not find out about until weeks later when she comes to the clinic. But now, the app is telling me in real time that she didn't take her pills and I get that message in clinic the next morning."
Brenda Jans Darling, 45, has another reason she sometimes fails to take her chemotherapy medication on time: She's a working mother of two juggling recovery with a daily onslaught of activities.
Darling, a sixth-grade math teacher in Minnesota, said she initially agreed to partake in a digital medicine pilot program about three months ago because she thought it sounded interesting. What she quickly discovered, she said, was that she enjoyed seeing her health data like her sleep and activity levels tracked on an app. She also found that information useful, comparing it to a "medical Fitbit."
What she finds most useful, however, is knowing she no longer needs to count her medication to ensure she hasn't missed a dose at the end of a busy two-week medication cycle.
"I'll be having breakfast and one kid needs a permission slip signed and the other one needs something else and — while all that's going on — I'm trying to make sure I took the right number of pills at the right time of day," Darling said. "Life happens, right, and it can be tough to keep track of everything."
Now when she forgets, she added, she gets an alert on her iPhone that says, "We have not received data from your 8 a.m. dose," which she finds extremely useful.
Darling said she isn't worried about her doctors having access to an ever-growing amount of personal health information and believes that data will be used to combat the disease she's spent the last two years battling.
"I'm going to do whatever I can to advance cancer research so we can find a cure to our modern-day plague," she said.

redpajama_set_name: RedPajamaCommonCrawl


Shop best items at best prices at DORCO now! Go check it and save your money now! Take advantages of the DORCO online shop!
Shop with these top sale goods at DORCO. Fantastic voucher code that you can find on DealJonesy.com. Amazing savings, time is limited.
Receive the shocking benefit of our DORCO voucher code. See details on DealJonesy.com. Hurry before it's gone, it can not wait for anymore.
Spend much less on the wonderful products with our DORCO promo. Do not miss it at DORCO. The offer has limited time, it can not wait any longer.
Shop at DORCO and enjoy excellent savings! Great offers cannot wait for you! It can not last any longer at DORCO!
Wonderful products with wonderful price at DORCO promo. Best voucher begin and come to save. Excellent savings, do not miss it.
Enjoy the excellent benefit of our DORCO discount code. Purchase now for a discount price at DORCO. It only last for a limited time!
Shop with the products and get more savings from DORCO for a discount price. Great items at discounts. Do not miss them at DORCO.
Enjoy yourself to use this big discount code at DORCO. Offers are only for a limited time, use them quickly. It only last for a limited time!
Take extra money with this wonderful coupon provided by DORCO! Do not miss the amazing chance when you shop at DORCO. Check it out before it expires.
Excellent products with excellent price at DORCO promo code. Come and save now! A promo code like this, is hard to miss.
DORCO is offering a special sale today. Do not miss the special chance when you shop at DORCO. Take advantages of the DORCO online shop!
Receive wonderful benefit of this wonderful promo and shop with a discounted price from DORCO. You can enjoy your shopping at DealJonesy.com. Begin to save!
Receive wonderful reductions at DORCO. Shop the wonderful goods at DORCO. Exclusive promo are waiting for you.
Get the items you want at a great price with our DORCO promotional code. Test them now at DORCO! Exclusions may apply.
Don't pass this chance to get extra money with DORCO promotional code. Come and begin your shopping at DORCO! Be the first to use it!
A great voucher code from DORCO which you can't pass it! The more you shop at DORCO, the more you can save. Exclusions may apply.
Spend much less on the shocking products with our Orthotic Shop voucher. Do not miss the shocking chance when you shop at Orthotic Shop. The more you shop at Orthotic Shop the more savings are waiting for you!
Enjoy yourself to use this big discount code at Revant Optics. Offers are only for a limited time, use them quickly. It only last for a limited time!
Shop with this great promo from ice chips for a great price. Offer for a limited time only. You can benefit a lot from ice chips.
Get special reductions with our special Demco promo code. Voucher ends quickly! The promotional code is useful only for a limited time.
With our Citizen Goods promo code, you can get special discounts. Add the items to your cart and save instantly at Citizen Goods. Promotional code ends soon, make your order today!
It is a fantastic chance to receive extra money when you use Pixel Eyewear coupon code. These voucher code may no longer exist tomorrow. Limited time are waiting for you at Pixel Eyewear.
Enjoy the excellent benefit of our Better Eyewear discount code. Purchase now for a discount price at Better Eyewear. It only last for a limited time!
Receive big benefit of this big promo and shop with a discounted price from Bauer Nutrition. These deals are only last for a limited time. Do not wait for anymore.
You have come to the right place if you are looking for the newest and valid DORCO voucher code. Take 45% off with our lasted vouchers. Never spend more than you have to and use our DORCO coupon & sales. Simply click on one of the 14 available Dorco Free Shipping Code & enjoy shopping. Get everything you need with the new range now launched at DORCO!


text:
package savers

/****************
 * Date: 30/05/16
 * Time: 5:48 PM
 ***************/

import (
	"encoding/gob"
	"github.com/twinj/uuid"
	"log"
	"os"
	"path"
	"time"
)

var _ uuid.Saver = &FileSystemSaver{}

// This implements the Saver interface for UUIDs
type FileSystemSaver struct {
	// A file to save the state to
	// Used gob format on uuid.State entity
	file *os.File

	// Preferred location for the store
	Path string

	// Whether to log each save
	Report bool

	// The amount of time between each save call
	time.Duration

	// The next time to save
	uuid.Timestamp
}

func (o *FileSystemSaver) Save(pStore *uuid.Store) {

	if pStore.Timestamp >= o.Timestamp {
		err := o.openAndDo(o.encode, pStore)
		if err == nil {
			if o.Report {
				log.Printf("UUID Saved State Storage: %s", pStore)
			}
		}
		o.Timestamp = pStore.Add(o.Duration)
	}
}

func (o *FileSystemSaver) Read() (err error, store uuid.Store) {
	store = uuid.Store{}
	gob.Register(&uuid.Store{})

	if _, err = os.Stat(o.Path); os.IsNotExist(err) {
		dir, file := path.Split(o.Path)
		if dir == "" || dir == "/" {
			dir = os.TempDir()
		}
		o.Path = path.Join(dir, file)

		err = os.MkdirAll(dir, os.ModeDir|0755)
		if err == nil {
			// If new encode blank store
			err = o.openAndDo(o.encode, &store)
			if err == nil {
				log.Println("uuid.FileSystemSaver created", o.Path)
			}
		}
		log.Println("uuid.FileSystemSaver.Read: error will autogenerate", err)
		return
	}

	o.openAndDo(o.decode, &store)
	return
}

func (o *FileSystemSaver) openAndDo(fDo func(*uuid.Store), pStore *uuid.Store) (err error) {
	o.file, err = os.OpenFile(o.Path, os.O_RDWR|os.O_CREATE, os.ModeExclusive)
	defer o.file.Close()
	if err == nil {
		fDo(pStore)
	} else {
		log.Println("uuid.FileSystemSaver.openAndDo error:", err)
	}
	return
}

func (o *FileSystemSaver) encode(pStore *uuid.Store) {
	// ensure reader state is ready for use
	enc := gob.NewEncoder(o.file)
	// swallow error for encode as its only for cyclic pointers
	enc.Encode(pStore)
}

func (o *FileSystemSaver) decode(pStore *uuid.Store) {
	// ensure reader state is ready for use
	o.file.Seek(0, 0)
	dec := gob.NewDecoder(o.file)
	// swallow error for encode as its only for cyclic pointers
	dec.Decode(pStore)
}

redpajama_set_name: RedPajamaGithub


You will definitely want to be thinking sentimental for gift ideas!!
over the past couple of years!
photo taking of the two of them only happens every so often!
Rag & Bone Bindery has an endless amount of ways to preserve your precious memories!
It is a family run business out of Pawtucket, Rhode Island.
are glued and wrapped by skilled bookbinders.
using traditional bookbinding tools and techniques.
Rag & Bone Bindery is offering one reader a $50 gift certificate!!
Love these...mom is going to love it!
I made sure she got it before she saw it posted here, and her exact words in text were "Well you got me all crying. LOVE, LOVE the picture book with a bunch of smileys & hearts. Thank you so much!..." #bestdaughteraward :) haha j/k love you sis!
can you believe that we're talking Mother's Day already!? Ahh this year is flying by! Thanks for the cute idea, I love giving gifts for people but sometimes struggle with what to buy!!


text:
In the Spotlight: Tina Turner
Posted by Patrick Robbins at 10:44 am Add comments
In the Spotlight showcases a cross-section of an artist's cover work. View past installments, then post suggestions for future picks in the comments!
Tina Turner's place in rock 'n' roll history has long been assured. In fact, those words could have been spoken forty years ago, thanks to her Hall of Fame-worthy career with Ike Turner. But it's what she's done since then that really puts her over the top. She overcame a textbook case of a hellish marriage, turning tragedy into triumph. Her solo work has become her signature work, something no Beatle or Supreme could ever say. And for the past twenty years, she's been in the Guinness Book of Records for selling more concert tickets than any solo artist before or since. Make no mistake: her second trip into the Rock and Roll Hall of Fame is even more deserving than her first.
In Turner's honor, we've selected five covers of songs we associate with her. They've all been reinvented to some degree, and that's only appropriate, as Turner reinvented herself from a struggling supper-club singer to, simply, the best.
Rose Cousins – What's Love Got To Do With It (Tina Turner cover)
It took a while for "What's Love Got To Do With It" to find a home. Cliff Richard turned it down, Donna Summer never got around to it, and Turner herself didn't like the song at all, saying, "I can't sing that wimpy crap." Fortunately, she figured out how to keep it from being either, and took it to number one, her first and last song to top the Billboard charts. Maybe she was afraid that if she sang it in a gentler way she would sound the way Rose Cousins sounds in her cover. But as Cousins proves, the song sounds just fine in that style as well.
Peaches – Private Dancer (Tina Turner cover)
Mark Knopfler wrote "Private Dancer," and while Dire Straits recorded the track for their Love Over Gold album, Knopfler decided the song needed to be sung by a woman to work, and passed it along to Turner. It was the right move – "I put the old soulful touch on it," Turner said, and it became the title track and a top-ten hit. Peaches' cover takes it back to voice and piano, and while she didn't pull it off to her personal satisfaction, she still delivers with a vocal less sung than wrought.
Northern Kings – We Don't Need Another Hero (Tina Turner cover)
Of all the Finnish symphonic metal cover supergroups out there, Northern Kings was one of the best. Their first single, a cover of "We Don't Need Another Hero" from the third Mad Max movie, was a big hit in Finland, and the album it came from, 2007's Reborn, went gold with 17,000 copies sold. Sadly, they just couldn't get beyond Thunderdome, and by 2011 they were no more.
Joanne Kelly – Better Be Good To Me (Spider cover)
"Better Be Good To Me" is a song that could easily have been in our That's a Cover? category. It was originally released in 1981 by Spider, whose keyboardist Holly Knight co-wrote the song; the band folded in 1984, the same year Turner took the song into the top five. This cover comes from Going the Distance, a 2004 road-trip movie aimed at the adolescents of the Great White North. Canadian Pie, if you will. We see Joanne Kelly, one of the hitchhikers, tear through "Better Be Good To Me" in a roadside bar, and our protagonist realizes that hey, this girl might be something special.
Noah Reid – Simply the Best (Bonnie Tyler cover)
Here's another one where Turner made the world forget this was a cover, totally eclipsing the hard work of original artist Bonnie Tyler (Holly Knight co-wrote this one too). Her version of "The Best" is so huge that it's easy to miss the romance in the words and think it's about the winner of the big game. But Dan Levy didn't make that mistake, and when writing the "Open Mic" episode of Schitt's Creek, he entrusted costar Noah Reid with making the lyrics matter. Reid came through in spades, turning the volume way down and creating the most unforgettable moment of the show. The studio version can't be beat, either, as you'll see/hear.
In Memoriam: Otis Redding
Good, Better, Best: "Love is a Battlefield" (Pat Benatar)
Five Good Covers: "Keep On Running" (Spencer Davis Group/Jackie Edwards)
Feature, In the Spotlight Tagged with: Joanne Kelly, Noah Reid, Northern Kings, Peaches, Rock and Roll Hall of Fame, Rose Cousins, Tina Turner
Depeche Mode's Dave Gahan Shares a Soul-Stirring Cat Power Cover Susanna Hoffs Takes on Badfinger with Aimee Mann on New Covers Record


text:
No. 12 Wisconsin knocks off Loyola
By Steve Megaree
MADISON, Wis. — Wisconsin is at its best when hitting 3-pointers. That meant Loyola of Chicago picked the wrong night to play the Badgers.
Brad Davison and D'Mitrik Trice each scored 17 points as No. 12 Wisconsin capitalized on its sizzling outside shooting to beat Loyola of Chicago 77-63 on Tuesday.
Wisconsin (5-1) went 10 of 18 from beyond the arc and made four of those 3s during a 19-2 run that helped the Badgers build an 18-point advantage in the second half.
"We like to play with five guys on the floor who can all shoot the 3," said Wisconsin's Micah Potter, who went 3 of 4 from long range to finish with 13 points and seven rebounds. "That really spaces the floor, opens up driving lanes, opens up the post, all that kind of stuff."
The game was set up just two days earlier after both teams had scheduled Wednesday matchups scrapped due to the pandemic. Wisconsin had planned on hosting Northern Iowa, and Loyola was supposed to face Duquesne in Indianapolis.
Instead, these two midwestern programs squared off for the first time since 1998.
Jonathan Davis added 12 points as Wisconsin earned its 13th straight home victory. Cameron Krutwig scored 19 and Lucas Williamson had 15 for Loyola (3-1).
Krutwig and Williamson both played on Loyola's 2018 Final Four team. Ramblers coach Porter Moser said this game showed how his squad must build more depth to make a postseason run this season.
"I thought Lucas and Krutwig competed their tails off, really competed at a high level," Moser said. "We've got to get more of that."
Northwestern 100, Quincy 48. Miller Kopp scored 23 points and Northwestern made 14 3-pointers in a 100-48 victory over Quincy on Tuesday. Kopp scored seven of Northwestern's nine points over the final three minutes of the first half, and he scored 12 of their opening 16 points after the break for a 40-point lead. Boo Buie finished with 15 points, Chase Audige added 14 and Pete Nance 12 for Northwestern (3-1). Tanner Stuckman scored 10 points for Quincy.
Illinois State 91, Chicago State 62. DJ Horne had 22 points, Josiah Strong had 17 points, six rebounds and six assists for Illinois State (3-3). Dedric Boyd added 13 points and Antonio Reeves had 11 points. Xavier Johnson had 22 points for the Cougars (0-7). Carlo Marble added 17 points and Ke'Sean Davis had 12 rebounds.
Southern Illinois vs. UT Martin. Canceled.


Personal Views: We have come a long way in building a strong community, so let's stay the course to complete what we have started. My theme has always been "Think Bonney Lake First", and everything I have tried to do is with the best interests of the community in mind. Since becoming your Mayor in November 2005, we have been successful in securing $3.0 million in grants for the City.
During this time, I have worked with the council to create a long-term vision and an atmosphere of "team". The city continues to achieve stability among staff in all departments keeping costs down, which allows us to focus on long term planning, planning for growth, continued efficient and effective public safety, improvement of our local lakes, parks, and implementing the sidewalk, streetlight and downtown plans.
It is time we think of Bonney Lake first.


text: Q: UAC setting ask Credentials for for Admin Taks but not for signing Up I want to set UAC Setting on my Windows 10 PC such that when I log in through Admin User, It should not ask for Password but when I do Some Admin Stuff (Installing aplications, etc), It should give a prompt asking for Credentials.

redpajama_set_name: RedPajamaStackExchange


text:
Q: is there anything special i need to do on my server to run jQuery.getJSON on a bookmarklet? I am making a bookmarklet that requires some user details.
I understood my best options for cross domain communication is either jQuery.getJSON or appending a form and iframe to the document and setting the form target to the iframe.
I've decided to start with getJSON as it seems alot less messy however seems i have some problem with it.
Here is what i use for the json call :
$.getJSON('http://mydomain.com/bookmarklet_getjson_functions.php', function(data) {
        alert(data);
});

here is the php page:
echo "{";
echo "item1: ", json_encode($item1), "\n";
echo "item2: ", json_encode($item2), "\n";
echo "item3: ", json_encode($item3), "\n";
echo "}"

When i click on the bookmarklet link while Firebug is opened i get:
GET http://mydomain.com/bookmarklet_getjson_functions.php 200 OK 116ms
But the alert wont happened.
Any clue ?

A: Alright - found out what i did wrong after testing several other example - problem was i did not returned callback that jquery produce - it is used as sort of a "password"
So the solution is basically:
$.getJSON('http://mydomain.com/bookmarklet_getjson_functions.php?jsoncallback=?', function(data) {
    alert(data);
});

The last ? in ?jsoncallback=? get replaced automatically with a random function name that is later on used as "password"
so now in my php file i first echo the callback function then () and inside parentheses i echo the variable in json format:
$toecho = array(
"firstname" => "someone",
"lastname" => "somoene's lastname"

);
//return in JSON format
echo $_GET['jsoncallback'].'('.json_encode($toecho).')';

This will echo something like:
jquery_37328327837823232_3236363({"firstname":"somsomeonee","lastname":"omoene's lastname"})
Hope that helps someone out there.

A: That "php page" is a bit strange. You can try this:
header('Content-Type: text/javascript');
echo json_encode(array(
    'item1' => $item1,
    'item2' => $item2,
    'item3' => $item3,
));
exit;

You should set the content type for the response and properly encode the data you want to send.
(...testing the response is better with console.log() instead alert() :)