Kershaw closures for linear transport equations in slab geometry II: high-order realizability-preserving discontinuous-Galerkin schemes

Kershaw closures for linear transport equations in slab geometry II: high-order realizability-preserving discontinuous-Galerkin schemes

Florian Schneider Fachbereich Mathematik, TU Kaiserslautern, Erwin-Schrödinger-Str., 67663 Kaiserslautern, Germany,

This paper provides a generalization of the realizability-preserving discontinuous-Galerkin scheme given in Schneider2015a () to general full-moment 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 high-order methods is demonstrated using numerical convergence tests and non-smooth benchmark problems.

moment models, minimum entropy, Kershaw closures, kinetic transport equation, realizability-preserving, discontinuous-Galerkin scheme
[2010] 35L40, 47B35, 65M08, 65M60, 65M70
journal: arXiv.orgjournal: Journal of Computational Physics

1 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 phase-space 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 so-called 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 Lewis-Miller-1984 (), where is the degree of the highest-order 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 non-negative. Thus solutions can contain negative values for the local densities of particles, rendering the solution physically meaningless.

Entropy-based 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,111 Positivity is actually not gained for every entropy-based 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 space-time 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 higher-order 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 so-called 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 high-order numerical solutions (in space and time) to the Euler equations, which indeed are an entropy-based 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 high-order method for hyperbolic systems is the Runge-Kutta 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 “positivity-preserving”) limiter Zhang2010 (), but so far these have been implemented for low-order moment systems (that is or ) Olbrant2012 () because here one can rely on the simplicity of the structure of the realizable set for low-order moments. For moments of large order , the realizable set has complex nonlinear boundaries: when the velocity domain is one-dimensional, the realizable set is characterized by the positive-definiteness of Hankel matrices Shohat1943 (); Curto1991 (); in higher dimensions, the realizable set is not well understood. In Schneider2015a (), using that a quadrature-based 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 minimum-entropy ansatz a new hierarchy of full-moment 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 minimum-entropy 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 discontinuous-Galerkin scheme is given with the necessary extensions to obtain a realizability-preserving 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 high-order space-time approximations. Finally, conclusions and an outlook on future work is given.

2 Modelling

In slab geometry, the transport equation under consideration has the form


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


where is non-negative, symmetric in both arguments and normalized to . In this paper the special case of the BGK-type isotropic-scattering operator with is used for the simulations.

(2.1) is supplemented by initial and boundary conditions:


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 high-dimensional equation into a system of lower-dimensional 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 so-called moments of a given distribution function with respect to are then defined by


where the integration is performed componentwise.

Assuming for simplicity , the quantity is called local particle density. Furthermore, normalized moments are defined as


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


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 infinite-dimensional 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 moment-space given by . This is the so-called 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 short-hand notation is neglected for notational simplicity.

In this paper, the full-moment 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 non-negative, a moment vector only makes sense physically if it can be associated with a non-negative 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 semi-definite. In particular denotes that is positive semi-definite.

For the full-moment 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


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


where the interpolation constant


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 first-order system of balance laws




4 Realizability-preserving discontinuous-Galerkin scheme

Recent numerical experiments have shown that high-order schemes (in space and time) outperform highly-resolved first-order methods, comparing degrees of freedom and running time versus approximation quality. This has been investigated in the case of minimum-entropy 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, higher-order 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 minimum-entropy 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 high-order kinetic scheme presented in Schneider2015b () is not obvious. This has been observed before in Vikas2011 () for quadrature-based moment methods. Therefore this paper focuses on the discontinuous-Galerkin scheme presented in Schneider2015a (). While there only quadrature-based minimum-entropy models have been investigated, the following sections will show how to generalize the scheme and its realizability limiter to the general case of full-moment 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


is the finite-element space of piecewise polynomials of degree .

The discontinuous-Galerkin 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 finite-element 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


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 cell-interfaces, the fluxes at the points of discontinuity are both replaced by a numerical flux , thus coupling the elements with their neighbours Toro2009 (). Several well-known examples for such a numerical flux exist in literature. The simplest example is the global Lax-Friedrichs flux


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 one222The results in Schneider2015 () prove this for but there is no general proof of this fact for arbitrary Kershaw closures yet..

The local Lax-Friedrichs flux could be used instead. This requires computing the eigenvalues of the Jacobian in every space-time cell to adjust the value of the numerical viscosity constant but possibly decreases the overall diffusivity of the scheme. However, since high-order space-time approximations are considered, the decrease in diffusivity achieved by switching to the local Lax-Friedrichs flux should be negligible.

