Stellarator FCI Validation And Reduced SOL Dynamics

Figures are release-hosted

The campaign figures referenced below live on the validation-artifacts-2026-04-28 GitHub release. Because this repository is private, release assets cannot render inline for docs readers; the links require repository access, or restore them locally with python scripts/fetch_example_artifacts.py.

Figures are release-hosted

The campaign figures referenced below live on the validation-artifacts-2026-04-28 GitHub release. Because this repository is private, release assets cannot render inline for docs readers; the links require repository access, or restore them locally with python scripts/fetch_example_artifacts.py.

This page documents the first native non-axisymmetric 3D scrape-off-layer lane in drbx. The lane is deliberately staged: it validates geometry and field-line-map operators first, then runs a compact reduced dynamics benchmark on the same geometry. It is not yet a real-device prediction.

The generated artifacts live in docs/data/stellarator_fci_validation_artifacts/ and are regenerated by:

PYTHONPATH=src .venv/bin/python examples/geometry-3D/stellarator-fci/validation_campaign.py

For a step-by-step, user-editable version of the same workflow, start with Stellarator examples. That page maps each README figure to a runnable script, links the source code functions used by the examples, and separates the linear mode, nonlinear reduced-turbulence, and geometry plotting cases. For the distinction between synthetic connection-length proxies, imported wall-hit lengths, adjacent-plane arc lengths, and effective \(L_\parallel\) closures, see Connection Length.

The focused regression gate is:

.venv/bin/pytest -q tests/test_geometry_fci_maps.py tests/test_validation_stellarator_fci_campaigns.py

The current suite includes three analytic 3D configurations so the native lane does not validate against a single geometry only: a baseline island SOL, a stronger island-map displacement case, and a higher mirror/shear case.

Literature Pattern

The validation pattern follows current stellarator edge/SOL literature rather than treating the first 3D movie as sufficient evidence. Recent work emphasizes field-line-map geometry, island and open-field connection length, curvature drive, fluctuation statistics, spectra, radial flux, and metric/operator verification before making device-scale turbulence claims.

The main public references used to define the validation checklist are:

  • Coelho et al. 2022, for global flux-driven stellarator island-divertor analysis and coherent mode/transport diagnostics.
  • Shanahan et al. 2024, for SOL fluctuation amplitude, skewness, island-width-scale structures, and connection between curvature drive and island-region transport.
  • The 2024 Journal of Plasma Physics article, for methods-section structure: model equations, geometry, grid, source, boundary conditions, and then steady/fluctuating results.
  • The 2025 non-axisymmetric edge-fluid solver paper, for the expected standard of metric verification, island profile tests, target-aware geometry, 3D visualization, and continuous verification.

The local reference-workflow QA added one practical plotting constraint: the most readable stellarator SOL media use R-Z cross-sections at several toroidal angles, common color limits, explicit time/angle labels, opened traced surfaces, radial cuts, and field-line overlays. The current figures now follow that pattern instead of showing only a radial-poloidal heatmap or a single pale outer torus surface.

Implemented Objects

