Data-driven linear time advance operators for the acceleration of plasma physics simulation

The simulation of magnetized plasmas in fusion devices is complicated by the presence of timescales spanning more than 10 orders of magnitude, that is, from those associated with the microscopic gyration of charged particles about magnetic field lines to those of the macroscopic evolution of the plasma (C.f. Fig. 16 in Ref. 1). Slow timescale dynamics can be controlled by the fast timescale physics, for example, the relaxation of a tokamak plasma to equilibrium partially governed by neoclassical and turbulent transport, such that computational methods for tractably predicting this coupled behavior are required. Here, the accuracy of the simulation is determined not only by the stability of the integration method but also by the fidelity of the equations representing the system dynamics.

Direct numerical simulation of microturbulence by gyrokinetic codes on timescales long enough to predict the resulting plasma equilibrium is prohibited by the excessive computational cost. Contemporary efforts utilizing extreme scale computing rely on nearly full utilization of the Titan supercomputer at ORNL for the investigation of heat-flux width on the divertor of a reactor class (ITER) device.^2–4 Still, these and other studies employ frugal expenditure of the available resources to perform attainable calculations. Reference 5 achieved full-f gyrokinetic simulation up to the confinement time required to ignite a burning plasma over a reduced computational domain spanning only a “wedge” of the toroidal volume, which allows for efficient computation with acceptable loss in the fidelity of the turbulent spectrum.⁶ 

The traditional alternative to direct numerical simulation is the separation and subsequent coupling of multiple timescales. Reference 7 treats the problem of obtaining a self-consistent temperature profile across the core and to the edge of a tokamak as the coupling between gyrokinetic turbulence and equilibrium timescale transport solvers. Reference 8 follows a similar methodology to couple the 3D BOUT++ fluid turbulence code over its own timescale of $∼ 10 − 6$ s to the multi-fluid plasma and kinetic neutrals 2D transport SOLPS-ITER (Scrape-off Layer Plasma Simulator - International Thermonuclear Experimental Reactor) code at timescales of $10 − 3 – 10 − 2$ s. The SOLPS-ITER code suite⁹ has also been independently targeted as a candidate for accelerating the convergence to a steady-state solution for use in design and scenario studies. For ITER simulation, SOLPS-ITER can take months of wall-clock time to calculate the interaction of multiple fluid plasma ion species with a Monte Carlo model for kinetic neutrals. Significant effort has been applied to relaxing restrictions on time step size,¹⁰ for example, approximating the contribution of kinetic effects in hybrid modifications to the fluid portion of the code through enforced consistency with kinetic transport rates and boundary conditions;¹¹ the use of different computational meshes in the collisionless kinetic dominated vacuum and highly collisional plasma dominated regions;¹² and direct kinetic corrections to the simulated plasma domain.¹³ 

Another class of methods attempts to accelerate multiple timescale simulation by parallelizing the time dimension, thereby further exploiting the utilization of large computing resources. A “parareal” algorithm was applied to the SOLPS code,^14,15 which parallelizes the simulation over temporal segments and iterations between a “fine” and “coarse” solver. The parareal method corrects coarse initializations of the physical state with limited parallel-in-time application of the fine solver. It was shown that both the fluid neutral configuration of SOLPS and a reduced-grid size kinetic neutral configuration could be used successfully as the coarse solver. A challenge for implementing this predictor–corrector approach is the sensitivity to the selection of parameters in the fine corrections that facilitate convergence of the coarse solution for optimal computational gain on the order of $10 ×$ (defined as the ratio of serial to parallel compute time).

The works of Refs. 16–20 propose an alternative serial approach that retains the physical fidelity of an inner-integrator (microscopic solver) by leveraging in serial outer computations that bridges gaps in the timescales toward latent dynamics (macroscopic solution). This method is “equation-free” in the sense that a model for the slow timescale behavior of the system is never explicitly prescribed. These efforts enforce “projective integration” that “reduces” the fast microscopic system dynamics down to a slow macroscopic system state whose approximate temporal derivative can be “projected” forward in time. The future state of the slow macroscopic system is then “lifted” to reinitialize the fast microscopic simulation, where the procedure can be continued in cycle. This equation-free projective integration framework has been applied to several problems in plasma physics.

Reference 21 first demonstrated projection of ion acoustic wave dynamics, achieving a runtime gain of $13 ×$ in terms of the total PIC (particle-in-cell) simulation cost over the equation-free implementation. The authors demonstrated a linear scaling on the gain with respect to system size on a 1x–1v domain, but found that this behavior required coarsening of the macroscopic gridscale to allow for a stable time step size under first-order explicit integration. Following this work, Ref. 22 adapted the methodology to the kinetic simulation of plasma slab expansion in a vacuum with a wavelet representation in velocity space to suppress particle noise by a sparse basis. Reference 23 generalized these results for gyrokinetic PIC simulation of neoclassical ion heat transport by using a functional representation of the distribution function and developing a polynomial-based reinitialization (lifting) scheme. This transformation eliminated the excitation of spurious transients due to self-consistent relaxation of the system state when large errors are incurred by even stable projective time steps. These works have shown that robust acceleration of plasma physics simulation with projective integration is challenged by ensuring macroscopic time step stability and consistently mitigating error in the lifting from macroscopic to microscopic scales, even in physical problems without multiple timescales.

In this paper, we continue progress in the deployment of projective integration for plasma physics and contribute a scheme that addresses the inherent limitations of stability and accuracy variability in the reduction and projection operations. To directly reveal the stability limits of the projection operator, we apply the data-driven dynamic mode decomposition (DMD) method.²⁴ DMD yields a linear operator approximating the dynamics embedded in a given discretized dataset—here the output of a plasma physics simulation to be accelerated. Like proper orthogonal decomposition (POD),²⁵ DMD gives a low-dimensional representation of the system, but, in contrast to POD, does not require knowledge of the governing equations. DMD further characterizes the model time-dependent state by the frequencies and growth rates of that reduced basis.²⁶ These features of DMD allow for both the determination of the stability limits of each timescale within the extracted linear operator, and to customize that operator by including only those frequencies of interest, for example, the stable slow timescale modes. Here, we utilize the linear operator produced by the basic DMD method²⁷ in both explicit and implicit time advance schemes. There have been many modifications to the standard DMD procedure that addresses its sensitivity with noise²⁸ and transient effects, allowing for extensions of the method to a variety of practical purposes as summarized by Ref. 29 and references therein. A number of these efforts introduce methods with built-in constraints that optimize the extraction of physically meaningful modes.^30–32 Recent work by Ref. 33, in particular, has further adopted statistical techniques to allow uncertainty quantification of both the regression and dynamics forecast from the data. We also note that several authors have presented algorithmic improvements tailored toward applications on multiscale data,^34–36 which we may utilize in future investigations.

