European effort towards edge fluid modelling tools for self-consistent reactor-relevant simulations: main outcomes and remaining gaps (Turbulence codes)

• Patrick Tamain (CEA)

Room: HOW room Predictive modelling of plasma edge and divertor physics is a key challenge for ITER and future fusion power plants, where narrow operational margins require self-consistent treatment of turbulence, neutrals, magnetic geometry, and plasma-wall interactions. The TSVV3 project was established to advance European edge fluid turbulence modelling tools towards reactor-relevant predictive capability. We report here on its main outcomes and discuss identified remaining gaps. The TSVV3 project coordinated developments across the main European edge turbulence codes (GBS, GRILLIX, FELTOR and SOLEDGE3X), combining advances in numerical methods, physical models and high-performance computing. The project addressed four major challenges: 1- realistic magnetic and wall geometries, 2- self-consistent neutral dynamics in detached divertor conditions, 3- turbulence modelling in high-performance plasma regimes, and 4- scalability towards reactor-sized simulations. Code developments were accompanied by validation against experiments on several European tokamaks and by applications to both tokamak and stellarator configurations. The project enabled arbitrary axisymmetric magnetic geometries in all contributing European edge turbulence codes and extended simulations to three-dimensional configurations, including stellarators. Self-consistent neutral models based on fluid, kinetic and coupled Monte Carlo approaches were implemented, allowing first turbulence simulations in detached divertor regimes and demonstrating significant turbulence-induced scrape-off layer broadening in line with experimental results. Electromagnetic turbulence and improved collisional closures extended the applicability of fluid models towards high-β and low-collisionality plasmas, eventually enabling first self-consistent edge turbulence simulations in H-mode conditions. In parallel, substantial advances in numerical algorithms and GPU acceleration improved computational performance. Validation activities across multiple experimental devices have demonstrated the growing flexibility and predictive capability of these models. While TSVV3 has significantly expanded the physics fidelity and applicability of European edge fluid turbulence codes, significant challenges remain on the path to self-consistent reactor-scale simulations. These include improved wall geometry treatment, scalable kinetic neutral modelling, impurity physics in three-dimensional turbulence simulations, and numerical strategies capable of addressing the spatial and temporal scales of fusion power plants. We will browse through some of these issues and highlight possible solutions currently explored in the TSVV-B project which took over TSVV3.

Numerical challenges in SOLPS-ITER modelling for STEP divertor design (Mean-field codes)

• Ryoko Osawa (UKAEA)

Room: HOW room We have been using the SOLPS-ITER code to design and optimise the divertor for STEP (Spherical Tokamak for Energy Production). While this modelling has yielded important physical insights — such as the impacts of fuelling locations and dome shape — it has also highlighted significant numerical challenges. These issues introduce additional uncertainties on top of existing variances driven by assumed input parameters and reduced physics models. This presentation focuses on the two most significant numerical challenges encountered: unexpected up-down asymmetries under nominally symmetric conditions, and particle balance errors. We will detail the possible causes of these numerical issues and discuss actionable strategies to achieve more robust, reliable SOL/Divertor modelling within the timescales demanded by STEP.

Global fluid turbulence simulations in a general frame (Turbulence codes)

• Brendan Shanahan (MIT Plasma Science and Fusion Center)

Room: HOW room Toroidal magnetic confinement fusion experiments require a rotational transform to achieve adequate confinement. This rotational transform can be achieved via any combination of internal plasma currents, torsion of the magnetic axis, and the toroidal rotation of the plasma cross section. Internal currents can lead to current-driven instabilities, and therefore the first stellarators employed magnetic axis torsion in a figure-8 form to provide their rotational transform [1]. These early configurations were unstable, however, and the figure-8 stellarator was abandoned in favor of the toroidally-shaped tokamaks and classical stellarators. As a result, nearly every simulation and analysis framework for magnetic fusion experiments has been written assuming a toroidal geometry. Recently, however, advancements in near-axis theory [2] and equilibrium reconstruction codes [3] have renewed the interest in configurations unrealizable by traditional simulation and analysis frameworks [4]. Here we demonstrate the capability of the BSTING project [5, 6] to simulate global fluid transport and turbulence a generalized frame. We utilize a hot-ion turbulence model within the Hermes-3 family of multifluid models [7], which provides an inherent flexibility of model fidelity. As a proof-of-principle, transport characteristics of an optimized quasi-isodynamic figure-8 stellarator [4] are presented. It is determined that the figure-8 indicates strong ballooning transport in high-curvature regions, with radial heat and particle flux bands extending into zero-curvature regions. A comparison of the turbulent flux with known turbulence scalings, as well as how the flux compares to that in conventional stellarator geometries is provided.

Dynamic system identification simulations of MAST-U using SOLPS-ITER (Kinetic neutrals)

• Stijn Kobussen (DIFFER, Dutch Institute for Fundamental Energy Research)

Room: HOW room The safe and controlled exhaust of heat and particles from magnetically confined fusion plasmas requires having dynamic information about the system. On present-day reactors, this information is obtained using system identification experiments, which involve perturbing the gas injection rate and observing the Scrape-Off-Layer (SOL) response in the frequency domain. However, this strategy is cumbersome for future reactors, as accidental disruptions or reattachment of the SOL plasma could lead to significant damage to the device. Therefore, numerical dynamic models of the exhaust plasma, validated against experiments, are needed. To this end, we present a multi-sine gas injection perturbation approach that uses time-dependent SOLPS-ITER simulations, with a grid extending to the vessel wall, to quantify the plasma response. Comparison to an existing set of experimental data on MAST-U is done to validate this approach. The simulation grid was generated using the Grid Optimization and Adaptation Toolbox (GOAT), and both the Advanced Fluid Neutral (AFN) model and the kinetic neutral model Eirene available in SOLPS-ITER are evaluated. In steady state, it is possible to simulate divertor states ranging from attached to deeply detached using both fluid and kinetic neutrals. The comparison in the frequency domain shows that the AFN SOLPS-ITER simulations predict faster response times to divertor fuelling perturbations than the experiments at low frequencies, likely in part due to the fluid approximation used for the neutrals. Simulations employing time-dependent Eirene for kinetic treatment of neutrals show good agreement with the experimental data throughout the measured frequency range of around 9 to 30 Hz. Simulations with Eirene using the Quasi-Steady-State (QSS) approximation result in different dynamic behaviour due to the non-linear coupling between the plasma solver B2.5 and Eirene. Proper time-dependent kinetic treatments of the neutrals is therefore essential for accurately simulating plasma evolution in response to changes in gas injection.

Hermes-3 simulation of plasmas confined by a levitated dipole magnet (Turbulence codes)

• Yichen Fu (Columbia University)

Room: HOW room Laboratory dipole confinement was proposed by Hasegawa in 1987 as a possible route to fusion-relevant plasmas and has been explored in levitated-dipole experiments such as LDX, RT-1, and the recent Junior device. A dipole magnetic field is characterized by its strong magnetic-flux expansion, i.e., the magnetic field strength varies strongly along a field line, which poses challenges for numerical simulations. In this work, we use Hermes-3 to study three-dimensional plasma dynamics near the edge of a dipole-confined plasma. Fluctuations in the closed-field-line region in dipole-confined plasmas are often believed to be flute-like with $k_\parallel\approx 0$. We find that this picture can break down near a marginal stationary profile. In slightly supercritical profiles, three-dimensional simulations reveal that non-flute fluctuations with finite parallel structure naturally emerge due to the strong variation of the local magnetic geometry along the field line. Furthermore, these fluctuations produce a differential radial transport: along the same field line, particles are pinched inward near the outer midplane while being driven outward near the inner midplane. This poloidal-angle-dependent transport is not captured by purely flute-like, flux-tube-averaged descriptions. This result is further supported by eigenmode solutions from a simplified drift-reduced fluid model, whose quasilinear flux reproduces the differential radial transport. Future work will extend the simulation domain to include the X-point and open-field-line regions in a dipole field.

