Example scripts

Complete pysimple scripts from examples/, run by test/python/test_examples.py. Run them from examples/; they read test/test_data/wout.nc, the equilibrium of Quick start.

Surface trace

examples/simple_api.py

#!/usr/bin/env python3
"""Simple particle tracing example using pysimple API."""

import pysimple

vmec_file = '../test/test_data/wout.nc'

pysimple.init(vmec_file, deterministic=True, trace_time=5e-5, ntestpart=32)

particles = pysimple.sample_surface(32, s=0.3)

results = pysimple.trace_parallel(particles, integrator='midpoint')

n_lost = ((results['loss_times'] > 0.0) & (results['loss_times'] < 5e-5)).sum()
n_confined = (results['loss_times'] >= 5e-5).sum()
n_skipped = (results['loss_times'] < 0.0).sum()
print(f"Traced 32 particles: {n_confined} confined, {n_lost} lost, {n_skipped} skipped")

Fast classification

examples/classify_fast.py

#!/usr/bin/env python3
"""Fast classification using jpar and topology."""

import pysimple

vmec_file = '../test/test_data/wout.nc'

pysimple.init(vmec_file, deterministic=True, trace_time=1e-4, ntestpart=16)

particles = pysimple.sample_surface(16, s=0.4)

results = pysimple.classify_fast(particles, nturns=8)

n_passing = results['passing'].sum()
n_trapped = (~results['passing']).sum()

jpar_regular = (results['jpar'] == 1).sum()
jpar_stochastic = (results['jpar'] == 2).sum()

topology_ideal = (results['topology'] == 1).sum()
topology_non_ideal = (results['topology'] == 2).sum()

print(f"Fast classification of 16 particles: {n_passing} passing, {n_trapped} trapped")
print(f"  J_parallel: {jpar_regular} regular, {jpar_stochastic} stochastic")
print(f"  Topology: {topology_ideal} ideal, {topology_non_ideal} non-ideal")

Fractal classification

examples/classify_fractal.py

#!/usr/bin/env python3
"""Fractal classification using fractal dimension."""

import pysimple

vmec_file = '../test/test_data/wout.nc'

pysimple.init(vmec_file, deterministic=True, trace_time=1e-4, ntestpart=16)

particles = pysimple.sample_surface(16, s=0.4)

results = pysimple.classify_fractal(particles, tcut=0.1)

n_passing = results['passing'].sum()
n_trapped = (~results['passing']).sum()
n_lost = results['lost'].sum()

fractal_regular = (results['fractal'] == 1).sum()
fractal_chaotic = (results['fractal'] == 2).sum()

print(f"Fractal classification of {16} particles:")
print(f"  Passing: {n_passing}, Trapped: {n_trapped}, Lost: {n_lost}")
print(f"  Fractal: {fractal_regular} regular, {fractal_chaotic} chaotic")

Macrostep orbit trace

examples/orbit_macro.py

#!/usr/bin/env python3
"""Trace orbit with trajectory output."""

import pysimple

vmec_file = '../test/test_data/wout.nc'

pysimple.init(vmec_file, deterministic=True, trace_time=1e-4, ntestpart=1)

particles = pysimple.sample_surface(1, s=0.35)
particle = particles[:, 0]

result = pysimple.trace_orbit(particle, integrator='midpoint', return_trajectory=True)

n_timesteps = result['times'].shape[0]
print(f"Traced 1 particle with {n_timesteps} timesteps")
print(f"Loss time: {result['loss_time']:.3e}s")