First-order continuous- and discontinuous-Galerkin moment models for a linear kinetic equation: realizability-preserving splitting scheme and numerical analysis

TL;DR

This paper develops a second-order, realizability-preserving numerical scheme for first-order continuous and discontinuous Galerkin moment models of linear kinetic equations, demonstrating competitive performance against classical full-moment models.

Contribution

It introduces a novel second-order splitting scheme that preserves realizability for first-order Galerkin moment models, with extensive numerical validation.

Findings

The new scheme maintains realizability in simulations.

The models can outperform classical full-moment models in test cases.

Numerical analysis confirms the effectiveness of the approach.

Abstract

We derive a second-order realizability-preserving scheme for moment models for linear kinetic equations. We apply this scheme to the first-order continuous and discontinuous models in slab and three-dimensional geometry derived in a previous paper as well as the classical full-moment $M_N$ models. We provide extensive numerical analysis as well as our code to show that the new class of models can compete or even outperform the full-moment models in reasonable test cases.

Tables

Table 1. Table 1 : Parameter choice for the different aspects of the simulation. Notation for the Newton algorithm as in [ 3 ] .

Newton algorithm
$k_0$	$k_{max}$	${\epsilon }_{\gamma }$	$ϵ$	$\chi$	$\xi$	$\tau$	$\left\{r_l\right\}$
$500$	$1000$	${10}^{-2}$	$2^{-52}$	$1/2$	${10}^{-3}$	${10}^{-9}$	$\{0,{10}^{-8},{10}^{-6},{10}^{-4},{10}^{-3}$
							${10}^{-2},0.05,0.1,0.5,1\}$
		Realizability limiter			Minima
		${\epsilon }_{ℛ}$	$\tilde{\epsilon }$		${\rho }_{vac}$	${\psi }_{\text{vac}}$
		${10}^{-11}$	${10}^{-11}$		${10}^{-8}$

Equations

$${\partial }_t\psi +\mathbf{\Omega }\cdot {\nabla }_{\mathbf{x}}\psi +{\sigma }_a\psi ={\sigma }_s\mathcal{C}\left(\psi \right)+Q,$$

$$\mathcal{C}\left(\psi \right)=\underset{{\mathcal{S}}^2}{\int }K(\mathbf{\Omega },{\mathbf{\Omega }}^{\prime })\psi (t,\mathbf{x},{\mathbf{\Omega }}^{\prime })\text{leavevmode}d{\mathbf{\Omega }}^{\prime }-\underset{{\mathcal{S}}^2}{\int }K({\mathbf{\Omega }}^{\prime },\mathbf{\Omega })\psi (t,\mathbf{x},\mathbf{\Omega })\text{leavevmode}d{\mathbf{\Omega }}^{\prime }.$$

$$\psi (0,\mathbf{x},\mathbf{\Omega })$$

$$\psi (t,\mathbf{x},\mathbf{\Omega })$$

$$\mathbf{\Omega }={\left(\sqrt{1-{\mu }^2}\cos (\phi ),\sqrt{1-{\mu }^2}\sin (\phi ),\mu \right)}^T=:{\left({\Omega }_x,{\Omega }_y,{\Omega }_z\right)}^T$$

$$\mathbf{u}=\underset{{\mathcal{S}}^2}{\int }\mathbf{b}\psi \text{leavevmode}d\mathbf{\Omega }=:⟨\mathbf{b}\psi ⟩$$

$${\partial }_t\mathbf{u}+{\nabla }_{\mathbf{x}}\cdot ⟨\mathbf{\Omega }\mathbf{b}\psi ⟩+{\sigma }_a\mathbf{u}={\sigma }_s⟨\mathbf{b}\mathcal{C}\left(\psi \right)⟩+⟨\mathbf{b}Q⟩.$$

$$\mathcal{H}(\psi )=⟨\eta (\psi )⟩\text{ under the moment constraints }⟨\mathbf{b}\psi ⟩=\mathbf{u}.$$

$$\mathbf{\alpha }(\mathbf{u}):=\underset{\tilde{\mathbf{\alpha }}\in {\mathbb{R}}^n}{\mathrm{argmin}}⟨{\eta }_{\ast }(\mathbf{b}\cdot \tilde{\mathbf{\alpha }})⟩-\mathbf{u}\cdot \tilde{\mathbf{\alpha }},$$

$${\hat{\psi }}_{\mathbf{u}}={\eta }_{\ast }^{\prime }\left(\mathbf{b}\cdot \mathbf{\alpha }(\mathbf{u})\right)$$

$$\eta (\psi )=\psi \log (\psi )-\psi$$

$${\partial }_t\mathbf{u}+{\partial }_x⟨{\Omega }_x\mathbf{b}{\hat{\psi }}_{\mathbf{u}}⟩+{\partial }_y⟨{\Omega }_y\mathbf{b}{\hat{\psi }}_{\mathbf{u}}⟩+{\partial }_z⟨{\Omega }_z\mathbf{b}{\hat{\psi }}_{\mathbf{u}}⟩+{\sigma }_a\mathbf{u}={\sigma }_s⟨\mathbf{b}\mathcal{C}\left({\hat{\psi }}_{\mathbf{u}}\right)⟩+⟨\mathbf{b}Q⟩.$$

$${\hat{\psi }}_{\mathbf{u}}=\mathbf{b}\cdot \mathbf{\alpha }(\mathbf{u})$$

$${\partial }_t\mathbf{u}+{\partial }_x{\mathbf{f}}_1\left(\mathbf{u}\right)+{\partial }_y{\mathbf{f}}_2\left(\mathbf{u}\right)+{\partial }_z{\mathbf{f}}_3\left(\mathbf{u}\right)=\mathbf{s}\left(\mathbf{u}\right),$$

$${\mathbf{f}}_1\left(\mathbf{u}\right)$$

$$\mathbf{s}\left(\mathbf{x},\mathbf{u}\right)$$

$${\partial }_t\psi +\mu {\partial }_z\psi +{\sigma }_a\psi ={\sigma }_s\mathcal{C}\left(\psi \right)+Q,t\in T,z\in X,\mu \in [-1,1].$$

$${\partial }_t\mathbf{u}+{\partial }_z⟨\mu \mathbf{b}{\hat{\psi }}_{\mathbf{u}}⟩+{\sigma }_a\mathbf{u}={\sigma }_s⟨\mathbf{b}\mathcal{C}\left({\hat{\psi }}_{\mathbf{u}}\right)⟩+⟨\mathbf{b}Q⟩.$$

$${\mathbf{f}}_N$$

$$h_l(\mu )=\text{mathbbm}{1}_{I_{l-1}}\frac{\mu -{\mu }_{l-1}}{{\mu }_l-{\mu }_{l-1}}+\text{mathbbm}{1}_{I_l}\frac{\mu -{\mu }_{l+1}}{{\mu }_l-{\mu }_{l+1}},$$

$${\mathbf{p}}_{I_i}=\text{mathbbm}{1}_{I_i}{\left(1,\mu \right)}^T,$$

$${\mathbf{f}}_N={\left(S_l^m(\mu ,\phi );l=0,\dots ,N,\text{leavevmode}m=-l,\dots ,l\right)}^T,$$

$${\mathbf{h}}_{n_v}=\left(h_0,\dots ,h_{n_v-1}\right),$$

$${\mathbf{p}}_n=\left({\mathbf{p}}_{\text{vbox}\text{leavevmode}\text{resizebox}4.51207pt0.0pt>\text{nointerlineskip}\text{set@color}K};\text{vbox}\text{leavevmode}\text{resizebox}9.20833pt0.0pt>\text{nointerlineskip}\text{set@color}K\in {\mathcal{T}}_h\right)=\left(\left(\text{mathbbm}{1}_{\text{vbox}\text{leavevmode}\text{resizebox}4.51207pt0.0pt>\text{nointerlineskip}\text{set@color}K},\text{mathbbm}{1}_{\text{vbox}\text{leavevmode}\text{resizebox}4.51207pt0.0pt>\text{nointerlineskip}\text{set@color}K}\mathbf{\Omega }\right);\text{vbox}\text{leavevmode}\text{resizebox}9.20833pt0.0pt>\text{nointerlineskip}\text{set@color}K\in {\mathcal{T}}_h\right),$$

$${\mathcal{R}}_{\mathbf{b}}=\left\{\mathbf{u}\in {\mathbb{R}}^n\text{leavevmode}:\text{leavevmode}\exists \psi (\mathbf{\Omega })\ge 0,\rho =⟨\psi ⟩>0,\text{ such that }\mathbf{u}=⟨\mathbf{b}\psi ⟩\right\}.$$

$$u_{2i}\ge 0\text{ and, for }{u}_{2i}>0,\frac{u_{2i+1}}{u_{2i}}\in I_i=[{\mu }_i,{\mu }_{i+1}]$$

$${\partial }_t\mathbf{u}$$

$${\partial }_t{\overline{\mathbf{u}}}_i=\frac{1}{\Delta z}\underset{z_{i-\frac{1}{2}}}{\overset{z_{i+\frac{1}{2}}}{\int }}\mathbf{s}\left(z,\mathbf{u}\right)\text{leavevmode}dz.$$

$${\overline{\mathbf{s}\left(z,\mathbf{u}\right)}}_i=\frac{1}{\Delta z}\underset{z_{i-\frac{1}{2}}}{\overset{z_{i+\frac{1}{2}}}{\int }}\mathbf{s}\left(z,\mathbf{u}\right)\text{leavevmode}dz=\mathbf{s}\left(z_i,{\overline{\mathbf{u}}}_i\right)+\mathcal{O}(\Delta z^2).$$

$${\partial }_t{\overline{\mathbf{u}}}_i=\mathbf{s}\left(z_i,{\overline{\mathbf{u}}}_i\right),$$

$$\mathcal{C}\left(\psi \right)=\frac{⟨\psi ⟩}{⟨1⟩}-\psi .$$

Full text

††journal: arXiv.org

1 Introduction

We consider moment closures, which are a type of (non-linear) Galerkin projection, in the context of kinetic transport equations. Here, moments are defined by taking velocity- or phase-space averages with respect to some (truncated) basis of the velocity space. Unfortunately, the truncation inevitably comes at the cost that information is required from the basis elements 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 $\text{P}_{N}$ closure [36], where $N$ is the degree of the highest-order moments in the model. The $\text{P}_{N}$ 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, $\text{P}_{N}$ solutions can contain negative values for the local densities of particles, rendering the solution physically meaningless. Entropy-based moment closures, typically denoted by $\text{M}_{N}$ models in the context of radiative transport [41, 18], have (for physically relevant entropies) all the properties one would desire in a moment method, namely positivity of the underlying kinetic density, hyperbolicity of the closed system of equations, and entropy dissipation [35]. These models are usually comparatively expensive as they require the numerical solution of an optimization problem at every point on the space-time grid. Practical interest in such models increased recently due to their inherent parallelizability [25]. While the cost of solving the local nonlinear problems in the $\text{M}_{N}$ model scales strongly with the number of moments $n$ (since one has to solve square problems of size $n$), the desired spectral convergence with respect to the moment order $N$ is only achieved for smooth test cases, which rarely occur in reality. This means that the gain in efficiency by increasing the order of approximation will become rather insignificant.

