Thea Energy, Inc. (United States)

GX is a code designed to solve the nonlinear gyrokinetic system for low-frequency turbulence in magnetized plasmas, particularly tokamaks and stellarators. In GX, our primary motivation and target is a fast gyrokinetic solver that can be used for fusion reactor design and optimization along with wide-ranging physics exploration. This has led to several code and algorithm design decisions, specifically chosen to prioritize time to solution. First, we have used a discretization algorithm that is pseudospectral in the entire phase space, including a Laguerre–Hermite pseudospectral formulation of velocity space, which allows for smooth interpolation between coarse gyrofluid-like resolutions and finer conventional gyrokinetic resolutions and efficient evaluation of a model collision operator. Additionally, we have built GX to natively target graphics processors (GPUs), which are among the fastest computational platforms available today. Finally, we have taken advantage of the reactor-relevant limit of small $\rho _*$ by using the radially local flux-tube approach. In this paper we present details about the gyrokinetic system and the numerical algorithms used in GX to solve the system. We then present several numerical benchmarks against established gyrokinetic codes in both tokamak and stellarator magnetic geometries to verify that GX correctly simulates gyrokinetic turbulence in the small $\rho _*$ limit. Moreover, we show that the convergence properties of the Laguerre–Hermite spectral velocity formulation are quite favourable for nonlinear problems of interest. Coupled with GPU acceleration, which we also investigate with scaling studies, this enables GX to be able to produce useful turbulence simulations in minutes on one (or a few) GPUs and higher fidelity results in a few hours using several GPUs. GX is open-source software that is ready for fusion reactor design studies.

Abstract This report presents an initial analysis of the fast ion confinement and losses within quasi-axisymmetric stellarator equilibria in consideration by Thea Energy. The equilibria have not yet been explicitly optimized for fast particle confinement and require validation. Modeling with the ASCOT5 code is used to directly examine the fast ion transport. The particle tracking simulations are purely (neo)classical in nature and simply contain the supplied equilibrium and collisions (pitch-angle, energy slowing, and velocity diffusion) from supplied thermal profiles. Uniform marker deposition is used to probe the general confinement properties of the equilibria while a realistic beam-born population is provided from the BEAMS3D code and an alpha particle population is calculated from a fusion source integrator. A first wall is included and defines the loss boundary. Analysis for NBI ions within Thea Energy’s conceptual Eos neutron source are presented along with alpha particles in an enlarged DT-plasma. The fast ions are assessed in regards to their confinement time, pitch, energy, and spatial coordinates. For each population, the impact of collisions and orbit drifts are discussed. It is found that NBI born ions in Eos are strongly confined until slowing-down, owing largely to the tangential injection geometry, while 22% of DT-born alpha energy is lost in the scaled device, indicating that any fusion pilot plant design optimization should include metrics for fast ion confinement.

To better understand the dependence of the magnetic field structure in the plasma edge on the plasma boundary shape, in the context of X-point and island divertor designs, we define and develop a class of stellarators called umbilic stellarators. These equilibria are characterised by a single continuous high-curvature edge on the plasma boundary that goes around multiple times toroidally before meeting itself. We develop a technique that allows us to simultaneously optimise the plasma boundary along with a curve lying on the boundary on which we impose a high curvature while imposing omnigenity – a property of the magnetic field that ensures trapped particle confinement throughout the plasma volume. We find that umbilic stellarators naturally tend to favour piecewise omnigenity instead of omnigenity with a specific helicity. After generating omnigenous umbilic stellarators, we design coil sets for some of them and explore the fieldline structure in the edge and its sensitivity to small fluctuations in the plasma. Finally, using single-stage optimisation, we simultaneously modify the plasma and coil shape and propose an experiment to modify an existing tokamak to a finite- $\beta$ stellarator using this technique and explore a potentially simpler way to convert a limited tokamak into a diverted stellarator.

We propose to build a Flexible Stellarator Physics Facility to explore promising regions of the vast parameter space of disruption-free stellarator solutions for Fusion Pilot Plants (FPPs). The FESAC Long Range Plan recognized the quasi-symmetric stellarator as “the leading US approach to developing disruption-free, low-recirculating-power fusion configurations”. To deliver the ambitious Decadal Vision for Commercial Fusion Energy, we must establish a persuasive stellarator program in parallel to the tokamak one: the stellarator will prove to be a better path to a reactor if theoretical predictions are confirmed and novel optimization techniques and strategies work as desired. Since the release of the Long Range Plan in 2020, stellarators have arguably made the most significant advances of all fusion concepts. Groundbreaking results from W7-X demonstrated low neoclassical transport and the successful operation of the island divertor. Advances in theory and modeling now allow us to minimize turbulent transport, to achieve equilibria with precise quasisymmetry, to reduce neoclassical transport and fast ion loss to levels far below what has been previously achieved, and to minimize the effect of coil manufacturing errors. If realized, these advances will lead to cost-effective stellarator designs with confinement comparable to tokamaks but without the fundamental challenges of disruptions and current drive.

