Resolution, Step and Solver Guidance

What to set, in what order, and how to check it. Values are starting points for Cyclone-like cases; the check at the end of each section is the test. Measurements are in Numerics, Solvers and Examples.

Defaults at a glance

What you get when you change nothing, and when to change it.

Choice	Default	Change it when
Time stepper	no dt (deck): nonlinear rk4-pi; nonlinear kinetic electrons at \(\beta \le 3\times10^{-4}\) imex-ark3 (adaptive, falls back to rk4-pi with a warning if its factor check refuses); linear and bare TimeConfig() rk4 + CFL. A dt gives rk2 at that step	a published step (rk3 at fixed dt); the speedups are in the table below
Linear growth rates	solver auto: banded shift-invert, else adaptive, certified	one mode of a hard branch (shift-invert)
ky layout	half (1 + Ny//2 rows)	full only to read a two-sided array from an older run
Velocity space	linear (12, 24), nonlinear (4, 8); with kinetic electrons per species: linear ions (12, 32), electrons (6, 48), nonlinear (4, 16) for both species	a deck or call that names Nl/Nm; TEM, KBM: see the sections below
Collisions	linearized Landau; hypercollisions on the last Hermite moments	nu = 0 collisionless runs keep the hypercollisions
End damping	one rate, damp_ends_rate, on every route	never with the step size; see Numerics
Precision	float32	float64 for eigensolves, conditioning and residual checks
Threads and devices	XLA’s pool; Eigen pool off below four CPUs; one device	nothing on CPU; GPUs: XLA_PYTHON_CLIENT_PREALLOCATE=false when shared

Default stepper: measured cost

Deck with no dt and no method against rk4 + CFL, one A4000, float32, wall time including setup, heat flux over the saturated window (transient dropped; error is \(\sigma/\sqrt{n_\mathrm{eff}}\) with the Sokal correlation time, at least 20 correlation times).

Run	Default	Wall, default / rk4 + CFL	\(Q\) default / rk4 + CFL
Cyclone, adiabatic, 64x64, (4, 8), t = 1000	rk4-pi	103 s / 135 s (1.3x)	8.51 +- 0.14 / 8.87 +- 0.15
QHS, adiabatic, 48x48, (4, 8), t = 600	rk4-pi	111 s / 107 s (none)	9.82 +- 0.10 / 9.91 +- 0.09
Cyclone, kinetic electrons, \(\beta\) 1e-4, 32x32, (4, 8), t = 400	imex-ark3	450 s / 1434 s (3.2x)	1129 +- 14 / 1068 +- 18 (rk4-pi 1162 +- 11 in 1123 s)
Kinetic electrons, \(\beta\) 1e-3, t = 60	rk4-pi; imex-ark3 measured 188 s against 129 s (0.7x)	(imex) 188 s / 129 s	(imex) 905 +- 50 / 919 +- 72
Kinetic electrons, \(\beta\) 0.015, t = 60	rk4-pi	182 s / 310 s (1.7x)	18900 +- 800 / 19500 +- 1000

Setup (compile and factor build, 30 to 130 s) is a fixed part of each wall time, so longer runs gain more. The linear eigensolver auto is the certified route but not faster than the time route at (4, 8) or (12, 24) on one GPU (23 s against 10 s and 44 s against 10 s, Cyclone ky 0.3); with kinetic electrons the time route without a dt does not yet integrate (non-finite), so use the eigensolver or give a dt.

Choose the time stepper

GKX advances \(\partial_t g = L g + N(g)\). Write \(L = S + R\), where \(S\) holds the stiff terms (parallel streaming with its dissipation, and the field response) and \(R\) the drifts, mirror, drive and collisions. \(N\) is the \(E\times B\) bracket.

Scheme	Use it for	Limit
rk4 + CFL (linear default without dt)	any deck; the controller sets the step from the stiffest explicit term	step follows electron streaming when electrons are kinetic
rk4-pi (nonlinear default without dt)	RK4 with an embedded error estimate and an I controller	1 to 2x fewer steps than rk4 + CFL at the same flux
imex-ark3 (adaptive, nonlinear kinetic electrons at low beta)	chunk-level adaptive step, no electron-streaming limit	refuses decks whose factor check fails (some stellarators); no gain at \(\beta \gtrsim 10^{-3}\)
rk2/rk3 at fixed dt	reproducing a published step, parity runs	warns when dt exceeds the CFL estimate
imex, imex2	linear runs with stiff collisions or hypercollisions	diagonal damping only; streaming stays explicit
imex-ars3	nonlinear kinetic electrons at 32x32 and up	fixed dt, chosen by hand; factor memory and build time
solver = "krylov"/"adaptive"	linear growth rates: no time step at all	one mode per solve

Accuracy of the adaptive IMEX step. imex-ark3 keeps the embedded ARK3(2) error of each step below rtol = 1e-3 of the state norm. Against rk4-pi (Cyclone, kinetic electrons, \(m_e/m_i\) = 1/1836, 32x32x16, (4, 16), t = 120, A4000, heat-flux window mean over t = 30-120 with the Sokal standard error of each run, three random seeds):

rtol	Steps	Wall	\(Q_i\) against rk4-pi	\(Q_e\) against rk4-pi
rk4-pi	37200	334-342 s	503, 507, 506	709, 724, 749
1e-2 (before)	960-1090	165-201 s	+34, +32, +30 (+6%)	+47, +39, +17 (+5%)
1e-3	2020-2140	204-222 s	+3, -3, +10 (+0.7%)	+12, +8, -33 (-0.6%)
1e-4 (one seed)	4830	349 s	-20 (-4%)	+11 (+2%)

Each window mean carries a standard error of 5-25, so one seed does not resolve the 1e-2 offset; all three seeds have it, and the 1e-3 offsets scatter about zero. Wall times include the 100-200 s of setup. The first turbulent burst is a trajectory comparison, which only holds until the peak: rk4 runs at CFL 0.9, 0.45 and 0.225 (16x16x16, 1/100) agree to 0.1% through the onset and peak, then differ by 10-30% by t = 26-30 (chaotic divergence), so errors of 40% after the peak are not a property of the step. Up to the peak, against rk4 CFL 0.45 (window 12-22) or rk4-pi (32x32x16, window 10-20), the error of \(Q(t)\) is:

Case	1e-2	1e-3	1e-4	3e-5
16x16x16, 1/100: error, steps to t = 30	6.7%, 384	4.7%, 544	3.9%, 1088	0.2%, 1760
16x16x16, 1/1836: error, steps to t = 25	4.2%, 256	8.4%, 288	0.4%, 608	n/a
16x16x16, 1/3670: error, steps to t = 25	3.0%, 288	4.0%, 320	1.0%, 672	n/a
32x32x16, 1/1836: error, steps to t = 22	5.2%, 224	5.2%, 256	2.6%, 416	0.4%, 704

The onset time (\(Q\) = 1e2 to 2e4) is within 0.05 of the reference (t = 12-20) at every tolerance, and a QA stellarator (16x16x32, t = 60) agrees over its first burst to 0.3% at 1e-2 and 0.1% at 1e-4. The peak error plateaus for 1e-2 to 1e-3 and drops only at 1e-4 and below, where the saturated phase takes 5x the steps of 1e-2, so 1e-3 is the default: it is the loosest tolerance whose seed-mean saturated flux agrees with rk4-pi. The bracket constant (dt bound 17 / w_nl) does not change the first-burst error: 2.5, 4 and 8.5 at 1e-2 give 5.3%, 5.2% and 5.2% (320, 256, 224 steps; 224 at 17), and 8.5 and 34 at 1e-4 repeat the result of 17 (2.6%, 416 steps). A bound of 4 costs 2.6x the steps in the saturated phase (one seed) for about the flux offset of 1e-3, so it stays 17. The error is not float32: the same 16x16x16, 1/100 runs in float64 give 14.7% against 14.5% at 1e-2 (window 12-23). With a fixed step imex-ars3 at 0.05 is 0.2% off at t = 20 and at 0.02 it is within 0.1% to t = 21 (0.6% at 22), while 0.1 and 0.15 go non-finite at the burst, so a fixed step needs the margin the adaptive controller gives.

Why IMEX helps. With kinetic electrons the explicit step is set by electron streaming, \(\Delta t \lesssim 1/(v_{te}k_{z,\max})\), about \(\sqrt{m_i/m_e}\approx 60\) times shorter than the ion scale. In imex-ars3 the operator \(S\) is inverted exactly per twist-shift chain (block Thomas plus a Woodbury correction for the fields), so \(\Delta t\) is set by \(R\) and \(N\). On Cyclone with kinetic electrons at 64x64x24 this gave 10x shorter time to solution than RK3 at the same heat flux. At a 16x16 tutorial grid or a small QA deck it does not pay: the factor build dominates.

Choosing the IMEX step. dt is not found by a controller. 0.05 is stable at 32x32x16 and non-finite at 64x64x24, where 0.035 and 0.025 work. Rerun at half the step and compare the heat flux. Set XLA_PYTHON_CLIENT_PREALLOCATE=false on a GPU.

Collisions in the IMEX step. Explicit ARK3 is stable to \(\rho(C)\,\mathrm{d}t \approx 2.5\), with \(\rho(C)\) the spectral radius of the Landau matrix: \(\nu_e\) times 2.9, 8.7, 16 and 24 at \((N_l, N_m)\) = (2, 4), (4, 8), (4, 16) and (8, 16). With collision_scheme = "implicit" (the default) GKX measures \(\rho(C)\) at startup and, only if \(\rho(C)\,\mathrm{d}t > 2\), takes collisions out of the RHS and advances them by a backward-Euler sub-step after each step (dense per point, field-particle part by a Woodbury correction, density, momentum and energy conserved to rounding). Otherwise nothing changes: the run is the explicit one. At dt 0.05 the switch is \(\nu_e \approx 4.6\) at (4, 8), \(\nu_i \approx 0.08\) for \(m_e/m_i = 2.7\times10^{-4}\). Below it the sub-step costs 10% per step on small grids and buys nothing; above it explicit goes non-finite at fixed dt, and an adaptive run is held to \(\mathrm{d}t \lesssim 2.5/\rho(C)\). The sub-step is first order (Lie): its splitting error is outside the embedded estimate, so with fixed_dt = false \(\mathrm{d}t\) is also capped at \(5/\rho(C)\). Measured on 8x10x16 flux tubes, \(\beta\) 1e-4 unless stated, (4, 8), ARK3, 15 time units; error is the heat flux at the end against adaptive rk4 (CFL 0.3), the electrostatic row against explicit at dt 0.0125:

Deck	\(\nu_e\)	dt	\(\rho(C)\mathrm{d}t\)	\(Q\) error, explicit	\(Q\) error, implicit	wall, explicit / implicit
Cyclone, linked	0.6	0.2	1.0	0.0%	0.2%	41 / 46 s
Cyclone, linked	3	0.05	1.3	0.0%	0.7%	74 / 75 s
Cyclone, linked	3	0.2	5.3	non-finite	1.7%	– / 45 s
Cyclone, linked	12	0.05	5.3	non-finite	2.4%	– / 76 s
Cyclone, linked	12	0.2	21	non-finite	4.9%	– / 56 s
Cyclone, \(\beta\) 1e-2	12	0.05, 0.2	5.3, 21	non-finite	2.1%, 5.1%	– / 65, 38 s
Cyclone, electrostatic	12	0.05, 0.2	5.3, 21	non-finite	1.7%, 2.9%	– / 68, 47 s
Cyclone, periodic	12	0.05, 0.2	5.3, 21	non-finite	0.1%, 0.4%	– / 70, 55 s
QA stellarator	12	0.05, 0.2	5.3, 21	non-finite	3.2%, 8.7%	– / 88, 62 s
W7-X	0.6	0.2	1.0	1.3%	3.0%	60 / 67 s
W7-X	12	0.0125, 0.025	1.3, 2.7	0.1%, –	2.7%, 5.7%	257 / 276, 144 s
W7-X	12	0.05, 0.2	5.3, 21	non-finite	10%, 26%	– / 97, 61 s

The growth rate of the Cyclone decks is within 0.35% at every row. The error grows with \(\rho(C)\,\mathrm{d}t\), as \(\sqrt{\mathrm{d}t}\) on the growing Cyclone, QA and periodic modes and about linearly on the transient W7-X flux, so a flux that matters at a few percent at high \(\nu_e\) wants dt near the explicit bound or collision_scheme = "explicit". Adaptive nonlinear Cyclone, \(\beta\) 1e-4, 16x16x16, (4, 8), t = 24 on two cores: \(\nu_e\) = 0.6 explicit 411 s and implicit 366 s, first-burst flux within 0.3% (the default is the explicit run); \(\nu_e\) = 12 explicit 979 s at dt 0.025, default (implicit, dt 0.05 at the cap) 490 s and 2.2% in the first burst. Adiabatic-electron decks are not affected: the chain IMEX factor check refuses them with or without the sub-step.

Resolution

• Velocity space. Linear runs default to \((N_l, N_m) = (12, 24)\), nonlinear runs to (4, 8). Vary \(N_m\) first. Kinetic electrons and TEM need more: TEM growth rates still change at (16, 48).

• Per species. A kinetic-electron deck that names no resolution gives the ions and the electrons their own counts; see Velocity space per species.

• Perpendicular. Nonlinear production runs use 64-96 points in \(x\) and \(y\); gkx --estimate prints the grid the geometry asks for, with the reason.

• Parallel. \(N_z\) = 24-64 per field-period segment, and nperiod >= 2 for any linear run with linked boundaries.

• Check. The heat flux or growth rate changes less than your tolerance when \(N_m\) doubles.

Velocity space per species

Ions and electrons converge at different moment counts, and one count for both pays the larger for each. [[species]] Nl and Nm (inputs) give each species its own, and the right-hand side runs each at that size, so the cost is the sum over species of \(N_l N_m\) rather than \(N_s \max N_l \max N_m\).

What each species needs (Cyclone s-alpha, \(m_e/m_i = 2.7\times10^{-4}\), \(\beta = 10^{-4}\), \(N_z = 32\), certified Krylov eigenvalue at \(k_y\rho_i = 0.3\); the other species held at (12, 48); reference uniform (12, 64), \(\gamma = 0.1394\)):

Electron \(N_l\) (\(N_m = 48\))	1	2	4	6	12
\(\gamma\) against reference	-26%	-13%	-1.1%	+0.2%	+0.3%
Electron \(N_m\) (\(N_l = 4\))	16	24	32	48
\(\gamma\) against reference	+5.1%	+4.6%	+1.9%	-1.1%
Ion \(N_m\) (\(N_l = 12\))	8	12	16	24	32
\(\gamma\) against reference	+11%	-3.4%	-4.7%	+1.8%	+1.7%
Ion \(N_l\) (\(N_m = 48\))	4	6	8	10	12
\(\gamma\) against reference	-22%	+7%	-2%	-5%	+0.3%

Trapped-electron physics needs a few Laguerre moments and many Hermite ones; the ions need Laguerre moments for the finite Larmor radius and fewer Hermite ones. This branch keeps a scatter of a few percent for ion \(N_l\) between 6 and 10 even at (12, 48): ITG and trapped-electron branches sit close at \(k_y\rho_i = 0.3\), so judge a choice by the table below, not by one entry.

Fallback and check. Ions (12, 32), electrons (6, 48) against the uniform (12, 64) reference (uniform (12, 48) is within 0.3% on the ITG point):

Deck	Reference \(\gamma\)	Per species	Difference
Cyclone ITG, \(k_y\rho_i = 0.3\)	0.13940	0.14161	+1.6%
Cyclone ITG, \(k_y\rho_i = 0.2\)	0.0807 (12, 48)	0.0798	-1.1%
TEM, \(a/L_{n}\) 2, \(a/L_{Te}\) 4, \(k_y\rho_i = 0.4\)	0.4492 (12, 48)	0.4495	+0.1%
ETG-range, \(k_y\rho_i = 30\)	19.381	19.514	+0.7%
Landreman-Paul QA stellarator, \(k_y\rho_i = 0.3\)	0.22193	0.22217	+0.1%

The reference moves by 5% between (8, 32) and (12, 48) on the TEM deck, so use no smaller counts there. The cost is \(12\cdot32 + 6\cdot48 = 672\) moments per grid point against 1152 for uniform (12, 48) and 1536 for (12, 64). Smaller counts (ions (8, 24), electrons (4, 32)) are within 1% on the ITG and stellarator decks and 12% low on the TEM deck.

Nonlinear: (4, 16) for both species is the fallback. Against uniform (8, 32) the heat-flux window mean (t = 50-100, Cyclone with kinetic electrons, 32x32x16, imex-ars3, one realization) is:

Counts	\(Q_i\)	\(Q_e\)
uniform (4, 8) (the old fallback)	+12.5%	+9.5%
uniform (4, 16)	-0.4%	+2.4%
ions (4, 8), electrons (4, 16)	+5.6%	+4.2%
ions (8, 24), electrons (4, 32)	-0.8%	+2.0%
ions (8, 16), electrons (4, 32)	-0.2%	+3.1%

One window mean moves by 5-10% with the window start, so differences below that are not resolved; the old (4, 8) is outside it.

Cost of per-species counts. One A4000, float32, 64x64x24, explicit rk3 at a fixed step (time per step from a long and a short run), peak device memory:

Counts	ms per step	Peak memory
uniform (4, 16)	34.9	0.89 GB
ions (4, 8), electrons (4, 16)	31.3 (1.12x)	0.67 GB
uniform (8, 32)	161.3	3.29 GB
ions (8, 24), electrons (4, 32)	110.3 (1.46x)	2.32 GB
ions (8, 16), electrons (4, 32)	84.4 (1.91x)	2.26 GB

One right-hand side at (12, 48): 147.6 ms and 5.5 GB uniform, 86.5 ms and 3.1 GB with ions (12, 32) and electrons (6, 48). Small states are launch-bound (two groups launch twice the kernels), so the gain grows with the state. The state arrays keep the [run] extents, which sets the floor on memory.

Implicit (IMEX) factor. imex-ars3/imex-ark3 build the block-Thomas factor, its field coupling and the collision sub-step per species group, at the counts each species keeps (zero rows above them stay the identity), so memory and solve time follow the counts. One A4000, float32, state at the [run] counts, electrons (4, 32), factor build plus one solve and one collision sub-step:

Grid, ions	Factor, padded / per group	Peak memory	Solve, ms	Build, s
32x32x16, uniform (4, 32)	0.93 / 0.93 GB	1.29 / 1.29 GB	24.5 / 24.0	63 / 62
32x32x16, ions (4, 16)	0.93 / 0.85 GB	1.30 / 1.22 GB	25.5 / 21.0	65 / 63
32x32x16, ions (8, 16)	3.22 / 0.96 GB	5.56 / 1.61 GB	48.2 / 22.2	75 / 68
64x64x24, ions (4, 16)	8.99 / 8.16 GB	10.65 / 9.98 GB	137 / 117	183 / 136
64x64x24, ions (8, 16)	out of memory (16 GB)	12.33 GB	283	166

The electron mirror blocks (Laguerre-coupled) dominate the factor, so the gain is large when the state’s Nl is set by the ions. Uniform counts give the same factor to round-off (solve 2e-14 relative, collision sub-step bitwise). On a loaded 10-core CPU the 32x32x16 ions (8, 16) build took 254 s against 3819 s padded (16 against 9 GB resident). The state and restart arrays keep the [run] extents (zeros above a species’ count). Hermite-sharded routes, a custom collision operator and the per-term diagnostics refuse unequal counts with NotImplementedError.

Kinetic electrons

• Leave fixed_dt = false with an explicit scheme; read the bound in diagnostics.cfl_scales.

• Use a small finite \(\beta\) with \(A_\parallel\) (\(\beta = 10^{-4}\)). At \(\beta = 0\) the electrostatic \(\omega_H\) mode cuts the step a further 6 times.

• For growth rates use solver = "krylov" or "adaptive".

• Time: allow \(10/\gamma\) for the linear phase, saturation by \(t\approx 60\)-100, then 100-200 for the average, or run_to = "saturation".

• Check. Rerun with damp_ends_rate 3x up and down and with nperiod + 1 (below).

Electromagnetic and KBM

• \(\beta\) in the deck sets \(A_\parallel\); check the mode against the next \(N_m\).

• The linear KBM branch moves with the end-damping rate and with nperiod. At nperiod = 1 the layer is a quarter of the line: use 2 or more.

• The Miller KBM needs more points along the line than the reference deck’s ntheta = 32. At \(\beta=0.015\), \(k_y=0.3\), nperiod = 2, (16, 48), \(\gamma\) is 0.3128, 0.2999, 0.2953 and 0.2920 at ntheta 32, 48, 64 and 96, falling as \(1/n_\theta^2\); use 64 or more and rerun at 1.5x. Nm 96 moves it by 0.5%. The comparison with GS2 is in Validation (X-kbm-beta-gs2).

• A time run that stops without a stationary growth window falls back to a log-linear window of at least 30% of the samples and warns. Read the slope of the late trace before trusting a rate from such a run; an electron-only start (init_electrons_only) shows a faster transient for the first 10 to 20 time units.

• Explicit steps must keep \(\nu_{end}\Delta t < 0.3\). An nperiod = 1 KBM deck (rate 598) blew up at \(\Delta t = 10^{-3}\) and ran at \(5\times10^{-4}\).

• Eigenvalues: auto (banded, else adaptive) certifies KBM; a shift-invert with a poor shift can certify a non-dominant branch.

End damping

One rate, [time] damp_ends_rate, applies on every route and no time step enters it. The default is 30 e-folds per thermal transit of the fastest kinetic species across the end layer, and every run prints it. Below about 1 e-fold the kinetic-electron and KBM rates are not converged. Rerun at 3x above and below the printed rate and at nperiod + 1; a converged rate moves by less than its tolerance in all three (Parallel boundary and end damping).

Collisions

The default operator is the linearized Landau operator (Operators And Terms). It needs the velocity resolution of the case; a factor is built once (0.5 s at (4, 8), tens of seconds at (16, 48)) and cached. Set nu per species; lenard_bernstein is cheaper and not accurate at high \(\nu\) (21% low on a Cyclone scan at \(\nu_i = 0.1\)).

Eigensolver

auto (default) runs the Hermite-banded direct shift-invert where its blocks fit in memory and adaptive otherwise; adaptive is Arnoldi on the time-advance map, as GS2 does with SLEPc. Both certify the pair against the full operator. Use matrix-free shift-invert only for a chosen interior branch (Methods and Decisions).