First-order quarter- and mixed-moment realizability theory and Kershaw closures for a Fokker-Planck equation in two space dimensions

First-order quarter- and mixed-moment realizability theory and Kershaw closures for a Fokker-Planck equation in two space dimensions

Florian Schneider Jochen Kall Andreas Roth Fachbereich Mathematik, TU Kaiserslautern, Erwin-Schrödinger-Str., 67663 Kaiserslautern, Germany, schneider@mathematik.uni-kl.de Fachbereich Mathematik, TU Kaiserslautern, Erwin-Schrödinger-Str., 67663 Kaiserslautern, Germany, kall@mathematik.uni-kl.de Fachbereich Mathematik, TU Kaiserslautern, Erwin-Schrödinger-Str., 67663 Kaiserslautern, Germany, roth@mathematik.uni-kl.de
Abstract

Mixed-moment models, introduced in Frank07 (); Schneider2014 () for one space dimension, are a modification of the method of moments applied to a (linear) kinetic equation, by choosing mixtures of different partial moments. They are well-suited to handle such equations where collisions of particles are modelled with a Laplace-Beltrami operator. We generalize the concept of mixed moments to two dimension. The resulting hyperbolic system of equations has desirable properties, removing some drawbacks of the well-known model. We furthermore provide a realizability theory for a first-order system of mixed moments by linking it to the corresponding quarter-moment theory. Additionally, we derive a type of Kershaw closures for mixed- and quarter-moment models, giving an efficient closure (compared to minimum-entropy models). The derived closures are investigated for different benchmark problems.

keywords:
radiation transport, moment models, realizability, mixed moments, Laplace-Beltrami operator
Msc:
[2010] 35L40, 35Q84, 65M08, 65M70
journal: arXiv.org

1 Introduction

The full discretization of kinetic transport equations like the Fokker-Planck equation is in general very expensive since the discretized variable resides in where and denotes the unit sphere in . Thus, the solution of the Fokker-Planck equation is very high-dimensional.

A common approach to reduce the dimensionality is given by the method of moments Eddington (); Lev96 (). One chooses a set of angular basis functions, tests the Fokker-Planck equation with it and integrates over the angular variable, removing the angular dependence while getting a system of equations in and . Assuming now a specific form of the underlying distribution function, different approximate systems arise. Typical examples are the well-known spherical harmonics or models Jea17 (); Eddington (); Brunner2005 () and their simplifications, the Gel61 () methods. These models are computationally inexpensive since they form an analytically closed system of hyperbolic differential equations. However, they suffer from severe drawbacks: The methods are generated by closing the balance equations with a distribution function which is a polynomial in the angular variable. This implies that this distribution function might be negative resulting in non-physical values like a negative particle density. Additionally, in many cases a very high number of moments is needed for a reasonable approximation of the transport solution. This is in particular true in beam cases, where the exact transport solution forms a Dirac delta. The entropy minimization -models Min78 (); DubFeu99 (); BruHol01 (); Monreal2008 (); AllHau12 () are expected to overcome this problem since their closure functions are always positive. In many situations these models perform very well. Still they produce unphysical steady-state shocks due to a zero netflux problem.

To improve this situation, half or partial moment models have been introduced in DubKla02 (); FraDubKla04 (). These models work especially well in one space dimension, because they capture the potential discontinuity of the probability density in the angular variable which in 1D is well-located. If however, a Laplace-Beltrami operator is used instead of the standard integral scattering operator, i.e. scattering is extremely forward-peaked Pom92 (), these half moment approximations fail Frank07 (); Schneider2014 ().

To improve this situation a new model with mixed moments was proposed in Frank07 (); Schneider2014 () which is able to avoid this problem. Instead of choosing full or half moments, a mixture of both is used. Contrary to a typical half moment approximation, the lowest order moment (density) is kept as a full moment while all higher moments are averaged over half-spaces. This ensures the continuity of the underlying distribution function.

Since in one spatial dimension mixed moments perform very well, it seems reasonable to extend them to multi-. It turns out, that here a mixture of full, half and quarter moments is necessary to derive the correct mixed-moment ansatz.

