Kershaw closures for linear transport equations in slab geometry II: highorder realizabilitypreserving discontinuousGalerkin schemes
Abstract
This paper provides a generalization of the realizabilitypreserving discontinuousGalerkin scheme given in Schneider2015a () to general fullmoment models that can be closed analytically. It is applied to the class of Kershaw closures, which are able to provide a cheap closure of the moment problem. This results in an efficient algorithm for the underlying linear transport equation. The efficiency of highorder methods is demonstrated using numerical convergence tests and nonsmooth benchmark problems.
keywords:
moment models, minimum entropy, Kershaw closures, kinetic transport equation, realizabilitypreserving, discontinuousGalerkin schemeMsc:
[2010] 35L40, 47B35, 65M08, 65M60, 65M701 Introduction
Moment closures are a class of spectral methods used in the context of kinetic transport equations. An infinite set of moment equations is defined by taking velocity or phasespace averages with respect to some basis of the velocity space. A reduced description of the kinetic density is then achieved by truncating this hierarchy of equations at some finite order. The remaining equations however inevitably require information from the equations which were removed. The specification of this information, the socalled moment closure problem, distinguishes different moment methods. In the context of linear radiative transport, the standard spectral method is commonly referred to as the closure LewisMiller1984 (), where is the degree of the highestorder moments in the model. The method is powerful and simple to implement, but does not take into account the fact that the original function to be approximated, the kinetic density, must be nonnegative. Thus solutions can contain negative values for the local densities of particles, rendering the solution physically meaningless.
Entropybased moment closures, referred to as models in the context of radiative transport Min78 (); DubFeu99 (), have all the properties one would desire in a moment method, namely positivity of the underlying kinetic density,^{1}^{1}1 Positivity is actually not gained for every entropybased moment closure but is indeed a property of those models derived from important, physically relevant entropies. hyperbolicity of the closed system of equations, and entropy dissipation Lev96 (). Practical implementation of these models has been traditionally considered too expensive because they require the numerical solution of an optimization problem at every point on the spacetime grid, but recently there has been renewed interest in the models due to their inherent parallelizability Hauck2010 (). However, while their parallelizability goes a long way in making models computationally competitive, in order to make these methods truly competitive with more basic discretizations, the gains in efficiency that come from higherorder methods (in space and time) will likely be necessary. Here the issue of realizability becomes a stumbling block.
The property of positivity implies that the system of moment equations only evolves on the set of socalled realizable moments. Realizable moments are simply those moments associated with positive densities, and the set of these moments forms a convex cone which is a strict subset of all moment vectors. This property, while indeed desirable since it is consistent with the original kinetic density, can cause problems for numerical methods. Standard highorder numerical solutions (in space and time) to the Euler equations, which indeed are an entropybased moment closure, have been observed to have negative local densities and pressures Zhang2010 (). Similar effects have been reported in the context of elastic flow Schar1996 (). This is exactly loss of realizability.
A recently popular highorder method for hyperbolic systems is the RungeKutta discontinuous Galerkin (RKDG) method Cockburn1989 (); Cockburn1989a (); Cockburn1991 (). An RKDG method for moment closures can handle the loss of realizability through the use of a realizability (or “positivitypreserving”) limiter Zhang2010 (), but so far these have been implemented for loworder moment systems (that is or ) Olbrant2012 () because here one can rely on the simplicity of the structure of the realizable set for loworder moments. For moments of large order , the realizable set has complex nonlinear boundaries: when the velocity domain is onedimensional, the realizable set is characterized by the positivedefiniteness of Hankel matrices Shohat1943 (); Curto1991 (); in higher dimensions, the realizable set is not well understood. In Schneider2015a (), using that a quadraturebased approximation of the realizable set is a convex polytope Alldredge2014 (), the realizability limiters of Zhang2010 (); Olbrant2012 () has been generalized for moment systems of (in principle) arbitrary order.
To avoid the expensive minimumentropy ansatz a new hierarchy of fullmoment models has been derived in Schneider2015 (), the class of Kershaw closures, based on the findings in Ker76 (). It provides a reasonably simple closure relation, closely related to minimumentropy models while being cheap to evaluate. This paper aims at generalizing the scheme given in Schneider2015a () to this class of models for (in principle) arbitrary moment order .
This paper is organized as follows. First, the transport equation and its moment approximations are given. Then, the available realizability theory is shortly reviewed, followed by a brief summary of the class of Kershaw closures. The discontinuousGalerkin scheme is given with the necessary extensions to obtain a realizabilitypreserving scheme. Numerical convergence of this scheme up to seventh order against an analytical solution is shown and the Kershaw closures are submitted to a set of benchmark tests investigating the effect of highorder spacetime approximations. Finally, conclusions and an outlook on future work is given.
2 Modelling
In slab geometry, the transport equation under consideration has the form
(2.1) 
The physical parameters are the absorption and scattering coefficient , respectively, and the emitting source . Furthermore, , and .
The shorthand notation denotes integration over .
Assumption 2.1.
Following Levermore1996 (), the collision operator is assumed to have the following properties.
A typical example for is the linear integral operator
(2.3) 
where is nonnegative, symmetric in both arguments and normalized to . In this paper the special case of the BGKtype isotropicscattering operator with is used for the simulations.
(2.1) is supplemented by initial and boundary conditions:
(2.4a)  
(2.4b)  
(2.4c) 
3 Moment models and realizability
In general, solving equation (2.1) is very expensive in two and three dimensions due to the high dimensionality of the state space.
For this reason it is convenient to use some type of spectral or Galerkin method to transform the highdimensional equation into a system of lowerdimensional equations. Typically, one chooses to reduce the dimensionality by representing the angular dependence of in terms of some basis .
Definition 3.1.
The vector of functions consisting of basis functions , of maximal order is called an angular basis.
The socalled moments of a given distribution function with respect to are then defined by
(3.1) 
where the integration is performed componentwise.
Assuming for simplicity , the quantity is called local particle density. Furthermore, normalized moments are defined as
(3.2) 
where is the remainder of the basis .
To obtain a set of equations for , (2.1) has to be multiplied through by and integrated over , giving
Collecting known terms, and interchanging integrals and differentiation where possible, the moment system has the form
(3.3) 
The solution of (3.3) is equivalent to the one of (2.1) if is a basis of .
Since it is impractical to work with an infinitedimensional system, only a finite number of basis functions of order can be considered. Unfortunately, there always exists an index such that the components of are not in the linear span of . Therefore, the flux term cannot be expressed in terms of without additional information. Furthermore, the same might be true for the projection of the scattering operator onto the momentspace given by . This is the socalled closure problem. One usually prescribes some ansatz distribution to calculate the unknown quantities in (3.3). Note that the dependence on the angular basis in the shorthand notation is neglected for notational simplicity.
In this paper, the fullmoment monomial basis is considered. However, it is in principle possible to extend the derived concepts to other bases like half DubKla02 (); DubFraKlaTho03 () or mixed moments Frank07 (); Schneider2014 ().
The rest of this section is a brief summary of the corresponding parts in Schneider2015 (). All details and further discussions can be found therein.
3.1 Realizability
Since the underlying kinetic density to be approximated is nonnegative, a moment vector only makes sense physically if it can be associated with a nonnegative distribution function. In this case the moment vector is called realizable.
Definition 3.2.
The realizable set is
If , then is called realizable. Any such that is called a representing density. If is additionally a linear combination of Dirac deltas Hassani2009 (); Mathematics2011 (); Kuo2006 (), it is called atomic Curto1991 ().
Definition 3.3.
Let be Hermitian matrices. The partial ordering on such matrices is defined by if and only if is positive semidefinite. In particular denotes that is positive semidefinite.
For the fullmoment basis the question of finding practical characterizations of the realizable set has been completely solved in Curto1991 (). See Schneider2015 () for more details. The following characterizations of the above realizable set holds.
Lemma 3.4.
Define the Hankel matrices
Then the realizable set satisfies
Due to the structure of the used Hankel matrices (the highest moment always appears exactly once in the entries of the matrices) it is always possible to rearrange the conditions involving this highest moment in Theorem 3.4 in such a way that
for functions and . Whenever is realizable and , is said to be on the lower order realizability boundary. Similarly, if , is said to be on the upper order realizability boundary.
The functions and can be specified using the pseudoinverses of combinations of Hankel matrices. To simplify notation later, the following corollary is written in terms of instead of .
Corollary 3.5.
The functions and satisfying
(3.4) 
are given by
where in the odd case
and in the even case
Remark 3.6.
By convention, if .
3.2 Kershaw closures
With the previous realizability theory it is now possible to develop the closure strategy which is called Kershaw closure. This class of moment models is defined by convexly combining upper and lower moments of order in such a way that the isotropic point is correctly reproduced.
Corollary 3.7.
The Kershaw closure of order is given by
(3.5) 
where the interpolation constant
(3.6) 
is defined via the functions and as given in Corollary 3.5 and .
For convenience, (3.3) using the Kershaw closure can be written in the form of a usual firstorder system of balance laws
(3.7) 
where
(3.8a)  
(3.8b) 
4 Realizabilitypreserving discontinuousGalerkin scheme
Recent numerical experiments have shown that highorder schemes (in space and time) outperform highlyresolved firstorder methods, comparing degrees of freedom and running time versus approximation quality. This has been investigated in the case of minimumentropy moment models in Schneider2015a (); Schneider2015b () for two different types of schemes. The most challenging part is to preserve realizability during the simulation since otherwise the closure cannot be evaluated. Unfortunately, higherorder schemes typically cannot guarantee this property on their own, as has been observed in the context of the compressible Euler equations (which are indeed in the hierarchy of minimumentropy models) in Zhang2010 () and for the model in Olbrant2012 ().
Due to the lack of smoothness in the underlying distribution of the Kershaw models (since it is atomic) the application of the highorder kinetic scheme presented in Schneider2015b () is not obvious. This has been observed before in Vikas2011 () for quadraturebased moment methods. Therefore this paper focuses on the discontinuousGalerkin scheme presented in Schneider2015a (). While there only quadraturebased minimumentropy models have been investigated, the following sections will show how to generalize the scheme and its realizability limiter to the general case of fullmoment models.
In the following, the spatial domain is divided into (for notational simplicity) (equidistant) cells , where the cell interfaces are given by for cell centres , and .
Furthermore, is the set of polynomials of degree at most on the interval , and
(4.1) 
is the finiteelement space of piecewise polynomials of degree .
The discontinuousGalerkin method for the general hyperbolic system (3.7), as outlined in Cockburn1989 (); Cockburn1989a (); Cockburn1991 (), can be briefly described as follows.
For each , seek an approximate solution whose components live in the finiteelement space as defined in (4.1).
Then follow the Galerkin approach: replace in (3.7) by a solution of the form , multiply the resulting equation by basis functions of and integrate over cell to obtain
(4.2a)  
(4.2b) 
where and again denote the limits from left and right, respectively, and is the projection of the initial distribution to the moment space. In order to approximately solve the Riemann problem at the cellinterfaces, the fluxes at the points of discontinuity are both replaced by a numerical flux , thus coupling the elements with their neighbours Toro2009 (). Several wellknown examples for such a numerical flux exist in literature. The simplest example is the global LaxFriedrichs flux
(4.3) 
The numerical viscosity constant is taken as the global estimate of the absolute value of the largest eigenvalue of the Jacobian . Following Schneider2015 (), the viscosity constant can be set to , because for the moment systems used here it can be shown that the largest eigenvalue is bounded in absolute value by one^{2}^{2}2The results in Schneider2015 () prove this for but there is no general proof of this fact for arbitrary Kershaw closures yet..
The local LaxFriedrichs flux could be used instead.
This requires computing the eigenvalues of the Jacobian in every
spacetime cell to adjust the value of the numerical viscosity constant
but possibly decreases the overall diffusivity of the scheme.
However, since highorder spacetime approximations are considered, the decrease in
diffusivity achieved by switching to the local LaxFriedrichs flux should be
negligible.
The usual approach is to expand the approximate solution on each interval as
(4.4) 
where denote a basis for with respect to the standard scalar product on the reference cell . It is convenient to choose an orthogonal basis like the Legendre polynomials scaled to the interval , denoted by
(4.5) 
With an orthogonal basis the cell means are easily available from the expansion coefficients , since
(4.6) 
Collecting the coefficients into the matrix
(4.7) 
equation (4.2) can be written in compact form as the coupled system of ordinary differential equations
(4.8) 
with initial condition (4.2b) and an appropriate choice of the local differential operator Schneider2015a ().
The incorporation of boundary conditions for moment systems is nontrivial. Here, an oftenused approach is taken that incorporates boundary conditions via ‘ghost cells’. First assume that it is possible to smoothly extend in to for (note that while moments are defined using integrals over all , the boundary conditions in (2.4b)–(2.4c) are only defined for corresponding to incoming data).
Then the moment approximations in the ghost cells at and simply take the form
(4.9a)  
(4.9b) 
Note, however, that the validity of this approach, due to its inconsistency with the original boundary conditions (2.4b)–(2.4c), is not entirely noncontroversial, but the question of appropriate boundary conditions for moment models is an open problem pomraning1964variational (); Larsen1991 (); Rulko1991 (); Struchtrup2000 (); levermore2009boundary () which is not explored here.
For Dirichletboundary conditions, the simplest approach is taken. The ghostcell moments are chosen to be the constant functions
with and defined as in (4.9).
For periodic boundary conditions, the obvious choice is
All that remains to obtain a highorder scheme in space and time is a suitable time integrator for (4.8). Such a class of integrators is given by the strongstabilitypreserving (SSP) methods, as used for example in Zhang2010 (); AllHau12 (). The
stages and steps of these type of methods are convex combinations of forwardEuler steps.
Since the realizable set is convex, the analysis of a forwardEuler step
then suffices to prove realizability preservation of the full method.
When possible, SSPRungeKutta (SSPRK) methods are used, but unfortunately they only exist up to order four Ruuth2004 (); Gottlieb2005 (). For orders the socalled twostep RungeKutta (TSRK) SSP methods Ketcheson2011 () as well as their generalizations, the multistep RungeKutta (MSRK) SSP methods Bresten2013 () can be applied. They combine RungeKutta schemes with positive weights and highorder multistep methods to achieve a total order higher than four while maintaining the important SSP property.
See Schneider2015b () for more information about the SSPschemes used in the actual implementation. Note that they differ from those used in Schneider2015a () where only discretizations up to third order were used, in contrast to the methods of order one to seven given in Schneider2015b ().
The rest of the methodology follows closely Schneider2015a (). Here, the standard TVBM corrected minmod limiter proposed in Cockburn1989a () is used.
Assuming that the major part of the spurious oscillations is generated in the linear part of the underlying polynomial, whose slope in the reference cell is simply , a limiter can be defined as
(4.10) 
for the cell and the case , that is piecewise quadratic or higherdegree polynomials, so that the final rows of zeros in the first case indicates that the coefficients for the higherorder spatial basis functions are set to zero for each moment component. The absolute value and the inequality are applied componentwise. The label “scalar” is used because the limiter is directly applied to each scalar component of . The function is the standard minmod function applied componentwise, defined by
The constant is a problemdependent estimate of the second derivative, though it has to be noted that in Cockburn1989a () the authors did not find the solutions very sensitive to the value chosen for this parameter.
However, it has been found that applying the limiter to the components themselves may introduce nonphysical oscillations around an otherwise monotonic solution Cockburn1989 (). Instead, the limiter is applied to the local characteristic fields of the solution. The flux Jacobian is obtained numerically using finite differences. This can be achieved cheaply since only the last equation of the flux has to be considered, all other components are trivial.
4.1 Realizability preservation
In order to evaluate the fluxterm at the spatial quadrature nodes in the cell, at least for each node is necessary^{3}^{3}3Although intuition expects for all , having realizable point values only indeed suffices to preserve realizability of the updated cell means..
To prove Theorem 4.2 the following rather strong assumption has to be made.
Assumption 4.1.
For every satisfying there exists a such that the moments of the collision operator with respect to the same angular basis can be written as
(4.11) 
This assumption is fulfilled by the integral collision operator (2.3).
While firstorder schemes (like the LaxFriedrichs method) automatically preserve realizability of the cell means Schneider2015a (), higherorder schemes () typically cannot guarantee this property on their own, as
has been
observed in the context of the compressible Euler equations (which are
indeed in the hierarchy of minimumentropy models) in Zhang2010 () and for the model in Olbrant2012 ().
It is, however, possible to show that, when the moments at the quadrature nodes are realizable, the presented schemes preserve realizability of the cell means under a CFLtype condition. With realizable cell means available, a pointwiserealizable polynomial representation can be obtained by applying a linear scaling limiter which pushes towards the cell mean and thus into the realizable set for each quadrature node .
Following the arguments in Zhang2010a (); Zhang2011 (), this limiter does not destroy the accuracy of the scheme in case of smooth solutions if is not on the boundary of the realizable set. This is verified numerically in Section 5.1. For convenience, the main result of Schneider2015a (); Schneider2015b () is summarized in the following theorem. Note that, since SSP time integrators are used, it suffices to investigate forward Euler steps in time, which are then convexly combined to obtain the designed order of time integration in the SSP scheme.
Theorem 4.2 (Schneider2015a (); Schneider2015b ()).
Assume that