This patch release updates package dependencies to resolve critical incompatibilities between jax-related libraries. DESC versions 0.16.0 and 0.17.0 contain a known bug when used alongside jax-finufft v1.3.0, and we recommend that users upgrade to DESC version 0.17.1 to ensure environment stability. If an immediate upgrade is not feasible, the issue can be mitigated by manually downgrading jax-finufft to version 1.2.0. Minor Changes Fixes incorrect units in the documentation of some curvature variables. Full Changelog: https://github.com/PlasmaControl/DESC/compare/v0.17.0...v0.17.1

New Features Adds particle tracing capabilities in desc.particles module. Particle tracing is done via desc.particles.trace_particles function. Particles can be initialized in couple different ways: ManualParticleInitializerLab : Initializes particles at given positions in lab coordinates. ManualParticleInitializerFlux : Initializes particles at given positions in flux coordinates. CurveParticleInitializer : Initializes N particles on a given curve. SurfaceParticleInitializer : Initializes N particles on a given surface. Implemented particle trajectory models are: VacuumGuidingCenterTrajectory : Integrates the particle motion by vacuum guiding center ODEs, conserving energy and mu. Particle trajectories can be plotted with desc.plotting.plot_particle_trajectories function. Adds new option for loss_function, "sum", which changes an objective to target the sum of the values computed. Adds utility functions desc.external.paraview.export_surface_to_paraview, desc.external.paraview.export_volume_to_paraview and desc.external.paraview.export_coils_to_paraview to export Paraview files for surfaces, volume and coils. These functions use an optional dependency pyvista which is not automatically installed. The x_scale option for eq.optimize and related functions can now be given as a dictionary mapping individual parameter names to their relevant scales, or if multiple things are being optimized, x_scale can be a list of dict, one for each optimizable thing. Adds new option x_scale='ess' to use exponential spectral scaling from (Jang 2025) which has been shown to improve performance and robustness as an alternative to fourier continuation methods. Adds x_scale='ess' option for the OmigenousField class. Adds "scipy-l-bfgs-b" optimizer option as a wrapper to scipy's "l-bfgs-b" method. The x_scale parameter can now be used with stochastic gradient descent type optimizers. Adds wrappers for optax optimizers. They can be used by prepending 'optax-' to the name of the optimizer (i.e. optax-adam). Additional arguments to the optimizer such as learning_rate can be pass via options = {'optax-options': {'learning_rate': 0.01}}. Even a custom optax optimizer can be used by specifying the method as 'optax-custom' and passing the optax optimizer via the 'update-rule' key of optax-options in the options dictionary. See the docstring of the optax-custom for details. Adds check_intersection flag to desc.magnetic_fields.FourierCurrentPotentialField.to_Coilset, to allow the choice of checking the resulting coilset for intersections or not. Changes the import paths for desc.external to require reference to the sub-modules. Adds a differentiable utility for finding constant offset toroidal surfaces inside of optimizations. See PR for more details. Add support for Python 3.14 Bug Fixes No longer uses the full Hessian to compute the scale when x_scale="auto" and using a scipy optimizer that approximates the hessian (e.g. if using "scipy-bfgs", no longer attempts the Hessian computation to get the x_scale). SplineMagneticField.from_field() correctly uses the NFP input when given. Also adds this as a similar input option to MagneticField.save_mgrid(). Fixes some bugs that hampered robustness of desc.geometry.FourierRZToroidalSurface.constant_offset_surface, particularly when the given grid had stellarator symmetry or when NFP=1. Fixes possible bug in computing normalizations when both kinetic and pressure profiles are assigned. Also adds warnings whenever an pressure is added to a kinetic-constrained equilibrium and vice-versa to alert user to ambiguous equilibrium setups. Adds error in MercierStability to guard against situation where if a grid with a point at rho=0 were used, NaN would be computed, asMercierStability is undefined on-axis. Performance Improvements ProximalProjection.grad uses a single VJP on the objective instead of multiple JVP followed by a manual VJP. This should be more efficient for expensive objectives. Deprecations sgd optimizer is deprecated in favor of optax-sgd, and will be removed in a future release. To achieve the same behavior with optimizer = Optimizer('sgd') and options={'alpha': ..., 'beta': ...} when the optimizer is removed, one can use optimizer = Optimizer('optax-sgd') and options={'optax-options': {'learning_rate': alpha, 'momentum': beta, 'nesterov': True}}. Removes FiniteDiffDerivative from the public API. This class was no longer actually usable with the current versions of DESC's optimization framework, as JAX is now required for running any equilibrium or optimization solves. New Contributors @byoungj made their first contribution in https://github.com/PlasmaControl/DESC/pull/1736 Full Changelog: https://github.com/PlasmaControl/DESC/compare/v0.16.0...v0.17.0