Implementation and Analysis of Transport Simulations with External 3D Fields in Hermes-3 (Mean-field codes)

• Sidney Williams (University of California San Diego)

Room: HOW room Despite tokamaks being primarily treated as axisymmetric beasts, many aspects of actually running such a machine are inherently nonaxisymmetric. This includes error fields, magnetic ripple, and even Edge Localized Mode-suppressing (ELM-suppressing) Resonant Magnetic Perturbations (RMPs). Each of these effects are extremely important to consider while the fusion field is moving towards next generation machines. Specifically, recent SPARC studies have addressed error field mitigation [1], and magnetic ripple-driven alpha particle losses [2], additionally many ITER studies have explored the implementation of RMPs to suppress machine-killing ELMs [3]. Upon adding 3D magnetic fields, transport is strongly affected. Loss of axisymmetry breaks integrability of field line trajectories, resulting in a chaotic equilibrium magnetic field. This chaotic field changes the structure of the differential operators within the governing transport equations, and provides additional transport in and of itself [4]. In this work, these effects are studied through a transport simulation of an RMP-applied DIII-D shot, although the methodology could nominally be applied to any other machine, and any other mode of externally-applied asymmetry. The implementation methodology is walked through, highlighting difficulties, and possible pitfalls as well as interesting physics effects which become necessary when considering a chaotic equilibrium field. Finally, the results of the study, such as the effect of 3D fields on anomalous transport coefficients, and divertor heat flux patterns are presented. [1] S. Munaretto et.al., Nuclear Fusion, (2025). [2] S. D. Scott et.al., Journal of Plasma Physics, (2020). [3] L. Zhou et.al., Plasma Phys. Control., (2016). [4] A. B. Rechester and M. N. Rosenbluth, Physical Review Letters, (1978).

Progress of the drift and turbulence capabilities of Hermes-3 for non-axisymmetric geometries (Mean-field codes)

• Tobias Tork (IPP Greifswald)

Room: HOW room Particle drifts play a crucial role in understanding of the edge physics of current and future fusion devices. While tokamak simulations routinely include drift effects, they have been largely neglected in stellarator simulations until now. Advancements in BOUT++ enable the efficient simulation of stellarator geometries. Hermes-3 is a plasma fluid model for the edge of fusion devices and covers a wide range of levels of fidelity. In this work, we report on the progress of the Hermes3 implementation for non-axisymmetric geometries using the flux-coordinates-independent (FCI) approach. We elucidate the implementation for steady-state transport simulations with drifts and show a first example simulation of W7-X. Additionally, we demonstrate the turbulence capabilities by simulating a W7-AS-like stellarator geometry. In addition to recent progress, we show some of the remaining challenges such as neutral physics or the creation of closed divertor grids.

Immersed Boundary Extension to BOUT++/Hermes-3 (Turbulence codes)

• Stefan Tirkas (Lawrence Livermore National Laboratory)

Room: HOW room Realizing commercially viable fusion power plants (FPPs) demands whole-device modeling efforts which integrate high‐fidelity plasma simulations with engineering systems and design capabilities. Among the most pressing challenges are the accurate prediction of heat and particle exhaust in the divertor region and the self‐consistent treatment of plasma–neutral and plasma–wall interactions in the edge and scrape‐off layer. In this work, the BOUT++/Hermes-3 [1] fluid framework is extended with a Cartesian flux-coordinate-independent (FCI) discretization [2], allowing for full-device transport and turbulence studies from the core to the first wall. To fully characterize the plasma, the device wall is represented as an immersed boundary [3], from which boundary conditions are imposed on the perpendicular anomalous diffusion and ExB advection terms. For terms employing a finite volume method, plasma cells intersected by the boundary are further treated using a cut-cell discretization [4]. The parallel discretization and associated boundary conditions are presently inherited from the quasi-FCI implementation BSTING [5]. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LLNL-ABS-2020022 [1] B. Dudson, M. Kryjak, H. Muhammed, P. Hill, and J. Omotani Comp. Phys. Comm. 296 108991 (2024) [2] F. Hariri and M. Ottaviani Comp. Phys. Comm. 184 2419-2429 (2013) [3] R. Ghias, R. Mittal, and H. Dong Journal of Computational Physics 225 528–553 (2007) [4] H. Johansen and P. Colella Journal of Computational Physics 147, 60-85 (1998) [5] B. Shanahan, B. Dudson, and P. Hill Plasma Phys. Control. Fusion 61 025007 (2019)

On the challenge of kinetic neutrals modelling in reactor relevant conditions: characterization of the issues and possible solutions (Kinetic neutrals)

• Yannick Marandet (Aix-Marseille University / CNRS)

Room: HOW room The particle, momentum and energy exchanges between the plasma and the neutral gas from recycling are key aspects of divertor physics. Atoms and molecules are often not in a collisional regime and a kinetic description is required to get transport right. In detached regimes, where the temperature in the divertor drops below a few eVs, the complexity of the reaction channels for neutrals increases. The Monte Carlo approach implemented in EIRENE allows to solve the kinetic problem on a 3D grid, and can accommodate both the geometrical complexity and the relevant species and reactions channels. However, in large machines such as ITER, in which the mean free path of neutrals can be very short compared to the size of the divertor, regions where neutrals become collisional can appear. Such regions make the Monte Carlo approach computationally much more costly, since atoms or molecules may undergo tens of thousands of collisions before being ionized or pumped. It also makes fluid approaches accurate at least in these regions. This presentation will focus on EIRENE simulations in ITER conditions, showing that the computational cost increase is related to physics and that sacrificing the physics is not an option if the goal is to assess e.g. peak heat fluxes in the divertor. Advanced fluid models, fully consistent with the underlying kinetic model used in EIRENE, have been developed for atoms and shown promising results. However, they are generally not valid everywhere in the simulation domain and hybrid kinetic-fluid models are shown to provide a way to combine accuracy and efficiency.

The KINetic Deterministic NEutral Solver for turbulence simulations in the boundary of magnetic confinement devices (Kinetic neutrals)

• Davide Mancini (EPFL)

Room: HOW room We present KINDNES (KINetic Deterministic NEutral Solver), a tool for solving the kinetic Boltzmann equation applied to neutral particle dynamics in the boundary plasma of tokamaks and stellarators. Originally developed as the neutral model within the GBS turbulence code, KINDNES is now developed as a standalone library. Neutral particles play a central role in edge plasma physics, governing key processes such as recycling, momentum and energy exchange, and the onset of plasma detachment. While existing neutral models typically rely either on fluid approximations or on Monte Carlo kinetic methods, KINDNES provides a fully kinetic, deterministic framework, free of statistical noise. The model discretizes the Boltzmann equation for each neutral species and integrates it along its characteristics, reducing the problem to the inversion of a discrete linear system. The fundamental physical mechanisms and numerical implementation of KINDNES are first presented, followed by its first benchmark against the Monte Carlo code EIRENE. Using a deuterium plasma background derived from a SOLPS simulation, the two codes are compared in an attached divertor configuration, with EIRENE configured to match the physical assumptions currently available in KINDNES. Within this common physics framework, the two codes show good agreement in atomic deuterium density, temperature and velocities, representing an encouraging validation of KINDNES as a deterministic alternative for neutral particle modeling. Building on this validation, the most recent developments of KINDNES are then presented: the extension to molecular deuterium dynamics, enabling the simulation of detached divertor conditions; the introduction of flexible wall geometry support, broadening the range of device configurations that can be modeled; the development of a full 3D solver, that includes the proper tratment of toroidal velocities; the relaxation of the adiabatic approximation, rendering the solver aware of the plasma evolution in time and the adoption of a Hierarchical Matrices method to substantially accelerate the solution of the linear system in terms of both memory and computational cost.

