An adjoint formulation leveraging a physics-informed neural network (PINN) is employed to advance the density moment of a runaway electron (RE) distribution forward in time. A distinguishing feature of this approach is that once the adjoint problem is solved, its solution can be used to project the RE density forward in time for an arbitrary initial momentum space distribution of REs. Furthermore, by employing a PINN, a parametric solution to the adjoint problem can be learned. Thus, once trained, this adjoint-deep learning framework is able to efficiently project the RE density forward in time across various plasma conditions while still including a fully kinetic description of RE dynamics. As an example application, the temporal evolution of the density of primary electrons is studied, with particular emphasis on evaluating the decay of a RE population when below threshold. Predictions from the adjoint-deep learning framework are found to be in good agreement with a traditional relativistic electron Fokker–Planck solver, for several distinct initial conditions, and across an array of physics parameters. Once trained, the PINN thus provides a means of generating RE density time histories with exceptionally low online execution time.
I. INTRODUCTION
The treatment of runaway electron (RE) dynamics is a critical component in the integrated description of tokamak disruptions1,2 and startup scenarios.3,4 While first-principle kinetic solvers are available to evaluate the transient RE dynamics during these phases of a tokamak discharge,5–9 such approaches are computationally intensive, thus, sharply limiting their application. This has led to the development of several reduced models of RE formation that may be conveniently coupled to magnetohydrodynamic (MHD) codes.10–12 These reduced models, however, are often based on the assumption of a slowly evolving background plasma and are, thus, unable to accurately capture highly transient phases of RE generation or decay. An exception to this trend includes reduced models of hot tail RE seed formation,13–15 where RE formation during a rapid thermal quench is estimated. Such models, however, have to date been developed specifically for the thermal quench, and, thus, are not applicable to other transient phases of a disruption, and do not self-consistently include the presence of large-angle collisions. The amplification of the RE seed by large-angle collisions, the so-called RE avalanche,16 is instead incorporated by a separate module17–19 that does not directly account for the transiently evolving momentum space distribution of REs.
The present work, together with Paper II,20 aims to address these limitations by employing an adjoint-deep learning framework that will allow the transient dynamics of REs to be evolved with an online execution time that is orders of magnitude less than a traditional relativistic Fokker–Planck solver. To achieve this, we will first focus on the solution to an adjoint problem,14,21–23 rather than directly solving the relativistic Fokker–Planck equation. A distinguishing feature of this approach is that a single solution of the adjoint problem can be used to evolve the RE density forward in time for an arbitrary initial momentum space distribution. Thus, once a solution to the adjoint problem has been identified, the impact of distinct initial momentum space distributions on RE density evolution can be rapidly evaluated for a given set of physics parameters. The second component of the framework will be the use of a physics-informed neural network (PINN)24 to identify a parametric solution to the adjoint of the relativistic Fokker–Planck equation. Physics-informed neural networks are a form of deep learning, whereby physical constraints are embedded in the training of a neural network (NN). While this approach has gained prominence due to its ability to integrate physical constraints with data,25 or treat inverse problems,26 for cases where a closed system of equations is available, PINNs can be used to learn the solution to a PDE across a broad range of physics parameters in the absence of data.27,28 Further noting that while the offline training time of a PINN is long, its online inference time is orders of magnitude faster than a traditional solver. The PINN will thus provide a rapid surrogate of the physical system of interest. Combining the above two components, the resulting adjoint-deep learning framework will be able to generate time histories of the RE number density at a fraction of the cost of a traditional solver, for an arbitrary initial momentum space distribution of REs and across a range of plasma conditions. The present paper, along with Paper II,20 will focus on developing the fundamental components of this adjoint-deep learning framework, while future works will focus on employing it for describing tokamak disruptions.
As specific example problems, we consider the decay of a primary RE distribution, while Paper II20 will treat the avalanche amplification of an initial RE seed. While large-angle collisions are essential to providing a description of RE generation during the current quench phase of a disruption, an accurate treatment of RE decay is essential when describing the decay of plasma current29 during the RE plateau. This decay rate will be inferred as a function of electric field strength, effective charge, and the strength of synchrotron radiation. The RE decay rate is shown to have a strong nonlinear dependence on the strength of the electric field, in contrast to commonly employed analytic expressions for RE growth and decay such as Ref. 17, which indicate a linear dependence. Furthermore, we find that characterizing the decay of REs using a simple decay rate, as is often done in current applications, is often inaccurate either due to the long time required for the RE distribution to settle into a well defined exponential decay or due to the RE density exhibiting temporal dynamics that cannot be fit by an exponential, even in the long time limit. While the present paper does not include large-angle collisions, Paper II20 finds that the decay of the overall RE distribution is well approximated by the decay of the primary population neglecting large-angle collisions. The present estimate, thus, provides an efficient means of evaluating the decay of a RE population when below threshold. As a further application, time histories of the RE density evolution are computed for distinct initial momentum space distributions of REs. It is shown that the adjoint-deep learning framework is able to accurately resolve these time histories for RE distributions with sharply different initial pitch distributions. The specific time histories predicted by the adjoint-deep learning framework are compared with those predicted by the particle-based RE solver RAMc,7 with good agreement evident across a range of initial momentum space distributions.
The remainder of this paper is organized as follows. Section II discusses the adjoint formulation of the relativistic Fokker–Planck equation. A brief description of PINNs is provided in Sec. III, along with a discussion of the tailored NN architecture employed in this paper. The temporal evolution of the number density of primary electrons for distinct initial momentum space distributions of REs is given in Sec. IV, along with evaluating RE decay for a range of plasma parameters. A brief discussion along with conclusions are given in Sec. V.
II. ADJOINT FRAMEWORK FOR RELATIVISTIC ELECTRON EVOLUTION
A. Adjoint of the relativistic Fokker–Planck equation
The physical interpretation of the RPF can be made clear by considering the orbits of several electrons initialized at different pitches as shown in Fig. 1. These orbits have been computed using characteristic equations derived from Eq. (1), where pitch-angle scattering was neglected leading to deterministic orbits. The parameters ( and ) were selected such that the system is slightly above the avalanche threshold, where effects from synchrotron radiation are most pronounced. After integrating forward in time until [see Fig. 1(a)], all of the electron orbits have slowed down, except an orbit with an initial pitch of . This last orbit was initialized with a pitch that put it in a region with (below the dashed black curve) allowing it to be accelerated. It is useful to note that the rate at which electrons in the region decelerate depends sensitively on their pitch, with electrons near decelerating the most rapidly due to synchrotron radiation being maximal in this region, where the presence of synchrotron radiation results in the region with contracting as p increases. At a later time [ Fig. 1(b)], the majority of electron orbits continue to slow down, but are also tilted toward the negative pitch region. This is due to the electric field pushing the electrons toward the axis. Finally, by , all of the electrons have either been accelerated through the surface, or decelerated to the bulk (i.e., since we have assumed a cold bulk plasma). Those that have been accelerated through the surface [which could be taken to be the upper momentum boundary] will be treated as REs according to Eq. (10). In contrast, two orbits in Fig. 1(c), the purple and red orbits initialized with have fallen back to the bulk, and are, thus, deemed to not run away. Thus, for the case where pitch-angle scattering is neglected, a separatrix will form delineating runaway and non-runaway orbits.
Example deterministic electron orbits. Here, electrons are initialized at the solid circle markers located at and evolved forward in time including electric field acceleration, collisional drag, and synchrotron radiation. Panel (a) indicates the orbits at , panel (b) for , and panel (c) is for . The dashed black curve indicates the contour , where regions below this curve will have . The parameters were taken to be and , where does not enter since pitch-angle scattering was neglected.
Example deterministic electron orbits. Here, electrons are initialized at the solid circle markers located at and evolved forward in time including electric field acceleration, collisional drag, and synchrotron radiation. Panel (a) indicates the orbits at , panel (b) for , and panel (c) is for . The dashed black curve indicates the contour , where regions below this curve will have . The parameters were taken to be and , where does not enter since pitch-angle scattering was neglected.
The primary impact of pitch-angle scattering on the RPF will be to broaden the sharp transition region and shift its location to higher energy. This is illustrated in Fig. 2, where we have computed the RPF from the steady state adjoint equation, but artificially multiplied the pitch-angle scattering term by factors of 0.01 [Fig. 2(a)], 0.1 [Fig. 2(b)], and 1 [Fig. 2(c)] for a plasma with , while keeping all other terms fixed, including the collisional drag. It is evident that for Fig. 2(a) (where the pitch-angle scattering term is multiplied by 0.01), a sharp transition region is present, where the region with closely aligns with the RE orbits shown in Fig. 1(c). However, as the strength of pitch-angle scattering is increased in Figs. 2(b) and 2(c), the transition region broadens, and the location of the contour shifts to higher energy, implying that RE generation will be less efficient as pitch-angle scattering is increased.
Steady state RPF with the pitch-angle scattering term multiplied by a factor of 0.01 [panel (a)], 0.1 [panel (b)], and 1 [panel (c)]. The parameters were taken to be , and . The RPF was evaluated by solving Eq. (5) at steady state using a physics informed neural network, whose specific implementation is described in Ref. 19.
Steady state RPF with the pitch-angle scattering term multiplied by a factor of 0.01 [panel (a)], 0.1 [panel (b)], and 1 [panel (c)]. The parameters were taken to be , and . The RPF was evaluated by solving Eq. (5) at steady state using a physics informed neural network, whose specific implementation is described in Ref. 19.
B. Runaway electron density evolution
Panel (a) indicates an adjoint problem with the terminal condition , and Eq. (5) is integrated backward to . The solution to this problem will be identical to a problem with the terminal condition at and integrating backward until .
Panel (a) indicates an adjoint problem with the terminal condition , and Eq. (5) is integrated backward to . The solution to this problem will be identical to a problem with the terminal condition at and integrating backward until .
III. PHYSICS-INFORMED NEURAL NETWORKS
A. Physics-constrained deep learning
The adjoint problem discussed in Sec. II allows the RE density to be evolved in time beginning with an arbitrary initial momentum space distribution of REs. Our aim in this section will be to extend the generalizability of the adjoint framework by utilizing a PINN to evaluate the parametric solution to Eq. (5). In this way, the resulting adjoint-deep learning framework will be able to evaluate the RE density for both an arbitrary initial electron distribution, and a range of plasma conditions. Once trained, the resulting framework will allow for nearly instantaneous projections of RE density incorporating fully kinetic physics.
For the present study, we will consider the limit where no data are available, such that we will rely exclusively on physical constraints for training the PINN. A key component when developing a robust PINN will be to enforce several properties of the solution as hard constraints in the NN architecture itself, rather than as soft constraints in the loss function. This prevents the PINN from making unphysical predictions while simultaneously, and perhaps more importantly, severely restricting the solutions that the optimizer searches across, leading to more robust convergence of the PINN. These hard constraints will be enforced by a customized “physics layer,” detailed in Sec. III B below.
B. Physics layer
IV. DECAY OF A PRIMARY RUNAWAY ELECTRON DISTRIBUTION
As a specific application of the above-mentioned framework, we will treat the decay of an initial primary RE distribution. Here, an initial population of REs will be initialized with a given momentum space distribution and the adjoint-deep learning framework described in Secs. II and III will be used to project the RE density forward in time. We will begin by describing representative solutions to the adjoint of the relativistic Fokker–Planck equation (Sec. IV A) and subsequently describe the decay of a primary RE population for a range of initial conditions (Sec. IV B). Finally, the predictions of the adjoint-deep learning framework will be compared with a traditional RE solver in Sec. IV C.
A. Temporal evolution of the runaway probability function
Using the NN architecture and loss function described in Sec. III B above, example solutions to the adjoint problem are shown in Figs. 5 and 7 with the test and training loss history shown in Fig. 4. Here, we have trained the PINN for electric fields in the range , effective charges , and synchrotron radiation . A much broader range of physics parameters will be considered in Ref. 20. We have also chosen a domain with a minimum energy of , a maximum energy of , and . The energy defining which electrons are considered to be REs is taken to be , which corresponds to an energy of roughly 1 MeV. The loss history of the PINN is shown in Fig. 4. Here, the first 15 000 epochs are performed using the ADAM optimizer, while L-BFGS is used thereafter. After periods of 50 000 epochs of training with L-BFGS, the training point distribution is resampled using the residual based adaptive distribution (RAD) method described in Ref. 41. Such a sampling routine redistributes roughly half of the training points using a uniform random distribution, with the remaining points selected with a probability proportional to the magnitude of the residual. In this way, the entire training domain is sampled, but with regions that have relatively large residuals receiving more training points [the final training point distribution is shown in Fig. 4(b)]. This training point resampling results in spikes in the training loss every 50 000 epochs as evident from Fig. 7(a). Despite these spikes in the training loss of the PDE, the test loss of the PDE monotonically decays to a value of , with the boundary loss term dropping to a lower value, suggesting that the PDE is well satisfied across the range of parameters under consideration.
(a) Training and test loss history. The solid blue curve indicates the training loss of the PDE, the solid red curve indicates the training loss for the boundary condition, and the “×” markers indicate the test loss (the test and training losses are the same for the boundary term). 1 000 000 training and testing points are used. The training points initially obey a Hammersley sequence, with the test distribution being uniform random. The training points are resampled using the residual based adaptive distribution described in Ref. 41 every 50 000 epochs of L-BFGS. (b) Projection of the final training point distribution onto the energy and pitch domain.
(a) Training and test loss history. The solid blue curve indicates the training loss of the PDE, the solid red curve indicates the training loss for the boundary condition, and the “×” markers indicate the test loss (the test and training losses are the same for the boundary term). 1 000 000 training and testing points are used. The training points initially obey a Hammersley sequence, with the test distribution being uniform random. The training points are resampled using the residual based adaptive distribution described in Ref. 41 every 50 000 epochs of L-BFGS. (b) Projection of the final training point distribution onto the energy and pitch domain.
Two example solutions are shown in Figs. 5 and 7 for electric fields below and above the avalanche threshold (the avalanche threshold is for these parameters). Considering the subthreshold case first ( , Fig. 5), the contour (black contour in the top row of Fig. 5) shifts to higher momentum and smaller pitch when moving from the to time slice. Physically, this implies that an electron requires an energy near 5 MeV and a pitch near at to have a better than fifty percent chance of remaining in the RE region or being accelerated through the high momentum boundary, by . This is due to the combination of drag and synchrotron losses exceeding electric field acceleration throughout the majority of momentum space, with only a narrow channel where (below the dashed black contour in the bottom row of Fig. 5). Even inside this narrow channel, the presence of pitch-angle scattering will result in electrons initially inside the region being scattered to a larger pitch-angle where , such that they are slowed down42–45 resulting in the RPF taking on a small value across nearly all of momentum space. More insight into the accuracy of the solution in regions where can be gained by considering the of both the RPF and the residual (see Fig. 6). Here, it is evident that the magnitude of the residual of the RPF is small throughout the momentum space, with the magnitude of the residual normalized to the RPF [Fig. 6(c)] being small across regions with significant values of the RPF. We note that at low energies this ratio can exceed one, however, the RPF has already decayed to a value of the order of at these energies, such that the contribution to the RE density will nearly always be negligible from this region.
Panels (a)–(c) indicate time slices of the RPF evolution where panels (d)–(f) indicate the associated residual for an electric field of and . The bottom row indicates the residual to Eq. (5a) multiplied by , where is defined by Eq. (20), such that the low energy divergence is removed, along with the potentially large residual near . The first column indicates the terminal condition , the second column and the last column . The other parameters are given by and . The solid black contour in the top row indicates , whereas the dashed black curve in the bottom row indicates the location of the contour.
Panels (a)–(c) indicate time slices of the RPF evolution where panels (d)–(f) indicate the associated residual for an electric field of and . The bottom row indicates the residual to Eq. (5a) multiplied by , where is defined by Eq. (20), such that the low energy divergence is removed, along with the potentially large residual near . The first column indicates the terminal condition , the second column and the last column . The other parameters are given by and . The solid black contour in the top row indicates , whereas the dashed black curve in the bottom row indicates the location of the contour.
(a) The RPF shown in Fig. 5(c) plotted on a scale. (b) The magnitude of the residual shown in Fig. 5(f) on a scale. (c) Magnitude of the ratio of the residual and RPF, with the result plotted on a scale.
While the magnitude of the residual in Fig. 5 is small for the three time slices shown, suggesting an accurate solution, the final time slice [ , Fig. 7(c)] suggests a limitation of the PINN approach employed. In particular, the high momentum boundary condition employed is when [see Eq. (9c)]. This is not satisfied at . Specifically, as evident from Fig. 5(f), a narrow region near is present where . However, the RPF is not unity at the high momentum boundary in this narrow channel. The reason for this inconsistency is due to pitch-angle scattering quickly diffusing electrons out of the channel where they are subsequently slowed down. Thus, the magnitude of the RPF will quickly drop below unity for even when . Since the present PINN implements the high momentum boundary condition as a soft constraint, deviation from the prescribed boundary condition in this narrow boundary region only leads to a modest penalty in the boundary loss term. We also note that the region where the RPF deviates from the boundary condition ( , and ) represents a corner in the training regime, it is, thus, likely that it is not directly sampled by the random distribution of training points. While this is indicative of a slight inconsistency in the PINN solution, we do not expect it to strongly impact predictions of the RE density since the PINN correctly captures the value of the RPF throughout the majority of the simulation domain, with only a small inconsistency present in the corner of the simulation domain when near marginality. This will be verified in Sec. IV B, where predictions of the PINN are compared with solutions from a traditional RE solver (see Fig. 13).
Same as Fig. 5 but with . Panels (a)–(c) indicate time slices of the RPF evolution where panels (d)–(f) indicate the associated residual.
Same as Fig. 5 but with . Panels (a)–(c) indicate time slices of the RPF evolution where panels (d)–(f) indicate the associated residual.
Considering now a case above threshold ( , Fig. 7), the channel with is substantially broadened, facilitating the acceleration of electrons. In particular, the location of the contour near shifts to lower momenta as time evolves, indicating more electrons will be accelerated by the electric field. Since the present case is only modestly above marginality, pitch-angle scattering introduces a substantial broadening of the transition region. In addition, the RPF reaches a near steady state condition after ten , such that there is only a modest difference between the RPF at compared to .
B. Temporal evolution of the runaway electron density
1. Time histories of RE density for different initial conditions
Runaway electron density evolution for a range of electric fields and for the initial distribution defined by Eq. (21) for [panel (a)], [panel (b)] and [panel (c)]. Panels (d)–(f) indicate the corresponding initial momentum space distribution. Each case took . The electric fields selected correspond to (red curves), (cyan curves), and (blue curves). The other parameters are given by and . The bottom row indicates the initial electron distributions.
Runaway electron density evolution for a range of electric fields and for the initial distribution defined by Eq. (21) for [panel (a)], [panel (b)] and [panel (c)]. Panels (d)–(f) indicate the corresponding initial momentum space distribution. Each case took . The electric fields selected correspond to (red curves), (cyan curves), and (blue curves). The other parameters are given by and . The bottom row indicates the initial electron distributions.
2. Runaway electron decay rate
The time evolution of the RE density can be used to estimate the decay rate of a primary electron population when below threshold. Using a RE decay rate to characterize RE evolution comes with two significant caveats. The first is that while the RE density exhibits an exponential decay when slightly below threshold in the long time limit,45 its temporal evolution cannot be described by an exponential when the system is well below threshold. As a concrete example, for an electric field , the entire population of electrons will slow down due to collisional drag exceeding electric field acceleration at all energies and pitches. This will result in the primary population of electrons slowing down to the bulk in a finite amount of time, a process that cannot be described by a decaying exponential. This physics is illustrated in Fig. 9(a), where the temporal evolution of the RE density for electric fields in the vicinity of for and is shown. It is evident that for (green curve) an exponential curve cannot be fit to the RE density, as evidenced by the poor fit between the dashed and solid green curves. Furthermore, when at threshold ( , red curves) the exponential fit curve only approximates the latter portion of the time history of the RE density. Hence, a second caveat is that it often takes a substantial amount of time, tens to hundreds of ,45 for the REs to settle into a distribution with a clear exponential decay rate. Further noting that for a plasma with and a Coulomb logarithm of , the relativistic collision time is given by , implying that the RE density often will not achieve an exponential decay for tens to hundreds of milliseconds even for . Considering a case with and [Fig. 9(b)], and, thus, a larger avalanche threshold of , a good exponential fit is not present for over twenty , with only a marginally better fit present for .
The time evolution of the RE density for different electric fields. Panel (a) is for and where as panel (b) is for and . The solid curves are predictions from the adjoint-deep learning framework, whereas the dashed curves are exponential fits. The last five were used in the exponential fits. The other parameters are given by , , , and .
The time evolution of the RE density for different electric fields. Panel (a) is for and where as panel (b) is for and . The solid curves are predictions from the adjoint-deep learning framework, whereas the dashed curves are exponential fits. The last five were used in the exponential fits. The other parameters are given by , , , and .
The time required for the RE distribution to settle into a well defined exponential decay depends sensitively on the initial momentum space distribution. Specifically, while the inferred RE decay rate should be independent of the initial momentum space distribution, assuming the initial electrons are at sufficiently high energy, there will inevitably be a transient period while the RE momentum space distribution settles into its saturated form. To assess this uncertainty, a scan of the initial momentum and pitch of the primary electron distributions is shown in Figs. 10(a) and 10(b), respectively. Here, it is evident that the early time evolution of the RE population varies substantially depending on both the initial momentum and pitch of the primary distribution, with resulting fitted exponentially implying a similar, though not identical, decay rate. The modest deviation in inferred decay rates is due to the finite time of the simulations, such that for initial conditions far from those that characterize the decaying primary distribution, twenty is not sufficient for the REs to settle into their saturated decay rate. A strength of the present approach is that by predicting the specific time history of the RE density, rather than only the decay rate, the accuracy of the exponential fit can be readily verified, and the specific time history can be used if required.
The time evolution of the RE density for initial momentum space distributions centered about different momentum [panel (a)] and pitch [panel (b)]. For panel (a), the blue curve is for , the cyan curve is for , the red curve is for , the green curve is for and the magenta curve is for . The pitch was taken to be for all cases. For panel (b), the blue curve is for , the cyan curve is for , the red curve is for , the green curve is for and the magenta curve is for . The momentum was taken to be for all cases. The other parameters are given by , and , with the width of the momentum and pitch distributions given by and , respectively.
The time evolution of the RE density for initial momentum space distributions centered about different momentum [panel (a)] and pitch [panel (b)]. For panel (a), the blue curve is for , the cyan curve is for , the red curve is for , the green curve is for and the magenta curve is for . The pitch was taken to be for all cases. For panel (b), the blue curve is for , the cyan curve is for , the red curve is for , the green curve is for and the magenta curve is for . The momentum was taken to be for all cases. The other parameters are given by , and , with the width of the momentum and pitch distributions given by and , respectively.
Runaway electron decay rate vs electric field for different effective charges [panel (a)] and synchrotron radiation [panel (b)]. Panel (a) took , whereas panel (b) took . The solid curves indicate predictions from the PINN, whereas the blue “×” markers are benchmark values taken from Ref. 45. Panels (c) and (d) indicate the MAPE for the exponential functions used to estimate the RE decay rates in the first row. The initial momentum space distribution assumed , , , and , respectively. The black dashed curve in panel (a) is Eq. (22) for , , and , which is consistent with a value of .
Runaway electron decay rate vs electric field for different effective charges [panel (a)] and synchrotron radiation [panel (b)]. Panel (a) took , whereas panel (b) took . The solid curves indicate predictions from the PINN, whereas the blue “×” markers are benchmark values taken from Ref. 45. Panels (c) and (d) indicate the MAPE for the exponential functions used to estimate the RE decay rates in the first row. The initial momentum space distribution assumed , , , and , respectively. The black dashed curve in panel (a) is Eq. (22) for , , and , which is consistent with a value of .
Insight into the accuracy of the exponential decay rate near threshold can be gained by comparing with RE decay rates inferred in Ref. 45. Here, the RE solver RAMc46 was used to evolve the primary RE distribution for , leading to well converged exponential decay rates. A comparison between the PINN and results from Ref. 45 are shown by the blue “×” markers in Figs. 11(a) and 12. The solid blue curve in both figures indicates results from a PINN trained with , whereas the dashed blue curve in Fig. 12 being an analogous PINN, but trained over a time period of . Twice as many training points were used in this latter case (2 000 000) resulting in an increase in offline training time. It is evident that both results are in good agreement with the benchmark results, though the PINN with the larger exhibits modestly improved agreement.
Comparison of decay rates predicted by a model with (solid blue curve) and (dashed blue curve) with benchmark results from Ref. 45 (“×” markers). The parameters were taken to be , and , , , and .
Comparison of decay rates predicted by a model with (solid blue curve) and (dashed blue curve) with benchmark results from Ref. 45 (“×” markers). The parameters were taken to be , and , , , and .
C. Verification of adjoint-deep learning framework
As a means of verifying predictions of the adjoint-deep learning framework we will compute several time histories of the RE density using the Monte Carlo code RAMc46 for several different initial momentum space distributions and physical parameters. While a comparison of inferred RE decay rates could readily be performed, due to the ambiguities that often emerge when inferring the RE decay rate from the RE time history (as described in Sec. IV B 2), a more stringent test will be to compare detailed time histories of the RE density. A comparison of several temporal histories of RE density with the Monte Carlo code RAMc46 is illustrated in Fig. 13. Here, we have used the initial momentum space distributions illustrated by the second row of Fig. 13, and the values of and are the same as the RPFs shown in Fig. 5 above. The form of these momentum space distributions were chosen to be consistent with the particle initialization scheme used in RAMc. Good agreement is evident for the three electric fields and three initial conditions considered. For the practically important case of an initial RE distribution approximately aligned with the magnetic field ( ) the adjoint-deep learning predictions are in excellent agreement with the RAMc code. For cases with electrons further from , progressively larger deviations are present, though good overall agreement is evident with the largest deviation appearing for the lowest electric field with electrons initialized with [see Fig. 13(c) (red curves)]. The origin of this deviation is due to the limited accuracy of the PINN's prediction of the RPF in regions where the RPF is nearly vanishing. In particular, the small number of REs evident in the solid red curve in Fig. 13(c) results from integrating the product of the RPF shown in Fig. 5(c) with the initial momentum space distribution shown in Fig. 13(f). Since the RPF is nearly zero in the region where the initial electron distribution is finite (i.e., at and ), this, thus, requires the RPF to be computed very accurately to resolve this near zero value. As evident from Fig. 13(c), the PINN is able to accurately track the initial density decay, but exhibits a quantitative discrepancy once only a trace RE population remains.
Panels (a)–(c) show a comparison of the PINN's prediction of the RE density evolution (solid curves) with the Monte Carlo code RAMc (dashed curves). Panels (d)–(f) indicate the corresponding initial primary distribution. The electric fields selected correspond to (red curves), (cyan curves), and (blue curves). The other parameters are given by and . The bottom row indicates the initial electron distributions.
Panels (a)–(c) show a comparison of the PINN's prediction of the RE density evolution (solid curves) with the Monte Carlo code RAMc (dashed curves). Panels (d)–(f) indicate the corresponding initial primary distribution. The electric fields selected correspond to (red curves), (cyan curves), and (blue curves). The other parameters are given by and . The bottom row indicates the initial electron distributions.
V. CONCLUSION
An adjoint-deep learning framework for evaluating the time evolution of the density moment of the RE distribution was derived. By solving the adjoint to the relativistic Fokker–Planck equation, this allows the RE density to be evolved from an arbitrary initial momentum space distribution for a given set of parameters. By utilizing a PINN to identify the parametric solution to the adjoint problem, this further generalizes the framework to treat a range of plasma conditions. While this framework only evolves the RE density moment, by solving the adjoint to the relativistic Fokker–Planck equation, this enables a fully kinetic treatment of relativistic electron dynamics. As an initial application, the evolution of a primary electron population was evolved from three distinct initial pitch distributions, for three different values of the electric field, yielding good agreement with a traditional RE solver. While the offline training time was substantially longer than a traditional relativistic Fokker–Planck solver, the online prediction time is substantially less, typically only requiring tens to hundreds of milliseconds. Furthermore, the decay rate of primary electrons was estimated as a function of electric field strength, effective charge, and the strength of synchrotron radiation. The estimated decay rate was found to vary nonlinearly with the electric field strength, in sharp contrast to commonly employed analytic theories17 that yield a linear dependence on electric field strength. This stiff electric field dependence implies that the electric field will need to remain near threshold during the RE plateau,47,48 to sustain the RE population. While substantial ambiguities in the use of a RE decay rate for characterizing RE evolution when below threshold were identified, by predicting specific time histories of the RE density, the present approach allows for such limitations to be quantified. In particular, aside from predicting the RE decay rate, the prediction of the RE time history enables the quality of the exponential fit to be readily evaluated, where the residual of the PDE provides further insight into the accuracy of the solution.
Despite these encouraging initial results, several limitations of the adjoint-deep learning framework are apparent. One limitation in the present approach is the use of collisional coefficients that are only valid in the completely screened limit. Weakly ionized impurities are typically present in the plasma during a tokamak disruption and modify the collisional processes between electrons.49 While accounting for these partial screening effects will require modifying the PINN, we anticipate a straightforward extension of the existing PINN and is ongoing work. A second limitation, evident from Fig. 10(c), is the modest accuracy that PINNs are achievable when solving a PDE in comparison to traditional Fokker–Planck solvers.5,6,8,44 In particular, for the case shown in Fig. 13(c) the initial RE population has decayed by several orders of magnitude, such that the predicted number of REs is determined by a region of momentum space where the RPF is nearly zero. Since such regions where yield small contributions to the loss function defined by Eq. (15), the PINN will not be able to accurately resolve the RPF in these regions without special care. A third limitation is that the range of parameters that the PINN was trained across in this paper is fairly narrow [ , , ], and thus, will not be able to describe a broad range of tokamak conditions. This narrow training regime was selected to provide a conceptually clear treatment of the adjoint-deep learning framework, while avoiding the technical complexities that emerge when training across a broad range of parameters. This limitation is addressed in Paper II,20 where an adaptive energy and time domain are introduced that enable a far broader parameter regime to be considered while using a comparable number of training points. In addition, while the present approach allows for a broad range of plasma parameters to be considered, it does not allow for a general time dependent variation of those parameters between and . As an initial application of the present results we anticipate using the saturated decay rate of the RE population for a given set of parameters as an estimate of RE decay in an integrated description of a disruption. Finally, an accurate description of RE evolution when above threshold requires the inclusion of large-angle collisions in order to describe the exponential growth of the RE population. This limitation is addressed in Paper II,20 where a Rosenbluth–Putvinski secondary source term is used to approximate large-angle collisions.
ACKNOWLEDGMENTS
This work was supported by the Department of Energy, Office of Fusion Energy Sciences at the University of Florida under Award Nos. DE-SC0024649 and DE-SC0024634 and at Los Alamos National Laboratory (LANL) under Contract No. 89233218CNA000001. The authors acknowledge the University of Florida Research Computing for providing computational resources that have contributed to the research results reported in this publication. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a Department of Energy Office of Science User Facility using NERSC Award No. FES-ERCAP0028155.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Christopher J. McDevitt: Conceptualization (lead); Data curation (equal); Formal analysis (lead); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal). Jonathan S. Arnaud: Conceptualization (equal); Data curation (lead); Software (equal); Visualization (equal); Writing – review & editing (equal). Xian-Zhu Tang: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.