Numerical tools for burning plasmas

The software stack under development within a European coordinated effort on tools for burning plasma modelling is presented. The project is organised as a Task (TSVV Task 10) under the new E-TASC initiative (Litaudon et al 2022 Plasma Phys. Control. Fusion 64 034005). This is a continued effort within the EUROfusion inheriting from the earlier European coordination projects as well as research projects based at various European laboratories. The ongoing work of the TSVV Tasks is supported by the Advanced Computing Hubs. Major projects requiring the high performance computing (HPC) resources are global gyrokinetic codes and global hybrid particle-magnetohydrodynamics (MHD) codes. Also applications using the integrated modelling tools, such as the Energetic-Particle Workflow, based on the ITER Integrated Modelling & Analysis Suite (IMAS), or the code package for modelling radio-frequency heating and fast-ion generation may require intensive computation and a substantial memory footprint. The continual development of these codes both on the physics side and on the HPC side allows us to tackle frontier problems, such as the interaction of turbulence with MHD-type modes in the presence of fast particles. One of the important mandated outcomes of the E-TASC project is the IMAS-enabling of EUROfusion codes and release of the software stack to the EUROfusion community.


Introduction
The ITER baseline scenario [1] will be characterised by plasma heating dominated by fusion-born alpha particles, and a future DEMO reactor will operate close to ignition, employing auxiliary power only for control purposes. In all present devices, generation of a significant fraction of alpha particles is difficult. This makes understanding of burning plasma physics an urgent point to be addressed through theory and simulations.
In this paper, we describe numerical codes under development for burning plasmas in the scope of the EUROfusion theory and advanced simulation coordination (E-TASC) programme [2]. This development is carried out within the Theory, Simulation, Verification and Validation (TSVV) Task 10, which is a part of the E-TASC initiative dedicated to theory and simulations of burning plasmas.
The TSVV projects, including the one described here, provide an increased level of coordination. They build on achievements and continue efforts of earlier European coordinating structures, such as the EFDA Integrated Tokamak Modelling [3,4], EUROfusion's Workpackage 'Code Development' [5], EUROfusion's Workpackage 'Infrastructure and Support Activities', the EFDA High Level Support Team (HLST) [6], and its continuation into the EUROfusion. Also, a number of past EUROfusion's Enabling Research projects (such as 'Multi-scale Energetic particle Transport in fusion devices', 'Nonlinear 3D simulations of plasma core instabilities beyond magnetohydrodynamics (MHD) in tokamak plasmas', etc) have played a fundamental role for the TSVV Task 10. Another key ingredient in the TSVV structure are the laboratory based research initiatives. For TSVV Task 10, the main contributions include the particle-in-cell (PIC) code development at the Max Planck Institute for Plasma Physics (IPP) and the Swiss Federal Institute of Technology Lausanne (EPFL); nonlinear MHD and hybrid MHD codes as well as theory from the French Alternative Energies and Atomic Energy Commission, French National Centre for Scientific Research, and Italian National Agency for New Technologies, Energy and Sustainable Economic Development; reduced energetic-particle and transport modelling from IPP and the Instituto Superior Técnico of the University of Lisbon. The ongoing work at the TSVV Tasks is supported by the Advanced Computing Hubs (ACHs) [2]. For the TSVV Task 10, the key contributors are the Hubs based at the EPFL and IPP (focusing on the High-Performance Computing aspects), as well as the Core Software Integration Hub at the Poznan Supercomputing and Networking Centre.
Fusion alphas and energetic particles in reactor-relevant plasmas have a unique and crucial role as mediators of crossscale couplings [7,8]. This makes transport in fusion devices a multi-spatiotemporal-scale process and introduces a coupling of various physical subsystems which cannot anymore be considered as isolated parts. A complicating factor here is that the physics interactions mediated by the energetic particles are volumetric rather than through simple local interfaces which makes the scale separation more complex (or impossible) compared to the most other areas of research. As a consequence, predictive analyses based on first principles computations becomes very challenging. Developing a detailed understanding of such complex nonlinear dynamics in fusionreactor plasmas in both tokamak and stellarator geometries is essential. In our TSVV Task 'Physics of Burning Plasmas' (TSVV Task 10), we address the following subsystems of burning plasma physics: (i) turbulence, turbulent structures (zonal flows, avalanches, profile corrugations), and transport (ii) energetic particles and Alfvénic modes; phase-space structures (holes, clumps, auto-resonant 'chirping' dynamics) (iii) MHD stability and profile evolution (sawteeth, tearing instabilities).
In our project, we focus on the core-plasma dynamics since it is the core of the plasma which is expected to burn. For the subsystem coupling, we note the following examples.
(i) Turbulence leads to anomalous transport but it is affected by emergent structures such as zonal flows and avalanches. Energetic particles can destabilize Alfvénic modes and be transported by these instabilities. (ii) The nonlinear dynamics of destabilized Alfvén Eigenmodes may result in a zonal flow generation [9,10] which can affect (e.g. stabilize [11]) turbulence. (iii) MHD phenomena, such as sawteeth may lead to strong particle redistribution (for the both bulk-plasma and energetic species) and hence affect the Alfvénic dynamics and fusion. On the other hand, the sawteeth can be affected by the energetic particles playing a stabilizing role on the underlying instabilities but also leading to the appearance of particularly strong 'monster' sawteeth [12].
Strong coupling between different subsystems of fusion plasma dynamics calls for a unified framework which includes all these pieces self-consistently on the same footing. Many features of the subsystems relevant for the core burning plasma are kinetic and global, and many links between these subsystems are kinetic and global, too. This implies that the global gyrokinetic method is the minimal inclusive approach which may take into account the relevant complexity and mutual dependencies. Note that the above list of items is not exhaustive. One aspect, not included in the list, concerns global kinetic turbulence self-organization, e.g. the transport barrier formation [13], with a strong profile evolution coupled to the neoclassical physics and controlled by the particle and energy sources. This kind of physics will be affected by the energetic ions in burning plasmas and shall be addressed using fluxdriven codes such as GYSELA [14] as well as ORB5 [15,16]. Presently, GYSELA is not a part of the TSVV Task 10 program. Up to now, this code has mainly been used in the electrostatic regime.
In our project, we employ the global electromagnetic gyrokinetic codes ORB5 [15,17] and EUTERPE [18] which can be used in tokamak and stellarator (in the case of EUTERPE) geometries. Using these gyrokinetic codes, we study electromagnetic turbulence in different regimes [19][20][21][22], consider meso-scale couplings of the turbulent dynamics to the MHD-like tearing and kink modes [22], and study interaction with the Alfvénic instabilities in the presence of energetic ions [23][24][25][26][27][28][29]. Note that ORB5 and EUTERPE are included in other TSVV Tasks as well. In our Task, the focus is on the electromagnetic applications of these codes.
Global gyrokinetic simulations are very expensive and time consuming. Therefore, developing reduced models, supporting and supported by the first-principle simulations, is indispensable. There are different degrees of reduction possible. To consider the nonlinear MHD time scale, such as the sawtooth cycle or a fishbone burst, a hybrid MHD approach can be used. Here, the energetic-particle dynamics is considered using the full-kinetic or the gyrokinetic method whereas the bulk plasma is described solving the MHD equations. For this purpose, we use the nonlinear-MHD code XTOR-K [30][31][32], which solves the full kinetic problem for the energetic particles. Using XTOR-K, the sawtooth dynamics can be considered in the presence of energetic ions. The latter can modify the stability properties of the precursor internal kink instability but can also be transported by the sawteeth, altering their couplings to the Alfvénic continuum and Eigenmodes, as it has been seen in recent JET experiments [33]. Also fishbones or energetic particle modes may lead to a non-diffusive transport of energetic particles [8,34]. In burning plasmas it can strongly affect the fusion reactivity and the overall system performance. For these instabilities, the main nonlinearity is on the energetic particles whereas the bulk plasma can be described with a good accuracy solving the linear full MHD equations, i.e. neglecting the wave-wave nonlinearity. This approach is realized in the HYMAGYC code [35] whereas the HMGC code [36] solves the nonlinear reduced-MHD equations to describe the bulk plasma. In stellarator geometry, the hybrid particle-MHD code CKA-EUTERPE [37] is used. An important feature of HYMAGYC is its almost complete ITER Integrated Modelling and Analysis Suite (IMAS)-enabling based on a machine-generic data dictionary [38]. IMAS [38] is a key standardization tool for E-TASC and EUROfusion data and code integration [2] based on the physics Data Dictionary and the Interface Data Structures (IDSs) [38] typically describing tokamak subsystems (such as equilibrium or heating). HYMAGYC can be run as an embedded IMAS actor [38] in an integrated environment. New IDS structures have been identified during the IMASenabling of HYMAGYC which should be included in the next IMAS releases. These new IDS structures can be used in future by other energetic particle and MHD stability codes. Note that most of the codes developed within the TSVV Task 10 implement unified interfaces based on the IDS structures. One exception is EUTERPE where IDS interface implementation is complicated by the stellarator geometry (IMAS is focused on tokamak applications).
Hybrid-MHD simulations are more affordable than the full gyrokinetic approach but they can also become prohibitively expensive on the full discharge time scale. Here, integrated modelling based on the IDS coupling of several codes becomes the only possible approach. The key part of this setup, further developed in our project, is the new Energetic-Particle Workflow. It is a pythonic framework which couples equilibrium codes, such as the CHEASE [39] or HELENA [40], transport codes, such as the european transport solver (ETS) [41], heating codes (e.g. NUBEAM [42]), and the LIGKA-HAGIS codes [43,44] which describe the nonlinear fast-ion dynamics on a short time scale. The resulting energetic-particle distribution function can be used on the next time step to advance the plasma profiles and the equilibrium, if needed. To speed-up the computations, pre-calculated phase-space zonals structures for a set of Alfvénic modes at fixed amplitudes (somewhat similar to the kick-model approach of [45]) can be employed to calculate phase-space dependent energetic-particle diffusion coefficients which can be used in the heating codes.
The new Energetic-Particle Workflow has only been developed for tokamaks so far. In future, this framework will be extended to stellarator geometries using the CONTI code [46] and the CKA-EUTERPE code [37] as the main ingredients. The CONTI code computes shear Alfvén continua in general toroidal geometry based on an Variational Moments Equilibrium Code (VMEC) [47] equilibrium. The CKA-EUTERPE code can be used to study the nonlinear energetic-particle dynamics in stellarator geometry interacting with a prescribed set of Alfvén Eigenmodes. In future, these tools can be combined in an automated manner resulting in a similar workflow for stellarators.
Although most of the heating in a burning plasma should be provided by fusion-born alpha particles, auxiliary heating will still be of importance before the plasma ignition is achieved and also later for discharge control. In our project, we employ the SCENIC framework [48], which describes ICRH plasma heating and energetic-ion generation via a combination of three codes: LEMan [49] solving the hotplasma wave-propagation equations, VENUS-LEVIS [50] tracing energetic-particle orbits and using Monte Carlo kick operators [51] for collisions and interaction with the wave field, and ANIMEC [52] solving for a plasma equilibrium with an anisotropic energetic-ion pressure tensor. The role of SCENIC in the TSVV Task 10 is to provide an energeticparticle distribution function needed as input by the stability codes.
In summary, our codes are employed to provide a selfconsistent description of the mutual interaction of energetic particles with MHD modes and turbulence, as well as their interplay with the kinetic plasma profiles in both tokamak and stellarator geometries. In future, they will be used to develop strategies optimizing deposition of the fusion alpha energy to the bulk plasma and improving the reactor performance.
This paper is organized as follows. In section 2, we detail the description of the codes developed by the TSVV Task 10 and address the main simulation results. In section 3, we summarize our conclusions.