Using SPARTA as a standalone DSMC test bed for divertor neutral transport and exhaust boundary conditions (Kinetic neutrals)

• Matteo Moscheni (Gauss Fusion GmbH, Parkring 29, 85748 Garching bei München, Germany)

Room: HOW room Kinetic neutral transport remains one of the least reducible components of fusion exhaust modelling. While coupled plasma-edge tools such as SOLPS-ITER/EIRENE remain the standard route for integrated divertor simulations, standalone kinetic-neutral calculations provide a complementary way to isolate rarefied-gas effects, test boundary-condition assumptions, and interpret neutral transport patterns without the full complexity of plasma feedback. In this contribution, we discuss the use of SPARTA, a direct simulation Monte Carlo code for rarefied-gas dynamics, as a standalone tool for divertor-relevant neutral transport studies. The emphasis is not on replacing established plasma-edge workflows, but on identifying where a DSMC treatment can provide useful diagnostic and modelling insight. We first outline the practical advantages and limitations of this approach relative to EIRENE-like neutral models, including the use of binary neutral--neutral collisions rather than BGK-type closures, the treatment of neutral--wall interaction models such as CLL and TRIM-informed reflection, and the construction of plasma-mimicking boundary conditions from prescribed neutral sources and fluxes. We then present a set of simplified private-flux-region and toroidal-plane test cases motivated by divertor exhaust applications. These cases show that neutrals need not behave as an isotropic diffusive background: ordered sources, wall interaction, geometry, and neutral--neutral collisions can generate structured flows, preferential capture, backflow, and non-trivial momentum transport. Such behaviour has direct implications for reduced exhaust models, including the common treatment of pump surfaces through fixed albedo or absorption coefficients. In particular, DSMC calculations allow pump capture, return flux, and effective exhaust conductance to be estimated from the local rarefied flow rather than prescribed as constants. The talk will focus on lessons for kinetic-neutral and exhaust-code developers: which assumptions are exposed by standalone DSMC, which quantities can be passed back to integrated edge codes, and where systematic comparison against EIRENE-like models is still needed.

A fluid model for the deuterium molecules in the plasma edge based on the AMJUEL database and benchmark with EIRENE (Mean-field codes)

• Wim Van Uytven (KU Leuven (EU))

Room: HOW room In plasma boundary modelling, fluid neutrals offer several practical disadvantages compared to their more accurate kinetic Monte Carlo counterparts. They are typically much faster, and their deterministic nature is a major benefit when doing optimization or when coupling to plasma turbulence codes. Besides missing kinetic effects, fluid neutral models are typically also lagging behind with regard to the physics processes they include. The vast majority of present-day fluid neutral models consider only atoms. The first commonly referenced fluid molecular models in the context of plasma edge models were developed in UEDGE. The UEDGE work focused primarily on the details of the molecular physics (Collisional Radiative Models) and their impact on detachment in DIII-D, while the benchmark with EIRENE was more restricted. In this contribution, we aim for a complementary approach, considering only the default molecular physics model from EIRENE, but performing for the first time a thorough benchmark between the fluid model and EIRENE for all macroscopic parameters (density, velocity and temperature) and all exchange terms with the plasma and atom populations (ion/atom particle sources, ion/atom parallel momentum sources, electron energy sources and ion/atom energy sources), on a plasma background representative of the detached regime. We show that the particle velocity distribution of the kinetic molecules approaches a Maxwellian in a large part of the relevant domain, due to the elastic collisions with the detached plasma background. We also show that the bulk of the molecular physics can be captured reasonably well with the presented fluid molecular model, showcasing it to be an excellent starting point for further research in this direction.

Propagator-based Multi-level Monte-Carlo for Kinetic Neutral Species (Kinetic neutrals)

• Gregory Parker (UC Berkeley)

Room: HOW room In recent work, we introduced a novel multi-level Monte-Carlo scheme for kinetic simulations of neutral particles in plasmas near the isotropic cross-section limit. Of particular interest is the case of edge plasmas in detatched or near-detached regimes. The multi-level scheme in this case is based on a propagator formalism enumerating discrete charge-exchange collision events. This scheme reproduces the results of known algorithms in the large particle number and fine mesh limit. More importantly, the new formalism shows promise for differentiable and deterministic particle-based models. At the end, we will discuss preliminary work implementing this method for fully-implicit coupled plasma-kinetic neutral models. This talk/poster is based on joint work with Maxim Umansky and Ben Dudson.

NESO-Particles: A Performance Portable Particle Framework for Fusion Plasma Simulation (Research Software Engineering)

• Will Saunders (UKAEA)

Room: HOW room Simulation of fusion plasma often requires representations of non-Maxwellian distributions. A common method is to represent the distribution function as a collection of individual particles where the state of each particle evolves over time. In this talk we describe NESO-Particles (NP) which is a performance portable framework for describing particle based algorithms. NP provides users with abstractions for describing particle data and looping operations for the fusion plasma use case. These abstractions form a separation of concerns between the description of particle based algorithms and the execution of these algorithms on modern HPC hardware. Our implementation is realised as a MPI+SYCL library that executes on both CPU and GPU hardware.

The kinetic neutral model in MHD code JOREK (Kinetic neutrals)

• Daniël Maris (DIFFER)

Room: HOW room A few years ago, the nonlinear extended MHD code JOREK [2] has been extended with a dedicated kinetic neutral model [3], opening up the possibility to transiently model the interplay between edge-plasma and plasma-wall interactions and the plasma processes for which JOREK is well known such as ELMs, disruptions and RE beams, where previously these worlds could not be self-consistently combined in any code. Since then, the kinetic neutral model has been further developed and used as exhaust code by a growing part of the JOREK community and has attracted increasingly more attention from the exhaust physics field. The MHD model in JOREK uses a flux surface aligned poloidal finite element grid up to the true first wall with the option to model variation in the toroidal direction using Fourier expansion. The kinetic neutral model is built on JOREK’s hybrid OMP-MPI particle-in-cell framework [4] and follows an array of variable weight particles through time, accumulating sources and sinks that are fed into the plasma equations. Both the MHD side and the kinetic side are inherently suited and optimized for time dependent modelling, making the model interesting for investigating dynamic physics and for control. Included in the neutral model physics are ionisation, charge exchange, wall – and volume recombination, line radiation and puffing, while molecular interactions are currently still missing. As part of recent developments, various kinetic applications in JOREK (neutrals, impurities, runaway electrons and fast particles) have been combined into a single flexible kinetic framework, with object oriented control of the simulation setup, such that the included particle species, the particle pushers used for each species and the kinetic-MHD coupling scheme can be selected. Additionally, the neutral model is now compatible with a two temperature plasma description; direct binary neutral self-collisions have been added; and capability to resolve simulation actions (such as neutral collisions) on timescales in between the MHD timestep and the kinetic particle timestep have been implemented. Furthermore, there is work underway for GPU acceleration and making the kinetic neutral model compatible with stellarator geometry. This talk will discuss JOREK’s kinetic neutral model, how it is implemented and what is and is not currently included in the physics, along with details on the recent upgrades, encountered issues, and a selection of application and validation simulations of small ELM burn-through, XPR formation and transition into a MARFE, QCE access, and dynamic detachment studies. [1] Hoelzl, Nucl. Fusion 66 (2026) 116006. [2] Hoelzl, Nucl. Fusion 61 (2021) 065001. [3] Korving, Phys. Plasmas 30 (2023) 042509. [4] van Vugt, PhD Thesis, TU Eindhoven, 978-90-386-4811-8 (2019).

