gcmotion.FrequencyAnalysis#

The FrequencyAnalysis class iterates through (μ, Pζ, Ε) values upon a given Tokamak, and finds the ωθ, ωζ frequencies and their ratio qkinetic by searching for contours.

class gcmotion.FrequencyAnalysis(tokamak: Tokamak, psilim: tuple, muspan: ndarray, Pzetaspan: ndarray, Espan: ndarray = None, **kwargs)#

Performs a Frequency Analysis on a given Tokamak, by calculating contours.

Parameters:
tokamak: Tokamak

The tokamak to perform the analysis upon.

psilimtuple(float, float)

The \(\psi\) limit to restrict the search for contours, relative to \(\psi_{wall}\). Avoid setting the lower limit exactly 0, but set it at a very low value (e.g. 1e-7) to avoid RuntimeWarnings when ψ becomes negative due to rounding errors.

muspannumpy.ndarray

The \(\mu\) span. Can be either 1D or 2D. See documentation for definitions

Pzetaspannumpy.ndarray

The \(P_\zeta\) span. Can be either 1D or 2D. See documentation for definitions

Espannumpy.ndarray

The Energy span. Can be either 1D or 2D. See documentation for definitions

Methods

scatter(x, y[, scatter_kw])

Draws a scatter plot of the calculated frequencies.

start([pbar])

Iterates through the triplets and finds the frequencies of the resulting found contours.

to_dataframe([extended])

Creates a pandas DataFrame with the resulting frequencies.
Other Parameters:
main_grid_density: int, optional

The Main Contour’s’ grid density. For analytical equilibria, a minimum of 500 is recommended, while numerical equilibria require at least 1000. This option mostly affects trapped orbits around O-points. Diminishing results are observed after 1800. Defaults to 1000.

local_grid_density: int, optional

The local contours’ grid density. 100 seems to work best for almost all cases. Defaults to 100.

theta_expansion, psi_expansion: float, optional

The bounding box expansion factor around which to create the local contour. Note that the expanded bounding box is still limited by the Main Contour’s limits. Default to 1.2.

pzeta_min_step: float, optional

The minimum \(P_\zeta\) step allowed. This step is usually smaller than the relative minimum step, but takes over for \(P_\zeta \approx 0\). Defaults to 2e-4.

passing_pzeta_rstep: float, optional

The relative \(P_\zeta\) step size for passing orbits. Defaults to 1e-3.

trapped_pzeta_rstep: float, optional

The relative \(P_\zeta\) step size for trapped orbits. Defaults to 1e-3.

energy_rstep: float, optional

The relative E step size. Defaults to 1e-3.

skip_trapped: bool, optional

Whether or not to skip calculations of trapped orbits. Defaults to False.

skip_passing: bool, optional

Whether or not to skip calculations of passing orbits. Defaults to False.

skip_copassing: bool, optional

Whether or not to skip calculations of co-passing orbits. Defaults to False.

skip_cupassing: bool, optional

Whether or not to skip calculations of cu-passing orbits. Defaults to False.

calculate_omega_theta: bool, optional

Whether or not to calculate each orbit’s \(\omega_\theta\). Defaults to True.

calculate_qkinetic: bool, optional

Whether or not to calculate each orbit’s \(q_{kinetic}\). Defaults to True.

cocu_classification: bool, optional

Whether or not to further classify passing orbits as Co-/Cu-passing. Defaults to True.

min_vertices_method_switch: int, optional

Minimum number of vertices below which the frequency calculations switch to double contour mode. Defaults to 40.

max_pzeta_method_switch: float, optional

Maximum \(P_\zeta\) (absolute) value, below which the frequency calculations switch to double contour mode. Defaults to 0.

relative_upper_E_factor: float, optional

Used only in dynamic minimum energy mode. Defines the upper Espan limit to search for trapped orbits, relative to the total minimum Energy. Defaults to 1.1.

logspace_len: int, optional