The DMD method has a recent history of successful application to plasma physics. It has been applied to data from the HIT-SI (Helicity Injected Torus with Steady Inductive helicity injection spheromak experiment) experiment to aid in the interpretation of nonlinear resistive magnetohydrodynamics (MHD) simulation by the reduction of NIMROD (Non-Ideal Magnetohydrodynamics with Rotation, Open Discussion) simulation output.³⁷ The authors of Ref. 38 also applied DMD to MHD simulation interpretation, and further compared several DMD variants, to identifying unstable modes measured by diagnostic probes. In fluid dynamics, DMD has been shown to accelerate the convergence of iterative methods for solving Navier–Stokes problems when trained on the background mean flow.³⁹ DMD is also closely related to the techniques of biorthogonal decomposition^40,41 and subspace system identification⁴² thoroughly utilized in understanding the stability of MHD modes in magnetic fluctuation data. Diagnostic measurements from experimental field coils have enabled the identification of 3D structures and the real-time evaluation of linear models for plasma control.^43–48 Here, we primarily combine the ideas of using DMD to both reduce and accelerate computations in order to demonstrate a DMD reduction/projection operator for data-driven projective integration of plasma physics simulation, and investigate its numerical properties. We benchmark the constructed operators on several cases of linear fluid simulation and then proceed with a practical application to nonlinear SOLPS-ITER simulation data with varying degrees of numerical noise.

In terms of the stability and controllable accuracy of a DMD-based reduction/projection operator (i.e., does accuracy improve at a predictable rate as the projection time step is reduced), we find both that the slow timescale explicit time step stability limits are now well-defined and that the accuracy of the projection step is indeed controllable when the operator is selectively applied to time intervals satisfying several online error metrics for integration as detailed in this report. We find that DMD linear operators for time advance are capable of predicting the long-term state of nonlinear SOLPS-ITER simulations at an $8 ×$ projected speedup within a range of errors similar to those acceptable to users of that code, that is, of the order 10%.^14,15 This report demonstrates an extension of equation-free projective integration with a data-driven reduced modeling approach to inform the acceleration of plasma physics simulation. The results presented here lay the groundwork for speeding up predictions of SOLPS-ITER time-dependent simulations with DMD, a technique which can also be leveraged in the feedback control of the tokamak plasma boundary where fast models are necessary. In addition, the successful application of DMD to coupled multi-fluid plasma and kinetic neutrals transport promotes further research on this methodology and the degree of utility it can achieve in gyrokinetic and kinetic projective integration on inherently multiple timescale phenomena.

This paper presents a comprehensive analysis of an operator scheme to advance projective integration of fusion plasma physics problems. The remainder of this work is organized as follows: in Sec. II, we introduce a reduction and projection operator for projective integration based on the dynamic mode decomposition (DMD) approximation. In Sec. III, we benchmark the constructed linear time advance operator using simulation data from 1D diffusion and 1D advection–diffusion fluid equations. Section IV demonstrates the application of this method to 1D profiles of electron density and temperature produced by SOLPS-ITER simulation on the inboard and outboard tokamak divertor targets, which indicates the need to limit projections to data intervals with minimal sensitivity to the degree of validity of the linear approximation. Section V summarizes the work and concludes that data-driven linear time advance operators for projective integration allow for the robust calculation of integration time step stability limits and consistent application to both steady-state linear fluid and nonlinear plasma dynamics.

In this section, we propose a data-driven technique for the extrapolation of latent timescales under projective integration. The equation-free approach to single timescale plasma physics simulation has been demonstrated for a posteriori application of reduction and projecting operations (2) and (3) separately. We introduce the formulation of linear time advance operators constructed via dynamic mode decomposition (DMD) to combine these two stages. We will show that the conceived operators provide an accurate least squares approximation of microscopic simulation data. By adopting the DMD framework for reduced modeling, projective integration can be informed by an a priori analysis of the extracted operator's macroscopic time step stability and applicability over nonlinear dynamics. Ensuring the viability of these characteristics for robust and consistent implementation is a critical challenge to the successful acceleration of codes for multiple timescale simulation.

Dynamic mode decomposition (DMD) extracts low-dimensional linear dynamics through matrix factorization of high-dimensional experimental or simulated data. Through the singular value decomposition, DMD identifies spatiotemporal features that are coherent in rank given the covariance structure of the data and separable in time given the eigendecomposition. We follow Algorithm 2, for “exact” DMD, as derived by Ref. 27 in order to construct a reduced model for projective integration. From the obtained set of DMD modes, $ϕ$⁠, and eigenvalues, λ, of an m × n data matrix (spatial coordinates vs time steps), we represent a linear solution of the dynamics up to rank r as

where the temporal coefficients $ω j = log ( λ ) / δ t$ are scaled by the discrete time step and the spatial mode amplitudes b_j are determined by setting an initial condition of the form $b = Φ + x ( 0 )$⁠. Equation (1) solves the dynamic system given by

This linear ODE is completely determined by the matrix exponential, such that for an initial value problem of $x ̃ ( t 0 ) = x ̃ 0$⁠, we write $x ̃ ( t ) = x ̃ ( t 0 ) exp ( T t )$⁠. In this manner, the discretized simulation of a physical system can be re-expressed on the best-fit least squares projected basis $x ̃ 0 = b T Φ + x ( t 0 )$ of the linear approximation to the dynamic timescales such that $x ( t ) = Φ b x ̃ ( t )$ and

where the operator $T adv$ advances each DMD mode independently according to a complex value $ω j = 2 π ν j + i 2 π f j$ with respect to multistep integration methods. The terms, ν and f, specify the growth/decay rates and frequencies of the solution basis, respectively. Thus, not only can unstable modes be characterized and eliminated from consideration but also persistent mode timescales can be identified directly and separated into fast and slow regimes.

The constructed linear time advance operator is subject to the exact same requirements for stable integration as standard finite-difference time domain discretizations of a partial differential equation (PDE). However, the operator is reduced in rank to lower dimensionality as well as separated into an orthogonal spectrum of spatial modes and corresponding dominant timescales. These features enable a robust application to physical systems in which a linear approximation is sufficient to capture the dynamics and improved performance over the computational cost of traditional techniques. Our aim is to demonstrate that the analysis of the extracted operator's properties can inform when such an application is consistently appropriate on a case of scrape-off layer plasma dynamics.

For a one-dimensional linear ordinary differential equation of the form

where u(x, t) is the time-varying solution and $L x$ is the linear operator defining its evolution, we may express the continuous problem in discrete form to solve the corresponding system of algebraic equations

Taking $L x → T adv$ for a reduced rank projected basis over $ϕ ∈ [ 1 , r ]$ instead of $x ∈ ( 0 , M δ x ]$⁠, we may apply the linear multistep integration methods listed in Table I. The first-order explicit forward Euler and implicit backward Euler schemes have stability constraints dependent on the eigenvalues, ω, of the diagonal operator $T adv$ according to

and

respectively, where $h = Δ t$ is an unspecified time step size. The matrix exponential scheme is unconditionally stable and exact, but provides a geometric progression on $T adv$ that will converge to the dominant eigenvalue over a large number of iterations, $N Δ t$⁠.

