Semi-implicit Hybrid Discrete \left(\text{H}^T_N\right) Approximation of Thermal Radiative Transfer
mmanuscript No. (will be inserted by the editor)
Semi-implicit Hybrid Discrete (cid:16) H ππ (cid:17) Approximation ofThermal Radiative Transfer
Ryan G. McClarren Β· James A. Rossmanith Β· Minwoo Shin β Received: date / Accepted: date
Abstract
The thermal radiative transfer (TRT) equations form an integro-diο¬erentialsystem that describes the propagation and collisional interactions of photons. Com-puting accurate and eο¬cient numerical solutions to TRT is challenging for severalreasons, the ο¬rst of which is that TRT is deο¬ned on a high-dimensional phase spacethat includes the independent variables of time, space, and velocity. In order to reducethe dimensionality of the phase space, classical approaches such as the P π (sphericalharmonics) or the S π (discrete ordinates) ansatz are often used in the literature. In thiswork, we introduce a novel approach: the hybrid discrete (H ππ ) approximation to theradiative thermal transfer equations. This approach acquires desirable properties ofboth P π and S π , and indeed reduces to each of these approximations in various limits:H π β‘ P π and H π β‘ S π . We prove that H ππ results in a system of hyperbolic partialdiο¬erential equations for all π β₯ π β₯
0. Another challenge in solving the TRTsystem is the inherent stiο¬ness due to the large timescale separation between propa-gation and collisions, especially in the diο¬usive (i.e., highly collisional) regime. Thisstiο¬ness challenge can be partially overcome via implicit time integration, although β Corresponding AuthorRyan G. McClarrenUniversity of Notre DameDepartment of Aerospace and Mechanical EngineeringNotre Dame, IN 46556, USAE-mail: [email protected] A. RossmanithIowa State UniversityDepartment of MathematicsAmes, IA 50011-2104, USAE-mail: [email protected] ShinUniversity of Notre DameDepartment of Aerospace and Mechanical EngineeringNotre Dame, IN 46556, USAE-mail: [email protected] a r X i v : . [ m a t h . NA ] F e b R.G. McClarren, J.A. Rossmanith, and M. Shin fully implicit methods may become computationally expensive due to the strong non-linearity and system size. On the other hand, explicit time-stepping schemes that arenot also asymptotic-preserving in the highly collisional limit require resolving themean-free path between collisions, making such schemes prohibitively expensive. Inthis work we develop a numerical method that is based on a nodal discontinuousGalerkin discretization in space, coupled with a semi-implicit discretization in time.In particular, we make use of a second order explicit Runge-Kutta scheme for thestreaming term and an implicit Euler scheme for the material coupling term. Further-more, in order to solve the material energy equation implicitly after each predictorand corrector step, we linearize the temperature term using a Taylor expansion; thisavoids the need for an iterative procedure, and therefore improves eο¬ciency. In or-der to reduce unphysical oscillation, we apply a slope limiter after each time step.Finally, we conduct several numerical experiments to verify the accuracy, eο¬ciency,and robustness of the H ππ ansatz and the numerical discretizations. Keywords thermal radiative transfer Β· discontinuous Galerkin Β· discrete hybridapproximation Β· semi-implicit time integration Β· asymptotic-preserving Mathematics Subject Classiο¬cation (2010) Β· The thermal radiative transfer (TRT) equations describe the interaction of matter andthermal radiation, and are used in a wide range of applications including remotesensing, glass manufacturing, and combustion [31]. The TRT equations model thepropagation, absorption, and emission of photons, as well as the coupling of thesedynamics to background media via a material-energy equation. Accurate and eο¬cientnumerical solutions of the TRT equations are challenging for several reasons:1. The TRT equations form an integro-diο¬erential system that is posed on a high-dimensional phase space: the time dimension, up to three spatial dimensions, twovelocity direction dimensions, and a frequency dimension;2. Optically thick media introduces stiο¬ness due to the large timescale separationbetween propagation and collisions;3. The material-energy coupling equation introduces strong non-linearity.1.1
Discretization of the phase space
Two broad classes of solution approaches have been used to tackle the TRT equation(1) stochastic methods and (2) deterministic methods; each approach oο¬ers their ownadvantages and disadvantages. The main workhorse of the stochastic approach is theMonte Carlo (MC) method, which is considered one of the most reliable methods inthe radiation community [8,24,36]. Monte Carlo methods produce statistical noisedue to under-sampling of the phase space; and therefore, to compensate for this under-sampling, many MC histories are required, which in turns makes the MC calculationexpensive. As a result, high performance computing approaches are needed. In the emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 3 realm of deterministic approaches the most popular methods are the discrete ordinates(S π ) method [3,13,14,15,16,29] and the spherical harmonics (P π ) method [18,21,25,26,28]. The advantages of these approaches are that, relative to Monte Carlo, theyare less computationally expensive. On the other hand they each suο¬er from their ownspeciο¬c disadvantages.The S π approximation scheme creates a system of equations along a discrete setof angular directions that are taken from a speciο¬c quadrature set [4,10,17,33], andangular integrals are calculated via the given quadrature set. The S π method has beenthe subject of intense research, and many large-scale eο¬cient solution techniques havebeen developed, including approaches that have been shown to scale on leadership-class computers [1]. The S π methodβs main drawback is the phenomenon known asray eο¬ects [14,15], which arises due to the fact that particles move only along certaindirections in the quadrature set. These eο¬ects are conspicuously observed in opticallythin materials with localized sources or sharp material discontinuities.On the other hand, the P π method approximates the solution with a sphericalharmonics expansion that by construction is rotationally invariant and converges inthe πΏ sense [9]. The main disadvantage of the P π approach is that it is plagued byGibbs phenomena (i.e., unphysical oscillations) due to the fact that it is a truncatedspectral method and can thereby produce the negative particle concentrations [9,22],which is physically undesirable. Additionally, in the P π method, when π is small, thewave speeds are reduced because the system eigenvalues can be far from unity.1.2 Stiο¬ness
Another challenging aspect in accurately and eο¬ciently solving the thermal radiativetransfer (TRT) equations is that the system is inherently stiο¬ . This stiο¬ness is especiallypervasive in the diο¬usive regime (i.e., highly collisional), where there is a largetimescale separation between propagation and collisions. The reduction of the kineticTRT system to a lower dimensional system of partial diο¬erential equations (PDEs)via an ansatz such as discrete ordinate (S π ) or spherical harmonics (P π ) does notalleviate this problem.From the perspective of numerical methods, the stiο¬ness challenge can be partiallyovercome via implicit time integration, although fully implicit methods may becomecomputationally expensive due to the strong nonlinearity and system size. On theother hand, explicit time-stepping schemes that are not also asymptotic-preserving inthe highly collisional limit require resolving the mean-free path between collisions,making such schemes prohibitively expensive.An eο¬cient alternative to both fully implicit and fully explicit schemes are semi-implicit schemes, where, at least roughly-speaking, the transport is handled explicitly,while the material coupling term (i.e., collisions) are handled implicitly. For example,Klar et al. [11,12] developed and analyzed an operator splitting approach that wasshown to be asymptotically-preserving while at the same time relatively eο¬cient in allregimes. A diο¬erent semi-implicit approach for TRT was subsequently developed byMcClarren et al. [21], in which a two-stage second-order Runge-Kutta (RK) schemewas used for the streaming term and the backward Euler scheme was used for the R.G. McClarren, J.A. Rossmanith, and M. Shin material coupling term. This scheme was also shown to be asymptotic-preserving inthe diο¬usive limit.1.3
Strong nonlinearity in the material coupling term
In the absence of material coupling, the thermal radiative transfer equations and itsvarious reductions (e.g., P π or S π ) represent a system of linear constant-coeο¬cientpartial diο¬erential equations. The presence of the material coupling term introducesstrong nonlinearity. For explicit time-stepping methods this nonlinearity does not posea direct challenge, although, as mentioned above, explicit methods will then suο¬erfrom small time-steps. For fully implicit time-stepping methods the nonlinearity re-sults in large nonlinear algebraic equations that must be inverted in time-step, resultingin signiο¬cant computational expense. For semi-implicit time-stepping approaches, anonlinear material coupling term only introduces a local nonlinear algebraic equation(i.e., local on each element), although even this can add computational expense. Fortu-nately, McClarren, Evans, and Lowrie [21] showed that the local nonlinear problemspresented by a semi-implicit approach can be linearized via a simple Taylor seriesargument. It was shown that this linearization reduces the computational complexityand does not adversely aο¬ect overall accuracy of the method.1.4 Scope of this work
In order to address some of the shortcomings of classical deterministic methods,we introduce in this work the hybrid discrete (H ππ ) approximation, which hybridizesaspects of the P π and S π methods. The H ππ approximation is equivalent to sphericalharmonics (P π ) approximation when the number of the discretized velocity space( π ) is one and only depends on the order of spherical harmonics basis functions ( π ).Also, it is equivalent to the discrete ordinate (S π ) method for a certain quadratureset. By hybridizing P π and S π , the H ππ approximation is able to acquire beneο¬cialproperties of both of these classical approaches. We prove that H ππ results in a systemof hyperbolic partial diο¬erential equations for all π β₯ π β₯ ππ approximateform, we will discretize the resulting PDE system via the same semi-implicit timeintegration scheme introduced by McClarren et al. [21]. In particular, this approachuses a second order Runge-Kutta explicit time discretization schemes for the streamingterm and a backward Euler scheme for the material energy term; this allows us toresolve the stiο¬ness of the TRT systems and to preserve the asymptotic diο¬usionlimit. Furthermore, in order to solve the material energy equation implicitly aftereach predictor and corrector step, we linearize the temperature term using a Taylorexpansion; this avoids the need for an iterative procedure, and therefore improveseο¬ciency. In order to reduce unphysical oscillation, we apply a slope limiter aftereach time step.In the present work we will consider only grey, i.e., frequency-averaged, radia-tion transport equations in slab geometry. Non-grey and multi-dimensional radiative emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 5 transfer will be studied in future work. For the grey TRT system in slab geometry weconsider several standard test cases in order to validate the accuracy, eο¬ciency, androbustness of scheme, as well as to highlight beneο¬ts over P π and S π solutions.The remainder of this paper is organized as follows. In section 2 the novel hybriddiscrete (H ππ ) approximations in slab geometry are derived ands discussed. In section 3we develop a semi-implicit nodal discontinuous Galerkin ο¬nite element scheme for theresulting H ππ systems. In section 4 we provide numerical results of H ππ approximationson various benchmark problems to show its robustness and asymptotic preservingproperty. Conclusions are presented in section 5. In this section, we formulate H ππ approximations to the radiative transfer in one-dimension. We begin with the thermal radiative transfer in slab geometry, and thenderive the hybrid discrete approximation. We prove that the H ππ system is alwayshyperbolic.2.1 Frequency-independent grey thermal transfer
Consider the 1D scattering-free thermal radiative transfer equation (e.g., see [21]):1 π β (cid:98) πΌ β π‘ + π β (cid:98) πΌ β π§ + (cid:98) π (cid:98) πΌ = (cid:98) π (cid:98) π΅ + (cid:98) π , (2.1)where (cid:98) πΌ ( π‘, π§, π, π ) : R β₯ Γ R Γ [β , ] Γ R β¦β R is the speciο¬c intensity, π‘ β R β₯ is time, π§ β R is the spatial variable, π β [β , ] is the angular variable, i.e., π = cos π with polarangle π β [ , π ] , π β R is the frequency of the photon, (cid:98) π ( π‘, π§, π ) : R β₯ Γ R Γ R β¦β R is anexternal source term, (cid:98) π ( π§, π ) : R Γ R β¦β R β₯ is the opacity, (cid:98) π΅ ( π, π ) : R Γ R β₯ β¦β R β₯ is the Planck function that satisο¬es: π΅ ( π ) : = β« π (cid:98) π΅ ( π, π ) d π = πππ , (2.2)where π ( π‘, π§ ) : R β₯ Γ R β¦β R β₯ is the material temperature in keV, π = Γ cm s β is the speed of light, and π = . Γ ergs cm β keV β is the radiation constant.The notation (cid:98) Β· is used to signify quantities that depend explicitly on the frequency π .Equation (2.1) couples to the material-energy equation: β π β π‘ = β« π β« π (cid:98) π ( π§, π ) (cid:16)(cid:98) πΌ ( π‘, π§, π, π ) β (cid:98) π΅ ( π, π ) (cid:17) d π d π, (2.3)where π ( π, π ) : R β₯ Γ R β₯ β¦β R β₯ is the material energy per volume and π ( π‘, π§ ) : R β₯ Γ R β¦β R β₯ is the material density. R.G. McClarren, J.A. Rossmanith, and M. Shin
For the remainder of this paper we consider only the case where the opacity is fre-quency independent: (cid:98) π ( π§, π ) β π ( π§ ) . In this approximation, we deο¬ne the followingfrequency-integrated quantities (which now removes the (cid:98) Β· notation): πΌ ( π‘, π§, π ) : = β« π (cid:98) πΌ ( π‘, π§, π, π ) d π, π ( π‘, π§ ) : = β« π (cid:98) π ( π‘, π§, π ) d π, (2.4)and equation (2.2). We also deο¬ne the angular moment of the radiation intensity (cid:98) πΌ asfollows: πΈ ( π‘, π§ ) : = π β« π β« π (cid:98) πΌ ( π‘, π§, π, π ) d π d π = π β« π πΌ ( π‘, π§, π ) d π. (2.5)Furthermore, we will assume throughout this work that the material density, π , isconstant, which results in the following: β π β π‘ = β π β π β π β π‘ + β π β π β π β π‘ = πΆ π£ β π β π‘ , (2.6)where πΆ π£ = . Γ erg / cm / keV is the heat capacity at constant volume.Under the assumption of a frequency independent opacity and a constant materialdensity, we can integrate (2.1) over the frequency π and arrive at the following greytransport equation and material-energy equation:1 π β πΌ β π‘ + π β πΌ β π§ + ππΌ = ππππ + π , (2.7) πΆ π£ β π β π‘ = π (cid:18)β« β πΌ ( π‘, π§, π ) d π β πππ (cid:19) . (2.8)respectively. These equations are deο¬ned on π‘ > π§ β ( π§ πΏ , π§ π ) , and π β [β , ] , andmust be equipped with initial conditions at π‘ = π§ = π§ πΏ and π§ = π§ π .2.2 Formulation of the 1D H ππ approximation Equations (2.7)β(2.8) represent a integro-diο¬erential equation for the intensity, πΌ ( π‘, π§, π ) ,and the material temperature, π ( π‘, π§ ) . There are two main classes of techniques forapproximating such systems: (1) direct Monte-Carlo simulations and (2) dimensionreduction via ansatz in π . Standard techniques for reducing the dimensionality ofthese equations to a system deο¬ned only over ( π‘, π§ ) include the P π (i.e., sphericalharmonics) (e.g., see [2]) and S π (i.e., discrete ordinates) [3] methods. Spherical har-monics (P π ) suο¬er from producing negative (and therefore unphysical) intensities,while discrete ordinates (S π ) suο¬er from ray eο¬ects. In this work we consider analternative approach that was ο¬rst developed in Shin [30], namely the H ππ (i.e., hybriddiscrete) approach. The main advantage of H ππ is that it allows us to work with anapproximation that combines desirable aspects of both P π and S π . In particular H ππ reduces to each of these approximations in various limits: H π β‘ P π and H π β‘ S π . emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 7(a) π§ πΏ π§ π π§π β π π π π π β π π β π π π§ πΏ π§ π π§π β π π π π π β π π β π π π§ π§ π§ π β π§ π Fig. 2.1: The H ππ approximation. Panel (a) shows the discrete velocity bands, eachof which is centered at π π β (β , ) for π = , , . . . , π . Panel (b) shows the discretevelocity bands with the physical π§ -mesh superimposed.The H ππ approximation of Shin [30] begins by constructing a mesh in velocityspace. In 1D this means a mesh in the angular variables π β [β , ] : π (cid:216) π = (cid:20) π π β Ξ π , π π + Ξ π (cid:21) , π π = β + (cid:18) π β (cid:19) Ξ π, and Ξ π = π , (2.9)which we depict in Figure 2.1(a). With this velocity mesh, the H ππ approach deο¬neson the velocity band centered at π π the following band-localized intensity: πΌ ( π‘, π§, π ( πΌ )) (cid:12)(cid:12)(cid:12)(cid:12) π β (cid:104) π π β Ξ π ,π π + Ξ π (cid:105) β πΌ ( π ) ( π‘, π§, πΌ ) : = π βοΈ π = (cid:18) π + (cid:19) π’ ( π ) π ( π‘, π§ ) π π ( πΌ ) , (2.10)for π = , , . . . ,π , where πΌ β [β , ] is a local variable in each velocity band centeredat π π with thickness Ξ π : π ( πΌ ) = π π + πΌ (cid:18) Ξ π (cid:19) , (2.11)and π π ( πΌ ) is the π th order Legendre polynomial: π π ( πΌ ) = (cid:18) π β π (cid:19) πΌπ π β ( πΌ ) β (cid:18) π β π (cid:19) π π β ( πΌ ) , π ( πΌ ) = , π ( πΌ ) = πΌ, and β« β π π ( πΌ ) π π ( πΌ ) d πΌ = (cid:18) π + (cid:19) πΏ ππ . (2.12)The moments in (2.10) are deο¬ned as π’ ( π ) π ( π‘, π§ ) = β« β πΌ ( π ) ( π‘, π§, πΌ ) π π ( πΌ ) d πΌ, (2.13) R.G. McClarren, J.A. Rossmanith, and M. Shin for π = , , Β· Β· Β· ,π . Note that in the H ππ approximation π is the number of π -cells,while π is the number of basis functions used in ansatz (2.10).Using the H ππ approximation, on each interval π β (cid:104) π π β Ξ π , π π + Ξ π (cid:105) , equation(2.7) can be written as follows (where we have also temporarily assumed that π β‘ π β πΌ ( π ) β π‘ + (cid:16) π π + πΌ Ξ π (cid:17) β πΌ ( π ) β π§ + ππΌ ( π ) = ππππ , (2.14)for π = , , Β· Β· Β· ,π and πΌ β [β , ] . Plugging ansatz (2.10) into (2.14), multiplyingthe equation by the Legendre polynomial π β , integrating over [β , ] , and using theorthogonality of Legendre polynomials (2.12), gives the following equations:1 π β π’ ( π ) β β π‘ + π π β π’ ( π ) β β π§ + Ξ π π βοΈ π = π βπ β π’ ( π ) π β π§ + ππ’ ( π ) β = ππππ πΏ β , (2.15)for β = , , Β· Β· Β· , π , where via the three-term recurrence relationship and orthogonalityfrom (2.12): π βπ = (cid:18) π + (cid:19) β« β πΌπ π ( πΌ ) π β ( πΌ ) d πΌ = (cid:18) π (cid:19) β« β π π β ( πΌ ) π β ( πΌ ) d πΌ + (cid:18) π + (cid:19) β« β π π + ( πΌ ) π β ( πΌ ) d πΌ =  π π β , if β = π β , π + π + , if β = π + , , otherwise . (2.16)We can write the resulting system of equations as follows:1 π β u ( π ) β π‘ + A ( π ) β u ( π ) β π§ = q (cid:16) u ( π ) , π (cid:17) , (2.17)for π = , , Β· Β· Β· ,π , where u ( π ) = (cid:16) π’ ( π ) π’ ( π ) . . . π’ ( π ) π (cid:17) π β R π + , A ( π ) = π π I + Ξ π (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) . . .. . . . . . π β π β π β π β π π β π π + (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) , (2.18) q (cid:16) u ( π ) , π (cid:17) = β π u ( π ) + ππππ e , (2.19)where e = ( , , . . . , ) π β R ( π + ) , I , A ( π ) β R ( π + )Γ( π + ) , and q (cid:0) u ( π ) , π (cid:1) : R ( π + ) Γ R β₯ β¦β R ( π + ) . Finally, we obtain the linear system:1 π β u β π‘ + A β u β π§ = Q ( u , π ) , (2.20) emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 9 where A = (cid:169)(cid:173)(cid:173)(cid:171) A ( ) . . . A ( π ) (cid:170)(cid:174)(cid:174)(cid:174)(cid:172) , Q ( u , π ) = (cid:169)(cid:173)(cid:173)(cid:171) q (cid:0) u ( ) , π (cid:1) ... q (cid:0) u ( π ) , π (cid:1)(cid:170)(cid:174)(cid:174)(cid:174)(cid:172) , and u = (cid:169)(cid:173)(cid:173)(cid:171) u ( ) ... u ( π ) (cid:170)(cid:174)(cid:174)(cid:174)(cid:172) , (2.21)where A β R ( π + ) π Γ( π + ) π , Q ( u , π ) : R ( π + ) π Γ R β₯ β¦β R ( π + ) π , and u β R ( π + ) π . Proposition 1
Equation (2.20) β (2.21) is a system of hyperbolic partial diο¬erentialequation.Proof To prove hyperbolicity we must show that the matrix A β R ( π + ) π Γ( π + ) π in(2.20)β(2.21) is diagonalizable with only real eigenvalues. Showing this is equivalentto showing that each block matrix A ( π ) β R ( π + )Γ( π + ) is diagonalizable with onlyreal eigenvalues, since A is simply a block diagonal matrix with blocks A ( π ) for π = , , . . . ,π .To show that A ( π ) is diagonalizable with only real eigenvalues, we write thetridiagonal matrix A ( π ) as follows and deο¬ne a diagonal matrix D : A ( π ) = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) π π π π π π π . . .. . . . . . π π β π π β π π π π π π (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) , D = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) βοΈ π π βοΈ π π . . . βοΈ π π π π (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) , (2.22)where π = π π and for all π = , , Β· Β· Β· , π : π π = π /( π β ) β π π = ( π + )/( π + ) β
0. It can then easily be shown that P ( π ) = D β A ( π ) D , where P ( π ) = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) π β π π β π π π β π π β π π π . . .. . . . . . β π π β π π β β π π β π π β π β π π π π β π π π π π (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) . (2.23)Now since the matrix P ( π ) is real and symmetric, it has only real eigenvalues. Also, allthe eigenvalues of P ( π ) are distinct since all oο¬-diagonal elements are nonzero [27].Hence, A ( π ) has only real eigenvalues and all the eigenvalues are distinct by similarityto P ( π ) . (cid:117)(cid:116) Proposition 2
The eigenvalues of the matrix A β R ( π + ) π Γ( π + ) π in (2.20) β (2.21) are real, distinct, and given by π π +( π β ) ( π β ) = π π + (cid:18) Ξ π (cid:19) π π for π = , . . . , π + , π = , . . . ,π, (2.24) where π π β (β , ) for π = , , . . . , π + are the real distinct roots of the degree ( π + ) Legendre polynomial. If we assume that the π + Legendre roots are ordered so that π π + is the largest root (i.e., closest to +1), then the largest eigenvalue (in absolutevalue) of A is given by π ( A ) = β π = π ( π + ) π = π π + (cid:18) Ξ π (cid:19) π π + . (2.25) Proof
It can be shown that the eigenvalues of A ( π ) β R ( π + )Γ( π + ) from (2.17) aregiven by the following formula (e.g., see Cohen [5]): π ( π ) π = π π + (cid:18) Ξ π (cid:19) π π , for π = , , . . . , π + , (2.26)where π π β (β , ) for π = , , . . . , π + ( π + ) Legendrepolynomial. Since A is just a block diagonal matrix with A ( π ) β R ( π + )Γ( π + ) as theblocks for π = , , . . . ,π (as shown in (2.21)), it follows that the eigenvalues of A areof the form: π π +( π β ) ( π + ) = π π + (cid:18) Ξ π (cid:19) π π for π = , . . . , π + , π = , . . . ,π . (2.27)The most negative and most positive eigenvalues in (2.27) are π ( π = π = π ( π + ) π ( π = π + π = π ), respectively. This follows directly from the assumedordering of π π , as well as the ordering assumed in the deο¬nition of π π (see equation(2.9)). In fact, these two eigenvalues have the same magnitude, since by deο¬nition: π = β π π + and π = β π π . This gives the desired result: equation (2.25). Proposition 3 π ( A ) < for any π and π of H ππ .Proof A simple calculation, using the fact that π π + <
1, shows that π ( A ) = π π + (cid:18) Ξ π (cid:19) π π + = β (cid:18) Ξ π (cid:19) + (cid:18) Ξ π (cid:19) π π + < β (cid:18) Ξ π (cid:19) + (cid:18) Ξ π (cid:19) = , = β π ( A ) < . Furthermore, we note that π ( A ) = β O ( Ξ π ) as Ξ π β (cid:117)(cid:116) Remark 2.1
We have shown via the above propositions that all of the eigenvalues of A are real, distinct, and strictly between β A not have a zero eigenvalue; this additional requirementis useful especially in the computation of steady-state solutions and in the impositionof inο¬ow/outο¬ow boundary conditions. We can always ensure that A does not have azero eigenvalue if either of the following conditions are satisο¬ed:1. π is even; or2. π is odd and π is odd.These assertions follow directly from deο¬nition (2.27). emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 11 The radiation energy density
Using the H ππ asantz (2.10), the angular moment of the radiation density (2.5) can bewritten as follows: πΈ ( π‘, π§ ) : = π β« β πΌ ( π‘, π§, π ) d π = π Ξ π π βοΈ π = β« β πΌ ( π ) ( π‘, π§, πΌ ) d πΌ = Ξ π π π βοΈ π = β« β π βοΈ π = (cid:18) π + (cid:19) π’ ( π ) π ( π‘, π§ ) π π ( πΌ ) d πΌ = Ξ π π π βοΈ π = π’ ( π ) ( π‘, π§ ) . (2.28)Using this result, we can rewrite material-energy equation (2.8) as follows: πΆ π£ β π β π‘ = π Ξ π π βοΈ π = π’ ( π ) ( π‘, π§ ) β ππππ . (2.29) In this section we develop a discontinuous Galerkin (DG) ο¬nite element method witha semi-implicit time discretization for solving the 1D H ππ system (2.20)β(2.21) alongwith the material-energy equation (2.29). The strategy that we employ is similar to themethod developed by McClarren et al. [21] for the P π equations, but here we extendthat work to the H ππ approximation as described above.3.1 Discontinuous Galerkin spatial discretization
The H ππ approximation as described in the previous sections divides the phase spaceinto discrete velocity bands, inside of which we approximate the speciο¬c intensity, πΌ ( π‘, π₯, π ) , by a ο¬nite polynomial ansatz; this is illustrated in Figure Figure 2.1(a). Inorder to spatially discretize the resulting H ππ equations: (2.20)β(2.21) and (2.29), weadditionally introduce a mesh in the π§ -coordinate: [ π§ πΏ , π§ π ] = π π§ (cid:216) π = T π , where T π = (cid:20) π§ π β Ξ π§ , π§ π + Ξ π§ (cid:21) , (3.1) π§ π = π§ πΏ + (cid:18) π β (cid:19) Ξ π§, and Ξ π§ = π§ π β π§ πΏ π π§ . (3.2)The full π§ β π phase space mesh is illustrated in Figure Figure 2.1(b). On each meshelement we deο¬ne a local coordinate as follows: π§ (cid:12)(cid:12)(cid:12) T π = π§ π + (cid:18) Ξ π§ (cid:19) π, where π β [β , ] . (3.3) The test and trial functions for the discontinuous Galerkin scheme will come fromfollowing broken ο¬nite element spaces: W Ξ π§π eqn : = (cid:26) w Ξ π§ β (cid:104) πΏ β [ π§ πΏ , π§ π ] (cid:105) π eqn : w Ξ π§ (cid:12)(cid:12)(cid:12) T π β (cid:2) P (cid:0) π deg (cid:1)(cid:3) π eqn βT π (cid:27) , (3.4)where π eqn is the number of equations and P (cid:0) π deg (cid:1) is the set of all polynomialswith maximum polynomial order π deg . On each mesh element we deο¬ne the π deg + π π β [β , ] for π = , . . . , π deg + Ξ¦ π ( π ) = π deg + (cid:214) π = π β π ( π β π π ) (cid:0) π π β π π (cid:1) , s.t. Ξ¦ π ( π β ) = πΏ πβ for π, β = , . . . , π deg + . (3.5)The approximate solution on each element can then be written as (cid:110) u Ξ π§ ( π‘, π§ ( π )) , π Ξ π§ ( π‘, π§ ( π )) (cid:111)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T π = π deg + βοΈ π = (cid:110) U π π ( π‘ ) , Ξ π π ( π‘ ) (cid:111) Ξ¦ π ( π ) , (3.6)where u Ξ π§ ( π‘, π§ ) : R β₯ Γ R β¦β W Ξ π§ ( π + ) π , π Ξ π§ ( π‘, π§ ) : R β₯ Γ R β¦β W Ξ π§ , Ξ¦ π ( π ) : [β , ] β¦β R , and U π ( π‘ ) : R β₯ β¦β R ( π + ) π . Similarly, we write the approximate source as Q Ξ π§ (cid:16) u Ξ π§ , π Ξ π§ (cid:17) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T π = π deg + βοΈ π = Q (cid:0) U π π ( π‘ ) , Ξ π π ( π‘ ) (cid:1) Ξ¦ π ( π ) . (3.7)To obtain the spatially discretized version of H ππ system (2.20) on each element T π , we replace the exact solution by (3.6), the exact source by (3.7), multiply (2.20) bya test function Ξ¦ π ( π ) , integrate over the element, and apply integration-by-parts in π :1 π π deg + βοΈ π = d U π π ( π‘ ) d π‘ π π π + (cid:18) Ξ π§ (cid:19) (cid:104) Ξ¦ π ( ) F π + ( π‘ ) β Ξ¦ π (β ) F π β ( π‘ ) (cid:105) β (cid:18) Ξ π§ (cid:19) A π deg + βοΈ π = U π π ( π‘ ) π π π = π deg + βοΈ π = Q (cid:0) U π π ( π‘ ) , Ξ π π ( π‘ ) (cid:1) π π π , (3.8)for each π = , . . . , π deg , where we used the fact that β π§ = ( / Ξ π§ ) β π , and where π π π = β« β Ξ¦ π ( π ) Ξ¦ π ( π ) d π and π π π = β« β Ξ¦ π ( π ) Ξ¦ (cid:48) π ( π ) d π. (3.9)The numerical ο¬ux on each element face is deο¬ned as follows: F π β ( π‘ ) = A π deg + βοΈ π = (cid:16) U π β π ( π‘ ) Ξ¦ π ( ) + U π π ( π‘ ) Ξ¦ π (β ))+ | A | π deg + βοΈ π = (cid:16) U π β π ( π‘ ) Ξ¦ π ( ) β U π π ( π‘ ) Ξ¦ π (β ) (cid:17) , (3.10) emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 13 for π = , . . . , π π§ +
1. In the above expression, | A | is deο¬ned through the eigenvalues of A : A = V π² V β = β | A | = V | π² | V β , (3.11)where V is the matrix of right eigenvectors of A , π² = diag (cid:0) π , . . . , π ( π + ) π (cid:1) is thediagonal matrix of eigenvalues of A , and | π² | = diag (cid:0) | π | , . . . , | π ( π + ) π | (cid:1) .In order to close this semi-discrete system, we also need to semi-discretizematerial-energy equation (2.29). Following all of the above outlined procedures, thisresults in the following equation: πΆ π£ π deg + βοΈ π = d Ξ π π ( π‘ ) d π‘ π π π = π deg + βοΈ π = (cid:32) π Ξ π π βοΈ β = π ( β ) π π ( ) ( π‘ ) β πππ (cid:0) Ξ π π ( π‘ ) (cid:1) (cid:33) π π π , (3.12)where π ( β ) π π ( ) refers to the ( + ( β β )( π + )) th component of U π π (i.e., the ο¬rstcomponent of U π π in the β th velocity band).For the remainder of this paper we consider the case π deg = Ξ¦ ( π ) = ( β π ) and Ξ¦ ( π ) = ( + π ) . (3.13)In this case, all of the above expressions simplify greatly. For example, the valuesdeο¬ned in (3.8) reduce to the following: M = (cid:20) (cid:21) and N = (cid:20) β β
11 1 (cid:21) . (3.14)After some simple algebra, we arrive at the following semi-discrete system:1 π d U π d π‘ β F π + + F π β β A (cid:0) U π + U π (cid:1) Ξ π§ = Q (cid:16) U π , Ξ π (cid:17) , (3.15)1 π d U π d π‘ + F π + + F π β β A (cid:0) U π + U π (cid:1) Ξ π§ = Q (cid:16) U π , Ξ π (cid:17) , (3.16) πΆ π£ d Ξ π π d π‘ = π Ξ π π βοΈ β = π ( β ) π π ( ) β πππ Ξ π π , for π = , , (3.17)for π = , . . . , ( π + ) π . The numerical ο¬uxes in the π deg = F π β ( π‘ ) = A (cid:0) U π ( π‘ ) + U π β ( π‘ ) (cid:1) β | A | (cid:0) U π ( π‘ ) β U π β ( π‘ ) (cid:1) , (3.18)for π = , . . . , π π§ + Semi-implicit time scheme: nonlinear version
For the time integration we adopt the semi-implicit scheme by McClarren et al. [21],which is a two-stage Runge-Kutta method. The ο¬rst stage (i.e., the predictor step) canbe written as follows:1 π U π + π β U ππ Ξ π‘ / = + F ππ + + F ππ β β A (cid:0) U ππ + U ππ (cid:1) Ξ π§ + Q (cid:16) U π + π , Ξ π + π (cid:17) , (3.19)1 π U π + π β U ππ Ξ π‘ / = β F ππ + + F ππ β β A (cid:0) U ππ + U ππ (cid:1) Ξ π§ + Q (cid:16) U π + π , Ξ π + π (cid:17) , (3.20) πΆ π£ Ξ π + π π β Ξ ππ π Ξ π‘ / = π Ξ π π βοΈ β = π ( β ) π + π π ( ) β πππ (cid:18) Ξ π + π π (cid:19) , for π = , . (3.21)Similarly, the second stage (i.e., the correction step) can be written as follows:1 π U π + π β U ππ Ξ π‘ = + F π + π + + F π + π β β A (cid:18) U π + π + U π + π (cid:19) Ξ π§ + Q (cid:16) U π + π , Ξ π + π (cid:17) , (3.22)1 π U π + π β U ππ Ξ π‘ = β F π + π + + F π + π β β A (cid:18) U π + π + U π + π (cid:19) Ξ π§ + Q (cid:16) U π + π , Ξ π + π (cid:17) , (3.23) πΆ π£ Ξ π + π π β Ξ ππ π Ξ π‘ = π Ξ π π βοΈ β = π ( β ) π + π π ( ) β πππ (cid:16) Ξ π + π π (cid:17) , for π = , . (3.24)The numerical ο¬ux in both the ο¬rst and second stages is of the following form: F β π β = A (cid:0) U β π + U β π β (cid:1) β | A | (cid:0) U β π β U β π β (cid:1) , (3.25)where β β (cid:8) π, π + (cid:9) .The time-stepping scheme described above is semi-implicit in that it is expliciton the wave propagation terms and implicit on the collision terms. Since in the TRTsystem the collision source is a nonlinear function of the temperature, the result is thatin each stage a nonlinear algebraic equation must be solved. In order to avoid this, weshow in the next section how to linearize the source.3.3 Semi-implicit time scheme: linearized version
The scheme we propose in this work for solving the 1D H ππ system (2.20)β(2.21)along with the material-energy equation (2.29) is a variant of the semi-implicit schemedescribed by (3.19)β(3.24), but with the additional feature that the source is linearized.The proposed scheme is a H ππ extension of the scheme developed for the P π systemin [21]. The details of this scheme are provided in this section. emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 15 We begin this discussion by recalling that source, Q , in (3.19)β(3.20) and (3.22)β(3.23) can be written as follows: Q (cid:16) U β π π , Ξ β π π (cid:17) = β π U β π π + πππ (cid:16) Ξ β π π (cid:17) Λe , (3.26)where β β (cid:8) π + , π + (cid:9) and Λe β R ( π + ) π is a vector with the following comoponents:Λ π π = (cid:40) π = + ( π β )( π + ) for π = , , . . . ,π, . (3.27)These sources are clearly nonlinear in the temperature Ξ . In order to linearize Q weinvoke the following two Taylor expansions in temperature: (cid:18) Ξ π + π π (cid:19) = (cid:16) Ξ ππ π (cid:17) + (cid:16) Ξ ππ π (cid:17) (cid:18) Ξ π + π π β Ξ ππ π (cid:19) + O (cid:32)(cid:18) Ξ π + π π β Ξ ππ π (cid:19) (cid:33) , (3.28) (cid:16) Ξ π + π π (cid:17) = (cid:16) Ξ ππ π (cid:17) + (cid:16) Ξ ππ π (cid:17) (cid:16) Ξ π + π π β Ξ ππ π (cid:17) + O (cid:18)(cid:16) Ξ π + π π β Ξ ππ π (cid:17) (cid:19) . (3.29)Since our overall method is only accurate to second order, we can safely disregard thesecond order corrections in the above expressions. Furthermore, the linear diο¬erencesin the Taylor expansions can be replaced via the update formulas (3.21) and (3.24),respectively, yielding: (cid:18) Ξ π + π π (cid:19) β (cid:16) Ξ ππ π (cid:17) + (cid:16) Ξ ππ π (cid:17) Ξ π‘ πΆ π£ (cid:32) π Ξ π π βοΈ β = π ( β ) π + π π ( ) β πππ (cid:18) Ξ π + π π (cid:19) (cid:33) , (3.30) (cid:16) Ξ π + π π (cid:17) β (cid:16) Ξ ππ π (cid:17) + (cid:16) Ξ ππ π (cid:17) Ξ π‘πΆ π£ (cid:32) π Ξ π π βοΈ β = π ( β ) π + π π ( ) β πππ (cid:16) Ξ π + π π (cid:17) (cid:33) . (3.31)Treating the above approximations as equalities and solving (3.30) and (3.31) for thefourth power of Ξ π + π π and Ξ π + π π , respectively, yields: (cid:18)(cid:101) Ξ π + π π (cid:19) : = (cid:16) Ξ ππ π (cid:17) (cid:32) πΆ π£ Ξ ππ π + Ξ π‘ Ξ ππ π βοΈ β = π ( β ) π + π π ( ) (cid:33) πΆ π£ + Ξ π‘πππ (cid:16) Ξ ππ π (cid:17) , (3.32) (cid:16)(cid:101) Ξ π + π π (cid:17) : = (cid:16) Ξ ππ π (cid:17) (cid:32) πΆ π£ Ξ ππ π + Ξ π‘ Ξ ππ π βοΈ β = π ( β ) π + π π ( ) (cid:33) πΆ π£ + Ξ π‘πππ (cid:16) Ξ ππ π (cid:17) . (3.33)Using these versions of the fourth power of the temperature successfully linearizesthe source terms in (3.19)-(3.20) and (3.22)-(3.23).In order to complete the linearization of the source terms, we now turn our attentionto (3.24). In particular, we replace the fourth power of Ξ π + π π in (3.24) by (3.33). After some simple algebra, we now arrive at the ο¬nal semi-implicit discontinuous Galerkinscheme that is advocated in this work. The ο¬rst stage is1 π U π + π β U ππ Ξ π‘ / = + F ππ + + F ππ β β A (cid:0) U ππ + U ππ (cid:1) Ξ π§ β π U π + π + πππ (cid:18)(cid:101) Ξ π + π (cid:19) Λe , (3.34)1 π U π + π β U ππ Ξ π‘ / = β F ππ + + F ππ β β A (cid:0) U ππ + U ππ (cid:1) Ξ π§ β π U π + π + πππ (cid:18)(cid:101) Ξ π + π (cid:19) Λe , (3.35)where (cid:101) Ξ is deο¬ned by (3.32) and Λe is deο¬ned by (3.27). The second stage is1 π U π + π β U ππ Ξ π‘ = + F π + π + + F π + π β β A (cid:18) U π + π + U π + π (cid:19) Ξ π§ β π U π + π + πππ (cid:16)(cid:101) Ξ π + π (cid:17) Λe , (3.36)1 π U π + π β U ππ Ξ π‘ = β F π + π + + F π + π β β A (cid:18) U π + π + U π + π (cid:19) Ξ π§ β π U π + π + πππ (cid:16)(cid:101) Ξ π + π (cid:17) Λe , (3.37) Ξ π + π π = Ξ ππ π + π Ξ π‘ (cid:34) Ξ π π βοΈ β = π ( β ) π + π π ( ) β ππ (cid:16) Ξ ππ π (cid:17) (cid:35) πΆ π£ + Ξ π‘πππ (cid:16) Ξ ππ π (cid:17) , for π = , , (3.38)where (cid:101) Ξ is deο¬ned by (3.33) and Λe is deο¬ned by (3.27). Note that by the time wereach (3.38), U π + π π is already known, meaning that this step has the computationalcomplexity of an explicit update. Again, the numerical ο¬uxes in both the ο¬rst andsecond stages are of the form (3.25).3.4 Slope limiter
In order to remove unphysical oscillations from the numerical method describedabove, we include a slope limiter. McClarren and Lowrie [23] pointed out that the double minmod slope limiter, also known as the monotonized central slope limiter,is asymptotic-preserving for hyperbolic systems with stiο¬ relaxation terms while theminmod limiter does not preserve the asymptotic limit due to discontinuities at thecell edge. Therefore, in this work we use the same double minmod limiter to preserveasymptotic limit. In particular, after each predictor and corrector step we compute thecell average: U π = U π + U π , (3.39) emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 17 and then modify the original nodal values as U π : = U π β s π U π : = U π + s π , (3.40)where π π ( β ) = mm (cid:16) π π ( β ) β π π ( β ) , πΌ (cid:16) π π ( β ) β π π β ( β ) (cid:17) , πΌ (cid:16) π π + ( β ) β π π ( β ) (cid:17)(cid:17) , (3.41)where πΌ β [ , ] , π = , Β· Β· Β· , π π§ , π = , Β· Β· Β· ,π ( π + ) , and the minmod function is deο¬nedas follows:mm ( π, π, π ) : = (cid:40) sign ( π ) min (| π | , | π | , | π |) , if sign ( π ) = sign ( π ) = sign ( π ) , , otherwise . (3.42)As explained in [21], πΌ = πΌ = πΌ = πΌ = Boundary conditions
To complete the numerical methods section, we brieο¬y explain how boundary con-ditions are implemented. The three types of boundary conditions considered in thiswork in various examples are reο¬ective, Dirichlet, and vaccuum conditions. In allcases we prescribe the intensities on the left and right boundaries via the followingexpressions: πΌ πΏ ( π‘, π ) = (cid:40) πΌ out ( π‘, π§ πΏ , π ) if π > ,πΌ ( π‘, π§ πΏ , π ) if π < , πΌ π ( π‘, π ) = (cid:40) πΌ ( π‘, π§ π , π ) if π > ,πΌ out ( π‘, π§ π , π ) if π < , (3.43)respectively, where πΌ out ( π‘, π§, π ) =  πΌ ( π‘, π§, β π ) if reο¬ective BC ,πΌ ( π‘, π§, π ) if Dirichlet BC , . (3.44)Note that π§ πΏ and π§ π denote the left and right boundaries, respectively. In this section, numerical results for six standard benchmark problems for the TRTsystem are provided, including examples in the optically thin and thick regimes: (1)bilateral inο¬ow (Β§4.1), (2) streaming in a vaccuum (Β§4.2), (3) Su-Olson problem(Β§4.3), (4) diο¬usive Marshak wave (Β§4.4), (5) Marshak wave in thin medium (Β§4.5),and (6) smooth Marshak wave problem (Β§4.6). The material-energy coupling equation(2.29) is only required for in the examples described in Β§4.3βΒ§4.6.
Unless otherwise stated, we choose a CFL condition for all simulations as follows:CFL : = π π ( A ) Ξ π‘ Ξ π₯ β€ . , (4.1)where π ( A ) β A as deο¬ned in equation (2.25). Fornumerical examples in which the diο¬usion dominates, for example problems like thediο¬usive Marshak wave problem, we are able to achieve very relaxed CFL numberbetween 2 and 3.4.1 Bilateral inο¬ow
This problem is used to test the H ππ scheme without the opacity π and the externalsource π ; this shows how well the hybrid discrete approximation of the free-streamingkinetic operator and its corresponding numerical discretization captures the correctwave speeds and resolves discontinuities. In this setting, equation (2.7) reduces to thefollowing: 1 π β πΌ β π‘ + π β πΌ β π§ = , (4.2)for which the analytic solution can easily be computed from the method of character-istics (e.g., see [7]): πΌ ( π‘, π§, π ) = πΌ ( π§ β πππ‘, π ) . (4.3)We choose the initial condition as πΌ ( π§, π ) =  πππΏ ( π β ) if π§ β€ . , . < π§ < . , . ππ if π§ > . , (4.4)where πΏ ( π β ) is a Dirac delta centered at π = π = Γ cm s β is the speed oflight, and π = . Γ ergs cm β keV β is the radiation constant. The exact angularmoment of the radiation intensity, πΈ ( π‘, π§ ) , is given by the following if 0 < ππ‘ < . π β πΈ ( π‘, π§ ) =  π§ β€ . + ππ‘, . + ππ‘ < π§ β€ . β ππ‘, ( π§ β . + ππ‘ ) /( ππ‘ ) if 0 . β ππ‘ β€ π§ β€ . + ππ‘, π§ β₯ . + ππ‘. (4.5)We run the code until ππ‘ end = .
1, in the physical domain π§ β [ , ] with inο¬ow/outο¬owboundary conditions based on the exact solution (4.3)β(4.4) . This example is chal-lenging for moment closure methods such as P π and H ππ due to the fact that there isa delta function in velocity, as well as discontinuities in both π§ and π . Furthermore,in this example there are no collisions to help smooth out the solution.The scaled angular moment of the radiation intensity, πΈ ( π‘, π§ )/ π , for various P π and H ππ approximations are shown in Figure 4.1. Each panel in Figure 4.1 showssolutions with models that have the same DOFs: (a) 4 moments: P and H , (b) 6 emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 19 moments: P , H , and H , (c) 8 moments: P , H , and H , and (d) 24 moments: P ,H , and H . In each case, the H ππ model with π > π model with the same DOFs.In all the simulations shown in Figure 4.1 we have used the double minmod limiterto control unphysical oscillations and to remove negative density values; without thedouble minmod limiters active, both P π and H ππ solutions suο¬er from negativedensities near the discontinuity. In order to show how the double minmod limiteraο¬ects the solution we also provide Figure 4.2, in which we show a direct comparisonof the unlimited and limited H solutions. (a) Bilateral inο¬ow: P vs H ππ z -0.500.511.52 E / a H P ref (b) Bilateral inο¬ow: P vs H ππ z -0.500.511.52 E / a H H P ref (c) Bilateral inο¬ow: P vs H ππ z -0.500.511.52 E / a H H P ref (d) Bilateral inο¬ow: P vs H ππ z -0.500.511.52 E / a H H P ref Fig. 4.1: (Β§4.1: Bilateral inο¬ow) Comparisons of P π and H ππ solutions for the bilateralinο¬ow problem at ππ‘ = . π π§ = = . Streaming in a vacuum
In this section, we test our scheme on the propagation of photons in a vacuum. Wechoose a zero initial condition, πΌ ( π§, π ) =
0, on the computational domain π§ β [ , ] z -0.500.511.52 E / a H w/o limiterH w/ limiter Fig. 4.2: (Β§4.1: Bilateral inο¬ow) H solution with and without the double minmodlimiter for the bilateral inο¬ow problem at ππ‘ = . π π§ = = . πΌ ( π‘, π§ = , π ) = ππ and πΌ ( π‘, π§ = , π ) = . (4.6)The analytic solution for π‘ > πΌ ( π‘, π§, π ) = (cid:40) ππ if π§ /( ππ‘ ) < π β€ , β β€ π β€ π§ /( ππ‘ ) . (4.7)The exact solution for the angular moment for π‘ > π β πΈ ( π‘, π§ ) = β« π§ /( ππ‘ ) d π = β π§ /( ππ‘ ) . (4.8)In Figure 4.3 and Table 4.1 we compute the absolute πΏ error in the scaled angularmoments, πΈ ( π‘, π§ )/ π , for various H ππ approximations: πΏ error : = (cid:34) ππ π§ π π§ βοΈ π = (cid:16) πΈ ππ ( π‘ end , π§ π ) β πΈ exact ( π‘ end , π§ π ) (cid:17) (cid:35) / , (4.9) πΏ β error : = max β€ π β€ π π§ (cid:12)(cid:12)(cid:12)(cid:12) π (cid:16) πΈ ππ ( π‘ end , π§ π ) β πΈ exact ( π‘ end , π§ π ) (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) , (4.10)where πΈ ππ is the angular moment solution calculated using the H ππ approximation. Inparticular, we show in Figure 4.3a the πΏ convergence of H ππ with increasing π andin Figure 4.3b the πΏ convergence of H ππ with increasing π . Figure 4.3a shows rapidconvergence in terms of π for the H π and H π solutions, however P π =(H π ) shows amuch slower convergence rate. Meanwhile, Figure 4.3b shows rapid convergence witha much steeper slope than Figure 4.3a, which suggests that increasing π is a better emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 21 way to achieve the desired accuracy than increasing π when a discontinuity exists inthe underlying intensity πΌ ( π‘, π§, π ) . This observation agrees with the numerical valuesshown in Table 4.1, where the πΏ and πΏ β errors are shown for various methods thatall have the same degrees of freedom (DOF = πΈ ( π‘, π§ ) for various methods in Figure 4.4.Each H ππ solution in Figure 4.4 shows the propagation of multiple waves dependingon the number of distinct eigenvalues of the matrix A deο¬ned in (2.18) and (2.21).In particular, in Figure 4.4a and Figure 4.4b we compare H π and H π solutions withvarious π , respectively. Analogously, in Figure 4.4c and 4.4d we compare variousH π and H π solutions with diο¬erent π , respectively. Finally, in Figure 4.4e through4.4h we plot in each panel a diο¬erent H ππ method with DOF =
16: (e) H , (f) H ,(f) H , and (g) H . Again, we conclude from these simulations that in the case whenthe intensity is discontinuous, increasing π is more eο¬ective than increasing π . Wealso demonstrate in these panels that when π is odd and π is even (i.e., A has a zeroeigenvalue β see Remark 2.1), the solution degrades signiο¬cantly: Figure 4.4a (H andH ) and in Figure 4.4d (H and H ). (a) H π (=P π ) vs H π vs H π N -2 -1 L - e rr o r H H H (b) H π vs H π vs H π T -2 -1 L - e rr o r H T1 H T3 H T5 Fig. 4.3: (Β§4.2: Streaming in a vacuum) Convergence study of H ππ solutions for thevacuum propagation problem at π‘ = . Γ β s with π π§ =
100 and CFL = . π and (b) π . error H H H H ( = P ) πΏ Γ β Γ β Γ β Γ β πΏ β Γ β Γ β Γ β Γ β Table 4.1: (Β§4.2: Streaming in a vacuum) πΏ and πΏ β errors for the vacuum propagationproblem for various methods that all have the same degrees of freedom (DOF = π z E / a H H H H (b) Vacuum propagation: H π z E / a H H H H (c) Vacuum propagation: H π z E / a H H H H (d) Vacuum propagation: H π z E / a H H H H (e) Vacuum propagation: H z E / a H exact (f) Vacuum propagation: H z E / a H exact (g) Vacuum propagation: H z E / a H exact (h) Vacuum propagation: H z E / a H exact Fig. 4.4: (Β§4.2: Streaming in a vacuum) The H ππ solutions for the vacuum propagationproblem at π‘ = . Γ β with π π§ =
100 and CFL = . emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 23 Su-Olson problem
The next benchmark problem we solve is the non-equilibrium Su-Olson problem withmaterial coupling [32]. In order to compare the numerical solutions to the semi-analytic solution, we follow the conditions used in [32], i.e., π =
1, and the externalsource term is given by π ( π‘, π§ ) = (cid:40) ππ if β . β€ π§ β€ . , π§ β (cid:2) β ππ‘ end β , ππ‘ end + (cid:3) with periodic boundaryconditions. For this problem, both P π = H π and H ππ perform very well even withsmall π .Figure 4.5 shows the results from three diο¬erent simulations with ππ‘ end =1.0, 3.16,and 10.0, respectively. The solutions to these three cases are shown on a linear scalein Figure 4.5a (P and H ) and 4.5b (P and H ), and on a log-log scale in Figure 4.5c(P and H ) and 4.5d (P and H ). The reference solutions are obtained from [32].4.4 Diο¬usive Marshak-wave problem
In this section, a diο¬usive Marshak-wave problem is investigated; this problem is astandard test case in the literature [6,19,21,34]. This problem consists of a semi-inο¬nite medium of material with the opacity π = / π . (4.12)The computational domain is π§ β [ . , . ] and the initial conditions are given by thefollowing with π = β keV: πΌ ( π‘ = , π§, π ) = πππ and π ( π‘ = , π§ ) = π . (4.13)We use the isotropic incoming boundary condition corresponding to a 1 keVtemperature source on the left boundary π§ πΏ =
0, and no incoming radiation on theright boundary at π§ π = . πΌ ( π‘, π§ πΏ , π > ) = ππ and πΌ ( π‘, π§ π , π < ) = . (4.14)We compute the material temperature, π ( π‘, π§ ) at various times: π‘ = β s, 5 Γ β s,and 10 β s. In our numerical tests for this problem all H ππ solutions look similar, thus,we only present H solutions in Figure 4.6. In this test, we use π π§ =
15 with themesh size Ξ π§ = .
04. Despite the fact that we use a coarse mesh, the H ππ solution isable to adequately capture the wave propagation front. Due to the fact that diο¬usiondominates, we are able to achieve stable results with a relaxed Courant number:CFL = .
7. The reference solution shown in this plot is the semi-analytic equilibrium-diο¬usion solution (e.g., see [21]). vs H z E / a H P exact (b) Su-Olson: P vs H z E / a H P exact (c) Su-Olson: P vs H (log-log) -1 z -2 -1 l og ( E / a ) H P exact (d) Su-Olson: P vs H (log-log) -1 z -2 -1 l og ( E / a ) H P exact Fig. 4.5: (Β§4.3: Su-Olson problem) Three diο¬erent simulations using H ππ for theSu-Olson problem with π π§ = = .
3, and ππ‘ end = . , . ,
10, respectively.Panels (a) and (b) show P , H , P , and H on a linear scale, while panels (c) and (d)show those same solutions on a log-log scale.4.5 Marshak-wave in thin medium
Here we apply our scheme to a Marshak-wave problem in an optically thin mediumwith an opacity given by π = / π . (4.15)The radiation temperature is given by π rad ( π‘, π§ ) : = βοΈ πΈ ( π‘, π§ ) π , (4.16)and the computational domain is π§ β [ , . ] . We use the isotropic incoming boundarycondition corresponding to a 1 keV temperature source on the left boundary π§ πΏ = π§ π = . πΌ ( π‘, π§ πΏ , π > ) = ππ and πΌ ( π‘, π§ π , π < ) = . (4.17) emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 25 z ( k e V ) ref H Fig. 4.6: (Β§4.4: Diο¬usive Marshak-wave problem) The H solutions for the materialtemperature, π ( π‘, π§ ) , for the diο¬usive Marshak-wave problem at various times: π‘ = β s, 5 Γ β s, and 10 β s, respectively, with π π§ = π = / π , π = β keV,and CFL = .
7. The reference solution is the semi-analytic equilibrium-diο¬usionsolution.The initial conditions are given by the following with π = β keV: πΌ ( π‘ = , π§, π ) = πππ and π ( π‘ = , π§ ) = π . (4.18)In Figure 4.7a we show the H solution with π π§ =
400 for the material temperature, π , and the radiation temperature, π rad (4.16). In Figure 4.7b we show the materialtemperature error for H for various π π§ . The error is computed by comparing the H solutions with various π π§ to the P solution with π π§ = π , H π ,H π and H π with π π§ =
512 as a function of π . In particular, we compute the followingapproximate πΏ error: πΏ error : = (cid:34) . π π§ π π§ βοΈ π = (cid:16) π π π β ( π‘ end , π§ π ) β π π ( π‘ end , π§ π ) (cid:17) (cid:35) / , (4.19)for π = , , . . . ,
9, where π ππ represents the material temperature as calculated withthe H ππ model. Due to the non-smoothness of the solution we note a fairly rapidconvergence as a function of π with ο¬xed π , and slower convergence as a function of π with ο¬xed π .4.6 Smooth Marshak-wave
Finally, a smooth Marshak-wave problem is considered in this section to observethe convergence rate of our numerical scheme on a smooth solution in an optically solutions with π π§ = z ( k e V ) radiation material (b) H convergence rate in position space -4 -3 -2 -1 z -3 -2 -1 L - e rr o r H slope 1 Fig. 4.7: (Β§4.5: Marshak-wave in thin medium) The H solutions and convergence ratefor the thin Marshak-wave problem at π‘ = β s with π = / π , π = β keV, andCFL = . N -8 -6 -4 -2 L - e rr o r O(N )O(N )O(N ) H H H H Fig. 4.8: (Β§4.5: Marshak-wave in thin medium) Error in various H ππ approximationswith π π§ =
512 and CFL = . π§ β [ , . ] andthe smooth initial conditions are given by: πΌ ( π‘ = , π§, π ) = ππ (cid:20) β . (cid:16) + tanh (cid:2) ( π§ β . ) (cid:3) (cid:17)(cid:21) , (4.20) π ( π‘ = , π§ ) = (cid:18) πΈ ( π‘ = , π§ ) π (cid:19) . (4.21) emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 27 The boundary conditions are given by πΌ ( π‘, π§ πΏ , π > ) = πΌ ( π‘ = , π§ πΏ , π ) and πΌ ( π‘, π§ π , π < ) = , (4.22)where π§ πΏ = π§ π = . π ( π‘, π§ ) , at varioustimes: (a)The initial condition for the material temperature, π ( π‘, π§ ) , is shown in Panel (a) ofFigure 4.9. In Figure 4.9(b) we show a direct comparison of the material temperature ascomputed with π and π» , both methods have DOF =
6, for π π§ = π and π» withvarious π π§ , respectively. (a) The initial condition z ( k e V ) initial T (b) P vs H with π π§ = z ( k e V ) P H (c) P solutions with various π π§ z ( k e V ) N z =8N z =16N z =32N z =64N z =128 (d) H solutions with various π π§ z ( k e V ) N z =8N z =16N z =32N z =64N z =128 Fig. 4.9: (Β§4.6: Smooth Marshak-wave) Comparison of the P and H π solutions forthe material temperature for the smooth Marshak wave problem with CFL = . π» withincreasing mesh resolution π π§ . In this ο¬gure, the dashed-line indicates a slope oftwo on a log-log scale. In these convergence experiments the approximate πΏ error is calculated via the formula: πΏ ( π β π ) : = (cid:118)(cid:117)(cid:116) . π π§ π π§ βοΈ π = (cid:16) π π β π /( π ) π (cid:17) , (4.23)where π /( π + ) π = (cid:16) π /( π ) ( π β )+ + π /( π ) ( π β )+ (cid:17) , (4.24)for π = , , Β· Β· Β· ,
7. Here π represents the reference temperature solution on meshwith π π§ = π stands for the π th grid cell. The ideaencapsulated in formulas (4.23) and (4.24) is that we project the reference solution π onto coarser mesh, i.e., π π§ = /( π ) , for π = , , Β· Β· Β· , π , π , Β· Β· Β· , π ; once we have projectedthis solution down to the mesh on which π is deο¬ned, we can directly compute the πΏ distance. In Figure 4.10 we have used this strategy with π π§ =
8, 16, 32, 64, 128, 256,512, and 1024 to show the convergence rate. -3 -2 -1 z -6 -5 -4 -3 -2 -1 L - e rr o r H slope 2 Fig. 4.10: (Β§4.6: Smooth Marshak-wave) Convergence rate of the material temperature π for the smooth Marshak wave problem. In this work we have developed the hybrid discrete (H ππ ) approximation methodfor the thermal radiative transfer (TRT) equations, and implemented a numericaldiscretization of these equations using a second order discontinuous Galerkin ο¬niteelement method in conjunction with a semi-implicit time-stepping scheme. The H ππ approach acquires desirable properties of two classical methods for the TRT equations,namely P π (spherical harmonics) and S π (discrete ordinates), and indeed reducesto each of these approximations in various limits: H π β‘ P π and H π β‘ S π . Weproved that the H ππ approximation results in a system of hyperbolic partial diο¬erentialequations for all π β₯ π β₯
0. In particular, in one spatial dimension, the H ππ emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 29 scheme is essentially a collection of P π approximations localized in each velocityband. Because of this structure, H ππ , just like P π , can exhibit negative densities.However, because H ππ has band-localized structure, we are able to control unphysicalnumerical oscillations by increasing the number of discrete regions in velocity space,thereby blending properties of P π and S π .Once the H ππ approximation was developed for TRT, we introduced a semi-implicit numerical method that is based on a second order explicit Runge-Kuttascheme for the streaming term and an implicit Euler scheme for the material couplingterm. Furthermore, in order to solve the material energy equation implicitly aftereach predictor and corrector step, we linearized the temperature term using a Taylorexpansion; this avoided the need for an iterative procedure, and therefore improvedeο¬ciency. In order to reduce unphysical oscillation, we applied a slope limiter aftereach time step.In the numerical results section we compared the solutions of the H ππ and P π schemes for the various benchmark problems. We demonstrated for a variety ofproblems that for a ο¬xed total number of moments, we are able to achieve accuracygains over the P π β‘ H π approximation by balancing π and π . With a more balancedchoice of π and π , H ππ shows less oscillation than P π , especially in the presence ofdiscontinuities.In future work, we will develop extensions of the H ππ method for multi-energygroup models or frequency-dependent equations in multiple dimensions. In particular,just as we have done in this work, we will investigate various choices of π and π toachieve accurate and eο¬cient moment closure in the multidimensional setting. Wewill also investigate adaptive strategies for selecting π and π in the presence of someappropriate error indicator. Acknowledgements
JAR was supported in part by NSF Grants DMSβ1620128 and DMSβ2012699.
References
1. Adams, M.P., Adams, M.L., Hawkins, W.D., Smith, T., Rauchwerger, L., Amato,N.M., Bailey, T.S., Falgout, R.D., Kunen, A., Brown, P.: Provably optimal paralleltransport sweeps on semi-structured grids. J. Comput. Phys. , 109β234 (2020)2. Brunner, T.A., Holloway, J.P.: Two-dimensional time-dependent Riemann solversfor neutron transport. J. Comput. Phys. (1), 386β399 (2005)3. Carlson, B.: Solution of the transport equation by the S π method. Los AlamosNational Laboratory (1955)4. Carlson, B.: Tables of symmetric equal weight quadrature EQn over the unitsphere. Los Alamos National Laboratory (1971)5. Cohen, A.: An algebraic approach to certain diο¬erential eigenvalue problems.Linear Algebra Appl. , 183β198 (1996)6. Evans, T.M., Urbatsch, T.J., Lichtenstein, H., Morel, J.E.: A residual Monte Carlomethod for discrete thermal radiative diο¬usion. J. Comput. Phys. (2), 539β556(2003)
7. Fan, Y.W., Li, R., Zheng, L.C.: A Nonlinear Moment Model for Radiative TransferEquation in Slab Geometry. J. Comput. Phys. , 109β128 (2020)8. Fleck, J.A., Jr., Cummings, J.D.: An implicit Monte Carlo scheme for calculatingtime and frequency dependent nonlinear radiation transport. J. Comput. Phys. (3), 313β342 (1971)9. Hauck, C.D., McClarren, R.G.: Positive P π closures. SIAM J. Sci. Comput. (5), 2603β2626 (2010)10. Jarrell, J., Adams, M.: Discrete-ordinates quadrature sets based on linear dis-continuous ο¬nite elements. In: Proc. International Conference on Mathematicsand Computational Methods applied to Nuclear Science and Engineering, Rio deJaneiro, Brazil (2011)11. Klar, A.: An asymptotic-induced scheme for nonstationary transport equations inthe diο¬usive limit. SIAM J. Numer. Anal. (3), 1073β1094 (1998)12. Klar, A., Unterreiter, A.: Uniform stability of a ο¬nite diο¬erence scheme for trans-port equations in diο¬usive regimes. SIAM J. Numer. Anal. (3), 891β913 (2002)13. Koch, R., Krebs, W., Wittig, S., Viskanta, R.: The discrete ordinate quadratureschemes for multidimensional radiative transfer. J. Quant. Spectrosc. Ra. (4),353β372 (1995)14. Lathrop, K.D.: Ray eο¬ects in discrete ordinates equations. Nucl. Sci. Eng. (3),357β369 (1968)15. Lathrop, K.D.: Remedies for ray eο¬ects. Nucl. Sci. Eng. (3), 255β268 (1971)16. Lathrop, K.D., Carlson, B.G.: Discrete ordinates angular quadrature of the neutrontransport equation. Los Alamos Scientiο¬c Laboratory Report 3186 (1965)17. Lau, C., Adams, M.: Discrete Ordinates Quadratures Based on Linear andQuadratic Discontinuous Finite Elements over Spherical Quadrilaterals. Nucl.Sci. Eng. (1), 36β52 (2017)18. Lewis, E.E., Miller, W.F.: Computational Methods of Neutron Transport. JohnWiley & Sons, DeKalb (1994)19. Li, W., Liu, C., Zhu, Y., Zhang, J., Xu, K.: Uniο¬ed gas-kinetic wave-particlemethods III: Multiscale photon transport. J. Comput. Phys. , 109β280 (2020)20. Lowrie, R.B.: A Comparison of Implicit Time Integration Methods for NonlinearRelaxation and Diο¬usion. J. Comput. Phys. (2), 566β590 (2004)21. McClarren, R.G., Evans, T.M., Lowrie, R.B., Densmore, J.D.: Semi-implicit timeintegration for P π thermal radiative transfer. J. Comput. Phys. (16), 7561β7586 (2008)22. McClarren, R.G., Holloway, J.P., Brunner, T.A.: On solutions to the pn equationsfor thermal radiative transfer. Journal of Computational Physics (5), 2864β2885 (2008)23. Mcclarren, R.G., Lowrie, R.B.: The eο¬ects of slope limiting on asymptotic-preserving numerical methods for hyperbolic conservation laws. J. Comput.Phys. (23), 9711β9726 (2008)24. McClarren, R.G., Urbatsch, T.J.: A modiο¬ed implicit monte carlo method fortime-dependent radiative transfer with adaptive material coupling. Journal ofComputational Physics (16), 5669β5686 (2009)25. Olson, G.: Second-order time evolution of P N equations for radiation transport.J. Comput. Phys. (8), 3072β3083 (2009) emi-implicit Hybrid Discrete Approximation of Thermal Radiative Transfer 31
26. Olson, G.L., Auer, L.H., Hall, M.L.: Diο¬usion, p , and other approximate formsof radiation transport. J. Quant. Spectrosc. Ra. (6), 619β634 (2000)27. Parlett, B.N.: The Symmetric Eigenvalue Problem. Prentice?Hall Series in Com-put. Math. (7), 277β348 (1981)28. Pomraning, G.C.: Variational boundary conditions for the spherical harmonicsapproximation to the neutron transport equation. Ann. Phys. , 193β215 (1964)29. Pomraning, G.C.: The Equations of Radiation Hydrodynamics. Pergamon Press,Oxford, U.K. (1973)30. Shin, M.: Hybrid discrete (H ππ ) approximations to the equation of radiative trans-fer. Ph.D. thesis, Iowa State University, Ames, IA (2019)31. Siegel, R., Howell, J.R.: Thermal radiation heat transfer, vol. 3. HemispherePublishing Corp., Washington, United States (1972)32. Su, B., Olson, G.L.: An analytical benchmark for non-equilibrium radiative trans-fer in an isotropically scattering medium. Ann. Nucl. Energy (13), 1035β1055(1997)33. Thurgood, C.P., Pollard, A., Becker, H.A.: The T π Quadrature Set for the DiscreteOrdinates Method. J. Heat Transfer (4), 1068β1070 (1995)34. Vikas, V., Hauck, C., Wang, Z., Fox, R.: Radiation transport modeling usingextended quadrature method of moments. J. Comput. Phys. (1), 221β241(2013)35. Wikipedia: Gaussian quadrature (2021). URL https://en.wikipedia.org/wiki/Gaussian_quadrature
36. Wollaber, A.B.: Four decades of implicit monte carlo. Journal of Computationaland Theoretical Transport45