VANTAGE-Reactions – performance-portable neutral reactive physics for exhaust modelling (Kinetic neutrals)

• Stefan Mijin (UKAEA)

Room: HOW room UKAEA (United Kingdom Atomic Energy Authority), Culham Campus, Abingdon, Oxfordshire, OX14 3DB, UK. The accurate simulation of neutral particles in the tokamak exhaust has been a topic of research across many decades, with the primary numerical approach being Markov Chain Monte Carlo methods, widely used in linear transport problems[1,2]. Historically, neutral-neutral collisions have been treated using the BGK method in order to fit in with the linear transport constraints[3], with reduced models recently gaining traction[4]. Recent work at UKAEA has been focused on performance-portable particle libraries built on top of NESO-Particles[5], with VANTAGE-Reactions[6] providing abstractions and general implementations for plasma-neutral, particle-surface, and neutral-neutral interactions. It is actively being integrated into finite volume (Hermes-3[7]) and finite element (PENKNIFE[8]) codes. Following a less common path to kinetic neutral modelling for tokamaks[9,10,11], we build on methods from rarefied gas dynamics, specifically the Direct Simulation Monte Carlo method[12], with a focus on performance-portability and flexibility in both steady state and time-dependent scenarios. In this talk, we will present the general weighted particle approach used by VANTAGE-Reactions, benchmarking results so far, as well as early progress in development of binary collision support through an event-splitting stochastic weighted particle method[13,14]. This work has been part-funded by the EPSRC Fusion Grant 2022/27 [grant number EP/W006839/1]. References: [1] D. Reiter, Journal of Nuclear Materials, 196–198 80–89 (1992) [2] D.P. Stotler et al. “DEGAS 2 neutral transport modeling of high density, low temperature plasmas.” (1997). [3] V. Kotov, et al., Plasma Physics and Controlled Fusion, 50 10 (2008) [4] D. V. Borodin et al., Nuclear Fusion, 62 8 (2022) [5] https://github.com/ExCALIBUR-NEPTUNE/NESO-Particles [6] https://github.com/UKAEA-Edge-Code/VANTAGE-Reactions [7] B. Dudson et al. Computer Physics Communications 296 2024 [8] https://github.com/ExCALIBUR-NEPTUNE/PENKNIFE [9] S. Varoutis et al., Fusion Engineering and Design, 121 13–21 (2017) [10] S. Q. Korving et al, Physics of Plasmas, 30 4 (2023) [11] K. Kvist et al., Physics of Plasmas, 31 3 (2024) [12] G. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Oxford University Press (1994) [13] S. Rjasanow et al., Journal of Computational Physics, 124. 2 243–253 (1996) [14] G. Oblapenko et al. Journal of Computational Physics 466 (2022)

Hermes-3 performance and neutral model update (Mean-field codes)

• Mike Kryjak

Room: HOW room Hermes-3 is an open-source fluid edge plasma code written using the BOUT++ framework. Building on the features of SD1D, Hermes-2 and STORM, it enables self-consistent simulation of scrape-off layer plasmas in 1D, 2D and 3D for both steady-state and unsteady/turbulent regimes. In this work, we present an update on several activities contributing towards improved neutral model capability in the code. Hermes-3 is a fully implicit code. In 2D mean-field, it relies on the SNES solver in PETSc with backward Euler time integration and an LU-type preconditioner. The implementation of separate flux limiters for cross-field neutral advection, conduction and viscosity led to severe performance degradation. Upon investigation, the limiter implementation was found to lead to several numerical problems due to instances of non-differentiability and insufficient regularisation. Resolving these has led to a marked speedup and improvement in robustness across a range of levels of limiter saturation. Following this, implementing a fully differentiable slope limiter has led to reduced solver failures, steadier timesteps and improved performance. A further significant speed enhancement was found by selectively lagging terms in the neutral flux limiter formulation, revealing their large impact on the stiffness of the system. We also present an update on coupling Hermes-3 to VANTAGE (Versatile, Accurate Neutral Transport and Analysis for GPU Execution), a new performance-portable kinetic neutral model. Mass conservation is demonstrated in a simple slab simulation with an ionization-recombination balance between the fluid plasma and kinetic neutrals. This work has been funded by the Fusion Futures Programme. As announced by the UK Government in October 2023, Fusion Futures aims to provide holistic support for the development of the fusion sector. Performed in part under the auspices of the U.S. DOE by LLNL under contract DE-AC52-07NA27344.

SOLEDGE-HDG: a novel unstructured approach for transport simulations in tokamak-relevant conditions (Mean-field codes)

• Frederic Schwander (Aix-Marseille University / M2P2)
• Eric Serre (Aix-Marseille University / M2P2)

Room: HOW room Fluid transport simulations based on mean-field quantities, with the averaged out turbulence fluctuations being modelled, remain a standard in the international community when considering engineering studies of tokamak relevant configurations (for a review of state-of-the-art transport codes, see Schwander et al. 2024). In the perspective of enhancing computational efficiency and expanding codes capabilities, we started developing SOLEDGE-HDG, a high-order, finite-element code within the SOLEDGE code suite approximately ten years ago [Giorgiani et al. 2018 ; Capasso et al. 2025]. The hybridization of the DG method enables a special choice of the numerical traces (or numerical fluxes) for spatial discretization, which makes the HDG stand out thanks to its stability features, reduced number of degrees of freedom, and super-convergence properties (the error converges faster than the degree of the polynomials used to represent approximate the solution) [Giorgiani et al. 2018]. The innovation of this HDG method for tokamak simulations lies in its mesh flexibility which eliminates the need for alignment with magnetic field lines or flux surfaces. This enables accurate discretization of complex plasma-facing components and singularities (e.g., X-points) as well as the simulation of non-steady phases (e.g., start-up) without costly remeshing, addressing a gap left by current codes. Relieving the constraint of aligment with magnetic field comes at the price of additional spurious numerical diffusion in the perpendicular direction that can be bounded as long as high-order interpolations are used [Giorgiani et al. 2020]. The HDG scheme also promotes implicit time integration that relies on rapid Newton-Raphson convergence, thus allowing large time steps and fast computational access to steady states, leading to superior performance of SOLEDGE-HDG in achieving 2D transport equilibria when compared to semi-implicit codes. In this presentation, we will detail the algorithm, highlighting both its attractive properties and its current limitations. Over recent years, the code has been successfully applied to simulate tokamak-relevant configurations [Scotto et al. 2022, Kudashev et al. 2026]. Selected results will be presented to illustrate the code’s new capabilities and its strong potential for addressing transport and heat exhaust challenges in existing devices, as well as in next-generation machines such as ITER [Scotto et al. 2024] or SPARC.

Advancing BIT1 Towards Exascale: Multi-GPU Hybrid Particle-in-Cell Monte Carlo Simulations at Scale (Research Software Engineering)

• Jeremy J. Williams (KTH Royal Institute of Technology)