To increase the accuracy of the $\text{M}_{N}$ models while maintaining the lower cost for small moment order $N$, a partition of the velocity space while keeping the moment order fixed is useful, similar to some h-refinement for, e.g., finite element approximations [6]. We focus on the continuous and discontinuous piece-wise linear bases derived in [56], which aim to be a generalization of the special cases provided in [20, 19, 53, 44, 57] in slab geometry and the fully three-dimensional case.

Besides their inherent parallelizability, 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) are 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, even though desirable due to its consistency with the original kinetic distribution, can cause problems in numerical simulations. 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 [62]. Similar effects have been reported in the context of elastic flow [46]. This is exactly loss of realizability.

We propose a second-order realizability-preserving scheme, that is based on a splitting technique and analytic solutions of the stiff part, combined with a realizability-preserving reconstruction scheme. It turns out that this scheme is very effective for (medium) smooth and non-smooth test cases, which can also occur in practice. The realizability-preserving property is achieved using the realizability limiter proposed in [2, 54, 51, 15]. This limiter requires information about the set of realizable moments, which turns out to be very simple in the context of our first-order models [56]. Again, this additionally makes the implementation of such models faster (and easier) compared to standard $\text{M}_{N}$ models.

This paper is organized as follows. First, we shortly recall the transport equation, its moment approximations and the relevant results from [56] (Sections 2 and 3). Then, we propose our second-order realizability-preserving scheme and investigate all the required properties that it should fulfill (Section 4). In Section 5, we discuss some implementation details of our scheme. Finally, in Section 6, we give a comprehensive numerical investigation of our models and the $\text{M}_{N}$ models in slab geometry and three dimension, to show that our models can indeed compete with or even outperform the full-moment models.

2 Modeling

This section closely follows the corresponding part in [56]. We consider the linear transport equation

[TABLE]

The equation is supplemented with initial condition and Dirichlet boundary conditions:

[TABLE]

where $\mathbf{n}$ is the outward unit normal vector in $\mathbf{x}\in\partial X$.

Parameterizing $\boldsymbol{\Omega}$ in spherical coordinates we obtain

[TABLE]

where $\varphi\in[0,2\pi]$ is the azimuthal and $\mu\in[-1,1]$ the cosine of the polar angle.

Definition 2.1**.**

The vector of functions ${\mathbf{b}}:\mathcal{S}^{2}\to\mathbb{R}^{n}$ consisting of $n$ basis functions $b_{i}$, $l=0,\ldots n-1$ of maximal order $N$ (in $\boldsymbol{\Omega}$) is called an angular basis.

The so-called moments $\mathbf{u}=\left(u_{0},\ldots,u_{n-1}\right)^{T}$ of a given distribution function $\psi$ are then defined by

[TABLE]

where the integration is performed component-wise.

Furthermore, the quantity $\rho=\rho(\mathbf{u})\coloneqq\left<\psi\right>$ is called the local particle density. Additionally, $\mathbf{u}_{\text{iso}}=\left<{\mathbf{b}}\right>$ is called the isotropic moment.

Equations for $\mathbf{u}$ can then be obtained by multiplying (2.1) with ${\mathbf{b}}$ and integration over $\mathcal{S}^{2}$, resulting in

[TABLE]

Depending on the choice of ${\mathbf{b}}$ the terms $\left<\Omega_{x}{\mathbf{b}}\psi\right>$, $\left<\Omega_{y}{\mathbf{b}}\psi\right>$, $\left<\Omega_{z}{\mathbf{b}}\psi\right>$, and in some cases even $\left<{\mathbf{b}}\mathcal{C}\left(\psi\right)\right>$, cannot be given explicitly in terms of $\mathbf{u}$. Therefore an ansatz $\hat{\psi}$ has to be made for $\psi$ closing the unknown terms. This is called the moment-closure problem.

In this paper the ansatz density $\hat{\psi}$ is reconstructed from the moments $\mathbf{u}$ by minimizing the entropy-functional

[TABLE]

The kinetic entropy density $\eta\colon\mathbb{R}\to\mathbb{R}$ is strictly convex and twice continuously differentiable and the minimum is simply taken over all functions $\psi=\psi(\boldsymbol{\Omega})$ such that $\mathcal{H}(\psi)$ is well defined. This problem, which must be solved over the space-time mesh, is typically solved through its strictly convex finite-dimensional dual,

[TABLE]

where ${\eta}_{*}$ is the Legendre dual of $\eta$. The first-order necessary conditions for the multipliers $\boldsymbol{\alpha}(\mathbf{u})$ show that the solution to (2.5), if it exists, has the form

[TABLE]

This approach is called the minimum-entropy closure [35]. The resulting model has many desirable properties: symmetric hyperbolicity, bounded eigenvalues of the directional flux Jacobian and the direct existence of an entropy-entropy flux pair (compare [35, 52]).

The kinetic entropy density $\eta$ can be chosen according to the physics being modelled. As in [35, 25], Maxwell-Boltzmann entropy

[TABLE]

is used, thus ${\eta}_{*}(p)={\eta}_{*}^{\prime}(p)=\exp(p)$. This entropy is used for non-interacting particles as in an ideal gas.

Substituting $\psi$ in (2.4) with $\hat{\psi}_{\mathbf{u}}$ yields a closed system of equations for $\mathbf{u}$:

[TABLE]

Remark 2.2**.**

Note that using the entropy $\eta(\psi)=\frac{1}{2}\psi^{2}$ the linear ansatz

[TABLE]

is obtained, leading to standard continuous/discontinuous-Galerkin approaches. If the angular basis is chosen as spherical harmonics of order $N$, (2.9) turns into the classical $\text{P}_{N}$ model [9, 11, 58].

For convenience, we write (2.9) in the standard form of a non-linear hyperbolic system of partial differential equations:

[TABLE]

where

[TABLE]

For ease of visibility, we also consider our models in slab geometry, which is a projection of the sphere onto the $z$-axis [58]. The transport equation under consideration then has the form

[TABLE]

The shorthand notation $\left<\cdot\right>=\int\limits_{-1}^{1}\cdot\leavevmode\nobreak\ d\mu$ then denotes integration over $[-1,1]$ instead of $\mathcal{S}^{2}{}$. Finally, the moment system is given by

[TABLE]

3 Angular bases

We shortly recall the angular bases under consideration. For a detailed derivation and further information, we refer the reader to [56] .

3.1 Slab geometry

• •

Full-moment basis

[TABLE]

with the monomials or the Legendre polynomials $P_{l}^{0}$, $l=0,\ldots,N$.

• •

Piecewise-linear angular basis (hat functions, continuous-Galerkin ansatz) $\mathbf{h}_{n}=\left(h_{0},\ldots,h_{n-1}\right)^{T}$

[TABLE]

where $-1=\mu_{0}<\mu_{1}<\ldots<\mu_{n-2}<\mu_{n-1}=1$ are some angular “grid” points and $\mathbbm{1}_{I_{{i}}}(\mu)$ is the indicator function on the interval $I_{{i}}=[\mu_{i},\mu_{i+1}]$ (with $\mathbbm{1}_{I_{{-1}}}\equiv\mathbbm{1}_{I_{{n-1}}}\equiv 0$).

• •

Partial moments (discontinuous-Galerkin ansatz) $\mathbf{p}=\left(\mathbf{p}_{I_{0}},\ldots\mathbf{p}_{I_{k-1}}\right)$

[TABLE]

where $k$ is the number of intervals.

Definition 3.1**.**

The resulting linear (compare (2.10)) and nonlinear models (compare (2.8)) will be called $\text{P}_{N}$/$\text{M}_{N}$ (full moment basis), $\text{HFP}_{n}$/$\text{HFM}_{n}$ (hat functions basis) and $\text{PMP}_{n}$/$\text{PMM}_{n}$ (partial moment basis), respectively.

3.2 Angular bases in three dimensions

Albeit both approaches are not limited to this, we consider moments on spherical triangles. To that end, let $\mathcal{T}_{h}$ be a spherical triangulation of $\mathcal{S}^{2}$ and \mathchoice{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\displaystyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{9.20833pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\textstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{6.44583pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptstyle K}}}{\vbox{ \hbox{\leavevmode\resizebox{4.60416pt}{0.0pt}{{\char 62\relax}}} \nointerlineskip\hbox{\set@color\scriptscriptstyle K}}}\in\mathcal{T}_{h} be a spherical triangle. In this paper, the triangulation $\mathcal{P}$ will be obtained by dyadic refinement of the octants of the sphere $V=\mathcal{S}^{2}$, i.e. the coarsest triangulation contains the eight spherical triangles obtained by projecting the octahedron with vertices $\{{(\pm 1,0,0)}^{T}$, ${(0,\pm 1,0)}^{T}$, ${(0,0,\pm 1)}^{T}\}$ to the sphere and finer partitions are obtained by iteratively subdividing each spherical triangle into four new ones, adding vertices at the midpoints of the triangle edges. After $r$ refinements, we thus obtain $n_{v}(r)=4^{r+1}+2$ vertices and $n_{t}(r)=2\cdot 4^{r+1}$ spherical triangles.

The bases that we use are the following.

• •

Full-moment basis

[TABLE]

where ${S_{{l}}^{{m}}}$ are the real-valued spherical harmonics on the unit sphere [9, 58].

• •

Barycentric-coordinate basis functions

[TABLE]

where $n_{v}$ is the number of vertices of the triangulation and $h_{l}$ is the basis function defined using spherical barycentric coordinates on the $l$-th vertex as in [14, 30, 45].

• •

Partial moments on the unit sphere

[TABLE]

where $n=4\cdot\left|\mathcal{T}_{h}\right|$ is the number of moments.

Naming of the models will be analogous to the slab-geometry case, compare Definition 3.1.

3.3 Realizability

The minimum-entropy moment problem (2.5) has a solution if and only if the moment vector is realizable.

Definition 3.2**.**

The realizable set $\mathcal{R}_{{\mathbf{b}}}$ is

[TABLE]

If $\mathbf{u}\in\mathcal{R}_{{\mathbf{b}}}$, then $\mathbf{u}$ is called realizable. Any $\psi$ such that $\mathbf{u}=\left<{\mathbf{b}}\psi\right>$ is called a representing density.

Unfortunately, checking whether a moment vector is realizable is not trivial for general bases. However, for the piecewise linear moment models, the realizability conditions are particularly simple (see [56]).

