A discrete-ordinate weak Galerkin method for radiative transfer equation

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.


Introduction
The propagation, absorption, and scattering of particles within the participating media are simulated using the radiative transfer 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.The RTE is an integro-differential equation that depends on a set of optical parameters (index of refraction, absorption, scattering, and scattering function) that describe the medium, further information in [1].Schuster first devised the radiative transfer theory to analyse light propagation in foggy atmospheres [2].Chandrasekhar [3] contributed greatly to the theory of radiative transfer and one of his main contributions was the creation of the vector equation of transfer for polarised light, focusing on Rayleigh scattering.In recent years, this theory has received significant attention for its applications in combustion system [4], remote sensing [5,6], optical tomography [7,8], photo-bioreactors [9] and solar energy harvesting [10,11].A comprehensive overview of the analytical and numerical solution of RTE in different participating mediums can be found in [12] and references therein.
where the spatial boundaries Ω ± s are where n(x) the unit outward normal for x ∈ ∂Ω.
The physical meaning of the medium parameters that determine the RTE is explained in the following few lines.
• Absorption: The absolute amount of absorption is directly proportional to the magnitude of the incident energy and the distance the beam travels through the medium.
• Scattering: or "out-scattering" (away from the direction s under consideration), is very similar to absorption, i.e., a part of the incoming intensity is removed from the direction of propagation, s ′ .
• Total attenuation: The total attenuation of the intensity in a pencil of rays is summed up by the contribution of both absorption and scattering.It is also known as "extinction".
The equation, σ t = σ a + σ s , describes how the overall attenuation σ t is influenced by both absorption σ a and scattering σ s .In simpler terms, it breaks down the total impact on the medium into the portions caused by absorption and scattering.
The functions f and u in are the source and external radiative intensity.The integral operator K in (1.1) is given by where Φ is a non-negative normalized phase function, which describes the scattering property of the participating medium.The precise form of the phase function is usually unknown for applications.In this study, we consider a benchmark choice the Henyey-Greenstein (H-G) phase function which is given by where the anisotropy factor η ∈ (−1, 1).
We assume that data of the model problem (1.1) satisfies the following conditions: where C(S) is a continuous space over S. Next, we denote the space H 1 2 (Ω × S) by with the norm .
By using (1.3), it can be shown that the model problem has a unique solution u ∈ H 1 2 (Ω × S).For more details, one may see [13,Theorem 1].
It is challenging to address the numerical approximation of the model problems because of the higherdimensional RTE.In literature, the numerical methods for RTE can be categorized into the stochastic and deterministic approaches.In the stochastic approach, one of the major computational techniques is the Monte Carlo simulations.There are various articles, that implemented the Monte-Carlo simulation, for example, [14,15,16,17].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 the expansion of radiation intensity in angles using spherical harmonics basis functions.Detailed analysis of these methods for RTE can be found in [18,19,20,13] and references therein.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 semidiscrete first-order hyperbolic system.The detailed description of this numerical approach can be found in [21,22,23].
Owing to RTE's hyperbolic nature, one of the most widely used numerical methods for spatial discretization in the context of finite element schemes is the discontinuous Galerkin (DG) approach.DG methods for differential mathematical models have been extensively studied in the last few decades, see e.g.[24,25,26] and references therein.Recently, the local discontinuous Galerkin (LDG) finite element method for the fractional Allen-Cahn equation is studied in [27].In [28], for the numerical simulation of the nonlinear fractional hyperbolic wave model, a second-order time discrete algorithm with a fast time two-mesh mixed finite element scheme is discussed.DG methods for spatial discretization of RTE based on the DOM method are discussed in [29] and also error analysis is presented.Later, a stabilized DG method along with DOM is used for numerical simulation in [30].A more sophisticated numerical method comprises a multigrid preconditioner for S N -DG approximation is presented in [31].Recently, the upwind DG method for scaled RTE is discussed in [32].In [33], a low-rank numerical method for time-dependent RTE is established.A sparse grid DOM discontinuous Galerkin method for the RTE is investigated in [34].In [35], 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 [36].
Discrete-ordinate DG type finite element methods for RTE are developed in [29,30].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 [37,38,39] and references therein.In the past, the WG finite element method has been used to solve several models.For example, the Cahn-Hilliard model is solved in [40], and Biot's consolidation model is investigated in [41].Recently, an explicit WG method has been used for numerical simulation of first-order hyperbolic systems in [42].Inspired by these works, we explore the possibility of using the WG method for the RTE model, where we will use a modified version of the weak Galerkin method introduced in [43,44] for spatial discretization.This is the first article that focuses on discretizing the spatial direction of RTE utilizing the WG methods.
The proposed numerical method combines DOM for the angular variable and the WG method for the spatial variable.The resulting technique is called the discrete-ordinate weak Galerkin method or DOWG method for short.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.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 discusses 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 [45] throughout this article.And C, sometimes with a subscript, denotes a generic positive constant, the value of which can be different in different occurrences.