Used only in dynamic minimum energy mode. The number of Energies for which to search trapped orbits, up until relative_upper_E_factor. Defaults to 50.

trapped_min_num: int, optional

Used only in dynamic minimum energy mode. Defines the minimum number of trapped orbits, after which calculation stops. Defaults to 1.

Examples

>>> import numpy as np
>>> import gcmotion as gcm
>>>
>>> Rnum = 1.65
>>> anum = 0.5
>>> B0num = 1
>>> species = "p"
>>> Q = gcm.QuantityConstructor(R=Rnum, a=anum, B0=B0num, species=species)
>>>
>>> B0 = Q(B0num, "Tesla")
>>> i = Q(0, "NUPlasma_current")
>>> g = Q(1, "NUPlasma_current")
>>>
>>> tokamak = gcm.Tokamak(
...     R=Q(Rnum, "meters"),
...     a=Q(anum, "meters"),
...     qfactor=gcm.qfactor.Unity(),
...     bfield=gcm.bfield.LAR(B0=B0, i=i, g=g),
...     efield=gcm.efield.Nofield(),
... )
>>>
>>> # Cartesian mode span
>>> muspan = np.array([1e-5])
>>> Pzetaspan = np.linspace(-0.04, -0.005, 10)
>>> Espan = np.logspace(7.8e-6, 4e-5, 20)
>>>
>>> freq = gcm.FrequencyAnalysis(
...     tokamak=tokamak,
...     psilim=(1e-7, 1.2),
...     muspan=muspan,
...     Pzetaspan=Pzetaspan,
...     Espan=Espan,
...     main_contour_density=800,
... )
>>>
>>> freq.start()
>>> print(freq)
>>> freq.scatter(x="Pzeta", y="qkinetic")
>>> df = freq.to_dataframe()
start(pbar: bool = True)#

Iterates through the triplets and finds the frequencies of the resulting found contours.

Parameters:
pbar: bool, optional

Whether or not to display a progress bar. Defaults to True.

to_dataframe(extended: bool = False)#

Creates a pandas DataFrame with the resulting frequencies.

Parameters:
extendedbool

Whether or not to add extra information for every orbit.

scatter(x: str, y: str, scatter_kw: dict = {})#

Draws a scatter plot of the calculated frequencies.

Parameters:
x: {“mu”, “Pzeta”, “E”, “qkinetic”, “omega_theta”, “omega_zeta”}

The x coordinate, as a column name of the Dataframe.

y: {“mu”, “Pzeta”, “E”, “qkinetic”, “omega_theta”, “omega_zeta”}

The y coordinate, as a column name of the Dataframe.

scatter_kw: dict, optional

Further arguements to be passed to matplotlib’s scatter method. Defaults to {}.

Notes on Input shapes#

The algorithm supports 3 modes, which are activated depending on the shape of the input arrays:

  1. Cartesian Mode: Activated by passing 3 1D arrays.

    The algorithm takes all combinations (cartesian product) of every array entry. If you want to iterate through only 1 COM value, use np.array([<value>]).

  2. Matrix Mode: Activated by passing 3 2D arrays, with the same shape.

    The algorithm creates triplets of COMs by stacking the 3 arrays. Each triplet is defined as (muspan[i,j], Pzetaspan[i,j], Espan[i,j]), where 0<=i<nrows and 0<=j<ncols. Useful when the grid is not orthogonal, for example when analysing a certain \(P_\zeta - E\) domain of the parabolas.

    Note: If we know that our grid is orthogonal, it’s much faster to use a row from each span array and use cartesian mode, since its significantly faster, especially when iterating through a lot of energies.

  3. Dynamic minimum energy Mode: Activated by only passing muspan and Pzetspan as 1D arrays.

    The algorithm finds the minimum vaule of the energy grid for every (\(\mu\), \(P_\zeta\)) pair, which is always found at an O-point, and slowly increments it until it finds 1 trapped orbit. This orbit’s frequency is (very close to) the O-point frequency, which is always the highest frequency of this specific family of trapped orbits. This frequency defines the maximum frequency with which the particles resonate, with 0 being the minimum (separatrix). This mode can only find the O-point frequency on the total minumum energy. If more O-points are present, then we must use method 4.

