 # Low-memory, discrete ordinates, discontinuous Galerkin methods for radiative transport

The discrete ordinates discontinuous Galerkin (S_N-DG) method is a well-established and practical approach for solving the radiative transport equation. In this paper, we study a low-memory variation of the upwind S_N-DG method. The proposed method uses a smaller finite element space that is constructed by coupling spatial unknowns across collocation angles, thereby yielding an approximation with fewer degrees of freedom than the standard method. Like the original S_N-DG method, the low memory variation still preserves the asymptotic diffusion limit and maintains the characteristic structure needed for mesh sweeping algorithms. While we observe second-order convergence in scattering dominated, diffusive regime, the low-memory method is in general only first-order accurate. To address this issue, we use upwind reconstruction to recover second-order accuracy. For both methods, numerical procedures based on upwind sweeps are proposed to reduce the system dimension in the underlying Krylov solver strategy.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Radiative transport equations [2, 9, 27, 25, 26, 12, 8] describe the flows of particles, such as photons, neutrons, and electrons, as they pass through and interact with a background medium. These equations are used in various applications, including astrophysics and nuclear reactor analysis.

In this paper, we consider the scaled, steady-state, linear transport equation equationparentequation

 Ω⋅∇Ψ(Ω,x)+(σs(x)ε+εσa(x))Ψ(Ω,x) =σs(x)ε¯¯¯¯Ψ(x)+εq(x), (Ω,x)∈S×D, (1.1a) Ψ(Ω,x) =α(Ω,x), (Ω,x)∈Γ−. (1.1b)

Here is an open, bounded, and Lipschitz domain; is the projection of the unit sphere in into (the interval for and unit disk for ); and , where

is the outward unit normal vector at any point

where the boundary is .

The angular flux is the flux of particles at the location moving with unit speed in the direction , and the scalar flux is the average of over .111Often the quantity is referred to as the scalar flux. The difference is simply a normalization factor from integration of the sphere. Here, we borrow the convention used in . The functions and are (known) non-dimensionalized scattering and absorption cross-sections, respectively, and is a (known) non-dimensionalized source. The function is the (known) incoming flux at moving in the direction . The constant is a scaling parameter which characterizes the relative strength of scattering.

Designing effective numerical methods for (1.1) is a serious challenge, and the intent of this paper is to address two of the main issues. Firstly, for a three-dimensional problem, the unknown intensity is a function of three spatial and two angular variables; the discretization of this five-dimensional phase space usually requires significant computational resources. Secondly, when the parameter is small, is nearly independent of and can be approximated by the solution of a diffusion equation in the variable only [16, 5, 6]. That is, away from the boundary, as , where satisfies

 −∇⋅(13σs∇Ψ(0)(x))+σaΨ(0)(x)=q(x),x∈D, (1.2)

along with appropriate boundary conditions. A numerical method for (1.1) should preserve this asymptotic limit without having to resolve the length scales associated with . In other words, in the limit , a discretization of the transport equation (1.1) should become a consistent and stable discretization of the diffusion equation (1.2). Otherwise a highly refined mesh is needed to approximate the solution accurately .222This issue is also known as “locking” in the elliptic literature .

Classical approaches for discretizing (1.1) often involve separate treatment of the angular and spatial variables, and a variety of options are available. Among them, the -DG method [17, 21, 1] has received significant attention due to it’s robustness, computational efficiency, and convenient implementation. The method (see for a substantial review and additional references) is a collocation method in which the angular variable is discretized into a finite number of directions and a quadrature rule is used to evaluate . The discretization preserves non-negativity of and can incorporate the boundary conditions from (1.1) in a straightforward way. It also preserves the characteristic structure of the advection operator in (1.1), which allows for the use of fast sweeping techniques for inverting the discrete form of the operator on the left-hand side of (1.1).

