Several applications in the field of physics require an accurate solution of the radiative transfer equation. This equation describes the evolution of the angular flux of particles moving through a material medium. Examples include nuclear engineering [duderstadt1976nuclear, henry1977nuclear], high-energy astrophysics [lowrie1999coupling, mcclarren2008manufactured], supernovae [fryer2006snsph, swesty2009numerical], and fusion [matzen2005pulsed, marinak2001three]. A major challenge when solving the radiative transfer equation numerically is the high-dimensional phase space on which it is defined. There are three spatial dimensions, two directional (angular) parameters, velocity, and time. In many applications, there is additional frequency or energy dependence. Hence, numerical methods to approximate the solution require a carefully chosen phase space discretization.
There are several strategies to discretize the angular variables, and they all have certain strengths and weaknesses [brunner2002forms]. The spherical harmonics (P) method [case1967linear, pomraning1973equations, lewis1984computational]
is a spectral Galerkin discretization of the radiative transfer equation. It uses the spherical harmonics basis functions to represent the solution in terms of angular variables with finitely many expansion coefficients, called moments. The Pmethod preserves rotational symmetry and shows spectral convergence for smooth solutions. However, like most spectral methods, P yields oscillatory solution approximations in non-smooth regimes, which can lead to negative, non-physical angular flux values. Filtering of the expansion coefficients has shown to mitigate oscillations [mcclarren2010robust] and a modified equation analysis has shown that filtering adds an artificial forward-peaked scattering operator to the equation if a certain scaling strength of the filter is chosen.
The discrete ordinates (S) method [lewis1984computational] approximates the radiative transfer equation on a set of discrete angular directions. The S discretization preserves positivity of the angular flux while yielding an efficient and straight forward implementation of time-implicit methods. However, the method is plagued by numerical artifacts, know as ray effects, when there are not enough ordinates to resolve the angular flux. Because increasing the number of ordinates significantly increases numerical costs of simulations, a major task to improve the solution accuracy of S methods is to mitigate these ray effects [lathrop1968ray, morel2003analysis, mathews1999propagation] without simply adding more ordinates.
Various strategies to mitigate ray effects at affordable costs have been developed. In [abu2001angular], a biased quadrature set, which reflects the importance of certain ordinates is used. Furthermore, [lathrop1971remedies] presents a method combining the P with the S method. Further studies for this method can be found in [jung1972discrete, reed1972spherical, miller1977ray], which show a reduction of ray effects. In [morel2003analysis], a comparison of these methods can be found. In [tencer2016ray], computing the angular flux for differently oriented quadrature sets and averaging over different solutions has been proposed to reduce ray artifacts. In [camminady2019ray], a rotated S method has been developed, which rotates the quadrature set after every time iteration. Consequently, particles can move on a heavily increased set of directions of travel, leading to a reduction of ray effects. Analytic results show that rotating the quadrature set plays the role of an angular diffusion operator, which smears out artifacts that stem from the finite number of ordinates. Unfortunately, this method does not allow a straight forward implementation of sweeping, complicating the use of implicit methods.
The idea of this work is to add angular diffusion directly with the help of a forward-peaked artificial scattering operator. We choose this operator so that the effect of artificial scattering vanishes in the limit of infinitely many ordinates, but at finite order adds angular diffusion in such a way that it mitigates ray effects. Unlike the rotated S method in [camminady2019ray], the current approach allows for a straight forward implementation of sweeping, which we use to implement an implicit method.
2 Main idea
In this section, we summarize the relevant mathematical background and introduce notation. We illustrate the problem of ray effects that occurs when discretizing the transport equation in angle and how artificial scattering can be used to mitigate these ray effects. We demonstrate that artificial scattering behaves like a Fokker-Planck operator in the appropriate limit.
2.1 Radiative transfer equation
The radiative transfer equation describes the evolution of the angular flux via
where denotes time, is the spatial variable, and represents the direction. The total cross section is . In the case of scattering, the in-scattering kernel operator describes the gain of particles that were previously traveling along direction and changed to direction . It is given by
is the probability of transitioning from directioninto direction or vice versa. We assume—for simplicity—that the source is isotropic.
2.2 Ray effects
As previously explained, the S method preserves positivity but suffers from ray effects. An example of these artifacts is demonstrated for the line-source benchmark in Fig.1. While the true scalar flux is radially symmetric, the numerical solution has artifacts in the form of oscillations. We will discuss the line-source problem in more detail in Section 4.1.
2.3 Artificial scattering
We propose to address the problem of ray effects by adding artificial scattering to the right-hand side of (1), in the form of an anisotropic scattering operator. The modified system is then
Here, can be any Dirac-like sequence222 While the idea of artificial scattering works with any Dirac-sequence, the asymptotic analysis that is performed later imposes stronger requirements to obtain a Fokker-Planck operator in the respective limit.
While the idea of artificial scattering works with any Dirac-sequence, the asymptotic analysis that is performed later imposes stronger requirements to obtain a Fokker-Planck operator in the respective limit., i.e.
for any sufficiently smooth function as . In our experiments, we choose
where the error function satisfies as . The proposed method, which we call artificial scattering-S (as-S), has the following effects:
Similar to the artificial viscosity used to stabilize spatial discretizations of hyperbolic operators [leveque1992numerical, Chapter 16.1]
, the artificial scattering adds an angular diffusion term to the radiative transfer equations. This term should vanish when the discretization is refined. Therefore, the variance of the artificial scattering kernel should be chosen to vanish in the limit. We choose the variance to be , where is a constant, user-determined parameter. This choice ensures that scales the average width of quadrature points, meaning that the domain of influence includes roughly the same number of ordinates when increases. In the limit, the as-S solution converges to the classical S solution.
The total number of particles is preserved by the artificial scattering term. Higher order moments are however dampened by the magnitude of artificial scattering. A beam of particles inside a void will be subject to scattering by the as-S method, however artifacts that result from the standard S method dominate the overall error, unless the beam is aligned with the quadrature set.
as-S has similarities to the P-equivalent S method [miller1977ray]: To mitigate ray effect, this method adds a fictitious source to the radiative transfer equation. This source, though derived by a different strategy, requires similar modifications of the standard S implementation. The main difference is that the artificial scattering kernel of as-S is forward peaked, which can be used to design an efficient numerical treatment.
as-S can be compared to filtered P [mcclarren2010robust, hauck2019filtered], since the artificial scattering acts as a filter on the moment level.
The as-S equation (3) can—with appropriate boundary and initial conditions—be solved in a straight-forward manner using common S implementations. When discussing one such implementation, we will focus on implicit discretization techniques and derive an efficient algorithm to treat the artificial scattering term.
2.4 Artificial scattering kernel
To better distinguish between the two types of scattering, we will call the naturally occurring scattering of (1) physical scattering and the scattering in (4) artificial scattering. The way the artificial scattering is written in (4), it includes in-scattering and out-scattering. We can split this further into
2.5 Modified equation analysis
According to Pomraning [pomraning1992fokker], the Fokker-Planck operator can be a legitimate description of highly peaked scattering. This is true if (i) the scattering kernel is a Dirac sequence, and (ii) the transport coefficients are of order . The resulting modified equation then reads
where is the Laplace operator in spherical coordinates. We have already shown (i). To verify (ii), let . Then
where and denote the gamma function and the upper incomplete gamma function, respectively. Since (ii) implies that , the operator vanishes if we let .
We set in the discrete case so that the angular diffusion vanishes if the number of ordinates tends to infinity. This analysis shows that the product controls the strength of the added angular diffusion. Section 4.1 confirms this numerically.
In the following section, we will discuss discretization and implementation of the presented as-S method, laying the focus on how to incorporate artificial scattering into existing S codes.
3.1 S discretization
For sake of completeness, we briefly summarize the S method. Given a finite number of ordinates and defining the S solution the S method solves the semi-discretized system of equations
Here, are quadrature weights, chosen such that
To compute numerical solutions, we still need to discretize (15) in space and time. In this work, we will investigate solutions for both implicit and explicit time discretizations. The explicit code uses Heun’s method as well as a minmod slope limiter. It is based on [camminady2019ray, garrett2013comparison], which provide a detailed description of the chosen methods. The implicit discretization and an efficient strategy to integrate artificial scattering in a given implicit code framework will be discussed in Sections 3.4 and 3.5.
3.2 Adding artificial scattering to the S equations
Here, and are normalization factors. While on the continuous level, these factors are the same for every direction, we obtain a dependency on the chosen ordinate due to the non-uniform discretization in angle. These normalization factors are needed to obtain a simple expression for the out-scattering terms. Moving these terms to the left-hand side of (17) stabilizes the source iteration used in the implicit method.
It remains to pick an adequate quadrature set. When applying artificial scattering, the solution smears out along the imprint of the discrete set of ordinates. To ensure an evenly spread artificial scattering effect, a quadrature with a highly uniform ordinate distribution should be chosen. Commonly, the construction of a quadrature set starts at a chosen planar geometry, which is discretized and then mapped onto the surface of a sphere. The mapped nodes of the previously chosen discretization are then taken to be the quadrature points while the weights are determined by the area associated with these points. An even node distribution is achieved by using an icosahedron as initial planar geometry, which one can find in Figure 2. Each face of the icosahedron is triangulated to generate the nodes which will be mapped onto the surface of the sphere. There are different strategies to perform this triangulation and we choose an equidistant spacing of points on each line of the triangle as well as the corresponding points inside the triangle. The corresponding weight is the area of the hexagon which lies around the given node and is defined by connecting the midpoints of the neighboring triangles. For more details on the icosahedron quadrature, see [icosahedron]. From now on, we will exclusively use this quadrature in all S and as-S computations.
3.4 Implicit time discretization
Implicit time discretization methods provide stability for large time steps, which are crucial in applications involving different time scales. However, when discretizing the radiative transfer equation, they require a matrix solve in every time step, which is commonly performed by a Krylov solver [faber1988look, ashby1991preconditioned]. We start with an implicit Euler discretization, where we, in an abuse of notation, denote the flux at the new time step by and at the old time step by . The equivalent as-S system is
Defining the streaming operator as well as the modified source , we can put this into more compact notation
First, let us numerically treat the artificial scattering in the same way as commonly done for physical scattering. The physical in-scattering kernel can be written as
where carries the respective expansion coefficients of the scattering kernel, maps from the ordinates to the moments and from the moments back to the angular space. Making use of this strategy to represent the artificial scattering, we get
When denoting the moments as , equation (19) becomes
Inverting and applying to both sides yields the fixed point equations
Note that with , this is the standard equation to which a Krylov solver is applied. Choosing a non-zero artificial scattering strength can result in significantly increased numerical costs when solving (22) with a Krylov method: To show this, let us move to the discrete level, i.e. discretizing the directional domain, which requires picking a finite number of moments. In this case becomes a diagonal matrix with entries falling rapidly to zero (in the case of isotropic scattering, only the first entry is non-zero). Hence, few moments are required to capture the effects of physical scattering. However, since the artificial scattering kernel is strongly forward-peaked, the entries of the diagonal matrix do not fall to zero quickly, meaning that the method requires a large number of moments to include artificial scattering, which results in a heavily increased run time [hauck2019filtered]. The slow decay of the Legendre moments for is visualized in Fig. 3.
In order to be able to choose the reduced number of moments required to resolve physical scattering, we move the artificial scattering into the sweeping step. Hence, going back to equation (19), we only perform the moment decomposition on the physical scattering to obtain
Moving the operator to the right hand side and taking moments yields
The Krylov solver is then applied to this fixed-point iteration. In contrast to (22), the physical scattering term dictates the number of moments. The computation of is performed by a source iteration, where the general equation is solved by iterating on
This iteration is expected to converge fast since effects of artificial scattering will be small in comparison to physical scattering.
3.5 Implementation details
At this point, we choose a finite number of ordinates and moments, i.e. the flux
is now a vector with dimensionand the moments have finite dimension . Consequently, operators applied to the directional space become matrices. For better readability, we abuse notation and reuse the same symbols as before.
We observed that a second-order spatial scheme is required to capture the behavior of the test cases used in this work. To ensure an efficient sweeping step, we use a second-order upwind stencil without a limiter. Let us denote the operator discretized in space and direction by . For ease of presentation, we assume a slab geometry. i.e. we have the spatial variable and the directional variable . In the following, we split the directional variable into and . An extension to arbitrary dimension is straight forward. Now with and , we can write the discretized streaming operator as
The numerical flux for is then given by
and for by
This scheme is L stable, which we show in Appendix B. Let us now discuss the implementation of the implicit method in more detail. As mentioned earlier, a source iteration (25) is required to invert the operator . For an initial guess and an arbitrary right hand side , this iteration is given by Alg. 1. Note that the discrete artificial in-scattering is a sparse matrix, which guarantees an efficient evaluation of the matrix vector product in (25). Furthermore, the inverse of can be computed by a sweeping procedure.
In order to get a good error estimator, we set the constantin Algorithm 1, where is an estimate of the Lipschitz constant and is a user-determined parameter. Our implementation solves the linear system of equations
using a GMRES solver. The solver requires the evaluation of the left-hand side for a given with an initial guess , which is given by Alg. 2.
The main time stepping scheme is then given by Alg. 3. After initializing and , the right hand side to (29a) is set up in line 4. Line 5 then solves the linear system (29a) and Line 6 determines the time-updated flux from the moments .
There exist several ways to modify the presented algorithm to achieve higher performance. For example, one can modify the presented method by not fully converging the source iteration in Alg. 1. Instead, only a single iteration can be performed to drive the moments and the respective angular flux to their corresponding fixed points simultaneously. In numerical tests, we observe that this will significantly speed up the calculation. However, since we do not focus on runtime optimization, we do not further discuss this idea and leave it to future work.
In the following, we evaluate the proposed method within the scope of two numerical test cases: (i) the line-source problem is used as it is inherently prone to ray-effects when using the S method, and (ii) the lattice test case models—in a very simplified way—neutrons in a fission reactor with a source and heterogeneous materials. For both problems, we present results for the explicit and implicit methods, respectively.
Both test cases are computed on a two-dimensional regular grid for the spatial variable. We project the for onto the --plane.
The code used to compute the numerical results is published under the MIT license in a public repository at https://github.com/camminady/SN.
4.1 Line-source test case
The goal of this test case is to numerically compute the Green’s function for an initial isotropic Dirac-mass at the origin, i.e. , which is realized as a narrow Gaussian in space with and . We choose . The spatial discretization varies from for a coarse grid to points on the domain for a fine grid. There exists a semi-analytical solution to the full transport equation for this problem due to Ganapol et al. [ganapol2001homogeneous]. The exact solution consists of a circular front moving away from the origin as well as a tail of particles which have been scattered or not emitted perpendicularly from the center. We chose the line-source because it is a test case that lays bare almost any artifact an angular discretization might suffer from.
The parameters of the artificial scattering have been set to with as the number of quadrature points. Obviously, choosing these parameters requires some experience. However, as in the case of filtering for [mcclarren2010robust], both parameters can be adjusted for coarse angular and spatial grids, and are expected to be valid for finer grids, which seems to be the case for the line-source problem as well. The CFL number, i.e. the ratio of the time step and the spatial cell size is for the explicit calculations and for the implicit calculation. For the implicit discretization, the tolerance for the GMRES solver was set to , and we considered the inner source iteration to be converged at an estimated error of .
An overview of the analytics that we perform is given in Fig. 4. We evaluate the scalar flux at the final time step. We have performed an explicit S computation with and . In both rows of Fig. 4, the left column shows the scalar flux at the final time step. The first row shows the solution along cuts through the domain on the right with the respective cuts on the left in white. The second row shows the solution along circles with different radii on the right and the respective circles on the left. Strong oscillations are visible due to ray-effects. For the first row, the analytical solution is given in green in the right column image. In the lower row, the analytical solution is constant along a circle with a certain radius, visualized for , , and in green.
In Fig. 5 and Fig. 6 we see the same summary of results, now for an explicit and implicit computation, respectively. In both computations, ray-effects have been reduced significantly when comparing the results with Fig. 4 despite the same number of quadrature points. The implicit calculation looks slightly more diffusive. However, the line-source problem is not a problem that would be computed implicitly in the first place and we only use it to illustrate the expected behavior for implicit computations.
The values for and in the fine calculations are determined from a parameter study using coarse spatial and angular grids. The results of this parameter study are given in Fig. 7 for the explicit algorithm and in Fig. 8 for the implicit algorithm.
A single simulation for the coarse configuration takes times the time of a single computation for the fine configuration. Consequently, the full parameter study with all configurations can be performed for less than the costs of a single fine computation. For the optimal parameter configuration, the error decreases down to for the explicit case and down to for the implicit case.
In both cases, implicit and explicit, we observe a region of parameters that yield similarly good results. This behavior mostly matches the predicted relation from the asymptotic analysis, i.e. when is small, controls the effect of artificial scattering.
We also investigate the performance of the as-S method when measured in runtime and in memory consumption. Consider therefore the results presented in Fig. 9 and Fig. 10. Both figures summarize the results for the line-source test case computed with the S and as-S methods for different values of . Fig. 9 measures the error between the numerical solution and the analytical solution in the L norm, called . Fig. 10 considers the semi-norm.
We observe an increase in runtime when activating artificial scattering, but a decrease in the errors and . On the right, the errors are plotted against the number of ordinates which ultimately dictates the memory consumption. For example, an S takes about as long as an as-S computation and yields a similar error. However, the number of ordinates can be reduced from 492 to 162. For both, and , the effect of artificial scattering vanishes in the limit of .
4.2 Lattice test case
We also investigate the lattice test case [brunner2002forms, brunner2005two], depicted in Fig. 11. A constant, isotropic source is placed in the center of the domain in the orange square. In the white cells, the material is purely scattering, whereas the orange and black squares are purely absorbing. The boundary conditions are vacuum. All test case parameters are listed in Table 1.
In Fig. 12 we see the as-S solution to the lattice problem on the left, the S solution in the center, and the S solution on the right. Here, S uses ordinates while S and as-S use ordinates.
We take the S solution with as our reference solution. When comparing the as-S solution with the S solution, we see an improvement in the solution quality. Ray-effects are better mitigated in regions where the scalar flux is small. The number of ordinates is kept constant.
Additionally, Fig. 13 puts the as-S solution and the S solution side-by-side in the center frame. Minor ray-effects are visible when looking at the white isoline. However, the number of ordinates has been reduced by a factor of .
Similar to the line-source test case, we set and for explicit calculation, and and for the implicit calculation.
We also perform simulations for the lattice problem with the implicit time discretization. However, since the chosen scheme is only L-stable, the solution becomes negative for the lattice test case as illustrated in Fig. 14. Nevertheless, Fig. 15 demonstrates the inherent advantage when performing implicit computations: We are able to use a very large CFL number, thus reducing the number of time steps and the overall computational costs drastically. Note that the scheme preserves positivity for the chosen CFL numbers.
5 Conclusion & Outlook
We have presented a new ray effect mitigation technique that relies on an additional, artificial scattering operator introduced into the radiative transfer equation. When the number of ordinates tends to infinity, the artificial scattering vanishes and the modified equation reduces to the original transport equation. In this case, when choosing the product of the scattering strength and the variance constant, the term tending to zero with the slowest rate is the Fokker-Planck operator.
The artificial scattering operator can be integrated into standard S codes. Solution algorithms, both for the explicit and implicit case, have been presented and rigorously analyzed in the non-standard, implicit case. To avoid using a large number of moments for the Krylov solver in the implicit case, we propose to invert the artificial scattering operator by a source iteration.
We have presented numerical results for the line-source and lattice test case. The results demonstrate that artificial scattering yields the same accuracy as S , but for a reduced number of ordinates.
For the second-order implicit computations the solutions might turn negative since L stability does not guarantee positivity of the solution. However, when choosing a sufficiently large CFL number, the solution values in our numerical experiments remain positive. A rigorous investigation of this effect and possibly the derivation of a CFL number ensuring positivity is left to future work.
Note that our test cases chose a constant value for the artificial scattering strength, however it seems plausible to make this strength spatially dependent to ensure that artificial scattering is only turned on when required. It remains to demonstrate the feasibility of the as-S method in real-world applications using large-scale, highly parallelizable codes.
The authors wish to thank Ryan G. McClarren (University of Notre Dame) for many fruitful discussions.
Appendix A Icosahedron quadrature
The quadrature points and weights for the quadrature in Section 3.3 are given for order 2 (12 quadrature points). Every line contains four entries: the , , and position, as well as the quadrature weight. The quadrature weights sum to . All entries are in double precision. The quadratures for oder 2, order 3 (42 quadrature points), order 4 (92 quadrature points), and order 5 (162 quadrature points) can be downloaded as .txt files from a public repository at github.com/camminady/IcosahedronQuadrature.
Appendix B Implicit second order upwind scheme
In the following, we show that the chosen numerical flux is stable. For simplicity, we look at the one-dimensional advection equation
with . A finite volume discretization is given by
where we use . A second order, implicit numerical flux is given by
Let us check if the scheme dissipates the entropy . For this we multiply our scheme (31) with , i.e. we obtain
Now one needs to remove the cross term which can be done by reversing the binomial formula
Plugging this formulation for the cross term into (34) and making use of the definition of the square entropy gives
This shows that in order to achieve entropy dissipation, i.e.
Note that the first term of , which essentially comes from the implicit time discretization is always positive. It remains to show that is positive as well. Let us rewrite this term for all spatial cells as a matrix vector product. I.e. when collecting the solution at time step for all spatial cells in a vector , this term becomes
where is a lower triangular matrix. This product can be symmetrized with , meaning that we have . For our stencil, the matrix has entries on the diagonal and as on the lower and upper diagonals. Positivity of and thereby of the entropy dissipation term in (B) is guaranteed if
is positive definite, i.e. has positive eigenvalues. The eigenvalues forhave been computed numerically to verify positivity in Fig. 16.