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
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 pointwhere 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
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.
|Unknowns per cell||Triangles ()||Rectangles ()|
|Memory cost ratio as|
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 -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.
Consider a quadrature rule with points and positive weights such that
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,
The method approximates the angular flux at the quadrature points by a vector-valued function whose components satisfy a coupled system with equations
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
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
The semi-norm induced by jumps at the cell interfaces is given by
To construct the -DG method, define the local operators equationparentequation
where is the upwind trace at , and is defined as zero when the limit is taken outside of . Then set
The -DG method is then: find such that
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
can be assembled into a linear system (detailed in Appendix A)
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
for the vector . This equation is derived by applying and then to (2.13). In a second step is recovered from the relation
The following theorem is proven in Appendix B.
The matrix is invertible.
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,
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
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
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
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
based on a proper subspace that maintains the diffusion limit, but requires fewer unknowns for a given mesh .
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
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),
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:
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),
where is the tensor product norm of in , the bound (ii) implies further that is uniformly bounded and hence converges subsequentially to a limit .
We combine and , using the fact that and invoking (2.2). This gives
Finally, the right-hand side of (3.1) is (for )
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
. 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
. 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.