Todo

  1. Dynamic O-point minimum energy Mode:

Algorithm#

Each contour represents a specific family of orbits represented by the same 3 Constants of Motion, and differing only in their initial conditions (however, a single triplet may correspond to more that 1 contour). By exploiting the fact that our poloidal angle is in fact the boozer theta, the area contained within the contour is equal to \(2πJ_\theta\), where \(J_\theta\) the corresponding action variable. We then use the definitions:

\[\begin{split}\omega_\theta = \dfrac{dE}{dJ_\theta} \\ q_{kin} = - \dfrac{dJ_\theta}{dJ_\zeta} = -\dfrac{dJ_\theta}{dP_\zeta} \\ \omega_\zeta = q_{kin} * \omega_\theta \\\end{split}\]

The algorithm follows these steps, regardless of the “scanning” method (cartesian, matrix, dymanic energy minimum, …), since the only thing they change is the way the triplets (μ, Ε, Ρζ) are created:

The area is found using the shoelace formula

  1. For a triplet, iterate first through all the energies, then Ρζ and then μ:

    >>> for mu in muspan:
>>>     ...
>>>     for pzeta in pzetaspan:
>>>         ...
>>>         create_main_contour()
>>>         ...
>>>         for energy in energyspan:
>>>             ...

    Since all orbits with the same μ and Ρζ share the same slice, we can calculate the Main Contour once and use it for all energies. As for the μ and Ρζ loops, they only update the Profile’s μ and Ρζ and are completely symmetrical.

    The Main Contour is always plotted with θ from -2π to 2π, while ψ is defined by the user. This ensures that all orbits are present, even trapped orbits around θ=π. These orbits also appear twice, which leads to calculating the same thing two times, but the performance impact is not worth the added complexity needed to avoid it…

    For matrix mode, we must create the main contour for every energy as well, since the grid is generally not orthogonal.

    Note that we do not use matplotlib’s contour methods, but contourpy’s ContourGenerator. Not only this is much faster and more memory efficient, but it also gives more control when extracting contour lines.

2. For every energy now, we extract all contour lines the generator found for that energy. The number of contour lines can be 0 up to 4-5, depending on the equilibrium and energy level.

The contour lines are returned as a list of (N, 2) numpy arrays.

3. For every line we create a ContourOrbit object. This object represents a single family orbit. It contains methods for validating itself and classifying the orbit type, and also stores all calculated quantities and frequencies.

4. We calculate the bounding box of the orbit (e.g. the smallest rectangle fully containing the orbit who’s sides are parallel to the x and y axes).

5. We classify the orbit as trapped or passing. This is immediately calculated by checking whether the bounding box touches both left and right walls (passing) or not (trapped).

If the orbit is passing, we further classify it as co- or cu-passing. For this we use the approximation that if rho>0 on all vertices then the orbit is co-passing, else cu-passing.

  1. We validate the orbit, by checking 2 things:

    • The orbit is fully contained within the contour limits and doesn’t get cutoff. This is True if its bounding box does not touch any of the upper or lower contour walls.

    • The orbit is not cutoff-trapped, e.g. a trapped orbit that gets cut off by the left and right contour walls. The orbit is cutoff-trapped if its bounding box touches the left or the right wall, but not both.

    These 2 simple conditions are enough to fully validate the orbit. If one of those checks fails, we discard the orbit.