Discontinuous Galerkin (DG) methods are a class of finite element methods that construct numerical solutions using piecewise polynomial spaces. The DG approach was introduced in  for the express purpose of solving equations like (1.1), followed shortly thereafter by a rigorous analysis in . Since then, DG methods have been applied to nonlinear hyperbolic conservation laws and convection-dominated problems , elliptic problems , and equations with higher-order derivatives [32, 31]. When used with upwind fluxes, DG methods preserve the characteristic structure of (1.1) that enables sweeps. Moreover, if the approximation space can support globally continuous linear polynomials, then DG methods with upwind fluxes will yield accurate numerical solutions for without the need to resolve with the spatial mesh [21, 1, 13]. However, this condition on the approximation space means that at least elements must be used for a triangular mesh and elements for a rectangular mesh.333This condition can be circumvented for non-upwind methods. In , the authors made the piecewise constant DG method asymptotic preserving with parameters adjusting numerical fluxes under different regimes. Similar techniques were introduced in finite volume contexts  as well and were recently used in  to develop a positive, asymptotic preserving method.

In order to reduce memory costs in the upwind -DG method, while still preserving the asymptotic diffusion limit and maintaining the characteristic structure needed for sweeps, we propose in this paper to couple the finite element spaces between different collocation angles in the discrete ordinate approximation. Since the solution becomes isotropic in the diffusion limit (), we hypothesize that only a (for triangles) or

(for rectangles) approximation of the angular average is necessary. Thus, instead of using a tensor product finite element space for the

-DG system, we seek the solution in a proper subspace, in which all the elements have isotropic slopes. This choice of finite element space yields a significant reduction in memory per spatial cell, as illustrated in Table 1.1.

In the diffusion limit, the low-memory approach typically displays second-order accuracy. However, because the finite element representation of each ordinate is coupled to all the other ordinates, the overall accuracy of the low-memory approach for fixed is only first-order. To address this drawback, we propose a modification of the low-memory scheme that uses local reconstruction to improve accuracy. As long as the reconstruction uses upwind information, the resulting transport operator can still be inverted with sweeps. While rigorous theoretic properties of this modified scheme are still under investigation, we observe numerically that it recovers second-order accuracy for arbitrary fixed and captures the asymptotic diffusion limit. However, the method does generate some small numerical artifacts at the discontiuity of the cross section, which we point out in the numerical results of Section 4.

The rest of the paper is organized as follows. In Section 2, we introduce the background and revisit the -DG method. Low-memory methods, including the original first-order approach and the second-order reconstructed scheme, are detailed in Section 3. Numerical tests are provided in Section 4 to illustrate the behavior of both approaches. Finally, conclusions and future work are discussed in Section 5.

## 2 The Sn-DG method

In this section, we review the -DG scheme and discuss its asymptotic properties and implementation. Throughout the paper, we consider the case and , unless otherwise stated. In general, the well-posedness of (1.1) also holds for . In some places, we will also assume that the cross-section is piecewise constant, either to simplify the exposition or to make connections between first- and second-order forms of the diffusion limit. In the numerics, we often consider nonzero boundary conditions. However, in proofs we often assume that . When is nonzero but isotropic, many of the results still hold. However, when is anisotropic, the diffusion equation requires a boundary layer correction in order to be uniformly accurate . At the discrete level, this situation requires more sophisticated analysis [21, 1, 13] than is presented here.

### 2.1 Formulation

Consider a quadrature rule with points and positive weights such that

 1|S|∫Sf(Ω)dΩ≈nΩ∑j=1wjf(Ωj),∀f∈C(S). (2.1)

We assume the quadrature is exact for polynomials in up to degree two444Level symmetric quadratures of moderate size will satisfy these properties. See, e.g.,  and references therein.; that is,

 (i) nΩ∑j=1wj=1,(ii) nΩ∑j=1wjΩj=0,and(iii) nΩ∑j=1wjΩj⊗Ωj=13Id. (2.2)

The method approximates the angular flux at the quadrature points by a vector-valued function whose components satisfy a coupled system with equations

 Ωj⋅∇ψj(x)+(σsε+εσa)ψj(x)=σsε¯¯¯¯ψ(x)+εq(x),¯¯¯¯ψ(x)=nΩ∑j=1wjψ(Ωj,x). (2.3)