Lemma 3.3**.**

For the hat function basis in one or three dimensions, $\mathbf{u}\in\operatorname{cl}\left(\mathcal{R}_{\mathbf{h}_{n}}\right)$ if and only if $u_{l}\geq 0$ for all $l=0,\ldots,n-1$.

Lemma 3.4**.**

For the partial moment basis in one dimension (slab geometry), $\mathbf{u}\in\operatorname{cl}\left(\mathcal{R}_{\mathbf{p}_{n}}\right)$ if and only if

[TABLE]

for all $i=0,\ldots,\frac{n}{2}-1$.

For more details on the realizability of the regarded models, see [56].

4 Second-order realizability-preserving splitting scheme

As already mentioned before, the minimum-entropy moment problem (2.5) has a solution if and only if the moment vector is realizable. This implies that it is mandatory to maintain realizability during the numerical simulation (since otherwise the flux function cannot be evaluated). Explicit high-order schemes have been developed in [2, 54]. Unfortunately, the physical parameters $\sigma_{s}$ and $\sigma_{a}$ directly influence the CFL condition, resulting in very small time steps for large scattering/absorption.

This can be overcome by using a first-order implicit-explicit time stepping scheme [51, 48, 49], treating the transport part explicit while implicitly solving the (time-)critical source term. Unfortunately, using higher-order IMEX schemes again results in a CFL condition of the same magnitude as for the fully explicit schemes.

We are interested in a second-order scheme for (2.11). This can be achieved by doing a Strang splitting for

[TABLE]

A second-order realizability preserving scheme will be obtained if both subsystems are solved with a (at least) second-order accurate and realizability-preserving scheme. For notational simplicity, we show the full scheme for one spatial dimension only. A generalization to structured meshes in higher dimensions is straightforward.

4.1 Source system

Let us start with the stiff part (4.1b) whose finite-volume form is given by

[TABLE]

Fortunately, using the midpoint rule, it holds that

[TABLE]

To obtain a second-order accurate solution of (4.2), it is thus sufficient to solve the system

[TABLE]

which is purely an ODE (in every cell). As mentioned above, we restrict ourselves to isotropic scattering, where we have $K(\mu,\mu^{\prime})=\frac{1}{\left<1\right>}=\frac{1}{2}$, i.e.

[TABLE]

The source term now becomes

[TABLE]

where $\mathbf{G}=\frac{\mathbf{u}_{\text{iso}}{(\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}})}^{T}}{\left<1\right>}$ is the matrix mapping the moment vector $\mathbf{u}$ to the isotropic moment vector with the same density $\mathbf{u}_{\text{iso}}(\mathbf{u})=\mathbf{u}_{\text{iso}}\cdot\frac{\rho(\mathbf{u})}{\left<1\right>}$. Here we assumed that there exists a vector $\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}}$ such that $\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}}\cdot{\mathbf{b}}\equiv 1$ (true for all regarded bases: $\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}}=(1,0,\ldots,0)^{T}$ for Legendre Polynomials, $\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}}=(\sqrt{4\pi},0,\ldots,0)^{T}$ for real spherical harmonics, $\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}}=(1,\ldots,1)^{T}$ for the hat functions basis, $\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}}=(1,0,1,0,\ldots)^{T}$ for the partial moments in slab geometry and $(1,0,0,0,1,0,0,0,\ldots)^{T}$ for the partial moment basis in three dimensions).

Since in this case, (4.4) is linear and the parameters $\sigma_{s}$, $\sigma_{a}$, $Q$ are time-independent, we solve it explicitly using matrix exponentials, trivially obtaining a realizable second-order accurate solution of (4.2).

Remark 4.1**.**

Note that in this specific situation, the solution of this sub-step does not depend on the moment closure used in the flux system.

Using the matrix exponential and the variation of constants formula, the solution to (4.4) is

[TABLE]

As $\mathbf{G}$ and $\mathbf{I}$ commute, we have

[TABLE]

It remains to compute the matrix exponential of ${\sigma_{s}}t\mathbf{G}$. As $\mathbf{G}\mathbf{u}_{\text{iso}}(\mathbf{u})=\mathbf{u}_{\text{iso}}(\mathbf{u})$, we have that $\mathbf{G}^{k}=\mathbf{G}$ for all $k\geq 1$. It follows

[TABLE]

Inserting (4.9) in (4.8), we get

[TABLE]

Plugging (4.10) into (4.7), we finally get

[TABLE]

If the source is also isotropic then $\mathbf{G}\left<{\mathbf{b}}Q\right>=\left<{\mathbf{b}}Q\right>=\left<{\mathbf{b}}\right>Q$ and (4.11) simplifies to

[TABLE]

which can easily be calculated without explicit calculation of $\mathbf{G}$ or any matrix operations.

4.2 Flux system

Let us now consider the non-stiff part (4.1a). This can be solved using standard realizability-preserving methods [54, 2, 52, 15], which will be summarized in the following.

The standard finite-volume scheme in semi-discrete form for (4.1a) looks like

[TABLE]

where $\mathbf{g}$ is a numerical flux function. The simplest example is the global Lax-Friedrichs flux

[TABLE]

The numerical viscosity constant $C$ is taken as the global estimate of the absolute value of the largest eigenvalue of the Jacobian $\mathbf{f}^{\prime}$. In our case, the viscosity constant can be set to $C=1$, because for the moment systems used here the largest eigenvalue is bounded in absolute value by one [2, 52, 42].

Another possible choice is the kinetic flux [54, 25, 22, 20]

[TABLE]

where $\left<\cdot\right>_{+}$ and $\left<\cdot\right>_{-}$ denote integration over the positive and negative half intervals $[-1,0]$ and $[0,1]$, respectively. The kinetic flux is less diffusive than the (global) Lax-Friedrichs flux and admits a more consistent implementation of kinetic boundary conditions (see [52] and Section 5.5). For this reason, we will use (4.15) in all our computations.

4.2.1 Polynomial reconstruction

The value $\mathbf{u}_{i+\frac{1}{2}}$ is the evaluation of a suitable linear reconstruction of $\mathbf{u}$ at the cell interface $z_{i+\frac{1}{2}}$. In one dimension, it can be obtained from a minmod reconstruction222Other second-order accurate reconstructions like WENO [28, 16] are also possible.

[TABLE]

where $\operatorname{minmod}\left(\cdot\right)$ is the minmod function

[TABLE]

applied componentwise. We then set $\mathbf{u}_{i+\frac{1}{2}}^{-}=\mathbf{u}_{i}(z_{i+\frac{1}{2}})$ and $\mathbf{u}_{i+\frac{1}{2}}^{+}=\mathbf{u}_{i+1}(z_{i+\frac{1}{2}})$.

To avoid spurious oscillations, the reconstruction has to be performed in characteristic variables. They are found by transforming the moment vector $\mathbf{u}$ using the matrix $\mathbf{V}_{i}$, whose columns hold the eigenvectors of the Jacobian $\mathbf{f}^{\prime}(\overline{\mathbf{u}}_{i})$ evaluated at the cell mean $\overline{\mathbf{u}}_{i}$. This leads to

[TABLE]

For details on the eigenvalue computation see Section 5.2. In several dimension, we perform a dimension-by-dimension reconstruction as in [60] using the minmod reconstruction in characteristic variables in each one-dimensional reconstruction step.

4.2.2 Realizability-preservation

While this already gives us a second-order scheme, we do not have the realizability-preserving property yet. To achieve this, we need to apply a realizability limiter, ensuring that $\mathbf{u}_{i}(z)$ is point-wise realizable at the interface nodes $z\in\{z_{i-\frac{1}{2}},z_{i+\frac{1}{2}}\}$. We follow the construction from [2].

We replace $\mathbf{u}_{i}$ with the limited version

[TABLE]

The limiter variable $\theta\in[0,1]$ dampens the reconstruction from unlimited ($\theta=0$) to first-order ($\theta=1$). Assuming that $\overline{\mathbf{u}}_{i}\in\mathcal{R}_{{\mathbf{b}}}$ is realizable, there exists at least one $\theta$ (namely $\theta=1$) such that $\mathbf{u}_{i}^{\theta}(z)\in\operatorname{cl}\left(\mathcal{R}_{{\mathbf{b}}}\right)$ for every $z$ in the set of quadrature nodes (where $\operatorname{cl}\left(\mathcal{R}_{{\mathbf{b}}}\right)$ is the closure of $\mathcal{R}_{{\mathbf{b}}}$). Since the realizable set is a convex cone, and by continuity, it is guaranteed that there exists a minimal $\theta$ satisfying this assumption. We are thus searching for the solution of the minimization problem

[TABLE]

In practice, given some interface node $z$, we search for the intersection of the line $\mathbf{u}_{i}^{\theta}(z)$ (wrt. $\theta$) with the boundary of the realizable set, check if the value is in $[0,1]$ and store it in the case that it is.

For the presented first-order moment models, the solution of the above limiter problem can often be computed explicitly (see Section 5.3 for more details).

If we discretize (4.13) with a second-order SSP Runge-Kutta (RK) scheme, e.g. Heun’s method or the general $s$ stage SSP ERK2 [23, 29], a realizability-preserving scheme is obtained under a CFL-like condition if reconstruction and limiting is performed in every stage of the RK method, see Lemma 4.5.

4.2.3 Solving the optimization problem

For the minimum-entropy models, in each stage of the time stepping scheme for (4.1a), we have to solve the optimization problem (2.5) once in each cell (to compute the Jacobians) and twice at each interface of the computational mesh (one optimization problem for the left and right reconstructed value at the interface, respectively). This usually accounts for the majority of computation time which makes it mandatory to pay special attention to the implementation of the optimization algorithm. In this section, we will focus on the stopping criteria for the optimization algorithm. For details on the implementation, see Section 5.1.

Recall that the objective function in the dual problem (2.6) is

[TABLE]

The gradient and Hessian of $p$ are given by

[TABLE]

and

[TABLE]

respectively. Note that ${\eta}_{*}^{\prime\prime}>0$ since we assumed that $\eta$ (and thus also ${\eta}_{*}$) is strictly convex, and remember that the basis functions contained in ${\mathbf{b}}$ are linearly independent. As a consequence, the Hessian $\mathbf{H}$ is symmetric positive definite.

To find a minimizer of $p$, we are searching for a root of the gradient $\mathbf{g}$ using Newton’s method. For simplicity, we will restrict ourselves to Maxwell-Boltzmann entropy (2.8) such that ${\eta}_{*}^{\prime}(\boldsymbol{\alpha}\cdot{\mathbf{b}})=\exp(\boldsymbol{\alpha}\cdot{\mathbf{b}})$. We use the Newton algorithm from [52, 2, 55] with some adaptions. Before entering the algorithm for the moment vector $\mathbf{u}$ we rescale it to