Realizability is the fact that a vector of moments is physically relevant, i.e. that it is the moment of a non-negative distribution function. While in 1D realizability theory for full moments Curto1991 () and mixed moments Schneider2014 () is completely solved, it remains an open problem in higher dimensions. We will give a realizability theory for quarter and mixed moments of order and derive, as in 1D, a corresponding Kershaw closure which provides an analytically closed system of equations, in contrast to minimum-entropy models which requires the solution of a nonlinear system of equations.

The paper is organized as follows: First, we will give a short introduction to the method of moments and its consequences for the Fokker-Planck equation. After setting up the basic notations in Section 2, we provide necessary and sufficient conditions of order for realizability in the case of quarter moments and mixed moments in Section 3. Section 4 deals with the construction of Kershaw closures. Finally, we present some numerical tests for the derived models in Section 5 as well as conclusions and outlook in Section 6.

2 Macroscopic Models

We consider the Fokker-Planck equation

(2.1a)
which describes the densities of particles with speed at position and time under the events of scattering (proportional to ), absorption (proportional to ) and emission (proportional to ). The equation is supplemented with initial condition and Dirichlet boundary conditions:
(2.1b)
(2.1c)
where is the outward unit normal vector in .

Similar to Seibold2012 () we assume that geometry, initial and boundary conditions are independent of the -direction, resulting in a solution which is also -independent. Therefore (2.1) can be reduced to . Parametrizing in spherical coordinates and taking the symmetry reduction into account we obtain

(2.2)

where is the azimuthal and the cosine of the polar angle. Then the Laplace-Beltrami operator on the unit sphere can be written as

(2.3)
Definition 2.1.

The vector of functions consisting of basis functions , of maximal order (in ) is called an angular basis.

The so-called moments of a given distribution function are then defined by

(2.4)

where the integration is performed component-wise.

Equations for can then be obtained by multiplying (2.1) with and integration over :

Collecting known terms, and interchanging integration and differentiation where possible, the moment system has the form

(2.5)

Depending on the choice of the terms , and in some cases even cannot be given explicitly in terms of . Therefore an ansatz has to be made for closing the unknown terms. This is called the moment-closure problem.

In this paper the ansatz density is reconstructed from the moments by minimizing the entropy-functional

(2.6)

under the moment constraints

(2.7)

The kinetic entropy density is strictly convex and twice continuously differentiable and the minimum is simply taken over all functions such that is well defined. The obtained ansatz , solving this constrained optimization problem, is given by

(2.8)

This problem, which must be solved over the space-time mesh, is typically solved through its strictly convex finite-dimensional dual,

(2.9)

where is the Legendre dual of . The first-order necessary conditions for the multipliers show that the solution to (2.8) has the form

(2.10)

where is the derivative of .

This approach is called the minimum-entropy closure Levermore1996 ().

The kinetic entropy density can be chosen according to the physics being modelled. As in Levermore1996 (); Hauck2010 (), Maxwell-Boltzmann entropy

(2.11)

is used, thus . This entropy is used for non-interacting particles as in an ideal gas or an ensemble of photons.

A closed system of equations for remains after substituting in (2.5) with :

(2.12)

Also note that using the entropy the linear ansatz

(2.13)

remains. If the angular basis is chosen as spherical harmonics of order , (2.12) turns into the classical model Blanco1997 (); Brunner2005 (); Seibold2012 ().

2.1 Angular bases

This moment approach strongly depends on the choice of the ansatz and the angular basis . In the following sections we will shortly derive the different angular bases for the models presented here. These bases will be generally collected in the basis-vector . If we need to further distinguish between the models, the corresponding symbols defined in the following sections will be used. If a result is independent on the choice of the basis we will just use as symbol.

2.1.1 Full moments

The full-moment basis of order consists of the tensorial powers of , i.e. (by abusing notation)

(2.14)

The corresponding tensorial moment of order is given by

(2.15)

which consists of the scalar moments

(2.16)

Note that an equivalent system can be obtained by using the corresponding real-valued spherical harmonics of order as angular basis Blanco1997 (); Brunner2005 (); Seibold2012 ().

2.1.2 Quarter moments

We follow the approach in Frank2006 () where general partial as well as the special case of quarter moments in two space-dimensions are treated. The main idea is not to integrate over the whole sphere but over subsets of it.

Definition 2.2.

For and we define its tensorial moment by

(2.17)

As for full moments the corresponding components of the tensorial moments are given by

(2.18)

In this paper, will be one of the following quarterspaces

