# A priori error analysis of discrete-ordinate weak Galerkin method for
radiative transfer equation

Maneesh Kumar Singh
Department of Mathematics
Imperial College London
SW7 2AZ, London, UK
<EMAIL_ADDRESS>
###### Abstract

This research article discusses a numerical solution of the radiative transfer
equation based on the weak Galerkin finite element method. We discretize the
angular variable by means of the discrete-ordinate method. Then the resulting
semi-discrete hyperbolic system is approximated using the weak Galerkin
method. The stability result for the proposed numerical method is devised. A
_priori_ error analysis is established under the suitable norm. In order to
examine the theoretical results, numerical experiments are carried out.

_Keywords_ Radiative transfer equation $\cdot$ Weak Galerkin method $\cdot$
Discrete-ordinate method $\cdot$ Stability results $\cdot$ Error analysis

##  1 Introduction

The propagation, absorption, and scattering of particles within the
participating media are simulated using the radiative transport equation
(RTE). The application of RTE takes place in numerous scientific fields,
including atmospheric physics, biomedical optics, heat transfer, neutron
transport, optical imaging, radiative transport, etc. Due to the higher-
dimensional RTE, it is difficult to discuss the numerical approximation of the
model problems. More details on the analytical and numerical solution of RTE
in different participating mediums can be found in Chandrasekhar (1960);
Modest (2013) and references therein.

###  1.1 The model problem

This research explores a stationary radiative transfer equation that is
modulated by