Gyrokinetic codes
ORB5 and EUTERPE are global codes solving the nonlinear electromagnetic gyrokinetic Vlasov-Maxwell system of equations [53] in tokamak and stellarator geometries using a PIC scheme. The codes assume multiple gyrokinetic ion species and drift-kinetic electrons. The choice of the ion species mix is completely general and can include energetic ions (e.g. alpha particles), impurities (Ar, Be, etc), isotope fuels (Deuterium and Tritium), or Helium ash. The same freedom of choice applies also for the leptonic component. For example, ORB5 has been used to simulated electronpositron pair plasmas [54]. All the species are treated on the same footing. The gyrokinetic ions imply that the codes can be used at all wavelengths ranging from the global to the ion gyro-radius scale and below. Several types of collision operators are available in the codes ranging from a simple Lorentz collision operator [55] to more general operators based on the Rosenbluth potentials [56]. A Fourier filter in the toroidal and poloidal directions as well as various control variates and noise reduction techniques [15] enable simulations with good signal-to-noise ratio at a limited numerical cost. For electromagnetic simulations, the control-variate and the pullback mitigation techniques [57,58] are available. As a subset, ORB5 and EUTERPE include fluid-electron and fluid-bulk kinetic-energetic-particle models [59] as well as the CKA-EUTERPE model [37]. ORB5 runs on both CPUs and GPUs [17]. EUTERPE has been extended to GPUs recently. Code testing in a production environment is ongoing. The codes are well benchmarked [60][61][62][63][64] against other similar codes and analytical predictions and show good scalability (see [15,17,65] addressing the performance of the codes on many-and multi-core architectures). EUTERPE is specialized for simulations in three-dimensional (stellarator) geometry using a plasma equilibrium calculated with the VMEC [47] or GVEC [66] codes. Two coordinate systems are used in EUTERPE: a system of magnetic coordinates (PEST) for the field solver and cylindrical coordinates for the equilibrium and to push the particles. The magnetic coordinate system employed for the field solver facilitates handling of complex shaped geometries, typical for stellarator plasmas. The change between the coordinate systems is done using linear interpolation. The quality of the interpolation (the corresponding grid size) should be good enough for a numerical accuracy of the codes. However, the interpolation is needed in the real space only and appears to be less problematic than in semi-Lagrangian codes [67] where the distribution function has to be interpolated in the phase space. The equations for the perturbed field are solved using the PETSc [68] library. The integration of the Vlasov equation is usually done using an explicit fourth order Runge-Kutta scheme. ORB5 can use realistic CHEASE [39] equilibrium.
Global gyrokinetic simulations of electromagnetic turbulence have been performed using ORB5 and EUTERPE in presence of energetic particles and global MHD-like (tearing) modes at different beta values in realistically shaped tokamak and stellarator geometries [21,22]. A transition from electromagnetic ITG to KBM regimes has been observed in ASDEX-Upgrade and ITER using ORB5. For stellarators, such a transition between the different turbulence regimes has been considered using EUTERPE in various configurations of W7-X as well as new optimized stellarator geometries. First ORB5 results (Toroidal Alfvén Eigenmodes, TAE modes) are available for JET plasmas, see figure 1 for an example. Non-adiabatic chirping dynamics of toroidal Alfvén eigenmodes (TAEs) and energetic particle modes (EPMs) has been studied [70] using HMGC (non-perturbative hybrid particle-MHD approach) and ORB5 (fully gyrokinetic for all particle species, including the electrons). The results compare well to each other and demonstrate numerically the importance of the 'Trap-Release-Amplify' mechanism [71] for the phasespace evolution of chirping instabilities. This mechanism has been studied analytically [72] both in the context of burning fusion and space plasmas, where a similar theoretical framework based on the Schwinger-Dyson equation can be applied to describe the evolution of phase space zonal structures [8,34]. These structures can be understood as the evolution of the system through neighbouring nonlinear equilibria [73]. The Hamiltonian-mapping diagnostics has been applied in HMGC to study the phase-space dynamics of the particles in the chirping modes which have also been studied with the fully gyrokinetic EUTERPE and the reduced CKA-EUTERPE model [37]. For ORB5, a dedicated finite element representation of phase space zonal structures has been introduced [74] which can be used for various purposes, including chirping dynamics. Fully numerical energetic-particle distribution functions (e.g. generated by the NBI) have been applied in ORB5 [26]. Validation of ORB5 using realistic distribution functions for ASDEX-Upgrade discharge #31213 has been extended to the nonlinear regime [25]. It has been shown that realistic energetic-particle distribution functions are essential to reproduce quantitatively the experimental results. Multi-scale analysis of global electromagnetic instabilities has been carried out using ORB5 in ITER Pre-Fusion Power Operation plasmas as a part of an international effort under the EP-ITPA umbrella [28]. It has been shown that Alfvénic instabilities can be excited by the thermal ions in ITER plasmas even in the absence of energetic ions providing the toroidal mode numbers (the diamagnetic frequency) are high enough [28].

