Architecture

This page maps the code: how a run flows through the modules, where fields, coordinates and integrators live, and what to change to add a field source or an integrator. Implementation detail of coordinates and fields is in Coordinate Systems and Magnetic Fields in SIMPLE; the equations are in Model and methods.

Run sequence

simple.x calls simple_main::main (src/simple_main.f90):

1. load_config reads &config into module params and &profiles into simple_profiles.

2. init_run loads the field (init_field), derives the normalisation (params_init: \(v_0\), \(\rho_0\), step \(\Delta\tau\)) and initialises collisions.

3. init_sampling_field and sample_start fill zstart according to startmode (src/samplers.f90).

4. init_tracing_field selects the integration field.

5. trace_parallel runs trace_orbit for each marker in an OpenMP loop (schedule(dynamic, 1)); with gpu_mode = 'production', trace_gpu_production replaces it.

6. write_output writes the text files and results.nc.

simple_api (src/simple_api.f90) exposes the same steps as separate calls (simple_init, simple_sample, simple_trace, …) plus single-orbit diagnostics. Public states are \((s, \vartheta, \varphi, v/v_0, \lambda)\) in reference coordinates; trace_orbit converts them with ref_to_integ before integration.

Module map

Module (file)	Role
params (src/params.f90)	&config keys, run state, exit codes, per-marker random streams
simple_main, simple_api	run sequence, per-orbit loop, wall and LCFS events, recovery, output
simple (src/simple.f90)	tracer_t, normalisation, encode_symplectic_state
samplers, fusion_source	start positions; fusion-weighted births
field (src/field.F90)	field_from_file: field object from field_input
src/field/field_*.f90	field classes (coils, Boozer chartmap, splined) and canonical fields (field_can_flux, _boozer, _meiss, _albert, _spectre, _test)
field_can_mod (src/field_can.f90)	canonical field dispatch, integ_to_ref / ref_to_integ
magfie_sub (src/magfie.f90)	non-canonical field evaluation for RK45 and sampling
reference_coordinates, src/coordinates/	reference coordinate system ref_coords, radial scaling
orbit_symplectic_base, orbit_symplectic, orbit_symplectic_quasi, orbit_symplectic_euler1	symplectic and quasi-symplectic steppers, LCFS event location
alpha_lifetime_sub (src/sub_alpha_lifetime_can.f90)	RK45 guiding-centre equations (velo_can) and recovery steps
orbit_fo_boris, orbit_fo_field	full-orbit Boris pusher
spectre_orbit, spectre_sympl_orbit, spectre_fo_field, spectre_fo_hybrid, field_line_spectre	SPECTRE tracing per volume, full orbit near interfaces, field lines
collis_alp (src/collis_alphas.f90)	Coulomb collision operator (stost)
classification, detect_oneline_mod, cut_detector	orbit classification
stl_wall_intersection, src/wall/stl_wall_cgal.cpp	STL wall segment test through CGAL
orbit_frequencies, orbit_invariants	bounce and toroidal frequencies; invariants to starting states
netcdf_results_output, netcdf_orbit_output	results.nc, orbits.nc
simple_gpu, simple_cuda_native, *.cu	OpenACC and native CUDA tracing for Boozer fields
python/pysimple/	Python layer over simple_api and params

Libraries: libneo supplies VMEC splines, Boozer conversion, coordinate systems (libneo_coordinates), the field base class magnetic_field_t, Biot–Savart, SPECTRE fields and the ODE and spline routines; fortplot draws diagnostic plots; fortnum provides the root finder and the GPU RK5(4) integrator. Versions are pinned in cmake/dependency_versions.json.

Fields and coordinates

Three layers separate field data from the coordinates used for integration.

1. Field objects extend libneo magnetic_field_t and implement evaluate(x, Acov, hcov, Bmod) in their own coordinates: vmec_field_t (libneo), coils_field_t, boozer_chartmap_field_t, splined_field_t (any field resampled on reference coordinates) and spectre_field_t (libneo). field_from_file chooses the class from the file.

2. Reference coordinates (ref_coords, a libneo coordinate_system_t) are VMEC \((s,\vartheta,\varphi)\) or a chartmap \((\rho,\vartheta,\zeta)\), set by coord_input.

3. Integration coordinates are chosen by integ_coords (Field sources and coordinates). Canonical choices fill a field_can_t with \(A_\vartheta, A_\varphi, h_\vartheta, h_\varphi, B\) and their first and second derivatives, then add \(H, p_\vartheta, v_\parallel\) and their derivatives in \((r,\vartheta,\varphi,p_\varphi)\). field_can_from_name sets the procedure pointers evaluate, integ_to_ref and ref_to_integ. RK45 uses magfie, a procedure pointer set by init_magfie(id), which returns \(B\), \(\sqrt g\), \(\nabla\ln B\), \(h_i\), \(h^i\) and \((\nabla\times\mathbf h)^i\).

Integrators

tracer_t (src/simple.f90) holds one marker’s field_can_t, symplectic_integrator_t, multistage_integrator_t (Gauss and Lobatto stages) and full-orbit state. orbit_sympl_init maps integmode (Running SIMPLE) to a step routine raw_timestep_sympl and its quasi-symplectic fallback; orbit_timestep_sympl wraps it with the LCFS event search. RK45 (integmode = 0) integrates velo_can with libneo odeint_allroutines. The GPU paths do not use procedure pointers: they evaluate the Boozer field locally and reuse the step algebra of orbit_symplectic_euler1 and orbit_rk_core.