or halfspaces

Note that for we have that .

For pure quarter moments, we will have for . The corresponding basis for quadrant is then given by

where should be understood as multiplication with every component. Consequently the complete set of basis functions is .

2.1.3 Mixed moments

As will be shown below the quarter-moment basis exhibits undesired properties when applied to the Laplace-Beltrami operator. This has inspired the works in Frank07 (); Schneider2014 () where so-called mixed moments were developed. The main problem of the quarter-moment basis is that the ansatz is not continuous in which is necessary for the solution of (2.1) in one space-dimension. There, mixed moments where constructed by starting from a half-moment ansatz and demanding continuity of the ansatz (2.10) with respect to this basis. There, it suffices to choose a full zeroth-order moment and half-moments for all higher order moments Schneider2014 ().

The construction of mixed moments in two dimensions works in the same spirit. We start with the general quarter-moment basis and demand continuity of the ansatz .

Having (2.10) in mind, we obtain multipliers for every quadrant . For example, for we have

At the poles of the sphere () we have , which implies that is continuous only if for all . Similarly it holds that along the quarter-space boundaries (i.e. , ) some of the multipliers have to be the same (exactly those whose component of does not vanish on the boundary, e.g. since at we have ).

Accordingly we obtain the moments

(2.19)
(2.20)
(2.21)

for , , , and .

In contrast to the one-dimensional setting one full moment, half moments for the basis functions contributing to either or direction, and quarter moments for the basis functions which contribute to both directions occur.

Note that in the fully three-dimensional setting the decomposition in -direction has to be taken into account, leading to octants instead of quadrants.

To embed this in the framework we choose our basis function of order as

(2.22)

Formulas to efficiently calculate the appearing integrals in case of a linear ansatz can be found in Appendix A.

2.2 Moments of the Laplace-Beltrami operator

All that remains to obtain a closed set of equations in (2.12) is to correctly evaluate . This can be done using the formal self-adjointness of the Laplace-Beltrami operator, using

Due to this, the calculation of these integrals does not a priori depend on the choice of the ansatz but is true for every .

2.2.1 Full moments -basis

Straight-forward calculations show that for , it holds that

(2.23)

Consequently, the corresponding moments are given by

2.2.2 Quarter moments -basis

We observe the following relations on the quadrant for , :

(2.24)
(2.25)

Similarly, the corresponding quantities in the other quarter-spaces can be obtained.

The moments of include the evaluation of the microscopic values and at the quarter-sphere boundaries. Note that, similar to the one-dimensional case, the mass-conservation property of the Laplace-Beltrami operator (i.e. ) is usually violated. This can be easily seen by summing up (2.24) over all quadrants, observing that at the angles , remains. Also note that this quantity is not rigorously defined for minimum-entropy closures due to the discontinuity in the ansatz .

2.2.3 Mixed moments -basis

We observe for the following:

The calculations for are equivalent to those of the quarter-moment basis.

Note that all these calculations are closure-independent. We therefore need to calculate , , , , and the semi-microscopic quantities for .

3 Realizability

In this section we will define precisely the concept of realizability. Furthermore we give necessary and sufficient conditions for first order quarter and mixed moment models.

Definition 3.1 (Realizability).

A moment vector is said to be realizable with respect to basis if there exists a non-negative distribution such that . The realizable set is defined as

(3.1)

is then called a realizing distribution.

Note that (2.9) is solvable if and only if .

A standard example for realizability are the realizability conditions of first order in the full-moment setting.

Example 3.2.

Let and . Then if and only if Ker76 ()

(3.2)

3.1 Quarter moments

In this section first order necessary and sufficient conditions for realizability of a quarter-moment vector will be given. This will be used later to derive the corresponding conditions for mixed moments.

Lemma 3.3.

For a vector of moments it is necessary and sufficient for the existence of a non-negative measure which realizes with respect to that

(3.3)

and the normalized first moment

(3.4)

satisfies .

Proof.

We only prove the statement for . The proof for the other quadrants works similarly.
Assume that in . Since we obtain

showing the necessity of (3.3). Since we obtain implying that . Together with (3.3) we have .

To show the sufficiency of (3.3) we give a realizing distribution function. A possible (but not necessarily unique) candidate is given by

(3.5)

