Energetic particle marginal stability profile for HL-2M integrated simulation based on neural network module

A critical gradient model is employed to develop a module of energetic particle (EP) marginal stability profiles in OMFIT integrated simulations for studying EP transport. Currently, each iteration of transport evolution is approximately 10 min in the integrated simulation, whereas, the EP marginal stability profile, which serves as an input in the integrated simulation could take much longer; the reason being a combination of the TGLFEP and EPtran codes is employed in our previous investigation. To reduce the simulation time, the critical gradient is predicted by a neural network instead of the TGLFEP code, and the EPtran code is revised with parallel computing, so that the running time of this module can be controlled to within 5 min. The predictions are in good agreement with previous approaches. The integrated simulation of HL-2M with Alfven eigenmodes transported by neutral beam EP profiles indicates that EP transport reduces the total pressure and current as expected, but could also under some conditions raise the safety factor in the core, which is favorable for reversed magnetic shear and high-performance plasmas.


Introduction
Prediction of energetic particle (EP) distributions in the presence of Alfven eigenmodes (AEs), error field ripples and toroidal geometry induced loss cones is a key issue for burning plasma investigation.In current experimental devices, EPs are mostly generated by neutral beam injection (NBI), and a Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
small fraction of EPs are generated by wave heating, such as ion cyclotron radio frequency.The ratio of the EP pressure to thermal pressure is ∼1/5 [1][2][3], and in some extreme discharges, the ratio can reach 1/1 [4,5].For future fusion reactors, such as ITER, the beta of alpha particle exceeds 0.01 onaxis [6].Thus, the error of EP prediction seriously affects the total pressure, which plays an important role in equilibrium reconstruction and is essential for diagnostics.Accurate EP profiles are also critical for predicting neutral beam (NB) and alpha heating as well as fusion power production.
In the current integrated simulation within the OMFIT framework [7], alpha particle distribution is calculated using the ONETWO code [8] and NUBEAM [9] is employed to calculate the NB deposition.Specifically, in the two codes, EP is calculated using the classical slowing down model [10,11], which typically does not agree with the EP distribution inferred by experiments [2,12] if instability, such as AE, is excited.In DIII-D experiments, the EP pressure profile is obtained by subtracting the thermal pressure from the motional-Stark-effect constrained equilibrium pressure profile [13].Unlike experimental analysis, to design HL-2M operations, we need a method to reliably predict the EP pressure without the benefit of experimental data.The method has to accommodate these features: (i) no parameters from AE or EP diagnostic data are used (we need thermal particle information, geometric parameters and equilibrium profiles, such as q.These parameters can be obtained from the OMFIT integrated simulation platform.);(ii) EP transport induced by instabilities (e.g.AE and turbulence) is included; (iii) the running time is sufficiently short so that multiple iterations of the integrated simulation can be employed to optimize the design.
Firstly, we select a critical gradient model (CGM) as the theoretical basis to establish a module for EP distribution prediction.Theoretically, EP spatial gradient is a primary driver for destabilizing AEs, and there is a threshold of EP gradient for AE marginal stability.When the EP gradient exceeds the threshold, the AE is excited and enhances the EP transport.Since CGM is a linear theory, mode-mode interaction is not considered, namely, the threshold is treated as independent of each toroidal mode number n.Generally, to consider multiple n, we calculate the threshold of each n separately, and select the lowest threshold for the local flux surface.In this paper, we aim to demonstrate that the CGM approach can be accelerated by the neural network (NN), thus enabling more efficient integrated simulations, therefore only n = 3 is considered for HL-2M simulation.Note, CGM is a physics-based reduced model, in which a single unstable AE can cause EP transport due to the breaking of adiabatic invariants by microturbulence scattering [14].Alternatively, multiple unstable AEs at different spatial locations could overlap to cause stochastic diffusion [15].Furthermore, a single AE can also cause orbit overlaps in phase space [16], but the stochastic diffusion requires sufficiently large amplitude of AE.It is still an on-going research topic as to whether EP redistribution is due to a single or multiple unstable AEs.What seems clear so far is that when a threshold is exceeded, the AE induced transport mechanism (whether it is a single mode or multiple modes) will quickly relax the EP profile to a marginal stability state [17].This paper focuses on the relaxation mechanism due to a single mode according to the CGM.
The CGM can only simulate the transport induced by EP driven instabilities, such as AE.Various MHD instabilities (kink mode, tearing mode, etc), which are driven by background plasma pressure or current, cannot be considered in the CGM, even though these modes are sometimes invoked in hybrid scenarios as responsible for the flattening of the central q profile.This occurs with q on-axis close to one.Hybrid modes can also exist with q > 1 using external q profile control.In section 2, for training the NN, some scenarios are generated with q < 1 since equilibria are generated automatically by EFIT code.It is difficult to control the value of q min because EFIT generates an equilibrium based on the input pressure and current profile.Since in our application of the NN module, EP redistribution is only induced by AE while neglecting MHD stability, it is important to recognize that the prediction of CGM is not accurate for scenarios with q < 1. Selectively, we focus on high performance scenarios, such as those in CFETR, so that the q-profile is raised above 1 to avoid unstable MHD such as sawtooth and tearing modes, which could trigger major or minor disruptions.
In our previous work [17], the TGLFEP code was used to calculate the critical gradient and the ORBIT code was used to calculate EP loss cones.The critical gradient and loss cones are input into the EPtran code to calculate EP redistribution.Since the TGLFEP code is a gyro Landau fluid code, the EP is treated as an equivalent Maxwellian distribution.The critical gradient predicted by the Maxwellian distribution is indicated to be a bit larger than that using the slowing down distribution [18].However, the EPtran code employs an anisotropic slowing down distribution to describe EP so that the CGM can simulate EP transport in phase-space.This simulation method, which combines the three codes, considers EP transport induced by AE and turbulence without diagnostic data, but the running time does not satisfy feature (iii).Thus, to reduce the simulation time, we explore a NN to replace the TGLFEP code and select analytical equations of [19] to calculate the particle loss boundary instead of the ORBIT code.In addition, we rewrite the EPtran code with parallel computing.
Zou et al [17] indicates that the coefficients k 1 and k 2 determine the local critical gradient evolution through the function a/L th nEP = k 1 /n EP + k 2 , where k 1 represents Landau damping of the bulk plasma, k 2 is the Landau damping of the EP itself as well as the temperature gradient drive of EP, L th nEP is critical characteristic length of the EP density, a is minor radius, and n EP is EP density.In our previous work, TGLFEP was used to estimate AE stability and calculate the two coefficients.The two quantities are evaluated for each profile iteration step, hence are quite time consuming.A possible approach to improve the efficiency is to develop an NN model, instead of TGLFEP, to first estimate if AE can be excited and second to predict the two coefficients for unstable AE threshold.To confirm that TGLFEP can be used for AE simulation, we have used the NOVA-K code for the verification of the TAE position, frequency and growth rate [20], and the MEGA code for the verification of threshold calculation [17].The input of the NN contains 18 variables, which provide a description of the plasma properties, geometry, magnetic field, and so on for each local flux surface.The details of the input variables are introduced in section 2.1.Focusing on HL-2M device [21], magnetic field (B) minor and major radius (a and R) are fixed.According to the parameter range of HL-2M operation, 150 equilibria are generated by the EFIT code [22], and 8 flux surfaces are selected randomly from each equilibrium, so that the database contains 1200 samples.We use 80% of the samples to train the NN, leaving the rest for validation.The trained model is employed to estimate AE instability and predict k 1 and k 2 for AE unstable region, hence the critical gradient profiles for three operating scenarios with monotonic, weak and strong Y. Zou et al reverse shear q-profile, respectively.The EP profiles predicted by NN are compared with the outputs of TGLFEP.The loss function (the criterion for convergence) is a key parameter for NN construction, and the standard mean square error (MSE) is employed in this paper.Since the selection of the loss function determines the accuracy of the NN prediction, we need to double check that the trained NN is adequate for predicting the coefficients k 1 and k 2 with MSE.Thus, considering the feature of the inverse proportional function governing the marginal stability threshold, we further define a more physical loss function (described in section 2.3) to estimate the sufficiency of MSE.
Finally, we establish a module of energetic-particle (MOE), which contains all critical EP physics and is implemented in the OMFIT framework for studying EP transport in integrated simulations.MOE reads the equilibrium from EFIT, bulk plasma density and temperature from TGYRO [23], and NB source from ONETWO.For a burning plasma, the alpha particle source can be calculated by D-T reaction cross section, which is related to the density and temperature.The output of the MOE is the EP distribution f ( r, E, v ∥ /v ) , and integral quantities such as EP pressure and current, which serve as input into the iteration of the integrated simulation.The equilibrium is updated through changes in MOE, until the AE activities, the EP profile and the transport evolution are all selfconsistent.
This paper is organized as follows: section 2 contains the database, training and prediction of NN.Section 3 describes the method for evaluating EP physics in MOE, and the coupling of MOE with integrated modeling using the OMFIT platform.Finally, a comparison of NN predictions with previous approaches, interesting new findings, discussions and conclusions are presented in section 4.

