Superconvergence Extraction of Upwind Discontinuous Galerkin Method Solving the Radiative Transfer Equation

TL;DR

This paper analyzes the superconvergence properties of the upwind discontinuous Galerkin method for radiative transfer equations and demonstrates how SIAC filters significantly improve accuracy and computational efficiency.

Contribution

It provides a theoretical proof of superconvergence orders and shows how SIAC filters enhance accuracy and efficiency in solving radiative transfer equations.

Findings

Superconvergence orders of (2k+2) and (2k+1/2) proven for steady-state and time-dependent problems.

SIAC filters improve accuracy by approximately 2.22 times for time-dependent and 4-9 times for steady-state problems.

Filtering reduces computational time by nearly 20 times for time-dependent problems without accuracy loss.

Abstract

We theoretically analyze the superconvergence of the upwind discontinuous Galerkin (DG) method for both the steady-state and time-dependent radiative transfer equation (RTE), and apply the Smooth-Increasing Accuracy-Conserving (SIAC) filters to enhance the accuracy order. Direct application of SIAC filters on low-dimensional macroscopic moments, often the quantities of practical interest, can effectively improve the approximation accuracy with marginal computational overhead. Using piecewise $k$-th order polynomials for the approximation and assuming constant cross sections, we prove $(2k+2)$-th order superconvergence for the steady-state problem at Radau points on each element and $(2k+1/2)$-th order superconvergence for the global $L^2$ and negative-order Sobolev norms for the time-dependent problem. Numerical experiments confirm the efficacy of the filtering, demonstrating post-filter convergence orders of $2k+2$ for steady-state and $2k+1$ for time-dependent problems. More significantly, the SIAC filter delivers substantial gains in computational efficiency. For a time-dependent problem, we observed an approximately $2.22 \times$ accuracy improvement and a $19.94 \times$ reduction in computational time. For the steady-state problems, the filter achieved a $4$--$9 \times$ acceleration without any loss of accuracy.

Equations

$${\partial }_t\psi (\mathbf{x},\mathbf{\Omega },t)+\mathbf{\Omega }\cdot {\nabla }_{\mathbf{x}}\psi (\mathbf{x},\mathbf{\Omega },t)={\sigma }_s(\overline{\psi }(\mathbf{x},t)-\psi (\mathbf{x},\mathbf{\Omega },t))-{\sigma }_a\psi (\mathbf{x},\mathbf{\Omega },t),$$

$$\overline{\psi }(\mathbf{x},t)=⟨\psi (\mathbf{x},\mathbf{\Omega },t)⟩=\frac{1}{4\pi }{\int }_{{\mathbb{S}}^2}\psi (\mathbf{x},\mathbf{\Omega },t)d\mathbf{\Omega },$$

$$v{\partial }_x\psi (x,v,t)={\sigma }_s(\overline{\psi }(\mathbf{x},t)-\psi (x,v,t))-{\sigma }_a\psi (x,v,t),$$

$$\overline{\psi }(x,t)=\frac{1}{2}{\int }_{-1}^1\psi (x,v)dv.$$

$$\text{“}\psi :=e^{{\sigma }_{a}t}\psi \text{”}.$$

$$\begin{array}{rl}& {\partial }_t\psi (\mathbf{x},\mathbf{\Omega },t)+\mathbf{\Omega }\cdot {\nabla }_{\mathbf{x}}\psi (\mathbf{x},\mathbf{\Omega },t)=\sigma (\overline{\psi }(\mathbf{x},t)-\psi (\mathbf{x},\mathbf{\Omega },t)),\\ & \overline{\psi }(\mathbf{x},t)=\frac{1}{4\pi }{\int }_{\mathcal{S}}\psi (\mathbf{x},\mathbf{\Omega },t)d\mathbf{\Omega }\end{array}$$

$$D^{|\alpha |}f(\mathbf{x})={\partial }_{x_1}^{{\alpha }_1}{\partial }_{x_2}^{{\alpha }_2}\cdots {\partial }_{x_d}^{{\alpha }_d}f(\mathbf{x}),$$

$$\Vert f{\Vert }_{L^p(\mathbf{X})}={\left({\int }_{\mathbf{X}}|f{|}^{p}d\mathbf{x}\right)}^{1/p},$$

$$W^{k,p}(\mathbf{X}):=\{v\in L^p(\mathbf{X}):D^{|\alpha |}f\in L^p(\mathbf{X}),|\alpha |\le k\},$$

$$\Vert f{\Vert }_{W^{k,p}(\mathcal{D})}=\sum _{|\alpha |\le k}\Vert f{\Vert }_{L^n(\mathcal{D})}.$$

$$\Vert f{\Vert }_{k,\mathbf{X}}=\Vert f{\Vert }_{H^k(\mathbf{X})}$$

$$(f,g{)}_{\mathbf{X}}={\int }_{\mathbf{X}}fgd\mathbf{x}.$$

$$\Vert F{\Vert }_{-\ell ,\mathbf{X}}=\underset{\Phi \in (C_0^{\infty }(\mathbf{X}){)}^N}{\sup }\left[\sum _{j=1}^N{\omega }_j(f^j,{\varphi }^j{)}_{\mathbf{X}}\right]/{\left(\sum _{j=1}^N{\omega }_j\Vert {\varphi }^j{\Vert }_{\ell ,\mathbf{X}}^2\right)}^{1/2}.$$

$$\begin{array}{rl}& {\partial }_t{\phi }^j+{\mathbf{\Omega }}_j\cdot {\nabla }_x{\phi }^j+\sigma (\overline{\phi }-{\phi }^j)=0,\\ & \overline{\phi }=\frac{1}{m(\mathcal{S})}\sum _{j=1}^{N_{N_{\mathbf{\Omega }}}}{\omega }_j{\phi }^j.\end{array}$$

$$\frac{d}{dt}\sum _{j=1}^{N_{\mathbf{\Omega }}}{\omega }_j({\psi }^j,{\phi }^j{)}_{\mathbf{X}}=0.$$

$$L_0(\xi )=1,L_1(\xi )=\xi ,L_{k+1}(\xi )=\frac{2k+1}{k+1}{L}_k(\xi )-\frac{k}{k+1}{L}_{k-1}(\xi ),k\ge 1.$$

