Guiding Center Dynamics in 4D
Guiding centre dynamics is a reduced version of charged particle dynamics, where the motion of the particle in a strong magnetic field $B$ is reduced to the motion of the guiding centre, that is the centre of the gyro motion of the particle about a magnetic field line. The dynamics of the guiding centre can be described in terms of only four coordinates (as compared to six for the full motion of the charged particle), the position of the guiding centre $r = (x,y,z)$ and the parallel velocity $u$, where parallel refers to the direction of the magnetic field.
Lagrangian Formulation
The simplest form of the guiding centre equations can be obtained from the Lagrangian
\[\begin{aligned} L &= \left( A (r) + u b (r) \right) \cdot \dot{r} - H(r,u) , & H &= \tfrac{1}{2} u^{2} + \mu \vert B (r) \vert , \end{aligned}\]
where $b = B / \vert B \vert$ is the unit vector of the magnetic field $B = \nabla \times A$ with $A$ the magnetic vector potential and $\mu$ is the magnetic moment. The Euler-Lagrange equations are computed as
\[\begin{aligned} \nabla \vartheta^{T} ( r(t), u(t) ) \cdot \dot{r} (t) - \dot{\vartheta} ( r(t), u(t) ) &= \nabla H ( r(t), u(t) ) , \\ b ( r(t) ) \cdot \dot{r} (t) &= u (t) , \end{aligned}\]
with $\vartheta(r,u) = A (r) + u \, b (r)$ and the gradient $\nabla$ denoting the derivative with respect to $r$.
Hamiltonian Formulation
Computing the time derivative of $\vartheta$, the Euler-Lagrange equations can be rewritten in an explicit form as
\[\begin{aligned} \dot{r} (t) &= \dfrac{u (t) \, \beta (r (t))}{b (r (t)) \cdot \beta (r (t))} + \dfrac{B (r (t))}{B (r (t)) \cdot \beta (r (t))} \times \nabla H (r (t),u (t)) , \\ \dot{u} (t) &= - \dfrac{\beta (r (t))}{b (r (t)) \cdot \beta (r (t))} \cdot \nabla H (r (t), u (t)) , \end{aligned}\]
where $\beta = \nabla \times \vartheta$. This constitutes a noncanonical Hamiltonian system of the form
\[\dot{q} = \Omega^{-T} (q) \nabla H(q) ,\]
with $q = (x,y,z,u)$ and the symplectic matrix $\Omega$ given by
\[\Omega_{ij} = \dfrac{\partial \vartheta_{i}}{\partial q^{j}} - \dfrac{\partial \vartheta_{j}}{\partial q^{i}} ,\]
which is the convention used throughout this package and shared with GeometricProblems.jl. With it the Euler-Lagrange equations of $L = \vartheta \cdot \dot{q} - H$ read $\Omega (q) \, \dot{q} = - \nabla H (q)$.
Modules
ChargedParticleDynamics.GuidingCenter4d — Module
The guiding centre dynamics in four dimensions: the gyro-averaged motion in the guiding centre position and the parallel velocity, $q = (X, u)$, with the magnetic moment $\mu$ a parameter.
A noncanonical, degenerate system built from the one-form $\vartheta = A + u b$. See the Guiding Center Dynamics in 4D page.
The equations are written once, in GuidingCenter4d, and read the field from params.field. They are documented below. Each equilibrium is its own module, holding its field as FIELD, its initial conditions and problem constructors that default to them, and it includes the diagnostics — plus guiding_center_4d_loop.jl and guiding_center_4d_surface.jl where it defines a Poincaré loop or surface. Those are documented once, under TokamakSmallCylindrical, and hold verbatim for all eleven. What differs between the modules is the chart, the equilibrium parameters and the initial conditions, which is what these docstrings record:
ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCylindrical — Module
Analytic axisymmetric small tokamak equilibrium in cylindrical coordinates.
ChargedParticleDynamics.GuidingCenter4d.SymmetricField — Module
First and second Poincaré invariant for a guiding center particle in an axisymmetric magnetic field of the form $B(x,y,z) = B_0 (1 + x^2 + y^2) e_z$.
The loop for the first invariant is initialized by
\[q (\tau) = \begin{pmatrix} r_x \cos (2\pi \tau) \\ r_y \sin (2\pi \tau) \\ z_0 + z_1 \sin (2\pi \tau) \\ u_0 + u_1 \cos (2\pi \tau) \\ \end{pmatrix}\]
with parameters
\[B_0 = 1, \quad r_x = 0.5, \quad r_y = 0.3, \quad z_0 = 0.0, \quad z_1 = 0.1, \quad u_0 = 0.5, \quad u_1 = 0.05, \quad \mu = 0.01 .\]
The surface for the second invariant is initialized by
\[q (\tau) = \begin{pmatrix} r_0 (\sigma - 0.5) \\ r_0 (\tau - 0.5) \\ z_0 + z_1 \cos (2\pi \sigma) \cos (2\pi \tau) \\ u_0 + u_1 \sin (2\pi \sigma) \sin (2\pi \tau) \\ \end{pmatrix}\]
with parameters
\[B_0 = 1, \quad r_0 = 0.5, \quad z_0 = 0.0, \quad z_1 = 0.1, \quad u_0 = 0.5, \quad u_1 = 0.01, \quad \mu = 0.01 .\]
ChargedParticleDynamics.GuidingCenter4d.ThetaPinchField — Module
First Poincaré invariant for a guiding center particle in an θ-pinch magnetic field of the form $B(x,y,z) = B_0 \, e_z$.
The loop for the first Poincaré invariant is initialized by
\[q (\tau) = \begin{pmatrix} r_x \cos (2\pi \tau) \\ y_0 \\ r_z \sin (2\pi \tau) \\ \end{pmatrix}\]
with parameters
\[B_0 = 1, \quad r_x = 0.5, \quad r_z = 0.3, \quad y_0 = 0.0, \quad u_0 = 0.5, \quad \mu = 2.5 \times 10^{-6} .\]
ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCartesian — Module
Analytic axisymmetric small tokamak equilibrium in cartesian coordinates.
ChargedParticleDynamics.GuidingCenter4d.TokamakSmallToroidal — Module
Analytic axisymmetric small tokamak equilibrium in circular coordinates.
ChargedParticleDynamics.GuidingCenter4d.TokamakMediumCartesian — Module
Analytic axisymmetric medium-size tokamak equilibrium in cartesian coordinates.
ChargedParticleDynamics.GuidingCenter4d.TokamakMediumCylindrical — Module
Analytic axisymmetric medium-size tokamak equilibrium in cylindrical coordinates.
ChargedParticleDynamics.GuidingCenter4d.TokamakIterCylindrical — Module
Analytic ITER-like Solov'ev equilibrium with X-point.
ChargedParticleDynamics.GuidingCenter4d.SolovevSymmetricField — Module
Analytic, quadratic Solov'ev equilibrium.
ChargedParticleDynamics.GuidingCenter4d.SolovevIter — Module
Analytic ITER-like Solov'ev equilibrium.
ChargedParticleDynamics.GuidingCenter4d.SolovevIterXpoint — Module
Analytic ITER-like Solov'ev equilibrium with X-point.
ChargedParticleDynamics.GuidingCenter4d.hamiltonian — Method
The guiding centre Hamiltonian,
\[H = \tfrac{1}{2} u^{2} + \mu \vert B \vert + \varphi ,\]
the parallel kinetic energy, the magnetic moment term and the electrostatic potential. The potential is the same one GuidingCenter3d carries; it is zero for every equilibrium this package ships, so adding it changed no result, but the model is now the one both reference papers write.
ChargedParticleDynamics.GuidingCenter4d.iodeproblem — Method
iodeproblem(qᵢ; kwargs...)
iodeproblem(ics::NamedTuple; kwargs...)The same dynamics as odeproblem in implicit form, as an IODEProblem built from the one-form $\vartheta = A + u b$, so that the variational and projection integrators can be applied to it. The initial momentum is $\vartheta(q_{i})$.
ChargedParticleDynamics.GuidingCenter4d.iodeproblem_λ — Method
iodeproblem_λ(qᵢ; kwargs...)iodeproblem with the Lagrange multiplier initialised explicitly, for the integrators that need a starting value for it rather than taking zero.
ChargedParticleDynamics.GuidingCenter4d.lodeproblem — Method
lodeproblem(qᵢ; kwargs...)
lodeproblem(ics::NamedTuple; kwargs...)The 4D guiding centre dynamics as an LODEProblem, from the degenerate phasespace Lagrangian $L = \vartheta \cdot \dot{q} - H$ together with the two-form $\omega$. Same equations as iodeproblem; the Lagrangian and the two-form are carried in addition, for the variational integrators that want them.
ChargedParticleDynamics.GuidingCenter4d.lodeproblem_formal_lagrangian — Method
lodeproblem_formal_lagrangian(qᵢ; kwargs...)This is currently a verbatim copy of lodeproblem — the same LODEProblem from the same ϑ, f, g, ω and lagrangian. It should instead be the formal Lagrangian of the equations of motion in noncanonical Hamiltonian form: writing them as $F_{i}(z, \dot{z}) = \Omega_{ij}(z) \dot{z}^{j} + \partial_{i} H(z) = 0$, the formal Lagrangian is $L(z, y, \dot{z}) = y^{i} F_{i}(z, \dot{z})$ on the doubled state $(z, y)$, whose Euler-Lagrange equations return the system on variation of $y$ and its adjoint on variation of $z$. That is an eight-dimensional problem, not the four-dimensional one built here. See TODO.md.
ChargedParticleDynamics.GuidingCenter4d.odeproblem — Method
odeproblem(qᵢ; kwargs...)
odeproblem(ics::NamedTuple; kwargs...)The 4D guiding centre dynamics as an explicit ODEProblem in the state $q = (X, u)$, the guiding centre position and the parallel velocity. The vector field is obtained from $\Omega \dot{q} = - \nabla H$.
The second form takes the named tuple that every initial_conditions_* of this module returns, whose params carry that condition's own magnetic moment μ.
ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCylindrical.compute_energy — Method
compute_energy(sol)The Hamiltonian along the solution, as a DataSeries.
ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCylindrical.compute_energy_error — Method
compute_energy_error(sol)Like every compute_*_error here, this forwards to GeometricSolutions.compute_invariant_error and so returns two series, not one: the invariant itself and its relative error. Destructure it — h, e = compute_energy_error(sol) — rather than passing the result somewhere a single series is expected.
ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCylindrical.compute_toroidal_momentum — Method
compute_toroidal_momentum(sol)The canonical toroidal momentum ϑ₃ along the solution, as a DataSeries. It is conserved in an axisymmetric equilibrium.
ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCylindrical.compute_toroidal_momentum_error — Method
compute_toroidal_momentum_error(sol)Returns a (value, error) pair; see compute_energy_error.
ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCylindrical.default_parameters — Method
The field, and the magnetic moment μ this equilibrium is set up for.
ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCylindrical.loop_ensemble — Method
loop_ensemble(prob, pinv)Sample the loop f_loop of this equilibrium with pinv and turn each point into the initial condition of one member of a GeometricEquations.EnsembleProblem built from prob. Integrate the result with GeometricIntegrators.integrate and pass the solution to PoincareInvariants.compute!.
ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCylindrical.poincare_invariant_1st — Method
poincare_invariant_1st(N; DT = Float64)Set up the computation of the first Poincaré integral invariant
\[I_{1} = \oint_{\gamma} \vartheta_{i} (q) \, dq^{i}\]
of the four-dimensional guiding centre dynamics, sampling the loop f_loop of this equilibrium at N points.
The guiding centre one-form $\vartheta = A + u b$ is state-dependent, so this is a noncanonical FirstPoincareInvariant built from ϑ, not the canonical variant. Pair it with one of loop_odeproblem, loop_iodeproblem or loop_lodeproblem through loop_ensemble:
pinv = poincare_invariant_1st(1000)
prob = loop_iodeproblem(; timespan = (0.0, 1E3), timestep = 1.0)
sol = integrate(loop_ensemble(prob, pinv), VPRKGauss(2))
I₁ = compute!(pinv, sol, parameters(prob))ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCylindrical.poincare_invariant_2nd — Method
poincare_invariant_2nd(N; DT = Float64, plan = SecondChebyshevPlan)Set up the computation of the second Poincaré integral invariant
\[I_{2} = \int_{S} \omega_{ij} (q) \, dq^{i} \wedge dq^{j}\]
of the four-dimensional guiding centre dynamics, sampling the surface f_surface of this equilibrium with the point specification N.
As for the first invariant, the guiding centre two-form $\omega = d\vartheta$ is state-dependent, so this is a noncanonical SecondPoincareInvariant built from ω. Pass N as a tuple (nx, ny) together with plan = SecondFinDiffPlan for a grid layout, which is what PoincareInvariants.plot_surface needs in order to draw the advected surface.
pinv = poincare_invariant_2nd(2000)
prob = surface_iodeproblem(; timespan = (0.0, 1E3), timestep = 1.0)
sol = integrate(surface_ensemble(prob, pinv), VPRKGauss(2))
I₂ = compute!(pinv, sol, parameters(prob))ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCylindrical.surface_ensemble — Method
surface_ensemble(prob, pinv)Sample the surface f_surface of this equilibrium with pinv and turn each point into the initial condition of one member of a GeometricEquations.EnsembleProblem built from prob. Integrate the result with GeometricIntegrators.integrate and pass the solution to PoincareInvariants.compute!.