Python binding

python/CMakeLists.txt wraps simple_api.f90, params.f90, orbit_symplectic_base.f90, orbit_frequencies.f90 and orbit_invariants.f90 with f90wrap into the extension _simple_backend. pysimple adds array handling and plotting only; physics stays in Fortran so that simple.x and pysimple agree. A Fortran program can link the CMake target simple and use simple_api.

Adding a field source

1. Write type, extends(magnetic_field_t) in src/field/ with evaluate returning covariant \(A_i\), \(h_i\) and \(B\) in Gaussian units (G, cm) in its native coordinates; convert units on load as create_coils_field does.

2. Detect the file in field_from_file (src/field.F90). If the native coordinates are not the reference coordinates, wrap the field with create_splined_field as for coils; it then runs with Meiss, Albert or RK45 on REFCOORDS.

3. A field that needs its own canonical coordinates also needs a field_can_<name>.f90 with evaluate, integ_to_ref and ref_to_integ, a case in field_can_from_name and an id in magfie_sub.

4. Add a roundtrip test and an orbit comparison against an existing source on the same equilibrium, and update Coordinate Systems and Magnetic Fields in SIMPLE.

Adding an integrator

1. Add a constant in orbit_symplectic_base and a step routine in orbit_symplectic.f90 that advances si%z from a field_can_t and returns a SYMPLECTIC_STEP_* status; reuse the Newton and LCFS event helpers.

2. Add the case to orbit_sympl_init and, if it has stages, size multistage_integrator_t.

3. Document the new integmode value in examples/simple_full.in.

4. Test the convergence order on the analytic field (integ_coords = -1), the boundedness of \(H\) and \(p_\varphi\) and a short loss comparison against an existing scheme.

IMAS actor

python/pysimple/imas (extra imas; usage: IMAS workflow) is Python on top of pysimple; physics stays in the Fortran library.

Module	Role
_common.py	DD 3/4 check, exact-time slice selection, array and scalar validation, output provenance (code.name, version)
markers.py	Markers; read_markers/write_markers by coordinate identifier for distribution_sources and distributions
equilibrium.py	equilibrium IDS to COCOS 3 G-EQDSK, then libneo convert_eqdsk_to_chartmap in a child Python process
plasma.py	core_profiles to a (6, n) table and collision settings
actor.py	run, run_entry and the simple-imas entry point

A wall IDS is rejected by run and run_entry before any tracing: the equilibrium chartmap ends at the LCFS, so a wall would never be queried. Wall tracing for IDS input needs a chartmap beyond rho_lcfs first.

run converts the equilibrium once, then per population calls pysimple.init with the actor-owned keys and plasma_profiles=table, cylindrical_to_reference, trace_parallel, reference_to_cylindrical and fills one distributions.distribution entry. The native side adds to src/simple_api.f90: simple_set_profile_table (set_profile_table in src/profiles.f90, linear in \(s\), replaces &profiles until reset_profiles), simple_cylindrical_to_reference, simple_reference_to_cylindrical, simple_reference_info (coordinate semantics: Coordinate Systems and Magnetic Fields in SIMPLE) and simple_wall_supported (compile-time CGAL flag stl_wall_available).

Keep when changing it:

• libneo’s field_eq shifts the G-EQDSK flux by its axis value and stops the surface scan at psimax in the axis-to-edge direction, so equilibrium.py passes the signed psimax \(=(\psi_\text{edge}-\psi_\text{axis})\cdot10^8\) G cm² (Wb/rad to G cm²). An unsigned value ends reversed-flux scans at the first surface (millimetre-sized chartmap). The signed scan needs libneo with itpplasma/libneo#434, pinned at 151f58b in pyproject.toml.

• The converter runs in a child process: a Fortran STOP or module state in efit_to_boozer must not reach the caller.

• Native exit codes other than 0, 1, 2 and invalid final states raise; they are never written as losses. Weights are never converted between rate and count implicitly.

• Every output IDS passes validate(); the coordinate API keeps its length scale after the temporary chartmap is deleted.

Tests: test/python/test_imas.py (CTest test_imas, labels imas, slow; needs pysimple[imas] with libneo’s _efit_to_boozer) covers NetCDF round trips in DD 3.42.2 and 4.1.1, identifier reordering, COCOS signs against an analytic circular field, chartmap extent for all flux and field signs, the coordinate round trip, particle and energy totals, rejection of invalid native states, the DBEntry path with core_profiles, and both wall forms; test/tests/test_profiles.f90 covers the profile table. A new IDS mapping gets a reader that validates through _common, selects time with select_slice and a test with real IMAS-Python IDS of both DD versions.

Invariants to keep

• Gaussian CGS inside the code; SI and eV only at input and output.

• Angles in radians, never wrapped during integration.

• The public state is \((s, \vartheta, \varphi, v/v_0, \lambda)\) in reference coordinates; \(s\) stays normalised toroidal flux at every interface, also for chartmaps where the coordinate is \(\rho\).

• Fields and canonical maps are built before the parallel region and read only inside it; per-marker state lives in tracer_t or threadprivate variables.

• No procedure pointers or class() dispatch inside GPU kernels; the CPU path is unchanged when acceleration is off.

Subsystem invariants for tests are listed in Testing. Known design limitations are tracked in issues #288 (VMEC inversion), #290 (sqgBctr), #291 (radial scaling), #292 (global magfie) and #293 (module-level canonical coordinates).