[TABLE]

such that $\rho(\boldsymbol{\phi})=1$. Let $\varrho(\boldsymbol{\alpha})=\left<\exp({\mathbf{b}}\cdot\boldsymbol{\alpha})\right>$ be the mapping $\boldsymbol{\alpha}\mapsto\rho(\mathbf{u}(\boldsymbol{\alpha}))$ which maps a set of multipliers $\boldsymbol{\alpha}$ to the density of its associated moment. If the optimization algorithm for $\boldsymbol{\phi}$ stops at an iterate $\boldsymbol{\beta}$, we return

[TABLE]

where $\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}}$ is the multiplier with the property that $\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}}\cdot{\mathbf{b}}\equiv 1$ (see Section 4.1). This ensures that the local particle density is preserved exactly:

[TABLE]

Given $\tau\in\mathbb{R}^{+}$, ${\varepsilon_{\gamma}}\in(0,1)$, we will stop the Newton iteration at iterate $\boldsymbol{\beta}$ if

[TABLE]

where $\widetilde{\boldsymbol{\alpha}}$ is obtained from $\boldsymbol{\beta}$ by (4.23) and, as always, $n$ is the number of moments.

In the following, we will explain the rationale behind these stopping criteria.

The first criterion guarantees that the gradient of the objective function is sufficiently small.

Lemma 4.2**.**

Let $\tau\in\mathbb{R}^{+}$. If (4.25a) is fulfilled, we have that $\left\lVert\mathbf{g}_{\mathbf{u}}(\widetilde{\boldsymbol{\alpha}})\right\rVert_{2}\leq\tau$.

Proof.

First note that, by its definition (4.20), the gradient can be written as

[TABLE]

where $\mathbf{u}(\widetilde{\boldsymbol{\alpha}})$ is the moment vector corresponding to the multipliers $\widetilde{\boldsymbol{\alpha}}$.

Let $\boldsymbol{\beta}^{\prime}\coloneqq\boldsymbol{\beta}-\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}}\log(\varrho(\boldsymbol{\beta}))$. Then it follows that

[TABLE]

Consequently, we have

[TABLE]

Moreover,

[TABLE]

where in the last step we used the explicit forms of $\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}}$ for the different bases (see Section 4.1). Since

[TABLE]

it follows from (4.28) and (4.25a) that

[TABLE]

for partial moments and hat functions, and

[TABLE]

for full moments, which directly gives

[TABLE]

respectively. Inserting these bounds in (4.27), we finally obtain

[TABLE]

for partial moments and hat functions, and similarly for full moments, removing $\sqrt{n}$ accordingly. ∎

The second criterion (4.25b) ensures that the ansatz density (2.7) corresponding to the multiplier $\widetilde{\boldsymbol{\alpha}}$ obtained from the Newton iteration is close enough to a representing density for the moments $\mathbf{u}$.

Lemma 4.3**.**

Let $\mathbf{u}\in\mathcal{R}^{+}_{{\mathbf{b}}}{}$ and let ${\varepsilon_{\gamma}}\in(0,1)$, $\widetilde{\boldsymbol{\alpha}}\in\mathbb{R}^{n}$ such that the second stopping criterion (4.25b) holds. Then there exists a representing distribution $\psi$ for $\mathbf{u}$, i.e. $\mathbf{u}=\left<{\mathbf{b}}\psi\right>$, such that

[TABLE]

Proof.

If (4.25b) is satisfied, there exists a positive distribution $\psi^{\varepsilon_{\gamma}}$ such that

[TABLE]

Then

[TABLE]

is a positive distribution representing $\mathbf{u}$ and satisfying (4.30). ∎

In Section 4.2.4 we will use Lemma 4.3 to show that the scheme is realizability-preserving although the optimization problems are only solved approximately.

Remark 4.4**.**

Note that

[TABLE]

is realizable for all ${\varepsilon_{\gamma}}>0$. Due to the openness of $\mathcal{R}^{+}_{{\mathbf{b}}}{}$, there exists a $\delta>0$ s.t. $\tilde{\mathbf{u}}\in\mathcal{R}^{+}_{{\mathbf{b}}}$ for all $\tilde{\mathbf{u}}$ with $\left\lVert\tilde{\mathbf{u}}-{\varepsilon_{\gamma}}\mathbf{u}\right\rVert_{2}<\delta$. Note further that

[TABLE]

so (4.25b) is fulfilled if $\left\lVert\mathbf{g}_{\mathbf{u}}(\widetilde{\boldsymbol{\alpha}})\right\rVert_{2}\leq\frac{\delta}{1-{\varepsilon_{\gamma}}}$, i.e. if our numerical solution to the approximation problem is close enough to the exact solution. For moments $\mathbf{u}$ that are very close to the realizable boundary (so $\delta$ is very small and in addition $\mathbf{H}$ may be very badly conditioned), we might not be able to achieve such an accuracy. In that case, we either use a regularized version of $\mathbf{u}$ (see (5.6)) or disable linear reconstruction (see Item 4 in Section 5). Choosing ${\varepsilon_{\gamma}}$ closer to $1$ makes it easier to fulfil (4.25b) at the expense of smaller time steps (see Lemma 4.5). In our computations, we used the value ${\varepsilon_{\gamma}}=0.01$ which worked well in practice.

4.2.4 Time-step restriction

Now we are able to put all the things together to show that one forward-Euler step of our scheme (4.13) is indeed realizability-preserving.

Lemma 4.5**.**

The finite volume scheme (4.13), using the kinetic flux (4.15) and the stopping criteria from Section 4.2.3, on a rectangular grid in $d$ dimensions preserves realizability under the CFL-like condition

[TABLE]

Proof.

Adapted from [52, Theorem 3.19]. As we are using time stepping schemes that consist of a convex combination of Euler forward steps, it is enough to show realizability preservation in a single Euler forward step. Consider the one-dimensional ($d=1$) case first. The update formula in one step is

[TABLE]

where $\psi_{j}^{(\kappa)}$ is an arbitrary representing density for $\mathbf{u}_{j}^{\left(\kappa\right)}$ and $\hat{\psi}_{j-\frac{1}{2}}^{+,\tau}$ is the ansatz distribution obtained from the approximate solution of the optimization problem. To preserve realizability, we have to ensure that $\psi_{j}^{(\kappa+1)}\geq 0$ for all $\mu\in[-1,1]$ and all cells $j$.

For $\mu>0$, after stripping away positive terms and using $\mu\leq 1$, we have

[TABLE]

where $\psi_{j+\frac{1}{2}}^{-}$ is the distribution from (4.30).

We have that

[TABLE]

where $\mathbf{u}_{i}^{\prime}$ is the (limited) slope on cell $i$. Thus we have

[TABLE]

and therefore a representing density for $\mathbf{u}_{j}^{\left(\kappa\right)}$ is $\frac{\psi_{j+\frac{1}{2}}^{-}+\psi_{j-\frac{1}{2}}^{+}}{2}$. Inserting this in (4.34) gives

[TABLE]

This is positive under the time step restriction

[TABLE]

The case $\mu\leq 0$ follows in a similar way.

In $d$ dimensions, the update formula changes to

[TABLE]

where $\mathbf{j}=(1,\ldots,d)^{T}$ is an index tuple. As in one dimension, we define the representing density $\psi_{\mathbf{j}}^{(\kappa+1)}$ and only regard the case $\boldsymbol{\Omega}_{l}>0\,\forall l$, the other cases follow similarly. After stripping away positive terms we are left with

[TABLE]

where $\Delta x=\min_{l}\Delta x_{l}$. We proceed as in one dimension and note that

[TABLE]

is a representing density for $\mathbf{u}_{\mathbf{j}}^{\left(\kappa\right)}$ for any partition of unity $\sum_{l=1}^{d}w_{l}=1$. Inserting this ansatz in (4.40) gives

[TABLE]

This is positive if

[TABLE]

So for given $\boldsymbol{\Omega}$ we have to find a partition of unity $\mathbf{w}$ such that the right-hand side of (4.43) is maximal, i.e., we want to find

[TABLE]

Obviously, the maximum is attained if $\frac{w_{l_{1}}}{\boldsymbol{\Omega}_{l_{1}}}=\frac{w_{l_{2}}}{\boldsymbol{\Omega}_{l_{2}}}$ for all $l_{1},l_{2}$ (otherwise we could increase the $w_{l}$ which belongs to the minimum and decrease the other ones a little). Taking the partition of unity property into account, we thus have to choose $w_{l}=\frac{\boldsymbol{\Omega}_{l}}{\left\lVert\boldsymbol{\Omega}\right\rVert_{1}}$. Inserting this in (4.44) gives

[TABLE]

Using this in (4.43), we end up with the time-step restriction

[TABLE]

∎

5 Implementation details

We implemented the whole scheme in the generic C**++** framework DUNE [8, 7], more specifically in the DUNE generic discretization toolbox dune-gdt [47] and the dune-xt-modules [39, 40].

As mentioned above, we advance the flux system in time using Heun’s method, which is a second-order strong-stability preserving Runge-Kutta scheme [23]. In each stage of the Runge-Kutta scheme, we perform the following steps:

Solve the optimization problem for the cell means $\overline{\mathbf{u}}_{i}$ in each grid cell (see Section 5.1). If regularization is needed, replace $\overline{\mathbf{u}}_{i}$ by its regularized version333This formally destroys the consistency of the scheme. However, since regularization rarely occurs (and only near the realizability boundary), this effect can usually be neglected in practice. (see Section 5.1.2). 2. 2.

Reconstruct the values at the cell interfaces using linear reconstruction in characteristic variables (see Section 4.2.1), using the solution of the optimization problems from step 1 to calculate the Jacobians. 3. 3.

Perform the realizability limiting (see Section 5.3). 4. 4.

Solve the optimization problem for all reconstructed values $\mathbf{u}_{i\pm\frac{1}{2}}$. If the solver fails for a reconstructed value, disable the linear reconstruction in that cell. 5. 5.

Evaluate the kinetic flux (4.15) and update the stage values according to (4.13).

In the following, we will give some details on the implementation of these steps.

5.1 Implementation of the minimum-entropy solver

Our solver for the optimization problem is based on the algorithm from [3]. It uses a Newton-type algorithm with Armijo line search, i.e. to find a minimizer of the objective function $p$ (see (4.19)), we are searching for a root of the gradient $\mathbf{g}$ (see (4.20)) in the Newton direction $\mathbf{d}(\boldsymbol{\alpha})$ which solves

[TABLE]

and then update the multipliers as

[TABLE]

where $\zeta_{k}$ is determined by a backtracking line search such that

[TABLE]

with $\xi\in(0,1)$.

We stop the optimization if the new iterate $\boldsymbol{\alpha}_{k+1}$ satisfies the stopping criteria (4.25), except that we use

