Tutorial 1: Loss fraction from a VMEC equilibrium¶

This tutorial computes the collisionless alpha-particle loss fraction of the quasi-axisymmetric (QA) and quasi-helically symmetric (QH) reactor-scale configurations of Landreman and Paul (2022), and its statistical error. Both VMEC equilibria are public simsopt test files with a mean field of about 5.9 T. Runtime: one to two minutes per configuration on 16 threads.

Step 1: Input data¶

mkdir -p tutorial1/qa tutorial1/qh
cd tutorial1
wget https://github.com/hiddenSymmetries/simsopt/raw/master/tests/test_files/wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc \
    -O qa/wout.nc
wget https://github.com/hiddenSymmetries/simsopt/raw/master/tests/test_files/wout_LandremanPaul2021_QH_reactorScale_lowres_reference.nc \
    -O qh/wout.nc

Step 2: Input file¶

Write tutorial1/simple.in and copy it to both directories (cp simple.in qa/ && cp simple.in qh/):

&config
netcdffile = 'wout.nc'        ! VMEC equilibrium
isw_field_type = 2            ! Boozer canonical coordinates
integmode = 1                 ! symplectic Euler
trace_time = 1d-1             ! s
sbeg = 0.6d0                  ! start surface, normalised toroidal flux
ntestpart = 256               ! number of alpha particles
deterministic = .true.        ! reproducible sampling
output_results_netcdf = .true.
/

trace_time = 0.1: Without collisions (the default) the result is the collisionless loss fraction up to this time. Losses of an optimised configuration occur mostly within the first 10 to 100 ms; the slowing-down time of 3.5 MeV alphas in a reactor is 0.1 to 1 s.
sbeg = 0.6: Markers start on the surface \(s = 0.6\) along one field line (startmode = 1, the default) with isotropic pitch. Larger sbeg probes the outer region, where losses are larger. Other start distributions: Running SIMPLE.
ntestpart = 256: Sets the statistical error (below). Production runs use several thousand markers.

Step 3: Run¶

(cd qa && ../../build/simple.x simple.in > simple.log)
(cd qh && ../../build/simple.x simple.in > simple.log)

The relative path assumes tutorial1 is in the repository root. OMP_NUM_THREADS sets the number of threads.

Step 4: Loss fraction¶

The last line of confined_fraction.dat holds time, the confined fractions of passing and trapped markers, and the number of numerically resolved markers:

tail -n 1 qa/confined_fraction.dat
tail -n 1 qh/confined_fraction.dat

0.10000000000000001  0.70703125000000000  0.17187500000000000  256
0.10000000000000001  0.67968750000000000  0.32031250000000000  256

In QA the confined fraction is \(0.707 + 0.172 = 0.879\): 12.1 % of the markers are lost within 0.1 s. In QH none is lost.

orbit_exit_code.dat gives the reason each orbit ended (second column; all codes in Outputs):

awk '{print $2}' qa/orbit_exit_code.dat | sort | uniq -c

0
1
3

Code 0 is traced to trace_time, code 1 crossed the last closed flux surface, code 3 was not traced: markers with trapping parameter below contr_pp (default \(-1\)) are strongly passing and counted as confined without tracing. Set contr_pp = -1d10 to trace all markers. Codes of 101 and above are numerically unresolved orbits; they are excluded from the confined fractions and reported in unresolved_fraction.dat.

Step 5: Statistical error¶

Each marker is lost or confined independently, so the number of losses \(k\) out of \(N\) resolved markers is binomial. The estimate \(\hat f = k/N\) has standard error

\[\sigma_f = \sqrt{\hat f (1-\hat f)/N}.\]

This plug-in error vanishes at \(k = 0\) and then says nothing. For zero observed losses the one-sided 95 % upper bound is

\[f < 1 - 0.05^{1/N} \simeq 3/N ,\]

which for \(N = 256\) is 1.2 %. A loss fraction quoted as zero means “below \(3/N\)”, not zero. For small \(k\) use the Wilson interval (Plotting results).

Save as tutorial1/plot_confined.py and run it in tutorial1:

import matplotlib.pyplot as plt
import numpy as np

for case in ("qa", "qh"):
    t, passing, trapped, resolved = np.loadtxt(
        f"{case}/confined_fraction.dat", unpack=True)
    loss = 1.0 - (passing + trapped)
    n = resolved[-1]
    error = np.sqrt(loss[-1] * (1.0 - loss[-1]) / n)
    bound = 1.0 - 0.05 ** (1.0 / n)
    print(f"{case}: loss {loss[-1]:.3f} +- {error:.3f}"
          f" (zero-loss 95% bound {bound:.3f})")
    plt.semilogx(t, loss, label=case.upper())

plt.xlabel("time (s)")
plt.ylabel("loss fraction")
plt.legend()
plt.savefig("loss_fraction.png", dpi=150)

qa: loss 0.121 +- 0.020 (zero-loss 95% bound 0.012)
qh: loss 0.000 +- 0.000 (zero-loss 95% bound 0.012)

The QA/QH difference (12.1 ± 2.0 % against < 1.2 %) is significant. The plot shows when the losses happen.

Step 6: Loss times of individual markers¶

times_lost.dat has one line per marker (columns in Outputs): index, loss time in s, trapping parameter \(\kappa\), start \(s\), … Confined markers have loss time trace_time, skipped markers -1. The five earliest losses in QA:

awk '$2 > 0 && $2 < 0.1' qa/times_lost.dat | LC_ALL=C sort -g -k2 | head -5

\(\kappa\) (third column) is negative for passing and in \((0,1)\) for trapped markers (Guiding-centre model). In this QA run 28 of the 31 lost markers are trapped, and the first losses occur after about 20 µs.

What next¶

Vary sbeg, ntestpart and trace_time and watch the loss fraction and its error.
Add collisions with swcoll = .true. and background profiles (Running SIMPLE).
Tutorial 2: Wall heat loads turns the QA results.nc into a heat-load pattern.