Abstract
We analyze the properties and compare the performance of several positivity limiters for spectral approximations with respect to the angular variable of linear transport equations. It is well-known that spectral methods suffer from the occurrence of (unphysical) negative spatial particle concentrations due to the fact that the underlying polynomial approximations are not always positive at the kinetic level. Positivity limiters address this defect by enforcing positivity of the polynomial approximation on a finite set of preselected points. With a proper PDE solver, they ensure positivity of the particle concentration at each step in a time integration scheme. We review several known positivity limiters proposed in other contexts and also introduce a modification for one of them. We give error estimates for the consistency of the positive approximations produced by these limiters and compare the theoretical estimates to numerical results. We then solve two benchmark problems with these limiters, make qualitative and quantitative observations about the solutions, and then compare the efficiency of the different limiters.





Similar content being viewed by others
Notes
In this paper, the term “concentration” refers the integral of the kinetic distribution over the momentum/angular space. The concentration is a function of position and time only.
In general, \(n = (N+1)^2\); however, in reduced geometries, the components of \(\mathbf {m}\) are not all linearly independent. In such cases, n will be smaller.
Note that \(\mathbf {m}\) depends on N and could be denoted as \(\mathbf {m}_N\). We suppress this notation here for the sake of simplicity.
Note that the results presented in this section only focus on the consistency properties of the limiters. A full convergence analysis for the P\(_N\) and FP\(_N\) equations with limiters is in the scope of future work.
Here we define the Sobolev spaces on \([-1,1]\) in (19) using weak derivatives and space interpolations as in [8, 41]. It is known that the Sobolev–Slobodeckij spaces used in [5] are equivalent to (19). On the other hand, we define the Sobolev spaces on \(\mathbb {S}^2\) in (20) via expansion coefficients as in [12, 20]. In [20, Section 8.1], it is shown that (20) is equivalent to a norm based on weak derivatives and space interpolations on \(\mathbb {S}^2\). We choose to use the form in (20) for simplicity.
This is primarily due to the fact that, on \(\mathbb {S}^2\), more points are required for quadratures to achieve precision \(2N+1\). This means that there is generally no polynomial interpolant of \(\varphi \) in \(\mathbb {P}_N\), thus the second equality in (31) does not hold on \(\mathbb {S}^2\).
For general problems, it may not be possible to take advantage of symmetries.
Note that the negative particle concentrations are colored in white.
Similar results are obtained for the \(L^1\) and \(L^\infty \) spatial errors.
The time step \(\varDelta t\) is also refined in such a way that the ratio \(\varDelta t/h\) stays fixed.
References
Alldredge, G.W., Hauck, C.D., Tits, A.L.: High-order entropy-based closures for linear transport in slab geometry II: a computational study of the optimization problem. SIAM J. Sci. Comput. 34(4), B361–B391 (2012)
Atkinson, K.: Numerical integration on the sphere. J. Aust. Math. Soc. Ser. B 23, 332–347 (1982)
Atkinson, K., Han, W.: Spherical Harmonics and Approximations on the Unit Sphere: An Introduction. Springer, Berlin (2012)
Bergh, J., Löfström, J.: Interpolation Spaces: An Introduction. Grundlehren der mathematischen Wissenschaften. Springer, Berlin (1976)
Bernardi, C., Maday, Y.: Polynomial interpolation results in Sobolev spaces. J. Comput. Appl. Math. 43(1), 53–80 (1992). https://doi.org/10.1016/0377-0427(92)90259-Z
Brunner, T.A.: Forms of approximate radiation transport. Technical Report SAND2002-1778, Sandia National Laboratories (2002)
Brunner, T.A., Holloway, J.P.: Two-dimensional time-dependent Riemann solvers for neutron transport. J. Comput. Phys. 210, 386–399 (2005)
Canuto, C., Quarteroni, A.: Approximation results for orthogonal polynomials in Sobolev spaces. Math. Comput. 38(157), 67–86 (1982)
Case, K., Zweifel, P.: Linear Transport Theory. Addison-Wesley, Reading, MA (1967)
Cercignani, C.: The Boltzmann Equation and its Applications, Applied Mathematical Sciences, vol. 67. Springer, New York (1988)
Cercignani, C., Illner, R., Pulvirenti, M.: The Mathematical Theory of Dilute Gases, Applied Mathematical Sciences, vol. 106. Springer, New York (1994)
Dai, F., Xu, Y.: Approximation Theory and Harmonic Analysis on Spheres and Balls. Springer Monographs in Mathematics. Springer, New York (2013). https://books.google.com/books?id=8SA_AAAAQBAJ. Accesses 2013
Dautray, R., Lions, J.L.: Mathematical Analysis and Numerical Methods for Science and Technology, Volume 6: Evolution Problems II. Spinger, Berlin (2000)
Deshpande, S.M.: Kinetic theory based new upwind methods for inviscid compressible flows. In: American Institute of Aeronautics and Astronautics, New York (1986). Paper 86-0275
Frank, M., Hauck, C., Küpper, K.: Convergence of filtered spherical harmonic equations for radiation transport. Commun. Math. Sci. 14(5), 1443–1465 (2016)
Ganapol, B.D.: Homogeneous infinite media time-dependent analytic benchmarks for X-TM transport methods development. Technical report, Los Alamos National Laboratory (1999)
Garrett, C.K., Hauck, C.D.: A comparison of moment closures for linear kinetic transport equations: the line source benchmark. Transp. Theory Stat. Phys. 42, 203–235 (2013)
Garrett, P.: Harmonic analysis on spheres, II (2011). http://www.math.umn.edu/~garrett/m/mfms/notes_c/spheres_II.pdf. Accesses 2011
Gottlieb, D., Gottlieb, S., Hesthaven, J.: Spectral Methods for Time-Dependent Problems. Cambridge University Press, New York (2007)
Guo, B.: Spectral Methods and Their Applications. World Scientific, Singapore (1998)
Harten, A., Lax, P.D., Leer, V.: On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev. 25, 35–61 (1983)
Hauck, C.D., McClarren, R.G.: Positive \({P_N}\) closures. SIAM J. Sci. Comput. 32(5), 2603–2626 (2010)
Laboure, V.M., McClarren, R.G., Hauck, C.D.: Implicit filtered P\(_N\) for high-energy density thermal radiation transport using discontinuous Galerkin finite elements. J. Comput. Phys. 321, 624–643 (2016). https://doi.org/10.1016/j.jcp.2016.05.046. http://www.sciencedirect.com/science/article/pii/S0021999116301917
Laiu, M.P.: Positive filtered P\(_n\) method for linear transport equations and the associated optimization algorithm. Ph.D. thesis, University of Maryland, College Park (2016)
Laiu, M.P., Hauck, C.D., McClarren, R.G., O’Leary, D.P., Tits, A.L.: Positive filtered P\(_{N}\) moment closures for linear kinetic equations. SIAM J. Numer. Anal. 54(6), 3214–3238 (2016)
Lebedev, V.: Quadratures on a sphere. Comput. Math. Math. Phys. 16, 10–24 (1976)
LeVeque, R.: Finite difference methods for ordinary and partial differential equations. Soc. Ind. Appl. Math. (2007). https://doi.org/10.1137/1.9780898717839
Lewis, E.E., Miller, W.F.J.: Computational Methods in Neutron Transport. Wiley, New York (1984)
Light, D., Durran, D.: Preserving nonnegativity in discontinuous Galerkin approximations to scalar transport via truncation and mass aware rescaling (TMAR). Mon. Weather Rev. 144(12), 4771–4786 (2016). https://doi.org/10.1175/MWR-D-16-0220.1
Liu, X.D., Osher, S.: Nonoscillatory high order accurate self-similar maximum principle satisfying shock capturing schemes I. SIAM J. Numer. Anal. 33(2), 760–779 (1996)
Liu, Y., Cheng, Y., Shu, C.W.: A simple bound-preserving sweeping technique for conservative numerical approximations. J. Sci. Comput. 73(2), 1028–1071 (2017). https://doi.org/10.1007/s10915-017-0395-x
Loubère, R., Staley, M., Wendroff, B.: The repair paradigm: new algorithms and applications to compressible flow. J. Comput. Phys. 211(2), 385–404 (2006). https://doi.org/10.1016/j.jcp.2005.05.010. www.sciencedirect.com/science/article/pii/S0021999105002706
Markowich, P.A., Ringhofer, C.A., Schmeiser, C.: Semiconductor Equations. Springer, New York (1990)
McClarren, R.G., Hauck, C.D.: Robust and accurate filtered spherical harmonics expansions for radiative transfer. J. Comput. Phys. 229(16), 5597–5614 (2010). https://doi.org/10.1016/j.jcp.2010.03.043. http://www.sciencedirect.com/science/article/pii/S0021999110001622
McClarren, R.G., Holloway, J.P., Brunner, T.A.: On solutions to the \({P_N}\) equations for thermal radiative transfer. J. Comput. Phys. 227(5), 2864–2885 (2008)
Nezza, E.D., Palatucci, G., Valdinoci, E.: Hitchhiker’s guide to the fractional Sobolev spaces. Bull. Sci. Math. 136(5), 521–573 (2012)
Olson, G.L.: Second-order time evolution of \({P_N}\) equations for radiation transport. J. Comput. Phys. 228(8), 3072–3083 (2009). https://doi.org/10.1016/j.jcp.2009.01.012
Perthame, B.: Boltzmann type schemes for gas dynamics and the entropy property. SIAM J. Numer. Anal. 27(6), 1405–1421 (1990)
Perthame, B.: Second-order Boltzmann schemes for compressible euler equations in one and two space dimensions. SIAM J. Numer. Anal. 29(1), 1–19 (1992). https://doi.org/10.1137/0729001
Pomraning, G.C.: Radiation Hydrodynamics. Pergamon Press, New York (1973)
Quarteroni, A.: Some results of Bernstein and Jackson type for polynomial approximation in \({L^p}\)-spaces. Jpn. J. Appl. Math. 1(1), 173–181 (1984). https://doi.org/10.1007/BF03167866
Radice, D., Abdikamalov, E., Rezzolla, L., Ott, C.D.: A new spherical harmonics scheme for multi-dimensional radiation transport I: static matter configurations. J. Comput. Phys. 242(0), 648–669 (2013). https://doi.org/10.1016/j.jcp.2013.01.048. http://www.sciencedirect.com/science/article/pii/S0021999113001125
Shashkov, M., Wendroff, B.: The repair paradigm and application to conservation laws. J. Comput. Phys. 198(1), 265–277 (2004). https://doi.org/10.1016/j.jcp.2004.01.014. www.sciencedirect.com/science/article/pii/S0021999104000270
Walters, W.: Use of the Chebyshev–Legendre quadrature set in discrete-ordinate codes. Technical report, LA-UR-87-3621, Los Alamos National Laboratory (1987)
Zhang, X., Shu, C.: On maximum-principle-satisfying high order schemes for scalar conservation laws. J. Comput. Phys. 229(9), 3091–3120 (2010). https://doi.org/10.1016/j.jcp.2009.12.030. http://www.sciencedirect.com/science/article/pii/S0021999109007165
Zhang, X., Shu, C.W.: Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 467(2134), 2752–2776 (2011). https://doi.org/10.1098/rspa.2011.0153
Author information
Authors and Affiliations
Corresponding author
Additional information
This manuscript has been authored, in part, by UT-Battelle, LLC, under Contract No. DE-AC0500OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for the United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan). M. Paul Laiu: Supported by the U.S. Department of Energy, under the SCGSR program administered by the Oak Ridge Institute for Science and Education under Contract No. DE-AC05-06OR23100. Supported by the U.S. National Science Foundation under Grant No. 1217170. Cory D. Hauck: This author’s research was sponsored by the Office of Advanced Scientific Computing Research and performed at the Oak Ridge National Laboratory, which is managed by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725. We thank Professor Xiangxiong Zhang for pointing out the sweeping limiter in [31] to the authors.
Appendix A: Implementation Details
Appendix A: Implementation Details
1.1 Appendix A.1: The Tensor Product Quadrature
We use a tensor product quadrature to define the positivity limiters and evaluate the numerical flux in the PDE solver. For functions g on \(\mathbb {S}^2\), we write \(g(\varOmega )\) in the polar coordinate as \(g(\mu ,\phi )\), where \(\mu =\varOmega _3\in [-1,1]\) and \(\phi =(0,2\pi ]\) is the azimuthal angle on the sphere. This quadrature rule integrates any integrable g by
Here \(\{\mu _k\}_{k=1}^{N_\mathcal {Q}}\) and \(\{{w}_k\}_{k=1}^{N_\mathcal {Q}}\) are the Gauss–Legendre abscissas and weights, \(\{\phi _j\}_{j=1}^{2N_\mathcal {Q}}\) are equally spaced points from 0 to \(2\pi \), and \(N_\mathcal {Q}\) denotes the degree of precision of \(\mathcal {Q}\). As discussed in Sect. 2.2, the quadrature is required to have degree of precision \(2N+1\) when the approximation is of order N. Thus, a grid of at least \(N+1\) (or \(\lceil (N+1)/2 \rceil \) for even functions on \(\mu \)) Gauss–Legendre points in the \(\mu \) direction and \(2(N+1)\) equally spaced points in the \(\phi \) direction is needed.
1.2 Appendix A.2: The Sweeping Permutation
As mentioned in Sect. 3.2.1, the swp approximation depends on the permutation \(\varPi \). In this appendix, we give the details on the choices of \(\varPi \) used in the numerical tests in Sects. 4.5 and 5.
For functions \(\varphi \) on \([-1,1]\) tested in Sect. 4.5, the regular swp limiter simply chooses \(\varPi \) to be the identity map, i.e., \(\varPi (i)=i\) for all \(i=1,\ldots ,|{Q_N}|\). This choice of \(\varPi \) implies that the Gauss–Legendre quadrature points are sorted as \(\mu _{\varPi (1)}\le \mu _{\varPi (2)}\le \cdots \le \mu _{\varPi (N+1)}\). For the swp-e limiter tested in Sect. 4.5, \(\varPi \) is chosen so that all negative nodal values are first visited in the forward sweeping procedure, followed by the nonnegative nodal values in the ascending order with respect to the associated quadrature weights. Specifically,
For functions \(\varPhi \) on \(\mathbb {S}^2\) in Sects. 4.5 and 5, the tensor-product quadrature points \(\{\varOmega _i=(\mu _k,\phi _j)\}\) are given in “Appendix A.1”. Suppose that \(\{\mu _k\}\) and \(\{\phi _j\}\) are both sorted in the ascending order, i.e., \(\mu _1\le \cdots \le \mu _{\lceil (N+1)/2 \rceil }\), \(\phi _1\le \cdots \le \phi _{2(N+1)}\), and \(\{\varOmega _i=(\mu _k,\phi _j)\}\) is ordered such that \(i=2(k-1)(N+1)+j\). Then the regular swp limiter chooses \(\varPi \) such that for \(\varOmega _i = (\mu _k,\phi _j)\), \(\varPi (i)=i=2(k-1)(N+1)+j\). Such choice of \(\varPi \) sorts \(\{\varOmega _i\}\) first in the ascending order of \(\mu _k\), and then in the ascending order of the \(\phi _j\). On the other hand, the swp-s limiter utilizes the rotational invariance of the line source solution and chooses \(\varPi \) based on the spatial position. Specifically, at position \((x_1,x_2)\), let
In the swp-s limiter, \(\varPi \) is chosen such that for \(\varOmega _i = (\mu _k,\phi _j)\), if \(x_1x_2\ge 0\),
and if \(x_1x_2<0\),
This choice of \(\varPi \) is designed to produce solutions that are symmetric with respect to the quadrants, as discussed in Sect. 5.1.
1.3 Appendix A.3: The Hierarchical Structure
The implementation of h-clp limiter requires a hierarchy of cells on the angular space \(\mathcal {S}\). When \(\mathcal {S}=[-1,1]\), we construct the hierarchy to satisfy (17) via uniform mesh refinements with \(K=3\). Specifically, we start with \(C_{0,1}:=[-1,1]\), and decompose \(C_{0,1}\) evenly into three cells \(C_{1,1}:=[-1, -\frac{1}{3})\), \(C_{1,2}:=[-\frac{1}{3}, \frac{1}{3})\), and \(C_{1,3}:=[\frac{1}{3}, 1]\). Then we continue the refinement by evenly decomposing \(C_{\ell ,k}\) into \(C_{\ell +1, 3k-2}\), \(C_{\ell +1, 3k-1}\), and \(C_{\ell +1, 3k}\), until there exists some \(C_{\ell ,k}\) that contains at most one quadrature point. Similarly, when \(\mathcal {S}=\mathbb {S}^2\), we use the tensor product mesh by applying the uniform mesh for \([-1,1]\) on both the \(\mu \) and \(\phi \) directions, and continue refining the tensor product mesh until some cells contain at most one quadrature point. Specifically, given \(C_{0,1}:=\mathbb {S}^2\), at level \(\ell =1\), \(C_{0,1}\) is decomposed into \(C_{1,1}:=[-1, -\frac{1}{3}) \times [0, \frac{2\pi }{3})\), \(C_{1,2}:=[-1, -\frac{1}{3}) \times [\frac{2\pi }{3}, \frac{4\pi }{3}), \ldots , C_{1,9}:=[\frac{1}{3}, 1] \times [\frac{4\pi }{3}, 2\pi ]\). For level \(\ell >1\), the cells are defined analogously via this refinement.
Rights and permissions
About this article
Cite this article
Laiu, M.P., Hauck, C.D. Positivity Limiters for Filtered Spectral Approximations of Linear Kinetic Transport Equations. J Sci Comput 78, 918–950 (2019). https://doi.org/10.1007/s10915-018-0790-y
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10915-018-0790-y