[TABLE]

as the first stopping criterion instead of simply using (4.25a). This avoids numerical difficulties for moments with small density, where $\tau^{\prime}$ is in the order of $1$ and thus some iterates $\boldsymbol{\beta}$ with very large (in absolute values) entries might fulfil the stopping criterion by chance. Moreover, checking the second stopping criterion (4.25b) might be quite expensive (depending on the basis ${\mathbf{b}}$). We therefore check this criterion only if additionally

[TABLE]

holds. This criterion approximately ensures (4.30) (see [52, 3]) but, in general, is much easier to evaluate than (4.25b). For the $\text{HFM}_{n}$ models, however, checking realizability is just checking positivity, so in that case we do not need to check (5.5) first.

To improve the performance and stability of the algorithm, we use several additional techniques which we will detail in the following. The values of the algorithms’ parameters that we use in all computations are given in Table 1.

5.1.1 Adaptive change of basis

Though the Hessian $\mathbf{H}(\boldsymbol{\alpha})$ is positive definite and thus invertible, it may be very badly conditioned, especially for multipliers $\boldsymbol{\alpha}$ corresponding to moments $\mathbf{u}(\boldsymbol{\alpha})$ close to the boundary of the realizable set. Moreover, in general, the integral in the definition (4.21) of $\mathbf{H}$ can only be calculated approximately using a numerical quadrature (see Section 5.4). If the quadrature is not sufficiently accurate, the approximate Hessian may have a significantly worse condition or may even be numerically singular.

To improve this situation, a change of basis can be performed after each Newton iteration such that the Hessian at the current iterate becomes the unit matrix in the new basis [3]. We use this procedure in our implementation for all bases except for the hat function bases $\mathbf{h}_{n}$.

For the hat function basis, all matrices and vectors required in the optimization algorithm are sparse and exploiting this fact in the implementation greatly speeds up the computations. Including the change of basis destroys the sparsity and thus harms performance. In theory, this could be compensated by faster convergence and thus less iterations of the algorithm due to the condition improvements. Further, the algorithm with change of basis might use regularization less frequently and thus introduce less errors in the solution, as shown for the full moments in [3]. We thus compared the algorithm with and without change of basis in several test problems. The differences in the results were negligible in all tests cases and the version without change of basis was significantly faster. We thus do not use the adaptive change of basis for the hat functions.

The first-order partial moments have a similarly simple structure as the hat functions, so the adaptive change of basis might also not be needed for these models. However, the change of basis does not have a significant performance impact in this case as the support of each basis function is restricted to a single interval or spherical triangle, and thus all matrix operations can be performed on the $2\times 2$ or $4\times 4$ submatrices corresponding to an interval in 1d and a spherical triangle in 3d, respectively. Similar, quadrature evaluations can be performed for each interval or spherical triangle separately. For this reason, we include the adaptive change of basis in our optimization algorithm for the partial moments.

For details on the change of basis algorithm see [3]. Note that the stopping criteria (4.25) are computed in the original basis such that knowledge of the change of basis algorithm is not required to understand the presentation in this paper.

5.1.2 Regularization

For tests with strong absorption, the local particle density may become very small in parts of the domain. As a consequence, also the entries of the Hessian $\mathbf{H}$ become very small which may cause numerical problems. We thus choose a “vacuum” density $\psi_{\text{vac}}$ with corresponding local particle density $\rho_{vac}=\left<\psi_{\text{vac}}\right>$. We then enforce a minimum local particle density of $\rho_{vac}$ by replacing moments $\mathbf{u}$ with local particle density $\rho(\mathbf{u})<\rho_{vac}$ by the isotropic moment with vacuum density $\rho_{vac}$. Obviously, this approach leads to a violation of the conservation properties of the scheme. However, since we only replace moments with very small local particle densities by moments with slightly larger but still very small densities, the effect should be negligible in practice.

Additionally, if the optimizer fails for a moment vector $\mathbf{u}$ (for example, by reaching a maximum number of iterations or being unable to solve for the Newton direction) we use the isotropic-regularization technique from [3], i.e. we replace $\mathbf{u}$ by the regularized moment vector

[TABLE]

and retry the optimization. If the optimizer still fails, we increase $r$ until the optimizer succeeds, which is guaranteed at least for $r=1$ where $\mathbf{u}^{r}$ is isotropic. As the regularized moment vector $\mathbf{u}^{r}$ always has the same local particle density as the original moment vector $\mathbf{u}$, this technique does not violate the mass conservation of the scheme but it may potentially completely alter the solution. In practice, regularization is only used rarely and if it is used, a small regularization parameters is usually sufficient.

5.1.3 Caching

We use two types of caching. First, for each grid cell we store the moment vector from the last time step and the corresponding multiplier obtained by entropy minimization. In this way we do not have to solve the optimization problem again if the moment vector in that grid cell did not change during the last time step. In addition, we store the last few solutions of the minimization problem with corresponding input moment vectors per thread of execution, so if several grid cells contain the same values, we only have to perform the optimization once and then use the cached values. If we encounter a moment vector that can not be found in the caches, we take the moment vector that is closest to the input vector (in one-norm) and use the corresponding multiplier as an initial guess.

5.1.4 Linear solvers

In each iteration of the Newton scheme, we have to apply the inverse of a positive definite Hessian matrix. We assemble the matrices using the quadratures described in Section 5.4. Inversion is then done by computing a Cholesky factorization of the assembled matrix. For the full moment models, the Hessian matrices are dense, so we use the LAPACK [5] routine dpotrf to compute the factorization and then use dtrsv to actually invert the linear systems. For the $\text{PMM}_{n}$ models, the Hessian is block-diagonal (each block corresponds to one interval/triangle of the partition) such that we can perform the Cholesky decomposition independently for each block. For the $\text{HFM}_{n}$ models in one dimension, the Hessian matrices are tridiagonal, so we can use the specialized LAPACK algorithms dpttrf and dpttrs. In three dimensions, the $\text{HFM}_{n}$ Hessians are not tridiagonal anymore but still sparse, so we use the sparse SimplicialLDLT solver from the Eigen library [24].

5.2 Solving the eigenvalue problems

To avoid spurious oscillations, the reconstruction has to be performed in characteristic coordinates (see Section 4.2.1). For that reason, we have to compute the eigenvectors of the flux Jacobians

[TABLE]

where

[TABLE]

and

[TABLE]

(compare Section 4.2.3). Note that, in general, the Jacobian (5.7) is not symmetric. However, since $\mathbf{J}$ is symmetric and $\mathbf{H}$ symmetric positive definite (see Section 4.2.3), we can see that $\mathbf{f}^{\prime}(\overline{\mathbf{u}}_{i})$ is similar to a symmetric matrix, i.e.

[TABLE]

and thus has real eigenvalues. In our implementation, we explicitly compute the matrix representation and then use an eigensolver for non-symmetric matrices (LAPACK’s dgeevx) to obtain the eigen decomposition. Unfortunately, though the Jacobian is a real matrix with real eigen values and thus also admits a set of real eigenvectors, the standard solvers for non-symmetric eigen problems (apart from dgeevx, we also tested the EigenSolver of the Eigen library [24]) often return complex eigenvectors. We thus add a step to compute real eigenvectors from the complex ones. Note that if

[TABLE]

is a set of linearly independent complex eigenvectors to the same eigenvalue $\lambda$ for the Jacobian $\mathbf{f}^{\prime}$, where $\mathrm{i}$ is the imaginary unit and $\mathbf{y}_{m}\in\mathbb{R}^{n}$ are real vectors, then

[TABLE]

is a set of $2k$ real eigenvectors for $\mathbf{f}^{\prime}$. Moreover, there are at least $k$ linearly independent vectors in this set. To see that, assume the opposite, i.e. that any $k$ vectors from the set (5.11) are linearly dependent. Without loss of generality, we assume that every vector in (5.11) can be written as a linear combination of the first $k-1$ vectors, i.e.

[TABLE]

with coefficients $a_{m,r}\in\mathbb{R}$. Then, the $k$ vectors $\mathbf{z}_{m}$ can also be written as (complex) linear combinations of these $k-1$ real vectors

[TABLE]

and thus cannot be linearly independent.

Consequently, to get real eigenvectors for $\mathbf{f}^{\prime}$ from the complex ones computed by the eigensolver, we first sort the eigenvectors into sets belonging to the same eigenvalue and then perform a Gram-Schmidt process with the real and imaginary parts for each of these sets.

Remark 5.1**.**

While this procedure works reasonably well, a better approach would probably be to use the structure of the Jacobian and, instead of solving the non-symmetric eigenvalue problem

[TABLE]

solve the symmetric generalized eigenvalue problem [37, 12]

[TABLE]

and then get the eigenvectors as $\mathbf{z}=\mathbf{H}\tilde{\mathbf{z}}$. Since the matrices $\mathbf{J}$ and $\mathbf{H}$ are both symmetric and $\mathbf{H}$ is positive definite, we can use a specialized algorithm like LAPACK’s dsygv and directly obtain real eigenvectors. Moreover, we can take advantage of the sparsity of these matrices and, e.g., use a generalized eigenvalue algorithm aimed at band matrices like dsbgv. In contrast, explicit assembly of the term $\mathbf{J}\mathbf{H}^{-1}$ might destroy the structure and result in a dense matrix even if the two factor matrices are sparse.

For the partial-moment models, the eigen decomposition can be done block-wise on the $2\times 2$ or $4\times 4$ matrix blocks which reduces the cubic complexity of the eigen decomposition [43] to a linear complexity for increasing number of moments and thus greatly accelerates computations for large problems.

5.3 Realizability limiting

The linear reconstruction process in the finite volume scheme does not guarantee preservation of realizability. Thus, we need an additional limiting step (4.17) to ensure that we are able to solve the optimization problem (2.5) for the reconstructed values. Since, in general, we cannot solve the integrals occurring in the optimization problem analytically and have to approximate them by a numerical quadrature $\mathcal{Q}$, the admissible moment vectors are further restricted to the numerically realizable set ($\mathcal{Q}$-realizable set)

[TABLE]

where for an integrable function $f$, $\left<f\right>_{\mathcal{Q}}=\sum_{i=0}^{{n_{\mathcal{Q}}}-1}w_{i}f(\boldsymbol{\Omega}_{i})\approx\left<f\right>$ is the approximation of the corresponding integral $\left<\cdot\right>$ with the quadrature rule $\mathcal{Q}$. In general, the numerically realizable set is a strict subset of the analytically realizable set.

The numerically realizable set can be described as the convex hull of the basis function values at the quadrature nodes (see [3] for the Legendre basis, the proof can be easily adapted for the other bases)

[TABLE]

If $\rho$ depends linearly on $\mathbf{u}$ it follows

[TABLE]