Discrete-ordinate WG finite element method
In this section, to solve the RTE, a discrete-ordinate weak Galerkin finite element method is presented (1.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.1) into a system of hyperbolic PDE.
Then, the resulting system of differential equations is further discretized by the weak Galerkin method.

Discrete-ordinate method
Before discussing the discrete-ordinate method, we briefly go through the quadrature rule, which is instrumental for the proposal of the discrete method.For any continuous function F , we can write due to the quadrature approximation For S ∈ R 2 , the quadrature rule (2.1) is expressed as where the angular variable is written in the spherical coordinate form.Afterwards, the standard quadrature formula is used.More details can be found in any standard numerical analysis book, for example, see [46].
In this paper, we have applied the composite trapezoidal quadrature formula in (2.2).To achieve this, we adopted the angular spacing parameter h θ = 2π/M with weights Similarly, for S ∈ R 3 , the following quadrature rule is used in (2.1) to obtain where {θ i } are chosen so that {cos θ i } and {ω i } are the Gauss-Legendre nodes and weights on [−1, 1].The discrete points {ψ j } are evenly spaced on [0, 2π] with a spacing of π/m.By following the quadrature error estimates [47], we have where C r is a positive constant depending only on r, and n denotes the degree of precision of the quadrature.
Here, norm • r,S in (2.4) where α is a multi-index and |∇ (k) . By using the normalized property of the scattering phase function with the estimate (2.4), for any fixed x ∈ Ω, we get Using the above numerical quadrature formulas, the integral operator K is approximated in the following unified way: (2.7)

Semi-discrete problem
By employing the approximation (2.7), the continuous problem (1.1) is discretized in each angular direction s m as follows: where and the incoming and outgoing boundaries are given by Γ ± m = {(x, s m ) : x ∈ ∂Ω ± s }.Note that, the exact solution u m = u m (x) of the semi-discrete system (2.8) is an approximation of the model solution u of the continuous problem (1.1) at points (x, s m ) for m = 0, 2, . . ., M .
Before going into the details of the error estimate, consider the following regularity conditions for the solution of the discrete method (2.8):In the 2D setting, we assume that the exact solution u of the continuous problem (1.1) and the phase function satisfy and sup And, in the 3D setting and sup where the Bochner spaces L 2 (Ω, C 2 (S)) and L 2 (Ω, H r (S)) are given by For simplification of the presentation, we introduce Then, from (1.3) and (2.6), one can deduce that This condition will be utilized while deriving the stability and convergence result of the numerical approximation.
In the next few lines, we discuss the error estimate using the discrete-ordinate method (2.8).
Theorem 2.1.Let u and u m are the solutions of the model problem (1.1) and the semi-discrete problem (2.8), respectively.In the 2D case, where h θ is sufficiently small.And, in the 3D case, we have where C is positive constant depending on r and the phase function Φ.

