gcmotion.Particle

class gcmotion.Particle(tokamak: Tokamak, init: InitialConditions)

Creates a specific Particle in a specific Profile and with specific InitialConditions.

A particle entity represents a fully-fledged particle inside a specific tokamak device, and defined initial conditions.

Parameters:

tokamakTokamak

Tokamak object containing information about the tokamak.

initInitialConditions

InitialConditions object containing the set of initial condtions of the particle.

Methods

quantities([which, everything])	Prints the pint Quantities of the object.
run([method, stop_after, events, info])	Calculates the particle's orbit.

Notes

To view the particle’s attributes, use its quantities() method.

Examples

Creating a Particle:

>>> import gcmotion as gcm
>>> import numpy as np
>>>
>>> # Quantity Constructor
>>> Rnum = 1.65
>>> anum = 0.5
>>> B0num = 1
>>> species = "p"
>>> Q = gcm.QuantityConstructor(R=Rnum, a=anum, B0=B0num, species=species)
>>>
>>> # Intermediate Quantities
>>> R = Q(Rnum, "meters")
>>> a = Q(anum, "meters")
>>> B0 = Q(B0num, "Tesla")
>>> i = Q(10, "NUPlasma_current")
>>> g = Q(1, "NUPlasma_current")
>>> Ea = Q(73500, "Volts/meter")
>>>
>>> # Construct a Tokamak
>>> tokamak = gcm.Tokamak(
...     R=R,
...     a=a,
...     qfactor=gcm.qfactor.Hypergeometric(a, B0, q0=1.1, q_wall=3.8, n=2),
...     bfield=gcm.bfield.LAR(B0=B0, i=i, g=g),
...     efield=gcm.efield.Radial(a, Ea, B0, peak=0.98, rw=1 / 50),
... )
>>>
>>> # Setup Initial Conditions for "RK45" method
>>> init = gcm.InitialConditions(
...     species="p",
...     muB=Q(0.5, "keV"),
...     Pzeta0=Q(-0.015, "NUCanonical_momentum"),
...     theta0=0,
...     zeta0=0,
...     psi0=Q(0.6, "psi_wall"),
...     t_eval=Q(np.linspace(0, 1e-3, 1000), "seconds"),
... )
>>>
>>> # Create the particle
>>> particle = gcm.Particle(tokamak=tokamak, init=init)

Running the particle with the RK45 method:

>>> particle.run()
>>>
>>> events = [gcm.events.when_theta(init.theta0, 3)]
>>> particle.run(events=events)

Running the particle with the NPeriods method:

>>> particle.run(method="NPeriods", stop_after=1)
>>> particle.run(method="NPeriods", stop_after=4)

Note

The InitialCondition’s t_eval arguement is ignored when the method NPeriods is used.

quantities(which: str = '', everything: bool = False)

Prints the pint Quantities of the object.

Parameters:

whichstr, optional

Options on which Quantities to print. Can include many options.

1. “init” :

   Prints all the Quantities defined upon the particle’s instanciation.

2. “NU” or “SI”:

   Print the respective subset of Quantities

Options can be chained together, for example “initNU”. Defaults to “” (prints all Quantites)

everythingbool, optional

Whether or not to print all particle’s attributes, and not just Quantities. Ignores “which” arguement if True. Defaults to False.

run(method: str = 'RK45', stop_after: int = None, events: list = [], info: bool = False)

Calculates the particle’s orbit. The results are stored in both SI and NU.

Parameters:

method: {“NPeriods”, “RK45”}, optional

The solving method, passed to SciPy’s solve_ivp. Options are “RK45” (Runge-Kutta 4th order), “NPeriods” (Stops the orbit after N full \(\psi-P_\theta\) periods). Events are ignored if the method is “NPeriods”. Defaults to “RK45”

stop_after: int, optional

After how many full periods to stop the solving, if the method is “NPeriods”. Ignored if the method is “RK45”. Defaults to None.

eventslist, optional

The list of events to be passed to the solver, if the method is “RK45”. Ignored if the method is “NPeriods”. Defaults to [].

infobool, optional

Whether or not to print an output message. Defaults to False.

Notes

