We consider the numerical solution of the radiative transfer equation in plane parallel geometry
where and , and denotes the thickness of the slab and
is the cosine of the polar angle of a unit vector. The functionmodels the equilibrium distribution of some quantity, like neutrons or photons Chandrasekhar60 ; CaseZweifel67 . The basic physical principles embodied in (1) are transport, which is modeled by the differential operator , attenuation with rate and scattering with . Internal sources are described by the function . In this work we will close the radiative transfer equation using inflow boundary conditions
Such transport problems arise when the full three-dimensional model posed on obeys certain symmetries CaseZweifel67 . It has been studied in many instance due to simpler structure compared to three dimensional problems without symmetries; while the methodology presented here directly carries over to the general case, see also the numerical examples presented below.
Classical deterministic discretization strategies are based on a semidiscretization in . One class of such strategies are the -approximations, which are spectral methods based on truncated spherical harmonics expansions Vladimirov61 ; LewisMiller84 , and we refer to Pitkaranta1977 ; ManResSta00 ; EggerSchlottbom12 for variational discretization strategies using this approximation. The major advantage of -approximations is that the scattering operator becomes diagonal. In addition, the matrix representation of the transport operator is sparse. The main drawbacks of the -method are that the variational incorporation of the inflow boundary condition introduces a dense coupling of the spherical harmonics expansion coefficients making standard -approximations quite expensive to solve. We note, however, that a modified variational formulation of the -equations has recently been derived that leads to sparse matrices EggerSchlottbom18 . In any case, the success of spectral approximation techniques depends on the smoothness of the solution. In general, the solution is not smooth for , which is related to the inflow boundary conditions (2). Hence, the -approximations will in general not converge spectrally.
A second class of semidiscretizations are discrete ordinates methods that use a quadrature rule for the discretization of the -variable CaseZweifel67 , with analysis provided in JohnsonPitkaranta83 ; Asadzadeh86 . Such methods are closely related to discontinuous Galerkin (DG) methods, see, e.g., ReedHill73 ; WareingMcGheeMorelPautz2001 ; RagusaGuermondKanschat12 ; KanschatRagusa14 ; KophaziLathouwers2015 . While allowing for local angular resolution, the main obstruction in the use of these methods is that the scattering operator leads to dense matrices, and a direct inversion of the resulting system is not possible in realistic applications. To overcome this issue, iterative solution techniques have been proposed. An often used iterative technique is Richardson iteration, i.e., the source iteration MarchukLebedev86 ; AdamsLarsen02 , but other Krylov space methods exist WarsaWareingMorel2004 . The key idea in the source iterations is to decouple scattering and transport, and to exploit that the transport part can be inverted efficiently. If , the convergence of these iterative methods is slow, and several preconditioning techniques have been proposed KanschatRagusa14 ; MarchukLebedev86 ; AdamsLarsen02 ; WarsaWareingMorel2004
. Among the most used and simple preconditioners is the diffusion synthetic acceleration method (DSA), in which a diffusion problem is solved in every iteration. This is well motivated by asymptotic analysisLarsenKeller ; EggerSchlottbom2014 . While Fourier analysis can be applied to special situations MarchukLebedev86 ; AdamsLarsen02 , the convergence analysis is mainly open for the general case, i.e., for jumping coefficients or non-periodic boundary conditions. A further complication in using the DSA preconditioner is that the resulting iterative scheme might diverge if the discretization of the diffusion problem and the discrete ordinates system is not consistent AdamsLarsen02 .
The contribution of this paper is to develop discretizations that allow for local resolution of the non-smoothness of the solution, and which lead to discrete problems that can be solved efficiently by diffusion synthetic accelerated source iterations. Our approach builds upon an even-parity formulation of the radiative transfer equations derived from the mixed variational framework given in EggerSchlottbom12 , where -approximations have been treated in detail. Using the framework of -methods MarchukLebedev86 , we show that an infinite dimensional DSA preconditioned source iteration converges already for the continuous problem. We present conforming -type approximations spaces for the semidiscretization in , and prove quasi-best approximation properties. In order to solve resulting linear systems, we employ a DSA preconditioned Richardson iteration, which is just the infinite dimensional iteration projected to the approximation spaces. In particular the finite dimensional iteration is guaranteed to converge for any discretization. Moreover, the inversion of the transport problem can be parallelized straightforward. In numerical experiments, even when employing low-order approximations, we observe that the developed method does not suffer from the ray effect, which is typically observed for discrete ordinates methods LewisMiller84 ; Brunner05 .
The outline of the paper is as follows: In Section 2 we introduce basic notation and recall the relevant function spaces. In Section 3 we introduce the even-parity equation for (1), show its well-posedness, and formulate an infinite dimensional DSA preconditioned source iteration, for which we show convergence. The approximation spaces are described in Section 4 and well-posedness of the Galerkin problems as well as quasi-best approximation results are presented. In Section 5 we discuss the efficient iterative solution of the resulting linear systems and provide a convergence proof for the discrete DSA scheme. Section 6 presents supporting numerical examples for slab geometry and for multi-dimensional problems that show the good approximation properties of the proposed methods as well as fast convergence of the iterative solver in multi-dimensional problems and in the diffusion limit. The paper ends with some conclusions in Section 7.
2 Function spaces and further preliminaries
Following Agoshkov98 we denote by with the usual Hilbert space of square integrable functions with inner product
and induced norm . Furthermore, we define the Hilbert space
of functions with square integrable weak derivatives with respect to the weighted differential operator endowed with the corresponding graph norm.
In order to deal with boundary data, let us introduce the Hilbert space that consists of measurable functions for which
is finite, and we denote by the corresponding inner product on . Similarly, denotes the space of outflow data. We have the following trace lemma Agoshkov98 .
If , then there exist traces and and
with a constant independent of and .
As a consequence of the trace lemma and the density of smooth functions in the following integration-by-parts formula is true
Throughout the manuscript we make the following basic assumption:
are non-negative and .
Let denote the -projection onto constants in , i.e.,
Since is strictly positive, we can define the norms on as follows
2.1 Even-odd splitting
and oddparts of a function are defined as
Even-odd decompositions are frequently used in transport theory Pitkaranta1977 ; Vladimirov61 . Since, as functions of , even and odd function are orthogonal in , we can decompose into orthogonal subspaces containing even and odd functions, respectively,
Similarly, we will write . As in EggerSchlottbom12 , we observe that for any , and for . It turns out that the natural space for our formulation is
3 Weak formulation of the slab problem
We follow the steps presented in EggerSchlottbom12 for multi-dimensional problems. The key idea is to rewrite the slab problem into a weak formulation for the even and odd parts of the solution. Multiplication of (1) with a test function and using orthogonality of even and odd functions gives that
Integration-by-parts (3) applied to the first term on the left-hand side yields that
Due to symmetries, we have that . Using (2), we have that on the inflow boundary, which leads to
Thus, for any , it holds that
Testing (1) with an odd test function , we obtain that
Let and . Find such that
where the bilinear form is given by
and the linear form is defined as
We endow with the norm induced by the bilinear form defined (7), i.e.,
Using the Cauchy-Schwarz inequality, we obtain the following result.
The linear form defined in (8) is bounded, i.e., for all it holds
Since the space endowed with the inner product induced by is a Hilbert space, the unique solvability of Problem 3.1 is a direct consequence of the Riesz representation theorem.
Let Assumption (A1) hold true. Then Problem 3.1 has a unique solution . Moreover, we have the bound
Setting and , one can show that , i.e., satisfies (1) in , cf. EggerSchlottbom12 . Using the trace lemma 2.1 and partial-integration (3), we further can show that satisfies the boundary conditions (2) in the sense of traces. Hence, Theorem 3.3, independently, leads to a well-posedness result as Lemma 2.2 for (1)-(2).
3.3 DSA preconditioned source iteration in second order form
As a preparation for the numerical solution of the discrete systems that will be described below, let us discuss an iterative scheme in infinite dimensions for solving the radiative transfer equation. The basic idea is a standard one and consists of decoupling scattering and transport in order to compute successive approximation, viz., the source iteration MarchukLebedev86 ; AdamsLarsen02 . Next to the basic iteration, we describe a preconditioner which resembles diffusion synthetic acceleration (DSA) schemes using the notation of AdamsLarsen02 or a scheme using the terminology of MarchukLebedev86 .
In order to formulate the method, we introduce the following bilinear forms
and denote the induced semi-norm and norm by and , respectively.
The iteration scheme is defined as follows: For given, compute as the unique solution to
The half-step error satisfies
The key idea is then to construct an easy-to-compute approximation to by Galerkin projection onto a suitable subspace . This approximation is then used to correct to obtain a more accurate approximation to . Define the following closed subspace of
i.e., consists of functions in that do not depend on . The correction is then computed by Galerkin projection of (11) to :
and the new iterate is defined as
If is a good approximation to , then is small. The convergence proof of the iteration , is based on the equivalence of the norms induced by and .
There exists a constant such that for all
From the definition of the bilinear forms and and positivity of , we deduce that
where we also used that . ∎
The proof follows (MarchukLebedev86, , X.9). Let us define the residual as the unique solution to
This implies for all , and hence
In addition, using Lemma 3.5, we have that
i.e., . Furthermore, by definition of and ,
for all , i.e, . Using these relations, we obtain that
Another application of Lemma 3.5 yields that
which proves the main assertion. The proof is concluded, by observing that, if and are constant, the constant in the latter inequality can be bounded by
Since is the orthogonal projection of to in the -inner product it holds that
The assertion then follows from Lemma 3.6. ∎
The convergence analysis presented in this section carries over verbatim to multi-dimensional problems without symmetries, and it can be extended immediately to more general (symmetric and positive) scattering operators.
Problem (13) is the weak formulation of the diffusion equation
with , complemented by Robin boundary conditions, which shows the close relationship to DSA schemes.
The convergence analysis for the iteration without preconditioning, can alternatively be based on the following estimates. These estimates are the only ones in this paper that exploit that the scattering operator is related to the
The convergence analysis for the iteration without preconditioning, can alternatively be based on the following estimates. These estimates are the only ones in this paper that exploit that the scattering operator is related to the-projector . First note that
and that . Setting , an the Cauchy-Schwarz inequality yields
for any . Choosing , shows that parts of the error are smoothened independently of , i.e.,
while the angular average is hardly damped. In any case, this shows that converges to zero. It remains open how to exploit such an estimate to improve the analysis of the DSA preconditioned scheme above as it seems difficult to relate the latter smoothing property to the best approximation error of in the -norm.
4 Galerkin approximations
In this section, we construct conforming approximation spaces in a two-step procedure. In a first step, we discretize the -variable using discontinuous ansatz functions. In a second step, we discretize the -variable by continuous finite elements. Before stating the particular approximation space, we provide some general results. Let us begin with the definition of the discrete problem.
Let , , and let and be defined as in Problem 3.1. Find such that for all there holds
Since the bilinear form induces the energy norm that we use in our analysis, we immediately obtain the following best approximation result.
Let be a closed subspace of . Then, there exists a unique solution of Problem 4.1 that satisfies the a-priori estimate
and the following best approximation error estimate
In the next sections, we discuss some particular discretizations. We note that these generalize the spherical harmonics approach presented in EggerSchlottbom12 .
4.1 semidiscretization in
Since we consider even functions, we require that the partition of the interval for the variable respects the point symmetry . For simplicity, we thus partition the interval only, and the partition of is implicitly defined by reflection.
Let be a positive integer, and define intervals , , such that and , and set and . Denote
the characteristic function of the interval. For , we define the piecewise functions
where denotes the th Legendre polynomial. Hence, is an -orthonormal basis for the space of polynomials of degree on each interval . For , we set , and for , we set . The semidiscretization of the even parts is then
If we partition the interval for the angular variable by a single element, we obtain truncated spherical harmonics approximations, see, e.g., LewisMiller84 ; EggerSchlottbom12 . Partitioning of into two symmetric intervals corresponds to the double method LewisMiller84 , which generalizes in multiple dimensions to half space moment methods
, which generalizes in multiple dimensions to half space moment methodsDubroca03 . The latter can resolve the non-smoothness of at , and, thus might yield spectral convergence on the intervals and .
4.2 Fully discrete scheme
In order to obtain a conforming discretization, we approximate the coefficient functions using -conforming elements. Let , and such that
be a partition of . Let and denote the space of polynomials of degree . The full approximation space is then defined by
We note that this space satisfies the compatibility condition , which makes this pair of approximation spaces suitable for a direct approximation of a corresponding mixed formulation, cf. EggerSchlottbom12 . The even-parity formulation that we consider here corresponds then to the Schur complement of the mixed problem, cf. EggerSchlottbom12 . The reader should note the different degrees in the polynomial approximations, e.g., if is piecewise constant in angle, then is piecewise linear.
5 Discrete preconditioned source iteration
In order to solve the discrete variational problem defined in Problem 4.1, we proceed as in Section 3.3 but with and replaced by and , respectively. We note that consists of functions in that do not dependent on .
The finite dimensional counterpart of the DSA preconditioned source iteration is then defined as follows: For given , compute as the unique solution to
The correction is defined by Galerkin projection of to :
and the new iterate is defined as
Using the same arguments as above, we obtain the following convergence result.
Similar to Remark 3.9, is the Galerkin projection to of the solution to
with -grid dependent diffusion coefficient and . If piecewise constant functions in angle are employed, then .
Once the scattering term in the right-hand side of (19) has been computed, the half-step iterate can be computed independently for each direction, and, thus, its computation can be parallelized.
6 Numerical examples
In this section, we report on the accuracy of the proposed discretization scheme and its efficient numerical solution using the DSA preconditioned source iteration of Section 5. We restrict our discussion to a low-order method that offers local resolution. To do so, we fix and in the definition of , while might be large. Hence, the approximation space consists of discontinuous, piecewise constant functions in the angular variable and continuous, piecewise linear functions in .
6.1 Manufactured solutions
To investigate the convergence behavior, we use manufactured solutions, i.e., the exact solution is
with parameters , and , and source terms defined accordingly. We computed the numerical solution using the DSA preconditioned iteration (19), (20), (21). We stopped the iteration using the a-posteriori stopping rule
where is a chosen tolerance. The approximations of the even and odd parts of the solutions are recovered as described in Remark 4.4.
For the chosen parameters, we have that , which leads to a predicted convergence rate of . Table 1 shows the errors for different discretization parameters and . As expected we observe a linear rate of convergence with respect to , cf. Theorem 4.2. In addition, we observe that the preconditioned source iteration converged within at most iterations, and the expression