The current lane is split into small, directly tested modules:

  • MetricTensor3D stores \(J\), \(|B|\), all six independent contravariant metric components, and all six independent covariant metric components.
  • build_metric_report checks finite values, positive \(J\), positive \(|B|\), and inverse consistency \(g^{ik}g_{kj}=\delta^i_j\).
  • FciMaps stores forward/backward traced plane intersections, boundary masks, and the toroidal plane spacing.
  • build_synthetic_stellarator_geometry creates an analytic five-period non-axisymmetric geometry with rotating elliptical cross-sections, helical mirror modulation, island-like map displacement, curvature proxy, and connection-length proxy.
  • interpolate_fci_plane, fci_yup, fci_ydown, grad_parallel_fci, laplace_parallel_fci, and laplace_perp_xz provide the first JAX-native field-line-map operators.
  • build_fci_target_masks, fci_bohm_sound_speed, and compute_fci_sheath_recycling provide the first JAX-native target closure on non-axisymmetric traced field lines: endpoint masks, normalized Bohm ion loss, zero-current electron particle reconstruction, sheath heat loss, and exact recycled-neutral source accounting.
  • conservative_parallel_diffusion_fci and conservative_perp_diffusion_xz provide the first conservative metric-weighted transport kernels, with constant-state and dissipation gates.
  • metric_weighted_scalar_laplacian_3d evaluates the full \(J^{-1}\partial_i(JK g^{ij}\partial_j f)\) scalar diffusion operator using all contravariant cross terms. The corresponding manufactured-solution campaign verifies the JAX-native stencil on an analytic identity metric and checks constant-state, dissipation, and cross-term activity on the synthetic stellarator metric.
  • compute_fci_neutral_reaction_diffusion adds neutral diffusion, ionisation, recombination, and charge-exchange source accounting on the same geometry.
  • apply_fci_vorticity_operator and solve_fci_vorticity_potential_cg add the first non-axisymmetric metric-weighted vorticity/potential inversion, radial \(E\times B\) proxy, and Boussinesq versus non-Boussinesq polarization gate. The solver default remains fixed-iteration CG; an opt-in preconditioner="jacobi" path applies a conservative diagonal approximation to the perpendicular operator and is gated by a manufactured residual-reduction test before it is used for larger imported-field movie residual probes.
  • FciDrbState and compute_fci_drb_rhs assemble the first transformable PyTree RHS surface combining target, neutral, and vorticity components. The compact RHS now threads the same Boussinesq/non-Boussinesq potential-solve switch used by the vorticity component gate and an opt-in potential-fed plasma \(E\times B\) advection path. That advection path acts on charged density, pressure, ion parallel momentum, and vorticity; neutral gas fields remain controlled by neutral diffusion and reaction closures. The focused regression suite checks that this combined state supports jax.jvp, jax.vmap, and multi-device execution when more than one local device is available.

The design rule is that geometry data and field-line maps are ordinary Python objects holding arrays that can be consumed by JAX kernels. The next refactor should make these objects explicit PyTrees so full reduced models can be batched, differentiated, and compiled without host-side orchestration.

Equations

Let \(x^1=x\), \(x^2=y\), \(x^3=z\), and let \(\mathbf{r}(x,y,z)\) be the physical position. The covariant basis and metric are

e_i = partial r / partial x^i
g_ij = e_i dot e_j
J = e_1 dot (e_2 cross e_3)
g^ij = inverse(g_ij)

The metric gate verifies:

min(J) > 0
min(|B|) > 0
all fields finite
max |g^ik g_kj - delta^i_j| < tolerance

For a scalar field \(f\), field-line-map interpolation defines:

f_up(i,j,k) = I[f(:,j+1,:)](x_fwd(i,j,k), z_fwd(i,j,k))
f_dn(i,j,k) = I[f(:,j-1,:)](x_bwd(i,j,k), z_bwd(i,j,k))

The first compact native operators are:

partial_parallel f ~= (f_up - f_dn) / (2 dphi)
L_parallel f ~= (f_up - 2 f + f_dn) / dphi^2
L_perp f ~= partial_xx f + partial_zz f

These are not the final production operators for an imported equilibrium. The production operator should be the metric-weighted conservative form:

div_parallel(K grad_parallel f)
  = J^(-1) partial_parallel(J K partial_parallel f)

with boundary-distance fallback where the traced line leaves the domain before a full stencil is available.

The full 3D metric gate now verifies the scalar form used for perpendicular and non-axisymmetric diffusion terms:

L_K f = J^(-1) partial_i(J K g^ij partial_j f)
F^i = J K g^ij partial_j f
L_K f = J^(-1) partial_i F^i