Machine-learning for critical gradient calculation
In [17], we employ the combination of TGLFEP and EPtran to simulate EP redistribution.The critical characteristic length of EP density, L th nEP , used in CGM is established as an inverse proportional function of EP density.
where k 1 represents Landau damping of the bulk plasma, k 2 is the Landau damping of EP itself as well as the temperature gradient drive of EP, a is minor radius, and n EP is EP density.By using TGLFEP, we adjust L nEP until AE marginal stable condition is satisfied (γ AE+ITG/TEM = γ ITG/TEM , where γ AE+ITG/TEM represents linear growth rate of AE by considering the background plasma, and γ ITG/TEM represents linear growth rate of turbulence) to find the L th nEP (ρ) with a fixed EP density for each radial location.Using the same method, a range of EP densities are scanned to obtain L th nEP (ρ, n EP ), which is fitted as equation (1) to get k 1 (ρ) and k 2 (ρ).Following that, L th nEP could be evolved during EP radial transport in the EPtran code with the input of k 1 and k 2 .However, to obtain L th nEP (ρ, n EP ), the TGLFEP code requires much more computational time than needed for obtaining n th EP (ρ) used in the previous simpler CGM [24,25].In general, the grid of n EP needs to be set at over 20 in order to obtain a fitting curve for equation (1) at each flux surface.In other words, the scanning time is 20 times that of previous simulations.With the radial grid of 50, it is equivalent to running TGLFEP 1000 times.In practice, the plasma profiles vary during the iteration in the integrated simulation to update the background profiles (see figure 10 below), which requires k 1 and k 2 to be recalculated.Obviously, it is prohibitive to employ TGLFEP in integrated simulation because of such a long running time.This motivates us to employ NN to predict the coefficients k 1 and k 2 as an attempt to save computational resources.