Bootstrap current plays a crucial role in the equilibrium of magnetically confined plasmas, particularly in quasi-symmetric stellarators and in tokamaks, where it can represent bulk of the electric current density. Accurate modeling of this current is essential for understanding the magnetohydrodynamic (MHD) equilibrium and stability of these configurations. This study expands the modeling capabilities of M3D-C1, an extended-MHD code, by implementing self-consistent physics models for bootstrap current. It employs two analytical frameworks: a generalized Sauter model (Sauter et al. 1999 Phys. Plasmas vol. 6 , no. 7, pp. 2834–2839), and a revised Sauter-like model (Redl et al. 2021 Phys. Plasmas vol. 28 , no. 2, pp. 022502). The isomorphism described by Landreman et al. (2022 Phys. Rev. Lett. vol. 128, pp. 035001) is employed to apply these models to quasi-symmetric stellarators. The implementation in M3D-C1 is benchmarked against neoclassical codes, including NEO, XGCa and SFINCS, showing excellent agreement. These improvements allow M3D-C1 to self-consistently calculate the neoclassical contributions to plasma current in axisymmetric and quasi-symmetric configurations, providing a more accurate representation of the plasma behavior in these configurations. A workflow for evaluating the neoclassical transport using SFINCS with arbitrary toroidal equilibria calculated using M3D-C1 is also presented. This workflow enables a quantitative evaluation of the error in the Sauter-like model in cases that deviate from axi- or quasi-symmetry (e.g. through the development of an MHD instability).

In the pursuit of the goal of commercial fusion as an abundant and safe source of energy, the stellarator is a leading concept with compelling attractiveness and demonstrated performance. A new mid-size stellarator is needed to retire risks and innovate towards a high performance, economically attractive, stellarator Fusion Pilot Plant. In this presentation we, as a community of US researchers from Universities, National Laboratories, and Private Industry, involved in studying the stellarator concept, lay out the programmatic and technical motivation for a new and modern mid-size stellarator research facility. A new mid-scale stellarator is needed to realize the potential predicted by a solid body of theory and simulation along with advances in computational tools for optimization and non-linear turbulence modeling. Notably, it is possible to combine the advantages of the stellarator (steady state, no current drive, no disruptions) with the good confinement regularly achieved in tokamaks. The top priorities for experimental work, and the motivation for a mid-scale stellarator experiment are: turbulence control, non-resonant divertor, MHD stability at large beta, confinement of fast particles, and coil simplification. A new mid-size quasi-symmetric stellarator, built as a user facility, would complement existing research at Wendelstein 7-X and LHD and strongly augment private industry. It would provide a program of innovative research, concept validation, theoretical advancement, and workforce development. Growing support and interest for stellarators by the fusion community and private industry affirms this rationale.