We do not want the limited moments to be too close to to the boundary of the numerically realizable set as we are not able to solve the optimization problem (2.5) in that case (see [4]). Moving the limited value away from the boundary can be done in several ways. A simple but often sufficient method can be employed for all limiters presented in this section. We simply add a small parameter $\tilde{\varepsilon}$ to the final limiter variable $\theta$ [52]. A problem with this approach is that the connecting line between $\mathbf{u}$ and $\overline{\mathbf{u}}$ might be almost parallel to the boundary which possibly results in a limited moment that is still too close to the boundary. Another approach is to require a fixed distance $\varepsilon_{\mathcal{R}}$ to the boundary of $\mathcal{R}_{{\mathbf{b}}}^{\mathcal{Q}}$, i.e., to limit to the $(\mathcal{Q},\varepsilon_{\mathcal{R}})$-realizable set

[TABLE]

where $d(\cdot,\partial\mathcal{R}_{{\mathbf{b}}}^{\mathcal{Q}})$ is the Euclidian distance to $\partial\mathcal{R}_{{\mathbf{b}}}^{\mathcal{Q}}$. Limiting to this set is possible whenever $\overline{\mathbf{u}}$ is farther than $\varepsilon_{\mathcal{R}}$ away from the boundary. If $\overline{\mathbf{u}}$ is already in the $\varepsilon_{\mathcal{R}}$-range of the boundary, we disable reconstruction in that cell.

Unfortunately, checking whether a reconstructed value lies within the numerically realizable set is not trivial in general. In the following, we detail the limiting procedure for the different models. For the remainder of this section, let $\overline{\mathbf{u}}$ be the moment vector before reconstruction and $\mathbf{u}$ a reconstructed moment vector 444In one dimension, there are always two reconstructed values per grid cell (one at each interface), so each of the limiters described in the following is applied to both values and the larger $\theta$ is used in (4.17). In several dimensions, both reconstruction and limiting are performed independently for each coordinate direction.. Let further $\overline{u}_{l}$ and $u_{l}$ be the $l$-th component of $\overline{\mathbf{u}}$ and $\mathbf{u}$, respectively.

5.3.1 $\text{M}_{N}$ models

In [2, 52], the half space representation for the convex hull (5.18) was explicitly calculated before starting the time stepping, yielding

[TABLE]

where $n_{f}$ is the number of facets of the convex hull. During the time stepping, the intersection of the connecting line between $\mathbf{u}$ and $\overline{\mathbf{u}}$ and each facet can then be computed efficiently by solving

[TABLE]

(compare (4.17)) for the limiter variable. We thus obtain

[TABLE]

If $\rho\left(\mathbf{u}\right)\geq 1$ or $\rho\left(\overline{\mathbf{u}}\right)\geq 1$, the moments can simply be rescaled before applying the limiter [2, 52]. Alternatively, we can ignore the facet corresponding to the condition $\rho=\boldsymbol{\alpha}_{{\mathbf{b}}}^{\mathbbm{1}}\cdot\mathbf{u}\leq 1$ which gives the half-space description for the full numerically realizable set (5.17). In that case, no rescaling is necessary and we can easily ensure a minimum distance of $\varepsilon_{\mathcal{R}}$ to the realizable boundary by moving each facet in normal direction before calculating the intersections, resulting in

[TABLE]

instead of $d_{i}$ in (5.21). As for the other limiters, we disable reconstruction if $\overline{\mathbf{u}}$ does not lie within the $\varepsilon_{\mathcal{R}}$-realizable set, i.e. if

[TABLE]

However, explicit calculation of the convex hull is only viable for a relatively small number of moments (such that the convex hull has to be calculated in a low-dimensional space) or very sparse quadratures (such that the convex hull has to be calculated from a small number of points). For a larger number of moments and a reasonable fine quadrature, the construction of the convex hull takes excessively long. Moreover, even when the convex hull is available, the performance of this approach might be unacceptable as the number of facets grows rapidly with both the number of moments and the number of quadrature points [52]. We thus use this approach only for the partial moments (see Section 5.3.3) where we only have to calculate low-dimensional convex hulls.

For the $\text{M}_{N}$ models, as proposed in [52, Section 3.62], we instead utilize the quadrature description (5.16) of the numerically realizable set and limit by solving the linear program (LP)

[TABLE]

where $\tilde{w_{i}}$ should not be confused with $w_{i}$ but rather represents $w_{i}\psi(\boldsymbol{\Omega}_{i})$ for the sought representing distribution $\psi$. This approach removes the prohibitively costly explicit calculation of the convex hull. However, the runtime cost during the time stepping algorithm might be considerably higher as a linear program has to be solved for every reconstructed value.

Instead of using a single limiter variable $\theta$, principally, we can limit each component of $\mathbf{u}$ independently. This has been done, e.g., in the context of the Euler equations in [62]. However, if the limiting is naively performed in ordinary coordinates, spurious oscillations may occur, as the limiting in ordinary coordinates may actually increase the slope in one of the characteristic components. In our implementation, we thus limit each of the characteristic components independently. Let $\mathbf{V}$ be the matrix of eigenvectors of the Jacobian $\mathbf{f}^{\prime}(\overline{\mathbf{u}})$ and let $\mathbf{u}^{\text{c}}=\mathbf{V}^{-1}\mathbf{u}$, be the respective moment vectors in characteristic coordinates. Then we can find limiter variables $\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{n})$ for each characteristic component by solving the LP

[TABLE]

where the matrix $\tilde{\mathbf{V}}$ is defined as $\tilde{V}_{ij}=V_{ij}\left(u^{\text{c}}_{j}-\overline{u}_{j}^{\text{c}}\right)$ and the $i$-th column of $\mathbf{B}$ is ${\mathbf{b}}\left(\boldsymbol{\Omega}_{i}\right)$.

For the LP-based limiter, it is not clear how to ensure a fixed distance $\varepsilon_{\mathcal{R}}$ to the boundary. We thus use the method of adding a small parameter $\tilde{\varepsilon}$ to the final limiter variable by replacing $\theta\geq 0$ by $\theta\geq-\tilde{\varepsilon}$ in (5.24c) and using $\theta+\tilde{\varepsilon}$ instead of $\theta$ if it is in the interval $[-\tilde{\varepsilon},1-\tilde{\varepsilon}]$.

For checking realizability of $\mathbf{u}$, which is needed for the stopping criterion (4.25b) in the optimization algorithm, we solve the simpler LP

[TABLE]

Note that the limiters presented thus far are not restricted to the $\text{M}_{N}$ models but could be used for any basis ${\mathbf{b}}$. However, for the $\text{HFM}_{n}$ and $\text{PMM}_{n}$ models, due to the simpler realizability conditions, limiters that are both faster and easier to implement can be used.

5.3.2 $\text{HFM}_{n}$ models

For the hat functions the numerically realizable set and the realizable set agree for suitable quadratures [56]. As a consequence, we can use a limiter based on the analytical realizability conditions which only require component-wise positivity (see Section 3.3). We thus calculate the limiter variable $\theta$ (limiting to $\mathcal{R}_{\mathbf{h}_{n}}^{\mathcal{Q},\varepsilon_{\mathcal{R}}})$ as

[TABLE]

5.3.3 $\text{PMM}_{n}$ models

In one dimension, $\mathcal{R}_{\mathbf{p}_{n}}^{\mathcal{Q}}=\mathcal{R}_{\mathbf{p}_{n}}$ for suitable quadratures, so a limiter based on the analytical realizability conditions (3.3) can be used. We use a limiter variable $\theta_{i}$ per interval $I_{i}=[\mu_{i},\mu_{i+1}]$. If we require a distance of at least $\varepsilon_{\mathcal{R}}$ to the boundary, the realizability conditions (3.3) become

[TABLE]

If $\overline{\mathbf{u}}_{I_{i}}$ is already not $\varepsilon_{\mathcal{R}}$-realizable, we disable reconstruction for that interval. This results in the following limiter for one-dimensional partial moments

[TABLE]

where

[TABLE]

For the partial moment basis in three dimensions, the analytical and numerical realizable set differ. However, note that (5.18) holds separately for each spherical triangle (see [56, Lemma 5.13]), so we can explicitly calculate the half space representation (5.20) for each spherical triangle. Instead of calculating a convex hull in $n$ dimensions, as would be needed for the full moment models, we only have to calculate $\frac{n}{4}$ convex hulls in $4$ dimensions, which is considerably faster and usually finished within a few seconds in our implementation (remember that this calculation has to be done only once before the time stepping).

5.4 Implementation of quadrature rules

In one dimension, we use Gauss-Lobatto quadratures on each interval. These quadratures include the endpoints of the interval, which ensures that the numerically realizable set (see (5.16)) equals the analytically realizable set for hat functions and partial moments, see [56]. To choose a suitable quadrature order, we solved some of our numerical test cases for different quadrature orders and calculated the errors with respect to the reference solution (see Section S2 in the supplementary materials). As suggested by this analysis, for the first-order models, we use a quadrature of order 15 per interval of the partition $\mathcal{P}$. For the full moment $\text{M}_{N}$ models, we split the domain in the two intervals $[-1,0]$ and $[0,1]$ that are needed for calculation of the kinetic flux and use a quadrature of order $2N+40$ on each interval.

In three dimensions, for partial moments and hatfunctions, we are using Fekete quadratures [59] (from the TRIANGLE_FEKETE_RULE library [13]) mapped to the spherical triangles. The library provides seven Fekete quadratures of order $3$, $6$, $9$, $12$, $12$, $15$ and $18$, using $10$, $28$, $55$, $91$, $91$, $136$, $190$ quadrature points, respectively. The second rule of order $12$ contains some negative quadrature weights, so we do not use that quadrature. If we want to improve the approximation, we subdivide each spherical triangle in several smaller ones as in [10] and use the mapped Fekete quadrature on each subtriangle. The Fekete rules correspond to Gauss-Lobatto rules on the triangle edges and thus also include the vertices of each triangle [59], which simplifies the realizability preservation (see [56]). As suggested by our quadrature sensitivity analysis (see Section S2), we use a quadrature order of $15$ for the $\text{HFM}_{6}$ and $\text{PMM}_{32}$ models and a quadrature order of $9$ for the other $\text{HFM}_{n}$ and $\text{PMM}_{n}$ models. For the $\text{M}_{N}$ models, we use tensor-product quadrature rules of order $2N+8$ on the octants of the sphere.

For the hat function basis in one dimension, we alternatively explicitly calculate all integrals needed in the Newton algorithm using the analytical formulas and Taylor expansion at the numerical singularities of the analytical formulas (see [32, Appendix A1] for the explicit formulas). The Taylor expansion is performed in a neighborhood of radius $0.1$ around the singularity, up to vanishing remainder or a maximal order of $200$. This completely removes the need for quadrature rules. Note that the same approach could be used for the partial moments in one dimension. However, as the quadrature-based adaptive-change-of-basis algorithm is very efficient for partial moments, we did not implement the analytical formulas for partial moments.