Hybrid particle-MHD approach
XTOR-K solves nonlinear 3D extended MHD equations coupled self-consistently with a full-f full-orbit or guiding centre PIC time advance for selected populations of kinetic particles. The implicit fluid part in the time advance is inverted using a pre-conditioned matrix-free GMRES solver [30]. A substantial effort has been undertaken in the last two years to parallelize the inversion of the physical pre-conditioner using a SPIKE LU solver which solved the problem of small time steps in the saturated kink regime. It has resulted in an overall factor 2.5 speed-up of the code. Another factor 2 speed-up has been achieved due to particle sorting which strongly reduces cache missing for moment depositions. Successful simulations have been performed for the TAE mode destabilized by the energetic particles (the 'ITPA benchmark' [75,76]). For a correct simulation of sawtooth oscillations, the model in the fluid part of the code had to be improved and benchmarked from resistive MHD (used in previously done fishbone studies [32]) to a more general two-fluid model. The model differs from the one in the old two-fluid MHD XTOR-2F because bulk ion and electron temperatures evolve separately in XTOR-K. Also, the components of the MHD equations parallel and perpendicular to the magnetic field are treated differently. All these physical model and performance improvements were necessary before starting sawtooth cycle simulations in the presence of fast ions.
HMGC describes the thermal plasma by nonlinear reduced O(ε 3 0 ) visco-resistive MHD equations, with ε 0 = a/R 0 being the inverse aspect ratio (with a and R 0 the minor and major radius of the torus, respectively), evolving the fluctuating electrostatic field ϕ and the perturbed vector potential component, parallel to the equilibrium magnetic field, A ∥ (low-β limit, with β being the ratio of the plasma pressure to the magnetic pressure), in the zero pressure limit. These equations are closed by a pressure term related to one or more particle species (e.g. fast ions), computed by solving the nonlinear gyrokinetic Vlasov equation, in the k ⊥ ρ H ≪ 1 limit (driftkinetic limit), k ⊥ being the perpendicular (to the equilibrium magnetic field) wave vector of perturbed fields. HYMAGYC studies energetic-particle driven Alfvénic modes in general high-β axisymmetric equilibria, with perturbed electromagnetic fields (electrostatic potential ϕ and vector potential A) fully accounted for. The thermal plasma is described by linear, resistive, full MHD equations and the calculations are carried out in flux coordinates (s, χ, φ), s = √ 1 − ψ/ψ a (with s = 0 in the centre and s = 1 at the plasma boundary), ψ is the poloidal flux function, χ is a poloidal angle and φ is the geometrical toroidal angle. The perturbed variables are then expanded using Fourier series in χ and φ, and finite elements in s. When considering two dimensional, axisymmetric equilibria, the solution for each toroidal mode number n is independent from the others, and can be solved one by one. After discretization, one obtains matrices of the coefficients that are block tridiagonal, with eight blocks for each flux-like radial coordinate point considered. The dimension of each elementary block scales with the square of the number of poloidal Fourier components considered (which, in turn, scales as N pol ∼ n∆q, with ∆q = q max − q min , q being the safety factor, and q min , q max the minimum and maximum values of the safety factor considered in the equilibrium, respectively) times the number of equations solved for each flux point. Moreover, also the number of flux points that have to be retained in the simulation, scales with the toroidal mode number n, being the characteristic radial extension of each poloidal Fourier component inversely proportional to n: thus, in turn, the memory requirement to store the matrices of coefficients for solving the linear MHD system scales as n 3 . The MHD field solver originates from the code MARS [77], which has been transformed from an eigenvalue solver to an initial-value one. A parallel version of this field solver has been implemented in the recent past, thanks to the EFDA-HLST projects PARFS [78,79] and PARFS2 [80,81], thus overcoming single-node memory limitations to the maximum toroidal mode number that can be considered in the simulations. The energetic-particle population is described by the nonlinear gyrokinetic Vlasov equation valid for k ⊥ ρ H ∼ 1 (FLR effects properly retained). HYMAGYC is fully interfaced with the European Integrated Modelling Framework [4] data structure and benchmarked against HMGC in the suitable limits [35]. HYMAGYC has also been compared with the hybrid kinetic-MHD MEGA code and fully-gyrokinetic ORB5 code using the NLED-AUG test case to study Alfvénic modes driven by energetic particles [62]. Both HMGC and HYMAGYC are written in Fortran90 and have been parallelized for distributed and shared memory architectures.
For the HYMAGYC code, a strong effort has been put into further IMAS-enabling. In the frame of the European integrated modelling effort [4] and, then, in IMAS [38], a platform for integrated tokamak modelling has been developed using the graphical user friendly Kepler software manager [82]. Within Kepler, interchangeable physics modules, named actors, are used to construct complex Workflows. The minimal HYMAGYCIMAS Kepler Workflow [38] works correctly on the EUROfusion Gateway machine. The UALInit actor is used to extract the required IDS data and pass them to the next HYMAGYCIMAS actor which contains the actual HYMAGYC code. The resulting output of the HYMAGYC simulations is written to the IMAS Database using the UALSliceCollector actor. An extensive optimization of the MHD_LINEAR IDS has been carried out for parallel job executions.