From Eqs. (6) and (7), we can implement any method(s) and calculate the maximum stable time step size permitted by the scheme. Note that the reduced model Eq. (2) differs from the time advance performed under traditional projective integration in that the latter approximates the macroscopic slow evolution via a first-order explicit extrapolation given by $U ( t + Δ t ) = U ( t ) + Δ t U ′ ( t )$⁠, where $U ′ ( t )$ is the finite difference derivative and is fixed to a discrete single timescale over the simulated spatial discretization x. In our approach, the linear operator $T adv$ is readily extendable to implicit integration and is separable into several timescales of ω over $ϕ$⁠. DMD provides both a reduction of the dynamics in terms of the singular value decomposition (SVD) rank and a projection of independent modes in terms of the exact matrix exponential or other high-order schemes. Sections III and IV will demonstrate the use of DMD as a best-fit least squares reduced model of complete linear dynamics, $u ( t )$⁠, under standard multistep time advance and as a projection operator applied to nonlinear simulation yielding insight into the practical integration, time step, and extrapolation error achievable.

The proposed data-driven extraction of linear time advance operators suitable for reduced-order modeling (within the sampled domain) and projective integration (at time step sizes larger than or outside of the sampled domain) relies on the tuning of several hyperparameters for successful application. Before benchmarking this technique, we outline the methodology for verifying the validity of our approach.

• Interval Size: the number of simulation time steps, N, used to construct an operator. Though the spatial discretization of the system state is fixed by the discretization scheme, the number of snapshots available in a time interval for application of DMD is subject to familiar restrictions. The duration of data $( N δ t )$ is shown to control the slowest frequencies captured by the reduced-order model in both the multi-resolution³⁴ and multi-scale³⁶ variants of DMD. This factor is comparable to the case of the discrete Fourier transform, where oscillatory components $exp ( i 2 π f j t )$ extracted by DMD as modes have amplitudes weighted by growing or diminishing rates $exp ( 2 π ν j t )$⁠. For our purposes, we investigated the integration properties of operators constructed over several fixed intervals of linear fluid dynamics as well as nonlinear plasma dynamics and choose to present a selection of the results for a single size.

• Rank Truncation: the number of singular value components, r, used in the DMD algorithm. For dynamical systems where there is a structured covariance between spatial measurements, the SVD often produces a singular value spectrum with distinct breaks in magnitude. This permits truncation of the orthonormal basis to a lower rank than the original high-dimensional data with an estimate of the least-squares error.^49,50 We specify that the minimum number of singular values kept for operator construction maintain at least 99% of the total energy over the discarded tail of the distribution and compare this truncation against the integration properties of an operator derived from the full SVD spectrum.

• Projection Limit: the number of integration time steps, L, taken by the multistep method. The eigendecomposition of the discrete flow map produces a distribution of frequencies and growth rates scaled by the time step size of the simulation. Due to the Nyquist–Shannon sampling theorem, the maximum fast frequency content that can be resolved by DMD is on the order of $f < 1 / ( 2 δ t )$ with respect to the simulation time step. Analyzing the multistep constraints of the constructed linear operator allows for larger stable $Δ t$ time steps of the slow frequency content than feasible for finite-difference time domain discretizations limited by the Courant–Friedrichs–Lewy condition. The accuracy of the time advance application can be improved by decreasing the time step size, increasing the order of the integration method, and limiting the amount of projection outside of the sampled data set. In this work, we show the full reconstruction of linear fluid dynamics and up to $8 ×$ projected speedup of nonlinear plasma dynamics from brief observations of the simulated dynamics.

• Ensemble Metrics: relative error, ε, from simulation snapshots or model predictions for a sequence of adjacent operator intervals. For nonlinear or complex dynamics where a linear approximation is insufficient to globally predict the solution forward in time, we develop criteria to identify when DMD is locally applicable. We utilize a “committee” approach where a running ensemble of O offset intervals in the set $E = { X | X = { x j , x j + 1 , … , x j + n − 1 } , and 1 ≤ j ≤ O }$ are processed in series and polled on the fidelity of the extracted collection of constructed linear operators. For each fixed interval of size N time steps, the interval reconstruction error and termination time step error within and leaving the sampled data can be calculated. In addition, the deviation of projected values at a fixed future time between each operator in the ensemble yields a measure of model consensus. When all three error metrics satisfy independently established ad hoc thresholds, the criteria for operator projection are met and we elect to project from the leading data interval in the ensemble. This is a simpler alternative to the Gaussian weighted combination of successive DMD models proposed by Ref. 36.

1. Speedup factor: Over a data interval of brief duration with subset size m × N, we construct a linear time advance operator for projection up to an imposed limit of L on the maximum number of integrated time steps that can be taken. For the cases of fluid dynamics presented here, we let N = 128 at the respective time step size necessary to produce simulation data from $t = 1 − N δ t$ and allow L = 1024, which captures the entire available time domain to determine the fidelity of the operators on linear data. For the cases of plasma dynamics, which include nonlinear features, we also consider N = 128 out of 3000 available time steps in the moving interval $t = ( 1 − N δ t ) + t k$ and choose to define the projection limit as

   $L limit = N S speedup$

    

   (8)

   such that the dominant persisting eigenvalue timescales are greater than this limit and the selected operators produce minimal error, ε, as defined below. We remark that this estimate is not normalized by the relative cost in terms of computation time of the preceding simulation over $N δ t$⁠, but it reflects the replacement factor, S, of the simulation by the reduced model.

2. Reconstructed and terminating error: The time-dependent relative error over m spatial observations between the extracted reduced models and simulation is calculated in this work via the discrete Euclidean (L₂) vector two-norm as

   $ε 2 ( t ) = ∑ i = 1 m ( u ̃ ( x i , t ) − u ( x i , t ) ) 2 ∑ i = 1 m u ( x i , t ) 2 .$

    

   (9)
   We define the reconstruction error over the sampled interval of N time steps as

   $ε R E 2 ( t k + N / 2 ) = ∑ i = 1 m ∑ j = 1 N ( u ̃ ( x i , t k + j δ t ) − u ( x i , t k + j δ t ) ) 2 ∑ i = 1 m ∑ j = 1 N u ( x i , t k + j δ t ) 2 ,$

    

   (10)

   where $t k + N / 2$ marks the median time of the recorded data interval. To obtain a measure of the outgoing error past the sampled data, we define the termination error as the discrepancy between the final snapshot in the sampled interval at $ε T M 2 ( t k + N / 2 ) = ε 2 ( t k + N δ t )$⁠.

3. Projected deviation: We devise a scheme to determine the correspondence between sequentially sampled intervals with projected values at future times when local nonlinearities might otherwise prohibit satisfactory global projection of the extracted linear models. For each N snapshot interval starting from t = t_k, we consider the preceding O = 25 intervals, each offset by a time step decrement, and compute projections up to $t p = t k + Δ t$ at a speedup of $S = Δ t / ( N δ t )$⁠. The deviation from the leading projection from the interval starting at time step t_k is defined as

   $d S ( x , t k + Δ t ) = { ( u ̃ i ( x , t p ) − u ̃ k ( x , t p ) ) 2 u ̃ k ( x , t p ) 2 ∈ i = 1 , … , k } .$

    

   (11)

   We take the $L ∞$-norm to represent the degree of agreement among predicted values within the entire operator ensemble, $ε P D , S ( x , t k + N / 2 ) = L ∞ ( d S ( x , t k + Δ t ) )$⁠, and calculate uncertainty bounds as the maximum and minimum values in the set ${ d S ( x , t p ) ∈ i = 1 , … , k w . r . t . x }$⁠.