6. Calculate the orbit’s frequencies. The idea is to generate 2 adjacent local contours by slightly increasing and decreasing one of the constants Pζ or E and calculating the same orbit in these slightly different contours. We can then calculate their Jθs. Their difference would be dJθ, while dPζ or dE is defined by the user. Using the above definitions, we can calculate each frequency using these values. This is basically calculating the derivative locally.

  • Since the Main Contour is calculated in the θ-ψ plane, we need to convert the ContourOrbit’s vertices’ y coordinate to Pθ to correctly calculate Jθ. No need to convert the whole ψ-grid to Ρθ-grid each time.

  • If the main orbit has enough vertices as to calculate its Jθ accurately enough, we can avoid creating a 2nd contour by using the main orbit itself as one of the adjacent contours. This is what happens with most orbits. However, for low-energy trapped orbits, the orbit dimensions are comparable to the contour’s grid spacing, making the contour line jagged. By creating 2 adjacent local contours we circumvent this problem, since the orbit in the adjacent contour has much more vertices, and the main orbit’s Jθ is not calculated at all.

  • If at any point of the calculation an invalid contour line is found, the calculation is aborted. The number of orbits lost by this restriction however is negligible.

Advantages#

1. The algorithm is very close to the theoretical formulation of the Action-Angle theory, and uses a purely geometrical way of calculating frequencies, without the need for particle tracking.

2. There is no restriction as to how close or how far the array spans need to be spaced, since each triplet’s frequencies are calculated independently. This makes it possible to hold one or two of the COMs constant and iterate through a range of the the third COM. Moreover, we can use a more dense array for areas we expect more trapped orbits and a less dense one for areas with more passing orbits, since the Hamiltonian tends to be “flatter” in areas around O-points.

Caveats#

1. The approximation of rho in the co-/cu-passing classification is usually correct, except for some orbits exceptionally close to the separatrix. At the presence of more O-points, new families of passing orbits may be created, which cannot be classified this way. Those orbits have the flag “undefined”.

2. When calculating the Jθ for a passing orbit, we must add the points [-2π, 0] and [0, 2π] to its vertices. Since the order (left-to-right or right-to-left) of the vertices returned by the ContourGenerator is effectively random (but is always one of the two), we must first check their direction and then add the points in the right order. We must also divide by 2, since the contour is calculated in [-2π, 2π].

3. When calculating same-level orbits in the adjacent contours, the ContourGenerator may return more than 1 orbits, in a random order. This happens more often than not, since for most passing orbits, there is a corresponding copassing orbit with the same COMs. It can also happen with trapped orbits when more than 1 O-points are present. This problem can be solved by comparing the distances of the orbits’ bounding boxes (which are already calculated at this point) from the main orbit’s bounding box, and picking the closest one. Specifically, we compare the distances of the bottom left corner.

  1. Does not work under the presence of perturbations.

ContourOrbit Object#

class gcmotion.scripts.frequency_analysis.contour_orbit.ContourOrbit(E: float, vertices: ndarray)#

Path-like object containing the vertices as well as calculated frequencies, flags and methods needed to classify the orbit.

The methods should be called in a specific order, which is done inside profile_triplet_analysis() since some extra parameters are needed.

Methods

ContourOrbit.validate(psilim)

Checks if the bbox of the contour line touches the upper or lower walls, which means the orbit gets cut off and must be discarded.

ContourOrbit.calculate_bbox()

Calculates the orbit's bounding box (smallest rectangle fully containing the orbit).

ContourOrbit.distance_from(bbox)

Returns a distance-like quantity of the origin point from self's origin point.

ContourOrbit.classify_as_tp()

Classifies orbit as trapped or passing.

ContourOrbit.close_segment()

If the segment is passing, figure out if the points are left-to-right or right-to-left and append the two bottom points coorectly.

ContourOrbit.convert_to_ptheta(findPtheta, Q)

Converts all ycoords of the vertices from ψ to Pθ.

ContourOrbit.calculate_Jtheta()

Calculates the action J.

ContourOrbit.classify_as_cocu(profile)

Classifies orbit as co-/counter-passing.

ContourOrbit.pick_color()

Sets the segment's color depending on its orbit type.

ContourOrbit.str_dump()

Generates a small string corresponding to the orbit's type.