Room: HOW room Particle-in-Cell (PIC) Monte Carlo (MC) simulations of the plasma edge play an important role in magnetic confinement fusion research, providing insights into plasma-surface interactions and plasma, impurity and neutral particle transport, relevant to both present and future fusion devices. As plasma simulations continue to increase in size and complexity, efficiently exploiting modern heterogeneous supercomputing systems presents significant challenges, including data movement, load imbalance, scalability, and resilience. This talk presents recent developments in advancing the Berkeley Innsbruck Tbilisi 1D3V (BIT1) PIC MC code **[1,2]** towards exascale computing through a portable hybrid MPI+OpenMP implementation capable of leveraging large-scale multi-GPU systems across both NVIDIA and AMD architectures. The work focuses on practical aspects of code development and optimisation, including GPU offloading strategies, overlapping communication and computation, memory management approaches, particle load balancing, and scalable checkpoint/restart capabilities. The integration of openPMD and ADIOS2 for standardised high-performance I/O is also discussed, enabling parallel I/O **[3]**, enhanced diagnostics **[4]**, in-memory data streaming **[5]**, and in-situ analysis and visualisation workflows **[4,5]**. Performance and scalability results from leading pre-exascale and exascale systems, including Frontier, LUMI-G and MareNostrum 5 ACC are presented, together with lessons learned from profiling **[6]**, porting **[7,8]**, and enabling portability **[9]**. The talk concludes with ongoing efforts to improve resilience and workflow efficiency **[10]**, together with future research extending hybrid BIT1 to Intel GPUs towards exascale, targeting Aurora and Europe's first exascale system, JUPITER Booster, to assess portability and performance across NVIDIA, AMD, and Intel GPUs and identify any remaining architecture-specific limitations at scale. ### References **[1]** D. Tskhakaya, et al., “Optimization of PIC codes by improved memory management,” Journal of Computational Physics, vol. 225, no. 1, pp. 829–839, 2007. [doi:10.1016/j.jcp.2007.01.002](https://doi.org/10.1016/j.jcp.2007.01.002) **[2]** D. Tskhakaya, et al., “PIC/MC code BIT1 for plasma simulations on hpc,” in 2010 18th Euromicro, pp. 476–481, IEEE, 2010. [doi:10.1109/PDP.2010.47](https://doi.org/10.1109/PDP.2010.47) **[3]** J.J. Williams, et al., "Enabling high-throughput parallel I/O in particle-in-cell Monte Carlo simulations with OpenPMD and darshan I/O monitoring." in 2024 IEEE International Conference on Cluster Computing Workshops (CLUSTER Workshops), IEEE, 2024. [doi:10.1109/CLUSTERWorkshops61563.2024.00022](https://doi.org/10.1109/CLUSTERWorkshops61563.2024.00022) **[4]** J.J. Williams, et al., “Understanding the Impact of OpenPMD on BIT1, a Particle-in-Cell Monte Carlo Code, Through Instrumentation, Monitoring, and In-Situ Analysis,” in European Conference on Parallel Processing, pp. 214–226, Springer Nature Switzerland, 2024. [doi:10.1007/978-3-031-90200-0_18](https://doi.org/10.1007/978-3-031-90200-0_18) **[5]** J.J. Williams, et al., “Integrating High Performance In-Memory Data Streaming and In-Situ Visualization in Hybrid MPI+ OpenMP PIC MC Simulations Towards Exascale,” The International Journal of High Performance Computing Applications, p. 10943420251409229, 2025. [doi:10.1177/10943420251409229](https://doi.org/10.1177/10943420251409229) **[6]** J.J. Williams, et al., “Leveraging HPC Profiling and Tracing Tools to Understand the Performance of Particle-in-Cell Monte Carlo Simulations,” in European Conference on Parallel Processing, pp. 123–134, Springer Nature Switzerland, 2023. [doi:10.1007/978-3-031-50684-0_10](https://doi.org/10.1007/978-3-031-50684-0_10) **[7]** J.J. Williams, et al., “Optimizing BIT1, a Particle-in-Cell Monte Carlo Code, with OpenMP/OpenACC and GPU Acceleration,” in International Conference on Computational Science, pp. 316–330, Springer Nature Switzerland, 2024. [doi:10.1007/978-3-031-63749-0_22](https://doi.org/10.1007/978-3-031-63749-0_22) **[8]** J.J. Williams, et al., “Accelerating Particle-in-Cell Monte Carlo Simulations with MPI, OpenMP/OpenACC and Asynchronous Multi-GPU Programming,” Journal of Computational Science, vol. 88, p. 102590, 2025. [doi:10.1016/j.jocs.2025.102590](https://doi.org/10.1016/j.jocs.2025.102590) **[9]** J.J. Williams, et al., “Multi-GPU Hybrid Particle-in-Cell Monte Carlo Simulations for Exascale Computing Systems,” in International Conference on Computational Science, pp. 32–47, Springer Nature Switzerland, 2026. [doi:10.1007/978-3-032-29921-5_3](https://doi.org/10.1007/978-3-032-29921-5_3) **[10]** J.J. Williams, et al., “High-Performance Resilient Multi-GPU Hybrid Particle-in-Cell Monte Carlo Simulations at Scale,” in Euro-Par 2026: Parallel Processing Workshops: Euro-Par 2026 International Workshops ([BIGHPC](https://bighpc2026.di.unipi.it/)), Pisa, Italy, August 24–28, Accepted For Publication, 2026 [doi:10.48550/arXiv.2606.28534](https://arxiv.org/pdf/2606.28534)

Drifts and Turbulence in Hermes-3 (Turbulence codes)

• Ben Dudson (LLNL)

Room: HOW room We report on progress in modeling of plasma drifts and 3D turbulence using Hermes-3, in tokamak and stellarator configurations. During the past year improvements have been made to models for the diamagnetic and polarization currents, ion viscosity and energy conservation. Tokamak applications include simulation of turbulence in NSTX and NSTX-U configurations, with comparison to experimental data. Stellarator boundary modeling is challenging due to the complex magnetic field topology and wall geometry, in addition to the many interactions between plasma, neutrals and surfaces present in modeling of tokamak divertors. The Flux-Coordinate Independent (FCI) method has emerged as a promising approach to plasma discretization [1,2] that has sufficient flexibility to model real 3D devices. Recent advances in BOUT++ include relaxation methods to find steady-state plasma potentials[3], support-operator method (SOM) based FCI operators that conserve fluxes to machine precision on complex meshes, and nonlinear solvers and preconditioners to accelerate the solution of steady-state transport calculations. Progress has been made to port BOUT++ and Hermes-3 to GPUs. Field operators have been rewritten to use template expressions to enable kernel fusion. The design and tradeoffs and future plans will be discussed.

Overview of solver, Kinetics for Transport (Research Software Engineering)

• Sarah Newton (UKAEA)

Room: HOW room The dynamics of the edge of tokamak plasmas is understood to be sensitive to the electric field profile. Meanwhile modelling the long-wavelength component of the radial electric field is challenging for core physics, but essential to understand the rotation profile. The moment-kinetics approach to evaluating transport was introduced under the Excalibur project ([report TN-11][1]). This aims to evolve a suitably normalised distribution function on the transport timescale consistently with momentum conservation, across the separatrix and in the presence of the strong changes in parameter values typical of detaching plasmas. We present the work to date on the Kinetics for Transport (KfT) solver, which implements the moment-kinetics approach. The solver is being developed in Julia, in a finite element framework with a Chebyshev basis, and supports various configurations for physics studies. Benchmarks and physics studies in 1D1V and 1D2V have been completed, while 2D1V simulations with adiabatic electrons demonstrated the presence of an instability relying on finite pitch angle of the magnetic field. We have focused on developing a numerical approach that can simulate the sheath boundary conditions with kinetic ions and electrons with scalable performance, to enable 2D2V simulations. The implementation of the moment-kinetics structure and approach to electron preconditioning will be discussed. [1]: https://excalibur-neptune.github.io/Documents/TN-11_2DDriftKineticsHelicalFieldSeparatrix.html