Under this methodology, linear time advance operators constructed from brief simulation data are subject to two constraints: (1) that the multistep integration method is stable to the extracted eigenvalues of the reduced rank model, and (2) projections are taken only when the in-sample error metrics are satisfied. We formulate the verification of these properties in Secs. III and IV by benchmarking operator construction on linear simulations of fluid dynamics with ASGarD and trialing projection operators for projective integration on nonlinear simulations of plasma dynamics with SOLPS-ITER.

We perform simulations of 1D fluid dynamics with the ASGarD discontinuous Galerkin (DG) sparse grid solver.⁵¹ ASGarD uses a wavelet basis implementation that scales efficiently for high dimensional problems of interest to plasma physics.^52–54 In general, sparse grid methods have the capability to reduce the overall degrees of freedom of DG discretizations to less than that required with standard finite element (FE) methods from $O ( h − d )$ to $O ( h − 1 | | log 2 h | | d − 1 )$ with respect to the dimensionality, d, and smallest mesh size, h. ASGarD calculates the numerical solution of a partial differential equation on a non-uniform mesh grid defined over a hierarchical basis set. Functional interpolation yields a discretized solution of high-order accuracy at the stability limits of standard finite-difference time-domain integration schemes.

In the simulations presented here (discretization and reduction parameters summarized in Tables II), a high-order Runge–Kutta (4) and (5) integration scheme with variable time step size is used to record the simulation at a specified temporal resolution.^55,56 We benchmark the integration performance of linear time advance operators on both diffusion and advection–diffusion type cases and show that projections of the final solution can be obtained at orders of magnitude improved accuracy and stability. Though the selected equations are linear in nature and should readily find expression in the DMD form as reduced-order models, they pose distinct challenges if attempted via the traditional projective integration method.

The 1D diffusion equation, Eq. (12), is a parabolic partial differential equation that describes the single timescale decay of an initial spatial profile subject to a boundary condition and viscosity, $α > 0$⁠. The solution has an exponential form in time, though small time steps are required for a stable explicit forward Euler method. Under a forward-in-time central-in-space (FTCS) discretization the first-order temporal and second-order spatial finite-difference derivatives yield a Von Neumann stability condition for explicit integration on the order of $α δ t / δ x 2 ≤ 1 / 2$⁠. It is desirable to implement a backward Euler implicit method, which is unconditionally stable for this problem. Much larger time steps can be taken to resolve the long-term stationary state of the solution, at somewhat diminished accuracy

From Eq. (12) with the viscosity set to $α = 0.01$ over a spatial domain of $x ∈ [ 0 , L = 2 ]$⁠, we specify an initial condition of

for $k = π$ and apply a Neumann boundary condition satisfying

The ASGarD solution is obtained at a conservative resolution of $δ x = 0.0417$ and $δ t = 0.01$⁠, indicating a maximum attainable step size of around $δ t max = 0.08$ for the explicit forward Euler method under the standard FTCS approach.

A snapshot matrix is assembled from a subset series of the simulation run for an interval size of $[ 49 × 128 ]$⁠, as shown by the vertical dashed white line of Fig. 1. The linear time advance operator representative of the full solution projected out to the final time at T = 10 is constructed from this brief interval following the DMD algorithm outlined in Sec. II.

A Von Neumann stability analysis of the extracted linear operator is straightforward, amounting to the identification of its complex eigenvalues along the diagonal by construction where $ω δ t = ln ( λ )$⁠. Figure 2 presents the distribution of these eigenvalues weighted by the simulation time step for a reduced rank-1 and full rank operator on the left and right panels, respectively. The stability region for explicit forward Euler integration is shown by the offset outlined gray unit circle. Eigenvalues with scaled positions beyond this region are unstable. Implicit backward Euler is unconditionally stable in quadrants I and III along the negative real axis and everywhere on the complex axis. In quadrants II and IV, the method is unstable within an offset unit circle centered on $δ t ω = 1$ on the positive real axis.

Eigenvalues with negative real components are damped with time, while the complex components contribute to oscillatory behavior in the linear combination of modes making up the DMD approximation of the solution [ $∼ exp ( ω t / δ t )$]. It is evident that rank truncation from a full rank linear time advance operator drastically increases the tolerance to large step size explicit integration. The maximum stable time step for the forward Euler method on the rank-1 operator is up to $Δ t = 20$⁠, orders of magnitude greater than the temporal resolution of the sampled data at $δ t = 0.01$⁠. The single eigenvalue is purely real with a value of $ω 1 = − 9.8657 × 10 − 2$⁠.

The presence of spurious modes in the full rank operator is low in singular value energy and does not contribute greatly to the approximated solution. Due to the resolved dynamics in the small interval size used to construct the operator, these modes have poor coherence for projections outside of the sampled domain. Higher rank features impart an instability to explicit integration except for restrictive time steps of the size $Δ t = 0.002$⁠, an order of magnitude less than the sampled simulation data. For these reasons, we only compare the reconstructed and projected accuracy of reduced rank linear time advance operators against the analytic solution

Figure 3 presents the integration of the rank-1 operator following explicit forward Euler (red), implicit backward Euler (blue), and the matrix exponential (black) methods at $δ t = 0.01$ in the left panel and at variable time step sizes, $Δ t$⁠, up to T = 10.24 in the right panel. As a reference, the high-order ASGarD simulation is also shown in these panels with open black circles. The total reconstructed L₂-norm relative error throughout the simulated dynamics is of the order of $10 − 3$ at $δ t = 0.01$ with the ASGarD baseline falling to a level of $10 − 4$ as shown by the left panel. Both the forward Euler and matrix exponential methods remain at an approximately constant level of accuracy as successive time steps are taken, whereas the backward Euler method increases to about a half order of magnitude poorer approximation. We note that under matrix exponential integration the extracted DMD solution closely approximates the specified viscosity used in the initial value problem, where factorization yields $ω 1 = − 0.0099 π 2$⁠.

Linear time advance operators decouple the spatial discretized derivatives by projecting onto an orthogonal basis of DMD spatial modes and being diagonal in form. It follows that a further test to the degree of error accumulation can be determined by varying only the simulation output time step size by factors for two from $8 δ t$ to $δ t / 2$⁠. Operators are constructed from each dataset over initial intervals spanning $t ∈ [ 0 , 1.28 ]$ and compared with the analytic solution at a fixed final integration time of T = 10.24. The right panel of Fig. 3 indicates that the approximated linear time advance maintains first-order convergence (referenced by the dashed black line) for time step sizes greater than or equal to $Δ t = 2 δ t$⁠, but begin to saturate with temporal error accumulation for the backward Euler method and achieve superconvergence for the forward Euler method at time step sizes of $Δ t = δ t$ (marked by the solid vertical black line). The matrix exponential method is predominantly tolerant to time step size error, but shows a slight increase for large number of iterations due to floating-point error magnified by the evaluation of a geometric power progression on the operator eigenvalues.

We next solve Eq. (12) for the initial profile given by

where $k 1 = π$ and $k 2 = 5 π / 2$ for the Neumann boundary condition of $∂ x u ( x = 0 , t ) = 0$⁠. In this case, the analytic solution is given by two exponentially decaying components and has the form