$\left\\{\begin{array}[]{ll}\displaystyle\mathbf{s}\cdot\nabla
u(\mathbf{x},\mathbf{s})+\sigma_{as}(\mathbf{x})u(\mathbf{x},\mathbf{s})-\sigma_{s}(\mathbf{x})\mathcal{K}u(\mathbf{x},\mathbf{s})=f(\mathbf{x},\mathbf{s}),\quad(\mathbf{x},\mathbf{s})\in\Omega\times\mathcal{S},\\\\[4.0pt]
u(\mathbf{x},\mathbf{s})=u_{\mbox{in}}(\mathbf{x},\mathbf{s}),\,\,\,(\mathbf{x},\mathbf{s})\in\Gamma^{-},\end{array}\right.$
(1)

where $\Omega$ is a bounded polytopal domain in $\mathbb{R}^{d}(d=1,2,3)$ with
a smooth boundary $\partial\Omega$ and $\mathcal{S}$ is the unit sphere in
$\mathbb{R}^{d}$. The incoming and outgoing boundaries are given by

$\Gamma^{\pm}=\\{(\mathbf{x},\mathbf{s}):\mathbf{x}\in\partial\Omega_{\mathbf{s}}^{\pm},\,\mathbf{s}\in\mathcal{S}\\},$

with

$\partial\Omega_{\mathbf{s}}^{-}=\\{\mathbf{x}\in\partial\Omega:\mathbf{s}\cdot\mathbf{n}(\mathbf{x})<0\\},\quad\partial\Omega_{\mathbf{s}}^{+}=\\{\mathbf{x}\in\partial\Omega:\mathbf{s}\cdot\mathbf{n}(\mathbf{x})\geq
0\\},$

where $\mathbf{n}(\mathbf{x})$ the unit outward normal for
$\mathbf{x}\in\partial\Omega$. Here $\sigma_{as}=\sigma_{a}+\sigma_{s}$,
$\sigma_{a}$ is the absorption coefficient, $\sigma_{s}(\sigma_{s}\geq 0)$ is
the scattering coefficient and $f$ is the source function with $u_{\mbox{in}}$
the external radiative intensity. Also, we assume that
$\sigma_{as},\,\sigma_{s},f$ and $u_{\mbox{in}}$ are smooth functions. And,
the integral operator $\mathcal{K}$ is given by

$\mathcal{K}u(\mathbf{x},\mathbf{s})=\int_{\mathcal{S}}\Phi(\mathbf{x},\mathbf{s},\mathbf{s}^{\prime})u(\mathbf{x},\mathbf{s}^{\prime})d\sigma(\mathbf{s}^{\prime}),\quad\int_{\mathcal{S}}\Phi(\mathbf{x},\mathbf{s},\mathbf{s}^{\prime})d\sigma(\mathbf{s}^{\prime})=1,\,\,\forall\,\mathbf{x}\in\Omega,\,\,\mathbf{s}\in\mathcal{S},$

where $\Phi$ is a non-negative normalized phase function. For example, we
consider Henyey-Greenstein (H-G) phase function which is given by

$\Phi(\mathbf{s},\mathbf{s}^{\prime})=\dfrac{1}{2^{d-1}\pi}\dfrac{1-\eta^{2}}{\left(1+\eta^{2}-2\eta(\mathbf{s}\cdot\mathbf{s}^{\prime})\right)^{3/2}},\,\,d=2,3.$
(2)

where the anisotropy factor $\eta\in(-1,1)$. We assume that data of the model
problem (1) satisfies the following conditions:

$\sigma_{as},\,\sigma_{s}\in
L^{\infty}(\Omega),\quad\sigma_{as}-\sigma_{s}\geq\sigma_{0}>0,\,\,\mbox{in}\,\,\Omega,\,\,\sigma_{0}\in\mathbb{R}.$
(3)

Denote the space $H^{1}_{2}(\Omega\times\mathcal{S})$ by

$H^{1}_{2}(\Omega\times\mathcal{S})=\\{u\in
L^{2}(\Omega\times\mathcal{S}),\,\mathbf{s}\cdot\nabla u\in
L^{2}(\Omega\times\mathcal{S})\\}.$

Using regularity of the data and positivity condition (3), it can be shown
that the model problem has a unique solution $u\in
H^{1}_{2}(\Omega\times\mathcal{S})$. For more details, one may see Agoshkov
(2012).

In literature, the numerical methods for RTE can be categorized into the
stochastic and deterministic approaches. In the stochastic approach, one of
the major numerical techniques is the Monte Carlo simulations. There are many
articles, which used Monte-Carlo simulation, e.g. Gentile (2001); Howell
(1998); Kong et al. (2008); McClarren and Urbatsch (2009). Because of its
iterative nature, the Monte-Carlo approach is not very time-efficient for the
numerical approximation of RTE. In the deterministic approach, some popular
numerical methods are the Spherical harmonics ($P_{N}$)-method and the
discrete-ordinate method (DOM), also called the $S_{N}$-method. The spherical
harmonics method is based on an expansion of the radiation intensity in angles
using spherical harmonics basis functions. Detailed analysis of this method
for RTE is presented in Egger and Schlottbom (2012, 2016, 2019); Sheng and
Wang (2021) and references therein. Whereas in DOM, the basic idea is to use a
finite number of discrete directions to discretize the angular variables of
RTE, which is also our area of interest. DOM reduces the original RTE into a
semi-discrete first-order hyperbolic system. The detailed description of this
numerical approach can be found in Koch et al. (1995); Lewis and Miller
(1984); Kanschat et al. (2008).

For spatial discretization, due to the hyperbolic nature of RTE, the
discontinuous Galerkin (DG) methods are one of the popular numerical
techniques. DG methods for differential equations have been extensively
studied in the last few decades, see e.g. Arnold et al. (2002); Brezzi et al.
(2004); Houston et al. (2002) and references therein. DG methods for spatial
discretization of RTE based on the DOM method are discussed in Han et al.
(2010) and also error analysis is presented. Later, a stabilized DG method
along with DOM is used for numerical simulation in Wang et al. (2016). A more
sophisticated numerical method comprises a multigrid preconditioner for
$S_{N}$-DG approximation is presented in Kanschat and Ragusa (2014). Recently,
the upwind DG method for scaled RTE is discussed in Sheng and Hauck (2021). In
Peng et al. (2020), a low-rank numerical method for time-dependent RTE is
established. A sparse grid DOM discontinuous Galerkin method for the RTE is
investigated in Huang and Yu (2021). In Ganesan and Singh (2021), we have
proposed an operator-splitting based finite element method for the time-
dependent RTE, where the key idea is to split the RTE model problem concerning
the internal (angular) and the external (spatial) directions, resulting in a
transient transport problem and a time-dependent integro-differential
equation. A numerical simulation based on A reduced basis method for the
stationary RTE model is proposed in Peng et al. (2022).

Discrete-ordinate DG type finite element methods for RTE are developed in Han
et al. (2010); Wang et al. (2016). In this article, we discuss the weak
Galerkin (WG) method for the spatial discretization of the resulting semi-
discrete system. One of the highlights of the WG method is that it has the
advantage of the general mesh, which provides a high degree of versatility in
numerical simulation, especially in the complex geometry domain. Unlike the DG
method, the WG method is parameter-free while it has the same degrees of
freedom as DG possesses. More details of this method can be found in Lin et
al. (2018); Mu et al. (2015); Wang and Ye (2014) and references therein. In
recent times, WG finite element method is used for different models, like the
Cahn-Hilliard model is solved in Wang et al. (2019a) and Biot’s consolidation
model is analyzed in Hu et al. (2018). In this work, we will use a modified
version of the weak Galerkin method introduced in Wang et al. (2014); Gao et
al. (2020) for spatial discretization. Recently, an explicit WG method is used
for numerical simulation of first-order hyperbolic systems in Zhang and Zhang
(2022). Therefore, it is interesting to present a numerical simulation of the
RTE model via the WG method.

This article is the first one, which explores the implementation of the WG
method for discretizing the spatial direction of RTE. The proposed numerical
method is the combination of DOM for the angular variable and the WG method
for the spatial variable. We call the resulting method, the discrete-ordinate
weak Galerkin (DOWG) method. The stability result of the DOWG method is
derived. A _priori_ error estimate is obtained by combining the projection and
discretization errors.

The rest of the research is framed as follows: In Section 2, we define the DOM
for the model problem (1). Then, we introduce the WG method for the resulting
semi-discrete hyperbolic system and this section ends with the stability
result of the proposed DOWG method. Section 3 discuss the error estimate of
the DOWG method by combining the error estimate associated with both the
numerical methods, _i.e_ , DOM, and WG method. The numerical experiments that
support our theoretical findings are covered in Section 4 and the article is
summarized in Section 5.

Notations: We have used the standard notations for Sobolev space and norm/semi
norm Brenner and Scott (2008) thorough out this article. And $C$, sometimes
with a subscript, denotes a generic positive constant, the value of which can
be different in different occurrences.

##  2 Discrete-ordinate WG finite element method

In this section, to solve the RTE, a discrete-ordinate weak Galerkin finite
element method is presented (1). The numerical simulations are developed in
two phases: First, the discrete-ordinate approach is applied to discretize the
angular variable, which deduces the model problem (1) into a system of
hyperbolic PDE. Then the resulting system of differential equations is further
discretized by the weak Galerkin method.

###  2.1 Discrete-ordinate method

Before deriving the discrete-ordinate method, we briefly discuss the
quadrature rule which is instrumental for the proposal of the discrete method.
For any continuous function $F$, we have

$\int_{\mathcal{S}}F(\mathbf{s})d\sigma(\mathbf{s})\approx\sum_{m=1}^{M}\omega_{m}F(\mathbf{s}_{m}),\quad\omega_{m}>0,\,\,\
\mathbf{s}_{m}\in\mathcal{S},\,1\leq m\leq M.$ (4)

For $\mathcal{S}\in\mathbb{R}^{2}$, the quadrature rule (4) is expressed as

$\int_{\mathcal{S}}F(\mathbf{s})d\sigma(\mathbf{s})=\int_{0}^{2\pi}\overline{F}(\theta)d\theta\approx\sum_{m=1}^{M}\omega_{m}\overline{F}(\theta_{m}),\quad\omega_{m}>0,\,\,\
\theta_{m}\in[0,2\pi],\,1\leq m\leq M,$ (5)

where the angular variable is expressed in the spherical coordinate system.
Afterwards, any standard quadrature formula can be used in (5). More details
can be found in any standard book of numerical analysis, for example, see Süli
and Mayers (2003). In particular, we have applied the composite trapezoidal
rule in (5). Hence, we use $h_{\theta}=2\pi/M$ and weight
$\omega_{0}=\omega_{M}=\dfrac{h_{\theta}}{2}$ and $\omega_{m}=h_{\theta}$ for
$1\leq m\leq M-1$.

For $\mathcal{S}\in\mathbb{R}^{3}$, the quadrature rule (4) is expressed as

$\begin{array}[]{ll}\displaystyle\int_{\mathcal{S}}F(\mathbf{s})d\sigma(\mathbf{s})&\displaystyle=\int_{0}^{2\pi}\int_{0}^{\pi}\overline{F}(\theta,\psi)\sin\theta
d\theta d\psi,\quad(\theta,\psi)\in[0,\pi]\times[0,2\pi]\\\\[12.0pt]
&\displaystyle\approx\dfrac{\pi}{m}\sum_{j=1}^{2m}\sum_{i=1}^{m}\omega_{i}\overline{F}(\theta_{i},\psi_{j}),\quad\omega_{i}>0,\end{array}$
(6)

where $\\{\theta_{i}\\}$ are chosen so that $\\{\cos\theta_{i}\\}$ and
$\\{\omega_{i}\\}$ are the Gauss-Legendre nodes and weights on $[-1,1]$. The
points $\\{\psi_{j}\\}$ are evenly spaced on $[0,2\pi]$ with a spacing of
$\pi/m$.

By following the quadrature error estimates Hesse and Sloan (2006) for $3D$
angular discretization, we have

$\left|\int_{\mathcal{S}}F(\mathbf{s})d\sigma(\mathbf{s})-\sum_{m=1}^{M}\omega_{m}F(\mathbf{s}_{m})\right|\leq
C_{r}n^{-r}\|F\|_{r,\mathcal{S}},\,\,F\in H^{r}(\mathcal{S}),\,\,r\geq 1,$ (7)

where $C_{r}$ is a positive constant depending only on $r$, and $n$ denotes
the degree of precision of the quadrature. And, the respective norm
$\|\cdot\|_{r,\mathcal{S}}$ is given by

$\|v\|^{2}_{r,\mathcal{S}}:=\int_{\mathcal{S}}\sum_{k=0}^{r}\left|\nabla^{(k)}_{\mathbf{s}}v\right|^{2},\qquad\nabla^{(k)}_{\mathbf{s}}v:=\left\\{\nabla^{(\alpha)}_{\mathbf{s}}v\big{|}|\alpha|=k\right\\},$
(8)

where $\alpha$ is a multi-index and $|\nabla^{(k)}_{\mathbf{s}}v|$ is given by

$|\nabla^{(k)}_{\mathbf{s}}v|=\sum_{|\alpha|=k}\left(|\nabla^{(\alpha)}_{\mathbf{s}}v|^{2}\right)^{1/2}.$

By using the standard numerical quadrature formula, the integral operator
$\mathcal{K}$ is approximated as follows:

$\mathcal{K}_{d}u(\mathbf{x},\mathbf{s}_{m})=\sum_{m=1}^{M}\omega_{m}\Phi(\mathbf{x},\mathbf{s}\cdot\mathbf{s}_{m})u(\mathbf{x},\mathbf{s}_{m}),\quad\omega_{m}>0,\,\,\mbox{and}\,\,\mathbf{s}_{m}\in\mathcal{S},\,\,0\leq
m\leq M.$ (9)

By employing the approximated operator (9) and discrete-ordinate method, the
continuous problem (1) is deduced as follows:

$\left\\{\begin{array}[]{ll}\displaystyle\mathbf{s}_{m}\cdot\nabla
u^{m}(\mathbf{x})+\sigma_{as}(\mathbf{x})u^{m}(\mathbf{x})-\sigma_{s}\sum_{\ell=1}^{M}\omega_{\ell}\Phi(\mathbf{x},\mathbf{s}_{m}\cdot\mathbf{s}_{\ell})u^{\ell}(\mathbf{x})=f^{m}(\mathbf{x}),\quad\mbox{in}\,\,\Omega,\\\\[12.0pt]
u^{m}(\mathbf{x})=u_{\mbox{in}}^{m}(\mathbf{x}),\quad\mbox{on}\,\,\partial\Omega_{m}^{-}\equiv\partial\Omega_{\mathbf{s}_{m}}^{-},\qquad\qquad\qquad
1\leq m\leq M,\end{array}\right.$ (10)

where
$f^{m}(\mathbf{x})=f(\mathbf{x},\mathbf{s}_{m}),\,u_{\mbox{in}}^{m}=u_{\mbox{in}}(\mathbf{x},\mathbf{s}_{m})$.

Note that, the exact solution $u^{m}=u^{m}(\mathbf{x})$ of the discrete system
(10) is an approximation of the exact solution $u$ of the continuous problem
(1) at points $(\mathbf{x},\mathbf{s}_{m})$ for $m=1,2,\ldots,M$.

Before going into the details of the error estimate, we assume the following
regularity conditions for the solution of the discrete method (10): In the
$2D$ setting, we assume that exact solution $u$ of the continuous problem (1)
and the phase function satisfy following regularity property:

$\left\\{\begin{array}[]{ll}u^{m}\in H^{p}(\Omega),\quad\mbox{for }\,\,p\geq
1,\quad 1\leq m\leq M,\\\\[6.0pt] \displaystyle u\in
L^{2}(\Omega,\mathcal{C}^{2}(\mathcal{S}))\qquad\mbox{and}\quad\sup_{(\mathbf{x},\mathbf{s})\in\Omega\times\mathcal{S}}\|\Phi^{\prime\prime}(\mathbf{x},\mathbf{s}\cdot)\|_{0,\infty,\mathcal{S}}<C_{\Phi},\quad\mbox{for
}\,\,r\geq
1,\quad\mbox{where}\,\,\Phi^{\prime\prime}(\mathbf{x},t)=\dfrac{d^{2}\Phi}{dt^{2}}.\end{array}\right.$
(11)

And, in the $3D$ setting, we consider

$\left\\{\begin{array}[]{ll}u^{m}\in H^{p}(\Omega),\quad\mbox{for }\,\,p\geq
1,\quad 1\leq m\leq M,\\\\[6.0pt] \displaystyle u\in
L^{2}(\Omega,H^{r}(\mathcal{S}))\qquad\mbox{and}\quad\sup_{(\mathbf{x},\mathbf{s})\in\Omega\times\mathcal{S}}\|\Phi(\mathbf{x},\mathbf{s}\cdot)\|_{r,\mathcal{S}}<C_{\Phi_{*}},\quad\mbox{for
}\,\,r\geq 1.\end{array}\right.$ (12)

where the Bochner spaces $L^{2}(\Omega,\mathcal{C}^{2}(\mathcal{S}))$ and
$L^{2}(\Omega,H^{r}(\mathcal{S}))$ are given by

$\displaystyle
L^{2}(\Omega,\mathcal{C}^{2}(\mathcal{S}))=\left\\{v:\Omega\to\mathcal{C}^{2}(\mathcal{S})):\int_{\Omega}\|v\|_{\mathcal{C}^{2}(\mathcal{S}))}<\infty\right\\}$
$\displaystyle L^{2}(\Omega,H^{r}(\mathcal{S})=\left\\{v:\Omega\to
H^{r}(\mathcal{S}):\int_{\Omega}\|v\|_{r,\mathcal{S}}<\infty\right\\}.$

For simplification of error analysis, denote

$b(\mathbf{x})=\max_{1\leq\ell\leq
M}\sum_{\ell=1}^{M}\omega_{m}\Phi(\mathbf{x},\mathbf{s}_{\ell}\cdot\mathbf{s}_{m}).$

By using the normalized property of the scattering phase function with the
estimate (7), for any fixed $\mathbf{x}\in\Omega$, one can write

$\left|1-\sum_{\ell=1}^{M}\omega_{m}\Phi(\mathbf{s}_{m}\cdot\mathbf{s}_{\ell})\right|\leq
C_{r}n^{-r}\|\Phi(\mathbf{s}_{m}\cdot)\|_{r,\mathcal{S}},\quad r\geq 1,$ (13)

Form (3) and (13), it is clear that

$\sigma_{as}(\mathbf{x})-b(\mathbf{x})\sigma_{s}(\mathbf{x})\geq
C_{*}\quad\mbox{in}\,\,\Omega,\quad C_{*}\,\,\mbox{a positive constant}.$ (14)

This condition has much significance while deriving the stability result of
the numerical approximation.

In the next few lines, we discuss the error estimate by means of the discrete
ordinate method (10).

###### Theorem 2.1

Let $u$ and $u^{m}$ be the solutions of (1) and (10), respectively.

In the $2D$ case,

$\left(\sum_{\ell=1}^{M}\omega_{\ell}\sum_{T\in\mathcal{T}_{h}}\|u^{m}(\mathbf{x})-u(\mathbf{x},\mathbf{s}_{m})\|^{2}_{T}\right)^{1/2}=\mathcal{O}(h_{\theta}^{2}),$
(15)

where $h_{\theta}$ is sufficiently small. In the $3D$ case,

$\left(\sum_{\ell=1}^{M}\omega_{\ell}\sum_{T\in\mathcal{T}_{h}}\|u^{m}(\mathbf{x})-u(\mathbf{x},\mathbf{s}_{m})\|^{2}_{T}\right)^{1/2}\leq
C\,n^{-r}\left(\int_{\Omega}\|u\|^{2}_{r,\mathcal{S}}\,\mbox{d}\mathbf{x}\right)^{1/2},$
(16)

where $C$ is positive constant depending on $r$ and the phase function $\Phi$.

Proof. By using the regularity properties (11)-(12), the proof of this theorem
is followed from (Wang et al., 2016, Theorem 4.3, Theorem 4.4).

###  2.2 Weak Galerkin finite element method

Let $\mathcal{T}_{h}$ be a mesh partition of the spatial domain $\Omega$
consisting of polyhedrons in three dimensions, satisfying a set of shape-
regular conditions Mu et al. (2015). Let $\mathcal{E}_{h}$ is the set of all
edges in $\mathcal{T}_{h}$ with
$\mathcal{E}^{0}_{h}=\mathcal{E}_{h}\setminus\partial\Omega_{h}$ is the set of
all interior edges. And $\partial\Omega_{h}$ be the subset of
$\mathcal{E}_{h}$ of all edges on $\partial\Omega$. For every element
$T\in\mathcal{T}_{h}$, we denote by $h_{T}$ the diameter of $T$ and the mesh
size $h=\max_{T\in\mathcal{T}_{h}}h_{T}$ for $\mathcal{T}_{h}$.

Let $T_{1}$ and $T_{2}$ are the mesh elements sharing a common edge $e$. We
define the average of the scalar valued $v$ and vector valued functions
$\mathbf{v}$ on edge $e$ as follows:

$\\{v\\}_{e}=\left\\{\begin{array}[]{ll}\dfrac{1}{2}(v|_{T_{1}}+v|_{T_{2}}),&e\in\mathcal{E}^{0}_{h},\\\\[8.0pt]
v,&e\in\partial\Omega,\end{array}\right.\quad\\{\mathbf{v}\\}_{e}=\left\\{\begin{array}[]{ll}\dfrac{1}{2}(\mathbf{v}|_{T_{1}}+\mathbf{v}|_{T_{2}}),&e\in\mathcal{E}^{0}_{h},\\\\[8.0pt]
\mathbf{v},&e\in\partial\Omega.\end{array}\right.$

For a fixed angular direction $\mathbf{s}_{m}$, the incoming and outgoing
boundaries of the mesh element $T\in\mathcal{T}_{h}$ is given by

$\partial^{-}T=\\{\mathbf{x}\in\partial
T:\mathbf{s}_{m}\cdot\mathbf{n}(\mathbf{x})<0\\},\qquad\partial^{+}T=\\{\mathbf{x}\in\partial
T:\mathbf{s}_{m}\cdot\mathbf{n}(\mathbf{x})\geq 0\\}.$

For convenience, we introduce the following notations:

$\displaystyle(\phi,\psi)_{\mathcal{T}_{h}}=\sum_{T\in\mathcal{T}_{h}}\int_{T}\phi\psi\mbox{d}\mathbf{x},\,\,\langle\phi,\psi\rangle_{\partial\mathcal{T}_{h}}=\sum_{T\in\mathcal{T}_{h}}\int_{\partial
T}\phi\psi\,\mbox{ds},\,\,\langle\phi,\psi\rangle_{\partial\Omega_{h}}=\sum_{T\in\mathcal{T}_{h}}\int_{\partial
T\cap\partial\Omega_{h}}\phi\psi\,\mbox{ds}.$

Given an integer $k\geq 1$, the finite element space is defined as follows

$V_{h}=\\{v_{h}\in
L^{2}(\Omega):v_{h}|_{T}\in\mathbb{P}_{k}(T),\,\,\forall\,T\in\mathcal{T}_{h}\\},$

where $\mathbb{P}_{k}(T)$ is the space of polynomials on set $T$ with degree
no more than $k$.

###### Definition 2.2

For any function $v^{m}_{h}\in V_{h}$, on each element $T$, we define the weak
gradient operator $\nabla_{w}v_{h}^{m}|_{T}\in[\mathbb{P}_{k-1}(T)]^{d}$ as
follows:

$(\nabla_{w}v_{h}^{m},\mathbf{q}_{h})_{T}=-(v^{m}_{h},\nabla\cdot\mathbf{q}_{h})_{T}+\big{\langle}\\{v^{m}_{h}\\},\mathbf{q}_{h}\cdot\mathbf{n}\big{\rangle}_{\partial
T},\quad\forall\mathbf{q}_{h}\in[\mathbb{P}_{k-1}(T)]^{d}.$

And we define the weak divergence
$\mathbf{s}_{m}\cdot\nabla_{w}v^{m}_{h}\in\mathbb{P}_{k}(T)$ associated with
$\mathbf{s}_{m}\cdot\nabla v_{h}^{m}$ on $T$, as the unique polynomial
satisfying

$(\mathbf{s}_{m}\cdot\nabla_{w}v^{m}_{h},w_{h})=-(v^{m}_{h},\mathbf{s}_{m}\cdot\nabla
w_{h})+\big{\langle}\\{v^{m}_{h}\\},\mathbf{s}_{m}\cdot\mathbf{n}\,w_{h}\big{\rangle}_{\partial
T},\quad\forall w_{h}\in\mathbb{P}_{k}(T).$

where $\mathbf{n}$ is the unit outward normal vector of $\partial T$.

Next, we introduce the projection operators, which will be used in the error
analysis of the weak Galerkin method. For each element $T\in\mathcal{T}_{h}$,
let $\mathcal{Q}_{T}:L^{2}(T)\to\mathbb{P}_{k}(T)$ and
$\Pi_{T}:[L^{2}(T)]^{d}\to[\mathbb{P}_{k-1}(T)]^{d}$ be the
$L^{2}$-projections on the associated local polynomial spaces respectively.
For any mesh element $T$, we define the global projections $\mathcal{Q}_{h}$
and $\Pi_{h}$ element-wise by

$\begin{array}[]{ll}(\mathcal{Q}_{h}v)|_{T}=\mathcal{Q}_{T}(v|_{T}),&\forall\,v\in
L^{2}(T),\\\\[4.0pt]
(\Pi_{h}\mathbf{q})|_{T}=\Pi_{T}(\mathbf{q}|_{T}),&\forall\,\mathbf{q}\in[L^{2}(T)]^{d}.\end{array}$

The finite dimensional space $\mathbf{V}_{h}$ is defined as
$\mathbf{V}_{h}=(V_{h})^{M}$ for the semi-discrete problem (10) and an element
$\mathbf{v}_{h}\in\mathbf{V}_{h}$ is given by
$\mathbf{v}_{h}:=\\{v_{h}^{m}\\}_{m=1}^{M}$.

The WG finite element method for the semi-discrete problem (10) reads: Find
$u_{h}^{m}\in V_{h}$ such that

$\displaystyle\mathcal{A}_{c}(u_{h}^{m},v_{h}^{m})+\mathcal{A}_{st}(u_{h}^{m},v_{h}^{m})=\mathcal{F}^{\mathbf{s}}_{c}(v_{h}^{m}),\quad\forall\,v_{h}^{m}\in
V_{h},$ (17)

where $\mathcal{A}_{c}(\cdot,\cdot),\,\mathcal{A}_{st}(\cdot,\cdot)$ and
$\mathcal{F}^{\mathbf{s}}_{c}(\cdot)$ are given by

$\begin{array}[]{ll}\displaystyle\mathcal{A}_{c}(u_{h}^{m},v_{h}^{m})=\left(\mathbf{s}_{m}\cdot\nabla
u_{h}^{m},v_{h}^{m}\right)+\bigg{(}\sigma_{as}u_{h}^{m}-\sigma_{s}\sum_{\ell=1}^{M}\omega_{\ell}\Phi(\mathbf{x},\mathbf{s}_{m}\cdot\mathbf{s}_{\ell})u_{h}^{\ell},v_{h}^{m}\bigg{)}-\langle\mathbf{s}_{m}\cdot\mathbf{n}\,u_{h}^{m},v_{h}^{m}\rangle_{\partial\Omega^{-}_{h}},\\\\[16.0pt]
\mathcal{A}_{st}(u_{h}^{m},v_{h}^{m})=\displaystyle\sum_{T\in\mathcal{T}_{h}}\left\langle\mathbf{s}_{m}\cdot\mathbf{n}\big{(}u^{m}_{h}-\\{u^{m}_{h}\\}\big{)},v^{m}_{h}-\\{v^{m}_{h}\\}\right\rangle_{\partial^{+}T},\\\\[16.0pt]
\mathcal{F}^{\mathbf{s}}_{c}(v_{h}^{m})=(f^{m},v_{h}^{m})_{\mathcal{T}_{h}}-\langle\mathbf{s}_{m}\cdot\mathbf{n}\,u_{\mbox{in}}^{m},v_{h}^{m}\rangle_{\partial\Omega^{-}_{h}}.\end{array}$

Note that, we have imposed the boundary condition in the weak sense.

The global formulation of the DOWG method (17) is expressed as: Find
$\mathbf{u}_{h}=\\{u_{h}^{m}\\}\in\mathbf{V}_{h}$ such that

$\begin{array}[]{ll}\displaystyle\sum_{m=1}^{M}\omega_{m}\,\mathcal{A}_{c}(u_{h}^{m},v^{m}_{h})+\sum_{m=1}^{M}\omega_{m}\,\mathcal{A}_{st}(u_{h}^{m},v^{m}_{h})=\displaystyle\sum_{m=1}^{M}\omega_{m}\mathcal{F}^{\mathbf{s}}_{c}(v_{h}^{m}).\end{array}$

In the simplified version, the DOWG method is presented as follows:

$\left\\{\begin{array}[]{ll}\mbox{Find}\,\,\mathbf{u}_{h}\in\mathbf{V}_{h}\,\,\mbox{such
that}\\\\[6.0pt]
\mathcal{A}_{WG}(\mathbf{u}_{h},\mathbf{v}_{h})=\mathcal{F}_{c}(\mathbf{v}_{h}),\,\,\forall\,\mathbf{v}_{h}\in\mathbf{V}_{h},\end{array}\right.$
(18)

where the bilinear form $\mathcal{A}_{WG}(\cdot,\cdot)$ and the linear
functional $\mathcal{F}^{\mathbf{s}}_{c}(\cdot)$ are given by

$\mathcal{A}_{WG}(\mathbf{u}_{h},\mathbf{v}_{h})=\sum_{m=1}^{M}\omega_{m}\big{(}\mathcal{A}_{c}(u_{h}^{m},v^{m}_{h})+\mathcal{A}_{st}(u_{h}^{m},v_{h}^{m})\big{)},\quad\mathcal{F}_{c}(\mathbf{v}_{h})=\sum_{m=1}^{M}\omega_{m}\mathcal{F}^{\mathbf{s}}_{c}(v_{h}^{m}).$

For any $\mathbf{v}_{h}\in\mathbf{V}_{h}$ , the norm
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
is defined by

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}=\sum_{m=1}^{M}\omega_{m}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v^{m}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}^{2}$
(19)

where the norm term
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}$
is given by

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}^{2}=\sum_{T\in\mathcal{T}_{h}}\left(\|v^{m}_{h}\|^{2}_{T}+\left\||\mathbf{s}_{m}\cdot\mathbf{n}|^{1/2}\big{(}v_{h}^{m}-\\{v_{h}^{m}\\}\big{)}\right\|^{2}_{\partial
T}\right)+\||\mathbf{s}_{m}\cdot\mathbf{n}|^{1/2}v_{h}^{m}\|^{2}_{\partial\Omega_{h}}.$

##  3 Error analysis

In this section, we discuss the error estimate of the DOWG method. Firstly, we
start with the stability result of the proposed numerical method (18) and
using this stability result, we derive the error estimate of the numerical
solution.

###  3.1 Stability result

Now, we discuss the coercivity of the bilinear form
$\mathcal{A}_{WG}(\cdot,\cdot)$ in the following lemma.

###### Lemma 3.1

For sufficiently small $h$, we have

$\alpha{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}\leq\,\mathcal{A}_{WG}(\mathbf{v}_{h},\mathbf{v}_{h}),\,\,\forall\,\mathbf{v}_{h}\in\mathbf{V}_{h},$

where $\alpha$ is a positive constant.

Proof. By using (Gao et al., 2020, Lemma 3.2), it can be easily deduce that

$\left(\mathbf{s}_{m}\cdot\nabla_{w}v^{m}_{h},v^{m}_{h}\right)_{\mathcal{T}_{h}}=-\dfrac{1}{2}\left\langle\mathbf{s}_{m}\cdot\mathbf{n}\big{(}v_{h}^{m}-\\{v_{h}^{m}\\}\big{)},v_{h}^{m}-\\{v_{h}^{m}\\}\right\rangle_{\partial\mathcal{T}_{h}}+\dfrac{1}{2}\langle\mathbf{s}_{m}\cdot\mathbf{n}\,v_{h}^{m},v_{h}^{m}\rangle_{\partial\Omega_{h}}.$

Then, we have

$\begin{array}[]{ll}\mathcal{A}_{c}(v_{h}^{m},v_{h}^{m})+\mathcal{A}_{st}(v_{h}^{m},v_{h}^{m})&=\displaystyle\sum_{T\in\mathcal{T}_{h}}\left(\sigma_{as}v^{m}_{h},v^{m}_{h}\right)_{T}+I+\dfrac{1}{2}\langle\mathbf{s}_{m}\cdot\mathbf{n}\,v_{h}^{m},v_{h}^{m}\rangle_{\partial\Omega_{h}}\\\\[16.0pt]
&\quad-\dfrac{1}{2}\left\langle\mathbf{s}_{m}\cdot\mathbf{n}\big{(}v_{h}^{m}-\\{v_{h}^{m}\\}\big{)},v_{h}^{m}-\\{v_{h}^{m}\\}\right\rangle_{\partial\mathcal{T}_{h}}\\\\[16.0pt]
&\displaystyle\quad+\sum_{T\in\mathcal{T}_{h}}\left\langle\mathbf{s}_{m}\cdot\mathbf{n}\big{(}u^{m}_{h}-\\{u^{m}_{h}\\}\big{)},v^{m}_{h}-\\{v^{m}_{h}\\}\right\rangle_{\partial^{+}T},\end{array}$

where $I$ is given by

$I=-\sum_{T\in\mathcal{T}_{h}}\left(\sigma_{s}\sum_{\ell=1}^{M}\omega_{\ell}\Phi(\mathbf{x},\mathbf{s}_{m}\cdot\mathbf{s}_{\ell})v^{\ell}_{h},v^{m}_{h}\right)_{T},$

and it satisfies

$|I|\leq\sum_{T\in\mathcal{T}_{h}}\left(b\sigma_{s}v^{m}_{h},v^{m}_{h}\right)_{T}.$

By using the assumption on the absorption and scattering coefficients (14), we
obtain that

$\begin{array}[]{ll}\mathcal{A}_{c}(v_{h}^{m},v_{h}^{m})+\mathcal{A}_{st}(v_{h}^{m},v_{h}^{m})&\geq\displaystyle\left((\sigma_{as}-b\sigma_{s})v^{m}_{h},v^{m}_{h}\right)+\dfrac{1}{2}\langle\mathbf{s}_{m}\cdot\mathbf{n}\,v_{h}^{m},v_{h}^{m}\rangle_{\partial\Omega_{h}}\\\\[14.0pt]
&\quad+\dfrac{1}{2}\left\langle\mathbf{s}_{m}\cdot\mathbf{n}\big{(}v_{h}^{m}-\\{v_{h}^{m}\\}\big{)},v_{h}^{m}-\\{v_{h}^{m}\\}\right\rangle_{\partial\mathcal{T}_{h}}\\\\[14.0pt]
&\geq\left(C_{*}v^{m}_{h},v^{m}_{h}\right)+\dfrac{1}{2}\langle\mathbf{s}_{m}\cdot\mathbf{n}\,v_{h}^{m},v_{h}^{m}\rangle_{\partial\Omega_{h}}\\\\[14.0pt]
&\quad+\dfrac{1}{2}\left\langle\mathbf{s}_{m}\cdot\mathbf{n}\big{(}v_{h}^{m}-\\{v_{h}^{m}\\}\big{)},v_{h}^{m}-\\{v_{h}^{m}\\}\right\rangle_{\partial\mathcal{T}_{h}}.\end{array}$

Therefore, the bilinear form satisfies

$\begin{array}[]{ll}\mathcal{A}_{WG}(\mathbf{v}_{h},\mathbf{v}_{h})&\displaystyle\geq\alpha\sum_{m=1}^{M}\omega_{m}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v^{m}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}^{2}\\\\[12.0pt]
&=\alpha{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2},\qquad\alpha=\min\\{C_{*},\frac{1}{2}\\}.\end{array}$

This completes the proof.

The direct consequence of the above lemma is the uniqueness result.

###### Theorem 3.2

There exist a unique solution of the DOWG method (18).

Next, we state the trace inequality: For any function $\phi\in H^{1}(T)$, the
following trace inequality holds true (see Wang and Ye (2014)):

$\|\phi\|^{2}_{\partial T}\leq
C\left(h_{T}^{-1}\|\phi\|^{2}_{T}+h_{T}\|\nabla\phi\|^{2}_{T}\right).$ (20)

###  3.2 Error estimate

Next, we discuss the error estimate of the DOWG method. Firstly, we establish
an error equation associated with the discrete problem (18).

Define the components $\mathcal{L}_{j}(\cdot,\cdot)$ as follows:

$\displaystyle\mathcal{L}_{1}(u^{m},v_{h}^{m})$
$\displaystyle=(u^{m}-\mathcal{Q}_{h}u^{m},\mathbf{s}_{m}\cdot\nabla
v_{h}^{m})_{\mathcal{T}_{h}},$ $\displaystyle\mathcal{L}_{2}(u^{m},v_{h}^{m})$
$\displaystyle=-\big{\langle}u^{m}-\\{\mathcal{Q}_{h}u^{m}\\},\mathbf{s}_{m}\cdot\mathbf{n}(v_{h}^{m}-\\{v_{h}^{m}\\})\big{\rangle}_{\partial\mathcal{T}_{h}},$
(21) $\displaystyle\mathcal{L}_{3}(u^{m},v_{h}^{m})$
$\displaystyle=\displaystyle-\left(\sigma_{as}(u^{m}-\mathcal{Q}_{h}u^{m})-\sigma_{s}\sum_{\ell=1}^{M}\omega_{\ell}\Phi(\mathbf{x},\mathbf{s}_{m}\cdot\mathbf{s}_{\ell})(u^{\ell}-\mathcal{Q}_{h}u^{\ell}),v_{h}^{m}\right)_{\mathcal{T}_{h}},$
$\displaystyle\mathcal{L}_{4}(u^{m},v_{h}^{m})$
$\displaystyle=\big{\langle}\mathbf{s}_{m}\cdot\mathbf{n}(u^{m}-\mathcal{Q}_{h}u^{m}),v_{h}^{m}\big{\rangle}_{\partial\Omega^{-}_{h}}.$

###### Lemma 3.3

The error term $\mathbf{e}^{m}=(\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h})$
satisfies:

$\mathcal{A}_{WG}(\mathbf{e}^{m},\mathbf{v}_{h})=\mathbf{L}(\mathbf{u},\mathbf{v}_{h}),$
(22)

with

$\mathbf{L}(\mathbf{u},\mathbf{v}_{h})=\sum_{\ell=1}^{M}\omega_{m}\mathcal{L}(u^{m},v_{h}^{m}),\qquad\mathcal{L}(u^{m},v^{m}_{h})=\sum_{j=1}^{4}\mathcal{L}_{j}(u^{m},v_{h}^{m})+\mathcal{A}_{st}(\mathcal{Q}_{h}u^{m},v_{h}^{m}),$

where $\mathcal{L}_{j}(\cdot,\cdot)$ are given in (3.2).

Proof. By using the integration by parts and definition (2.2), we deduce that

$\begin{array}[]{ll}\left(\mathbf{s}_{m}\cdot\nabla
u^{m},v^{m}_{h}\right)_{\mathcal{T}_{h}}&=-(u^{m},\mathbf{s}_{m}\cdot\nabla
v_{h}^{m})_{\mathcal{T}_{h}}+\langle
u^{m},\mathbf{s}_{m}\cdot\mathbf{n}v_{h}^{m}\rangle_{\partial\mathcal{T}_{h}}\\\\[8.0pt]
&=-(\mathcal{Q}_{h}u^{m},\mathbf{s}_{m}\cdot\nabla
v_{h}^{m})_{\mathcal{T}_{h}}-(u^{m}-\mathcal{Q}_{h}u^{m},\mathbf{s}_{m}\cdot\nabla
v_{h}^{m})_{\mathcal{T}_{h}}\\\\[8.0pt]
&\quad+\langle\\{\mathcal{Q}_{h}u^{m}\\},\mathbf{s}_{m}\cdot\mathbf{n}v_{h}^{m}\rangle_{\partial\mathcal{T}_{h}}+\left\langle
u^{m}-\\{\mathcal{Q}_{h}u^{m}\\},\mathbf{s}_{m}\cdot\mathbf{n}(v_{h}^{m}-\\{v_{h}^{m}\\})\right\rangle_{\partial\mathcal{T}_{h}}\\\\[8.0pt]
&=\left(\mathbf{s}_{m}\cdot\nabla_{w}(\mathcal{Q}_{h}u^{m}),v^{m}_{h}\right)_{\mathcal{T}_{h}}-\mathcal{L}_{1}(u^{m},v_{h}^{m})-\mathcal{L}_{2}(u^{m},v_{h}^{m}).\end{array}$
(23)