Integrated modelling
The fully IMASified Energetic-Particle Workflow has been developed and implemented for discharge modelling in ASDEX-Upgrade, JET, and TCV. Predictive studies for JT-60SA and ITER have been performed. The Energetic-Particle Workflow is a pythonic framework coupling transport codes (such as ETS), heating or NBI codes (such as NUBEAM), and self-consistent energetic-particle codes (such as LIGKA/HA-GIS) via interfaces based on the IDS data structures. LIGKA [43] is a linear gyrokinetic eigenvalue solver. HAGIS [44] follows the particle guiding-centre orbits providing a selfconsistent model for nonlinear wave-particle interaction. The framework has been applied, for example, to the slow LH transition in the ASDEX-Upgrade discharge #39681 in presence of TAE instabilities where both bursty and steady-state phases have been observed. The Energetic-Particle Workflow gains an increasing popularity in EUROfusion's experimental community [83,84].
As a first step towards a fully integrated simulation of burning plasmas, part of the Energetic-Particle Workflow has been integrated into the transport solver ETS [41] by adding the equilibrium (HELENA) and MHD (LIGKA) Kepler actors in the MHD composite actor. The linear MHD spectrum can now be computed consistently in the convergence loop. To prepare for the integration of more demanding hybrid-kinetic-MHD or reduced model codes for MHD stability assessment in the presence of energetic particles, the IMASified NBI suite of codes BBNBI/ASCOT/AFSI were revised and upgraded to output in IMAS all the energetic particle distributions relevant for DT burning with sufficient resolution. A first implementation in the kick model limit [45] of the general phase space zonal structure transport equations [85] based on elements of the Energetic-Particle Workflow has been finalised and is currently benchmarked.
In stellarator plasmas, numerical tools for computing ICRH heating and the resulting energetic-particle distribution function have been further developed. The SCENIC code package [48] integrates the ANIMEC [52] (3D MHD equilibrium solver), LEMan [49] (linear full-wave global solver), and the fast-particle code VENUS-LEVIS [50]. Predictive simulations of radio-frequency heating and fast-ion generation have been carried out [48,86,87] in Wendelstein 7-X using the SCENIC framework. Finally, a numerical code describing the formation of shear Alfvén continua in stellarator plasmas in presence of magnetic islands has been developed [88]. Continuum calculations play an important role in the reduced modelling of energetic-particle interaction with Alfvénic instabilities and transport.