The ASGarD simulation of these dynamics is presented in Fig. 4 at resolution of $δ t = 0.01$ for a domain size of $[ 49 × 1024 ]$⁠. It is evident that the higher spatial scale oscillations, at k₂, damp more quickly than the lower spatial scale oscillations, at k₁. The construction of a linear time advance operator is limited to the initial interval from $t = 0 − 1.28$ as a test of the successful extraction of these features.

The application of the DMD algorithm over this interval size yields an operator with rank greater than the number of features in order to capture the dynamics with sufficient accuracy. Figure 5 shows the distribution of eigenvalues stable to explicit forward Euler integration for the truncated rank-6 and full rank operators in the left and right panels, respectively. For the reduced rank operator, the leading eigenvalues have magnitudes well characterizing the timescale of diffusion at $ω 1 = − 9.8699 × 10 − 2 = − 0.01 k 1 2$ and $ω 2 = − 6.1689 × 10 − 1 = − 0.01 k 2 2$⁠. However, the additional modes have diminished stability properties that contribute to a modest maximum forward Euler stable time step size of $Δ t = 0.028$ around a factor of 3 greater than the simulation at $δ t = 0.001$⁠.

The integration error displayed in the left panel of Fig. 6 follows the ASGarD-simulated solution for the matrix exponential method and is indistinguishable between the less accurate first-order methods through most of the time domain. The time step error shown in the right panel reproduces similar performance to the first case with a single perturbation. The backward Euler method resolves the solution to first order, while the forward Euler method is also first order when stable. Here, the inclusion of additional modes at a fixed higher rank leads to instability in the explicit integration for dynamics resolved at time steps larger than $Δ t = 2 δ t$⁠. Finally, the matrix exponential total relative error is nearly uniform with time step size and on the order of the simulation accuracy.

The 1D advection–diffusion equation (18) is a hyperbolic partial differential equation that describes the multiple timescale transport of an initial spatial profile subject to propagation at constant velocity c, and viscosity $α > 0$⁠. With periodic conditions, the solution is also periodic in time and contains exponential decay. Taking α = 0, the equation is unconditionally unstable to the FTCS discretization scheme. An alternative upwinding method, which is only first-order in both spatial and temporal finite-difference derivatives, is skewed in the direction of the flow and yields a von Neumann stability limit for explicit integration on the order of the Courant–Friedrichs–Lewy (CFL) condition $c ≤ δ x / δ t$⁠. The inclusion of $α > 0$ leads to further restrictive time steps for stable integration, and in general, numerical analysis of 18 is also inhibited by the introduction of artificial dispersion that depends on the finite-difference scheme used and its resolution

In the simple case of pure advection, from Eq. (18) with the viscosity set to α = 0 and a velocity of c = 40 over a spatial domain of $x ∈ [ − 1 , 3 ]$⁠, we specify an initial condition of

where $k = π / 2$ and apply a periodic Dirichlet boundary condition satisfying

The ASGarD solution is obtained at a resolution of $δ x = 0.0816$ and $δ t = 0.0004$⁠. Figure 7 presents the simulation output with respect to an interval size for operator construction that captures only a half cycle of the translating dynamics.

The eigenvalue distributions of the reduced and full rank linear time advance operators are shown in Fig. 8. A set of two dominant DMD modes with complex eigenvalues, $ω 1 , 2 = − 0.0010 ± 62.8321 i$ shown in the left panel, is sufficient to account for the prescribed dynamics. In this case, the maximum forward Euler stable time step size attainable by the operator with explicit integration is $Δ t = 0.032$ s, about two orders of magnitude greater than the sampled data. On the right panel, the higher rank representation clearly produces a collection of spurious and unstable modes.

Figure 9 shows the integration error of the rank-2 operator in the left panel. The matrix exponential preserves the solution at the fidelity of the high-order ASGarD simulation and avoids the large accumulation of numerical dispersion introduced by the first-order methods. Note that the analytic solution is $u ( x , t ) = sin ( k ( x − c t ) )$⁠, which indicates that a significant portion of the total relative error is due to the real part of the extracted eigenvalues contributing minor diffusion at approximately $α = 0.0004$⁠. The right panel of Fig. 8 confirms these results across a range of time step sizes and shows that the matrix exponential integration achieves comparable performance to the ASGarD simulation.

We now consider the case of diffusion-dominated advection for the remaining demonstration of this modeling methodology on linear test problems. Figure 10 presents the ASGarD simulation of Eq. (18) with α = 2 and c = 40 at a time step size of $δ t = 0.0004$⁠. The strength of the diffusive term results in a strongly diminished initial profile that nearly vanishes by the end of the time domain. Again, the interval size for operator construction only captures a half cycle of the periodic translation.

In Fig. 11, the eigenvalue distributions obtained from timescales present in the constructed operators are shown. The right panel indicates that the rank-2 operator retains the same improved forward Euler explicit stability limit of $Δ t = 0.032$ as the pure advection case. These identified timescales have form $ω 1 , 2 = − 4.957 ± i 62.8355$⁠, consistent with the specified system given $Re ( ω 1 , 2 ) = − 2.009 k 2$⁠. Note that the analytic solution in this case has the form $u ( x , t ) = exp ( − α k 2 t ) sin ( k ( x − c t ) )$⁠. The full rank operator in the right panel contains spurious unstable eigenvalues further displaced along the negative real axis than the other cases.

Figure 12 displays the reconstructed integration error in the left panel and highlights that the extracted operator under matrix exponential achieves the same level as accuracy as the high-order ASGarD simulation at a maximum relative error of $10 − 2$⁠. In this case, the first-order schemes are just greater than an order of magnitude increased error over the matrix exponential. The overall accuracy is improved here due to the physical real component of the solution captured by the eigenvalues as opposed to the artificial component introduced in the pure advection case. The right panel of Fig. 12 confirms the rate of convergence with time step size for the respective schemes. When utilized as an exact method, the linear time advance operators demonstrate sufficient reduced modeling of linear diffusion and advection–diffusion fluid test problems that can now be trialed as a projection operator for the integration of nonlinear plasma dynamics.

We perform simulations of the tokamak plasma boundary for dynamics of 1D profiles at the divertor targets with the SOLPS-ITER coupled plasma and neutral transport package.⁵⁷ SOLPS-ITER employs the 2D multi-fluid plasma B2.5 solver and 3D kinetic neutrals EIRENE Markov chain Monte Carlo (MCMC) code. The B2.5 plasma equations for fluid continuity, energy, and momentum of the plasma species are advanced implicitly on a field-aligned grid and account for plasma–neutral/plasma–surface interactions through parallelized MCMC computation of the neutral trajectories. The B2.5 code is not capable of effective parallelization due to the iterative procedure used to calculate quantities in each equation and numerical instabilities preventing convergence.^10,14,15,58 EIRENE, however, utilizes MPI processing on multiple CPU cores for an allotted total time per internal iteration.