Next, one can obtain that

$\begin{array}[]{ll}\displaystyle\left(\sigma_{as}u^{m}-\sigma_{s}\sum_{\ell=1}^{M}\omega_{\ell}\Phi(\mathbf{x},\mathbf{s}_{m}\cdot\mathbf{s}_{\ell})u^{\ell},v^{m}_{h}\right)_{\mathcal{T}_{h}}&=\displaystyle\left(\sigma_{as}\mathcal{Q}_{h}u^{m}-\sigma_{s}\sum_{\ell=1}^{M}\omega_{\ell}\Phi(\mathbf{x},\mathbf{s}_{m}\cdot\mathbf{s}_{\ell})\mathcal{Q}_{h}u^{\ell},v^{m}_{h}\right)_{\mathcal{T}_{h}}\\\\[16.0pt]
&\qquad-\mathcal{L}_{3}(u^{m},v_{h}^{m}).\end{array}$ (24)

Also, it can be easily deduced that

$-\langle\mathbf{s}_{m}\cdot\mathbf{n}\,u^{m},v_{h}^{m}\rangle_{\partial\Omega^{-}_{h}}=-\langle\mathbf{s}_{m}\cdot\mathbf{n}\,\mathcal{Q}_{h}u^{m},v_{h}^{m}\rangle_{\partial\Omega^{-}_{h}}-\mathcal{L}_{4}(u^{m},v_{h}^{m}).$
(25)

