On a convergent DSA preconditioned source iteration for a DGFEM method for radiative transfer
OOn a convergent DSA preconditioned source iteration for a DGFEMmethod for radiative transfer
Olena Palii, Matthias Schlottbom ∗ Department of Applied Mathematics, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands
Abstract
We consider the numerical approximation of the radiative transfer equation using discontinuous angular andcontinuous spatial approximations for the even parts of the solution. The even-parity equations are solvedusing a diffusion synthetic accelerated source iteration. We provide a convergence analysis for the infinite-dimensional iteration as well as for its discretized counterpart. The diffusion correction is computed by asubspace correction, which leads to a convergence behavior that is robust with respect to the discretization.The proven theoretical contraction rate deteriorates for scattering dominated problems. We show numericallythat the preconditioned iteration is in practice robust in the diffusion limit. Moreover, computations for thelattice problem indicate that the presented discretization does not suffer from the ray effect. The theoreticalmethodology is presented for plane-parallel geometries with isotropic scattering, but the approach and proofsgeneralize to multi-dimensional problems and more general scattering operators verbatim.
Keywords: radiative transfer, discontinuous angular approximation, discrete ordinates method, diffusionsynthetic acceleration, convergence rates
1. Introduction
We consider the numerical solution of the radiative transfer equation in plane parallel geometry µ∂ z φ ( z, µ ) + σ t ( z ) φ ( z, µ ) = σ s ( z )2 (cid:90) − φ ( z, µ (cid:48) ) dµ (cid:48) + q ( z, µ ) , (1)where 0 < z < Z and − < µ <
1, and Z denotes the thickness of the slab and µ is the cosine of thepolar angle of a unit vector. The function φ ( z, µ ) models the equilibrium distribution of some quantity, likeneutrons or photons [1, 2]. The basic physical principles embodied in (1) are transport, which is modeled bythe differential operator µ∂ z , attenuation with rate σ t and scattering with σ s . Internal sources are describedby the function q . In this work we will close the radiative transfer equation using inflow boundary conditions φ (0 , µ ) = g ( µ ) µ > , and φ ( Z, µ ) = g Z ( µ ) µ < . (2)Such transport problems arise when the full three-dimensional model posed on R × (0 , Z ) × S obeyscertain symmetries [2]. It has been studied in many instance due to simpler structure compared to threedimensional problems without symmetries; while the methodology presented here directly carries over tothe general case, see also the numerical examples presented below.Classical deterministic discretization strategies are based on a semidiscretization in µ . One class of suchstrategies are the P N -approximations, which are spectral methods based on truncated spherical harmonics ∗ Corresponding author
Email addresses: [email protected] (Olena Palii), [email protected] (Matthias Schlottbom)
Preprint submitted to Computers and Mathematics with Applications September 19, 2019 a r X i v : . [ m a t h . NA ] S e p xpansions [3, 4], and we refer to [5, 6, 7] for variational discretization strategies using this approximation.The major advantage of P N -approximations is that the scattering operator becomes diagonal. In addition,the matrix representation of the transport operator µ∂ z is sparse. The main drawbacks of the P N -methodare that the variational incorporation of the inflow boundary condition introduces a dense coupling of thespherical harmonics expansion coefficients making standard P N -approximations quite expensive to solve.We note, however, that a modified variational formulation of the P N -equations has recently been derivedthat leads to sparse matrices [8]. In any case, the success of spectral approximation techniques depends onthe smoothness of the solution. In general, the solution φ is not smooth for µ = 0, which is related to theinflow boundary conditions (2). Hence, the P N -approximations will in general not converge spectrally.A second class of semidiscretizations are discrete ordinates methods that use a quadrature rule for thediscretization of the µ -variable [2], with analysis provided in [9, 10]. Such methods are closely related to dis-continuous Galerkin (DG) methods, see, e.g., [11, 12, 13, 14, 15]. While allowing for local angular resolution,the main obstruction in the use of these methods is that the scattering operator leads to dense matrices,and a direct inversion of the resulting system is not possible in realistic applications. To overcome this issue,iterative solution techniques have been proposed. An often used iterative technique is Richardson iteration,i.e., the source iteration [16, 17], but other Krylov space methods exist [18]. The key idea in the sourceiterations is to decouple scattering and transport, and to exploit that the transport part can be invertedefficiently. If σ s /σ t ≈
1, the convergence of these iterative methods is slow, and several preconditioningtechniques have been proposed [14, 16, 17, 18]. Among the most used and simple preconditioners is thediffusion synthetic acceleration method (DSA), in which a diffusion problem is solved in every iteration.This is well motivated by asymptotic analysis [19, 20]. While Fourier analysis can be applied to specialsituations [16, 17], the convergence analysis is mainly open for the general case, i.e., for jumping coeffi-cients or non-periodic boundary conditions. A further complication in using the DSA preconditioner is thatthe resulting iterative scheme might diverge if the discretization of the diffusion problem and the discreteordinates system is not consistent [17].The contribution of this paper is to develop discretizations that allow for local resolution of the non-smoothness of the solution, and which lead to discrete problems that can be solved efficiently by diffusionsynthetic accelerated source iterations. Our approach builds upon an even-parity formulation of the radia-tive transfer equations derived from the mixed variational framework given in [7], where P N -approximationshave been treated in detail. Using the framework of KP -methods [16], we show that an infinite dimensionalDSA preconditioned source iteration converges already for the continuous problem. We present conforming hp -type approximations spaces for the semidiscretization in µ , and prove quasi-best approximation prop-erties. In order to solve resulting linear systems, we employ a DSA preconditioned Richardson iteration,which is just the infinite dimensional iteration projected to the approximation spaces. In particular thefinite dimensional iteration is guaranteed to converge for any discretization. Moreover, the inversion of thetransport problem can be parallelized straightforward. In numerical experiments, even when employinglow-order approximations, we observe that the developed method does not suffer from the ray effect, whichis typically observed for discrete ordinates methods [4, 21].The outline of the paper is as follows: In Section 2 we introduce basic notation and recall the relevantfunction spaces. In Section 3 we introduce the even-parity equation for (1), show its well-posedness, andformulate an infinite dimensional DSA preconditioned source iteration, for which we show convergence. Theapproximation spaces are described in Section 4 and well-posedness of the Galerkin problems as well asquasi-best approximation results are presented. In Section 5 we discuss the efficient iterative solution of theresulting linear systems and provide a convergence proof for the discrete DSA scheme. Section 6 presentssupporting numerical examples for slab geometry and for multi-dimensional problems that show the goodapproximation properties of the proposed methods as well as fast convergence of the iterative solver inmulti-dimensional problems and in the diffusion limit. The paper ends with some conclusions in Section 7.2 . Function spaces and further preliminaries Following [22] we denote by L ( D ) with D = (0 , Z ) × ( − ,
1) the usual Hilbert space of square integrablefunctions with inner product ( φ, ψ ) = (cid:90) − (cid:90) Z φ ( z, µ ) ψ ( z, µ ) dzdµ and induced norm (cid:107) φ (cid:107) L ( D ) = ( φ, φ ) . Furthermore, we define the Hilbert space H ( D ) = { φ ∈ L ( D ) : µ∂ z φ ∈ L ( D ) } of functions with square integrable weak derivatives with respect to the weighted differential operator µ∂ z endowed with the corresponding graph norm.In order to deal with boundary data, let us introduce the Hilbert space L − that consists of measurablefunctions for which (cid:107) ψ (cid:107) L − = (cid:90) | ψ (0 , µ ) | | µ | dµ + (cid:90) − | ψ ( Z, µ ) | | µ | dµ is finite, and we denote by (cid:104) ψ, φ (cid:105) L − the corresponding inner product on L − . Similarly, L denotes the spaceof outflow data. We have the following trace lemma [22]. Lemma 2.1. If φ ∈ H ( D ) , then there exist traces φ | Γ − ∈ L − and φ | Γ + ∈ L and (cid:107) φ | Γ ± (cid:107) L ± ≤ C √ − e − Z (cid:107) φ (cid:107) H ( D ) with a constant C > independent of φ and Z . As a consequence of the trace lemma and the density of smooth functions in H ( D ) the followingintegration-by-parts formula is true( µ∂ z φ, ψ ) = − ( φ, µ∂ z ψ ) + (cid:104) φ, ψ (cid:105) L − (cid:104) φ, ψ (cid:105) L − . (3)Throughout the manuscript we make the following basic assumption:(A1) σ s , σ t ∈ L ∞ (0 , Z ) are non-negative and σ a = σ t − σ s ≥ γ > Lemma 2.2.
Let assumption (A1) hold, and let g − ∈ L − and q ∈ L ( D ) , then (1) - (2) has a unique solution φ ∈ H ( D ) that satisfies the a-priori bound (cid:107) φ (cid:107) H ( D ) ≤ C ( (cid:107) g − (cid:107) L − + (cid:107) q (cid:107) L ( D ) ) . Assumption (A1) is not required to prove well-posedness for bounded geometries [23], or for slab problemswith constant coefficients [24, Thm. 2.25].Let P : L ( D ) → L ( D ) denote the L -projection onto constants in µ , i.e.,( P ψ )( z, µ ) = 12 (cid:90) − ψ ( z, µ (cid:48) ) dµ (cid:48) . Since σ t ∈ L ∞ (0 , Z ) is strictly positive, we can define the norms on L ( D ) as follows (cid:107) ψ (cid:107) σ t = ( σ t ψ, ψ ) and (cid:107) ψ (cid:107) σt = ( 1 σ t ψ, ψ ) . (4)3 .1. Even-odd splitting The even φ + and odd φ − parts of a function φ ∈ L ( D ) are defined as φ + ( z, µ ) = 12 ( φ ( z, µ ) + φ ( z, − µ )) , φ − ( z, µ ) = 12 ( φ ( z, µ ) − φ ( z, − µ )) . Even-odd decompositions are frequently used in transport theory [5, 3]. Since, as functions of µ , even andodd function are orthogonal in L ( − , L ( D ) into orthogonal subspaces containingeven and odd functions, respectively, L ( D ) = L ( D ) + ⊕ L ( D ) − . Similarly, we will write H ( D ) ± = H ( D ) ∩ L ( D ) ± . As in [7], we observe that µ∂ z φ ± ∈ L ( D ) ∓ for any φ ∈ H ( D ), and P φ ± ∈ L ( D ) + for φ ∈ L ( D ). It turns out that the natural space for our formulation is W = H ( D ) + ⊕ L ( D ) − .
3. Weak formulation of the slab problem
We follow the steps presented in [7] for multi-dimensional problems. The key idea is to rewrite the slabproblem into a weak formulation for the even and odd parts of the solution. Multiplication of (1) with atest function ψ ∈ W + and using orthogonality of even and odd functions gives that( µ∂ z φ − , ψ + ) + (( σ t − σ s P ) φ + , ψ + ) = ( q + , ψ + ) . Integration-by-parts (3) applied to the first term on the left-hand side yields that( µ∂ z φ − , ψ + ) = − ( φ − , µ∂ z ψ + ) + (cid:104) φ − , ψ + (cid:105) L − (cid:104) φ − , ψ + (cid:105) L − . Due to symmetries, we have that (cid:104) φ − , ψ + (cid:105) L = −(cid:104) φ − , ψ + (cid:105) L − . Using (2), we have that φ − = φ − φ + = g − φ + on the inflow boundary, which leads to( µ∂ z φ − , ψ + ) = − ( φ − , µ∂ z ψ + ) + 2 (cid:104) φ + , ψ + (cid:105) L − − (cid:104) g, ψ + (cid:105) L − . Thus, for any ψ + ∈ W + , it holds that2 (cid:104) φ + , ψ + (cid:105) L − − ( φ − , µ∂ z ψ + ) + (( σ t − σ s P ) φ + , ψ + ) = ( q + , ψ + ) + 2 (cid:104) g − , ψ + (cid:105) L − . (5)Testing (1) with an odd test function ψ − ∈ W − , we obtain that( µ∂ z φ + , ψ − ) + ( σ t φ − , ψ − ) = ( q − , ψ − ) , which implies that φ − = σ t ( q − − µ∂ z φ + ) ∈ W − . Using this expression for φ − in (5), we deduce that u = φ + is a solution to the following problem; cf. [7]. Problem 3.1.
Let q ∈ L ( D ) and g ∈ L − . Find u ∈ W + such that a ( u, v ) = (cid:96) ( v ) for all v ∈ W + . (6) where the bilinear form a : W + × W + → R is given by a ( u, v ) = 2 (cid:104) u, v (cid:105) L − + ( 1 σ t µ∂ z u, µ∂ z v ) + (( σ t − σ s P ) u, v ) , (7) and the linear form (cid:96) : W + → R is defined as (cid:96) ( v ) = ( q + , v ) + 2 (cid:104) g − , ψ + (cid:105) L − + ( q − , σ t µ∂ z v ) . (8)4 .2. Well-posedness We endow W + with the norm induced by the bilinear form a defined (7), i.e., (cid:107) u (cid:107) W = (cid:107) u (cid:107) a = a ( u, u ) for u ∈ W + . (9)Using the Cauchy-Schwarz inequality, we obtain the following result. Lemma 3.2.
The linear form (cid:96) : W → R defined in (8) is bounded, i.e., for all v ∈ W + it holds (cid:96) ( v ) ≤ ( (cid:107) q + (cid:107) σa + (cid:107) q − (cid:107) σt + 2 (cid:107) g − (cid:107) L − ) (cid:107) v (cid:107) a . Since the space W + endowed with the inner product induced by a is a Hilbert space, the unique solvabilityof Problem 3.1 is a direct consequence of the Riesz representation theorem. Theorem 3.3.
Let Assumption (A1) hold true. Then Problem 3.1 has a unique solution u ∈ W + . Moreover,we have the bound (cid:107) u (cid:107) a ≤ ( (cid:107) q + (cid:107) σa + (cid:107) q − (cid:107) σt + 2 (cid:107) g − (cid:107) L − ) . Remark 3.4.
Setting φ + = u and φ − = σ t ( q − − µ∂ z u ) , one can show that µ∂ z φ − ∈ L ( D ) , i.e., φ = φ + + φ − ∈ H ( D ) satisfies (1) in L ( D ) , cf. [7]. Using the trace lemma 2.1 and partial-integration (3) , wefurther can show that φ satisfies the boundary conditions (2) in the sense of traces. Hence, Theorem 3.3,independently, leads to a well-posedness result as Lemma 2.2 for (1) - (2) .3.3. DSA preconditioned source iteration in second order form As a preparation for the numerical solution of the discrete systems that will be described below, letus discuss an iterative scheme in infinite dimensions for solving the radiative transfer equation. The basicidea is a standard one and consists of decoupling scattering and transport in order to compute successiveapproximation, viz., the source iteration [16, 17]. Next to the basic iteration, we describe a preconditionerwhich resembles diffusion synthetic acceleration (DSA) schemes using the notation of [17] or a KP schemeusing the terminology of [16].In order to formulate the method, we introduce the following bilinear forms k ( u, v ) = ( σ s P u, v ) and b ( u, v ) = a ( u, v ) + k ( u, v ) for u, v, ∈ W + , and denote the induced semi-norm and norm by (cid:107) u (cid:107) k and (cid:107) u (cid:107) b , respectively.The iteration scheme is defined as follows: For u k ∈ W + given, compute u k + ∈ W + as the uniquesolution to b ( u k + , v ) = k ( u k , v ) + (cid:96) ( v ) for all v ∈ W + . (10)The half-step error e k + = u − u k + satisfies a ( e k + , v ) = k ( u k + − u k , v ) for all v ∈ W + . (11)The key idea is then to construct an easy-to-compute approximation to e k + by Galerkin projection onto asuitable subspace W +1 ⊂ W + . This approximation is then used to correct u k + to obtain a more accurateapproximation u k +1 to u . Define the following closed subspace of W + W +1 = { u ∈ W + : u = P u } , (12)i.e., W +1 consists of functions in W + that do not depend on µ . The correction u k + D ∈ W +1 is then computedby Galerkin projection of (11) to W +1 : a ( u k + D , v ) = k ( u k + − u k , v ) for all v ∈ W +1 , (13)5nd the new iterate is defined as u k +1 = u k + + u k + D . (14)If u k + D is a good approximation to e k + , then e k +1 = e k + − u k + D is small. The convergence proof of theiteration W + → W + , u k (cid:55)→ u k +1 is based on spectral analysis, see, e.g., [25, 26]. Lemma 3.5.
The half-step error e k + = u − u k +1 / of the iteration (10) satisfies (cid:107) e k + (cid:107) a ≤ c (cid:107) e k (cid:107) a , with constant c = (cid:107) σ s /σ t (cid:107) ∞ .Proof. We endow W + with the inner product induced by the bilinear b in this proof, and define bounded,self-adjoint and positive operators A and K on W + by b ( Au, v ) = a ( u, v ) , b ( Ku, v ) = k ( u, v ) u, v ∈ W + . Using a = b − k , we obtain that A = I − K . Using u k +1 / − u k = e k − e k +1 / , (11) can be written as e k +1 / = Ke k . We thus have that (cid:107) e k +1 / (cid:107) a = b (( I − K ) Ke k , Ke k ) = b ( K ( I − K ) / e k , ( I − K ) / e k ) ≤ max σ ( K ) (cid:107) ( I − K ) / e k (cid:107) b ≤ c (cid:107) e k (cid:107) a . In the last step we have used the following bounds on the numerical range0 ≤ b ( Kv, v ) = k ( v, v ) ≤ (cid:107) σ s σ t (cid:107) ∞ ( σ t v, v ) ≤ cb ( v, v ) , which yields the spectral bounds σ ( K ) ⊂ [0 , c ]. Theorem 3.6.
Let Assumption (A1) hold, and let c = (cid:107) σ s /σ t (cid:107) ∞ < be as in Lemma 3.5. For any u ∈ W + , the iteration defined by (10) , (13) , (14) converges linearly to the solution u of Problem 3.1 with (cid:107) u − u k +1 (cid:107) a ≤ c (cid:107) u − u k (cid:107) a . Proof.
Since u k +1 D is the orthogonal projection of e k + to W +1 in the a -inner product it holds that (cid:107) e k +1 (cid:107) a = (cid:107) e k + − u k +1 D (cid:107) a = inf v ∈ W +1 (cid:107) e k + − v (cid:107) a ≤ (cid:107) e k + (cid:107) a . The assertion then follows from Lemma 3.5.
Remark 3.7.
The convergence analysis presented in this section carries over verbatim to multi-dimensionalproblems without symmetries, and it can be extended immediately to more general (symmetric and positive)scattering operators. Moreover, if a Poincar´e-Friedrichs inequality is available, cf., e.g., [6], then the case σ a ≥ can be treated similarly as long as σ t is uniformly bounded away from , and the source iterationconverges also in this situation. Remark 3.8.
Problem (13) is the weak formulation of the diffusion equation − ∂ z ( 13 σ t ∂ z u D ) + σ a u D = f in (0 , Z ) , with f = σ s P ( u k + − u k ) , complemented by Robin boundary conditions, which shows the close relationshipto DSA schemes. emark 3.9. The convergence analysis for the iteration without preconditioning, can alternatively be basedon the following estimates. These estimates are the only ones in this paper that exploit that the scatteringoperator is related to the L -projector P . First note that (cid:107) e k + (cid:107) b = (cid:107)P e k + (cid:107) σ t + (cid:107) ( I − P ) e k + (cid:107) σ t + (cid:107) e k + (cid:107) L − + (cid:107) µ∂ z e k + (cid:107) σt and that (cid:107) e k + (cid:107) b = k ( e k , e k + ) . Setting c = (cid:107) σ s /σ t (cid:107) ∞ , the Cauchy-Schwarz inequality yields (1 − ε (cid:107)P e k + (cid:107) σ t + (cid:107) ( I − P ) e k + (cid:107) σ t + (cid:107) e k + (cid:107) L − + (cid:107) µ∂ z e k + (cid:107) σt ≤ c ε (cid:107)P e k (cid:107) σ t for any ε ∈ (0 , . Choosing ε = 2 shows that parts of the error are smoothened independently of c , i.e., (cid:107) ( I − P ) e k + (cid:107) σ t + (cid:107) e k + (cid:107) L − + (cid:107) µ∂ z e k + (cid:107) σt ≤ c (cid:107)P e k (cid:107) σ t , while the angular average is hardly damped. In any case, this shows that P e k converges to zero. It remainsopen how to exploit such an estimate to improve the analysis of the DSA preconditioned scheme above as itseems difficult to relate the latter smoothing property to the best approximation error of e k + in the a -norm.
4. Galerkin approximations
In this section, we construct conforming approximation spaces W + h,N ⊂ W + in a two-step procedure. Ina first step, we discretize the µ -variable using discontinuous ansatz functions. In a second step, we discretizethe z -variable by continuous finite elements. Before stating the particular approximation space, we providesome general results. Let us begin with the definition of the discrete problem. Problem 4.1.
Let q ∈ L ( D ) , g ∈ L − , W + h,N ⊂ W + and let a and (cid:96) be defined as in Problem 3.1. Find u h ∈ W + h,N such that for all v ∈ W + h,N there holds a ( u h , v h ) = (cid:96) ( v h ) . Since the bilinear form a induces the energy norm that we use in our analysis, we immediately obtainthe following best approximation result. Theorem 4.2.
Let W + h,N be a closed subspace of W + . Then, there exists a unique solution u h ∈ W + h,N ofProblem 4.1 that satisfies the a-priori estimate (cid:107) u h (cid:107) a ≤ ( (cid:107) q + (cid:107) σa + (cid:107) q − (cid:107) σt + 2 (cid:107) g − (cid:107) L − ) , (15) and the following best approximation error estimate (cid:107) u − u h (cid:107) a ≤ inf v h ∈ W + h,N (cid:107) u − v h (cid:107) a . (16)In the next sections, we discuss some particular discretizations. We note that these generalize thespherical harmonics approach presented in [7]. hp semidiscretization in µ Since we consider even functions, we require that the partition of the interval [ − ,
1] for the µ variablerespects the point symmetry µ (cid:55)→ − µ . For simplicity, we thus partition the interval [0 ,
1] only, and thepartition of [ − ,
0] is implicitly defined by reflection. 7et N ∈ N be a positive integer, and define intervals ¯ µ n = ( µ n − , µ n + ), n = 1 , . . . , N , such that µ = 0and µ N + = 1, and set ∆ µ n = µ n + − µ n − and µ n = ( µ n + + µ n − ) /
2. Denote χ n ( µ ) the characteristicfunction of the interval ¯ µ n . For µ >
0, we define the piecewise functions Q n,l ( µ ) = (cid:114) l + 12 P l (2 µ − µ n − ∆ µ n − χ ¯ µ n ( µ ) , µ > , where P l denotes the l th Legendre polynomial. Hence, { Q n,l } Ll =0 is an L -orthonormal basis for the spaceof polynomials of degree L on each interval ¯ µ n . For µ >
0, we set Q ± n,l ( µ ) = Q n,l ( µ ), and for µ <
0, we set Q ± n,l ( µ ) = ± Q ± n,l ( − µ ). The semidiscretization of the even parts is then u ( z, µ ) ≈ u h ( z, µ ) = N (cid:88) n =1 L (cid:88) l =0 φ + n,l ( z ) Q + n,l ( µ ) . Remark 4.3.
If we partition the interval [ − , for the angular variable by a single element, we obtaintruncated spherical harmonics approximations, see, e.g., [4, 7]. Partitioning of [ − , into two symmetricintervals ( − , ∪ (0 , corresponds to the double P L method [4], which generalizes in multiple dimensionsto half space moment methods [27]. The latter can resolve the non-smoothness of φ at µ = 0 , and, thusmight yield spectral convergence on the intervals µ > and µ < .4.2. Fully discrete scheme In order to obtain a conforming discretization, we approximate the coefficient functions φ + n,l using H (0 , Z )-conforming elements. Let J ∈ N , and ¯ z j = ( z j − , z j ) such that[0 , Z ] = ∪ Jj =1 clos(¯ z j )be a partition of (0 , Z ). Let p ≥ P p the space of polynomials of degree p . The full approximationspace is then defined by W + h,N = { ψ + h ( z, µ ) = N (cid:88) n =1 L (cid:88) l =0 ψ + n,l ( z ) Q + n,l ( µ ) : ψ + n,l ( z ) ∈ H (0 , Z ) , ψ + n,l | ¯ z j ∈ P p } . (17)The choice (17) for the approximation space W + h,N corresponds to an hp finite element method, for whichthe assertion of Theorem 4.2 holds true. Remark 4.4.
If the solution u h ∈ W + h,N to Problem 4.1 is computed, we compute even and odd approxi-mations to the solution φ of (1) via φ + h = u h and the odd part φ − h ∈ W − h as the solution to the variationalproblem ( φ − h , ψ − h ) = ( σ t ( q − − µ∂ z u h ) , ψ − h ) for all ψ − h ∈ W − h , where W − h,N = { ψ − h ( z, µ ) = N (cid:88) n =1 L +1 (cid:88) l =0 ψ − n,l ( z ) Q − n,l ( µ ) : ψ − n,l | ¯ z j ∈ P p − } . We note that this space satisfies the compatibility condition µ∂ z W + h,N ⊂ W − h,N , which makes this pair ofapproximation spaces suitable for a direct approximation of a corresponding mixed formulation, cf. [7].The even-parity formulation that we consider here corresponds then to the Schur complement of the mixedproblem, cf. [7]. The reader should note the different degrees in the polynomial approximations, e.g., if φ + h is piecewise constant in angle, then φ − h is piecewise linear. . Discrete preconditioned source iteration In order to solve the discrete variational problem defined in Problem 4.1, we proceed as in Section 3.3but with W + and W +1 replaced by W + h,N and W + h, , respectively. We note that W + h, ⊂ W +1 consists offunctions in W + h,N that do not dependent on µ .The finite dimensional counterpart of the DSA preconditioned source iteration is then defined as follows:For given u kh ∈ W + h,N , compute u k + h ∈ W + h,N as the unique solution to b ( u k + h , v h ) = k ( u kh , v h ) + (cid:96) ( v h ) for all v h ∈ W + h,N . (18)The correction u k + h,D ∈ W + h, is defined by Galerkin projection of e k + h to W + h, : a ( u k + h,D , v h ) = k ( u k + h − u kh , v h ) for all v h ∈ W + h, , (19)and the new iterate is defined as u k +1 h = u k + h + u k + h,D . (20)Using the same arguments as above, we obtain the following convergence result. Theorem 5.1.
Let Assumption (A1) hold, and let c < be as in Lemma 3.5. For any u h ∈ W + h,N , theiteration defined by (18) , (19) , (20) converges linearly with (cid:107) u h − u k +1 h (cid:107) a ≤ c (cid:107) u h − u kh (cid:107) a . Remark 5.2.
Similar to Remark 3.8, u k + h,D is the Galerkin projection to W + h, of the solution to − ∂ z ( D ( z ) ∂ z u D ) + σ a u D = f, < z < Z, with µ -grid dependent diffusion coefficient D ( z ) and f = σ s P ( u k + h − u kh ) . If piecewise constant functionsin angle are employed, then D ( z ) = σ t ( z ) (1 + (cid:80) Nn =1 ∆ µ n ) . Remark 5.3.
Once the scattering term in the right-hand side of (18) has been computed, the half-step iterate u k + h can be computed independently for each direction, and, thus, its computation can be parallelized.
6. Numerical examples
In this section, we report on the accuracy of the proposed discretization scheme and its efficient numericalsolution using the DSA preconditioned source iteration of Section 5. We restrict our discussion to a low-order method that offers local resolution. To do so, we fix p = 1 and L = 0 in the definition of W + h,N ,while N might be large. Hence, the approximation space W + h,N consists of discontinuous, piecewise constantfunctions in the angular variable and continuous, piecewise linear functions in z . To investigate the convergence behavior, we use manufactured solutions, i.e., the exact solution is φ ( z, µ ) = | µ | e − µ e − z (1 − z ) , (21)with parameters σ a = 1 / σ s ( z ) = 2 + sin( πz ) / Z = 1, and source terms defined accordingly. Wecomputed the numerical solution u h using the DSA preconditioned iteration (18), (19), (20). We stoppedthe iteration using the a-posteriori stopping rule (cid:107) u kh − u k − h (cid:107) a ≤ ε, ε = 10 − is a chosen tolerance. The approximations of the even and odd parts of the solutions arerecovered as described in Remark 4.4.For the chosen parameters, we have that κ ≤ . N and J . As expected we observe a linearrate of convergence with respect to N , cf. Theorem 4.2. In addition, we observe that the preconditionedsource iteration converged within at most 15 iterations, and the expression (cid:107) u kh − u k − h (cid:107) a decreased by 0 .
21 ineach iteration. The convergence rate with respect to J is initially quadratic, which is better than predictedby Theorem 4.2, and then saturates for fixed N = 8192. As before, the preconditioned source iterationconverged within at most 15 iterations with a decrease by a factor of 0 .
21 of (cid:107) u kh − u k − h (cid:107) a .The observed convergence rate is thus as the one obtained by classical Fourier analysis for constantcoefficients and periodic boundary conditions, which is bounded by 0 . N E h rate512 1.61e-041024 8.07e-05 0.992048 4.04e-05 0.994096 2.04e-05 0.998192 1.06e-05 0.95 J E h rate16 7.88e-0432 1.99e-04 1.9864 5.14e-05 1.96128 1.63e-05 1.66256 1.06e-05 0.62 Table 1: Observed errors E h = (cid:107) φ − φ h (cid:107) L ( D ) between finite element solution φ h and the manufactured solution φ defined in(21) together with the rate of convergence of E h . Left: Convergence for different discretization parameters N , and J = 256.Right: Convergence for different discretization parameters J and N = 8192. In the next section, we present a more detailed study of spectrum of the preconditioned operator.
In this section, we consider the spectrum of the error propagation operator P e nh (cid:55)→ P e n +1 h , where P isthe L -projection onto constants in angle defined in Section 2 and e nh is the sequence of errors generatedby the DSA preconditioned source iteration, cf. Remark 3.9. The function P e nh depends only on z , and,assuming σ s >
0, we can measure the projected error using the norm induced by σ s , i.e., (cid:107)P e nh (cid:107) σ s = (cid:107) e nh (cid:107) k . We choose the following scattering and absorption parameters σ s ( z ) = (cid:40) πz ) , z ≤
102 + sin(2 πz ) , z > , σ a ( z ) = (cid:40) . , z ≤ . , z > . We notice that both parameters have huge jumps and that the predicted convergence rate is c = (cid:107) σ s /σ t (cid:107) ∞ ≈ . . N . The spectra for N ≥ . Although the theory has been presented for slab problems, our results carry over verbatim to otherproblems. In the following we report on the performance of the method for the lattice problem, cf. [21].10 igure 1: Spectra of the error propagation operator P e nh (cid:55)→ P e n +1 h for different spatial discretizations J = 16 , ,
512 (fromleft to right). Each plot contains the corresponding spectra for N = 2 i , i = 1 , . . . , N = 4 and N = 64 triangles. Right: Geometry of the latticeproblem. This test problem is a three-dimensional problem with certain symmetries. The radiative transfer equationhere writes as s · ∇ φ ( x, s ) + σ t ( x ) φ ( x, s ) = σ s ( x ) P φ + q ( x, s )where x ∈ R ⊂ R and s ∈ S ⊂ R . The computational domain is a square R = (0 , × (0 , σ a = 10 in the black regions shown in Figure 2, and σ a = 0 elsewhere. We set σ s = 1 in the grey and white regions and σ s = 0 elsewhere. The source is defined by q ( x, s ) = 1 in the whiteregion and q ( x, s ) = 0 elsewhere. Note that due to the availability of a Poincar´e-Friedrichs inequality [6], thecase σ a = 0 leads to a well-posed radiative transfer problem, and, since σ s + σ a ≥
1, the theory presented hereis applicable. Moreover, the constant c will depend on the constant from the Poincar´e-Friedrichs inequalityand c < σ a vanishes.We discretize L ( S ) by approximating S using a geodesic polyhedron consisting of flat triangles, seeFigure 2. The approximation space in angle then consists of standard discontinuous finite element spacesassociated to this triangulation. In the following, we focus on piecewise constant approximations, but higherorder approximations are straight-forward if the geometry approximation is also of higher order.In Figure 3 we show the angular averages of the computed solutions for two different grids with J = 9 801vertices in the spatial grid and N = 4 triangles on a half-sphere and for J = 78 961 and N = 64, respectively.We note that our solutions do not exhibit ray effects, cf. [4, 21]. The preconditioned source iterationconverged with 9 and 17 iterations for the coarse grid approximation and for the fine grid approximation,respectively. This amounts to an error reduction per iteration of at least 0 .
04 for the coarse grid discretizationand of at least 0 .
20 for the fine grid discretization.Next, let us investigate the behavior of the preconditioned source iteration for scaled parameters in more11 igure 3: Angular average of the computed solution in a log -scale for the lattice problem for J = 9 801 spatial vertices and N = 4 triangles on a half-sphere (left) and J = 78 961 spatial vertices and N = 64 triangles on a half-sphere (right). detail. Introducing a scale parameter δ >
0, a diffusion limit is obtained by the scaling [28]¯ σ s ( x ) δ , δ ¯ σ a ( x ) , δ ¯ q ( x ) , where both parameters ¯ σ s and ¯ σ a are bounded and strictly bounded away from zero. Since this is not thecase for the lattice problem, we consider the following parameters σ δs ( x ) = σ s ( x ) + 1 / δ , σ δa = δ ( σ a ( x ) + 1 / , q δ ( x, s ) = δq ( x ) . For δ →
0, the corresponding solution u δ will converge to the solution of a diffusion problem; for non-smooth coefficients see [20]. The parameter c defined in Lemma 3.5 is bounded by O (1 /δ ). Table 2 shows J = 9 801 J = 78 961 N = 4 N = 64 N = 4 N = 641 /δ k rate k rate k rate k rate1 9 0.04 15 0.16 9 0.04 15 0.1710 9 0.06 15 0.22 9 0.06 16 0.25100 8 0.06 13 0.22 9 0.07 15 0.271000 5 0.01 7 0.06 6 0.05 10 0.17 Table 2: Iteration counts k and minimal reduction rates for (cid:107) u kh − u k − h (cid:107) a for the lattice problem with scaled parameters σ δs , σ δa and q δ for different δ and discretizations with N triangles on a half-sphere and J vertices in the spatial mesh. the iteration counts and the minimum reduction of (cid:107) u kh − u k − h (cid:107) a during the iteration. We observe that thepreconditioned iteration is robust with respect to δ →
7. Conclusions
We investigated discontinuous angular and continuous spatial approximations of the even-parity for-mulation for the radiative transfer equation. Certain instances of these approximations are closely related12o classical discretizations, such as truncated spherical harmonics approximations, double P L -methods ordiscrete ordinates methods.We considered a diffusion accelerated preconditioned source iteration for the solution of the resultingvariational problems that has been formulated in infinite dimensions. Convergence rates of this iterationhave been proven. The Galerkin approach used for the discretization allowed to translate the results for theinfinite dimensional iteration directly to the discrete problems. Moreover, the discrete iteration convergesindependently of the chosen resolution.The theoretically proven convergence rate in Theorem 3.6 is not robust in the limit of large scattering,while numerical results show that in practice the preconditioned iteration converges robustly even in scat-tering dominated problems. One approach to obtain an improved convergence rates estimate is to estimatethe best approximation error inf v ∈ W +1 (cid:107) e k + − v (cid:107) a , which, however, seems rather difficult, and we postponea corresponding rigorous analysis to future research.The term diffusion acceleration is linked to the usage of the space W +1 in (12) that consists of functionsconstant in angle. This choice is, however, not essential, and other closed subspace can be employed, whichallows for the construction of multi-level schemes. The multi-level approach might lead to a feasible approachto estimate the best approximation error.In any case, the DSA preconditioner can be combined with a conjugate gradient method in order tofurther reduce the number of iterations. Moreover, unlike many discrete ordinates schemes, the numericalapproximation did not suffer from the ray effect in our numerical examples. Acknowledgments
The authors would like to thank Herbert Egger and Markus Bachmayr for stimulating discussions.
References [1] S. Chandrasekhar, Radiative Transfer, Dover Publications, Inc., 1960 (1960).[2] K. M. Case, P. F. Zweifel, Linear transport theory, Addison-Wesley, Reading, 1967 (1967).[3] V. S. Vladimirov, Mathematical problems in the one-velocity theory of particle transport, Tech. rep., Atomic Energy ofCanada Ltd. AECL-1661. translated from Transactions of the V.A. Steklov Mathematical Institute (61) (1961).[4] E. E. Lewis, W. F. Miller Jr., Computational Methods of Neutron Transport, John Wiley & Sons, Inc., New YorkChichester Brisbane Toronto Singapore, 1984 (1984).[5] J. Pitk¨aranta, Approximate solution of the transport equation by methods of Galerkin type, J. Math. Anal. and Appl. 60(1977) 186–210 (1977).[6] T. A. Manteuffel, K. J. Ressel, G. Starke, A boundary functional for the least-squares finite-element solution for neutrontransport problems, SIAM J. Numer. Anal. 2 (2000) 556–586 (2000).[7] H. Egger, M. Schlottbom, A mixed variational framework for the radiative transfer equation, Math. Mod. Meth. Appl.Sci. 22 (2012) 1150014 (2012).[8] H. Egger, M. Schlottbom, A perfectly matched layer approach for radiative transfer in highly scattering regimes., SIAMJ. Numer. Anal. (accepted) (2019).[9] C. Johnson, J. Pitk¨aranta, Convergence of a fully discrete scheme for two-dimensional neutron transport, SIAM J. Numer.Anal. 20 (1983) 951–966 (1983).[10] M. Asadzadeh, Analysis of a fully discrete scheme for neutron transport in two-dimensional geometry, SIAM J. Numer.Anal. 23 (1986) 543–561 (1986).[11] W. H. Reed, T. R. Hill, Triangular mesh methods for the neutron transport equation, Tech. rep., Los Alamos ScientificLaboratory of the University of California (1973).[12] T. Wareing, J. McGhee, J. Morel, S. Pautz, Discontinuous finite element sn methods on three-dimensional unstructuredgrids, Nuclear Science and Engineering 138 (2001) 1—13 (2001).[13] J. C. Ragusa, J.-L. Guermond, G. Kanschat, A robust S N -DG-approximation for radiation transport in optically thickand diffusive regimes, J. Comput. Phys. 231 (4) (2012) 1947–1962 (2012). doi:10.1016/j.jcp.2011.11.017 .[14] G. Kanschat, J.-C. Ragusa, A robust multigrid preconditioner for S N DG approximation of monochromatic, isotropicradiation transport problems, SIAM J. Sci. Comput. 36 (5) (2014) A2326–A2345 (2014). doi:10.1137/13091600X .[15] J. K´oph´azi, D. Lathouwers, A space-angle DGFEM approach for the Boltzmann radiation transport equation with localangular refinement, J. Comput. Phys. 297 (2015) 637–668 (2015). doi:10.1016/j.jcp.2015.05.031 .[16] G. I. Marchuk, V. I. Lebedev, Numerical Methods in the Theory of Neutron Transport, Harwood Academic Publishers,Chur, London, Paris, New York, 1986 (1986).[17] M. L. Adams, E. W. Larsen, Fast iterative methods for discrete-ordinates particle transport calculations, Progress inNuclear Energy 40 (1) (2002) 3–159 (2002).
18] J. S. Warsa, T. A. Wareing, J. E. Morel, Krylov iterative methods and the degraded effectiveness of diffusion synthetic ac-celeration for multidimensional sn calculations in problems with material discontinuities, Nuclear Science and Engineering147 (3) (2004) 218–248 (2004). doi:10.13182/NSE02-14 .[19] E. W. Larsen, J. B. Keller, Asymptotic solution of neutron transport problems for small mean free paths, J. MathematicalPhys. 15 (1974) 75–81 (1974). doi:10.1063/1.1666510 .[20] H. Egger, M. Schlottbom, Diffusion asymptotics for linear transport with low regularity, Asymptot. Anal. 89 (3-4) (2014)365–377 (2014).[21] T. A. Brunner, Forms of approximate radiation transport, in: Nuclear Mathematical and Computational Sciences: ACentury in Review, A Century Anew Gatlinburg, American Nuclear Society, LaGrange Park, IL, 2003, tennessee, April6-11, 2003 (2003).[22] V. Agoshkov, Boundary Value Problems for Transport Equations, Modeling and Simulation in Science, Engineering andTechnology, Birkh¨auser, Boston, 1998 (1998).[23] H. Egger, M. Schlottbom, An L p theory for stationary radiative transfer, Appl. Anal. 93 (6) (2014) 1283–1296 (2014). doi:10.1080/00036811.2013.826798 .[24] J. C. Blake, Domain decomposition methods for nuclear reactor modelling with diffusion acceleration, Ph.D. thesis,University of Bath (2016).[25] G. Helmberg, Introduction to spectral theory in Hilbert space, North-Holland Series in Applied Mathematics and Me-chanics, Vol. 6, North-Holland Publishing Co., Amsterdam-London; Wiley Interscience Division John Wiley & Sons, Inc.,New York, 1969 (1969).URL https://mathscinet.ams.org/mathscinet-getitem?mr=0243367 [26] D. Werner, Funktionalanalysis, extended Edition, Springer-Verlag, Berlin, 2000 (2000).URL https://mathscinet.ams.org/mathscinet-getitem?mr=1787146 [27] B. Dubroca, M. Frank, A. Klar, G. Th¨ommes, A half space moment approximation to the radiative heat transfer equations,ZAMM Z. Angew. Math. Mech. 83 (12) (2003) 853–858 (2003). doi:10.1002/zamm.200310055 .[28] R. Dautray, J. L. Lions, Mathematical Analysis and Numerical Methods for Science and Technology, Evolution ProblemsII, Vol. 6, Springer, Berlin, 1993 (1993)..[28] R. Dautray, J. L. Lions, Mathematical Analysis and Numerical Methods for Science and Technology, Evolution ProblemsII, Vol. 6, Springer, Berlin, 1993 (1993).