Database
The TGLFEP code calculates the function L th nEP (ρ, n EP ), which is fitted by equation ( 1) to obtain the coefficients k 1 (ρ) and k 2 (ρ).The EPtran code only reads the two k values as inputs from TGLFEP.Instead of calculating the coefficients k 1 and k 2 in each iteration by TGLFEP, our goal is to have them accessible as the output of the NN.
Next, we need to determine the input of the NN.For ease of training NN, the input variables should be set as few as possible.Thus, a series of conditions are assumed as follows.By ignoring impurity, i.e., Z eff = 1, the bulk plasma and EP are both deuterium; thus, we only need three species of particles, namely, electrons, bulk ions (D) and EP (D).Note that since TGLFEP is a local stability code, the gradient is treated as independent of the corresponding physical quantity (e.g.density gradient and density are independent) for the given flux surface.For electrons, four variables are needed, i.e., density (n e ), temperature (T e ), normalized density gradient (L ne ) and normalized temperature gradient (L Te ).Here, we temporarily assume n i = n e , L ni = L ne for the input of NN, but the quasi-neutral condition will be satisfied by reducing n i when EP is added to the TGLF simulation.For convenience, we assume that T i = T e and L Ti = L Te .The temperature and temperature gradient of EP, T EP and L TEP , need to be input as well.As mentioned, we scan n EP and L nEP to get k 1 and k 2 , therefore k 1 and k 2 are independent of the values of n EP and L nEP .They are, however, dependent on the background parameters and geometry, which come from the NN inputs.For geometric parameters, we also need to input minor radius (a), flux surface centroid minor and major radius (r/a and R/a), elongation (κ), triangularity (δ), the gradient of elongation (r∂κ/κ∂r) and triangularity (r∂δ/δ∂r), and Shafranov shift.
Here, the squareness of the flux surface is neglected.In addition, magnetic field (B), safety factor (q), magnetic shear and pressure gradient are also needed.Finally, we define 18 independent variables (table 1) for the NN input.Some of the variables, such as p_prime is normalized, to conform with those in TGLF.In this paper, we focus on the HL-2M scenario, so that the magnetic field, minor radius and aspect ratio are fixed (B = 2.2 T, a = 0.65 m, R/a = 1.78/0.65),i.e., parameter 11 and 18 is fixed.Note that the parameter 12 represents the major 18 Minor radius a radius for the local flux surface, which is not fixed due to the Shafranov shift.
Although the gradient is assumed to be independent of the corresponding physical quantities at each local flux surface, a random selection of these parameters without constraints will produce abnormal samples not typically observed in experiments, e.g., a flat density gradient with a large electron density at mid-radius.To avoid including these samples and ensure the satisfaction of the G-S equation, the EFIT code is employed to construct an equilibrium, from which a local flux surface is selected.The corresponding 18 parameters of this flux surface compose a sample.This method can avoid illogical samples that reduce the accuracy of machine learning.To generate a sample, the total pressure and current profiles are assumed according to the HL-2M baseline [26].On the axis, the density is assumed to be within 5-6 × 10 19 m −3 , and the temperature is within 3-4 keV.The density and temperature profiles are taken to be Gaussian functions in the core.To construct the equilibrium for an H-mode, we also add an edge pedestal in the profiles.These pedestals are not sensitive to AEs because AEs are always destabilized only in the core.EP density and temperature profiles are assumed as Gaussian function with n EP (0) = [0.2− 0.25] n e (0) and T EP (0) = [8 − 10] T e (0), so that we can obtain T i = T e , and n i = n e − n EP .The total pressure profile can be obtained by summing the density and temperature profiles of the three species of particles.The current profile is modeled by a polynomial function, and the total current lies within 0.6-1.2MA.With the selected pressure and current profiles, EFIT solves the G-S equation to generate an equilibrium.From the equilibrium, we select a local flux surface (variable 5 in table 1) to determine the other 17 variables.Since AEs are only destabilized in the core region, the radial location is selected with ρ = 0.1 − 0.9.Using this method, we produce a total of 150 equilibria, with figure 1 depicting a range of total pressure and safety factor profiles.The red curves are the upper and lower limits of these 150 equilibria, from which 10 examples are selected and depicted by the black curves.After the equilibria are constructed, we select eight locations in each equilibrium, which would amount to 1200 samples.The radial range [ρ = 0.1 − 0.9] is separated into eight parts with a grid thickness of 0.1, and we randomly select a location within each part to ensure even radial coverage of the samples.These samples are input into TGLFEP to estimate if AE is destabilized.The output equals 1 if AE is destabilized, otherwise the output equals 0. TGLFEP is then applied to calculate the coefficients k 1 and k 2 of equation ( 1) for AE unstable region.This provides a set of code-generated data for training and benchmarking the accuracy of NN generated data.To avoid misunderstanding, the outputs of TGLFEP are called 'target' and the results of NN are called 'prediction' in the rest of this paper.