By using (23)-(25), we get

$\mathcal{A}_{c}(\mathcal{Q}_{h}u^{m},v_{h}^{m})-\sum_{j=1}^{4}\mathcal{L}_{j}(u^{m},v_{h}^{m})=(f^{m},v_{h}^{m})_{\mathcal{T}_{h}}-\langle\mathbf{s}_{m}\cdot\mathbf{n}\,u_{\mbox{in}}^{m},v_{h}^{m}\rangle_{\partial\Omega^{-}_{h}}.$

Adding $\mathcal{A}_{st}(\mathcal{Q}_{h}u^{m},v_{h}^{m})$ both side of the
above equation, we get

$\mathcal{A}_{c}(\mathcal{Q}_{h}u^{m},v_{h}^{m})+\mathcal{A}_{st}(\mathcal{Q}_{h}u^{m},v_{h}^{m})=(f^{m},v_{h}^{m})_{\mathcal{T}_{h}}-\langle\mathbf{s}_{m}\cdot\mathbf{n}\,u_{\mbox{in}}^{m},v_{h}^{m}\rangle_{\partial\Omega^{-}_{h}}+\mathbf{L}(u^{m},v_{h}^{m}).$
(26)

Then, one can easily write the above equation into the global formulation as
follows:

$\begin{array}[]{ll}\mathcal{A}_{WG}(\mathcal{Q}_{h}\mathbf{u},\mathbf{v}_{h})=\mathcal{F}_{c}(\mathbf{v}_{h})+\mathbf{L}(\mathbf{u},\mathbf{v}_{h}).\end{array}$
(27)