See Particle documentation for examples

Solvers

2 Solvers are available for calculating a particle’s orbit, and each can be useful in different circumstances.

Runge-Kutta 5(4)

This method is SciPy’s RK45 method. It is the default method when calling Particle.run().

The Particles InitialConditions must have t_eval defined, since RK45 integrates over a specific timespan.

Events are also available for use with this solver, and more than one can be used at a time. The terminal parameter can be used to halt the integration after the event was trigger terminal times.

Important

Even though multiple events with the same or different terminal can be used, they trigger completely intdependantly from one another. This means that we can’t use events to halt the integration at a time where more than one events are triggered, and currently SciPy has no support for such method.

NPeriodSolver

The NPeriodSolver is a custom-made Solver basing SciPy’s OdeSolver, and was created to circumvent the above problem. This is extremely useful, because it gives us the ability to stop the integration after the condition

\[\theta = \theta_0 \quad \text{and} \quad \psi = \psi_0\]

has been met N times, which amounts to exactly N completed orbit periods. This is derived from our Hamiltonian theory.

Note

Even though the actual conjugate variables are \(\psi\) and \(P_\theta\), the condition for a full period still holds if and only if the \(\psi-P_\theta\) relation is 1-1, which is always true since \(\partial P_\theta/\partial\psi> 0\)

Solver

class gcmotion.scripts.orbits.nperiod_solver.NPeriodSolver(fun: Callable, t0: float, y0: ndarray, t_bound: float, first_step=None, rtol: float = 1e-12, atol: float = 1e-12, max_step: float = 40, **extraneous)

Bases: OdeSolver

Custom OdeSolver modeled after SciPy’s RK45 Solver.

It is essentially identical to SciPy’s RK45 method, with a hard-coded event, and an extra step: For every solving step, the event triggers when the corresponding variable returns to its initial value. It then checks if the current value of its conjugate variable is also close to its initial value. If so, the Solver treats this as a full period of the orbit.

Depending on the value of the stop_after arguement, the integration can halt after an arbitrary number of full periods.

The coefficients belong to the Dormand-Prince method, which is a member of the Runge-Kutta family.

Methods

last_step_recursion()	Take recursing steps with smaller and smaller steps to reach the exact end of the integration.
period_check()	Checks if the orbit is close to completing a full period, and calls last_step_recursion() to halt the integration after stop_after periods.

References

[1]

J. R. Dormand, P. J. Prince, “A family of embedded Runge-Kutta formulae”, Journal of Computational and Applied Mathematics, Vol. 6, No. 1, pp. 19-26, 1980.

[2]

L. W. Shampine, “Some Practical Runge-Kutta Formulas”, Mathematics of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.

_step_impl()

Identical to RK45, with an extra call to period_check().

period_check()

Checks if the orbit is close to completing a full period, and calls last_step_recursion() to halt the integration after stop_after periods.

This method is called inside every step() iteration, but only proceeds if the integration is ‘close’ to a trigger point, e.g. the event variable is closing in on its initial value, and calculates the exact time this happens by a root finding algorithm.

It then checks if the conjugate variable is also closing in on its initial value. If so, the current state of integration, (t, y) is stored as an event (in contrast to how normal SciPy events work). The periods_completed is incremented by 1.

If the periods_completed attribute has reach the stop_after parameter, it means that we are just before or just after the point we want to halt the integration, depending on where the last step landed with respect to the event point. If the solver overstepped, we take one step back and call last_step_recursion().

last_step_recursion()

Take recursing steps with smaller and smaller steps to reach the exact end of the integration.

This method is called only when the full number of periods have been completed and the solving is reaching its end. At this point, the period_check() method is deactivated, and step() works identically to the RK45 method, except that the step_size is pre-set.

It recursively calls step() with a step_size equal to half the distance between of the last point and the event point. This asymptotically leads the orbit to its closing point, since we basically keep taking smaller and smaller steps to our target without overshooting.

Note that this recursion happens within the solver itself, in contrast to step() which is normally called only in solve_ivp(). This means that the solver’s dense_output is not updated, so we set the Particle’s y_final attribute to the last values we calculated. Those are later appended to the total calculated timeseries.