Data-Driven Acceleration of Thermal Radiation Transfer Calculations with the Dynamic Mode Decomposition and a Sequential Singular Value Decomposition
DData-Driven Acceleration of Thermal RadiationTransfer Calculations with the Dynamic ModeDecomposition and a Sequential Singular ValueDecomposition
Ryan G. McClarren a, ∗ , Terry S. Haut b a University of Notre Dame, Dept. of Aerospace & Mechanical Engineering, 365 FitzpatrickHall, Notre Dame, Indiana, USA b Lawrence Livermore National Laboratory, Livermore, California, USA
1. Introduction
The dynamic mode decomposition (DMD) [1, 2] is a data-driven methodfor understanding the spectral properties of an operator. It relies solely on asequence of vectors generated by an operator and requires no knowledge of theoperator. It has enjoyed considerable acceptance in the computational fluiddynamics community for understanding the properties of flows [3] and for com-paring simulation and experiment [1]. For neutron transport problems it wasintroduced as a technique to estimate time-eigenvalues [4], for creating reducedorder models [5], to understand stability [6], and for accelerating power itera-tions for k -eigenvalue problems [7]. In this work we turn to the problem of ac-celerating discrete ordinates solutions to radiative transfer problems (primarilyx-ray radiative transfer for time-dependent high-energy density physics applica-tions).In radiative transfer problems positive solutions, by which we mean positiveradiation densities, are essential due to the coupling of the radiation transportequation to an equation for the material temperature. Negative densities can ∗ Corresponding author
Email addresses: [email protected] (Ryan G. McClarren), [email protected] (Terry S.Haut)
Preprint submitted to Elsevier September 25, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] S e p ead to negative material temperatures that are both nonphysical and causeinstabilities. Moreover, many numerical methods based on high-order repre-sentations of the solution [8] or different angular discretizations [9] can lead tonegative radiation densities.Methods to remove the negative solutions that can arise from these prob-lems have been presented over the years. The zero-and-rescale fix [10, 11] setsany negative values to zero and scales other unknowns to conserve particles.This method has been shown to be effective, but it does not preserve certainmoments of the transport equation. The consistent set-to-zero method (CSZ)[12] addresses this problem by solving a local nonlinear equation to remove non-linearities. Other attempts to address negative solutions are the exponentialdiscontinuous method [13] and the positive spherical harmonics method [14, 15].All of these methods to remove negative solutions (with the exception of theexponential discontinuous scheme) render the solution of radiation transportequation nonlinear. Positive source iteration (a form of nonlinear Richardsoniteration) can still be used, but can be arbitrarily slow to converge on diffusiveproblems [16]. Nevertheless, the nonlinear nature of the solution techniqueimplies that standard acceleration techniques based on linear problems such asdiffusion synthetic acceleration (DSA) [16] and preconditioned GMRES [17] canno longer be used. There have been attempts to derive acceleration methodsbased on Jacobian-free Newton-Krylov [18] and nonlinear acceleration througha quasi-diffusion approach [11].In this paper we propose to use a simple acceleration based on the dynamicmode decomposition. We use DMD to estimate the slowest decaying error modesin positive source iteration, and then estimate the solution. Because DMD is adata-driven method, it is simple to implement. DMD relies on the computationof a singular value decomposition (SVD) of a data matrix containing the solutionfor the scalar intensity (i.e., the scalar flux) at several iterations. To alleviatethe expense of this decomposition we employ a sequential algorithm to performthe SVD that estimates the SVD using rank-one updates. Additionally, becausewe use a sequential algorithm, we can automatically determine the number of2terations to include in the DMD update. The inclusion of a sequential SVD andthe automatic selection of the number of iterations are the two key improvementsover our preliminary work presented at a recent conference [19].While data-driven methods may seem outside the tessitura of particle trans-port research, it is worth noting that Krylov methods such as GMRES can bethought of as data-driven because they do not require the knowledge of thematrix, rather just the action of the matrix. We also believe that the explosionof data being generated in the computational sciences will be another rowel toinvestigate more of these kinds of methods.This paper is organized as follows. In section 2 we introduce the dynamicmode decomposition and its properties. We then in section 3 discuss the gray,discrete ordinates radiative transfer equations and the discontinuous Galerkindiscretization of those equations using Bernstein polynomials. Section 4 givesthe standard, unaccelerated positive source iteration method, before we presentthe DMD acceleration for that method in Section 5. Section 6 gives numericalresults followed by conclusions and future work in section 7.
2. The dynamic mode decomposition
Here we discuss the properties of the dynamic mode decomposition (DMD)for approximating an operator based on information from the action of theoperator. A thorough treatment of the theory of this decomposition can befound in [1, 2, 20, 3].We consider a sequence vectors y k that are related by the application of anoperator A : y k +1 = Ay k . (1)The vectors y k ∈ R N , A is an operator of size N × N , and k = 0 , . . . , K . Thevectors y k could come from a discretized PDE, experimental measurements,sensor readings, etc. As we will see, knowledge of A is not required; only the y k need to be known. 3o find the DMD decomposition, we append the vectors into a data matrices of size N × K as Y + = | | | y y . . . y K | | | Y − = | | | y y . . . y K − | | | . (2)With the data matrices can write Eq. (1) as Y + = AY − . (3)We then take the thin singular value decomposition (SVD) of Y − to write Y − = U Σ V T , (4)where U is a N × K orthogonal matrix, Σ is a diagonal K × K matrix withnon-negative entries on the diagonal, and V is a K × K orthogonal matrix.The matrix U has columns that form an orthonormal basis for the row space of Y − ⊂ R N . In the case when there are only r < K nonzero singular values, weuse the compact SVD where U is N × r , Σ is r × r , and V is r × K .We substitute the SVD of Y − into Eq. (3), to get Y + = AU Σ V T , (5)and then we use the orthonormality properties of V and U , and the fact that Σis a diagonal matrix with non-zero entries to write˜ A ≡ U T AU = U T Y + V Σ − . (6)The matrix ˜ A is an r × r matrix that is a rank r approximation to A , where r is the number of nonzero singular values in the SVD of Y − . Notice in Eq. (6)that ˜ A can be formed using only the data matrices and no knowledge of A .The dynamic modes of A are determined from the eigenvalues of ˜ A . Thisrequires solving an r × r eigenvalue problem. If ( λ, w ) are eigenvalue/eigenvectorpairs of ˜ A then ϕ = 1 λ Y + V Σ − w (7)4re the r dynamic modes of A . The mode with the largest value of λ is said tobe the dominant mode.One of the properties of DMD is that the dynamic modes found will dependon the modes excited by the data. For instance, if y is an eigenvector of A ,then only one mode will be excited. This property was used previously in time-eigenvalue problems in neutron transport to find eigenmodes important to theevolution of an experiment [4].Before moving on, we point out that despite the fact that DMD is derivedas a linear method, it has been shown that DMD can be applied to nonlinearoperators. In particular DMD will find an approximation to the Koopman op-erator for the nonlinear update [3]. This will allow us to use DMD on nonlinearsolution techniques for the radiative transfer equations.
3. The Gray, Discrete Ordinates Radiative Transfer Equations
We will apply DMD to accelerate the solution to gray, radiative transfer cal-culations using discrete ordinates (S N ). The S N equations of thermal radiativetransfer in slab geometry the high-energy density regime are given by1 c ∂I n ∂t + µ n ∂I n ∂x + σ a ( x, t, T ) I n = 12 σ a acT + 12 Q ( x, t ) , (8a) ∂e∂t = σ a ( x, t, T )( φ − acT ) , (8b) φ ( x, t ) = N (cid:88) n =1 w n I n . (8c)Here x [cm] is the spatial variable, t [ns] is the time variable, w n and µ n arethe weights and abscissas of a quadrature rule over the range ( − , I n ( x, t )[GJ/(cm · s · steradian)] is the specific intensity of radiation in the quadraturedirection n , φ ( x, t ) [GJ/(cm · s)] is the scalar intensity, T ( x, t ) [keV] is the ma-terial temperature, and e ( T ) [GJ] is the internal energy density of the materialrelated to T via a known equation of state. Additionally, c ≈
30 [cm/ns] is thespeed of light, a = 0 . · keV )], σ a ( x, t, T ) [cm − ] is the absorp-tion opacity, and Q ( x, t ) is a known, prescribed source. For quadrature rules weapply Gauss-Legendre quadrature rules of even order.5he boundary conditions for Eq. (8) prescribe an incoming intensity on theboundary: I n (0 , t ) = g n ( t ) µ n > , I n ( X, t ) = h n ( t ) µ n < , (9)where g n and h n are known functions of time and X is the right boundary ofthe problem domain. Initial conditions specify I n ( x,
0) throughout the problem.For time discretization we use the backward Euler method with a lineariza-tion of the nonlinear temperature term. We write the solution at time t = m ∆ t using the superscript m : I mn ( x ) = I ( x, m ∆ t ). The semi-discrete equations are[21] µ n ∂I m +1 n ∂x + σ ∗ I m +1 n = 12 (cid:0) σ m s φ m +1 + σ m a f ac ( T m ) (cid:1) + Q ∗ , (10a) e m +1 = e m + ∆ tσ m a ( φ m +1 − ac ( T m ) ) , (10b)where σ ∗ = σ m a + ( c ∆ t ) − , Q ∗ = Q m +1 / c ∆ t ) − , σ s = (1 − f ) σ m a is theeffective scattering term, and the factor f is defined as f ( x, t, T ) = 11 + βcσ a ∆ t , β = 4 aC v (11)with C v the heat capacity at constant volume for the material. It is also usefulto define a radiation temperature as T r = (cid:112) φ/ ( ac ).The system in (10) is a quasi-steady transport problem to which we applya discontinuous Galerkin finite element method in space using the Bernsteinpolynomials as a basis [11, 22]. The resulting equations are( µ n G + F n + M ∗ ) I m +1 n = 12 (cid:0) M s φ m +1 + M a ac ( T m ) (cid:1) + Q ∗ , (12a) M e ( e m +1 − e m ) = M a ( φ m +1 − ac ( T m ) ) , (12b)where the superscript m denotes a time level, µ n G + F n is the upwinded repre-sentation of the derivative term, M ∗ , M s , and M a are the mass matrices asso-ciated with the σ ∗ , σ s , and σ a terms, respectively. The vectors I mn , φ m , T m , and Q ∗ are vectors that contain the coefficients of the finite element representationsof the intensity, scalar intensity, temperature, and source.6he system in (12) can be advanced in time by solving Eq. (12a) and thenevaluating the material internal energy update in Eq. (12b). However, theaddition of the effective scattering term on the RHS of Eq. (12a) couples all ofthe N quadrature directions together.
4. Positive Source Iteration Method
The matrices on the LHS of Eq. (12a) can be written in block lower-triangularform [17]. Therefore, we can perform the following iterative procedure to find I m +1 n I m +1 n (cid:12)(cid:12) k +1 = ( µ n G + F n + M ∗ ) − (cid:20) (cid:0) M s φ m +1 (cid:12)(cid:12) k + M a ac ( T m ) (cid:1) + Q ∗ (cid:21) . (13)Here we denote the k th iteration of a quantity as ( · ) | k . The application of theinverse of the the lower triangular operator ( µ n G + F n + M ∗ ) is known as atransport sweep: it involves moving information for a particular direction n across the problem domain. Note that if we take the quadrature sum of bothsides of Eq. (13) we get an update in terms of the scalar intensity only: φ m +1 (cid:12)(cid:12) k +1 = D ( µ n G + F n + M ∗ ) − (cid:20) (cid:0) M s φ m +1 (cid:12)(cid:12) k + M a ac ( T m ) (cid:1) + Q ∗ (cid:21) , (14)where D represents the quadrature sum (cid:80) Nn =1 w n operator.The iteration scheme in Eq. (14) can be very slow to converge when f → t → ∞ . In this scenario the discrete equations have no absorption ofradiation, leading to the iterations having a spectral radius approaching unity[16]. It has been shown that the iterations can be accelerated by using a diffu-sion correction, called diffusion synthetic acceleration, and by “wrapping” theiterations in a Krylov solver and preconditioning the solver [17]. Physically, the specific intensity is a phase-space density, and as such itshould be non-negative. Nevertheless, it is known that solutions to discrete or-dinates problems can give negative solutions [11, 23]. This is particularly vexing7n radiative transfer problems because negative intensities can lead to negativetemperatures [9] that can cause issues with evaluating material properties.To address this issue we use the zero-and-rescale fix [11] to impose posi-tivity on the intensities in our calculations. This is a nonlinear method thatduring the transport sweep monitors the solution during the sweep. If one ofthe coefficients is negative, this implies that the finite element representationwill have negative values. Therefore, we zero out any negative coefficients andrescale the other coefficients local to a zone to conserve the total intensity ofthe solution locally. Using transport sweeps with the zero-and-rescale fix is aform of nonlinear Richardson iteration.The addition of this nonlinear fix renders acceleration techniques such asdiffusion synthetic acceleration and preconditioned GMRES impotent as thesetechniques require a linear iterative strategy. Recently, Yee, et al. showedthat this nonlinear fix could be accommodated in a nonlinear quasi-diffusioniteration [11]. Here we will show how DMD can be used to handle this type ofnonlinearity as well.
5. DMD Acceleration
In this section we show how DMD can be applied to source iteration usinga sequential SVD. To begin we write Eq. (14) in the following shorthand: y k +1 = Ay k + b (15)where y k +1 = φ m +1 (cid:12)(cid:12) k +1 , and A = 12 D ( µ n G + F n + M ∗ ) − M s , (16) b = D ( µ n G + F n + M ∗ ) − (cid:104) M a ac T m ) + Q ∗ (cid:105) . (17)By substituting in the converged solutions, we can see that Eq. (15) is aniterative procedure for solving ( I − A ) y = b, (18)8here I is the identity operator. Also, if we subtract successive iterations weget the following relationship for the difference between iterations: y k +1 − y k = A ( y k − y k − ) . (19)It is this relationship that we will use DMD to formulate an approximation to A . We define data matrices to contain the differences between iterations Y + = [ y − y , y − y , . . . , y K +1 − y K ] , (20) Y − = [ y − y , y − y , . . . , y K − y K − ] . (21)These are each N × K matrices, where N is the number of spatial degrees offreedom. As before we define an approximate A as the K × K matrix:˜ A = U T AU = U T Y + V Σ − . (22)We can use ˜ A to construct the operator ( I − ˜ A ) − and use this to approximatethe solution.Using Eq. (18) we can write the difference between the solution and the K thiteration as ( I − A )( y − y K ) = b − ( I − A ) y K = b − y K + ( y K +1 − b )= y K +1 − y K . (23)Next, we define ∆ z as the length K vector that satisfies y − y K = U ∆ z, (24)and substitute this into the LHS of Eq. (23) and left multiply by U T to get( I − ˜ A )∆ z = U T ( y K +1 − y K ) . (25)This is a linear system of size r ≤ K where r is the number of non-singularvalues in the SVD of Y − . We solve this system and approximate the solution as y ≈ y K + U ∆ z . 9his algorithm uses the changes between iterations to estimate the operator A that governs the iterative change. We then, in effect, use this approximatedoperator to extrapolate the solution to convergence. This update requires taking K +1 iterations of source iteration, the computation of an SVD, and the solutionof a small linear system. In the previous algorithm, we needed to compute K + 1 iterations in addi-tion to the SVD. However, choosing how many iterations is not obvious. Addi-tionally, the SVD will require O ( N K ) operations to compute where N is thenumber of spatial degrees of freedom. To address both of these problems we usea sequential SVD generated by rank-one updates.Recently, Choi et al. [24] presented an algorithm for taking the SVD of a datamatrix where the elements in the matrix were generated sequentially. We usethis approach to build up the SVD using successive source iterations and deter-mine, based on the results, when K is large enough. The function defined in [24] incrementalSVD takes as inputs a new column vector, u , a tolerance for lineardependence (cid:15) SVD , a minimum size for a singular value (cid:15) SV , the current singularvalue decomposition U, S, V and the column index of the vector, k . Thus, wewrite a call of the incremental SVD as incrementalSVD ( u, (cid:15) SVD , (cid:15) SV , U, S, V, k ). We specify the algorithm for applying automatic DMD with sequential SVDin Algorithm 1. This algorithm uses the incremental SVD function defined inAlgorithm 2 of [24]. Our automatic DMD acceleration takes successive sourceiterations to build up the data matrices defined in Eqs. (20) and (21). After twosource iterations, we have enough data to start applying the acceleration. Wecompute a new value of the estimated solution based on the approximation ˜ A as in Eq. (25). We continue making approximations until either the maximumnumber of iterations is reached or until we find that the two source iterationsdid not add to the rank of ˜ A . This stopping criteria is used because the dataindicates that further source iterations are not improving the approximation ˜ A .10here are other particularities of the automatic DMD acceleration that wepoint out here. Firstly, we remove small singular values from the SVD of Y − .This is done to remove singular values that are unimportant and could addnumerical noise to the update. Additionally, we do not compute an updatebased on DMD if there are eigenvalues of ˜ A with a magnitude larger than one(c.f. line 18 of the algorithm). This is because these large eigenvalues couldallow the solution to diverge in the update.The time update using automatic DMD acceleration is shown in Algorithm2. When we apply the automatic DMD acceleration to compute a time update,we add J additional source iterations outside the DMD acceleration. This isdone to damp any high-frequency errors introduced by the DMD acceleration.In practice we typically use J = 2 or 3. In Algorithm 2 we check for convergencein the source iterations outside the DMD acceleration step. In practice we alsocheck for convergence in the DMD acceleration function to save on iterations,but this detail is omitted from our listing for clarity in the algorithms.To compute a time step for the radiative transfer solver, the storage require-ments for the radiation variables are • Two angular flux vectors for the previous and current angular flux, • The data matrices, Y + and Y − , each of size N × k where N is the number ofspatial degrees of freedom, and k ≤ K is the number of iterations requiredin the DMD acceleration step.The number of iterations (transport sweeps) required for convergence will bethe sum of the iterations outside the DMD acceleration step and those requiredin the DMD update. For comparison with standard source iteration, we use thenumber of transport sweeps as our metric for efficiency.For the nonlinear zero-and-rescale fix, we apply that nonlinearity during thetransport sweep, i.e., in the application of D ( µ n G + F n + M ∗ ) − . This is notexplicitly called out in Algorithm 2, but will be understood in our results.11 lgorithm 1 Automatic DMD Acceleration[ φ ] = AutomaticDMDAcceleration( A , b , φ | , (cid:15) ) Input:
Sweep operators A and b , initial guess φ | , maximum iterations K ,current residual estimate (cid:15) Output:
Approximate solution φ tmpOld = φ | Y + = [], Y − = [] U = [], Σ = [], V = [] for k = 1 to K+1 do tmp= A · tmpOld+ b , i.e., Perform a sweep φ = tmp if k < K + 1 then Append ∆ k = (tmp − tmpOld) to Y − end if if k > then Append ∆ k = (tmp − tmpOld) to Y + [ U, Σ , V, r ] = incrementalSVD (∆ k , (cid:15) × − , (cid:15) × − , U , S , V , k − end if if k > then Remove Singular values from Σ smaller than (cid:15) × − of the trace of Σ Remove columns of U and V associated with the removed singularvalues ˜ A = U T Y + V Σ − Compute eigenvalues of ˜ A as λ k if max k ( | λ k | ) < then Solve ( I − ˜ A )∆ z = U T ∆ k for ∆ z φ = tmpOld + U ∆ z else Not enough iterations to estimate ˜ A , continue end if end if tmpOld = tmp if k > r + 2 then Exit For Loop end if end for return φ m +1 lgorithm 2 Radiative Transfer Time Step Update with DMD[ I n , φ , e , T ] = RadStep( φ m , I mn , T m , e m , J , K , (cid:15) , (cid:15) ∞ ,. . . ) Input:
Previous solutions φ m , I mn , T m , and e m , material properties, quadra-ture rule, number of extra iterations J and maximum DMD iterations K , (cid:15) and (cid:15) ∞ as the L and L ∞ tolerances Output:
Solutions at time level m + 1: φ , I n , T , and e . Compute b = D ( µ n G + F n + M ∗ ) − (cid:2) M a ac ( T m ) + Q ∗ (cid:3) φ | = φ m while Not Converged do { Apply Source Iteration J times } for j=1 to J do φ | j = A φ | j − + b change = (cid:107) φ | j − φ | j − (cid:107) if change < (cid:15) and (cid:107) φ | j − φ | j − (cid:107) ∞ < (cid:15) ∞ then { The iterations are converged } I n = ( µ n G + F n + M ∗ ) − (cid:104) (cid:16) M s φ | j + M a ac ( T m ) (cid:17) + Q ∗ (cid:105) e = e m + M − M a ( φ − ac ( T m ) ) Compute T by inverting the equation of state at e return I n , φ | j , e , and T . end if end for { Apply DMD Acceleration } φ | = AutomaticDMDAcceleration( A , b , φ | J , change) end while . Numerical Results To demonstrate the effectiveness of our acceleration strategy, we consider astandard, diffusive Marshak wave problem [25, 26, 21]. We use σ = 300 T − , andan equation of state given by e = C v T with C v = 0 . · cm − ); thereis no source in the problem. The initial conditions are T ( x,
0) = 0 .
001 keVand φ ( x,
0) = acT ( x, . The domain has an incoming boundary condition of g n = ac/ x = 0 and no incoming radiation at the right edge of the domain.There exists semi-analytic diffusion solutions for this problem that we cancompare with our method to assure that we are converging to the correct answer.We run the problem with different values of the FEM expansion order, numberof spatial zones, and with a time step size of ∆ t = 0 .
01 ns and S Gauss-Legendre quadrature. In Figure 1 results from the DMD solution with cubicelements of size 0 .
02 cm are compared with the semi-analytic diffusion solution.The source iteration solutions are identical on the scale of the figure to the DMDsolutions and are, therefore, not shown. In the figure we see that the S solutionagrees with the semi-analytic solution except near the wavefront, as has beenpreviously observed in comparisons with the diffusion solution [27, 28, 29].To compare the efficacy of our positive, automatic DMD acceleration withstandard positive source iteration we vary the time step size, number of zones,and order of the FEM expansion. In Figure 3 the average number of iterations,that is the number of transport sweeps, per time step is shown. We note thatfor this problem the positivity fix we utilize is needed as the solution near thewavefront can become negative when the fix is not applied.For this problem we inspect the three dominant dynamic modes in the updateof φ found from the first application of DMD acceleration in the step at t = 1and 10 ns. These dominant modes will be estimates of the slowest decayingerror modes from source iteration. In Figure 2 we plot modes as calculated byEq. (7); we normalize the modes by dividing by the maximum magnitude in themode. Note that these magnitude are only defined up to a factor of ±
1. From14 .0 0.1 0.2 0.3 0.4 0.5x (cm)0.00.20.40.60.81.0 T ( ke V ) MaterialRadiation
Figure 1: Comparison of S solutions obtained with DMD and semi-analytic diffusion solutionfor the Marshak wave problem with at 4 different times. The S solutions used ∆ t = 0 .
01 ns,zone sizes of 0 .
02 cm, and a cubic polynomial basis. .000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225x (cm)−1.00−0.75−0.50−0.250.000.250.500.751.00 n o r m a li z e d m ag n i t u d e t=1 ns t=10 ns Dyn. Mode 1Dyn. Mode 2Dyn. Mode 3
Figure 2: Three most dominant dynamic modes found during the time steps at times t = 1and 10 ns. Figure 2 we see that the dominant modes highlight the wavefront and the heatedregion behind it. This is expected because most of the change in the solutionis occurring at the wavefront. Also, it is here that positivity preservation isneeded. The fact that the three modes have nearly the same shape indicatesthat there are several, similar error modes that are slowly decaying.
For this Marshak wave problem, the time step size is a proxy for the scat-tering ratio, σ s /σ ∗ . Using the definition of these quantities and the materialproperties of this problem, we find the scattering ratio simplifies to σ s σ ∗ = (1 − f ) σ a σ a + c ∆ t (26) ≈ t t + 1 . For ∆ t = 0 . , . , and 0 .
02 ns the corresponding scattering ratios are 0.8889,0.9412, 0.9697, respectively.We solve the Marshak wave problem until a final time of t = 10 ns using avariety of spatial resolutions, time steps, and finite element expansion orders.We report the number of iterations (i.e., transport sweeps) required to solve16he problem to the final time in Figure 3. From the figure we see that theDMD-accelerated solutions require significantly fewer iterations than positivesource iteration, fewer than half as many iterations. The difference betweenthe required number of iterations gets larger as the scattering ratio increases toabout a 40% reduction when ∆ t = 0 .
02. For this problem we do not increasethe time step further as this causes overheating due to the large time step andactually makes the medium behave less diffusive.We also notice that as the order of the finite element expansion increasesthe number of iterations required also increases. This was consistent in thesource iteration and DMD-accelerated results. We also observe a slight increasein the number of iterations required in DMD as the number of zones increasesin the ∆ t = 0 .
005 and 0.01 ns cases. The number of iterations required in the∆ t = 0 .
02 ns case, however, decreases as the number of zones increases above40.
There also exist semi-analytic solutions for radiative transfer in an opticallythin problem, driven by a radiation source where the heat capacity of the ma-terial is proportional to the cube of the temperature. Transport and diffusionsolutions for this problem can be found in [30] and S solutions are given in[31]. We solve this problem to demonstrate that the DMD-accelerated solutionis not slower to converge than positive source iteration in optically thin media,where we would expect that no acceleration is need.For this problem we observe that the number of iterations required is almostidentical between DMD and source iteration. At most we see a 5% decrease inthe number of iterations per time step. We can also use this problem to verifyour method; the solutions obtained with DMD obtain reasonable agreementwith the analytic results with a mesh width 0 . t = 6 . × − ns and an expansion order of 3. A comparison of the numerical The solutions are given for the P equations, but S with Gauss quadrature is equivalentto P in slab geometry. I t e r a t i o n s p e r t i m e s t e p Δ t = (a) ∆ t = 0 .
20 40 60 80 100 120 140 160Number of zones2030405060 I t e r a t i o n s p e r t i m e s t e p Δ t = 0.001 order = 1order = 2order = 3order = 4 (b) ∆ t = 0 .
20 40 60 80 100 120 140 160Number of zones406080100120 I t e r a t i o n s p e r t i m e s t e p Δ t = 0.002order = 1order = 2order = 3order = 4 (c) ∆ t = 0 . t = 10 ns. Solid lines denote positive source iterationresults; dashed lines are for DMD-accelerated calculations. .0 0.5 1.0 1.5 2.0 2.5 3.0 3.5x (cm)0.00.51.01.52.0 R a d . En e r g y D e n s i t y ( G J / c m ) S S analyticS transport (a) Radiation T e m p e r a t u r e ( ke V ) S S analyticS transport (b) Material TemperatureFigure 4: Comparison of the numerical results from DMD-accelerated S N and the S /P and transport analytic solutions at times t = 0 . σ/c, σ/c, . σ/c, and 10 σ/c with σ = 1 cm − . and analytic solutions are shown in Figure 4. The final problem we solve is a radiative transfer problem inspired by ex-periments involving laser-driven shocks [32]. In these experiments a laser pulsestrikes a beryllium (Be) disk that is on the end of a xenon (Xe) filled tube. Thelaser launches a shock wave into the Be disk that eventually breaks out into theXe gas. The state of the system at a given time [33] is used to set up our testproblem. The radiative transfer in these shock experiments is complicated dueto the large thin sources that arise [34, 35, 36, 37]. This problem will test howthe DMD acceleration technique performs on a realistic problem with multiplematerials, large variations in density, and optically thin and thick regions.The density, temperature, and material (either Be or Xe) for the test problemare given in Table 1. This table gives the initial conditions for the temperatureand the intensity in equilibrium (i.e., I = acT / . γ = 5 / γ = 1 .
45 in beryllium as calibrated from experiment [38, 39] to19osition (mm) Density (g/cm ) Temperature (eV) Material0.0000 0.0168 40.1723 Be0.0462 0.1681 11.6676 Be0.1046 0.3429 5.5016 Be0.1080 0.5107 3.6279 Be0.1134 0.7383 1.9934 Be0.1241 0.1829 4.6652 Be0.1300 0.1384 6.2204 Be0.1302 0.7405 16.2775 Xe0.1322 0.0493 74.5469 Xe0.6000 0.0065 14.9126 Xe Table 1: Density, initial temperature, and material as a function of position for the radiatingshock problem. D e n s i t y ( g / c m ) T e m p e r a t u r e ( e V ) I n t e r n a l e n e r g y d e n s i t y ( G J / c m ) A b s o r p t i o n o p a c i t y σ a ( c m − ) Figure 5: The density and initial temperature, internal energy density, and σ a for the radiatingshock problem. C v (cid:20) GJkeV · cm (cid:21) = . . . (27)Additionally, we use an approximate bremsstrahlung opacity [40] as σ a ( T ) (cid:2) cm − (cid:3) = 0 . ρ Z T − , (28)for T in keV, ρ the density in g/cm , and Z is the atomic number of the material,4 for Be and 54 for Xe.In Figure 5 the density and the initial values for the temperature, internalenergy density, and σ a are shown as a function of position. At the Be/Xe inter-face there is a jump in the temperature due to the fact that the two materialshave different heat capacities. Between this interface and the other density max-imum at 0.1134 cm there is a region where the absorption opacity drops. Thisis where most of temperature change due to radiative transfer in this problemwill occur. Due to the stiffness of the problem from the large opacity and smallvalue for the heat capacity, we use the modified linearization from [41] with (cid:96) = 5. This has the effect of reducing the value of f in Eq. (11) and increasingthe scattering.To compare the efficiency of our DMD acceleration we solve the radiativetransfer problem for this shock profile over a time step of 0 .
01 ns. We consider aspatial domain extended from x = 0 to 0 .
25 mm with 500 spatial zones and useorder 3 finite elements. The DMD-accelerated solution required 42 iterationswhile positive source iteration required 827; the acceleration led to a speed upof nearly a factor of 20 times. In other words the accelerated solution couldcomplete 20 time steps for the cost of a single, unaccelerated time step. Thesolution for this problem after 100 time steps (i.e., t = 1 ns), is shown in Figure6.
7. Conclusions and Future Work
We have presented a novel method for accelerating the discrete ordinatessolution of radiative transfer problems that includes nonlinear positivity preser-21 .10 0.11 0.12 0.13 0.14 0.15 0.16 0.17x (mm)0.000.010.020.030.040.050.060.07 T e m p e r a t u r e ( ke V ) MaterialRadiationInitial
Figure 6: Solution from DMD-accelerated S for the radiating shock problem at t = 1 nsalong with the initial condition. vation. The acceleration technique, based on the dynamic mode decompositionand an incremental singular value decomposition, was found to be a factor ofnearly 20 faster than positive source iteration on a radiative shock problem andabove three times faster for a standard Marshak wave problem.There is clearly more research that should be performed on this method. Forinstance, we used a linearization of the radiative transfer equations and did notconsider nonlinear temperature updates or non-gray radiative transfer. UsingDMD as a part of multigroup, nonlinear elimination, as outlined in [42], couldbe a fruitful avenue of investigation. There are other radiative transfer prob-lems that could also benefit from a DMD approach. For instance, the iterativeimplicit Monte Carlo (IIMC) method of Gentile and Yee [43] requires a seriesof iterations similar to source iterations. It is possible that DMD accelerationcould also perform well on that formulation. ACKNOWLEDGMENTS
Lawrence Livermore National Laboratory is operated by Lawrence LivermoreNational Security, LLC, for the U.S. Department of Energy, National Nuclear22ecurity Administration under Contract DE-AC52-07NA27344. This document(LLNL-JRNL-1019986) was prepared as an account of work sponsored by anagency of the U.S. government. Neither the U.S. government nor LawrenceLivermore National Security, LLC, nor any of their employees makes any war-ranty, expressed or implied, or assumes any legal liability or responsibility forthe accuracy, completeness, or usefulness of any information, apparatus, prod-uct, or process disclosed, or represents that its use would not infringe privatelyowned rights. Reference herein to any specific commercial product, process, orservice by trade name, trademark, manufacturer, or otherwise does not neces-sarily constitute or imply its endorsement, recommendation, or favoring by theU.S. government or Lawrence Livermore National Security, LLC. The views andopinions of authors expressed herein do not necessarily state or reflect those ofthe U.S. government or Lawrence Livermore National Security, LLC, and shallnot be used for advertising or product endorsement purposes.
References [1] P. J. Schmid, Dynamic mode decomposition of numerical and experimentaldata, Journal of Fluid Mechanics 656 (2010) 5–28.[2] P. J. Schmid, L. Li, M. P. Juniper, O. Pust, Applications of the dynamicmode decomposition, Theoretical and computational fluid dynamics 25 (1-4) (2010) 249–259.[3] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, J. N. Kutz,On dynamic mode decomposition: Theory and applications, Journal ofComputational Dynamics 1 (2) (2014) 391–421.[4] R. G. McClarren, Calculating time eigenvalues of the neutron transportequation with dynamic mode decomposition, Nuclear Science and Engi-neering (2019) DOI:10.1080/00295639.2018.1565014.[5] Z. K. Hardy, J. E. Morel, C. Ahrens, Dynamic mode decomposition for23ubcritical metal systems, Nuclear Science and Engineering 193 (11) (2019)1173–1185.[6] A. Di Ronco, C. Introini, E. Cervi, S. Lorenzi, Y. S. Jeong, S. B. Seo,I. C. Bang, F. Giacobbo, A. Cammi, Dynamic mode decomposition for thestability analysis of the molten salt fast reactor core, Nuclear Engineeringand Design 362 (2020) 110529.[7] J. A. Roberts, L. Xu, R. Elzohery, M. Abdo, Acceleration of the powermethod with dynamic mode decomposition, Nuclear Science and Engineer-ing 193 (12) (2019) 1371–1378.[8] K. A. Mathews, On the propagation of rays in discrete ordinates, Nuclearscience and engineering 132 (2) (1999) 155–180.[9] R. G. McClarren, J. P. Holloway, T. A. Brunner, On solutions to the P n equations for thermal radiative transfer, J. Comput. Phys. 227 (5) (2008)2864–2885. doi:10.1016/j.jcp.2007.11.027 .[10] S. Hamilton, M. Benzi, J. Warsa, Negative flux fixups in discontinuousfinite element SN transport, in: International Conference on Mathematics,Computational Methods Reactor Physics, Saratoga Springs, New York,2009.[11] B. C. Yee, S. S. Olivier, T. S. Haut, M. Holec, V. Z. Tomov, P. G. Mag-inot, A quadratic programming flux correction method for high-order dgdiscretizations of sn transport, arXiv preprint arXiv:1910.02918 (2019).[12] P. G. Maginot, A nonlinear positive extension of the linear discontinu-ous spatial discretization of the transport equation, Master’s thesis, TexasA&M University (2010).[13] W. F. Walters, T. A. Wareing, An accurate, strictly-positive, nonlinearcharacteristic scheme for the discrete-ordinate equations, Transport Theoryand Statistical Physics 25 (2) (1996) 197–215.2414] C. D. Hauck, R. G. McClarren, Positive P N closures, SIAM Journal onScientific Computing, submitted for publication (2009).[15] M. P. Laiu, C. D. Hauck, R. G. McClarren, D. O. S. J. on, 2016, PositiveFiltered P Moment Closures for Linear Kinetic Equations, SIAM 54 (6)(2016) 3214–3238.[16] M. L. Adams, E. W. Larsen, Fast iterative methods for discrete-ordinatesparticle transport computations, Prog. Nucl. Energy 40 (1) (2001) 3–159.[17] J. S. Warsa, T. A. Wareing, J. E. Morel, Krylov iterative methods andthe degraded effectiveness of diffusion synthetic acceleration for multi-dimensional s-n calculations in problems with material discontinuities,Nucl. Sci. Eng. 147 (3) (2004) 218–248.[18] E. D. Fichtl, J. S. Warsa, J. D. Densmore, The Newton-Krylov MethodApplied to Negative-Flux Fixup in S NTransport Calculations, Nuclearscience and engineering 165 (3) (2017) 331–341.[19] R. G. McClarren, T. S. Haut, Acceleration of source iteration using thedynamic mode decomposition, arXiv preprint arXiv:1812.05241 (2018).[20] P. J. Schmid, Application of the dynamic mode decomposition to experi-mental data, Experiments in fluids 50 (4) (2011) 1123–1130.[21] R. G. McClarren, T. M. Evans, R. B. Lowrie, J. D. Densmore, Semi-implicittime integration for P N thermal radiative transfer, J. Comput. Phys.227 (16) (2008) 7561–7586. doi:http://dx.doi.org/10.1016/j.jcp.2008.04.029 .[22] T. Haut, P. Maginot, V. Tomov, B. Southworth, T. Brunner, T. Bailey,An efficient sweep-based solver for the sn equations on high-order meshes,Nuclear Science and Engineering 193 (7) (2019) 746–759.[23] K. A. Mathews, On the propagation of rays in discrete ordinates,Nucl. Sci. Eng. 132 (1999) 155. 2524] Y. Choi, P. Brown, B. Arrighi, R. Anderson, Space-time reduced ordermodel for large-scale linear dynamical systems with application to boltz-mann transport problems, arXiv preprint arXiv:1910.01260 (2019).[25] A. G. Petschek, R. E. Williamson, J. K. Wooten, Jr., Penetration of ra-diation with constant driving temperature, Tech. Rep. LAMS-2421, LosAlamos Scientific Laboratory (1960).[26] T. K. Lane, R. G. McClarren, New self-similar radiation-hydrodynamics so-lutions in the high-energy density, equilibrium diffusion limit, New Journalof Physics 15 (9) (2013) 095013.[27] J. E. Morel, B. T. Adams, T. Noh, J. M. McGhee, T. M. Evans, T. J. Ur-batsch, Spatial discretizations for self-adjoint forms of the radiative transferequations, Journal of Computational Physics 214 (1) (2006) 12–40.[28] R. G. McClarren, R. B. Lowrie, The effects of slope limiting on asymptotic-preserving numerical methods for hyperbolic conservation laws, Journal ofComputational Physics 227 (23) (2008) 9711–9726.[29] R. P. Smedley-Stevenson, R. G. M. J. o. C. Physics, 2015, Asymptotic dif-fusion limit of cell temperature discretisation schemes for thermal radiationtransport, Elsevier 286 (2015) 214–235.[30] B. Su, G. L. Olson, An analytical benchmark for non-equilibrium radiativetransfer in an isotropically scattering medium, Annals of Nuclear Energy24 (13) (1997) 1035 – 1055. doi:DOI:10.1016/S0306-4549(96)00100-4 .[31] R. G. McClarren, J. P. Holloway, T. A. Brunner, Analytic p1 solutions fortime-dependent, thermal radiative transfer in several geometries, Journal ofQuantitative Spectroscopy and Radiative Transfer 109 (3) (2008) 389–403.[32] J. P. Holloway, D. Bingham, C.-C. Chou, F. Doss, R. P. Drake, B. Fryx-ell, M. Grosskopf, B. Van der Holst, B. K. Mallick, R. McClarren, et al.,Predictive modeling of a radiative shock system, Reliability Engineering &System Safety 96 (9) (2011) 1184–1193.2633] R. G. McClarren, D. Ryu, R. P. Drake, M. Grosskopf, D. Bingham, C.-C.Chou, B. Fryxell, B. Van der Holst, J. P. Holloway, C. C. Kuranz, et al.,A physics informed emulator for laser-driven radiating shock simulations,Reliability Engineering & System Safety 96 (9) (2011) 1194–1207.[34] R. G. McClarren, R. P. Drake, Anti-diffusive radiation flow in the coolinglayer of a radiating shock, J. Quant. Spec. Rad. Transf. 111 (2010) 2095–2105.[35] R. G. McClarren, R. P. Drake, J. E. Morel, J. P. Holloway, Theory ofradiative shocks in the mixed, optically thick-thin case, Phys. Plasmas 17(2010).[36] R. P. Drake, Theory of radiative shocks in optically thick media, Phys. Plas-mas 14 (4) (2007) 043301.[37] A. Holgado, J. Ferguson, R. McClarren, Anti-diffusive-like-behavior insemi-analytic radiative shocks via multigroup sn transport with constantcross sections, High Energy Density Physics 17 (2015) 114–118.[38] H. Stripling, R. McClarren, C. Kuranz, M. Grosskopf, E. Rutter, B. Tor-ralva, A calibration and data assimilation method using the bayesian marsemulator, Annals of Nuclear Energy 52 (2013) 103–112.[39] A. Chakraborty, B. K. Mallick, R. G. Mcclarren, C. C. Kuranz, D. Bing-ham, M. J. Grosskopf, E. M. Rutter, H. F. Stripling, R. P. Drake, Spline-based emulators for radiative shock experiments with measurement error,Journal of the American Statistical Association 108 (502) (2013) 411–428.[40] Y. B. Zel’dovich, Y. P. Raizer, Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena, Dover, Mineola, New York, 2002.[41] R. G. McClarren, T. J. Urbatsch, A modified implicit Monte Carlo methodfor time-dependent radiative transfer with adaptive material coupling,J. Comput. Phys. 228 (16) (2009) 5669–5686. doi:10.1016/j.jcp.2009.04.028 . 2742] T. A. Brunner, T. S. Haut, P. F. Nowak, Nonlinear elimination appliedto radiation diffusion, Nuclear Science and Engineering 0 (0) (2020) 1–13. doi:10.1080/00295639.2020.1747262doi:10.1080/00295639.2020.1747262