By subtracting the discrete problem (18) from (27), we obtain the required
error equation.

In the next lemma, we derive bounds for each term of
$\mathbf{L}(\mathbf{u},\mathbf{v}_{h})$. Then, we use the coercivity of the
DOWG method to find the estimates of error term
$(\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h})$.

###### Lemma 3.4

Let $u^{m}$ be the solution of the discrete problem (10). For $v^{m}_{h}\in
V_{h}$, we have the following estimates:

$\begin{array}[]{ll}|\mathcal{L}_{1}(u^{m},v_{h}^{m})|&\leq
Ch^{k+1}\|u^{m}\|_{k+1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}},\\\\[8.0pt]
|\mathcal{L}_{2}(u^{m},v_{h}^{m})|&\leq
Ch^{k+1/2}\|u^{m}\|_{k+1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}},\\\\[8.0pt]
|\mathcal{L}_{3}(u^{m},v_{h}^{m})|&\leq
Ch^{k+1}\|u^{m}\|_{k+1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}},\\\\[8.0pt]
|\mathcal{L}_{4}(u^{m},v_{h}^{m})|&\leq
Ch^{k+1/2}\|u^{m}\|_{k+1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}},\\\\[8.0pt]
|\mathcal{A}_{st}(\mathcal{Q}_{h}u^{m},v_{h}^{m})|&\leq
Ch^{k+1/2}\|u\|_{k+1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}.\end{array}$

Proof. By using $L^{2}$-projection error estimates given in (Lin et al., 2018,
Lemma 3.6), we deduce that

$\begin{array}[]{ll}|\mathcal{L}_{1}(u^{m},v_{h}^{m})|=\displaystyle\left|\sum_{T\in\mathcal{T}_{h}}\left(u^{m}-\mathcal{Q}_{h}u^{m},\mathbf{s}_{m}\cdot\nabla
v_{h}^{m}\right)_{T}\right|\leq
Ch^{k+1}\|u^{m}\|_{k+1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}.\end{array}$

Using Cauchy-Schwarz inequality, we have

$\begin{array}[]{ll}|\mathcal{L}_{2}(u^{m},v_{h}^{m})|&=\displaystyle\left|\sum_{T\in\mathcal{T}_{h}}\left\langle
u^{m}-\\{\mathcal{Q}_{h}u^{m}\\},\mathbf{s}_{m}\cdot\mathbf{n}(v_{h}^{m}-\\{v_{h}^{m}\\})\right\rangle_{\partial
T}\right|\\\\[16.0pt]
&\leq\displaystyle\sum_{T\in\mathcal{T}_{h}}\|u^{m}-\\{\mathcal{Q}_{h}u^{m}\\}\|_{\partial
T}\||\mathbf{s}_{m}\cdot\mathbf{n}|^{1/2}(v_{h}^{m}-\\{v_{h}^{m}\\})\|_{\partial
T}\\\\[16.0pt] &\leq
Ch^{k+1/2}\|u^{m}\|_{k+1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}.\end{array}$

In a similar way $\mathcal{L}_{3}(u^{m},v_{h}^{m})$ is bounded by

$\begin{array}[]{ll}|\mathcal{L}_{3}(u^{m},v_{h}^{m})|&=\displaystyle\sum_{T\in\mathcal{T}_{h}}\left|\left(\sigma_{as}(u^{m}-\mathcal{Q}_{h}u^{m})-\sigma_{s}\sum_{\ell=1}^{M}\omega_{\ell}\Phi(\mathbf{x},\mathbf{s}_{m}\cdot\mathbf{s}_{\ell})(u^{\ell}-\mathcal{Q}_{h}u^{\ell}),v^{m}_{h}\right)_{T}\right|\\\\[16.0pt]
&\leq\displaystyle
C\sum_{T\in\mathcal{T}_{h}}\|u^{m}-\mathcal{Q}_{h}u^{m}\|_{T}\|v^{m}_{h}\|_{T}\leq
Ch^{k+1}\|u^{m}\|_{k+1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}},\end{array}$

where the assumption (14) is used. Next, we get

$\begin{array}[]{ll}|\mathcal{L}_{4}(u^{m},v_{h}^{m})|&\leq
C\|u^{m}-\mathcal{Q}_{h}u^{m}\|_{\partial\Omega^{-}_{h}}\||\mathbf{s}_{m}\cdot\mathbf{n}|^{1/2}v_{h}^{m}\|_{\partial\Omega^{-}_{h}}\\\\[8.0pt]
&\leq
Ch^{k+1/2}\|u^{m}\|_{k+1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}.\end{array}$