For the hat function basis in three dimensions, integrals cannot be evaluated analytically anymore. We experimented with an approach where the integrals are expanded in a Taylor series representation (see Section S1 in the supplementary materials). However, it turned out that the Taylor series had to be computed up to a prohibitively high order. For this reason, we dismissed that approach and use the Fekete quadrature approach described above also for the hat functions in three dimensions.

5.5 Implementation of initial and boundary conditions

The initial values for the finite volume scheme are computed by integration of the kinetic equation’s initial values (2.1c).

[TABLE]

Since the initial values in our test cases are isotropic (see Section 6), i.e. $\psi_{t=0}(\mathbf{x},\boldsymbol{\Omega})=\psi_{t=0}(\mathbf{x})$, we only have to compute the velocity integral of the basis $\left<{\mathbf{b}}\right>$. For this integral, we use the same quadratures as in Section 5.4 to ensure that the result is numerically realizable. Except for the plane-source and point-source tests, the initial values are constant in each grid cell, so we use the midpoint quadrature to evaluate the spatial integral. For the plane-source test, we always use an even number of grid cells and distribute the Dirac delta at $\mathbf{x}=0$ into the two adjacent grid cells, i.e. the initial value in these grid cells is set to the constant $\psi_{t=0}(\mathbf{x})=\psi_{\text{vac}}+\frac{1}{2\Delta x}$. For the point-source test, we use a Gauss-Legendre tensor product quadrature of order 20 to evaluate the spatial integrals for the initial values.

Boundary conditions for the moment equations are implemented by replacing the ansatz function $\hat{\psi}_{\mathbf{u}_{j}}$ belonging to a grid cell ${j}$ outside of the computational domain (such cells often called “ghost cells”) by the boundary condition $\psi_{b}$ of the kinetic equation (2.1d) in the computation of the kinetic flux $\eqref{eq:kineticFlux}$.

6 Numerical results

We want to apply our moment models to several test cases in the one- and three-dimensional setting. We follow the FAIR guiding principles for scientific research [61] and publish the code that generates the following results in [34]. As already mentioned (see Section 5), the scheme was implemented in the DUNE generic discretization toolbox dune-gdt [47]. The computations were done on a varying number of nodes of a distributed memory computer cluster555Each node encloses two Intel Intel Skylake Xeon Gold 6140 CPUs ($2\times 18$ cores) and $92$GB RAM.. Communication between nodes was done via MPI (Message Passing Interface) [38]. On each node, we used a work-stealing task-based shared-memory parallelization which was implemented using Intel TBB [27].

6.1 Slab geometry (1D)

6.1.1 Plane source

In this test case an isotropic distribution with all mass concentrated in the middle of an infinite domain $z\in(-\infty,\infty)$ is defined as initial condition, i.e.

[TABLE]

where the small parameter $\psi_{\text{vac}}=0.5\cdot 10^{-8}$ is used to approximate a vacuum. In practice, a bounded domain must be used which is large enough that the boundary should have only negligible effects on the solution. For the final time $t_{f}=1$, the domain is set to $X=[-1.2,1.2]$ (recall that for all presented models the maximal speed of propagation is bounded in absolute value by one).

At the boundary the vacuum approximation

[TABLE]

is used again. Furthermore, the physical coefficients are set to $\sigma_{s}\equiv 1$, $\sigma_{a}\equiv 0$ and $Q\equiv 0$.

All solutions are computed with an even number of cells, so the initial Dirac delta lies on a cell boundary. Therefore it is approximated by splitting it into the cells immediately to the left and right. In all figures below, only positive $z$ are shown since the solutions are always symmetric around $z=0$.

Note that since the method of moments is indeed a type of spectral method, it can be expected that due to the non-smoothness of the initial condition the convergence towards the kinetic solution of this test case is slow (note that $\psi_{t=0}(\cdot,\mu)\notin L^{p}$ for any $p$). Nevertheless, it is an often-used benchmark revealing many properties of a moment model (see e.g. [22]).

Some exemplary solutions at the final time are shown in Figures 1, 2 and 3. Remember that the full-moment models are indexed by the basis order $N$ while the piecewise linear bases are indexed by the number of moments $n$. Since, in one dimension, a basis of the space of polynomials up to order $N$ has $n=N+1$ elements, we compare the $\text{P}_{N}$ and $\text{M}_{N}$ models to the piecewise linear models with $N+1$ moments.

As expected, there are strong oscillations about the reference solution (the analytical solution from [21]) for all tested models. With increasing number of moments, the number of peaks increases while their height decreases. The oscillations are considerably stronger for the linear models than for the corresponding minimum-entropy models. The $\text{M}_{N}$ models are closest to the reference solution, particularly for the low-order models.

This is also reflected in the convergence results, which can be found in Figure 4. Depicted are the $L^{1}$ and $L^{\infty}$ errors between the local particle densities $\rho(\mathbf{u})$ of the moment models and the analytic reference solution from [21] at the final time $t_{f}$. As expected, overall convergence is slow. The $\text{HFP}_{n}$, $\text{PMP}_{n}$ and $\text{P}_{N}$ models show very similar $L^{1}$ errors at all orders. With respect to $L^{\infty}$ norm, the $\text{PMP}_{n}$ models are slightly better than the $\text{P}_{N}$ models. As observed before for the $\text{P}_{N}$ models [50, 52], $\text{HFP}_{n}$ models with odd $n$ show a higher $L^{\infty}$ error than models with even $n$ due to a zero eigenvalue of the flux jacobian. A similar but much less pronounced behaviour can be seen for the $\text{PMP}_{n}$ models (for odd and even number of intervals $\frac{n}{2}$). For odd $n$, $\text{HFP}_{n}$ $L^{\infty}$ errors are close to the corresponding $\text{P}_{N}$ error. For even $n$ the $\text{HFP}_{n}$ models perform better than the $\text{P}_{N}$ models and similar to the $\text{PMP}_{n}$ model.

The entropy-based $\text{HFM}_{n}$, $\text{PMM}_{n}$ and $\text{M}_{N}$ models show lower errors than their linear counterparts both in $L^{1}$ and $L^{\infty}$ norm. The $\text{PMM}_{n}$ models perform slightly better than $\text{HFM}_{n}$ models of the same order. The difference between odd and even orders/number of intervals is much more pronounced than for the $\text{HFP}_{n}$ and $\text{PMP}_{n}$ models. The $\text{M}_{N}$ models give the lowest errors of all tested models. However, the errors are still high and the rate of convergence is equally bad for all models although the convergence rate seems to improve with higher orders, especially for the minimum-entropy-based models.

6.1.2 Source beam

The discontinuous version of the source-beam problem from [26] is presented. The spatial domain is $X=[0,3]$, and

[TABLE]

with initial and boundary conditions

[TABLE]

The final time is $t_{f}=2.5$ and the same vacuum approximation $\psi_{\text{vac}}$ as in the plane-source problem is used.

Convergence results can be found in Figure 5. The reference solution for this test case is computed from a direct finite difference discretization of the kinetic equation on a grid with $21000\times 14000$ elements. As expected due to the higher regularity of the test case666$\psi$ is “only” a discontinuous function compared the distributional setting in the plane-source test., the convergence for all tested models is much better than in the plane-source test. As a consequence, most of the models are relatively close to the reference solution which is why we do not show exemplary solutions plots. The piecewise linear models and the $\text{M}_{N}$ models converge with second and first order in $L^{1}$ and $L^{\infty}$ norm, respectively. The $\text{P}_{N}$ models show spectral convergence which is reflected in comparatively high errors and a slow rate of convergence for the lower-order models and then a rapid convergence for the high-order models. Despite the eventual high rate of convergence, for the maximal moment number $n=50$ regarded here the $\text{P}_{N}$ models are mostly outperformed by the other models. Again, the minimum-entropy-based models perform better than their linear counterparts although the rate of convergence is the same. In $L^{\infty}$ norm, the $\text{M}_{N}$, partial-moment and $\text{HFP}_{n}$ models again show a zig-zag pattern where, e.g., the $\text{HFP}_{n}$ models using an odd number of intervals (even number of moments $n$) perform better than the ones using an even interval number (odd $n$). For the partial-moment models, an even number of intervals gives better results. Notably, the $\text{HFM}_{n}$ do not show such an alternating behaviour.

6.2 Three dimensions

We now consider numerical results in three spatial dimensions with velocities on the unit sphere.

6.2.1 Point source

The point-source test is the three-dimensional analogue of the plane-source test (Section 6.1.1) in slab geometry. Due to the limitations in the resolution we use a smoothed version of the initial Dirac delta:

[TABLE]

where $\sigma=0.03$, $\psi_{\text{vac}}=\frac{10^{-8}}{4\pi}$. As before, we choose $\sigma_{s}\equiv 1$, $\sigma_{a}\equiv 0$ and $Q\equiv 0$. All models are calculated on $X=[-1,1]^{3}$ to the final time $t_{f}=0.75$. The grid size is chosen to be $\Delta x=\Delta y=\Delta z=0.02.$ The point-source test is well-suited to demonstrate symmetries (or symmetry breaks) appearing in the solution. We show some selected models in Figures 6 and 7, where we use the endcap geometry $[0,1]\times[-1,1]\times[-1,1]$ for the isosurfaces.

The reference solution itself is rotationally symmetric and can be computed analytically using the formulas by Ganapol [21]. It can be observed that the hat functions have a preferred directions of propagation, directly related to the position of the vertices in the spherical triangulation (e.g., the octahedron that defines the $\text{HFM}_{6}$ basis can be easily identified in Figure 6). Similar effects occur for the partial moments. However, the discontinuity of their basis is also reflected in the peaks along the boundaries of the spherical triangles (compare $\text{PMM}_{32}$). In contrast to this, the full-moment models preserve the rotational symmetry (compare $\text{M}_{3}$, where small irregularities in the solution arise due to the spherical quadrature rule) but adding more waves to the solution.

Finally, we show error plots for our models in Figure 8. The models show the expected slow convergence in the $L^{1}$-norm, similar to the plane-source test (Section 6.1.1). All first-order models show roughly order $\frac{1}{2}$, whereas the full-moment models have varying convergence rates. In the $L^{\infty}$-norm, the first-order models show order $1$ convergence in the beginning, which then slows down to order $\frac{1}{2}$ as well. The full-moment models are showing no (or very slow) convergence, which is the well-known Gibbs phenomenon.

Note that the $\text{PMM}_{n}$ models clearly outperform the other methods. In particular, they are as good as or even slightly better than the $\text{M}_{N}$ models whose calculation is significantly more expensive for the same degrees of freedom (see Section 6.3).

6.2.2 Checkerboard