NN model
We construct an NN with the parameters of table 1, while keeping the magnetic field, minor and major radius fixed.In addition, we only focus on n = 3 AE in this paper.Since TGLFEP is an eigenvalue code, there is no nonlinear interaction between different toroidal mode numbers.Considering other n number modes requires repeating the database construction as described in section 2.1.
The flow chart of NN is depicted in figure 2(a) and the outputs of the NN include two steps.As described, each sample contains 18 inputs of table 1, and 3 'target' quantities calculated by the TGLFEP code, which includes (i) if AE can be excited (0 for stable and 1 for unstable), (ii) k 1 and (iii) k 2 .
The NN algorithm gathers input information from all samples to form a training database from which 3 corresponding output quantities are obtained after proceeding through the 2 steps described in the flow chart.In the first step, the NN predicts the probability of AE being excited.If the probability is larger than 0.5, AE is considered to be unstable; conversely, AE is stable with probability less than 0.5.The results calculated by the NN are called 'prediction'.Note that the output quantities (predictions) are in general different from the input quantities (targets) even for the same 18 input values of table 1.As a first check, we compare the NN stability prediction with the stability finding of the target set.
We shuffle the 1200 sample profiles obtained in section 2.1, using 80% of them for training the NN, and the rest will form a validation set, which is used to assess the trained NN.The architecture of NN adopts a sequential model with multiple hidden layers between input and output.An absence of hidden layer can only reflect a linear relation between input and output.Adding one layer introduces nonlinear relations important for prediction.Furthermore, two layers can obtain more accurate and smoother mapping relations.Going beyond two layers adds complexity that costs much more computational time, and with only 1200 samples does not yield significant  gain.Thus, we choose to use 2 hidden layers for NN construction.Each layer has 128 dimensions, which is depicted by figure 2(b).All layers of the NN use the rectified linear unit (ReLU) activation function [27], which has good performance and ease of training.ReLU has faster convergence than sigmoid and tanh activation functions, and is favorable to avoid the gradient vanishing.A loss function measures the difference between the prediction and the target.The form of the loss function can be set with some flexibility such as conforming to known physics.For a given loss function, the NN is trained to decrease the loss with the training set, and a mapping relation between input and output can be obtained.Next, the validation set is input into the relation to calculate the loss as well.The smaller loss of validation set means that the relation is more accurate.As a standard choice, we employ the MSE as a loss function in this step, and the MSE and the "Accuracy" with training epoch are depicted in figure 3. The loss with MSE can be expressed as: and 'Accuracy' measures the ratio of the correct cases (prediction is the same as the target) to total cases in validation set.Here, one epoch means the NN trains all samples once.The blue and red curve of figure 3(a) depict the MSE of training and validation set.We show that the loss of validation set decreases at first, but increases when epoch exceeds 200, which indicates that over-fitting [28], i.e. model does not generalize well from observed data to unseen data, exists in the training, which could cause prediction diverging from target.This is a common observation in NN training.Thus, In the second step, the NN is constructed to predict k 1 and k 2 , which are the coefficients of equation ( 1).We collect the samples with unstable AE to make a new database.Resembling the above, we select 80% data of this database for training and the rest data construct the validation set to assess the accuracy of the trained NN.The loss function adopts MSE, which can be written as where k 1 and k 2 represent the target values, k1 and k2 represents the predicted value, and N is the number of samples.The MSE decreases to 0.01 after 100 epochs, and reaches 0.005 at 500 epochs (figure 4 As indicated, the prediction is not acceptable with epoch = 50, while the difference between epoch = 500 and 5000 is not significant.Thus, we select epoch = 500 as appropriate for our model.

EP distribution predicted by NN
Using the method described above, we have obtained a model to estimate if AE can be excited and an accurate estimate of the coefficients k 1 and k 2 for local positions with unstable AE.In this section, the EPtran code is employed to calculate EP distributions with the coefficients predicted by NN, and k 1 and k 2 calculated by TGLFEP are directly used to calculate EP profiles for comparison.To avoid confusion, subscript 'N' expresses the coefficients predicted by NN, and subscript 'T' expresses the coefficients calculated by the TGLFEP.In this section, we pick three equilibria from section 2.1 to test the NN, and calculate the EP transport profile by using the predicted k 1 and k 2 .Figure 6 presents total pressure and safety factor for the three equilibria, which are within the range of our database (red curve).The three equilibria have monotonic (blue), weak reverse shear (green) and strong reverse shear (purple) q-profile, respectively.Different from section 2.2, we no longer use MSE as the only loss function to estimate the performance of NN in this section.This is because k 1 and k 2 are the coefficients of an  inverse proportional function, which is key for predicting the critical gradient profile.Since the accuracy of the predicted profile depends on a combination of the coefficients, it might not be sufficient to focus on the separate accuracy of the two coefficients.As an example, figure 7 illustrates two curves (blue and red) have the same MSE as the black curve.The three curves are assumed to be of the form of y = k 1 /x + k 2 .k 1 = k 2 = 1 for black curve, k 1 = 0.8, k 2 = 0.6 for blue curve, and k 1 = 0.6, k 2 = 0.8 for red curve.By equation ( 3), the MSE = 0.2 for both blue and red compared to the black curve, but the blue one is clearly closer to the target in the chosen range.This suggests that other loss functions should also be explored.Generally, n EP is of the order of O(10 −1 ) [10 19 m −3 ]; thus, we define another loss function as where x is selected from 0.1 to 1 with step of 0.1 to represent n EP .The coefficients of k 1 and k 2 are predicted by the new loss function and compared with TGLFEP results in figure 8, where the black curve depicts k 1T and k 2T , and the blue curve depicts k 1N and k 2N .The value of zero on the vertical axis represents AE is stable at the local position.For comparison, the NN results using MSE as loss function are also depicted by blue curves.It can be found that the NN predicted unstable AE region is slightly different from that calculated by the TGLFEP for each case, since the accuracy is ∼0.9 according to figure 3(b).For the three cases, the predictions of k 1 and k 2 show good consistence with TGLFEP calculations in the middle of the minor radius, but they are different at either end of the AE region.Obviously, the prediction error is the largest for the third case, which indicates that the trained NN is not sufficient for the scenario with strong reverse shear, although it could be improved by increasing the number of samples.Figure 8 also shows that, compared with the loss function of MSE, the predictions of k 1 k 2 by the custom loss function (equation ( 4)) are closer to the TGLFEP calculations at most locations.
The coefficients k 1 and k 2 obtained by TGLFEP and NN, EP density profiles (figure 9) are calculated by the EPtran code for the three cases.For each case, the green curve depicts the classical slowing down profile, and the black curve depicts the EP transported profile with the coefficients calculated by TGLFEP.The blue and red curves depict EP profiles using the coefficients predicted by the NN with the loss function of equations ( 3) and ( 4), respectively.The results indicate that both predictions with loss function of equations ( 3) and ( 4) are quite close to the target (black curve).Compared to the target, the mean error is <1.5%, and the difference between the blue and red curves are insignificant, suggesting that selecting MSE as the loss function is acceptable.Certainly, the step of x in equation ( 4) could affect the training result, but it is ignored since the prediction of EP profile is quite close to the target.Here, we only construct the NN to replace the TGLFEP, therefore the trained NN has the same disadvantage as the CGM, i.e., it cannot consider EP transport induced by MHD stability with q < 1.

Module of EP base on OMFIT
In this section, a MOE particle (MOE) is constructed using the above NN model and the EPtran code, which is based on CGM.MOE is added to the OMFIT platform for integrated simulations to study the effect of EP transport on plasma profiles and magnetic equilibrium.

EPtran upgrade
The main equation of EPtran [29] is written as equation ( 5), in which only the last term considers EP transport in pitch-angle space.Therefore, all terms except the last one can be allocated to different processors to calculate EP transport with a given pitch-angle distribution.For the last term, it only needs to receive the value of f EP from adjacent processors to update the pitch-angle modification.By this method, the code can calculate EP transport with different pitch-angles in parallel.For NBI, the source S 0 is obtained from NUBEAM output.But for α particles, S 0 is calculated using density, temperature and reaction cross section.Thus, the model can be employed for burning plasma investigation as well even though our examples are for NB heated plasmas only.
In addition, we previously employed the ORBIT code to calculate the loss cone region in physical and velocity space, which we input into EPtran for inclusion of finite orbit width effect, i.e., particles move across the flux surface due to gradient and curvature drift.Figure 6(b) of [17] indicates that considering the perturbation of AE will increase the size of the loss cone; however, the additional loss cone due to AE only reduces EP radial profile slightly, because the amplitude of AE is δB/B∼O (10 −4 ) in experiments [12,30].Thus, we neglect the enlargement of the loss cone by AE and adopt an analytical method to calculate the loss cone instead of the ORBIT code.Estimation of particle loss only needs equations (3.46) and (3.47) of [19], namely where is toroidal momentum, ψ w is the poloidal flux at the LCFS, B min and B max are magnetic field at ψ w with θ = 0 (low magnetic field side) and θ = π (high magnetic field side), θ is Boozer poloidal angle, 2π ḡ is poloidal current, µ is magnetic moment, and E is particle energy.The two equations determine the boundary of loss particles, including trapped and passing particles.With the new additions described in this paper, the inputs to EPtran no longer require the outputs of TGLFEP and ORBIT, instead, the coefficients of critical gradient is calculated by NN, and loss cone is calculated by analytical equations.The code itself is updated with parallel  3) and ( 4), respectively.computing so that the running time of the MOE is within 5 min.