At last, the stabilization term can be estimated by

$|\mathcal{A}_{st}(\mathcal{Q}_{h}u^{m},v^{m})|\leq
Ch^{k+1/2}\|u^{m}\|_{k+1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}.$

This completes the proof.

###### Theorem 3.5

Let $u^{m}$ and $u^{m}_{h}$ be the solutions to the problem (10) and the
discrete WG scheme (18), respectively. Then, the error term
$\mathbf{e}^{m}=\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h}$ satisfies:

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq
Ch^{k+1/2}\|\mathbf{u}\|_{k+1}.$ (28)

Proof. Using the stability result from Lemma 3.1, we have

$\alpha{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}\leq\mathcal{A}_{WG}(\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h},\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h}).$
(29)

By employing Lemma 3.4, one can obtain that

$\begin{array}[]{ll}\mathcal{A}_{WG}(\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h},\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h})&=\mathbf{L}(\mathbf{u},\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h})\\\\[6.0pt]
&\leq
Ch^{k+1/2}\|\mathbf{u}\|_{k+1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\\\\[6.0pt]
&\leq
Ch^{2k+1}\|\mathbf{u}\|^{2}_{k+1}+\dfrac{\alpha}{2}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathcal{Q}_{h}\mathbf{u}-\mathbf{u}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}.\end{array}$
(30)

A straightforward computation provides the desired result (28).

The projection error is discussed here. By using the trace inequality (20) and
result form (Lin et al., 2018, Lemma 3.6), it can be easily verified that

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u^{m}-\mathcal{Q}_{h}u^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}\leq
Ch^{k+1/2}\|u^{m}\|_{k+1}.$ (31)

Next, we combine the above results to deduce the error estimate associated
with weak Galerkin approximation in the following theorem.

###### Theorem 3.6

Let $\mathbf{u}_{h}=(u^{m}_{h})$ be the numerical approximation of
$\mathbf{u}=(u^{m})$. Then, the following estimates

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{u}-\mathbf{u}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq
Ch^{k+1/2}\|u^{m}\|_{k+1}.$ (32)

Proof. Combing the error estimates (28) and (31), we obtain that

$\begin{array}[]{ll}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u^{m}-u^{m}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}&={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u^{m}-\mathcal{Q}_{h}u^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathcal{Q}_{h}u^{m}-u^{m}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{T}_{h}}\\\\[8.0pt]
&\leq Ch^{k+1/2}\|u^{m}\|_{k+1}.\end{array}$ (33)

Then, use the norm definition to obtain (32). This completes the proof.

Finally, we discuss the estimate for the fully discrete error
$e^{m}_{h}=(u-u_{h}^{m})$. By adding estimates from Theorem 2.1 and Theorem
3.6, the fully discrete error term $e^{m}_{h}$ satisfies:

In the $2D$ model setting,

$\left(\sum_{m=1}^{M}\omega_{m}\sum_{T\in\mathcal{T}_{h}}\|e^{m}_{h}\|^{2}_{T}\right)^{1/2}\leq
Ch^{k+1/2}\|u^{m}\|_{k+1}+\mathcal{O}(h_{\theta}^{2}).$ (34)

In the $3D$ model setting,

$\left(\sum_{m=1}^{M}\omega_{m}\sum_{T\in\mathcal{T}_{h}}\|e^{m}_{h}\|^{2}_{T}\right)^{1/2}\leq
C_{1}h^{k+1/2}\|u^{m}\|_{k+1}+C_{2}\,n^{-r}\left(\int_{\Omega}\|u\|^{2}_{r,\mathcal{S}}\,\mbox{d}\mathbf{x}\right)^{1/2}.$
(35)

##  4 Numerical discussion

In this section, we consider a few test examples employing manufacturing
methods to verify the theoretical estimations provided in the earlier
sections. For the theoretical validation, we focus the numerical experiments
only on the two-dimension models to lower the computational cost. We have
followed the numerical construction from Wang et al. (2016) and the software
library deal.II Wang et al. (2019b) is used for error simulation.

For numerical validation, we discuss error $e_{h}$ in the following norm

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|e_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{1}=\left(\sum_{\ell=1}^{M}\omega_{m}\sum_{T\in\mathcal{T}_{h}}\|e^{m}_{h}\|^{2}_{T}\right)^{1/2},$

and convergence order $eoc_{h}$ is given by

$eoc_{h}=\log_{2}\left(\dfrac{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|e_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{1}}{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|e_{h/2}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{1}}\right).$

###  4.1 Numerical examples

Let $\Omega=(0,1)^{2}$. We consider the following examples of the radiative
transfer equation (1):

###### Example 4.1

The H-G phase function (2) with $\eta=0.5$ and the exact solution is

$u(\mathbf{x},\mathbf{s})=\sin(\pi x_{1})\sin(\pi x_{2})$

.

###### Example 4.2

The scattering phase function $\Phi(t)=\dfrac{(2+t)}{4\pi}$ and the exact
solution is

$u(\mathbf{x},\mathbf{s})=e^{-ax_{1}-bx_{2}}(1+c\cos(\theta)),\quad\theta\in[0,2\pi],$

where $a,b=\dfrac{\sigma_{a}}{3}$ and
$c=\dfrac{\sigma_{a}}{\sigma_{a}+6\sigma_{s}}$.

For both test problems, the source functions $f$ satisfy the radiative
transfer equation (1). Furthermore, the value of absorption and scattering
coefficients are taken as $\sigma_{as}=2,\,\,\sigma_{s}=1/2$ for all the test
problems.

In spatial discretization, the first level grid consists of four squares.
Then, each square is equally decomposed into $4$ squares, to get the next
level grid. For the numerical experiments, $Q_{1}$, the bilinear finite
elements on quadrilaterals and $Q_{2}$, the biquadratic finite elements on
quadrilaterals are chosen and error estimates with convergence orders are
derived. In both cases, the mesh width in angular direction is taken as
$h_{\theta}=\pi/10$.

In table 1, the error estimate, for example, 4.1 is presented whereas the
error estimate for Example 4.2 in table 2. Also, the order of convergence
($eoc_{h}$) is calculated for the error estimates. From tables 1, 2, one can
verify the theoretical error estimate of $\mathcal{O}(h^{k+1/2})$ with $Q_{1}$
and $Q_{2}$-finite elements. In particular, $eoc_{h}$ is close to $2.5$ rather
than $3$ for $Q_{2}$-finite elements. This reduction happens due to the
composite trapezoidal formula which is of second order.

Table 1:  Discretization error $e_{h}$ and order of convergence $eoc_{h}$ for Example 4.1.  Level |  $Q_{1}$-element |  |  $Q_{2}$-element |
---|---|---|---|---
| ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|e_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{1}$ | $eoc_{h}$ | ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|e_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{1}$ | $eoc_{h}$
3 | 8.5643e-02 | – | 6.4527e-02 | –
4 | 2.7626e-02 | 1.63 | 1.2936e-02 | 2.32
5 | 8.4227e-03 | 1.71 | 2.4672e-03 | 2.39
6 | 2.4413e-03 | 1.78 | 4.4057e-04 | 2.48
7 | 6.7628e-04 | 1.85 | 7.4519e-05 | 2.56
Table 2:  Discretization error $e_{h}$ and order of convergence $eoc_{h}$ for Example 4.2.  Level |  $Q_{1}$-element |  |  $Q_{2}$-element |
---|---|---|---|---
| ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|e_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{1}$ | $eoc_{h}$ | ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|e_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{1}$ | $eoc_{h}$
3 | 7.5954e-02 | – | 5.7826e-02 | –
4 | 2.5150e-02 | 1.59 | 1.1892e-02 | 2.28
5 | 8.0352e-03 | 1.64 | 2.3037-03 | 2.36
6 | 2.4497e-03 | 1.72 | 4.2676e-04 | 2.43
7 | 6.9794e-03 | 1.81 | 7.4895e-05 | 2.51

##  5 Concluding remarks

In this article, we have proposed a weak Galerkin finite element method for
the numerical solution of the radiative transfer equation. The discrete-
ordinate weak Galerkin method has been used to approximate the angular and
spatial variables. For the resulting numerical method, the stability result is
established. In the specific norm, the error estimate is devised for the
numerical solutions under the suitable regularity assumptions. In the end, the
numerical experiments are presented for theoretical verification.