$${\delta }_{jk}=\{\begin{array}{l}1,j=k\\ 0,j\ne k\end{array}.$$

$$R_k^{\pm }(\xi )=\{\begin{array}{ll}{L}_k(\xi )\pm L_{k-1}(\xi ),& \text{if}k\ge 1,\\ 0,& \text{if}k=0.\end{array}$$

$${\overline{f}}_h=\sum _{j=1}^N{\omega }_j{f}_j,$$

$$\left|\frac{1}{m(\mathcal{S})}{\int }_{\mathcal{S}}f(\mathbf{\Omega })d\mathbf{\Omega }-\sum _{j=1}^N{\omega }_{j}f({\mathbf{\Omega }}_j)\right|\le C_s{N}^{-s}\Vert f{\Vert }_{s,\mathcal{S}}.$$

$${\int }_{\mathbf{X}}d\mathbf{x}=\sum _{K\in {\mathcal{T}}_h}{\int }_{K}d\mathbf{x},{\int }_{\mathbf{s}}d{\mathcal{E}}_h=\sum _{e\in {\mathcal{E}}_h}{\int }_{\mathbf{s}}d\mathbf{s}.$$

$$V_{h}^{k}=\{v\in L^{2}(\boldsymbol{X}):v\bigg{|}_{K}\in\mathbb{Q}^{k}(K),\,K\in\mathcal{T}_{h}\}.$$

$$[[g{]]}_{\mathbf{s}}=g^{+}{\mathbf{n}}_{\mathbf{s}}^{+}-g^{-}{\mathbf{n}}_{\mathbf{s}}^{-},\{\{g\}{\}}_{\mathbf{s}}=\frac{1}{2}\left(g^{+}+g^{-}\right),$$

$${\partial }_t{\phi }^j+{\mathbf{A}}_j(\mathbf{x})\cdot {\nabla }_{\mathbf{x}}{\phi }^j+\sigma \left[\overline{\phi }-{\phi }^j\right]=0,1\le j\le N_{\mathbf{\Omega }}.$$

$$\sum _{j=1}^{N_v}{\omega }_j\Vert {\phi }^j(\cdot ,t){\Vert }_{\ell ,\mathbf{X}}^2\le C\sum _{j=1}^{N_v}{\omega }_j\Vert {\phi }^j(\cdot ,T){\Vert }_{\ell ,\mathbf{X}}^2.$$

$$\begin{array}{rl}\frac{1}{2}\frac{d}{dt}\sum _{j=1}^{N_v}{\omega }_j\Vert {\phi }^j(\cdot ,t){\Vert }_{\ell ,\mathbf{X}}^2=& -\sigma \sum _{j=1}^{N_v}{\omega }_j{\int }_{\mathbf{X}}(\overline{\phi }-{\phi }_j){\phi }_{j}d\mathbf{x}\\ & =-\sigma \sum _{j=1}^{N_v}{\omega }_j{\int }_{\mathbf{X}}(\overline{\phi }-{\phi }_j{)}^{2}d\mathbf{x}\ge 0.\end{array}$$

$$\Vert {\eta }_h{\Vert }_{L^2(K)}+h_K\Vert {\nabla }_x{\eta }_h{\Vert }_{L^2(K)}+h_K^{1/2}\Vert {\eta }_h{\Vert }_{L^2(\partial K)}\le C{h}_K^{k+1}\Vert \psi {\Vert }_{k+1,K},\forall K\in {\mathcal{T}}_h,$$

$$\Vert {\nabla }_{x}g{\Vert }_{L^2(K)}\le C{h}_K^{-1}\Vert g{\Vert }_{L^2(K)},\Vert g{\Vert }_{L^2(\partial K)}\le C{h}_K^{-1/2}\Vert g{\Vert }_{L^2(K)},$$

$$v{\partial }_x{\psi }^j(x,v_j,t)+{\sigma }_t{\psi }^j(x,v_j,t)={\sigma }_s\overline{\psi }(x,t)$$

$$\overline{\psi }(x,t)=\frac{1}{2}\sum _{j=1}^{N_v}{\omega }_j{\psi }^j(x).$$

$${\partial }_t{\psi }^j(\mathbf{x},t)+{\mathbf{\Omega }}_j\cdot {\nabla }_{\mathbf{x}}{\psi }^j(\mathbf{x},t)=\sigma (\overline{\psi }(\mathbf{x},t)-{\psi }^j(\mathbf{x},t))$$

$$\overline{\psi }(x)=\frac{1}{m(\mathcal{S})}\sum _{j=1}^{N_{\mathbf{\Omega }}}{\omega }_j{\psi }^j.$$

Full text

Using piecewise $k$-th order polynomials for the approximation and assuming constant cross sections, we prove $(2k+2)$-th order superconvergence for the steady-state problem at Radau points on each element and $(2k+\frac{1}{2})$-th order superconvergence for the global $L^{2}-$ and negative-order Sobolev norms for the time-dependent problem.

Numerical experiments confirm the efficacy of the filtering, demonstrating post-filter convergence orders of $2k+2$ for steady-state and $2k+1$ for time-dependent problems. More significantly, the SIAC filter delivers substantial gains in computational efficiency. For a time-dependent problem, we observed an approximately $2.22\times$ accuracy improvement and a $19.94\times$ reduction in computational time. For the steady-state problems, the filter achieved a $4–9\times$ acceleration without any loss of accuracy.

1 Introduction

The radiative transfer equation (RTE) ([40, 30, 8]) describes particles, e.g. neutrons or photons, propagating through and interacting with a background medium by means of a kinetic distribution that gives the density of particles with respect to the local phase space measure . It has a wide range of applications in nuclear engineering, medical imaging, astrophysics and remote sensing. Efficient numerical simulation tools for RTE are highly desired. In this article, we take advantage of the underlying superconvergent properties of the approximation to improve accuracy as well as computational performance.

In order to approximate solutions to RTE, the upwind discontinuous Galerkin (DG) method is used. DG has been a popular and powerful deterministic solver for RTE since its first introduced by Reed and Hill in simulating neutron transport [41]. A key strength of the high-order upwind DG method is its asymptotic preserving (AP) property, which ensures it correctly captures the diffusion limit of the RTE without resolving the small particle mean free path. This property was formally established through asymptotic analysis [29, 2], and later rigorously proved [22]. Uniform convergence with respect to the mean free path has been analyzed in [43, 44]. Besides upwind DG, DG methods have been actively developed for solving RTE. These developments include hp-adaptive hybridized DG approach [15], approaches leveraging micro-macro decomposition and implicit-explicit (IMEX) time integrator [26, 37, 39, 38, 1], semi-Lagrangian approach [7], and high-order low-order method (HOLO) [16].

Despite its success, the discontinuous Galerkin (DG) method for the radiative transfer equation (RTE) suffers from the curse of dimensionality, as the unknown distribution function is posed in a high-dimensional phase space. An attractive mitigation is post-processing techniques that can significantly enhance resolution of numerical solutions. Exploiting pointwise or Fourier space superconvergence, post-processing techniques typically enhance accuracy of a finite element solution through a convolution that performs effective local averaging to damp dominating high frequency errors. The development of post-processing dates back to Bramble and Schatz [6] in the context of continuous finite elements methods for elliptic problems, which exploits superconvergence in negative-order norms. Formally defined as weighted dual norms, negative-order norms can be seen as a bridge connecting information in the Fourier/signal space with the one in physical space. For time-dependent problems, the DG method is proven to exhibit higher-order accuracy in negative-order norms— $2k+1$ for dispersion errors and $2k+2$ for dissipative errors—compared to only $k+1/2$ in the $L^{2}$-norm. The pioneering work of Cockburn, Luskin, Shu, and Süli [12] established the framework for superconvergence analysis in negative norms and post-processing for time-dependent linear hyperbolic equations. This relies on a properly designed kernel that respects moment [35]. This technique was later extended to nonlinear scalar conservation laws [28, 27, 32], nonlinear symmetric systems [33], and, more recently, to the collisionless Vlasov–Maxwell system [21]. Besides exploiting negative-order norm, super-convergence for time-dependent problems can also be proved through a direct Fourier analysis [23]. Beyond signal space, superconvergence of DG methods for steady-state linear and nonlinear hyperbolic problems is analyzed through polynomial expansions of error equations [4, 5].

In this paper, we exploit post-processing based on Smooth-Increasing Accuracy-Conserving (SIAC) filters [46, 34, 14] to enhance the accuracy of both steady state and time-dependent RTE, which is a collisional kinetic equation. In practice, the desired quantity of interest is often low-dimensional macroscopic moments (e.g., the density) rather than the high-dimensional particle distributions. We theoretically establish the superconvergence results for both steady state and time-dependent RTE and numerically demonstrate that directly applying SIAC filters to these macroscopic moments significantly enhances their approximation accuracy. Since the effectiveness of SIAC filter relies on the superconvergence of the DG scheme, we also establish superconvergence results for both steady state and time-dependent RTE. The main contributions of this paper are as follows.

For the steady-state problem, we prove $k+2$ convergence order at Radau points within each element and $2k+2$-th order at the outflow edge. This result is established by extend the framework based on coefficient matching for a polynomial expansion of the error equation [4, 5] which proves $2k+2$-th order convergence at outflow flux and $k+2$-th order at interior Radau points for linear advection and nonlinear scalar hyperbolic problems. Instead of using a key argument inapplicable to RTE based on the theory for the Runge-Kutta method, we have a novel presentation of the analysis and show that the local polynomial is the same on each element and only weighted by the error coefficients. 2. 2.

For the time-dependent problem with constant cross-section, we establish, for the first time, an $L^{2}-$ error estimate for the DG-DO method. By carefully designing a dual problem, we prove a $2k+\tfrac{1}{2}$ superconvergence rate in a negative Sobolev norm and observe a $2k+1$ convergence rate numerically. This analysis, which addresses discretization in both space and angle for time-dependent problems, is new. Previous works focused either on steady-state problems [24, 45] or assumed continuity in physical space [17, 20, 9]. 3. 3.

Through a series of benchmark tests, we validate our theoretical results and showcase the SIAC filter’s capacity to drastically improve computational efficiency. By enhancing the accuracy achievable on coarse meshes, the filter reduces the time required to reach a target error tolerance in comparison to mesh refinement. Specifically, we observe an approximately $3.26\times$ accuracy boost along with a $19.65\times$ computational time reduction for a time-dependent problem, and $4\times$ to $9\times$ speedup without sacrificing accuracy for a steady state problem.

This paper is organized as follows. We introduce the model problem, notation, preliminaries and the discretization in Sec. 2, our superconvergence results are presented in Sec. 3, sketch the key steps for our proofs in Sec. 4. The SIAC filter is introduced in Sec. 5 and the numerical tests of our methods are presented in Sec. 6. Our conclusions are discussed in Sec. 7.

2 Background

We consider one-group, linear radiative transfer equation with isotropic scattering:

[TABLE]

equipped with proper boundary conditions. Here, $\psi(\boldsymbol{x},\boldsymbol{\Omega},t)$ is the particle distribution function (also known as intensity or angular flux) for spatial location $\boldsymbol{x}\in\boldsymbol{X}\subset\mathbb{R}^{3}$, angular direction on the unit sphere $\boldsymbol{\Omega}\in\mathbb{S}^{2}$ and time $t$. The macroscopic density (also known as scalar flux), as $\overline{\psi}(\boldsymbol{x},t)=\langle{\psi(\boldsymbol{x},\boldsymbol{\Omega},t)}\rangle=\frac{1}{4\pi}\int_{\mathbb{S}^{2}}\,\psi(\boldsymbol{x},\boldsymbol{\Omega},t)\,d\boldsymbol{\Omega}$ is the average of the distribution function with respect to angular direction, and $\langle{\cdot}\rangle$ represents the averaging operator over angular space. Here, $\sigma_{s}(\boldsymbol{x})\geq 0$ is an isotropic scattering cross section, while $\sigma_{a}(\boldsymbol{x})\geq 0$ is the absorption cross section. The total cross section is defined as $\sigma_{t}(\boldsymbol{x})=\sigma_{s}(\boldsymbol{x})+\sigma_{a}(\boldsymbol{x})$. In this work we will assume the cross-sections $\sigma_{t}$ and $\sigma_{a}$ are constant in physical space $\boldsymbol{x}$.

In our analysis for the steady-state problem, we focus on the one-dimensional slab geometry. For a one-dimensional slab geometry, the particle distribution $\psi(\boldsymbol{x},\boldsymbol{\Omega},t)$ only depends on a location $x\in\mathbb{R}$ and the cosine of angle between the angular direction and the $x$-axis $\boldsymbol{\Omega}\in[-1,1]$. With this assumption, the equation can be simplified as $v=\cos(\Omega)$:

[TABLE]

On the other hand, following [20], to help the error analysis in the time, since the cross-sections $\sigma_{a}$ and $\sigma_{t}$ are constant in space, we can reduce the analysis to Equation (1a) to a purely scattering system for the function (abusing notation)

[TABLE]

$\psi$ satisfies

[TABLE]

where $\sigma:=\sigma_{s}$. Henceforth, we restrict our attention to (3). The corresponding results for the case of nonzero absorption can then be recovered by reversing the transformation, which introduces exponential time decay when $\sigma_{a}>0$.

In order to define the numerical discretization and prove superconvergence, we first note a few preliminaries.

2.1 Notation and preliminaries

In this subsection, we discuss notation and estimates that will be used throughout this article.

2.1.1 Notation

Spaces and Norms. In this article, we denote the domain for angular space as $\mathcal{S}=\mathbb{S}^{d-1}$ and for Euclidean space as $\boldsymbol{X}\subset\mathbb{R}^{d}.$ $\boldsymbol{X}$ refers to a general domain – either in Euclidean or angular space. Hence, when $\boldsymbol{X}\subset\mathbb{R}^{3},\,\mathcal{S}=\mathcal{S}.$

We start by defining the multi-index $\alpha$, where $\alpha=(\alpha_{1},\dots,\alpha_{d}),\quad|\alpha|=\sum_{j=1}^{d}\,\alpha_{j}$ and

[TABLE]

where $f$ is some given function.

The usual $L^{p}-$norm is defined as

[TABLE]

where $\|{f}\|_{L^{\infty}({\mathcal{D}})}=\text{ess sup}_{\mathcal{D}}\,|f|$.

The general Sobolov space for $k\in\mathbb{N}_{+},\,1\leq p\leq\infty$ is defined as

[TABLE]

with norm

[TABLE]

We note that when $p=2,$ we denote $H^{k}=W^{k,2}$ we simplify notation as

[TABLE]

and define the $L^{2}-$inner product as

[TABLE]

Further, for purposes of error analysis contained in this article, we define a weighted negative-order Sobolev norm for $F=(f^{1},\cdots,f^{N})$ as

[TABLE]

The weights $\omega_{j}$ correspond to the quadrature weight ones in the $S_{N}$ discretization for the phase-angle space; a more detailed discussion will be provided below.

Dual problem. To obtain superconvergence for the time-dependent equation, we will need to bound the negative-order Sobolev norm of the divided differences of the error. This requires defining the dual problem. The dual problem for $j=1,\dots,N_{\boldsymbol{\Omega}}$ is defined as

[TABLE]

together with the final time condition $\varphi^{j}(x,T)=\varphi^{j}(x)\in C^{\infty}(X)$.

Note that for this dual problem satisfies

[TABLE]

This is an important property that allows us to phrase our error estimates in terms of the initial projection as well as the temporal growth in jumps across element interfaces.

Approximation polynomials. Throughout this article, we will utilize the Legendre polynomials defined on a reference interval $[-1,1]$, which are defined by the recurrence relation

[TABLE]

These satisfy $\int_{-1}^{1}\,L_{k}(\xi)L_{j}(\xi)d\xi=\frac{2}{2k+1}\delta_{jk}$ with $\delta_{jk}$ being the Kronecker delta

[TABLE]

The Radau polynomials will also be useful. The $k$-th order left Radau polynomial, $R_{k}^{+},$ and the right Radau polynomial, $R_{k}^{-},$ on the reference element $[-1,1]$ are defined as

[TABLE]

Note that $R_{k}^{+}(-1)=0$ and $R_{k}^{-}(1)=0$. In other words, the left (right) end point of the reference interval is a root of the left (right) Radau polynomial.

Quadrature. For the discretization in phase space-angle, we will need to define a quadrature that can integrate exactly polynomials of total degree less than or equal to $N$ ($|\alpha|=N$) [18] and denote it by:

[TABLE]

with $\omega_{j}$ denoting the quadrature weights and $f_{j}$ denoting the function $f$ evaluated at the quadrature points. The weights are defined such that $\sum_{j=1}^{N}\,\omega_{j}=1$.

For the quadrature, we have the following approximation result:

Lemma 1** (Quadrature accuracy bound ([25], Corollary)).**

Let $f\in H^{s}(\mathcal{S}),\,s>1$ (cf. [18])). Then

[TABLE]

$C_{s}$* is an universal constant depending only on $s$ and $m(\mathcal{S})$ represents the Lebesgue measure $\mathcal{S}$.*

Discretization. Let $\mathcal{T}_{h}$ be a spatial tessellation of the domain $\boldsymbol{X}$ and $K$ an element in the tessellation. For simplicity of the analysis, we assume that $K$ is a rectangular element, that is $K$ is defined as $K=I_{K,1}\otimes\dots\otimes I_{K,d}$, with $I_{K,j}$ being a one-dimensional interval. For simplicity, we denote $h_{K}$ to be the size of $K$. We also introduce the edges, $\mathbf{s}$, and $\mathcal{E}_{h}=\cup\,\mathbf{s}.$ For convenience, we introduce the shorthand notation

[TABLE]

Given an integer $k\geq 1$, a polynomial approximation space for the tensor product polynomials, $\mathbb{Q}^{k}=\mathbb{P}^{k}(I_{K,1})\otimes\dots\otimes\mathbb{P}^{k}(I_{K,d})$, is defined as

[TABLE]

On the edge, $\mathbf{s}$, of element $K$, $\mathbf{n}_{x}^{\pm}$ denotes the unit outward normal across $K$. This allows for defining the jumps and averages on edges. In order to do so, denote $g^{\pm}=\lim_{x\rightarrow\mathbf{s}^{\pm}}g$. The jumps and averages across $\mathbf{s}$ are then defined as

[TABLE]

respectively. Note that we will often drop the superscript for the normal vector and utilize $\mathbf{n}_{x}^{+}=-\mathbf{n}_{x}^{-}=-\mathbf{n}_{x}.$

2.1.2 Estimates

In this subsection we discuss some of the basic estimates that will be useful in the proof of superconvergence. We note that we do not distinguish between constants and denote them by $C$.

We begin by discussing the following regularity results for a linear PDE system, which will be useful for the proof in the time-dependent case.

Lemma 2** (Regularity of the dual problem).**

Consider the following system of equations with periodic boundary conditions in $\boldsymbol{x}$ for all $t\in[0,T]$:

[TABLE]

where the given functions $\mathbf{A}_{j}\in W^{l+1,\infty}$ are periodic, and satisfy the divergence free constraint $\nabla_{x}\cdot\mathbf{A}_{j}(x)=0$. Then for all $t\in[0,T]$ there exists a constant independent of time $t$ satisfying:

[TABLE]

Proof.

For the time dependent case we will assume $\sigma$ is constant. By testing the equation above by $\varphi^{j}$ for each $j$, and using the periodic and divergence free properties of the $\mathbf{A_{j}}$’s, we have the following,

[TABLE]

The conclusion follows. Now since we are considering the full Sobolev norm, we still need to estimate the $L^{2}$- weighted norms of the higher order spatial derivatives $\partial_{x_{j}}^{\ell}$, to do so we apply $\partial_{x_{j}}$ to the equation above and then we repeat the steps above. ∎

For the discretization, the following approximation properties for $V_{h}^{k}$, as well as some inverse inequalities [10] will be necessary. Here, we denote by $\Pi^{k}$ the usual $L^{2}-$projection onto $V_{h}^{k}$ and define $\eta_{h}^{j}=\Pi^{k}\psi^{j}-\psi^{j}$ as the projection error as $\xi^{j}_{h}=\Pi^{k}\psi^{j}-\psi^{j}_{h}$ as the error in space. Then, $e^{j}_{h}=\psi^{j}-\psi^{j}_{h}=\xi^{j}_{h}-\eta^{j}_{h}$ and we have the following estimate for the projection error:

Lemma 3** (Approximation properties[10]).**

There exists a constant $C>0$, such that for any $g\in H^{k+1}(\mathcal{D})$,

[TABLE]

where the constant $C$ is independent of the mesh size $h_{{K}}$, but depends on the polynomial degree, ${K}$, and the regularity of $\psi$.

We will also make use of the following inverse inequality:

Lemma 4** (Inverse inequality [10]).**

There exists a constant $C>0$, such that for any $\psi\in\mathbb{P}^{k}({K})$ the following holds:

[TABLE]

where the constant $C$ is independent of the mesh size $h_{K}$, but depends on polynomial degree $k$ and the shape regularity of the mesh.

2.2 Discretization in angular space

In angular space, we apply the discrete ordinates method $S_{N_{v}}$ for the 1D-steady-state case and $S_{N_{\boldsymbol{\Omega}}}$ for the time-dependent model. For the steady state case this is the Gauss-Legendre quadrature points $\{v_{j}\}_{j=1}^{N_{v}}$ and weights $\{\omega_{j}\}_{j=1}^{N_{v}}$ in $[-1,1]$. By approximating the macroscopic density with the numerical quadrature, we obtain the $S_{N_{v}}$ system which seeks $\psi(x,v_{j},t)\approx\psi^{j}(x,t)$ satisfying

[TABLE]

That is, equation (2) is solved at a set of quadrature points. Analogously for the time-dependent problem (3) using the $S_{N_{\boldsymbol{\Omega}}}$ nodes $\{\boldsymbol{\Omega}_{j}\}_{j=1}^{N_{v}}$ and weights $\{\omega_{j}\}_{j=1}^{N_{v}}$, we obtain the $S_{N_{\boldsymbol{\Omega}}}$ system which seeks $\psi(\boldsymbol{x},\boldsymbol{\Omega}_{j},t)\approx\psi^{j}(\boldsymbol{x},t)$ satisfying

[TABLE]

where $\boldsymbol{x}\in\boldsymbol{X},\;1\leq j\leq N_{\boldsymbol{\Omega}}.$

We note that the accuracy of the $S_{N}$ approximation is given by lemma 1 applied to angular space.

2.3 Spatial discretization via the DG method

We apply high order upwind discontinuous Galerkin (DG) method to the $S_{N}$ system (21). That is, we seek $\psi^{j}_{h}\in V_{h}^{k}$ such that for all $\tau_{h}\in V_{h}^{k}$ and $1\leq j\leq N_{\boldsymbol{\Omega}}$:

[TABLE]

Here, $\widehat{\left({\boldsymbol{\Omega}_{j}\psi^{j}_{h}}\right)}_{\partial K}$ denotes the upwind numerical flux defined by :

[TABLE]

where $\mathbf{n}$ is the outward normal direction with respect to the element $K$.

For the one-dimensional case, we consider the case of steady-state and denote the domain as $X=[x_{L},x_{R}]$. For simplicity, we assume a uniform mesh given by $x_{L}=x_{0}<x_{1}<\dots<x_{N_{x}}=x_{R}$ with mesh size $h=\frac{x_{R}-x_{L}}{N_{x}}$. Additionally, we set $v_{j}=\cos(\boldsymbol{\Omega}_{j}).$ Then, the upwind DG scheme is: Find $\psi^{j}_{h}\in V_{h}^{k},$ and $j=1,\dots,N_{v}$ such that:

[TABLE]

where $K_{n}=(x_{n-1},x_{n})$ represents the $n-$th element. Here, $x_{n}^{\pm}$ stands for the right and left limit of the function value. For the one-dimensional case, the upwind flux is simply

[TABLE]

Particularly, given inflow boundary conditions, $\psi^{j}_{h}(x^{-}_{0})=\psi^{j}(x_{0},v_{j})$ for $v_{j}>0$ and $\psi^{j}_{h}(x^{+}_{N})=\psi^{j}(x_{N},v_{j})$ for $v_{j}<0$.

We emphasize that the upwind DG discretization with polynomial order $k\geq 1$ is proven to be asymptotically preserving for the radiative transport equation (RTE) [3, 22, 43]. Hence, it can capture the correct diffusion limit of RTE without resolving the particle mean free path.

2.4 Error estimates for the approximation

For the error estimates, recall that we define $\eta^{j}_{h}=\Pi^{k}\psi^{j}-\psi^{j}$ to be the error between the projection of the exact solution and the exact solution and $\xi^{j}_{h}=\Pi^{k}\psi^{j}-\psi^{j}_{h}$ to be the error between the projection and the approximation. It is easy to see that $e^{j}_{h}=\psi^{j}-\psi^{j}_{h}=\xi^{j}_{h}-\eta^{j}_{h}$. From the approximation, Lemma 3, the following inequality holds:

[TABLE]

This will be used frequently.

We also have the following estimate:

Theorem 5**.**

Let $k\geq 0$ and let $\psi^{j},\,1\leq j\leq N_{\boldsymbol{\Omega}}$ be the exact solutions to (1) and assume that $\psi^{j}\in C([0,T];H^{k+1}(\boldsymbol{X})),\,1\leq j\leq N_{\boldsymbol{\Omega}}$. If $\psi^{j}_{h}$ is an approximation obtained via Equation (22) with the numerical initial conditions $\psi^{j}_{h}(\cdot,0)=\Pi^{k}\psi^{j}_{0}$, then

[TABLE]

Here, the constant $C$ depends on the upper bound of $\sum_{j=1}^{N_{\boldsymbol{\Omega}}}\,\omega_{j}{\|{\psi^{j}}\|_{k+1,\boldsymbol{X}}}$ .

Now we analyze the error coming from the angular part.

Theorem 6**.**

Let $N$ be the accuracy of the $S_{N}$ quadrature, if $\psi\in C([0,T];L^{2}(\boldsymbol{X};H^{s}(\mathcal{S})))$, then we have the following error estimate

[TABLE]

That depends on $c_{s}$ is an universal constant depending only on $s$ and the upper bounds of $\|{\psi}\|_{L^{2}(\boldsymbol{X};H^{s}(\mathcal{S}))}$.

For the proofs see the Appendix.

By combining Theorems 5 and theorem 6, as well as the triangle inequality, we can conclude the following estimate:

Theorem 7**.**

For polynomial degree $k\geq 1$, if $\psi\in C([0,T];L^{2}(\boldsymbol{X};H^{s}(\mathcal{S})))$ and each $\psi^{j}\in C([0,T];H^{k+1}(\boldsymbol{X})),\,j=1,\dots,N_{\boldsymbol{\Omega}}$, the semi-discrete discrete-ordinates approximation via the Discontinuous Galerkin method for the time-dependent radiation transport equation, eq. 1, has the following error estimate

[TABLE]

where $C$ and $D$ depend on upper bounds of $\sum_{j=1}^{N_{\boldsymbol{\Omega}}}\,\omega_{j}{\|{\psi^{j}}\|_{k+1,\boldsymbol{X}}},\,\|{\psi}\|_{L^{2}(\boldsymbol{X};H^{s}(\mathcal{S}))}$ respectively.

3 Superconvergence and extraction

In this section we discuss the superconvergence of the discontinuous Galerkin approximation to Equation eq. 1. We begin by introducint the superconvergent points in the one-dimensional steady-state case. We then proceed to discuss the multi-dimensional time-dependent case and the underlying superconvergence in the negative-order norm. In Section section 5 we show how to extract this information via the Smoothness-Increasing Accuracy-Conserving (SIAC) filter.

In our superconvergence analysis, we assume periodic or zero inflow boundary conditions as well as a constant scattering cross section, $\sigma_{s}$. Note that though our analysis is restricted to constant scattering cross sections, we numerically observe improved accuracy after post-processing for general cross sections.

3.1 Superconvergence of the steady state problem for a 1D slab geometry

In this section, we discuss the superconvergent points of the one-dimensional steady-state model (19).

The discontinuous Galerkin approximation given in section 2.3 in terms of the local coordinate mapping becomes: Find $\psi^{j}_{h}\in V_{h}^{k}$ such that

[TABLE]

The mapping to the reference element, $\mathcal{I}=[-1,1]$, is defined by $\zeta=\frac{2}{h}\left(x-\overline{x_{n}}\right)-1,$ where $\overline{x_{n}}=\frac{1}{2}(x_{n}+x_{n+1}).$ We then have the following theorem:

Theorem 8**.**

(Superconvergence at Radau points for steady state problem) Suppose $\psi^{j}$ is a Lipschitz continuous function belonging to $C^{\infty}([x_{L},x_{R}])$, with $\sigma_{s}$ and $\sigma_{t}$ being constant, and $\psi^{j}_{h}\in V_{h}^{k}$ being the approximation to the $S_{N}$ system (19) and its DG approximation (2.3). Define the error function as $e_{h}^{j}(x,t)=\psi^{j}(x,t)-\psi^{j}_{h}(x,t)$. When the mesh size is sufficiently small, the following statements hold.

When $v_{j}>0$, $e_{h}^{j}(x_{p}^{-},t)=\mathcal{O}(h^{2k+2})$ for $1\leq p\leq N_{v}$. When $v_{j}<0$, $e_{h}^{j}(x_{n-1}^{+},t)=\mathcal{O}(h^{2k+2})$ for $1\leq p\leq N_{v}$. 2. 2.

Let the roots of the $(k+1)$-th order right and left Radau polynomial $R^{\pm}_{k+1}$ be $\{\nu^{+}_{\ell}\}_{\ell=1}^{k}\bigcup\{1\}$ and $\{-1\}\bigcup\{\nu^{-}_{\ell}\}_{\ell=1}^{k}$. For the interior roots, when $v_{j}>0$,

[TABLE]

for $1\leq\ell\leq k$. And when $v_{j}<0$,

[TABLE]

for $1\leq\ell\leq k$.

Notice that when $\sigma_{s}=\sigma_{t}=0$, the $S_{N}$ system (2) degenerates to a $N_{v}$ decoupled linear advection equation. As a result, Theorem 8 implies that the upwind DG method has $2k+2$ order of accuracy at downwind edges. In other words, for linear advection, the global superconvergence order for the downwind edge of each element is $2k+2$ instead of $2k+1$ proved in [4]. Our proof, which will be outlined in Sec. 4, can be seen as an extension of [4].

3.2 Superconvergence for time dependent case

The time-dependent case considers multi-dimensional approximations to eq. 21. For this analysis, we will need to establish a bound on a weighted negative-order norm of the error as we utilized the $S_{N}$ system.

Theorem 9**.**

Let $k\geq 0$, Let $\psi^{j}$, $1\leq j\leq N_{\boldsymbol{\Omega}}$ be the exact solutions to (1) and assume that $\psi^{j}\in C([0,T];H^{k+2}(X))$, $1\leq j\leq N_{\boldsymbol{\Omega}}$. If $\psi^{j}_{h}$ is an approximation obtained via the DG formulation, eq. 22, with the numerical initial conditions $\psi^{j}_{h}(\cdot,0)=\Pi^{k}\psi^{j}_{0}(\cdot)$, then

[TABLE]

Where $\Psi=(\psi^{j})_{j=1}^{N_{\boldsymbol{\Omega}}}$ and $\Psi_{h}=(\psi_{h}^{j})_{j=1}^{N_{\boldsymbol{\Omega}}}$ Here, $\|{\cdot}\|_{-(k+1),\boldsymbol{X}}$ represents a weighted negative-order norm defined in eq. 6 and $C$ depends on upper bounds of $\|{\psi}\|_{k+2,\boldsymbol{X}}$

For the case with constant scattering, we note that this estimate also holds for the divided-differences of the error.

The proof of the negative order norm is given in Section 6.1.2. It relies on the following three estimates, which follow the same structure as in [11], with the additional consideration of the collision operator.

Lemma 10** (Projection Estimate).**

Assume that the same assumptions hold as in Theorem 9, then, defining $\Theta_{M}=-\sum_{j=1}^{N_{\boldsymbol{\Omega}}}\,\omega_{j}{(\eta^{j}_{h},\varphi(0))},$ the following estimate holds for the projection error

[TABLE]

where $C$ depends on $\left(\sum_{j=1}^{N_{\boldsymbol{\Omega}}}\,\omega_{j}{\|{\psi^{j}_{0}}\|_{k+1,\boldsymbol{X}}^{2}}\right)^{1/2}$.

Lemma 11** (Residual estimate).**

Let $\chi^{j}=\Pi^{k}\varphi^{j}$ for each $1\leq j\leq{N_{\boldsymbol{\Omega}}}$. Assume that the same assumptions hold as in Theorem 9, then, defining

[TABLE]

where $\chi^{j}$ in $V_{h}$, we have the following:

[TABLE]

where $C$ depends on the upper bounds of $\|{\psi^{j}_{h}}\|_{k+2,\boldsymbol{X}}$.

Lemma 12** (Consistency estimate).**

Assume that the same assumptions hold as in Theorem 9 and define

[TABLE]

One can show that $\Theta_{C}$ satisfies

[TABLE]

4 Superconvergence analysis

In this section, we outline main steps in the proofs of our superconvergence results.

4.1 Analysis of superconvergence for the one-dimensional steady-state problem

Here, we sketch key steps to prove Theorem 8. Throughout this subsection, we restrict the analysis of the DG solution $\psi^{j}_{h}$ and its error function to the $n-$th element, $K_{n}=[x_{n-1},x_{n}].$ The error on $K_{n}$ is denoted as e_{h}^{j}\bigg{|}_{K_{n}}=e_{n}^{j},\,j=1,\dots,{N_{v}} and $e_{K}=(e_{n}^{1},e_{n}^{2},\dots,e_{n}^{{N_{v}}})$.

We use the same framework for the superconvergence analysis of the linear advection equation as used by Adjerid et al. [4] and match the coefficients of the polynomial expansion of the error equation. A key step in [4] is to utilize the theory for the Runge-Kutta method by viewing the one-dimensional linear advection equation as an ODE. However, unlike the one-direction linear advection equation, the RTE equation involves both left- and right-going particles coupled through an integral term. As a result, the ODE argument in [4] is not applicable for kinetic equations. To bypass this difficulty, we utilize element-by-element mathematical induction.

Before proceeding with the analysis, we note that the exact solution on element $K_{n}$ can be expressed as a Taylor series centered around the element center, $\overline{x_{n}}$:

[TABLE]

where $\zeta$ is the local coordinate mapping. Noting the relation between monomials, $\zeta^{m}$, and Legendre polynomials,

[TABLE]

where

[TABLE]

This can be written as

[TABLE]

This relation can also be seen by using an $L^{2}-$projection of $\psi^{j}(x,v_{j})$ onto the Legendre polynomials:

[TABLE]

Derivation of the error equation on the reference element. Denote the error on element $n$ as

[TABLE]

Then, the error equation on the reference element is given by

[TABLE]

Using an upwind flux, the error equation then becomes

[TABLE]

Polynomial expansion of the error function on the reference element as a series of $\Delta x$. Using assumptions in Theorem 8, and the Taylor series expansion eq. 35, the error in terms of powers of $h$ is

[TABLE]

where $Q_{{m}}(\zeta)$ is a polynomial of degree $m.$ We will show that it is defined as

[TABLE]

Notice that $Q_{{m}}(\zeta)$ is written in terms of the local coordinates and does not depend on the element itself. Further, note that naively applying the expansion in eq. 36 would lead to

[TABLE]

We will see that this is not the case.

Grouping the error equation in terms of powers of $h$ leads to the system:

[TABLE]

for $\tau\in\mathbb{P}^{k}$. In terms of the coefficients, this becomes:

[TABLE]

for all $\tau\in\mathbb{P}^{k}$. Note that for the $S_{N}$ system

[TABLE]

Proof that $Q_{{k+1}}(\zeta)=\widetilde{\beta_{k+1}}R_{k+1}^{\pm}(\zeta)$. To prove that the leading order error term is a Radau polynomial, we set $\tau(\zeta)=L_{q}(\zeta),\,q=0,\dots,p$ in eq. 41. Then, the leading order error is:

[TABLE]

For $q=0$ this simplifies to

[TABLE]

This implies that, for $v_{j}>0$, either $\alpha^{j}_{{n},{k+1}}-\alpha^{j}_{{n-1},{k+1}}=0,\,n=1,\dots,N_{x},\,$ or $Q_{{k+1}}(1)=0$ and for $v_{j}<0$ either $\alpha^{j}_{{n+1},{k+1}}-\alpha^{j}_{{n},{k+1}}=0,\,n=1,\dots,N_{x},\,$ or $Q_{{k+1}}(-1)=0$. If the difference in the coefficients is zero, one can prove by induction that $\alpha^{j}_{{n},{m}}=0$ for all $m,\,n$ since we assume the inflow boundary is exact. This implies that the error is identically zero, which can only occur if our exact solution is a polynomial in the approximation space.

Next, consider $q=1,\dots,k:$

[TABLE]

Noting that

[TABLE]

leads to the conclusion that the quantity on the right is always zero as $q-(2s+1)\leq k-1.$ Hence $Q_{{k+1}}(\zeta)$ is orthogonal to polynomials of degree $\leq k-1$ and, combining the two results above, leads to

[TABLE]

Combining (i) the local Taylor series expansions for both the exact solution and the DG approximation; (ii) the relation between the monomials and the Legendre polynomials, eq. 36; and (iii) Galerkin orthogonality leads to the expression, for $m\geq k+2$,

[TABLE]

Thus proving the second part of our theorem.

Orthogonality of $Q_{{m}}$ for $m\geq k+2$. Consider the case where $m\geq k+2$ and $\tau(\zeta)=L_{q}(\zeta),\,q=0,\dots,k$. In this case, the error equation is given by

[TABLE]

Rearranging the to obtain a relation for the $Q_{m}(1)$ gives

[TABLE]

By Galerin orthogonality, the right side is zero for $q=0,1,\dots,k$. Hence the first non-vanishing term at $\zeta=1$ is $Q_{2(k+1)}$ and our theorem theorem 8 is proven.

4.2 Analysis of superconvergence for time dependent case

In this section, we outline the proof of Theorem theorem 9:

Proof.

Let $\Phi\in(C_{0}^{\infty}(\boldsymbol{X}))^{N}$, Then by (8) and the dual problem definition,

[TABLE]

Notice that for any $\chi^{j}$ in $V_{h}$,

[TABLE]

This allows for the error to be written as

[TABLE]

where

[TABLE]

$\Theta_{M}$, $\Theta_{N}$ and $\Theta_{C}$ are respectively the projection, residual and consistency terms. Using Lemmas 10, 11, and 12 together with the dual estimate, eq. 14 gives our desired estimate. ∎

5 Extracting superconvergence using the SIAC filter

Now that we have proven that higher-order accuracy exists in the negative-order norm, we show how to extract that information via the Smoothness-Increasing Accuracy-Conserving (SIAC) filter.

To illustrate the ability of SIAC to perform on a given data set, it is useful to outline how SIAC works for general data as well as through the error estimates.

Given $\psi^{j}_{h},\,j=1,\dots,N_{\boldsymbol{\Omega}},$ superconvergence can be extracted through convolving with a specially designed kernel, ${K_{H}\left(\cdot\right)}$:

[TABLE]

where $H$ represents the kernel scaling, in this case the uniform mesh size.

We show the reliance on the negative-order norm by decomposing the filtered error into a term that only depends on the number of moments and a term that relies on the error in the negative-order norm:

[TABLE]

where $\|\cdot\|$ is some norm. Here, the choice of $r$ is the number of moments the filter is designed to capture and $s$ is the order of accuracy of the approximation. In this article, the ability to bound the $L^{2}-$norm by the negative-order norm is utilized and $r=2k+1.$

The success of the filter relies on the following results.

Theorem 13** (Bramble and Schatz [6]).**

Let $k\geq 0$. For $T>0$, let $\Psi=(\psi^{n})_{n=1}^{N_{\boldsymbol{\Omega}}}$ be the exact solution to problem (21b), satisfying $\psi^{j}\in C([0,T];H^{k+1}(\boldsymbol{X}))$, $1\leq j\leq N_{\boldsymbol{\Omega}}$. Let $\Omega_{0}+2supp(K_{h}^{2(k+1),k+1}(\boldsymbol{x}))\subset\subset\boldsymbol{X}$ and $\Psi_{h}=(\psi^{j}_{h})_{j=1}^{N_{\boldsymbol{\Omega}}}$ then

[TABLE]

where $C_{\mathrm{P}}$ depends solely on $\Omega_{0}$, $\boldsymbol{X}$, $k$, and it is independent of $h$.

In (49), we used the notation of the divided differences, which are defined as

[TABLE]

where $\mathbf{e}_{i}$ is the unit multi-index whose $i$-th component is $1$ and all others [math].

For any multi-index $\lambda=(\lambda_{1},\lambda_{2},\lambda_{3})$ we set the $\lambda$-th order difference quotient to be

[TABLE]

5.1 SIAC formulation.

The Smoothness-Increasing Accuracy-Conserving (SIAC) kernel is comprised of $r+1$ (scaled) function translates of a given function,

[TABLE]

in this article, central B-Splines are used. ${B_{\text{\tiny T},n}}$ represents $n^{th}$-order central B-spline with knot sequence, T and smoothness $n-2$. The scaling, $H$, is generally tied to the mesh size. The central B-splines are defined through the relations

[TABLE]

where ${\bf T}_{n}$ represents a knot matrix for the $n^{th}$ order spline (i.e. B-spline breaks) [42]. For a symmetric kernel of $2k+1$ B-splines, the general form of the knot matrix is

[TABLE]

Each row of the knot matrix gives the B-Spline breaks of the $\gamma^{th}$ B-spline [42] ($\gamma=1,\dots,{2k+1}$). The $c_{\gamma}$ are weights of the B-splines, which are determined by ensuring that the kernel satisfies consistency plus $r=2k$ moments. We further note that where ${K_{H}\left(\cdot\right)}=\frac{1}{H}{K\left(\frac{\cdot}{H}\right)}$ can be viewed as a normalized probability density function. We note that Mock and Lax [35] introduced the importance of satisfying moment conditions and pre-processing data. This allows for recovering accuracy for discontinuous functions – away from any discontinuities. The pre-processing of data is important for methods not based on Galerkin orthogonality. Further, utilizing a linear combination of B-Splines allows for writing derivatives can be written as divided differences of lower order splines,

[TABLE]

where $\partial_{H}^{\alpha}$ represents the $\alpha^{th}$ divided difference. This ensures that when we pass to the negative order norm, the order of accuracy is not reduced. These ideas were introduced by Bramble and Schatz [6] and Thomeé [47].

Here, we note that, using equally spaced knots for the B-spline filter, the Fourier transform of the SIAC kernel is given by

[TABLE]

As can be seen, the smoothness chosen for the B-splines controls the amount of dissipation and the number of moments controls the accuracy.

There are a few choices for extending the filter to multi-dimensions. One method is via a tensor product:

[TABLE]

However, for computational efficiency, a rotated one-dimensional filter, the Line SIAC (LSIAC) kernel [13], is often used

[TABLE]

In [13], results in two- and three-dimensions are demonstrated. For two-dimensions, filtering is performed along the line $\Gamma(t)=(\overline{x},\overline{y})+\lambda(\cos(\theta),\ \sin(\theta))$, with an angle of rotation $\theta=\tan\left(\frac{\Delta y}{\Delta x}\right)$.

5.2 A note on computation

An illustration of the improved performance from using the Line SIAC kernel in two-dimensions can be seen by considering a kernel consisting of 5 B-splines of order 3 on a structured mesh. To filter one point, 196 two-dimensional integrals are required for the implementation of the tensor product filter, while only 21 one-dimensional integrals are required for LSIAC. Computing these integrals require quadratures that respect both B-Spline breaks and mesh breaks. Computing these breaks is the most costly aspect of the filter. We note that if polynomials of degree $k$ are used in the DG approximation and B-Splines of order $k+1$ are used, $\lceil\frac{2k+1}{2}\rceil$ quadrature points per region are required.

5.3 Post-processing error estimates

Here we present the main theorem for the $L^{2}$ error estimate for the postprocessed solution for the time dependent problem.

Theorem 14**.**

Let $\psi$ be the exact solution to (3), and let us assume it satisfies $\psi\in C([0,T];L^{2}(X;H^{s}(\mathbb{S}^{2})))$. $\psi^{j}\in C([0,T];H^{k+2}(X))$, $1\leq j\leq N_{\boldsymbol{\Omega}}$. if $\psi^{j}_{h}$ is the DG solution to (22) with numerical initial conditions $\psi^{j}_{h}(\cdot,0)=\Pi^{k}\psi^{j}_{0}(\cdot)$, then

[TABLE]

where $C_{\mathrm{T}}$ depends on the upper bound of $\left(\sum_{j=1}^{N_{\boldsymbol{\Omega}}}\,\omega_{j}{\|{\partial_{h}^{\lambda}\psi^{j}}\|_{k+2,X}}^{2}\right)^{1/2}$, for all $|\lambda|\leq k+1$ and the constant $C_{\mathrm{P}}$ in Theorem 13 and $D$ depends on upper bounds of $\|{\psi}\|_{L^{2}(X;H^{s}(\mathbb{S}^{2}))}$.

Proof.

A direct application of triangle inequality gives,

[TABLE]

The first term in the inequality comes from the proof of Theorem 7. To bound the second term, we use the fact that since the scattering $\sigma$ is constant, then for all multi-indices $\lambda$, $\partial_{h}^{\lambda}\psi$ satisfies equation (3). Then by Theorem 9

[TABLE]

with $C_{\lambda}$ depending on upper bounds of $\|{\partial_{h}^{\lambda}\psi^{j}}\|_{k+2,\boldsymbol{X}}$, $1\leq j\leq N_{\boldsymbol{\Omega}}$. Then the conclusion follows applying Theorem 13. ∎

6 Numerical results

Here, we demonstrate the performance of SIAC filter through a series of numerical examples. Due to its superior efficiency discussed in Sec 5.2, we apply the line SIAC filter [13] in all our numerical tests.

The linear system resulting from the discretization of the steady-state problem or implicit time-marching is solved through Source Iteration with Diffusion Synthetic Acceleration (SI-DSA) [4]. A partially consistent DSA strategy is applied as our preconditioner (see [48] and Appendix A of [36] for details). The diffusion equation inside our DSA preconditioner is solved by conjugate gradient method with algebraic multigrid (AMG) preconditioner. We set the stopping criteria of the source iteration as $||\overline{\psi}^{(\ell)}-\overline{\psi}^{(\ell-1)}||<\epsilon_{\textrm{SI}}$ with $\overline{\psi}^{(\ell)}$ be the density in the $\ell$-th iteration. We set $\epsilon_{\textrm{SI}}$ as $10^{-10}$ for $K\leq 2$ and $10^{-11}$ for $K=3$.

We implement our code in the Julia language, leveraging IterativeSolvers.jl package for Krylov solver, AlgebraMultigrid.jl package for AMG preconditioner and MSIAC.jl [14] for the SIAC filter. Numerical tests are performed on a MacBook with Apple M1 chip.

6.1 Accuracy test

We present a series of accuracy tests to demonstrate superconvergence and the computational saving gained to reach the same level of accuracy by applying the SIAC filter.

6.1.1 Steady-state problem

The first problem that we consider the computational domain $[-1,1]^{2}$ with vacuum boundary conditions. A uniform mesh with $N_{x}\times N_{y}$ rectangular elements is applied to partition the computational domain. The mesh size is defined to be $h=\min(1/N_{x},1/N_{y})$. We impose source terms so that the manufactured solution $f(x,y,v_{x},v_{y})=\sin(\pi x)\sin(\pi y)$ is satisfied. We consider two different material properties: (1) constant material, $\sigma_{s}(x,y)=1$ and $\sigma_{a}(x,y)=0$; and (2) variable material, $\sigma_{s}(x,y)=2+\sin(16\pi x)\sin(16\pi y)$ and $\sigma_{a}=0$.

Superconvergence. In Fig. 1, we present $L^{2}$ error of the numerical solution before and after post-processing with the SIAC filter. Before applying the filter, we observe the expected $(k+1)$-th order of accuracy. After applying the SIAC filter, we observe approximately $3.5$ order of accuracy for $Q_{1}$ elements, and $2(k+2)$-th order of accuracy for $k=2,3$. This observation matches the theories given in Sec. 3.

We also observe that the filtered solution is more accurate when the mesh resolution is sufficiently high, while it may not be more accurate on a coarse mesh.

Efficiency gain. In this example, directly applying the SIAC filter to $\overline{\psi}$ always leads to less than $10\%$ more computational time. In Fig. 2, we present the relation between the $L^{2}$ error and the computational time. We observe that the SIAC filter enables us to obtain significantly more accurate results with almost the same computational time as when a more refined mesh is used. Moreover, we observe that, after a break-even point, applying the SIAC filter takes less computational time to reach the desired level accuracy compared to $h$-refinement.

For constant scattering, to reach an error of pproximately $7\times 10^{-5}$ with $Q_{1}$ elements, it takes approximately $15.42$ seconds with $h=N_{y}=20$ when the SIAC filter is applied. It takes approximately $140.06$ seconds with $N_{x}=N_{y}=80$ without the help of SIAC filter to reach the same level of accuracy. In this case, SIAC filter leads to approximately $9$ times acceleration compared to computing on a refined mesh.

For the constant scattering, when $Q_{2}$ element are used, the SIAC filters enables us to obtain $6.63\times 10^{-8}$ in the $L^{2}$ error with approximately $160.71$ seconds and $N_{x}=N_{y}=40$, while only $1.35\times 10^{-7}$ for the $L^{2}$ error is achieved with $621.45$ seconds and $N_{x}=N_{y}=80$ without filtering. In this case, applying SIAC filter permits us to obtain $2$ times more accurate results with only $25\%$ computational time compared to refining the mesh.

Similar observations can be seen for the variable scattering case as well.

6.1.2 Time dependent problem

For the time-dependent case, we consider the computational domain $[-1,1]^{2}$ with vaccum boundary conditions. A uniform mesh with $N_{x}\times N_{y}$ rectangular elements is used with mesh size $h=\min(1/N_{x},1/N_{y})$. The material properties used in this example are $\sigma_{s}(x,y)=1$ and $\sigma_{a}(x,y)=0$. We impose source terms so that $f(t,x,y,v_{x},v_{y})=\exp(-t)\sin(\pi x)\sin(\pi y)$ is an exact solution. We run the simulation from time $t=0$ to $t=0.5$.

Superconvergence. For the time-dependent problem, we test both a BDF$2$ and BDF$3$ schemes. For these methods, the temporal error will become dominant when the SIAC filter is applied and, in order to observe spatial superconvergence, either a higher-order temporal discretization should be used, or the CFL should be modified accordingly. We consider the following discretizations (1) $Q1$-BDF$2$ with a time step size $\Delta t=h$, (2) $Q1$-BDF$3$ with a time step size $\Delta t=h$ and (3) $Q2$-BDF$3$ with $\Delta t=4h^{\frac{5}{3}}$.

We present the relation between the $L^{2}-$error and the mesh size in Fig. 3. From the top left plot, we observe that SIAC filter is able to improve the accuracy of $Q_{1}$-BDF$2$, but the accuracy order is still second order due to the dominance of the temporal error. From the top right result, we observe that the SIAC filter improves the accuracy order of $Q_{1}$-BDF$3$ from second order to third order. From the bottom left figure, slightly higher than $5$-th order accuracy is observed for $Q_{2}$-BDF$3$ with $\Delta t=4h^{\frac{5}{3}}$ when SIAC filter is applied. In summary, the SIAC filter is able to improve the order of accuracy to $2k+1$ with sufficiently accurate temporal discretization.

Efficiency gain. We compare the efficiency of $Q_{1}$-BDF$2$, $Q_{1}$-BDF$3$ without post-processing with the filtered $Q_{1}$-BDF$3$. The time step size is chosen as $\Delta t=h$ for all methods.

Though BDF$3$ involves more vector operations per time step than BDF$2$, both of them only require one linear solve per time step through SI-DSA. Since the linear solve SI-DSA takes significantly longer time than additional vector operations in BDF$3$, the overall computational efficiency of these two time integrators are comparable.

In the bottom right plot of Fig. 3, we present the relation between the computational time, namely $T_{\textrm{comp.}}$, and the $L^{2}-$error. We observe that the $L^{2}-$error scales roughly as $T_{\textrm{comp.}}^{-0.8}$ for $Q_{1}$-BDF$2$ and $Q_{1}$-BDF$3$ without filtering, while it scales as $T_{\textrm{comp.}}^{-1.3}$ for $Q_{1}$-BDF$3$ post-processed using the SIAC filter. Utilizing the same computational time, $Q_{1}$-BDF$3$ with post-processing achieves significantly more accurate results.

Specifically, with $N_{x}=N_{y}=96$, $Q_{1}$-BDF$2$ and $Q_{1}$-BDF$3$ before filtering achieves $2.52\times 10^{-5}$ and $1.71\times 10^{-5}$ $L^{2}-$error with $3004.53$ seconds and approximately $3049.48$ seconds computational time, respectively. By applying SIAC filter, $Q_{1}$-BDF$3$ with $N_{x}=N_{y}=32$ is able to achieve $7.70\times 10^{-6}$ error with only $152.94$ seconds. Consequently, using the SIAC filter and BDF$3$ achieves an approximately $3.26$ times more accurate result with $19.65$ times acceleration over $Q_{1}$-BDF$2$, and a $2.22$ times more accurate result with $19.94$ times acceleration over $Q_{2}$-BDF$3$.

Compared to the steady-state case, the time of post-processing using the SIAC filter is less significant – it is smaller than $5\%$ – since the total number of source iterations is significantly higher in the time-marching case.

6.2 Steady-state with variable scattering

We consider a steady-state problem on the computational domain $[-1,1]^{2}$ with vacuum boundary conditions and a Gaussian source $G(x,y)=\frac{10}{\pi}\exp(-100(x^{2}+y^{2}))$. There is no absorption, while the scattering cross section is defined as

[TABLE]

The configuration of $\sigma_{s}$ and a reference solution generated with $Q_{1}$ DG using $128\times 128$ rectangular mesh in space and CL($80,40$) quadrature rule is presented in Fig. 4.

As discussed in [31, 19], the regularity of RTE is low in the angular space, and the $S_{N}$ method may suffer from significant order reduction in the angular space. With low regularity in the angular space, we are not able to observe superconvergence when applying the SIAC filter due to the dominating angular error. However, as shown in Fig. 4, the SIAC filter still removes non-physical oscillations in the center of computational domain and improves the resolution of the solution on a coarse space-angle mesh with $32\times 32$ elements in the physical space and CL($20,10$).

7 Conclusions

In this paper, we theoretically and computationally present the benefit of applying SIAC filters to the upwind DG method for solving the steady-state and time-dependent RTE.

We have proven $(2k+2)$-th order accuracy for the steady state RTE at the outflow edge and $(k+2)$-th order accuracy at the interior roots of the Radau points on each spatial element, and $(2k+\frac{1}{2})$-th order accuracy for the time-dependent problem with respect to a weighted negative-order norm. Our numerical results validate $(2k+2)$-th order superconvergence for steady state problem and demonstrate $(2k+1)$-th order superconvergence for the time-dependent problem. 2. 2.

Furthermore, we numerically demonstrate that directly applying SIAC filters to low-dimensional macroscopic density can greatly reduce the computational time to reach a desired level of accuracy.

Potential future directions are as follows: (1) For varying cross sections, our current superconvergence proof can be extended by establishing divided difference results for the DG approximation. Though numerically observed superconvergnce, a complete superconvergence proof for spatial cross sections is still an open questions. Additionally, analysis for problems involving nonlinear thermal radiation, multi-energy groups and anisotropic scattering is challenging and worth investigation. (2) Our current analysis mainly focuses on improving convergence order in physical space. However, another important area of exploration is enhancing angular accuracy and mitigating ray effects through post-processing in angular space.

8 Declarations

Acknowledgments

A. Galindo-Olarte gratefully acknowledges support from the Oden Institute for Computational Engineering and Sciences. Z. Peng also acknowledges support from the Hong Kong Research Grants Council (grant ECS-26302724 and GRF-16306825). The work of J. Ryan was funded by the Swedish Research Council (VR grant 2022-03528). Z. Peng and J. Ryan were also funded by internal support from the HKUST-KTH Global Knowledge Network Awards.

Appendix A Analysis of spatial error, theorem 5

Proof.

Summing (22) over all $K\in\mathcal{T}_{h}$ we have that $\psi_{h}$ satisifies

[TABLE]

where

[TABLE]

It is clear that the exact solution satisfies the weak form above, and it is linear. Then the $e^{j}_{h}$ satisfies

[TABLE]

where $\overline{e_{h}}=\frac{1}{m(\mathcal{S})}\sum_{j=1}^{N_{\Omega}}\omega_{j}e^{j}_{h}$.

By setting the test function $\tau_{h}=\xi^{j}_{h}$

[TABLE]

We first the term on the left hand-side of Equation (62) which is equivalent to

[TABLE]

Now just focus on the integral terms

[TABLE]

then the left-hand side of the equality is given by:

[TABLE]

Now let us concentrate on the right hand side, since $\xi^{j}_{h}$ it is orthogonal to any piecewise polynomial of degree $k$:

[TABLE]

then the only two terms that survuve from the right hand side are

[TABLE]

For the first term in the above sum, notice that since the Time derivative commutes with $\Pi^{k}$, then by orthogonality and the fact that $\xi^{j}_{h}\in V_{h}^{k}$,

[TABLE]

for the second term, by the definition of the upwind flux (23) and Lemma 3 gives,

[TABLE]

Finally notice that, by the orthogonality of the $L^{2}$ projection, $\xi^{j}_{h}\bot\eta^{i}_{h}$, then

[TABLE]

we can easily see that

[TABLE]

then combining our computations for the lef-hand and right-hand sides together and computing the weighted sum,

[TABLE]

Hence,

[TABLE]

then by using (68),

[TABLE]

then we end up with the following Gronwall inequality

[TABLE]

Then using the fact that the initial error $\psi^{j}_{h}(x,0)=\Pi^{k}\psi(x,\boldsymbol{\Omega}_{j},0)$,

[TABLE]

This give us the final result. ∎

Appendix B Analysis of angular error, theorem 6

Proof.

Let us define the error in angle $e^{j}(x)=\psi(\boldsymbol{x},\mathbf{\Omega}_{j})-\psi^{j}(\boldsymbol{x})$ and $\mathbf{e}(x)=(e^{j}(\boldsymbol{x}))_{j=1}^{N_{\boldsymbol{\Omega}}}$ Then, $e^{j}$ satisfies

[TABLE]

where $\overline{\mathbf{e}}=\frac{1}{m(\mathcal{S})}\sum_{j=1}^{N_{\boldsymbol{\Omega}}}\omega_{j}e^{j}$.

[TABLE]

if we multiply (78) by $\omega_{j}e^{j}$ and integrate over $X$ and use the periodic boundary conditions in $\boldsymbol{x}$, we have the following error equation,

[TABLE]

Notice that since

[TABLE]

then

[TABLE]

from we can easily see that

[TABLE]

Thus that since the $\omega_{j}>0$, then a double application of Cauchy–Schwarz inequality, first for integrals and then for sums, gives

[TABLE]

An application of Lemma 1 gives

[TABLE]

The conclusion follows using Gronwalls inequality and the fact that $e^{j}(\boldsymbol{x},0)=0.$ ∎