Some aspects of the finite element method applied to edge plasma modelling (Turbulence codes)

• Ed Threlfall (UKAEA)

Room: HOW room The finite element method (FEM) concerns the numerical discretization of continuum partial differential equations. This talk briefly covers some applications of FEM to problems in plasma physics, including strongly orthotropic transport and plasma turbulence. We will emphasize the potential advantages of FEM in terms of numerical properties – stability, rapid convergence, and the suppression of numerical diffusivity (even without a fieldline-following computational geometry) – and geometrical fidelity – higher-order geometries, including faithful representations of curved surfaces. Finally, we will discuss how certain types of discrete representation, known as compatible - or mimetic- discretizations, preserve features of the continuum model.

Plasma Edge Numerics using Kinetic Neutrals Integrated into Finite Elements (Turbulence codes)

• James Edgeley

Room: HOW room UKAEA's edge code PENKNIFE is under development, incorporating the finite element framework Nektar++ and the reactions library VANTAGE. The code is designed to be used with modern hardware, exploiting the higher arithmetic intensity of finite element methods and the GPU parallelism of NESO-Particles for improved performance. A range of features have been implemented, encompassing a hierarchy of complexity from simple mean-field transport to 3D turbulence. Results from the latest simulations are presented, including comparisons with SOLPS.

KN1DPy: A new Python translation of the KN1D neutral transport code (Kinetic neutrals)

• Jamie Dunsmore (MIT)

Room: HOW room We present KN1DPy, a Python translation of KN1D, a spatially 1D kinetic neutral transport code that was originally written in IDL to model gas fuelling on the Alcator C-Mod tokamak. Given input profiles of electron density, $n_e$, and temperature, $T_e$, in the scrape-off-layer and on closed flux surfaces, KN1D solves the Boltzmann Equation to obtain the distribution functions for atomic and molecular neutral hydrogen. While its 1D geometry represents a step down in fidelity compared to 2D codes like EIRENE and DEGAS-2, this simplification also vastly reduces the computation time. Each KN1D run takes only a few seconds, and as a result KN1D can easily be included within integrated transport models to give an estimate for the edge ionisation source from gas fuelling. Furthermore, because KN1D returns the full distribution function for the neutrals, it also yields detailed information about the effective neutral temperature, the fraction of impinging neutrals reflected by the plasma, and the fraction of neutrals that have undergone charge-exchange. We show some examples from DIII-D and Alcator C-Mod, and include comparisons with results from 2D codes on these devices. We also demonstrate how KN1DPy can be used to provide a quick estimate for the neutral source in integrated transport models.

GUERNICA: A continuum kinetic code for neutral dynamics in fusion plasmas (Kinetic neutrals)

• Jack Gabriel (William & Mary)

Room: HOW room In this talk, GUERNICA, a new continuum kinetic neutral code, is presented for the simulation of atomic hydrogen transport in fusion edge plasmas, solving the Boltzmann equation with charge exchange, ionization and elastic neutral-neutral collisions. Neutral atoms strongly influence tokamak edge and divertor physics through fueling, momentum exchange and volumetric power loss [1]. Accurate neutral modeling is increasingly important for future devices, particularly for detachment which is central to managing power exhaust [2]. In contrast to Monte Carlo approaches, which suffer from noise and high computational cost in strongly collisional regimes, GUERNICA combines a discontinuous Galerkin discretization in configuration space with a discrete velocity formulation, enabling noise-free, high-fidelity solutions whose cost is independent of plasma-neutral collisionality. Verification and validation are demonstrated using analytic test problems and comparisons against DEGAS2 in one-dimensional configurations. Performance results show favorable scaling with velocity resolution and high throughput on single-GPU hardware. These results are extended to initial two-dimensional simulations in a MAST-U geometry driven by a static plasma background, underscoring the potential of deterministic kinetic neutral modeling for next-generation edge studies. [1] S. Mordijck, Nucl. Fusion 60, 082006 (2020) [2] U. Fantz et al, J. Nucl. Mater. 290–293, 367–73 (2001) Work supported by US DOE under DE-SC0007880.

Non-Local Transport of Vibrationally Excited D2 and its Impact on Divertor Detachment in MAST Upgrade Super-X Plasmas (Mean-field codes)

• Joseph Bryant (Tokamak Energy)

Room: HOW room Molecular processes are known to influence divertor detachment, but conventional effective-rate models assume vibrationally excited molecules remain in local equilibrium and neglect their transport [1, 2]. In this work, vibrationally resolved D2 simulations using the X1EXT molecular database are applied to both simplified divertor-leg and full-device MAST Upgrade Super-X SOLPS-ITER simulations to investigate the role of vibrational transport,electron cooling, and plasma-surface interactions during detachment [3, 4]. In isolated divertor-leg simulations, vibrationally excited molecules were found to survive for tens of microseconds and travel distances of several tens of centimetres before reacting. More than 55% of molecular trajectories exceeded the local electron-temperature gradient length, demonstrating strongly non-local behaviour. This transport enables excited molecules generated in warmer regions to penetrate detached plasmas, enhancing molecular charge exchange (MCX), molecular-activated recombination (MAR), and molecular-activated dissociation (MAD) far from their point of origin. Electron-impact excitation of D2 was also identified as an important power-loss channel, accounting for approximately 10% of total divertor power losses. These excitation losses (“electron quenching”) reduce electron temperature, promote low-temperature plasma-molecular interactions, and accelerate the onset of detachment. The impact of these mechanisms was assessed in full-device simulations of MAST Upgrade Super-X plasmas and compared with spectroscopic measurements. Relative to effective-rate approaches, the vibrationally resolved model predicts detachment rollover at lower upstream density through enhanced MCX-driven MAR and stronger divertor cooling. The simulations reproduces experimentally observed Dα emission profiles, Fulcher-band emission,ionisation-front migration, and target ion-flux rollover more accurately than either AMJUEL or vibrationally unresolved X1EXT simulations. Sensitivity studies show that plasma-surface interaction models have a strong impact on the vibrational state distributions of recycled molecules. The Eley-Rideal recycling model produces an overpopulation of high vibrational states and leads to unrealistically early detachment rollover, whereas the partial-thermalisation approach gives the closest agreement with experimental observations [5, 6]. These results highlight that predictive modelling of detached divertor plasmas requires a self-consistent treatment of vibrational state distributions, non-local molecular transport, and electron quenching, alongside plasma-surface interactions. [1] K Verhaegh, et.al. The role of plasma-atom and molecule interactions on power & particle balance during detachment on the MAST Upgrade Super-X divertor. Nuclear Fusion, 63(126023), 2023. [2] S. Kobussen, et.al. Collisional radiative modelling with improved cross sections to investigate plasma molecular interactions in divertor plasmas. Technical report, Masters Thesis, 2023. [3] J Bryant, et.al. Impact of yacora evaluated molecular effective rate coefficients on detached solps-iter simulations. Nuclear Fusion, 2 2025. [4] D Wunderlich, et.al. Yacora on the Web: Online collisional radiative models for plasmas containing H, H2 or He. Journal of Quantitative Spectroscopy and Radiative Transfer, 240(106695), 1 2020. [5] Anthony J.H.M. Meijer, et.al. Isotope effects in the formation of molecular hydrogen on a graphite surface via an eley-rideal mechanism. Journal of Physical Chemistry A, 106, 2002. [6] K Verhaegh, et.al. Divertor shaping with neutral baffling as a solution to the tokamak power exhaust challenge. Communications Physics, 8, 5 2025