New Features New basis vector and metric elements derivatives in PEST coordinates and quantities useful for a global MHD stability solver. Adds desc.external.TERPSICHORE objective for computing linear ideal MHD stability growth rates. This objective subclasses from ExternalObjective and requires access to the TERPSICHORE code, which is not included with DESC or its dependencies. Adds docs/dev_guids/external_objectives.rst as a tutorial for how to use external objectives, with TERPSICHORE as an example using parallel processing. Adds keyword argument normalize to plot_1d, plot_3d. normalize is a string to use for normalization. Adds new linear objective ShareParameters which can enforce that the chosen parameters of two or more objects of the same type remain identical during an optimization. Potentially useful for flexible stellarator optimization, where one has two coilsets with the same geometry but differing currents, and attempts to optimize for two different stellarator equilibria. Changes related to field_line_integrate, see #1839: field_line_integrate now returns diffrax.diffeqsolve.stats and diffrax.diffeqsolve.result if the flag return_aux is set to True. Renames maxsteps argument of field_line_integrate to max_steps. Now the argument has a consistent meaning with the diffrax package and specifies the maximum number of steps allowed for whole integration. Previously, it was used as maximum number of iterations between integration steps. field_line_integrate function doesn't accept additional keyword-arguments related to diffrax, if it is necessary, they must be given through options dictionary. poincare_plot and plot_field_lines functions can now plot partial results if the integration failed. Previously, user had to pass throw=False or change the integration parameters. Users can ignore the warnings that are caused by hitting the bounds (i.e. Terminating differential equation solve because an event occurred.). chunk_size argument is now used for chunking the number of field lines. For the chunking of Biot-Savart integration for the magnetic field, users can use bs_chunk_size instead. Bug Fixes Fixes straight field line equilibrium conversion, see #1880 desc.compat.rescale will now return ScaledProfile instances for most of its profiles, to fix a bug where improper scaling could occur for certain profile types. Now always use sym=False in the default grid for plot_fsa to ensure correct averages Fixes bug that could lead extra compilation of jit-compiled functions that include field_line_integrate. Fixes inaccurate normalizations scales that could be computed for certain equilibria which had m=1 n=0 R and m=-1 n=0 Z components much smaller than their actual average minor radius, see #1954 Fix bug in PlasmaCoilSetMinDistance that occured using a FourierRZToroidalSurface object without passing in an the evaluation grid, see #2013 Equilibrium profile assignments are now guaranteed to be consistent with the equilibrium resolution—automatically increasing lower-resolution profiles to match the equilibrium (while keeping higher-resolution profiles untouched)—meaning users who relied on lower-resolution profiles to implicitly restrict optimization must now explicitly use the FixParameters constraint. Allow desc.vmec.VMECIO.load to load wout files that lack lrfp__logical__, like those outputted by VMEC++. This change assumes that those output files don't have poloidal flux label. Fixes a bug that had prevented passing the legend kwarg to desc.plotting.plot_surfaces. Backend When using any of the "proximal-" optimization methods, the equilbrium is now always solved before beginning optimization to the specified tolerance (as determined, for example, by options={"solve_options":{"ftol"...}} passed to the desc.optimize.Optimizer.optimize call). This ensures the assumptions of the proximal projection method are enforced starting from the first step of the optimization. desc.continuation.solve_continuation_automatic now falls back to performing shape perturbations first and then pressure if the default pressure-then-shaping path fails, increasing robustness in arriving at the final equilibrium. To go directly to the path of applying shaping then pressure, pass the flag shaping_first=True. desc.equilibrium.Equilibrium.set_initial_guess now sets lambda to zero for most use cases, and the docstring has been updated to be more explicit on what is done in each case. Minimum JAX version bumped up to 0.5.0 Performance Improvements Coordinate mapping uses partial summation, see #1826 Non-uniform FFTs (NUFFTS) are now used by default for computing bounce integrals, see #1834. NUFFT functionality is added through jax-finufft package. If the GPU installation fails, users can fall back to the older (much slower) implementation by setting nufft_eps=0 in the computation of the related quantities. New Contributors @singh-jaydeep made their first contribution in https://github.com/PlasmaControl/DESC/pull/1905 Full Changelog: https://github.com/PlasmaControl/DESC/compare/v0.15.0...v0.16.0

To better understand the dependence of the magnetic field structure in the plasma edge on the plasma boundary shape, in the context of X-point and island divertor designs, we define and develop a class of stellarators called umbilic stellarators. These equilibria are characterised by a single continuous high-curvature edge on the plasma boundary that goes around multiple times toroidally before meeting itself. We develop a technique that allows us to simultaneously optimise the plasma boundary along with a curve lying on the boundary on which we impose a high curvature while imposing omnigenity – a property of the magnetic field that ensures trapped particle confinement throughout the plasma volume. We find that umbilic stellarators naturally tend to favour piecewise omnigenity instead of omnigenity with a specific helicity. After generating omnigenous umbilic stellarators, we design coil sets for some of them and explore the fieldline structure in the edge and its sensitivity to small fluctuations in the plasma. Finally, using single-stage optimisation, we simultaneously modify the plasma and coil shape and propose an experiment to modify an existing tokamak to a finite- $\beta$ stellarator using this technique and explore a potentially simpler way to convert a limited tokamak into a diverted stellarator.

In the pursuit of the goal of commercial fusion as an abundant and safe source of energy, the stellarator is a leading concept with compelling attractiveness and demonstrated performance. In this white paper we, as a community of US researchers from Universities, National Laboratories, and Private Industry, involved in studying the stellarator concept, lay out the programmatic and technical motivation for a new mid-size stellarator research facility in the US. This contribution is complementary to several other whitepapers authored by members of our community which address different mid-scale stellarator aspects. A community based technical facility proposal has been prepared by F. Parra, et al: Flexible Stellarator Physics Facility. Two private stellarator companies have submitted proposals supporting the development of a mid-scale stellarator: Thea Energy (C.P.S. Swanson, et al.), Type One Energy (W. Guttenfelder, et al.).