Field-line Tracing (pyna.flt)

The pyna.flt package provides field-line integration routines built on top of the abstract pyna.system hierarchy.

Backends:

• CPU/serial – pure-Python RK4 integrator

• CPU/parallel – multi-process or multi-threaded variants

• CUDA – optional GPU backend via CuPy (up to 118× speedup)

• OpenCL – experimental

---

Core Tracer

Field line tracing with parallel execution.

Parallelism strategy (Windows + Python 3.13 standard GIL build): - Primary: ThreadPoolExecutor — works with any callable, numpy releases GIL

during integration loops, so real speedup is achievable

• ProcessPoolExecutor: available but requires module-level picklable functions; use only when explicitly requested and field_func is picklable

• CUDA: CuPy backend for GPU-accelerated batch tracing

Progress reporting

All batch methods accept an optional progress= parameter. Pass any TraceProgressBase instance to receive per-task or aggregate progress updates:

from pyna.progress import TqdmProgress, LogFileProgress, CompositeProgress

# tqdm bar in CLI / notebook
tracer.trace_many(starts, t_max, progress=TqdmProgress())

# bar + log file simultaneously
prog = CompositeProgress([TqdmProgress(), LogFileProgress("run.jsonl")])
tracer.trace_many(starts, t_max, progress=prog)

Wall-hit handling (boundary field lines)

Near the plasma boundary, field lines may strike the first wall after only a few integration steps. FieldLineTracer.trace_many marks these as “wall-hit” (non-blocking) and returns the short trajectory. Use reseed_boundary_field_lines() to automatically add extra seed points around the boundary where the hit-fraction is high, ensuring the boundary topology is still captured.

Legacy API

已移除。原 bundle_tracing_with_t_as_DeltaPhi、save_Poincare_orbits、 load_Poincare_orbits 及 OdeSolution.mat_interp monkey-patch 已删除。 请使用 FieldLineTracer.trace_many 替代。

New API

FieldLineTracer

RK4 integrator with .trace() / .trace_many() (ThreadPoolExecutor). 当传入网格数据时底层由 cyna C++ 扩展驱动；若只有 callable 场函数， 则使用纯 Python RK4 fallback（向后兼容）。 若要完全使用 C++ 底层加速，请使用 pyna.toroidal.flt.trace_poincare_batch 并传入网格数组。

WallModel

Parametric / polygon wall geometry for wall-hit detection.

reseed_boundary_field_lines

Adaptively add more seed points where wall-hits are dense.

get_backend(mode, **kwargs)

Factory: ‘cpu’, ‘cuda’, ‘opencl’.

class pyna.flt.FieldLineTracer(field_func, dt=0.04, RZlimit=None, n_workers=None)

Bases: object

RK4 field-line integrator with parallel trace_many.

The field function signature:

f(rzphi: ndarray[3]) -> ndarray[3]   # (dR/dl, dZ/dl, dphi/dl)

where rzphi = [R, Z, phi] is the current position.

Parallelism uses ThreadPoolExecutor on all platforms. This is deliberately chosen over ProcessPoolExecutor because:

• Works with any callable — no pickle requirement (important on Windows).

• NumPy/SciPy release the GIL during heavy computation, so threads achieve real parallel speedup.

• No process-spawn overhead on Windows.

Parameters:

• field_func (callable) – Unit tangent vector field.

• dt (float) – Integration step size (arc length).

• RZlimit (tuple or None) – Optional (R_min, R_max, Z_min, Z_max) domain boundary. Integration stops when the trajectory leaves this box.

• n_workers (int or None) – Default thread-pool size for trace_many(). None → min(os.cpu_count(), 16).

trace(start_pt, t_max, wall=None)

Trace a single field line, optionally stopping at a wall.

Parameters:

• start_pt (array-like of length 3) – Starting point (R, Z, phi).

• t_max (float) – Maximum arc-length parameter.