To formulate the upwind DG discretization of the system (2.3), let be a quasi-uniform partition of the domain . We assume to avoid unnecessary technicalities. Let be the collection of cell interfaces and let be the collection of boundary faces. Given a cell , we denote by the outward normal on and for any , let and . Given a face , we denote by a prescribed normal (chosen by convention) and, for any , let . For convenience, we assume trace values are identically zero when evaluated outside of .

The standard -DG method uses the tensor-product finite element space

 Vh=nΩ∏j=1Vh,Vh={vj:vj|K∈Z1(K)}, (2.4)

where for triangular or tetrahedral meshes, is the space of linear polynomials on and for Cartesian meshes is the space of multilinear polynomials on . The space can be equipped with an inner product and associated norm given by

 (u,v)=∑K∈ThnΩ∑j=1wj∫Kujvjdxand∥v∥=√(v,v). (2.5)

The semi-norm induced by jumps at the cell interfaces is given by

 ⟦v⟧=⎛⎝∑F∈FhnΩ∑j=1wj∫F|Ωj⋅νF|(v−j−v+j)2dx⎞⎠1/2. (2.6)

To construct the -DG method, define the local operators equationparentequation

 Lj,K(u,v)= −∫KujΩj⋅∇vjdx+∫∂KˆujΩj⋅νKvintjdx (2.7a) +∫K(σsε+εσa)ujvjdx, Sj,K(u,v)= ∫Kσsε¯¯¯uvjdx,with ¯¯¯u=nΩ∑j=1wjuj, (2.7b) Qj,K,α(v)= ∫Kεqvjdx−∫∂K∩F∂hαΩj⋅νKvintjdx, (2.7c)

where is the upwind trace at , and is defined as zero when the limit is taken outside of . Then set

 B(u,v)=L(u,v)−S(u,v), (2.8)

where

 L(u,v)=∑K∈ThnΩ∑j=1wjLj,K(u,v)andS(u,v)=∑K∈ThnΩ∑j=1wjSj,K(u,v), (2.9)

and let

 Qα(v)=∑K∈ThnΩ∑j=1wjQj,K,α(v). (2.10)

The -DG method is then: find such that

 B(ψh,v)=Qα(v),∀v∈Vh. (2.11)

#### 2.1.1 Implementation

Recall that is the number of discrete ordinates in the discretization. Let be the number of mesh cells in and let be the dimension of . Then the dimension of is .

Let be a set of basis functions for , with locally supported on . Then the set , where () and is the Kronecker delta, gives a complete set of basis functions for . With this choice of basis functions, the variational formulation in (2.11), written as

 L(ψh,v)=S(ψh,v)+Qα(v),∀v∈Vh, (2.12)

can be assembled into a linear system (detailed in Appendix A)

 LΨ=MPΨ+Q. (2.13)

In the above equation, is an block diagonal matrix, where the -th block () corresponds to the discretization of the operator ; is an injective matrix, is an matrix; is an vector assembled from the source and the inflow boundary ; and is an vector such that .

If upwind values are used to evaluate the numerical trace , each block of can be inverted efficiently with a sweep algorithm. The system in (2.13) can be solved numerically with a Krylov method by first solving the reduce system

 Φ−PL−1MΦ=PL−1Q (2.14)

for the vector . This equation is derived by applying and then to (2.13). In a second step is recovered from the relation

 Ψ=L−1MΦ+L−1Q. (2.15)

The following theorem is proven in Appendix B.

###### Theorem

The matrix is invertible.

[Sherman–Morrison formula] According to the Sherman-Morrison formula (see for example [11, Section 2.1.3]): given invertible matrices and ,

 B−1=A−1−A−1U(I+VA−1U)−1VA−1. (2.16)

The direct application of (2.16) with , and , yields the formula in (2.15) with given by (2.14).

#### 2.1.2 Asymptotic scheme

As , the -DG scheme gives a consistent approximation to the asymptotic diffusion problem. For simplicity, we focus here on the zero inflow boundary condition . The analysis of more general boundary conditions can be found in [1, 13, 14, 21].