Weak Galerkin finite element method
Let T h be a mesh partition of the spatial domain Ω consisting of elements that are closed and simply connected polygons in R 2 (polyhedrons in R 3 ) satisfying a set of shape-regular conditions [38].Let E h is the set of all edges/faces in T h with E 0 h = E h \ ∂Ω h is the set of all interior boundaries (faces for d = 3 or edges for d = 2).And ∂Ω h be the subset of E h of all boundaries on ∂Ω.For every element T ∈ T h , we denote by h T the diameter of T and the mesh size h = max T ∈T h h T for T h .
For a fixed angular direction s m , the incoming and outgoing boundaries of the mesh element T ∈ T h is given by 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 v on edge e as follows: And, define the jump of v on e by For convenience, we introduce the following notations: Now, we introduce the finite element space V h .Given an integer k ≥ 1, the finite element space is defined as follows where 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 ∈ V h , on each element T , we define the weak gradient operator And we define the weak divergence s m • ∇ w v m h ∈ P k (T ) associated with s m • ∇v m h on T , as the unique polynomial satisfying Next, we introduce the projection operators, which will be used in the error analysis of the weak Galerkin method.For each element projections on the associated local polynomial spaces respectively.For any mesh element T , we define the global projections Q h and Π h element-wise by The finite dimensional space V h is defined as V h = (V h ) M for the semi-discrete problem (2.8) and an The WG finite element method for the semi-discrete problem (2.8) reads: Find where where, the incoming boundary ∂Ω − h = {e ∈ ∂Ω h : s • n(x) < 0, x ∈ e}.Note that, we have imposed the boundary condition in the weak sense.

Fully-discrete method
The global formulation of the DOWG method (2.16) is expressed as: Find In the simplified version, the DOWG method is presented as follows: where the bilinear form A W G (•, •) and the linear functional F s (•) are given by For any v h ∈ V h , the proposed norm |||•||| is defined as where |||•||| T h is given by

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 (2.17).By using this stability result, we derive the error estimate of the numerical solution.

Stability result
Now, we discuss the coercivity of the bilinear form A W G (•, •) in the following lemma.
Lemma 3.1.For sufficiently small h, we have where α is a positive constant.
Proof.By using the definition of weak divergence (2.15) and integral by parts, it can be easily deduced that Then, we have where I is given by and it satisfies By using condition (2.11), we obtain that Therefore, the bilinear form satisfies This completes the proof.
The direct consequence of the above lemma is the uniqueness result.
Next, we state the trace inequality: For any function φ ∈ H 1 (T ), the following trace inequality holds true (see [39]):

Error estimate
Next, we discuss the error estimate of the DOWG method.Firstly, we establish an error equation associated with the discrete problem (2.17).
Proof.By using the integration by parts and definition (2.15), we deduce that Next, one can obtain that and Finally, by combining (3.4)-(3.6),we get Adding A st (Q h u m , v m h ) both side of the above equation to obtain The above equation into the global formulation is given by By subtracting the discrete problem (2.17) from (3.8), we obtain the required error equation.
In the next lemma, we derive bounds for each term of L(u, v h ).Then, we use the coercivity of the DOWG method to find the estimates of error term e Q h .
Lemma 3.4.Let u m be the solution of the discrete problem (2.8).For v m h ∈ V h , we have the following estimates: Proof.By using L 2 -projection error estimates given in [37, Lemma 3.6], we deduce that Using Cauchy-Schwarz inequality and trace inequality (3.1), we have In a similar way L 3 (u m , v m h ) is bounded by where the assumption (2.11) is used.By the means of trace inequality (3.1), At last, the stabilization term can be estimated by This completes the proof.
Theorem 3.5.Let u m and u m h be the solutions to the problem (2.8) and the discrete WG scheme (2.17), respectively.Then, the error term e Q h ≡ Q h u − u h satisfies: (3.9) Proof.Using the stability result from Lemma 3.1, we have By employing Lemma 3.4, one can obtain that A straightforward computation provides the desired result (3.9).
The projection error is discussed here.By using the trace inequality (3.1) and result form [37, Lemma , it can be easily verified that Next, we combine the above results to deduce the error estimate associated with weak Galerkin approximation in the following theorem.
Theorem 3.6.Let u h = {u m h } be the numerical approximation of u = {u m }.Then, the following estimates Proof.Combing the error estimates (3.9) and (3.12), we obtain that Next, apply the norm definition from (2.18) to derive (3.13).This finalizes the proof.
Finally, we discuss the estimate for the fully discrete error e h := {e m h } M m=0 , where e m h = (u − u m h ).By adding estimates from Theorem 2.1 and Theorem 3.6, the fully discrete error term e h ≡ u − u h satisfies: In the 2D model setting, In the 3D model setting,

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-dimensional models to lower the computational cost.The numerical construction is followed from [30], which is described briefly here.
Firstly, we state that the spatial domain is Ω = (0, 1) 2 and the mesh width in the angular direction (parameterized form) is h θ = π/10 for all the numerical examples.Let T m h be the triangulation of Ω with respect to the angular direction s m .For example, the angular direction s 0 ≡ (cos 0, sin 0) = (1, 0) is associated with the triangulation T 0 h , which has the incoming boundary ∂Ω − s0 = {(0, y)|y = [0, 1]}.In this way, we can fix the spatial triangulation T m h associated with the angular direction s m , then the following algorithm is used for the numerical simulation.

Algorithm 1: Source iteration
Set parameters tol, δ h } on all elements T ∈ T m h by solving the discrete problem (2.17) end We complete the aforementioned process in a single iteration step for all directions, terminate the iteration if a stopping condition (for below test problems, tol = 10 −3 ) is fulfilled, and take u m,i h as u m h .We have employed the software library deal.II [48] for the same computation.All the numerical simulations were executed on an Intel Core i7-10510U CPU@1.80GHzmachine with 8 GB of RAM.For both test problems, the source functions f satisfy the radiative transfer equation (1.1).Furthermore, the absorption and scattering coefficients value is σ t = 2, σ s = 1/2 for all the test problems.
For numerical validation, we discuss error e h in the following norm , and convergence order eoc h is given by In spatial discretization, the coarse mesh (at first level) consists of four squares.Then, each square is equally decomposed into 4 squares, to get the next fine mesh.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 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 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.