• wall (WallModel or None) – Wall model for early-termination on wall hit. None uses the RZlimit box-boundary from the constructor.

Return type:

ndarray

Returns:

ndarray, shape (N, 3) – Trajectory points (R, Z, phi). For wall-hit cases N may be much smaller than int(t_max / dt); the trajectory is not padded.

trace_many(start_pts, t_max, n_workers=None, progress=None, wall=None, min_valid_steps=0)

Parallel field-line tracing with optional wall-hit handling.

Wall-hitting trajectories are returned as short arrays and do not block or slow down other traces in the same batch.

Parameters:

• start_pts (array-like, shape (N, 3)) – Starting points.

• t_max (float) – Maximum arc-length for each field line.

• n_workers (int or None) – Override the default worker count.

• progress (TraceProgressBase or None) – Optional progress reporter.

• wall (WallModel or None) – Wall model. When provided, trajectories that hit the wall are returned early (shorter arrays).

• min_valid_steps (int) – If > 0, trajectories with fewer than this many steps are flagged as “wall-hit” in the returned metadata. This does not discard the trajectory — callers can filter using the returned list.

Return type:

List[ndarray]

Returns:

list of ndarray – One trajectory per starting point.

class pyna.flt.WallModel(R_wall, Z_wall, min_clearance=0.0)

Bases: object

First-wall geometry for wall-hit detection during field-line tracing.

A wall is represented as a closed polygon in the (R, Z) poloidal cross- section. The same shape is assumed toroidally symmetric (axisymmetric wall). More complex 3-D first-wall descriptions can be approximated by providing an appropriate 2-D outline.

The model also provides a minimum clearance check so that trajectories that approach very close to the wall but haven’t yet crossed it can be flagged for early termination.

Parameters:

• R_wall (array_like of shape (N,)) – R-coordinates of the wall polygon vertices (m). The polygon is closed automatically; do not repeat the first point.

• Z_wall (array_like of shape (N,)) – Z-coordinates of the wall polygon vertices (m).

• min_clearance (float, optional) – Minimum distance from wall (m) before a point is considered a near-miss / wall-hit. 0 means only count true crossings. Default 0.

Examples

>>> import numpy as np
>>> theta = np.linspace(0, 2*np.pi, 64, endpoint=False)
>>> wall = WallModel(2.0 + 0.5*np.cos(theta), 0.5*np.sin(theta))
>>> wall.is_outside(np.array([2.6, 0.0]))
True

classmethod circular(R0, a, n_vertices=128, min_clearance=0.0)

Construct a circular (axisymmetric) wall.

Parameters:

• R0 (float) – Major radius and minor radius of the circular wall.

• a (float) – Major radius and minor radius of the circular wall.

• n_vertices (int) – Number of polygon vertices.

• min_clearance (float) – Proximity clearance (m).

Return type:

WallModel

Returns:

WallModel

is_outside(RZ)

Return True if the point (R, Z) is outside the wall polygon.

Uses the ray-casting (winding-number) algorithm.

Parameters:

RZ (array_like of shape (2,)) – Point [R, Z] to test (m).

Return type:

bool

Returns:

bool

pyna.flt.get_backend(mode='cpu', **kwargs)

Get a field-line tracer backend.

Parameters:

• mode (str) –

  'cpu' — FieldLineTracer with ThreadPoolExecutor

  (always available).

  'cuda' — FieldLineTracerCUDA

  (requires CuPy, only for analytic fields).

  'opencl' — reserved, not yet implemented.

• **kwargs – Passed to the backend constructor.

Returns:

Backend object with .trace_many(start_pts, t_max) method.

Examples

CPU:

tracer = get_backend('cpu', field_func=my_field, dt=0.02)
trajs  = tracer.trace_many(starts, t_max=100.0)

CUDA:

tracer = get_backend('cuda', R0=1.0, a=0.3, B0=1.0, q0=2.0)
trajs  = tracer.trace_many(starts, t_max=100.0)