We use an overline to represent isotropic subspaces. For example,

 ¯¯¯¯Vh={v=(v1,…,vnΩ)∈Vh:vi=¯¯¯v,∀i}. (2.17)

We further define to be the space of continuous functions in that vanish on . is used to represent the tensor product space of with an induced norm still denoted as . In particular, since and are isomorphic, we often identify with . To facilitate the discussion, we also define

 Jh=1εnΩ∑j=1wjΩjψh,j=nΩ∑j=1wjΩjψh,j−¯¯¯¯ψhε, (2.18)

which is a vector field in . The following result is proved in 555The result in  is actually stated for more generally. In particular it allows to be nonzero and possibly anisotropic.; see also  and Section 3.1 in this paper.

###### Theorem (Asymptotic scheme)

Suppose . Then as , and converge to and , respectively, that are the unique solution to the mixed problem: equationparentequation

 ∑K∈Th∫K(−J(0)h⋅∇φ+σaψ(0)hφ)dx=∫Dqφdx, (2.19a) ∑K∈Th∫K(13∇ψ(0)h+σsJ(0)h)⋅ζdx=0, (2.19b)

and .

## 3 Low-memory strategies

In this section, we generalize the statement of Section 2.1.2 slightly to allow for proper subspaces of in the finite element formulation. Based on the analysis, a first-order low-memory scheme is constructed. We then apply the reconstruction technique to lift the accuracy of the method to second-order.

### 3.1 Asymptotic schemes with subspaces of Vh

The results of Section 2.1.2 suggest that, rather than , it is the approximation of the integrated quantities and that play an important role in the diffusion limit. In particular, the continuity requirement on plays a crucial role. Indeed, as is well known , if the space is constructed from piecewise constants, then (2.19) implies that is a global constant and . This solution is clearly inconsistent with the diffusion limit. However, it is possible to construct a DG method: find such that

 B(ψh,v)=Qα(v),∀v∈Wh (3.1)

based on a proper subspace that maintains the diffusion limit, but requires fewer unknowns for a given mesh .

###### Theorem

For each and linear subspace , (3.16) has a unique solution. In particular, if

, the solution satisfies the energy estimate

 1ε∥σ12s(ψh−¯¯¯¯ψh)∥2+ε2∥σ12aψh∥2+12⟦ψh⟧2≤ε2δa∥q∥2. (3.2)

The proof is based on coercivity of and we refer to  and  for details. Here, is assumed for simplicity. Energy estimates with general inflow boundary condition can be found in [13, Lemma 4.2]. In , the case is studied and error estimates are derived using the coercivity with respect to a modified norm.

We next characterize sufficient conditions for . Define the spaces

 ¯¯¯¯¯¯¯¯¯¯¯¯ΩWh:={nΩ∑j=1wjΩjvj:v∈Wh}⊂¯¯¯¯VdhandΩ⋅¯¯¯¯¯¯¯¯¯¯¯¯ΩWh:={Ω⋅ζ:ζ∈¯¯¯¯¯¯¯¯¯¯¯¯ΩWh}⊂Vh, (3.3)

where . According to (2.18), . Section 2.1.2 can now be generalized to the space .

###### Theorem

Suppose . Suppose is a linear space such that . Then as , and converge to and , respectively, that are the unique solution to the mixed problem (2.19),

and .

###### Proof

Because the proof follows the arguments in [13, Section 4] closely, we provide only a brief outline, emphasizing where the condition on the space plays a role.

1. The stability estimate in (3.2) provides the following three bounds:

 (i) ∥ψh∥2≤1δ2a∥q∥2,(ii) ∥ψh−¯¯¯¯ψh∥2≤ε2δaδs∥q∥2,and(iii) ⟦ψh⟧2≤εδa∥q∥2. (3.4)

Bounds (i) and (ii) imply that converges (via a subsequence) to a function . Bound (iii) implies that .

2. Since, from the definition in (2.18),

 ∥Jh∥≤nΩ∑j=1wj∥ψh−¯¯¯¯ψh∥ε, (3.5)