where denotes the multi-dimensional Dirac-delta distribution111We assume for notational simplicity that has mass even on the boundary of integration.. If the distribution function is supported in . Thus

Therefore, is a realizing distribution for under the given assumptions. ∎

3.2 Mixed moments

With the knowledge of Lemma 3.3 we are able to provide the realizability conditions for mixed moments of order . ‚

Theorem 3.4 (First order necessary and sufficient conditions).

For a vector of moments

it is necessary and sufficient for the existence of a non-negative measure which realizes with respect to that

(3.6)

and .

Proof.

We start again with the necessity of (3.6):
Note that the vector of component-wise absolute values satisfies . Therefore for every unitvector , implying

This can be reformulated to

(3.7)

The left hand side will be extremal if is collinear to , i.e. . Thus (3.7) implies (3.6). The sign-constraints on the moments follow again from the signs of and in the corresponding halfspaces.

For sufficiency we give again a reproducing distribution. In the following we will make use that under (3.6) the reduced moment vector is located in . Similarly the quantities

(3.8)

will be in . Reformulating (3.6) in terms of the normalized moments collected in we see that

(3.9)

Also note that . We now embed the moments into the normalized mixed-moment space (which is a subset of ) in the following way. Define

(3.10a)
(3.10b)
(3.10c)
(3.10d)
which then fulfils .

Furthermore, any convex combination of those modified moment vectors along neighbouring quadrants (e.g. and ) lies on an isoline of , i.e. for we have

and analogously for the other half-space combinations. This is visualized in Figure 1.

Figure 1: Linear interpolation between two realizability boundaries in the projected three-space along an isoline of . The realizable set with respect to the quadrant is plotted in grey.

It can be shown that (and consequently ) can be written as a convex combination of the moments to . Indeed, defining

(3.11)

we see that

and finally

Since we see that can be realized (with respect to ) by a distribution function with support in , namely the quarter-moment distribution (realizing with respect to ) as given in equation (3.5) in Lemma 3.3.

Therefore, due to the linearity of the problem, a non-negative realizing distribution for is given by

(3.12)

4 Kershaw closures

A typical drawback of the minimum-entropy models defined by (2.10) is that the dual problem (2.9) cannot be solved analytically. The numerical solution, which has to be calculated at least once in every space-time cell, is challenging and expensive Hauck2010 (); Alldredge2014 ().

On the other hand, standard models may give physically irrelevant solutions since they do not ensure positivity of the underlying distribution function.

Due to this, Kershaw closures became recently a topic of increasing interest. They are constructed in such a way that they are automatically generated by a nonnegative distribution function . Therefore, the moment vector including the unknown highest moment is also realizable with respect to the basis of one order higher. Furthermore the flux function is chosen to be exact (i.e. where realizes the moment vector including the unknown highest moment) if is the isotropic moment. The last condition is also called the isotropic interpolation condition.

In one spatial dimension for a full-moment basis, the Kershaw closure is also exact on the realizability boundary, because there the realizing measure on the realizability boundary is unique. This property is no longer true in higher dimension (see e.g. Monreal ()) or for other models. However, it gives an idea of how to construct such a closure.

The name of the closure is dedicated to David Kershaw who first proposed such an idea in Ker76 (). For an introduction into Kershaw closures in one space dimension we refer to Schneider2014 (). The construction of Kershaw closures of first order as done in Schneider2014 () requires realizability information of second order. With this it is possible to linearly interpolate between different parts of the higher-order realizability boundaries (choosing the interpolation parameter in such a way that the isotropic moment is interpolated). The resulting model is then cheap to evaluate because it is analytically closed (in contrast to minimum-entropy models).

Unfortunately, we were not able to provide a closed second-order realizability theory for mixed or quarter moments yet which implies that we can’t identify the correct parts of this second-order realizability boundary for interpolation. However, it turns out that under some assumptions on the second-order realizability information from the half moments in one space dimension are sufficient to define the Kershaw closure for quarter moments in a similar fashion.

For mixed moments we choose a different approach to build the unknown second moment. Abusing the constructive procedure in Theorem 3.4 we are able to provide a closure for mixed moments by combining the quarter-moment closures accordingly.

4.1 Quarter moments

It was shown in Lev84 () that, by assuming that the distribution function is symmetric around a preferred direction, the second moment of the closure can be decomposed into