Module of EP within OMFIT framework
The original OMFIT integrated simulation workflow without pedestal evolution is depicted in figure 10 with blue border.Plasma profiles are evolved using the TGYRO code, in which the turbulent transport flux is calculated by TGLF [31] and neoclassical transport flux is calculated by NEO [32].Using these plasma profiles, the ONETWO code is employed to calculate particle and heating source profiles, including calling NUBEAM for NB heating.In addition, the EFIT code reads the total pressure gradient (P ′ ) and FF ′ (consistent with current) from ONETWO to reconstruct the equilibrium, which is input to TGYRO and ONETWO.On this basis, MOE is added for inclusion of EP transport.In MOE, the inputs of the upgraded EPtran code include (i) density and temperature profiles of ions and electrons, which can be calculated by the TGYRO code; (ii) the magnetic field, current and geometry parameters, which are recorded in the output of the EFIT code; (iii) the source of NBI; (iv) coefficients k 1 and k 2 predicted by the NN.The EPtran code calculates the EP distribution, which is integrated to obtain the transported EP pressure and current.Meanwhile, the thermal pressure and current are calculated using the ONETWO code.The two parts are added to the total pressure and current, which are input into EFIT for equilibrium reconstruction.In the present scheme, EP transport is added in the iteration (described by the workflow with the red boundary in figure 10).Here, we focus on the immediate impact of EP transport on the equilibrium without accounting for changes in the thermal pressure and current, which will be upgraded in the next iteration of the background plasma using ONETWO.Thus, the thermal pressure and current profiles, such as bootstrap current, are held fixed in the following analysis.
An HL-2M H-mode scenario is generated by previous OMFIT integrated simulations (blue border in figure 10) as an initial state.The flux surface, total pressure and safety factor are depicted by blue curves in figure 11.The q-profile is enlarged in the inset.For this hybrid scenario, q (0) is raised above 1 by the off-axis ECCD [33], which is a more controllable way to produce the internal flat magnetic shear profile.Based on this, MOE is employed to calculate EP transport, and the transported EP pressure profile is input into EFIT for equilibrium reconstruction.As mentioned, we do not use TGYRO for further iteration, so that the pressure of thermal particles do not change.We also first assumed that the current sources are unchanged.The simulation results are exhibited by green curves in figure 11.AE activity is localized within ρ < 0.4, in which EP gradient exceeds the critical gradient (indicated by the shaded region in figure 11(d)).Indeed, the EP gradient is close to the critical gradient within ρ = 0.3 − 0.4, so that EP transport is weak in this region.In contrast, EP transport is stronger near the axis.The mode structure calculated by NOVA (figure 11(e)) indicates the mode is a BAE and its location is consistent with the MOE result.The pressure profile varies only near the axis, resulting in a minimal change in the safety factor profile.To enhance the consistency, we next consider both pressure and NB driven current modification due to EP transport.The modified NB current profile is obtained using the modified NB parallel momentum acting on the background plasma, following the standard technique for calculating NBCD [34].The simulation results are highlighted by the red curves.It is shown that the safety factor slightly increases near the axis due to decreasing current, which is consistent with previous DIII-D experiment [5].
For comparison, an L-mode scenario of HL-2M is generated with a higher fraction of EP and strong off-axis AE activity (indicated by the shaded region), so that significant EP transport is expected and could induce more noticeable pressure and safety factor modifications.In figure 12, the initial state is depicted by the blue curve.Similarly, we consider only pressure modification (green), and both pressure and current modification (red) due to EP transport.Since the radial range of unstable AE is large and the AE activity is strong, the pressure profile varies significantly.Beam current decreases throughout the core, therefore the safety factor rises within ρ < 0.5.The comparison between EP gradient and critical gradient is exhibited in figure 12(d), and the unstable region is shaded.Figure 12(e) depicted the mode structure of TAE calculated by NOVA.Compared to the H-mode scenario, the modification of the safety factor profile is larger and extends over a broader range because the fraction of beam current is higher.An interesting observation is that if one could design an H-mode plasma with off-axis unstable AE activity such as the DIII-D discharges discussed in [17], we could in principle create a flat or weakly reversed magnetic shear core.This presents a self-organized phenomenon that could replace the use of off-axis ECCD to produce a high confinement hybrid mode plasma [35].In this paper, we shall note that only modifications of NB driven current is considered as the primary effect.Other current components, such as bootstrap current, are fixed in this initial step.Indeed, the modification of ohmic and bootstrap currents due to EP transport has been illustrated by TRANSP simulation [36].In OMFIT, we can consider these changes in subsequent iterations by updating plasma parameters with TGYRO, which is then used by ONETWO to update the new ohmic and bootstrap current profiles.