Here \(K\) is a scalar transport coefficient, \(g^{ij}\) includes the non-orthogonal cross terms \(g^{12}\), \(g^{13}\), and \(g^{23}\), and summation over repeated logical-coordinate indices is implied. The present implementation uses centered differences with periodic toroidal/binormal axes and one-sided radial boundaries by default. The manufactured-solution test uses fully periodic axes on an identity metric so the exact source is known, while the synthetic stellarator probe checks that constants are annihilated and that cross terms contribute nontrivially on a non-axisymmetric metric.

The reduced dynamics benchmark evolves a scalar fluctuation field:

partial_t n_tilde
  = chi_parallel L_parallel(n_tilde)
  + chi_perp L_perp(n_tilde)
  + alpha C(x,y,z) R[n_tilde]
  + N[n_tilde]
  - beta n_tilde
  - gamma n_tilde^3

where \(C\) is the curvature proxy, \(R\) is a radial-gradient drive, \(N\) is a conservative nonlinear transfer proxy, and the last two terms keep the compact benchmark bounded.

The non-axisymmetric sheath/recycling gate is the first closure promoted beyond geometry and scalar dynamics. On every forward/backward traced endpoint it evaluates:

c_s = sqrt((T_e + T_i) / m_i)
Gamma_i,target = N_endpoint n c_s
Gamma_e,target = Gamma_i,target
q_e,target = gamma_e Gamma_i,target T_e
q_i,target = gamma_i Gamma_i,target T_i
Gamma_n,recycle = R_recycle Gamma_i,target
Q_n,recycle = E_recycle Gamma_n,recycle

The current implementation keeps this as a local source-accounting gate rather than a full transient solve. That distinction matters: the figure below proves the traced endpoint masks and source identities, while the next solver step is to route these losses and sources through the full density, pressure, neutral, momentum, and vorticity RHS.

The compact neutral gate evaluates diffusion plus the local reaction block:

S_ion = k_ion n_n n_e sqrt(T_e)
S_rec = k_rec n_i n_e / sqrt(T_e)
S_cx = k_cx n_n n_i sqrt(T_n + T_i)

partial_t n_n |_react = -S_ion + S_rec
partial_t n_i |_react =  S_ion - S_rec
partial_t n_e |_react =  S_ion - S_rec

Charge exchange exchanges ion and neutral momentum and thermal energy while leaving particle number unchanged. The campaign therefore checks particle and momentum conservation separately from the metric-weighted diffusion integral.

The first vorticity gate solves the perpendicular polarization relation

Omega = - div_perp(K_pol grad_perp phi)
K_pol = <n / B^2>      (Boussinesq)
K_pol = n / B^2        (non-Boussinesq)

with a mean-free conjugate-gradient inversion on the metric-weighted perpendicular operator. The Boussinesq form freezes the polarization coefficient at the grid mean used by the implementation, while the non-Boussinesq form keeps the local density and magnetic-field modulation. The campaign now verifies both inversions against the same manufactured potential, checks that the variable-coefficient operator differs from the Boussinesq operator on a nonuniform density field, and checks that both forms collapse to the same discrete operator when \(n/B^2\) is constant. This is the numerical seam needed before replacing the reduced radial-flux proxy with a true \(E\times B\) transport diagnostic. The native nonlinear advection seam now exposes the logical perpendicular bracket used by the pedagogical vorticity example:

{phi, f}_{s,theta}
  = (partial_theta phi partial_s f - partial_s phi partial_theta f) / B.

The corresponding manufactured-solution test verifies convergence on a periodic identity metric, zero bracket for a constant advected field, and vanishing domain-mean bracket to roundoff. This is the physics-backed path for nonlinear density, pressure, and vorticity advection in the example suite. The older scalar nonlinear_transfer used in the showcase movie remains a reduced visualization gate; it should not be cited as a Braginskii vorticity/Poisson-bracket derivation.

Validation Metrics

The geometry campaign currently passes:

The multi-configuration geometry suite currently passes all three metric/map gates. The three cases span mirror ratios from about 0.45 to 0.61, mean connection-length proxies from about 23.8 to 24.7, and maximum radial map shifts from about 0.17 to 0.36 cells:

The operator campaign currently passes:

The full metric manufactured-solution campaign currently passes:

  • identity-metric MMS observed order: about 1.90;
  • RMS error decreases from about 2.00 at 16 x 16 x 32 to about 2.48e-1 at 48 x 48 x 96;
  • synthetic non-axisymmetric constant-state residual: below display precision;
  • synthetic metric energy monotone fraction: 1.0;
  • synthetic metric energy drop fraction: about 4.86e-3;

  • cross-term contribution fraction relative to the full operator: about 0.17.

  • Stellarator full metric MMS validation — release-hosted figure: docs__data__stellarator_fci_validation_artifacts__metric_mms__images__stellarator_metric_mms_campaign.png

The sheath/recycling campaign currently passes:

The neutral physics campaign currently passes:

The vorticity/potential campaign currently passes:

  • Boussinesq relative potential reconstruction error: about 1.30e-3;
  • non-Boussinesq relative potential reconstruction error: about 1.44e-3;
  • Boussinesq relative residual: about 2.04e-4;
  • non-Boussinesq relative residual: about 4.76e-4;
  • variable-coefficient operator difference: about 0.365 relative to the Boussinesq operator;
  • constant-\(n/B^2\) Boussinesq/non-Boussinesq operator mismatch: about 1.66e-10;
  • density-over-\(B^2\) contrast in the manufactured probe: about 3.11;

  • radial \(E\times B\) proxy RMS: about 1.04.

  • Stellarator vorticity validation — release-hosted figure: docs__data__stellarator_fci_validation_artifacts__vorticity__images__stellarator_vorticity_campaign.png

The PyTree DRB campaign currently passes:

  • short fixed-layout 3D transient on the combined ion/electron/neutral density, pressure, momentum, and vorticity state;
  • JVP-versus-finite-difference relative derivative error: about 6.4e-14;
  • non-Boussinesq objective JVP-versus-finite-difference relative derivative error: about 2.9e-9;
  • vmap-versus-serial objective mismatch: about 8.9e-16;
  • final compact potential residual: about 1.32;
  • Boussinesq versus non-Boussinesq potential relative difference on the same state: about 1.04e-1;
  • RHS source-state difference when only the potential model is switched: below display precision when potential feedback is disabled, preserving the old diagnostic-only default;
  • potential-fed ExB gate: the same Boussinesq/non-Boussinesq switch changes plasma RHS terms when plasma_exb_advection_strength > 0, leaves neutral ExB terms absent, and passes a JVP-versus-finite-difference objective gate;
  • density-over-\(B^2\) contrast in the compact PyTree state: about 4.74;
  • local CPU warmed campaign time: about 2.1e-2 seconds for the documentation grid;

  • remote two-device GPU smoke profile: passing, with pmap execution around 7.0e-2 seconds on the same compact grid.

  • Stellarator PyTree/JVP/scaling validation — release-hosted figure: docs__data__stellarator_fci_validation_artifacts__pytree_drb__images__stellarator_drb_pytree_campaign.png

The imported FCI campaign currently passes:

  • external field-line maps converted to fixed-shape FCI endpoints on an scaled VMEC Landreman-Paul QA shell;
  • forward and backward boundary fractions: about 0.67 and 0.75;
  • target endpoint fraction: about 0.90;
  • magnetic-field modulation: about 1.42;
  • scaled-surface non-axisymmetric major-radius RMS: about 0.116;

  • neutral particle and momentum balance errors: below display precision.

  • ESSOS imported FCI validation — release-hosted figure: docs__data__essos_imported_fci_artifacts__images__essos_imported_fci_campaign.png

The VMEC-coordinate imported FCI control closes the field-line map on VMEC surfaces, disables target endpoint masks, and preserves neutral source accounting without sheath heat load. The hybrid imported FCI control uses VMEC-coordinate interpolation locations with coil-derived endpoint masks, connection-length proxy, and \(|B|\), so it is the current open-field SOL bridge for non-axisymmetric sheath/recycling work.