for all cells it holds that , ;

the cell means at time step are realizable;

at the quadrature nodes of the point GaussLobatto rule on each cell , ,^{4}^{4}4Where is the ceiling function, that is, it returns smallest integer bigger than or equal to its argument. Since the GaussLobatto rule is exact for polynomials of degree this choice guarantees to exactly integrate the occurring polynomials of degree . the pointwise values of the moment approximation (componentwise in ) are realizable.
Then under the CFL condition
(4.12) 
the cell means after one forwardEuler step are realizable, where
(4.13) 
All that remains is to ensure that assumption (iii) in Theorem 4.2 is always fulfilled. Due to assumption (ii) and the convexity of the realizable set, this can be achieved using a linearscaling limiter, pushing the polynomial representation towards the (realizable) cell mean. This approach has been outlined in Zhang2010 (); Zhang2010a (); Zhang2011a () for the Euler equations and in Schneider2015a () for two classes of minimumentropy models.
For ease of notation, time indices are dropped.
Recall the definition (4.4) of , given by
Due to the convexity of the realizable set, if is realizable, then for each quadrature point there exists a such that
(4.14) 
is realizable. Indeed, by inserting the definition of from above, the limited moment vector can be written as
thus when limiting is necessary, the higherorder coefficients , ,
are damped while the cell mean remains unchanged.
The task of the limiter is now to choose for each the minimal value of such that is realizable at all quadrature nodes . This choice is optimal in the sense that the least information of the original polynomial is lost ( corresponds to no limiting while resembles limiting to first order).
Remark 4.3.
For readability reasons, the dependence on the cell index is dropped sometimes throughout the following examples.
Having the nonlinear structure of the full realizable set in mind, computing the smallest such that requires some effort.
Theorem 4.4.
The solution to the limiter problem
s.t.  
requires to calculate the roots of two polynomials of degree at most .
Proof.
Assume that . Define , , and to be the Hankel matrices associated with and , respectively. Then the Hankel matrices associated with are
By assumption (compare Lemma 3.4), . This implies that all eigenvalues of are nonnegative, and therefore . Being on the realizability boundary corresponds to having at least one zero eigenvalue, which is equivalent to a vanishing determinant of either or . Note that is a polynomial of degree in . Since the realizable set is convex, the maximal in that is a root of one of the two polynomials is the optimal limiter value.
The case works analogously. ∎
Example 4.5.
Realizability conditions for are very simple: . Plugging in from (4.14) gives
Solving these equations for equality (which is equivalent to finding roots of a polynomial of degree ) results in
Example 4.6.
For the realizability conditions are given through the Hankel matrices
and the conditions and . The matrices defining the limiter value are given by
The required polynomials are given by , i.e.
Let e.g. and . Then it follows that
which have roots , and . Since it follows that and , which indeed satisfies the secondorder realizability condition with equality. This example is visualized in Figure 1.
Remark 4.7.
To analytically obtain the coefficients for the resulting polynomials in is in general hard. However, this can be avoided using a simple trick. In the odd case evaluate the determinants of at distinct values of , e.g. at the linearlyspaced values in . These points uniquely define the desired polynomial. A similar approach can be done in the even case.
Remark 4.8.
It has been shown in Schneider2016 () that the slope limiter (4.10) (either evaluated in primitive or conserved variables) always has to be applied before the realizability limiter since the application of the slope limiter may destroy pointwise realizability (and thus Theorem 4.2 cannot be applied).
Both limiters have to be applied at every stage and step of the SSP time integrators.
5 Numerical experiments
This section contains numerical convergence results and some oftenused benchmark problems for moment models. They serve as a reference for the efficiency of Kershaw closures with a high number of moments combined with highorder spacetime discretizations.
The approximations of highest order in space and time are discretized with on a grid with cells, the medium order is represented by a solution on a grid with cells and the firstorder variant is calculated on a grid with cells. If not stated otherwise the TVB constant in the modified minmod limiter is set to for and for (compare Cockburn1989a (); Schneider2015a ()).
5.1 Convergence results
5.1.1 Manufactured solution
In general, obtaining analytical solutions for moment models is a hard task. In the case of Kershaw closures it is possible to provide a solution in some special cases. Consider the initial distribution
with some positive . Setting , the analytical solution of the transport equation (2.1) is given by
On the moment level this corresponds to a linear advection with transport speed since
Since can be reproduced exactly by the Kershaw closures Schneider2015 () the moments of the transport solution are also the moments of the Kershaw closure if . For this example the local mass is defined as on . The final time is set to and periodic boundary conditions are applied.
Errors are computed in the zeroth moment of the solution . Then  and errors for the zeroth moment (that is, the zeroth component of a numerical solution ) are defined as
(5.1) 
respectively. The integral in is approximated using a 100point GaussLobatto quadrature rule over each spatial cell , and is approximated by taking the maximum over these quadrature nodes. The observed convergence order is defined by
(5.2) 
where , , , is the error for the numerical solution using cell size .
A convergence table for orders is presented in Table 1.
10  7.  721e02  —  1.  418e04  —  4.  449e06  —  1.  235e07  —  2.  595e09  — 
20  2.  168e02  1.8  9.  004e06  4.0  1.  430e07  5.0  1.  898e09  6.0  2.  107e11  6.9 
40  1.  017e02  1.1  5.  788e07  4.0  4.  465e09  5.0  2.  951e11  6.0  1.  715e13  6.9 
80  2.  580e03  2.0  3.  667e08  4.0  1.  401e10  5.0  4.  629e13  6.0  5.  961e14  1.5 
160  6.  467e04  2.0  2.  296e09  4.0  4.  408e12  5.0  3.  028e14  3.9  1.  159e13  1.0 
10  5.  978e02  —  1.  864e04  —  6.  520e06  —  1.  750e07  —  4.  422e09  — 
20  1.  592e02  1.9  1.  155e05 