pyna.flt.reseed_boundary_field_lines(tracer, start_pts, trajs, t_max, wall=None, min_valid_fraction=0.5, n_reseed_factor=4, reseed_radius=0.02, n_workers=None, progress=None)

Adaptively reseed field lines near wall-hit starting points.

Near the plasma boundary, a fraction of seed points may produce very short trajectories (wall hits). This function:

1. Identifies seeds with trajectories shorter than min_valid_fraction * t_max / dt steps.

2. Surrounds each such seed with n_reseed_factor new seeds in a small circle of radius reseed_radius.

3. Traces the new seeds and returns all trajectories (original + reseeded).

Parameters:

• tracer (FieldLineTracer) – Tracer used for the original batch.

• start_pts (array, shape (N, 3)) – Original starting points [R, Z, phi].

• trajs (list of ndarray) – Trajectories from the original batch (one per start point).

• t_max (float) – Arc-length used in the original batch.

• wall (WallModel or None) – Wall model passed to the new traces.

• min_valid_fraction (float) – Fraction of t_max / dt steps below which a trajectory is considered a wall hit. Default 0.5.

• n_reseed_factor (int) – Number of new seeds to place around each wall-hit seed. Default 4.

• reseed_radius (float) – Radius (m) of the reseeding circle around each hit seed. Default 0.02.

• n_workers (int or None) – Worker count for the reseeded traces.

• progress (TraceProgressBase or None) – Progress reporter for the reseeded traces.

Return type:

List[ndarray]

Returns:

list of ndarray – Original trajectories plus new reseeded trajectories (appended).

---

Dynamical System Hierarchy

Core dynamical system abstractions for pyna.

This module defines the class hierarchy for dynamical systems, from general abstract base classes down to physically motivated special cases.

class pyna.system.AutonomousDynamicalSystem

Bases: DynamicalSystem

Autonomous dynamical system: dx/dt = f(x).

The right-hand side does not depend explicitly on time. Key consequences:

• Poincaré section maps are well-defined without a time reference.

• Invariant manifolds, fixed points, and KAM tori are time-independent geometric objects in state space.

class pyna.system.DynamicalSystem

Bases: ABC

Abstract base for all dynamical systems.

A dynamical system defines the evolution law for a state vector. Subclasses must implement state_dim() and __call__().

abstract property state_dim: int

Dimension of the state space.

class pyna.system.NonAutonomousDynamicalSystem

Bases: DynamicalSystem

Non-autonomous dynamical system: dx/dt = f(x, t).

Time appears explicitly on the right-hand side. Geometric structures such as Poincaré maps depend on the phase of the driving and are therefore more complex to analyse than their autonomous counterparts.

class pyna.system.VectorField

Bases: AutonomousDynamicalSystem

Autonomous dynamical system defined by a smooth vector field.

The vector field is the right-hand side f(x) of dx/dt = f(x). Subclasses specialise to a particular spatial dimension.

abstract property dim: int

Spatial dimension of the vector field.

property state_dim: int

Dimension of the state space.

class pyna.system.VectorField1D

Bases: VectorField

One-dimensional vector field: dx/dt = f(x).

property dim: int

Spatial dimension of the vector field.

class pyna.system.VectorField2D

Bases: VectorField

Two-dimensional vector field: dx/dt = f(x), x ∈ ℝ².

The 2-D case is where KAM theory is most transparent and area-preserving (Hamiltonian) maps are most widely studied.

property dim: int

Spatial dimension of the vector field.

class pyna.system.VectorField4D

Bases: VectorField

Four-dimensional vector field: dx/dt = f(x), x ∈ ℝ⁴.

Relevant for relativistic charged-particle dynamics (guiding-centre equations with energy as a fourth coordinate) and certain Hamiltonian systems with two degrees of freedom.

property dim: int

Spatial dimension of the vector field.