The imported PyTree/JVP campaign currently passes:

The ESSOS-imported QA-coil DRB movie campaign currently passes movie-QA checks, but its report keeps publication_ready=false while coil connection-length, grid-refinement, and time-refinement gates remain open:

The high-resolution Jacobi stationarity media is the preferred current imported-map showcase. It uses the hybrid VMEC/coil map source, the 16 x 96 x 48 physics grid, the Jacobi-preconditioned FCI potential solve, and the same fixed-layout DRB state with sheath, recycling, neutral, and vorticity closures. The paired report-only refinement summaries pass the grid/time gate, and the JSON-only stationarity report passes on the same settings. The media still represents a reduced transient demonstration rather than a final long-time turbulence validation, but it is no longer a renderer-only or coarse-grid movie.

The reduced dynamics campaign currently passes:

Interpretation

The current geometry plot is a metric and map-quality result. It shows that the analytic 3D geometry has nontrivial rotating cross-sections, finite \(|B|\), short- and long-connection-length regions, and nonzero radial/poloidal map displacements. This is the minimum evidence needed before any native operator claim on non-axisymmetric geometry.

The operator plot is the first numerical result. The interpolation curve is close to second order, as expected for bilinear interpolation on smooth fields. The parallel-gradient curve is lower because the field-line map, shifted sampling, and central difference are all active at once. The monotone diffusion energy curve confirms that the compact parallel diffusion probe is dissipative on the tested geometry.

The sheath/recycling figure is the first non-axisymmetric target closure result. The top row shows target particle loss, heat-weighted connection length, and target heat load in toroidal-poloidal space; the bottom row shows radial localization, toroidal modulation, and integrated source-balance bars. The important validation point is not the absolute normalized magnitude, but the fact that particle recycling, neutral energy recycling, and zero-current balance are closed to roundoff on non-axisymmetric endpoint masks.

The neutral figure exercises the next physics layer: recycling sources can now feed a neutral field whose diffusion and local ionisation/recombination/CX block are conservative at the source level. The vorticity figure exercises the potential-solve layer needed for actual \(E\times B\) transport. Together, these two gates turn the non-axisymmetric lane from a geometry/movie scaffold into a sequence of tested solver components.

The PyTree/JVP figure is the first end-to-end transformability result on this lane. It is a component-level transformability result rather than a full implicit transient solve, but it does prove that the fixed-layout component RHS can be compiled, batched, differentiated by forward-mode AD, and dispatched across multiple local devices without changing the numerical objective.

The dynamics plot and movie are a compact 3D showcase. The positive skewness, finite kurtosis, and radial localization metrics show that the field is not a static coloring of a torus. The benchmark is still a reduced scalar model, so the correct claim is that the geometry/operator/sheath/movie pipeline is now native and tested. The next broader physics claim requires carrying the potential-fed ExB path, source-driven steady state, neutral closures, and grid or timestep validation onto imported geometry.

Remaining Gates

The promoted synthetic stellarator lane now covers metric MMS, traced parallel-diffusion conservation, target sheath/recycling closure, neutral reaction-diffusion closure, vorticity/potential residuals, fixed-layout PyTree/JVP transformability, and compact reduced dynamics. The remaining work needed before a long-time stellarator turbulence claim is narrower and more physical:

  1. run longer nonlinear source-driven transients with grid and timestep sensitivity rather than relying on the current compact visualization window;
  2. replace synthetic endpoint masks with wall-resolved target geometry and connection-length metrics on imported maps;
  3. route the same non-axisymmetric closure arrays through a full production PyTree residual and implicit transient solve;
  4. refine imported-map agreement for coil, VMEC-coordinate, and hybrid maps against independent field-line diagnostics;
  5. extend autodiff, UQ, and inverse-design examples from compact geometry parameters to transport/source controls on the full 3D lane.