In the simulations presented here (discretization and reduction parameters summarized in Table III), we specify a plasma state corresponding to the magnetic equilibrium of DIII-D shot 174310 at 3500 ms with deuterium plasma species only in the fluid ions. We trial the performance of data-driven linear time advance for projections of 1D electron density and temperature profile dynamics and show that operator selection metrics enable a principled characterization of nearly order of magnitude potential speedup with bounded error gained from consistent least squares linear approximation with DMD. Though the plasma physics of the tokamak plasma boundary is inherently nonlinear, robust results can be maintained by taking an ensemble approach to the evaluation of candidate operators derived sequentially in time for projective integration.

Our baseline case resolves 1D profiles of electron density and temperature on the inboard and outboard divertor targets when the input power from a steady-state solution is doubled. The B2.5 simulation is discretized at 38 grid points across the nonuniform scrape-off layer magnetic field geometry and recorded at an internal time step size of $δ t = 10 − 5$ s. Two parameters distinguish the EIRENE iterations of this run from the following case. First, SOLPS-ITER utilizes 16 cores of the ORNL FusionT6 cluster, and second, the maximum time allotted to the CPU for each MCMC calculation is 50 s. A request on the total number of particles to launch can also be set, but is left unspecified for these trial cases. This setup took approximately 8 h of total run time and represents convergence to quasi-steady state.

Figure 13 displays the subset of data produced by the SOLPS-ITER simulation for the relevant target quantities. Linear time advance operators are successively constructed over intervals of I = 128 time steps shifted incrementally by $δ t = 10 − 5$ s across the time domain of 3000 total time steps. In each panel, the initial duration of this interval is indicated by a vertical white line marking $t i = I δ t = 1.3 × 10 − 3$ s. The top left panel shows vertical black lines marking relative projection times from this interval up to a speedup factor (fraction of replaced SOLPS-ITER time steps) of $2 × , 4 ×$⁠, and $8 ×$ for reference.

The top row of Fig. 13 presents the evolution of the inboard and outboard divertor target electron density profile with respect to distance from the separatrix. In the bottom row, the electron temperature is shown for the same regions. For the simulations considered here, asymmetries in the spatial distribution of the quantities as well as gradual temporal approach to steady-state conditions are evident. The linear DMD operators are constructed from finite intervals in time over these nonlinear dynamics and subject to a variable degree of fidelity. To permit consistent analysis, we proceed with a comprehensive examination of the extracted timescales for each variable separately and fix the decomposition rank to r = 12.

Candidate operators are evaluated following the criteria for projection as defined in Eqs. (10), (9), and (11) with respect to the reconstruction accuracy across a sampled interval, the approximation error at the terminating snapshot, and the maximum deviation between projections from an ensemble of adjacent operators to a fixed speedup factor. We impose ad hoc thresholds on each of these metrics, which are required to be satisfied simultaneously in order to select an interval for projection. In this case, the threshold values are set to $τ R E = 0.005 , τ T M = 0.003$⁠, and $τ P D = 0.2$ in terms of relative error. Actual projections are triggered when at least five such consecutive intervals are identified as a means of reducing false-positive detection. The leading operator in an ensemble is used as the starting interval from which the unconditionally stable matrix exponential is taken and integrated out to specified projection time to determine the corresponding error.

Figure 14 showcases the sequence of extracted operator timescales and fidelity of linear time advance operators constructed from an incrementally shifted interval of I = 128 snapshots up to the end of simulation at T = 0.03 s for the inboard electron density. Data are recorded at the central index, rounded down, for each corresponding time interval. The top panel presents the distribution of operator eigenvalue timescales, $Δ t = 1 / | ω i |$⁠, in terms of discrete simulation time steps, N, for $δ t = 10 − 5 s$⁠. Each data point is given a normalized grayscale shading within the operator eigenvalue distribution according to the time-dependent DMD mode amplitude growth/decay, $ψ i = b i exp ( Re ( ω i ) t ) / ψ max$⁠, up to a speedup factor of $8 ×$⁠. In descending order, the horizontal dashed red/blue/black lines indicate a total duration of $N = 8 I$⁠, 4I, $2 I = 1024$⁠, 512, 256 time steps, respectively. These reference annotations are used as a proxy for the projected speedups of $8 × , 4 × , 2 ×$ from the analyzed interval size, indicated by the horizontal solid line. To prevent extrapolation of spurious features, eigenvalues with timescales shorter than the sampled interval are dropped from operator projections past the terminating snapshot.

The overall trend in the top panel of Fig. 14 shows an abrupt period of highly varying timescales in the first 250 time steps of simulation. For the remaining duration, the extracted operators are confined to a narrow band of steadily increasing eigenvalue timescales over three orders of magnitude, from 10³ to 10⁵ time steps. The normalized amplitudes represent the relative contribution of each mode up to the final 8× projection, with a majority of the eigenvalues having timescales above this limit. We track the validity of the local linear DMD approximation in the central panel of Fig. 14. The reconstruction error is plotted in black, and the outgoing error at the terminating snapshot of the operator interval is plotted in gray. The two thresholds placed on these metrics are annotated by the horizontal lines, respectively. Only a few intervals surpass either of these error thresholds and are eliminated as candidates for projection. Most notably in the first series of intervals, the high variability in eigenvalue timescales (shown by the top panel) implies that the interval size is not suited to the rapid change in inboard electron density at the onset of simulation.

Overall, the reconstruction accuracy of the DMD operators is high in the approach to steady-state conditions. The bottom panel of Fig. 14 presents the actual projection error from leading intervals up to a speedup factor of $2 × , 4 ×$⁠, and $8 ×$ in black, red, and blue, respectively. These values have been filtered to only include operators that also satisfy the ensemble deviation metric (shown in the next figure). Between $t = 0.003 − 0.015$ s, both the $2 ×$ and $4 ×$ projections fall nearly an order of magnitude in relative error from 0.1 to 0.01. Though the viable range of $8 ×$ projections has less coverage, the projected error does not exceed a relative error of 0.3 and in the last half of the simulation reaches a mean value of 0.04. We note that the absence of values from the end of the time domain is primarily due to the duration spanned by $2 × , 4 × , 8 ×$ projections.

Figure 15 highlights the uncertainty estimate given by the ensemble deviation metric at two different operator intervals centered at t = 0.0132 and 0.0157 s The top panels show good agreement between operator projections, plotted in blue, and the simulation snapshot, plotted in red, at a speedup factor of $2 × , 4 ×$⁠, and $8 ×$ integration forward in time when the leading interval satisfies all three criteria for projection. The range of uncertainty, indicated by error bars, is computed as the maximum and minimum value at each spatial coordinate among the operator ensemble. In the bottom panel, we show an interval where the ensemble relative deviation exceeds the specified threshold at an $8 ×$ speedup factor. Though this example features a reasonable leading projection, the ensemble contains at least one invalid operator that extrapolates to values with potentially large maximum error. The contaminating timescales appear to have biased growth rates accounting for the several factors of error across the profile outside of the separatrix. Here, it is sufficient to lower the projected duration down to $4 ×$⁠, where the deviation threshold is again satisfied with a high degree of confidence for the same operator ensemble.