Comparison with DG methods
In this part of our article, we discuss a detailed comparison of our proposed DOWG method with the existing DG-type methods such as the discrete-ordinate discontinuous-Galerkin (DODG) method, the discrete- where with F s (v m h ) is given as in (2.16) and Here, c p > 0 is a penalty parameter and ûm h is the numerical trace, which is given by

DODSD method for RTE model (1.1)
DODSD method for the model problem (1.1) is given by where and u m h = u m in on Γ − m .The stabilization parameter δ satisfies δ = ch with some c > 0. Note that, the numerical solutions u h of both DG-type methods satisfy the error estimates (3.15) and (3.16).Now, we discuss the comparison of these methods with the proposed method (2.17) As, we can see from both numerical schemes (4.1) and (4.2), these methods are parameter-dependent.In contrast, the proposed method (2.17) is parameter-free, one of the major advantages of our numerical method.In the numerical section of [29], the authors discussed the impact of parameter c p on the accuracy of the numerical method.
To compare the accuracy of numerical solution with the exact solution of DG methods versus WG methods, we discuss the error estimate |||e h ||| 1 for the Q 1 finite elements.The penalty parameter c p = 0.1 and stabilization parameter δ = h is taken for numerical methods (4.1) and (4.2), respectively.We have presented error estimates in tables 3 and 4 for all three numerical methods, one can see that our proposed WG method performs better than other DG methods.The convergence rate of our method is also slightly better than both DG-type methods.We conclude this subsection with a few remarks that the DOWG method is a viable alternative to DODG/DODSD methods.Compared to DG-type methods, it is easier to use and performs better and no need for penalty/stabilization parameters.However, a complete discussion of this topic is beyond the scope of this paper.

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 the justification of theoretical findings.

s
Email address: maneesh-kumar.singh@imperial.ac.uk (Maneesh Kumar Singh )Preprint submitted to Elsevier February 13, 20241.1.The model problemThis research explores a stationary radiative transfer equation that is modulated by• ∇u(x, s) + σ t (x)u(x, s) − σ s (x)Ku(x, s) = f (x, s), (x, s) ∈ Ω × S, u(x, s) = u in (x, s), (x, s) ∈ Γ − , (1.1)where u is the radiative intensity for all x in the domain Ω, an open bounded polytopal domain inR d (d = 1, 2, 3) with a smooth boundary ∂Ω and for all directions s of the unit sphere S in R d .The incoming and outgoing boundaries are given by

Table 1 :
Discretization error e h and order of convergence eoc h for Example 4.1.

Table 2 :
Discretization error e h and order of convergence eoc h for Example 4.2.