CRFAX: Collisional-Radiative modelling Framework in JAX (Kinetic neutrals)

• Andrei Ludvig-Osipov (UKAEA)

Room: HOW room The plasma behaviour in the exhaust region of tokamaks is strongly influenced by atomic and molecular processes, as experimentally demonstrated on, e.g., MAST-U [1], JET, AUG [2], Alcator C-Mod [3]. Faithful representation of these processes is therefore necessary for robust exhaust simulation capability. However, fully-resolved atomic and molecular reaction systems typically include hundreds and even thousands of excitation and rovibration levels across ionisation stages [4], rendering direct inclusion into exhaust fluid or kinetic codes computationally unfeasible. At the same time, an asymptotic treatment is not applicable, as neither local-thermodynamic nor coronal equilibriums hold in the typical divertor conditions. This, combined with high stiffness due to extreme timescale ranges, warrants for reduced (effective) models accounting for both collisional and radiative processes. CRFAX (Collisional Radiative modelling Framework in JAX), being developed at UKAEA, is a Python package for modelling, analysis, reduction and post-processing of atomic and molecular systems in plasmas. A non-exhaustive list of present capabilities includes reduction methods for linearised reaction systems (a method by Greenland [5] and quasi-steady-state assumption), automatic reduced-model identification method [5], calculation of validity timescales and metrics of reduced models, solvers for linear and non-linear models, as well as calculation of rate contributions (emission spectrum and energy balance). Multi-dimensional scanning, automatic differentiation [6] and gradient-based fitting applied to parametrised reaction systems allow for applications such as tabulation of effective rate coefficients, sensitivity analysis and spectrum fitting. The provenance of CRFAX objects is tracked using W3C Prov data model [7]. The future development of CRFAX is aimed at reduction techniques for non-linear systems [8, 9] and implementation of uncertainty quantification. [1] K. Verhaegh et al., Nucl. Fusion 63, 016014 (2023). [2] M. Bernert et al., Nucl. Mater. Energy 12, 111-118 (2017). [3] A. Mathews et al., Rev. Sci. Instrum. 93, 063504 (2022). [4] K. Sawada, M. Goto, Atoms 4(4), 29 (2016). [5] P. T. Greenland, Proc. R. Soc. A: Math. Phys. Eng. Sci., 457, 1821-1839 (2001). [6] R. Frostig, M. J. Johnson, C. Leary, SysML Conf. 2018, hal-05188750, (2019). [7] L. Moreau et al., PROV-DM: The PROV data model, http://www.w3.org/TR/prov-dm/ (2013). [8] U. Maas, S.B. Pope, Proc. Combust. Inst. 24, 103-112 (1992). [9] S. H. Lam, D. A. Goussis, Int. J. Chem. Kinet. 26(4), 461-486 (1994).

Overview of the GRILLIX code (Turbulence codes)

• Andreas Stegmeir (Max Planck Institute for Plasma Physics)

Room: HOW room An overview of the GRILLIX [1] code is presented, a three-dimensional turbulence simulation code for the self-consistent modeling of the plasma edge and scrape-off layer. The code is based on a global drift-reduced fluid model coupled to a three-moment fluid neutral model. It employs the Flux-Coordinate Independent (FCI) approach, enabling efficient treatment of complex magnetic geometries, and has recently been extended to non-axisymmetric stellarator configurations [2]. Key scientific results obtained with GRILLIX are summarized, including simulations of advanced confinement regimes such as the X-point radiator [3] and Quasi-Continuous Exhaust regime [4], along with their experimental validation, studies of confinement transitions (L–H) [5], and first simulations of island divertor configurations in Wendelstein 7-X. We present the modeling, numerical, and technical developments required to obtain these results. On the modeling side, in particular the self-consistent treatment of electromagnetic effects and extensions to low-collisionality regimes proved essential. While the Flux-Coordinate Independent (FCI) approach provides an efficient framework for accurately capturing strongly anisotropic dynamics, it also introduces several challenges, including the treatment of boundary conditions, the coupling to neutral models, and the achievement of good parallel scalability. Furthermore, the HPC aspects of GRILLIX are discussed. Next-generation supercomputing systems increasingly require efficient execution on GPU-accelerated hardware. We therefore present the porting of the elliptic field solver to GPUs [6], which represents both the most algorithmically complex and the most computationally demanding component of GRILLIX, while also featuring a well-defined interface that enables modular porting. Significant performance improvements are achieved with the GPU-accelerated implementation across different architectures. The presentation concludes with an overview of ongoing work on open problems. In particular, the breakdown of the fluid assumptions may necessitate a kinetic description, for which the GENE-X code is being developed in close collaboration with GRILLIX. To address the challenges associated with boundary conditions, we present an approach that departs from the currently employed immersed boundary method and offers the potential for a more accurate and flexible implementation of boundary conditions within the FCI framework. Finally, lessons learned from the porting of the elliptic field solver have important implications for the GPU porting of the remaining GRILLIX components, which will be discussed. [1] A. Stegmeir, D. Coster, A. Ross et al., PPCF 60:035005 (2018). [2] A. Stegmeir, M. Finkbeiner, C.Pitzal et al., CPC 318:109874 (2026). [3] K. Eder, W. Zholobenko, A. Stegmeir et al., NF 65:096029 (2025). [4] K. Zhang, W. Zholobenko, A. Stegmeir et al., PRL submitted (2026). [5] W. Zholobenko, F. Jenko, K. Zhang et al., PRL 136:075101 (2026). [6] A. Stegmeir, C. Lalescu, M. Lin et al., PASC26 (accepted) (2026)

Porting Hermes-3 to GPUs: Current Status, Challenges, and Roadmap (Research Software Engineering)

• Jony Castagna (UKRI-STFC Hartree Centre)

Room: HOW room Hermes-3, a multi-fluid plasma simulation framework built on BOUT++, is currently being ported to GPUs using the RAJA performance portability layer and CUDA. This work builds upon the GPU infrastructure previously developed at Lawrence Livermore National Laboratory (LLNL) for BOUT++, enabling the reuse of established programming models and optimization strategies. The porting effort is being validated using the Annulus plasma turbulence model and three-dimensional edge plasma simulations of the TCV tokamak. We present the current status of GPU-enabled kernels and preliminary performance results, highlighting key challenges related to memory management, solver portability, and data movement. The roadmap toward full GPU capability, including further physics module acceleration and optimization for future exascale systems, is also discussed.

Old Solver, New Tricks: modernizing and developing the UEDGE codebase in the age of AI (Mean-field codes)

• Andreas Holm (Lawrence Livermore National Laboratory)