where is the tensor product norm of in , the bound (ii) implies further that is uniformly bounded and hence converges subsequentially to a limit .

3. The equation in (2.19a) is derived by testing (3.1) with and using the fact that is independent of and continuous in .

4. It is the derivation of (2.19b) which uses the condition . Specifically, if with , then this condition implies that . Therefore, we can test (3.1) with this choice of to find that

 L(ψh,Ω⋅ζ)−S(ψh,Ω⋅ζ)= −nΩ∑j=1wj∑K∈Th∫Kψh,jΩj⋅∇(Ωj⋅ζ)dx (3.6) +nΩ∑j=1wj∑K∈Th∫∂Kˆψh,j(Ωj⋅νK)(Ωj⋅ζint)dx =: I+II+III.

We combine and , using the fact that and invoking (2.2). This gives

 limε→0(I+II) =nΩ∑j=1wj(Ωj⊗Ωj):∑K∈Th(−∫K¯¯¯¯ψ(0)h∇ζdx+∫∂K¯¯¯¯ψ(0)hνK⊗ζintdx) (3.7) =13Id:∑K∈Th∫K∇¯¯¯¯ψ(0)h⊗ζdx=∑K∈Th∫K13∇¯¯¯¯ψ(0)h⋅ζdx.

Since ,

 limε→0III=limε→0∑K∈Th∫K(σsε+εσa)nΩ∑j=1(wjψh,jΩj)⋅ζdx=∫KσsJ(0)h⋅ζdx. (3.8)

Finally, the right-hand side of (3.1) is (for )

 Q0(v)=nΩ∑j=1wj∑K∈Th∫KΩj⋅ζqdx=0. (3.9)

Combining (3.7), (3.8), and (3.9) recovers (2.19b).

5. Uniqueness of the subsequential limits and follows from the uni-solvency of (2.19). Indeed if is the difference between any two solutions of (2.19), then

 3σs∥˜Jh∥2+σa∥˜ψh∥2=0. (3.10)

Since and are assumed positive, it follows that and are identically zero.

We then discuss the choice of and the corresponding space pair, and , in the diffusion limit. Let be the space spanned by constants on . Then we define the piecewise constant space and its orthogonal complement The isotropic subspace of is denoted by and the subsequent product space is denoted by , .

1. When or , we have , which implies and .

2. When , it can be shown that , 666Since , which forces . Here we have used (iii) in (2.2) for the last equality. and . The asymptotic scheme is the same as that of the original -DG method. If and are both piecewise constant, then the asymptotic scheme has the primal form: find , such that

 ∑K∈Th∫K(13σs∇ψ(0)h⋅∇φ+σaψ(0)hφ)dx=∫Dqφdx, (3.11)

. This is the classical continuous Galerkin approximation, which is stable and second-order accurate.

3. When , then , and . With elements and triangular meshes, the asymptotic scheme is essentially the scheme suggested by Egger and Schlottbom in  with . If elements and Cartesian meshes are used, the scheme yields the same variational form as that in , while the space pair no longer satisfies the condition .

From another point of view, suppose and are piecewise constant, the primal form is: find , such that

 ∑K∈Th∫K(13σsΠ0(∇ψ(0)h)⋅Π0(∇φ)+σaψ(0)hφ)dx=∑K∈Th∫Kqφdx, (3.12)

. For elements on triangular meshes, (3.12) is identical to (3.11). For elements on Cartesian meshes, one can show that (3.12) is unisolvent. Furthermore, , if . While the accuracy is hard to analyze under the finite element framework. Assume a uniform square mesh with cell length . Let and be globally constant. Then (3.12) can be rewritten as a finite difference scheme under the Lagrange basis functions.

 −ψi−1,j−1+ψi−1,j+1−4ψi,j+ψi+1,j−1+ψi+1,j+13σs⋅2h2+σaA[ψi,j]= A[qi,j], (3.13) A[ψj]:=136(ψi−1,j−1+ψi−1,j+1+ψi+1,j−1+ψi+1,j+1) (3.14) +19(ψi−1,j+ψi,j−1+