The projection characteristics for the inboard divertor electron temperature are similar, but the extracted operators contain a less coherent distribution of extracted timescales as evident in the top panel of Fig. 16. The initial steep rise and fall of eigenvalues is more prominent here than in the electron density over the first 250 time steps, and between 250 and 750 time steps, there is a minority branch of timescales near an 8× speedup factor above 0.3 normalized amplitude. After about 1250 time steps, during which the fixed interval has shifted by 0.0125 s, the operators settle into a narrow band of timescales nearing a discrete level of $4 × 10 4$⁠. In the central panel of Fig. 16, the relative reconstruction errors over the operator intervals are found to be increased by a factor of 3 overall but still fall below the fixed threshold for nearly the entire simulation time domain. The termination error is also offset and crosses its respective threshold frequently. Candidate operator dropouts triggered by the projection criteria are removed from consideration in the bottom panel of Fig. 16. Overall, the relative projection relative error falls from 0.05 (0.1) down to less than 0.01 (0.02) at $2 ×$ (⁠ $4 ×$⁠) speedup. For projections at an $8 ×$ speedup factor, the gap before t = 0.005 s is due to the minority branch of eigenvalues that introduce increased uncertainty while the rest of the selected projections for the remaining time reach a relative error of 0.5 by t = 0.02 s.

Figure 17 selects two different time intervals at T = 0.0125 and 0.0169 s where the projection deviation among the operator ensembles does and does not cross the relevant threshold for the inboard electron temperature. In the top panel of Fig. 17, the projection criteria are satisfied up to an $8 ×$ speedup. The largest deviation among the ensemble rises to $L ∞ = 0.09$ relative to the operator extracted from the leading interval, peaking just below the electron temperature at the separatrix. In the bottom panel of Fig. 17, only the $2 ×$ speedup projected duration is below this threshold. For the $4 ×$ speedup, the ensemble relative deviation reaches $L ∞ = 0.39$ at the separatrix, but the maximum and minimum uncertainty falls to a range of 0.2 and lower further from that point. This is not true for the $8 ×$ speedup projections where extrapolation pushes past the coherence of the extracted modes resulting in a deviation of $L ∞ = 22.56$⁠. Though the leading projection is satisfactory in both cases, the overall lower uncertainty of the $4 ×$ case shows that the ensemble threshold could be reasonably increased to allow for longer realistic projections.

For the outboard profiles on the divertor target, the estimated timescales are highly variable leading to fewer intervals satisfying all projection criteria for $8 ×$ speedup. Figure 18 shows the selection of operators from the electron density simulation. In the top panel of Fig. 18, the extracted eigenvalues are distributed across several orders of magnitude throughout the simulation time domain. For this variable, there is a broadening of the normalized amplitude that is most apparent between n = 250 to 500 time steps as well as continuously around the 8× speedup horizontal annotation. The central panel indicates the fidelity of the extracted operators within the sampled intervals and has an elevated relative error for the reconstructions with a mean relative error of 0.02 at a factor of 2 greater than the inboard case, but only crosses the threshold level once. The termination error exceeds its respective threshold several times in the latter half of the simulation time domain. In the bottom panel of Fig. 18, we see that the relative error of the $2 ×$ speedup projections is confined to within one order of magnitude with the $4 ×$ speedup having minor spread around a relative error of 0.06. In contrast to the inboard simulation, the projection criterion in terms of the ensemble deviation for the 8× speedup is only satisfied at a few short periods owing to the variability in extracted timescales between adjacent intervals.

Figure 19 showcases a selection of intervals from which the ensemble leading operator is projected forward at $2 × , 4 ×$⁠, and $8 ×$ speedup factors. In the top row of Fig. 19, the ensemble deviation falls below the established threshold for each projection. Even at $8 ×$ speedup the simulation snapshot at the projected time is within the ensemble uncertainty at less than $L ∞ = 0.16$ error for the asymmetric feature outside of the separatrix. The operator interval in the bottom row of Fig. 19 only satisfies the projection criterion up to $2 ×$ speedup. At $4 ×$ speedup, the ensemble deviation of $L ∞ = 0.25$ just exceeds the threshold, with the largest uncertainty at the transition between the peak of the separatrix and the adjacent index on positive side. For the $8 ×$ speedup projection, the ensemble maximum on the range of uncertainty greatly exceeds the simulation snapshot in the positive direction where the outboard electron density has significant content. Even the leading projection of the operator ensemble has a peak deviation from the simulated values above the required threshold at around 3 cm outside of the separatrix. This degree of error suggests that the ensemble time interval includes dynamics not adequately extrapolated up to $8 ×$ by the linear operator approximation.

Similar performance is still achieved with operators extracted from simulation of the outboard divertor target electron temperature where the dynamics are more challenging to isolate with linear analysis. Figure 20 demonstrates the high degree of variability of operator timescales in the top panel. The most pronounced features are seen to be multiple clusters of modes below the interval duration, annotated by the solid horizontal black line. Intermediate eigenvalues have reduced normalized amplitudes in the range of 0.2, about half the weight of the corresponding minority distribution in the outboard electron density case. From 1000 timesteps and onward, the dominant operator timescales occur past the $8 ×$ speedup factor up to a discrete value of 10⁴ except for the previously mentioned clustering. In the central panel of Fig. 20, the operator reconstruction relative error has a baseline on the order $10 − 3$⁠, but shows the largest peaks up to an order of magnitude greater than those from among the other simulation variables. These localized discrepancies, and the corresponding termination error peaks, are filtered out from consideration by the projection criteria thresholds. In the bottom panel of Fig. 20, the relative error of projections from intervals satisfying the required criteria is confined below a value of 0.3 for all speedup factors. As shown previously in the analysis of the outboard electron density case, the deviation among operator ensembles restricts $8 ×$ speedup projections to the intervals centered at t = 0.012, 0.017, and 0.02 s.

Figure 21 presents a series of comparisons between projections at several speedup factors and the corresponding simulation snapshots for the outboard electron temperature case. In the top row of panels, the deviation of ensemble projections from the leading interval operator is consistently under the established threshold at or below $L ∞ ≤ 0.2$⁠. Both the $2 ×$ and $4 ×$ projections show minimal disagreement with the recorded simulation snapshot. At the peak electron temperature of the separatrix, the maximum relative error does not exceed a value of 0.05. For the $8 ×$ speedup factor, this peak error is minimal, but the region on outside of the separatrix reaches a relative error of 0.15 well within the estimated ensemble uncertainty range at distances up to 5 cm. The panels in the bottom row consider intervals where the ensemble projection deviation threshold is not satisfied. In contrast to the close approximation of the dynamics by the leading interval operator at a $2 ×$ speedup factor, by $4 ×$ speedup the uncertainty at the peak temperature has increased substantially to about $L ∞ = 0.4$ even though the leading projection is still satisfactory. At an $8 ×$ speedup factor, this interval no longer produces a physically relevant approximation of the target profile due to the decoherence of modes contained in the operator. Here, the greatest deviation among the ensemble from this leading interval exceeds $L ∞ = 25$⁠. As the matrix exponential used in integrating these projections has an exact expression allowing for variable speedup factor, the $L ∞$-norm was chosen more to identify potential outliers within a sampled ensemble duration than to account for a bound on acceptable unknown error of the projection from the simulation. Thus, a conservative threshold of $L ∞ = 0.2$ on the projection deviation was used to avoid issues with greater disagreement found in the rightmost panel than the central panel of the bottom row of Fig. 21.