## References

  * Chandrasekhar [1960] S. Chandrasekhar.  _Radiative Transfer_.  Dover Publications, New York, 1960.
  * Modest [2013] M. F. Modest.  _Radiative Heat Transfer_.  Academic press, 2013.
  * Agoshkov [2012] V. Agoshkov.  _Boundary Value Problems for Transport Equations_.  Springer Science & Business Media, 2012.
  * Gentile [2001] N. A. Gentile.  Implicit monte carlo diffusion–an acceleration method for monte carlo time-dependent radiative transfer simulations.  _J. Comput. Phys._ , 172(2):543–571, 2001.
  * Howell [1998] J. R. Howell.  The monte-carlo method in radiative heat transfer.  _J. Heat Transfer_ , 120:547–560, 1998.
  * Kong et al. [2008] R. Kong, M. Ambrose, and J. Spanier.  Efficient, automated monte carlo methods for radiation transport.  _J. Comput. Phys._ , 227(22):9463–9476, 2008\.
  * McClarren and Urbatsch [2009] R. G. McClarren and T. J Urbatsch.  A modified implicit monte carlo method for time-dependent radiative transfer with adaptive material coupling.  _J. Comput. Phys._ , 228(16):5669–5686, 2009\.
  * Egger and Schlottbom [2012] Herbert Egger and Matthias Schlottbom.  A mixed variational framework for the radiative transfer equation.  _Math. Models Methods Appl. Sci._ , 22(03):1150014, 2012.
  * Egger and Schlottbom [2016] Herbert Egger and Matthias Schlottbom.  A class of Galerkin schemes for time-dependent radiative transfer.  _SIAM J. Numer. Anal._ , 54(6):3577–3599, 2016\.
  * Egger and Schlottbom [2019] Herbert Egger and Matthias Schlottbom.  A perfectly matched layer approach for $P_{N}$-approximations in radiative transfer.  _SIAM J. Numer. Anal._ , 57(5):2166–2188, 2019\.
  * Sheng and Wang [2021] Qiwei Sheng and Cheng Wang.  A spherical harmonic discontinuous Galerkin method for radiative transfer equations with vacuum boundary conditions.  _Journal of Scientific Computing_ , 88(1):1–35, 2021.
  * Koch et al. [1995] R Koch, W Krebs, S Wittig, and R Viskanta.  Discrete ordinates quadrature schemes for multidimensional radiative transfer.  _J. Quant. Spectrosc. Radiat. Transfer._ , 53(4):353–372, 1995.
  * Lewis and Miller [1984] E. E. Lewis and W. F. Miller.  _Computational Methods of Neutron Transport_.  John Wiley & Sons, 1984.
  * Kanschat et al. [2008] Guido Kanschat, Erik Meinköhn, Rolf Rannacher, and Rainer Wehrse.  _Numerical Methods in Multidimensional Radiative Transfer_.  Springer-Verlag, Berlin, 2008.
  * Arnold et al. [2002] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini.  Unified analysis of discontinuous Galerkin methods for elliptic problems.  _SIAM J. Numer. Anal._ , 39(5):1749–1779, 2002\.
  * Brezzi et al. [2004] F. Brezzi, L. D. Marini, and E. Süli.  Discontinuous Galerkin methods for first-order hyperbolic problems.  _Math. Models Methods Appl. Sci._ , 14(12):1893–1903, 2004.
  * Houston et al. [2002] P. Houston, C. Schwab, and E. Süli.  Discontinuous hp-finite element methods for advection-diffusion-reaction problems.  _SIAM J. Numer. Anal._ , 39(6):2133–2163, 2002\.
  * Han et al. [2010] W. Han, H. Jianguo, and A. E. Joseph.  Discrete-ordinate discontinuous Galerkin methods for solving the radiative transfer equation.  _SIAM J. Sci. Comput._ , 32(2):477–497, 2010\.
  * Wang et al. [2016] C. Wang, Q. Sheng, and W. Han.  A discrete-ordinate discontinuous-streamline diffusion method for the radiative transfer equation.  _Commun. Comput. Phys._ , 20:1443–1465, 2016.
  * Kanschat and Ragusa [2014] Guido Kanschat and J. C. Ragusa.  A robust multigrid preconditioner for $S_{N}$-DG approximation of monochromatic, isotropic radiation transport problems.  _SIAM J. Sci. Comput._ , 36(5):2326–2345, 2014\.
  * Sheng and Hauck [2021] Qiwei Sheng and Cory Hauck.  Uniform convergence of an upwind discontinuous Galerkin method for solving scaled discrete-ordinate radiative transfer equations with isotropic scattering.  _Mathematics of Computation_ , 90(332):2645–2669, 2021.
  * Peng et al. [2020] Zhuogang Peng, R. G. McClarren, and Martin Frank.  A low-rank method for two-dimensional time-dependent radiation transport calculations.  _J. Comput. Phys._ , page 109735, 2020.
  * Huang and Yu [2021] Jianguo Huang and Yue Yu.  A sparse grid discrete ordinate discontinuous Galerkin method for the radiative transfer equation.  _Communications in Computational Physics_ , 30(4):1009–1036, 2021.
  * Ganesan and Singh [2021] Sashikumaar Ganesan and Maneesh Kumar Singh.  An operator-splitting finite element method for the numerical solution of radiative transfer equation.  _arXiv preprint arXiv:2112.07949_ , 2021.
  * Peng et al. [2022] Zhichao Peng, Yanlai Chen, Yingda Cheng, and Fengyan Li.  A reduced basis method for radiative transfer equation.  _Journal of Scientific Computing_ , 91(1):1–27, 2022.
  * Lin et al. [2018] R. Lin, X. Ye, S. Zhang, and P. Zhu.  A weak Galerkin finite element method for singularly perturbed convection-diffusion-reaction problems.  _SIAM J. Numer. Anal._ , 56(3):1482–1497, 2018\.
  * Mu et al. [2015] L. Mu, J. Wang, and X. Ye.  Weak Galerkin finite element methods on polytopal meshes.  _International Journal of Numerical Analysis & Modeling_, 12(1), 2015.
  * Wang and Ye [2014] J. Wang and X. Ye.  A weak Galerkin mixed finite element method for second order elliptic problems.  _Math. Comput._ , 83(289):2101–2126, 2014.
  * Wang et al. [2019a] Junping Wang, Qilong Zhai, Ran Zhang, and Shangyou Zhang.  A weak Galerkin finite element scheme for the Cahn-Hilliard equation.  _Mathematics of Computation_ , 88(315):211–235, 2019a.
  * Hu et al. [2018] Xiaozhe Hu, Lin Mu, and Xiu Ye.  Weak Galerkin method for the Biot’s consolidation model.  _Computers & Mathematics with Applications_, 75(6):2017–2030, 2018.
  * Wang et al. [2014] X. Wang, N. S. Malluwawadu, F. Gao, and T.C. McMillan.  A modified weak Galerkin finite element method.  _J. Comput. Appl. Math._ , 271:319–327, 2014.
  * Gao et al. [2020] F. Gao, S. Zhang, and P. Zhu.  Modified weak Galerkin method with weakly imposed boundary condition for convection-dominated diffusion equations.  _App. Numer. Math._ , 157:490–504, 2020.
  * Zhang and Zhang [2022] Tie Zhang and Shangyou Zhang.  An explicit weak Galerkin method for solving the first order hyperbolic systems.  _Journal of Computational and Applied Mathematics_ , 412:114311, 2022.
  * Brenner and Scott [2008] S Brenner and R. Scott.  _The Mathematical Theory of Finite Element Methods_.  Texts in Applied Mathematics. Springer, 2008.
  * Süli and Mayers [2003] Endre Süli and David F Mayers.  _An Introduction to Numerical Analysis_.  Cambridge university press, 2003.
  * Hesse and Sloan [2006] Kerstin Hesse and Ian H Sloan.  Cubature over the sphere $\mathcal{S}^{2}$ in Sobolev spaces of arbitrary order.  _Journal of Approximation Theory_ , 141(2):118–133, 2006.
  * Wang et al. [2019b] Zhuoran Wang, Graham Harper, Patrick O’Leary, Jiangguo Liu, and Simon Tavener.  deal.II Implementation of a weak Galerkin finite element solver for darcy flow.  In _International Conference on Computational Science_ , pages 495–509. Springer, 2019b.