Room: HOW room The 2D multi-fluid edge-plasma code UEDGE employs the fully time-implicit NKSOL Jacobian-Free Newton-Krylov (JFNK) solver (P.N. Brown, Y. Saad, SIAM J. Sci. Comput. 11, 1990). UEDGE solves a set of coupled plasma-neutral equations delivering 10–100$\times$ speedup compared explicitly-coupled fluid-plasma/kinetic-neutral code and numerical convergence to steady-state when including magnetic and $\mathbf{E}\times\mathbf{B}$ drift flows. Maintaining compatibility with the NKSOL solver while extending the physics model has shaped the development of UEDGE over three decades. The codebase, tracing its lineage to Braam’s B2 code, has grown organically into a tightly integrated Fortran architecture with nested physics loops with a stencil extending beyond local cells due to higher-order terms and the staggered grid, and order-dependent sections. This heritage of hard-won institutional expertise, encoded in its numerics, presents a challenge to code refactoring efforts: modernization is in perpetual competition with the interdependencies that provides robustness to the solver. UEDGE’s time-tested but long-challenged fluid neutral model is one of its most notable features. The model has employed a 9-point stencil from its inception, as needed for non-orthogonal meshes: a choice since validated as critical for physics accuracy (Dekeyser et al., Nucl. Mater. Energy 18, 2019), dispelling much of the early criticism of fluid neutral models and paving the way for advanced fluid neutral (AFN) models. The focus of AFN model development has turned to capturing kinetic effects and the role of plasma chemistry involving minority species, such as hydrogenic molecules, and there are recent UEDGE development efforts addressing both topics. To capture molecule-assisted reaction pathways, a self-consistent model of diffusive fluid molecules, coupled to the plasma and atoms via effective collisional-radiative rates calculated with the bespoke CRUMPET code, was developed and implemented into UEDGE. To preserve compatibility with the JFNK solver, an expedient vacuum transport model using pre-computed DEGAS2 neutral trajectories through the halo plasma was developed and implemented, avoiding the significant refactoring associated with an unstructured, wall-conformal grid implementation. The inherent uncertainties carried by the effective dissociation rates, arising from the inevitable truncation of the reaction chains, and intricacies of the boundary conditions will be presented and discussed. Several recent numerical and computational developments have focused on improving solver performance. OpenMP parallelization of the preconditioning Jacobian assembly and right-hand side evaluation achieves 80% reduction in wall-clock computational time on 32 CPU threads. The optimal decomposition is non-trivial and Random Forest regression models have been applied to identify the dominant performance factors. A continuation solver, leveraging the NKSOL JFNK solver to trace a continuous path through parameter space, producing a dense set of steady-state solutions, with ~22 s computational time per solution for a benchmark scan, was developed and implemented, enabling the creation of large databases needed for data-driven applications. Regression models trained on solver timings failed to identify meaningful improvements to the solver hyperparameters, and the underlying challenges will be reviewed. The discussion will conclude with an overview of ongoing UEDGE developments. These include modernization efforts to decompose the monolithic Fortran source into modular, Python-level interfaces and potentially reimplement physics submodules in JAX; domain extensions to reach the vessel wall using semi-structured grids and API-like recursive calls; and a multi-layer agentic workflow for autonomous, accessible simulation. A demonstration of simulation capabilities will serve as an opportunity to introduce the Python-installable UEDGE and UETOOLS interface packages, both openly available on GitHub and PyPI. This work was supported by US DOE under contract nos. DE-AC52-07NA27344 and was supported by the LLNL-LDRD Program under Projects No. 25-ERD-014 and23-ERD-015. LLNL-ABS-2020769

A One-Dimensional Integrated Workflow for Tungsten Sputtering, Transport, and Redeposition in the Divertor Region (Kinetic neutrals)

• Tengfei Tang (Agency of Science, Technology and Research)

Room: HOW room Tungsten sputtering, redeposition, and transport in the divertor region are key processes governing impurity migration and contamination of the plasma core in fusion devices. Predictive modelling of these processes requires a coupled description of edge-plasma transport, sheath acceleration, plasma–surface interaction, sputtering, ionization, recombination, and impurity redeposition. In this work, we present a one-dimensional integrated workflow for modelling tungsten impurity production and transport in the divertor. The framework couples Hermes-3 for edge-plasma transport with a kinetic/analytical sheath model, a Monte Carlo impurity ion model, and a Direct Simulation Monte Carlo model for neutral tungsten atoms. Tungsten sputtering yields are calculated using a Monte Carlo ion tracer together with the Eckstein formula. Ionisation and recombination processes are included for tungsten charge states from neutral tungsten to $\text{W}^{20+}$. The model also supports multi-species incident ion populations contributing to tungsten sputtering. The developed workflow provides self-consistent estimates of tungsten redeposition rates, sputtering yields, and charge-state-resolved tungsten impurity distributions. In addition, the one-dimensional model can supply sputtering and redeposition data for two-dimensional impurity transport simulations. Development of a fully self-consistent two-dimensional model is currently underway.

2D transport mechanisms near the X-points in snowflake divertors (Turbulence codes)

• Dominic Power (LLNL / UKAEA)

Room: HOW room Snowflake divertor (SFD) experiments show enhanced radial transport near the X-points, more than can be explained by standard 2D edge transport modelling (i.e. with poloidally uniform diffusion coefficients and without drifts). While 3D effects may be relevant here, several 2D transport mechanisms have been proposed in the literature which could explain this phenomenon: 1) a region of enhanced ExB drifts driven by large poloidal gradients in plasma profiles close to the primary X-point [1], and 2) the ‘churning mode’ [2], a toroidally-symmetric plasma vortex in the region of weak poloidal magnetic field near the X-points. We have investigated these effects through two separate modelling approaches. Firstly, we have carried out interpretive modelling of recent SFD experiments at MAST-U using the code UEDGE with equilibrium drifts and currents included in the model. Secondly, we have developed a model of convective transport in the vicinity of the X-points which includes both the effect of equilibrium drifts and the churning mode. This latter model, implemented in BOUT++, allows us to directly compare the two mechanisms at conditions relevant to current and future SFD experiments. From UEDGE modelling with drifts, we have found that ExB drifts can drive transport across the X-point region leading to activation of secondary strike points (SPs) and reduced peak heat fluxes at the primary SPs, consistent with the predictions in [1]. Direct heat flux measurements at the SPs were not available on the experiments modelled, but the UEDGE results show good agreement in midplane $T_e$, $n_e$ profiles and radiated power profiles in the divertor region. From BOUT++ modelling, we find that the equilibrium drifts effect is dominant when the ratio of plasma pressure at the X-points to magnetic pressure at the midplane is low,$\beta_{pm} \lesssim 1\%$. Above this, electromagnetic effects become increasingly important in driving additional transport. Such conditions are attainable in future SFD experiments at MAST-U and NSTX-U, as well as during ELMs in current devices. This work was carried out under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC5207NA27344. [1] Canal et al., Nucl. Fusion 55 (2015) [2] Ryutov et al., Phys. Scr. 89 (2014)

A large EDGE2D-EIRENE MAST-U database and learned surrogate models for diverse tasks (Mean-field codes)

• George Holt

Room: HOW room Accurate scrape-off layer (SOL) modelling is essential for tokamak heat and particle exhaust studies, supporting experiment interpretation, scenario optimisation, and reactor design. Existing methods require a compromise between speed and fidelity: analytical models are fast but simplified, while comprehensive simulations capture plasma and neutral physics but are too costly for rapid-turnaround applications. This work presents a machine-learning framework for accelerating high-fidelity SOL modelling. A scalable automation workflow for EDGE2D-EIRENE was developed to generate tens of thousands of simulations spanning input power, fuelling and impurity seeding rates and locations, magnetic geometry, and transport coefficients. Neural networks, ensemble methods, tabular foundation models and neural operators were trained for tasks including two-dimensional plasma state reconstruction, derived quantity prediction (e.g., divertor peak heat flux), and boundary condition generation for core-edge integrated modelling. The resulting surrogate models provide millisecond inference while retaining the fidelity of the underlying simulations, enabling high-fidelity SOL modelling for experimental support and broad design-space exploration.