We estimate the maximum achievable speedup with reasonable relative error on projections of 1D electron density and temperature profiles on the inboard and outboard divertor targets. Extracting linear time advance operators constructed with dynamic mode decomposition and evaluated against a multi-threshold error criterion identifies intervals of simulation data where a projection can be applied. The previous results demonstrate the performance of these metrics and the importance of the operator ensemble on finding a linear approximation to the nonlinear dynamics inherent in SOLPS-ITER calculations that is consistent over a prescribed time duration. Taking all intervals where the criterion is satisfied from a speedup of $2 ×$ at integer multiples up to $8 ×$⁠, we demonstrate that the expected projection error converges to a linear in time trend. Figure 22 compares each selected ensemble's leading operator up to multiple projection factors against the relative error from the respective simulation snapshots. In the left panel of Fig. 22, the inboard projections of electron density, plotted as dotted blue lines, and electron temperature, plotted as dotted red lines, have mean error trends shown by the solid lines with a lower bound of 0.01 and an upper bound of 0.08 relative error for both variables. Overall, operators obtained from the electron temperature simulation have a projected relative error that is a factor of 2 greater than the operators of the electron density simulation. The dashed black line annotating each panel indicates a linear slope on the log –log axis for the range of speedup factors considered here. In the right panel of Fig. 22, there are substantially fewer applicable intervals shown for the outboard divertor target simulation. However, the $8 ×$ speedup factor limiting the criterion corresponds to a duration of 0.01 s, which is approximately one third of the entire simulated time domain in these cases allowing for fewer required projections. The mean trends in the relative error are again seen to perform just as well or better than the identified linear slope.

We analyze a modified case of SOLPS-ITER simulation to confirm the robustness of the extracted linear operators from the presence of numerical noise in nonlinear data. For this run, the B.2 transport equations are resolved over the same spatial and time domain as the previous cases except that the MC computation in EIRENE is only allotted 1 s per iteration and half (8) of the total amount of cores used previously on the ORNL FusionT6 cluster. The resulting dynamics of the 1D profiles of electron density and temperature on the inboard and outboard divertor targets were computed in approximately 3 h of total runtime, but at the expense of introducing deviations in convergence from the nominal solution. To account for the presence of random fluctuations contaminating the error in calculations between the reduced rank operators and this simulation run, the reconstruction and termination thresholds are raised an order of magnitude to $τ R E = 0.03$ and $τ T M = 0.02$⁠, respectively. The ensemble deviation threshold is left at $τ P D = 0.2$ and still yields satisfactory discrimination between projections, as the linear approximation allows for exact matrix exponentiation of all the extracted timescales, up to a rank of 12.

Figure 23 summarizes the expected application of linear time advance operators for projective integration of SOLPS-ITER simulations approaching steady-state conditions and demonstrates that the method is adequately robust to the presence of numerical noise considered here. The mean relative error of projections of electron density and temperature taken on the inboard side across the range of speedup factors from $2 ×$ to $8 ×$ follows a predominantly linear in time trend with a minimum on the order of 0.03 and maximum of 0.08. The overall performance of the operators is seen to be nearly identical to the baseline scenario except for a wider distribution of values. On the outboard side, there are far fewer recorded intervals that satisfy the projection criteria at a speedup of 8× for the mean relative error to be taken with the same level of trust. However, the trend still follows a linear slope with about a factor of 2 greater relative error starting from 0.05 and ending at 0.15. As the speedup factor is increased, less operator coverage within the bounds of expected error is required to perform successful projective integration. For the dynamics analyzed in this work, a speedup of SOLPS-ITER on the order of $8 ×$ and below 0.2 relative error appears to be feasible. Depending on the timescales present in the simulation a larger or variable speedup factor could be applied with reasonable success following the metric threshold methodology developed here.

We demonstrate the application of dynamic mode decomposition (DMD) to the projection stage for projective integration. Linear time advance operators for standard multistep integration schemes are conceived as locally reduced models approximating the global dynamic characteristics of the physical system. We show that for both 1D parabolic and hyperbolic linear fluid equations of diffusion and advection–diffusion, these operators robustly extract stable solutions. Utilizing rank truncation allows for orders of magnitude larger stable time steps than attainable via finite difference time domain approaches at the discretized resolution of these cases. Moreover, the operators maintain first-order and high-order error convergence with respect to time step size. These benchmarks indicate how DMD can be used to extend the projective integration framework to problems where a first-order forward Euler approximation of the time derivative is prohibited by the dynamics.

We further verify the consistent application of linear time advance operators on nonlinear simulations of the tokamak plasma boundary. The SOLPS-ITER coupled multi-fluid plasma and kinetic neutrals transport solver is used to obtain 1D profiles of inboard and outboard divertor target conditions with respect to electron density and temperature. An ensemble of operators constructed from a running series of fixed intervals is down selected to leading projections via three criteria: low reconstruction error, low terminating snapshot error, and low ensemble deviation. These in-sample metrics allow for physically realistic long time step projections up to a speedup factor of $8 ×$⁠. When taking all possible intervals that satisfy the projection criteria, we determine that, for the approach to quasi-steady-state simulations considered here, the model relative error increases to first order with speedup factor up to a maximum of 20% and mean error of 8%. In addition, we obtain the same performance even when this methodology was applied to a simulation run with low MCMC convergence that introduced a higher degree of numerical noise.

Though SOLPS-ITER cannot be fully parallelized, due to the semi-implicit iterative scheme used to evolve the fluid equations, the computation of dominant dynamic modes via DMD over higher dimensional data can be facilitated by efficient algorithms, such as parallelized bootstrap aggregation of the model ensemble, for better performance. The research presented here, therefore, promotes two primary areas of future investigation. First, DMD is shown to achieve satisfactory approximation of SOLPS-ITER plasma dynamics when conservatively applied. Such cases where long-term predictions from limited reduced models that are computationally fast may be useful include the feedback control of the tokamak scrape-off layer, where divertor target conditions must be adequately managed in order to preserve the material integrity of the plasma facing components. Second, the prospective speedup of SOLPS-ITER through projective integration can be used to improve the feasibility of parameter scans in simulations that are computationally expensive due to restrictive time step sizes, such as in the study of impurity seeding for detachment control and the impact of edge-localized modes (ELMs) on tokamak plasma boundary transport. Further development of the discussed computational techniques in this work is required to extend our approach for reduction and projection stages of projective integration with previous efforts toward the capability for accelerated plasma physics simulations of inherently multiple timescale phenomena.

This work was supported in part by the U.S. DOE under Contract No. DE-AC05-00OR22725. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is managed by UT-Battelle, LLC, for the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This research was also sponsored by the Laboratory Directed Research and Development Program of Oak Ridge National Laboratory.

This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy (DOE). The publisher acknowledges the U.S. government license to provide public access under the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

The authors have no conflicts to disclose.

Sebastian De Pascuale: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Writing – original draft (lead). David L. Green: Conceptualization (equal); Data curation (supporting); Funding acquisition (lead); Investigation (supporting); Methodology (supporting); Supervision (supporting). Jeremy D. Lore: Data curation (supporting); Methodology (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.