The usual approach is to expand the approximate solution on each interval as


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


With an orthogonal basis the cell means are easily available from the expansion coefficients , since


Collecting the coefficients into the matrix


equation (4.2) can be written in compact form as the coupled system of ordinary differential equations


with initial condition (4.2b) and an appropriate choice of the local differential operator Schneider2015a ().

The incorporation of boundary conditions for moment systems is non-trivial. Here, an often-used 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


Note, however, that the validity of this approach, due to its inconsistency with the original boundary conditions (2.4b)–(2.4c), is not entirely non-controversial, 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 Dirichlet-boundary conditions, the simplest approach is taken. The ghost-cell 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 high-order scheme in space and time is a suitable time integrator for (4.8). Such a class of integrators is given by the strong-stability-preserving (SSP) methods, as used for example in Zhang2010 (); AllHau12 (). The stages and steps of these type of methods are convex combinations of forward-Euler steps. Since the realizable set is convex, the analysis of a forward-Euler step then suffices to prove realizability preservation of the full method.

When possible, SSP-Runge-Kutta (SSP-RK) methods are used, but unfortunately they only exist up to order four Ruuth2004 (); Gottlieb2005 (). For orders the so-called two-step Runge-Kutta (TSRK) SSP methods Ketcheson2011 () as well as their generalizations, the multi-step Runge-Kutta (MSRK) SSP methods Bresten2013 () can be applied. They combine Runge-Kutta schemes with positive weights and high-order multistep methods to achieve a total order higher than four while maintaining the important SSP property.

See Schneider2015b () for more information about the SSP-schemes 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


for the cell and the case , that is piece-wise quadratic or higher-degree polynomials, so that the final rows of zeros in the first case indicates that the coefficients for the higher-order 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 problem-dependent 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 non-physical 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 flux-term at the spatial quadrature nodes in the cell, at least for each node is necessary333Although 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


This assumption is fulfilled by the integral collision operator (2.3).

While first-order schemes (like the Lax-Friedrichs method) automatically preserve realizability of the cell means Schneider2015a (), higher-order 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 minimum-entropy 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 CFL-type condition. With realizable cell means available, a point-wise-realizable 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

  1. for all cells it holds that , ;

  2. the cell means at time step are realizable;

  3. at the quadrature nodes of the -point Gauss-Lobatto rule on each cell , ,444Where is the ceiling function, that is, it returns smallest integer bigger than or equal to its argument. Since the Gauss-Lobatto rule is exact for polynomials of degree this choice guarantees to exactly integrate the occurring polynomials of degree . the point-wise values of the moment approximation (componentwise in ) are realizable.

Then under the CFL condition


the cell means after one forward-Euler step are realizable, where


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 linear-scaling 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 minimum-entropy 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


is realizable. Indeed, by inserting the definition of from above, the limited moment vector can be written as

thus when limiting is necessary, the higher-order 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 non-linear 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


requires to calculate the roots of two polynomials of degree at most .


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 non-negative, 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 second-order realizability condition with equality. This example is visualized in Figure 1.

Figure 1: Limiter example for the second-order basis with . The realizable set is plotted in grey.
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 linearly-spaced 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 point-wise 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 often-used benchmark problems for moment models. They serve as a reference for the efficiency of Kershaw closures with a high number of moments combined with high-order space-time 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 first-order 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


respectively. The integral in is approximated using a 100-point Gauss-Lobatto quadrature rule over each spatial cell , and is approximated by taking the maximum over these quadrature nodes. The observed convergence order is defined by


where , , , is the -error for the numerical solution using cell size .

A convergence table for orders is presented in Table 1.

10 7. 721e-02 1. 418e-04 4. 449e-06 1. 235e-07 2. 595e-09
20 2. 168e-02 1.8 9. 004e-06 4.0 1. 430e-07 5.0 1. 898e-09 6.0 2. 107e-11 6.9
40 1. 017e-02 1.1 5. 788e-07 4.0 4. 465e-09 5.0 2. 951e-11 6.0 1. 715e-13 6.9
80 2. 580e-03 2.0 3. 667e-08 4.0 1. 401e-10 5.0 4. 629e-13 6.0 5. 961e-14 1.5
160 6. 467e-04 2.0 2. 296e-09 4.0 4. 408e-12 5.0 3. 028e-14 3.9 1. 159e-13 -1.0
10 5. 978e-02 1. 864e-04 6. 520e-06 1. 750e-07 4. 422e-09
20 1. 592e-02 1.9 1. 155e-05