Conclusion
OMFIT integrated simulation is improved by including EP transport, which is predicted by CGM.To reduce the simulation time, an NN is employed for critical gradient calculation instead of using TGLFEP, and the EPtran code is upgraded by using parallel computing.With these improvements, the MOE particle can be run within 5 min, which is similar to a single iteration of current integrated simulations.
For the NN, 18 variables are defined by assuming the same profile for ions and electrons, while the impurity is neglected.150 equilibria are produced by the EFIT according to HL-2M parameters, and 1200 flux surfaces are selected from these equilibria to construct a database for NN training.The NN is built by a sequential model with a loss function of MSE, and only single n unstable AE is considered.For the estimation of AE excitation, the accuracy could reach 0.9 after 200 epochs, and increasing epoch could cause over-fitting.For the prediction of coefficients k 1 and for critical threshold, the training of the NN shows good convergence.Three equilibria with different safety factor profiles, i.e., monotonic, weak reverse shear and strong reverse shear q-profile, are employed to validate the trained NN.As a test of sensitivity, we apply a custom loss function to compare with MSE.The results indicate that there is some local distinction between the predicted k 1 and k 2 with the two loss functions, but the transported EP density profiles show only minor differences.The predicted EP density profile with k 1 and k 2 calculated by NN is quite close to the target with k 1 and k 2 calculated by TGLFEP.
Finally, an integrated simulation with the MOE-particle (MOE) is employed to reconstruct the equilibrium modified by AE activities for future HL-2M operation.We add the EP pressure (current) calculated by MOE with the thermal pressure (current) calculated by ONETWO to obtain the total pressure (current).The total pressure and current are input into EFIT to solve the G-S equation.In this paper, we only consider transported EP distribution for equilibrium reconstruction, but the updated equilibrium is not returned to TGYRO for iteration of the background.Two scenarios are used to interpret equilibrium reconstruction due to EP transport.Here, we focus on cases with q min > 1 to avoid low frequency MHD instabilities such as sawtooth and internal tearing modes, which could lead to major or minor disruptions for future burning plasmas.The simulations indicate that the degree of EP pressure reduction is different, which is related to the radial location and range of the AE activities.The safety factor increases in the core with decreasing EP current, especially in the L-mode scenario, suggesting that the reduction of EP current is favorable for producing local reverse shear, which may be suitable for hybrid mode advanced regime operation.The hybrid mode has exhibited good thermal confinement, and the self-organizing effect of AE activities and EP transport on energy confinement is an interesting topic for future studies using the new coupled MOE/OMFIT platform.In the above scenarios, the EP pressure can be of the same order as that from the thermal particles, so that EP transport can result in measurable equilibrium variations.
For the NN training, a narrow range of parameters in the database may cause large errors with extreme input, but the convergence property is not good with a wide range.To solve this, the convenient approach is increasing samples in the database.Besides, considering AEs with multiple n could further enhance EP transport, but it also needs more samples for NN.