The checkerboard test case is a lattice problem which is loosely based on a part of a reactor core [11]. We extend it in a straightforward manner to the three-dimensional case. The used geometry is shown in Figure 9. There are scattering (orange and green) and highly absorbing (black) regions. The parameters are chosen to be the following.

• •

Domain: $X=[0,7]^{3}$, subdivided into the three regimes

[TABLE]

• •

Final time: $t_{f}=3.2$,

• •

Parameters (compare Figure 9):

[TABLE]

• •

Initial condition: $\psi_{t=0}(\mathbf{x},\boldsymbol{\Omega})=\psi_{\text{vac}}:=\cfrac{10^{-8}}{4\pi}$ (approx. vacuum),

• •

Boundary conditions: $\psi_{b}(t,\mathbf{x},\boldsymbol{\Omega})=\psi_{\text{vac}}$.

Due to the discontinuous nature of the physical parameters, this test case is a challenging task for a numerical solver. We align our grid with the discontinuities of the parameters by using a multiple of $7$ (usually $70$) regularly spaced grid points in each direction.

Solution plots for selected models can be found in the supplementary materials (Figures S3.1–S3.3). The $\text{PMP}_{32}$ model and the $\text{P}_{N}$ models of order $N\in\{2,\ldots,9\}$ (only $N=9$ shown) have negative particle densities $\rho$. Surprisingly, the hat function basis $\text{HFP}_{n}$ has positive densities for all $n$ that we calculated.

We compare our models to a discrete ordinate implementation [31, 22] of second order. $L^{1}$- and $L^{\infty}$-errors of the local particle density $\rho$ can be found in Figure 10. All tested models converge with about first order both in $L^{1}$ and $L^{\infty}$ norm. Again, $\text{PMM}_{n}$ and $\text{HFM}_{n}$ models are comparable to the $\text{M}_{N}$ models.

6.2.3 Shadow

The shadow test case represents a particle stream that is partially blocked by an absorber, resulting in a shadowed region behind the absorber. The used geometry is shown in Figure 11.

The parameters are chosen to be the following.

• •

Domain: $X=[0,12]\times[0,4]\times[0,3]$

• •

Final time: $t_{f}=20$,

• •

Parameters:

[TABLE]

• •

Initial condition: $\psi_{t=0}(\mathbf{x},\boldsymbol{\Omega})=\psi_{\text{vac}}:=\cfrac{10^{-8}}{4\pi}$ (approx. vacuum),

• •

The isotropic particle stream with density $\rho=2$ enters the region via the boundary condition at $x=0$. At all other boundaries, vacuum boundary conditions are used.

[TABLE]

We show slices and isovalues of several models at the final time in the supplementary materials (Figures S.4–S.6). Again, several of the linear models (e.g. $\text{PMP}_{32}$ or $\text{P}_{22}$) show negative values (depicted in red). As in the previous test case, the partial moments perform very well. Compare, for example, the linear $\text{PMP}_{512}$ model to $\text{P}_{22}$ (which has roughly the same number of degrees of freedom). The partial-moment model approximates the reference much better, especially in the far field where also the small oscillations are captured accurately. A similar tendency can be observed for the hat function model $\text{HFP}_{258}$.

Both hat function and discontinuous minimum-entropy models show a good approximation of the absorber (compare $\text{PMM}_{32}$ and $\text{HFM}_{6}$). However, they are not able to provide a reasonable approximation in the far field. Further repartitioning of the sphere yields much better results in this case.

Investigating again the convergence towards the reference solution (see Figure 12), we see that the full-moment models are slightly superior in the beginning, but convergence slows down for higher $n$. Both $\text{HFM}_{n}$ as well as $\text{PMM}_{n}$ show a similar convergence behaviour. Again, taking running time into account, both models outperform the classical $\text{M}_{N}$ model in terms of efficiency (see Section 6.3). Note that we only computed $\text{M}_{N}$ models up to a order of $N=4$ since the higher-order models did not finish in the available computation time. The same is true for the $\text{HFM}_{n}$ models with $n>66$. Distributing the workload among more nodes of the distributed cluster did not improve computation times for these models, which is probably due to load-balancing issues. While the task-stealing algorithm (see above) ensures that the work is distributed evenly on each node, there is no load-balancing between nodes. In our implementation, the grid is distributed to the nodes using a decomposition of the domain in connected parts. Since the optimization problems are particularly challenging near the absorber (due to low particle densities and very anisotropic distributions), the node(s) containing the absorber often need much longer than the other nodes to solve the optimization problems. Here, an MPI-based load-balancing implementation may be required which would, however, significantly increase the communication overhead. As an alternative, in [33], we investigate a numerical scheme that avoids the non-linear optimization problems and thus does not show the same load-balancing issues.

In conclusion, moment models based on piecewise first-order continuous ($\text{HFP}_{n}$, $\text{HFM}_{n}$) or discontinuous ($\text{PMP}_{n}$, $\text{PMM}_{n}$) basis functions often approximate the true solution as good as or even better than the standard models ($\text{P}_{N}$, $\text{M}_{N}$) using polynomial bases. In contrast to the standard models, however, these models can be implemented very efficiently. This is especially true for the entropy-based models since the necessary realizability limiting can be based on the analytical realizability conditions (see Section 5.4). In addition, in case of the discontinuous models, all matrix operations can be performed on small matrix blocks which provides further performance advantages (also for the $\text{PMP}_{n}$ models, see Section 6.3).

6.3 Timings

Performance measurements can be found in Figure 13. The times were measured without parallelization. Displayed times are the minimum of three runs. Quadratures were chosen as described in Section 5.4. Measurements were done both for the first-order scheme without linear reconstruction and for the realizability-preserving second-order scheme (see Section 4). Profiling shows that the first-order scheme spends most of the time solving the optimization problems. For the second-order scheme, solving the eigen problems for the reconstruction in characteristic coordinates also has a large impact on the execution time. Here, computation times could probably be improved by using a generalized eigen solver which takes the structure of the Jacobians into account (see Section 5.2). Both the adaptive-change-of-basis scheme and the eigensolver have third-order complexity. We thus asymptotically expect third-order complexity in $n$ for both the first-order and the second-order scheme for the $\text{M}_{N}$ models. For the $\text{PMM}_{n}$ models, all operations (including the solution of the eigen problems) can be done block-wise, so we expect first-order complexity in that case. Regarding only the optimization problem, the same is true for $\text{HFM}_{n}$ models as all matrices involved are tridiagonal (in slab geometry) or very sparse (in three dimensions). However, the Jacobian of the flux function is not sparse in general for the $\text{HFM}_{n}$ models, so the eigen problems are currently solved with a standard third-order-complex eigensolver. These models would particularly benefit from an improved implementation of the eigensolver which exploits the fact that the Jacobians are products of two sparse symmetric matrices (see Section 5.2).

In slab geometry, we used a reduced version of the plane-source test case ($1000$ grid cells, final time $t_{f}=0.1$). For the $\text{HFM}_{n}$ models, two different implementations were tested: the backtracking Newton solver without change of basis (see Section 5.1.1) using quadratures to calculate the integrals and the same backtracking Newton solver where all needed integrals were solved using the analytical formulas (and Taylor expansion at the singularities, see Section 5.4). In three dimensions, we used a reduced version of the point-source test case ($I=10^{3}$ grid cells, single Runge-Kutta step).

As can be seen in Figure 13, for the first-order scheme, results are as expected except that the $\text{M}_{N}$ models show second-order complexity in slab geometry, probably because the matrices are relatively small here and thus the third-order matrix operations do not dominate the execution time. The $\text{HFM}_{n}$ implementation using analytic integrals is faster than the quadrature version but the difference is negligible in practice. Given that the implementation using analytic integrals is considerably more complex and that analytic realizability conditions can also be used for the quadrature-based version, we suggest to generally use a quadrature-based implementation also for the $\text{HFM}_{n}$ models.

For the second-order scheme, as expected, the $\text{PMM}_{n}$ models show first-order complexity both in slab geometry and in three dimension and thus are several orders of magnitudes faster than the other models. Curiously, the $\text{HFM}_{n}$ models are close to second-order complexity also in three dimensions. For even larger $n$, we expect the $\text{HFM}_{n}$ models to also increase with third-order due to the eigensolver but the results show that the $\text{HFM}_{n}$ models are much faster than the $\text{M}_{N}$ models for a long time.

Note that though the $\text{PMM}_{2}$ and $\text{M}_{2}$ models are equivalent, the measured times are different as the $\text{M}_{2}$ model uses the convex-hull based realizability limiter while the $\text{PMM}_{2}$ model uses the limiter based on the analytical realizability conditions to be consistent with the models with higher $n$.

7 Conclusions and outlook

We derived two classes of minimum-entropy moment models based on a continuous finite element basis as well as a discontinuous piece-wise linear basis. Both types of models are realizable, i.e., generated by a non-negative ansatz, such that important physical properties like positivity of mass are preserved. We demonstrated in various numerical tests in one and three dimensional geometry that those models are qualitatively competitive with the classical full-moment $\text{M}_{N}$ models of the same number of degrees of freedom if the solution of the kinetic equation only has a limited smoothness (since otherwise the $\text{M}_{N}$ models typically show spectral convergence). Additionally, the new models are much cheaper (with respect to running time) than the full moment models since the non-linear problems that have to be solved locally are much smaller and typically much easier to solve as well. Consequently, the new models are considerably more efficient in the sense that they reach the same approximation error with much less computation effort. In particular, the partial moments $\text{PMM}_{n}$ show a linear relation between wall time and number of moments, which is also true for the hat function basis if only a first-order scheme is used. If a higher-order discretization in space and time is required, the partial moments appear to be the model of choice. However, in some cases the discontinuity in the basis functions may lead to severe problems, for example when collision is modeled with the Laplace-Beltrami operator [53]. In such cases, $\text{HFM}_{n}$ might be favorable.

We provided a second-order realizability-preserving scheme by using a splitting technique and analytic solutions of the stiff part, combined with a realizability-preserving reconstruction scheme. Higher-order variants of this scheme can in principle be derived similarly, but we emphasize that we strictly focused on non-smooth problems, where the sense in applying schemes with (much) more than second order is questionable.

If the underlying problem admits more smoothness (especially in the velocity domain), higher-order moment models might be more appropriate to enhance the speed of convergence towards the kinetic solution. While this is rather straight-forward to define both in slab as well as three-dimensional geometry (partial moments can be constructed immediately while the hat-function basis can be extended to higher-order splines on the unit interval/unit sphere,respectively [1]), special care is required since the realizability conditions are needed in order to use our realizability-preserving scheme. Up to our knowledge, the corresponding realizability problems are only solved for partial moments (of arbitrary order) in slab geometry [17], while first approaches are given for second-order partial moments on quadrants/octants of the sphere [57].