Conclusions
In this paper, we have described the software stack under development within a European coordinated effort on tools for burning plasma modelling. The project is organised as a Task (TSVV Task 10) under the new E-TASC initiative [2]. It is a continued effort within the EUROfusion inheriting from the earlier European coordination projects as well as research projects based at various European laboratories. The ongoing work of the TSVV Tasks is supported by the ACHs. The major projects requiring the high performance computing resources are the global gyrokinetic codes ORB5 and EUTERPE. Nonlinear gyrokinetic simulations are performed to address the interaction of fast ions with turbulence and in presence of MHD-like perturbations. Also, we perform linear and nonlinear gyrokinetic simulations to study the mutual influence of MHD/EPM instabilities and energetic ions.
Interaction of energetic particles with the MHD/EPM instabilities can also be addressed using the hybrid particle-MHD codes such as HMGC, HYMAGYC, and XTOR-K. These codes are developed to investigate the role of a large population of fusion alphas and other suprathermal particles on the global MHD stability limits of the plasma including the impact on the global β-limit of the low-n kink modes stabilization and determination of the sawtooth period in presence of the kinetic effects.
Taking advantage of the IMAS standard, we integrate the hybrid particle-MHD code HYMAGYC into the EUROfusion software eco-system where it can readily be used to analyze 'off-line' some temporal snapshots of a transport simulation or even be directly coupled to a transport code, such as the ETS. Such code combination can be used to address selfconsistently the mutual interaction of MHD/EPM instabilities driven by fusion alphas and other suprathermal particles with the respective deposition profiles and, consequently, with the bulk density and temperature profiles evolution.
On the faster-execution lower-fidelity level, we develop integrated modelling tools which include the IMAS-based Energetic-Particle Workflow (based on HAGIS-LIGKA) and the SCENIC code package computing radio-frequency heating and fast-ion generation. The integrated modelling enables exploration of active strategies to optimize the deposition of alpha particle energy aiming at maximization of the fusion power yield. Furthermore, it enables modelling of burn control through auxiliary heating and fuelling strategies as well as prediction of current profile consistent with bootstrap contributions from pressure profile and fast particles. Reduced models of Alfvénic stability and nonlinear dynamics are being developed for use in the integrated modelling and systems codes addressing Tritium burn-up rates and core plasma Helium content.
In future, we will apply the software stack described in this paper at all fidelity levels to burning plasmas expected in ITER [1], SPARC [89], and also in stellarator-reactor configurations, such as the HELIAS reactor [90], as well as in a DEMO reactor [91]. The software stack developed at the TSVV Task 10 is being made available on EUROfusion's GitLab and on EUROfusion's Gateway which play the role of a central repository and a central 'device' for EUROfusion's modelling effort [3]. The training for the Energetic Particle Workflow is planned for EUROfusion's user community.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Acknowledgments
This work has been carried out within the framework of the EUROfusion Consortium, partially funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200-EUROfusion). The Swiss contribution to this work has been funded by the Swiss State Secretariat for Education, Research and Innovation (SERI). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union, the European Commission or SERI. Neither the European Union nor the European Commission nor SERI can be held responsible for them. This work was supported in part by the Swiss National Science Foundation. We acknowledge PRACE for awarding us access to Marconi100 at CINECA, Italy, and to Joliot-Curie at GENCI@CEA, France. Finally, we thank the organisers of the 3rd Fusion HPC Workshop and, in particular, Edilberto Sanchez, for suggesting this paper.