Disclaimer
This report was prepared as an account of work sponsored by an agency of the United States Government.Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights.Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof.The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Figure 1 .
Figure 1.Profiles of (a) total pressure and (b) safety factor.The red curves depict the upper and lower limit for the 1200 samples, from which 10 samples are selected (black curves) for example.

Figure 2 .
Figure 2. (a) Flow chart of the NN.The NN estimates if the AE can be excited at first, and calculates the two coefficients k 1 and k 2 for AE unstable location.(b) Schematic for the first step of the NN.

Figure 3 .
Figure 3. (a) Loss (mean square error) and (b) accuracy for AE stability estimation versus training epoch.After 500 epochs, the loss reduces to 0.05 and accuracy reaches 0.85.
(a)).Further training indicates MSE converges even for epoch = 5000, which indicates that there is no overfitting in the second step.The predicted k 1 and k 2 are compared with the target in figures 5(a) and (b), in which the blue, green and red points represent epoch = 50, 500 and 5000, respectively.With increasing epoch, the points approach closer to the diagonal, implying that the predicted value is closer to the target.These points are fitted by a linear regression, and the coefficients of determination (R 2 ) with epoch = [50, 500, 5000] are [0.257,0.532, 0.693] for k 1 and [0.548, 0.773, 0.850] for k 2 .

