XGC-Analysis package

xgc_analysis.AnalyticDiffusionProfiles module

class xgc_analysis.AnalyticDiffusionProfiles.AnalyticDiffusionProfiles(psi, time, S, Vol, grad_avg, psi_norm)

Bases: object

Creates diffusion profile data using analytic functions.

Attributes:

density : A 3D array (shape (n_species, n_samples, n_surf)) of the diffusion density. flow : A 3D array (shape (n_species, n_samples, n_surf) of the parallel mean flow. temp : A 3D array (shape (n_species, n_samples, n_surf)) of temperature (in eV). n_samples : An integer (n_samples) representing the number of time samples. n_species : An integer (default 2) representing the number of particle species. n_surf : An integer (n_surf) representing the number of ψ-surfaces. psi : A 1D array (shape (n_surf,)) containing the poloidal magnetic flux values. steps : A 1D array (shape (n_samples,)) of time step indices. time : A 1D array (shape (n_samples,)) of simulation times (in seconds).

xgc_analysis.DiffusionCoefficients module

class xgc_analysis.DiffusionCoefficients.DiffusionCoefficients(nspecies=2, npsi=100)

Bases: object

Stores the profiles of diffusion model coefficients defined as functions of the poloidal magnetic flux ψ.

The four diffusion coefficients are:

• ptl_diffusivity (m^2/s)

• momentum_diffusivity (m^2/s)

• heat_conductivity (m^2/s)

• ptl_pinch_velocity (m/s)

The data layout of each coefficient is (nspecies, npsi) where the second dimension (the ψ-grid) is contiguous in memory. Also the ψ–grid is stored as a 1D NumPy array of length npsi.

Parameters:

nspecies: Number of particle species (default: 2). npsi: Number of points in the ψ–grid.

write_to_file(filename='xgc.diffusion_coeff.bp')

Writes the diffusion coefficient profiles and the ψ–grid to an Adios BP file as a new step.

For each diffusion coefficient, the data for each species is split into an individual variable using a species suffix. The species suffixes (for species indices 0 through 6) are:

0: _elec 1: _ion 2: _imp1 3: _imp2 4: _imp3 5: _imp4 6: _imp5

In addition, the file will store:

• n_species: a scalar integer (int64).

• psi: a 1D array (of length npsi).

The method opens the specified file (filename) using the Adios2 FileWriter in append mode, writes all the individual diffusion profiles, and then ends the step.

xgc_analysis.DiffusionProfiles module

class xgc_analysis.DiffusionProfiles.DiffusionProfiles(filename='xgc.diffusion_profiles.bp')

Bases: object

Reads diffusion profile data from an Adios BP file (directory) “xgc.diffusion_profiles.bp”.

The file has the following structure:

double density n_steps*{n_species, n_samples, n_surf} double flow n_steps*{n_species, n_samples, n_surf} double temp n_steps*{n_species, n_samples, n_surf} int32_t n_samples n_steps*scalar int32_t n_species n_steps*scalar int32_t n_surf n_steps*scalar double psi n_steps*{n_surf} int32_t steps n_steps*{n_samples} double time n_steps*{n_samples}

That is, there are (for example) n_steps steps recorded in the BP file. The constructor of this class reads only the last step (i.e. the most recent one).

Attributes:

density : A 3D array (shape (n_species, n_samples, n_surf)) of the diffusion density. flow : A 3D array (shape (n_species, n_samples, n_surf) of the parallel mean flow. temp : A 3D array (shape (n_species, n_samples, n_surf)) of temperature (in eV). n_samples : An integer (n_samples) representing the number of time samples. n_species : An integer (default 2) representing the number of particle species. n_surf : An integer (n_surf) representing the number of ψ-surfaces. psi : A 1D array (shape (n_surf,)) containing the poloidal magnetic flux values. steps : A 1D array (shape (n_samples,)) of time step indices. time : A 1D array (shape (n_samples,)) of simulation times (in seconds).

xgc_analysis.constants module

Physical constants and CGS-to-SI conversion factors for plasma physics.

Values based on the 2019 NRL Plasma Formulary.

All constants are in SI units unless otherwise noted.

xgc_analysis.csr_matrix module

class xgc_analysis.csr_matrix.SparseMatrixCSR(reader_function)

Bases: object

A class for storing a sparse matrix in CSR format (using SciPy’s csr_matrix).

The constructor takes a reader function as input. This reader function should return the matrix data as:

(ncols, nrows, nelement, width, col_indices, values)

where:

• ncols (int): Number of columns.

• nrows (int): Number of rows.

• nelement (np.ndarray): A 1D array (length nrows) with the number of nonzeros in each row.

• width (int): The maximal number of nonzeros in any row.

• col_indices (np.ndarray): A 2D array (nrows x width) with 0-based column indices.

• values (np.ndarray): A 2D array (nrows x width) with the nonzero values.

The constructor reformats this data and stores the matrix in SciPy’s CSR format.

get_csr_matrix()

Returns:

csr_matrix: The sparse matrix in CSR format.

multiply_left(vector)

Multiplies the CSR matrix from the left by a 1D NumPy array (i.e. computes vector * A).

Args:

vector (np.ndarray): 1D array with length equal to the number of rows of the matrix.

Returns:

np.ndarray: 1D array of length equal to the number of columns of the matrix.

Raises:

ValueError: if the length of vector does not match the number of rows.

multiply_transpose_left(vector)

Multiplies the transpose of the CSR matrix from the left by a 1D NumPy array (i.e. computes vector * A^T).

Args:

vector (np.ndarray): 1D array with length equal to the number of columns of the matrix.

Returns:

np.ndarray: 1D array of length equal to the number of rows of the matrix.

Raises:

ValueError: if the length of vector does not match the number of columns.

xgc_analysis.csr_matrix.sparse_matrix_reader(filename, var_names)

Read the sparse‑matrix variables stored in an ADIOS2 BP directory.

Parameters:

• filename (str) – Path to the *.bp directory.

• var_names (list[str]) –

  The six variable names in this exact order:

  • ncols_name – number of columns

  • nrows_name – number of rows

  • width_name – max. non‑zeros per row

  • nelement_name – 1‑D array of non‑zeros per row

  • col_indices_name – 2‑D array of 1‑based column indices

  • values_name – 2‑D array of non‑zero values

Returns:

(ncols, nrows, nelement, width, col_indices, values), where col_indices has been converted to 0‑based indices.

Return type:

tuple

xgc_analysis.field_data module

class xgc_analysis.field_data.FieldData(mesh, work_dir='./', file_indices=None, is_axisymmetric=False, split_n0_turb=False)

Bases: object

get_array(var_name)

Returns the raw NumPy array for all steps of a given variable.

Parameters:

var_name (str) – Name of the variable to extract (e.g., “eden”, “time”, etc.).

Returns:

array –

A NumPy array of shape:

• (n_steps,) for scalar variables (e.g., time)

• (n_steps, n_planes, n_vertices) for mesh-based variables

Return type:

np.ndarray

xgc_analysis.fieldline_mapping module

xgc_analysis.fieldline_mapping.compute_fieldline_mapping(mesh, magnetic_field, *, rk4_substeps=2, n_int=250, direction='forward')

Build a magnetic‑field‑line mapping between adjacent planes.

Parameters:

• mesh (Mesh) – Instance of the Mesh defined in mesh.py.

• magnetic_field (MagneticField) – Instance of MagneticField (see magnetic_field.py).

• rk4_substeps (int, optional) – Number of RK4 sub‑steps per plane‑to‑plane advance (default 2).

• n_int (int, optional) – Interpolation‑grid resolution used when creating a matplotlib.tri.TriFinder (default 250).

• direction ({"forward", "backward"}, optional) – Direction of tracing; forward increases φ, backward decreases φ.

Returns:

triangle_index(n_vert, n_steps) int32

Triangle hosting each intersection.

bary_weights(n_vert, n_steps, 3) float64

Barycentric weights inside that triangle.

plane_index(n_steps,) int32

Index of the target plane within mesh.planes.

delta_phifloat

Toroidal separation between neighbouring planes (positive).

directionstr

Copy of the direction argument.

Return type:

dict

xgc_analysis.magnetic_field module

class xgc_analysis.magnetic_field.MagneticField(plane_instance, data_dir='.')

Bases: object

compute_background_field(R, Z)

Return the magnetic-field components (B_R, B_Z, B_t) at the given cylindrical coordinates.

Parameters:

• R (float or ndarray) – Cylindrical coordinates. Scalars and arbitrary-shaped NumPy arrays are both accepted; broadcasting is supported as long as R and Z have compatible shapes.

• Z (float or ndarray) – Cylindrical coordinates. Scalars and arbitrary-shaped NumPy arrays are both accepted; broadcasting is supported as long as R and Z have compatible shapes.

Returns:

B_R, B_Z, B_t – Radial, vertical (poloidal) and toroidal components. Their type/shape matches the broadcasted shape of R and Z.

Return type:

float or ndarray

xgc_analysis.mesh module

class xgc_analysis.mesh.Mesh(is_axisymmetric, data_dir='.')

Bases: object

A general mesh class that consists of an array of “plane” objects.

The mesh data for all planes comes from the same BP file (“xgc.mesh.bp”). The file also contains the variables “delta_phi” and “wedge_angle”. The number of wedges is computed as:

wedge_n = 2*pi / wedge_angle,

and the number of planes is computed as:

nphi = wedge_angle / delta_phi.

If the input flag is_axisymmetric is True, then the mesh is axisymmetric and only one plane is set up. Otherwise, nphi planes are created (all identical in this version).

apply_fieldline_mapping(meshdata, *, ref_plane: int = 0)

Pull the data of every plane onto the vertices of ref_plane along magnetic-field lines.

Returns:

If the stored field is scalar:

shape (n_steps, n_vert)

If the stored field is vector:

shape (n_steps, n_vert, n_components)

Axis-0 walks along the field line (φ = φ₀ + (s+1)Δφ); Axis-1 enumerates the vertices of the reference plane.

Return type:

np.ndarray

evaluate_grad_par(field_2d: ndarray, bt_sign: int) → ndarray

Compute b·∇(field) on an axisymmetric mesh.

Parameters:

• field_2d (np.ndarray, shape (nphi, n_n)) – Scalar field values on the full mesh (first index = poloidal plane, second = node).

• bt_sign (int, sign of the toroidal magnetic field ( +1 or -1 ).)

Returns:

Same shape as field_2d with the parallel derivative.

Return type:

np.ndarray

Notes

Uses the centred, 2-nd-order accurate stencil

b·∇f ≈ sgn * [ -l_r/(l_l l_tot) f_{i-1}

• (l_r-l_l)/(l_l l_r) f_i

• l_l/(l_r l_tot) f_{i+1} ]

where the forward / reverse projection matrices ff_1dp_fwd and ff_1dp_rev as well as the parallel arc-lengths dl_par_1dp_fwd and dl_par_1dp_rev are taken from planes[0] (valid only when self.is_axisymmetric is True).

evaluate_scalar_gradients(field_2d: ndarray) → ndarray

Apply per-plane gradient operators to a scalar field on the full 3-D mesh.

Parameters:

field_2d (np.ndarray, shape (nphi, n_n)) – field_2d[i_phi, i_node] holds the scalar value at node i_node on poloidal plane i_phi.

Returns:

grads – grads[i_phi] is the (n_n, 2) array returned by the corresponding plane’s evaluate_scalar_gradients.

Return type:

np.ndarray, shape (nphi, n_n, 2)

ff_map_real2ff(meshdata: ndarray) → ndarray

Maps mesh data from the regular representation of data on planes (points with the same index on different planes are connected toroidally) to field-following representation (vertices with the same index are connected along the magnetic field). The output format is list of MeshData objects of the same type, one for the left, one for the right plane.

flux_avg_from_surf(flux_avg: ndarray) → ndarray

Map flux‑surface averaged data back to vertices on every plane (producing a representative PlaneData).

flux_avg_to_from_surf(meshdata: ndarray) → ndarray

Compute flux-surface average on each plane and project back to the vertices φ.

flux_avg_to_surf(meshdata: ndarray) → ndarray

Flux‑surface average (⟨f⟩_ψ) of a 3‑D field, averaged over φ.

Returns:

Shape (nsurf,) for scalar, or (nsurf, n_components) for vector data.

Return type:

ndarray

flux_avg_to_surf_norm(meshdata: ndarray) → ndarray

Flux‑surface average of vector‑norm, averaged over the toroidal direction.

get_plane(index=0)

Returns the plane at the specified index. If the mesh is axisymmetric, index is ignored.

get_vertex_counts()

Returns a list of the number of vertices on each of the mesh’s planes.

For an axisymmetric mesh (is_axisymmetric=True), the same number of vertices (from the single plane stored) is repeated for each plane (nphi times).

For a non-axisymmetric mesh, the function returns the number of vertices for each plane.

setup_fieldline_mapping(magnetic_field, *, direction='forward', rk4_substeps: int = 2, n_int: int = 250)

Trace all mesh vertices along B and build one CSR propagator matrix for every stored toroidal step. The objects are cached as

self.fl_mapping – raw dict from compute_fieldline_mapping self.fl_csr_matrices – list[CSR] (length = n_steps)

Call this after construction if the mapping was not required earlier.

toroidal_average(meshdata: ndarray) → ndarray

Simple arithmetic mean over the toroidal direction (returns vertex‑shaped ndarray).

xgc_analysis.mesh_data module

class xgc_analysis.mesh_data.MeshData(mesh, data_array=None, field_type='scalar', n_components=1, dtype=<class 'numpy.float64'>, mesh_is_axisym=False)

Bases: object

Stores field data for a mesh as a list of PlaneData instances (one per plane).

The constructor accepts:

• vertex_counts: a list of integers indicating the number of vertices on each plane.

• field_type: ‘scalar’ or ‘vector’.

• n_components: for vector fields (2 or 3); for scalar fields this should be 1.

• dtype: the NumPy data type for the stored data.

• mesh_is_axisym (bool): True if the mesh is axisymmetric (i.e. all planes have the same number of vertices

  and the same (R,Z) coordinates), False otherwise.

The data for each plane is stored as a separate PlaneData instance.

apply_fieldline_mapping(*, ref_plane: int = 0)

Pull the data of every plane onto the vertices of ref_plane along magnetic-field lines.

Returns:

If the stored field is scalar:

shape (n_steps, n_vert)

If the stored field is vector:

shape (n_steps, n_vert, n_components)

Axis-0 walks along the field line (φ = φ₀ + (s+1)Δφ); Axis-1 enumerates the vertices of the reference plane.

Return type:

np.ndarray

evaluate_full_gradients(bt_sign: int) → MeshData

Return a 3-component gradient vector field as a new MeshData:

grads[…, 0] = ∂f/∂ψ grads[…, 1] = ∂f/∂θ grads[…, 2] = b·∇f

Parameters:

bt_sign (int) – Sign of the toroidal magnetic field Bᵗ (±1).

Return type:

MeshData # field_type = ‘vector’

evaluate_grad_par(bt_sign: int) → MeshData

Parallel derivative b·∇(f) returned as a new scalar MeshData.

Parameters:

bt_sign (int) – Sign of the toroidal magnetic field (±1).

Return type:

MeshData # field_type = ‘scalar’

evaluate_scalar_gradients() → MeshData

Return (R,Z)/(psi,theta)-gradients as a new MeshData object (vector field).

Return type:

MeshData # field_type = ‘vector’

extract_component(component_index)

For vector fields only.

Returns a new MeshData instance where each plane contains only the specified component of the original vector field (i.e. a scalar field).

ff_map_real2ff()

Maps mesh data from the regular representation of data on planes (points with the same index on different planes are connected toroidally) to field-following representation (vertices with the same index are connected along the magnetic field). The output format is list of MeshData objects of the same type, one for the left, one for the right plane.

Returns:

MeshData (vector-type): Field on a mesh mapped to field-aligned representation.

flux_avg_from_surf()

Project the flux-surface average data back to vertices and package into PlaneData.

flux_avg_to_from_surf()

Compute the flux-surface average of self.get_data() and project back to the plane in self.flux_avg_plane.

flux_avg_to_surf()

Compute the flux-surface average and store in self.flux_avg_psi.

flux_avg_to_surf_norm()

Flux‑surface average of vector‑norm, averaged over the toroidal direction.

get_data()

Returns the data from all planes as a 2D NumPy array with shape (nphi, n_vertices).

get_plane(plane_index)

Returns the PlaneData instance for the specified plane.

toroidal_average()

Return toroidal average as PlaneData (also cached on Mesh).

xgc_analysis.periodic_table module

xgc_analysis.periodic_table.get_element_name_by_mass(mass_au, tolerance=0.5)

Returns the name of the element with atomic mass closest to the given mass_au.

Parameters:

• mass_aufloat

  Mass in atomic mass units.

• tolerancefloat

  Maximum allowed deviation to consider a match.

Returns:

• str: Element symbol or ‘Unknown’ if no match is found.

xgc_analysis.plane module

class xgc_analysis.plane.Plane(filename='xgc.mesh.bp', data_dir='.', setup_ff2real_mapping=False)

Bases: object

A class to represent a planar unstructured triangle mesh (in cylindrical coordinates at constant toroidal angle) read from an Adios BP (BP version 5) file using the Adios2 Python FileReader API (v2.10.2+).

The BP file (e.g., ‘xgc.mesh.bp’) is assumed to contain the following variables:

Scalars:

• n_n : Number of mesh vertices.

• n_t : Number of triangles.

• nsurf : Number of psi contour surfaces.

• surf_maxlen : Maximum number of vertices on any psi surface.

Arrays:

• rz : (n_n, 2) coordinates (R, Z) of vertices.

• nd_connect_list : (n_t, 3) triangle connectivity (1-based indices).

• tr_area : (n_t,) triangle areas.

• node_vol : (n_n,) Voronoi volume (measure 1).

• node_vol_nearest: (n_n,) Voronoi volume (nearest neighbor measure).

• psi : (n_n,) poloidal magnetic flux.

• theta : (n_n,) straight field-line poloidal angle.

• region : (n_n,) region index of the vertices (1) core, (2) SOL, (3) privat flux, (100) wall

• surf_len : (nsurf,) number of vertices on each psi surface.

• surf_idx : (nsurf, surf_maxlen) 1-based vertex indices on each psi surface.

• rmin : (nsurf,) minor radius of each surface.

• rmaj : (nsurf,) major radius of each surface.

• epsilon : (nsurf,) inverse aspect ratio.

• m_max_surf : (nsurf,) Fourier space resolution limit.

• psi_surf : (nsurf,) psi value on each surface.

• qsafety : (nsurf,) safety factor.

• trapped : (nsurf,) trapped particle fraction.

Plotting functionality:

• triObjMatplotlib triangulation object for unstructured meshes.

  This is needed for plotting fields defined on the planar mesh

Matrices:

• cnv_to_surf : …

• cnv_from_surf : …

• grad_r_psi : …

• grad_z_theta : …

Diagnostic functions:

• plot_vertices_with_wall : …

• plot surfaces : …

• plot_triangles : …

calculate_flux_surface_area()

Calculates the surface area of each flux surface assuming that the 2D ψ-surfaces are rotated around the Z-axis by 2π.

For each ψ-surface, the routine:

• Loops over its vertices.

• For each segment between adjacent vertices, computes the midpoint (R_m, Z_m) and the distance l_seg.

• If the surface is closed (i.e. self.surf_region[i] == 1), the last vertex connects to the first; if open (region==2 or 3), the last vertex is skipped.

• The area of each segment is 2π * R_m * l_seg.

Returns:

A list of surface areas (one per ψ-surface).

Assumes that the Plane class has the following attributes:

• self.nsurf: total number of ψ-surfaces.

• self.surf_len: 1D array of length nsurf (number of vertices per surface).

• self.surf_idx: 2D array (nsurf × surf_maxlen) of 0-based vertex indices.

• self.rz: 2D array (n_n × 2) with (R,Z) coordinates.

• self.surf_region: 1D array (length nsurf) with region codes (1 for closed, 2 or 3 for open).

evaluate_scalar_gradients(field_1d: ndarray) → ndarray

Evaluate the radial and poloidal gradients of a scalar field. (Note that it depends on the value of self.basis whether the derivatives are (R,Z) or (psi,theta) directonal derivatives!)

Parameters:

field_1d (np.ndarray, shape (n_n,)) – Nodal values of the scalar field on this plane.

Returns:

grads – Column 0 → ∂f/∂ψ (gradient_r_psi ⋅ f) Column 1 → ∂f/∂θ (gradient_z_theta ⋅ f)

Return type:

np.ndarray, shape (n_n, 2)

flux_avg_from_surf(flux_avg: ndarray) → ndarray

Map a flux-surface averaged quantity <f>(psi) back to vertices.

flux_avg_to_from_surf(vertex_data: ndarray) → ndarray

Smooth vertex_data by averaging on ψ and mapping back to vertices.

flux_avg_to_surf(vertex_data: ndarray) → ndarray

Return the flux-surface average ⟨f⟩ for a vertex‑defined scalar or vector field f.

flux_avg_to_surf_norm(field: ndarray) → ndarray

Flux-surface average of the Euclidean norm of a vertex-centred field on this poloidal plane.

Parameters:

field (np.ndarray) –

• If shape == (n_n,) → treated as a scalar field.

• If shape == (n_n, n_comp) → treated as a vector field; the norm is computed along the last axis.

Returns:

1-D array of length n_psi (the number of flux-surfaces on the plane) containing the volume-weighted flux-surface averages.

Return type:

np.ndarray

Raises:

• RuntimeError – If the conversion matrix cnv_to_surf is unavailable because its backing .bp file was missing when the Plane was constructed.

• ValueError – If the input array does not have the required first dimension n_n.

gen_ff2real_mapping()

get_psi_surface_index_map(tol=0.0001)

Identify the ψ‑surfaces that are valid for subsequent analysis.

Strategy

For every surface i (whose length is surf_len[i]):

1. Single‑vertex surface

   • If the surface contains exactly one vertex and that vertex’s coordinates (R, Z) are within tol of the magnetic‑axis position (axis_r, axis_z), accept the surface.

2. Multi‑vertex surface

   • Restrict to vertices with R >= axis_r.

   • Among them, find • the vertex whose Z is largest below the axis • the vertex whose Z is smallest above the axis.

   • Accept the surface only if both vertices exist.

param tol:

Tolerance that defines “close to the magnetic axis”, measured in the same units as R and Z. Default is 1e‑4.

type tol:

float, optional

returns:

Zero‑based indices of the ψ‑surfaces that satisfy the above criteria.

rtype:

list[int]

get_surface_indices()

Returns:

np.ndarray: 2D array (shape: (nsurf, surf_maxlen)) containing 1-based vertex indices on each surface.

get_surface_lengths()

Returns:

np.ndarray: 1D array of length nsurf indicating the number of vertices on each surface.

get_surface_properties()

Returns:

dict: A dictionary of surface properties with keys:

• ‘rmin’ : Minor radius (np.ndarray of shape (nsurf,))

• ‘rmaj’ : Major radius (np.ndarray of shape (nsurf,))

• ‘epsilon’ : Inverse aspect ratio (np.ndarray of shape (nsurf,))

• ‘m_max_surf’ : Fourier space resolution limit (np.ndarray of shape (nsurf,))

• ‘psi_surf’ : Psi value on each surface (np.ndarray of shape (nsurf,))

• ‘qsafety’ : Safety factor (np.ndarray of shape (nsurf,))

• ‘trapped’ : Trapped particle fraction (np.ndarray of shape (nsurf,))

get_surface_vertex_indices(surface_index)

Returns the vertex indices for a given psi surface.

Args:

surface_index (int): The index of the surface (0-based).

Returns:

np.ndarray: An array of vertex indices (0-based) corresponding to the given surface.

get_triangle_areas()

Returns:

np.ndarray: Array of triangle areas with shape (n_t,).

get_triangle_connectivity()

Returns:

np.ndarray: Array of triangle connectivity indices (1-based) with shape (n_t, 3).

get_triangulation_data(n_int: int = 200, R_bounds: tuple[float, float] | None = None, Z_bounds: tuple[float, float] | None = None)

Build a regular (R,Z) interpolation grid and a Matplotlib Triangulation that matches the unstructured mesh.

Parameters:

• n_int (int, default 200) – Number of grid points in each dimension (R and Z).

• R_bounds ((min, max) tuple or None, optional) – Explicit coordinate limits for the interpolation grid. If None (default) the routine uses the min/max of the vertex coordinates.

• Z_bounds ((min, max) tuple or None, optional) – Explicit coordinate limits for the interpolation grid. If None (default) the routine uses the min/max of the vertex coordinates.

Returns:

Object containing (RI, ZI, triObj).

Return type:

TriangulationData

get_vertex_coordinates()

Returns:

np.ndarray: Array of vertex coordinates with shape (n_n, 2) representing (R,Z) pairs.

get_voronoi_volumes()

Returns:

tuple: Two np.ndarray objects representing the vertex Voronoi volumes

(node_vol and node_vol_nearest), each with shape (n_n,).

plot_contour(values_in: ndarray, *, n_int: int = 200, plot_norm: bool = False, i_comp: int = 0, var_name: str = '', title: str = '', levels: int | None = 60, lower: float | None = None, upper: float | None = None, R_bounds: tuple[float, float] | None = None, Z_bounds: tuple[float, float] | None = None, filename: str | None = None)

Plot the field stored in this PlaneData instance as a filled-contour map in the (R,Z) plane.

Parameters:

• values_in (np.ndarray) – 1-D array of length self.n_n (= number of vertices). If you pass a 2-D array, its Euclidean norm is taken or a single component (convenient for vector data), controlled by plot_norm and i_comp.

• n_int (int, default 200) – Grid resolution in both R and Z used for interpolation.

• plot_norm (bool, default False) – In case of vector data, whether to plot the norm | v | of the vector.

• i_comp (int, default 0) – In case of vector data, which component to plot.

• var_name (str, optional) – Label for the colour-bar (e.g. "nₑ (m⁻³)"). Defaults to "".

• title (str, optional) – Header for the plot window. Defaults to "".

• levels (int or sequence, optional) – If an int, the number of evenly-spaced contour levels (Matplotlib default is fine). If an explicit 1-D sequence, those values are used as contour levels. Defaults to 60.

• lower (float or None, optional) – Explicit lower/upper bounds for the colour-map. If either is None the routine uses the corresponding data extreme.

• upper (float or None, optional) – Explicit lower/upper bounds for the colour-map. If either is None the routine uses the corresponding data extreme.

• R_bounds ((min, max) tuple or None, optional) – Explicit limits for the interpolation grid. Passed straight through to plane.get_triangulation_data. None ⇒ auto-limits.

• Z_bounds ((min, max) tuple or None, optional) – Explicit limits for the interpolation grid. Passed straight through to plane.get_triangulation_data. None ⇒ auto-limits.

plot_fieldline_contour(data: ndarray, surf_idx: int, *, var_name: str = '', title: str = '', levels: int | ndarray | None = 60, lower: float | None = None, upper: float | None = None, filename: str | None = None)

Filled-contour plot of data along the magnetic field on the requested ψ-surface.

Parameters:

• data (ndarray) – Result of MeshData.apply_fieldline_mapping, shape (n_steps, n_vertices). Axis-0 = distance along B, axis-1 = vertex index on this plane.

• surf_idx (int) – Which ψ-surface to plot (index into self.surf_* arrays).

• var_name (str) – Labels for colour-bar and figure title.

• title (str) – Labels for colour-bar and figure title.

• levels (int or 1-D array or None) – Passed to contourf, like in plot_contour.

• lower (float or None) – Manual colour limits (None → data min/max).

• upper (float or None) – Manual colour limits (None → data min/max).

plot_surfaces()

Plot all psi surfaces.

plot_triangles()

Plots all the triangles of the mesh. For each triangle, the 3 vertex indices from nd_connect_list are used to get the corresponding (R, Z) coordinates from rz. The first vertex is appended to close the triangle, and the outline is plotted.

plot_vertices_with_wall()

Plots the mesh vertices as dots and overlays the wall as a closed line in the background. The wall is defined by the grid_wall_nodes variable.

class xgc_analysis.plane.TriangulationData(RI, ZI, triObj)

Bases: object

A simple container class to hold interpolation grid data and a Matplotlib Triangulation object.

Attributes:

RI (np.ndarray): 2D array of R coordinates of the interpolation grid. ZI (np.ndarray): 2D array of Z coordinates of the interpolation grid. triObj (Triangulation): A Matplotlib Triangulation object for the mesh.

xgc_analysis.plane_data module

class xgc_analysis.plane_data.PlaneData(plane, data_array=None, n_components=1, dtype=<class 'numpy.float64'>)

Bases: object

Stores field data for a single plane.

For a scalar field, data is a 1D NumPy array of shape (n,). For a vector field, data is a 2D NumPy array of shape (n, n_components), where for each vertex the vector components are stored contiguously in memory.

Attributes:

• n: number of vertices in this plane.

• field_type: ‘scalar’ or ‘vector’.

• n_components: number of components (1 for scalar, 2 or 3 for vector fields).

• dtype: data type for the array.

• data: the NumPy array storing the field data.

evaluate_scalar_gradients() → PlaneData

Compute ∇_{R,Z}f and return the result as a new PlaneData instance (two-component vector on the same plane).

Return type:

PlaneData # field_type = ‘vector’

extract_component(component_index)

For vector fields only.

Returns a new PlaneData instance (a scalar field) containing the specified component from the vector.

flux_avg_from_surf()

Back‑project flux_avg_psi and store result in self.flux_avg_plane.

flux_avg_to_from_surf()

Compute the flux-surface average of self.data and project back to the plane in self.flux_avg_plane.

flux_avg_to_surf()

Compute flux-surface-average and store it in self.flux_avg_psi.

flux_avg_to_surf_norm()

Compute flux-surface-average of the norm of a vector field and store it in self.flux_avg_norm_psi.

get_data()

Returns the underlying NumPy array.

plot_contour(plane: Plane, **kwargs)

Convenience wrapper that forwards to Plane.plot_contour using this instance’s data.

Any keyword arguments accepted by Plane.plot_contour can be supplied.

xgc_analysis.read_bp_file module

xgc_analysis.read_bp_file.ReadBPFile(filename, variables=None, step_range=None)

Reads data from an ADIOS2 BP file using the Stream API (ADIOS2 v2.10+).

Parameters:

• filename (str) – Path to the ADIOS2 BP file.

• variables (list of str, optional) – List of variable names to read. If None, all variables will be read.

• step_range (tuple of int, optional) – Tuple (start, end) indicating the range of steps to read. If None, only the last step is read.

Returns:

data – Dictionary with structure: {step_idx: {variable_name: np.ndarray}}.

Return type:

dict

xgc_analysis.simulation module

class xgc_analysis.simulation.Simulation(directories=['./'], is_stellarator=False, sim_is_axisymmetric=False)

Bases: object

convert_fortran_value(val)

xgc_analysis.species module

class xgc_analysis.species.Species(sim, index)

Bases: object

Class representing a plasma species with basic physical properties.

interpolate_profile(name: str, psi_new: ndarray, *, fill: str | float | None = 'edge') → ndarray

Interpolate one of this species’ initial profiles onto a new ψ grid.

Parameters:

• name (str) – Profile key in self.initial_profiles (“density”, “temperature”, “flow”, …).

• psi_new (ndarray) – Target ψ‑grid (1‑D NumPy array).

• fill ({"edge", None, float}, default "edge") –

  How to handle ψ outside the original range:

  • ”edge” – use nearest edge value (default)

  • None – raise ValueError if extrapolation needed

  • float – fill with that scalar

Returns:

Interpolated profile values at psi_new.

Return type:

ndarray

Module contents