Figure 6 .
Figure 6.(a) Density profiles of electron and EP.(b) Temperature profiles of electron and EP.(c) Safety factor profile.

Figure 7 .
Figure 7. Sample for three inverse proportional function.The red and blue curve have the same MSE with the black curve, but the blue curve is closer to the black curve in the given region.

Figure 8 .
Figure 8.The coefficients of k 1 and k 2 predicted by NN with loss function of equation (3) (blue) and equation (4) (red).For comparison, TGLFEP results are depicted by black curve.Panel (a) and (b) are for case 1, panel (c) and (d) are for case 2, panel (e) and (f ) are for case3.

Figure 9 .
Figure 9. EP profile comparison.In each panel, green curve depicts initial EP profile with classical slowing down distribution and black curve depicts EP profile with the critical gradient calculated by TGLFEP.The blue and red curves depict EP profiles by NN with loss function of equations (3) and (4), respectively.

Figure 10 .
Figure 10.OMFIT workflow with MOE.The typical iteration is in the blue border, and MOE is in the red border.

Figure 11 .
Figure 11.(a) Flux surface, (b) total pressure and (c) safety factor with classical slowing down (blue) and transported (green and red) EP distribution for HL-2M H-mode scenario.Green curve only considers pressure modification, and red curve considers both pressure and current modification.(d) Comparison of EP gradient and critical gradient.(e) Mode structure.

Figure 12 .
Figure 12.(a) Flux surface, (b) total pressure and (c) safety factor with classical slowing down (blue) and transported (green and red) EP distribution for HL-2M L-mode scenario.Green curve only considers pressure modification, and red curve considers both pressure and current modification.(d) Comparison of EP gradient and critical gradient.(e) Mode structure.

Y.
Zou et al and 12125502, and US DOE Office of Fusion Energy Science under Contract Number DE-SC0017992.The numerical simulations of this paper were carried out in Cluster 02 of SWIP and Tian He Supercomputing Center.

Table 1 .
Variables for NN input.