# Cubature rules for weakly and fully compressible off-lattice Boltzmann
methods

Dominik Wilde  Andreas Krämer  Mario Bedrunka  Dirk Reith  Holger Foysi  Chair
of Fluid Mechanics, University of Siegen, Paul-Bonatz-Straße 9-11, 57076
Siegen-Weidenau, Germany  Institute of Technology, Resource and Energy-
efficient Engineering (TREE),
Bonn-Rhein-Sieg University of Applied Sciences, Grantham-Allee 20, 53757 Sankt
Augustin, Germany  Department of Mathematics and Computer Science, Freie
Universität Berlin, Arnimallee 6, 14195 Berlin, Germany  Fraunhofer Institute
for Algorithms and Scientific Computing (SCAI), Schloss Birlinghoven, 53754
Sankt Augustin, Germany

###### Abstract

Off-lattice Boltzmann methods increase the flexibility and applicability of
lattice Boltzmann methods by decoupling the discretizations of time, space,
and particle velocities. However, the velocity sets that are mostly used in
off-lattice Boltzmann simulations were originally tailored to on-lattice
Boltzmann methods. In this contribution, we show how the accuracy and
efficiency of weakly and fully compressible semi-Lagrangian off-lattice
Boltzmann simulations is increased by velocity sets derived from cubature
rules, i.e. multivariate quadratures, which have not been produced by the
Gauß-product rule. In particular, simulations of 2D shock-vortex interactions
indicate that the cubature-derived degree-nine D2Q19 velocity set is capable
to replace the Gauß-product rule-derived D2Q25. Likewise, the degree-five
velocity sets D3Q13 and D3Q21, as well as a degree-seven D3V27 velocity set
were successfully tested for 3D Taylor-Green vortex flows to challenge and
surpass the quality of the customary D3Q27 velocity set. In compressible 3D
Taylor-Green vortex flows with Mach numbers
$\mathrm{Ma}={\\{0.5;1.0;1.5;2.0\\}}$ on-lattice simulations with velocity
sets D3Q103 and D3V107 showed only limited stability, while the off-lattice
degree-nine D3Q45 velocity set accurately reproduced the kinetic energy
provided by literature.

###### keywords:

Lattice Boltzmann method , Cubature , Semi-Lagrangian , Gauss-Hermite
quadrature , Compressible

††journal: Journal of Computational Science

##  1 Introduction

The lattice Boltzmann method (LBM) [1, 2, 3, 4, 5] is an efficient approach
for the numerical simulation of fluids. Compared to other methods discretizing
the Boltzmann equation, such as discrete-velocity models [6, 7] or (unified)
gas-kinetic schemes [8, 9], the LBM exhibits three key properties. First, the
distribution functions are integrated along characteristics in time, leading
to the well-known 0.5-shift in the relaxation time [10]. Second, the
equilibrium distribution function is expressed as Hermite series, mostly of
degree two, but higher for compressible [11] and thermal cases [12]. Third,
the velocity space is integrated via quadrature rules, in particular by Gauß-
Hermite quadratures [12, 13]. By the latter, the unbounded velocity space of
the Boltzmann equation is expressed by only a limited number of weighted
particle velocities, tied together as a velocity set. The type of numerical
integration with the underlying weight function
$\mathrm{exp}\left(-|\mathbf{x}|^{2}\right)$ does not only appear in the LBM,
but for instance also in Kalman filters [14]. Therefore, the literature on
numerical integration is a good starting point for the search of well-suited
velocity sets [15, 16, 17]. More specifically, the literature specifies two
main approaches to derive the respective rules: Gauß-product rules and non-
product rules [15]. The literature denotes multivariate non-product rules also
as _cubature rules_ [16, 17].

The simpler and widely-used approach of the above-mentioned ones is the Gauß-
product rule, which is calculated by the outer product of a 1D quadrature. For
example, the D1Q3 velocity set is essentially the 1D Gauß-Hermite quadrature
constructed by the roots of the third order Hermite polynomial [18]. By
applying the Gauß-product rule, the D1Q3 is turned into a D2Q9 in two
dimensions, and by applying the same rule once more the D3Q27 velocity set is
derived [18]. Since the number of particle velocities exceeds the number of
the encoded physical moments, pruning the D3Q27 velocity set leads to the
D3Q15 and D3Q19 sets with different high-order truncation errors [19, 20, 21].
A modification of the regular D3Q27 is the recently introduced
crystallographic LBM [22], which expresses the particle distributions in a
body centered cubic arrangement of grid points. All above mentioned velocity
sets fit into a regular grid, which is another key asset of the well-
established on-lattice Boltzmann method, with the vector of the particle
velocities ending on one of the neighbouring grid points. For weakly
compressible flows, this collection of velocity sets in combination with a
second-order polynomial expansion of the Maxwell-Boltzmann equilibrium is
widely used, although not perfect in terms of Galilean invariance, due to
well-known errors in the stress tensor [23]. However, when a higher degree of
precision for the numerical integration is needed [24, 25], in particular for
compressible [26] or dilute flows [27], the resulting velocity sets become
unfeasible, for two reasons. Firstly, Gauß-product rules suffer from the
”curse of dimensionality”, i.e. the number of abscissae rapidly increases
especially for high spatial dimensions with high degree of precision.
Secondly, when derived from Hermite polynomials the degree-nine one-
dimensional velocity set D1Q5 holds abscissae, which do not fit on a regular
grid, so that an even larger velocity set D1Q7 with equidistant abscissae must
be used to obtain a sufficiently high integration order [28]. These
multivariate lattices with equidistant nodes are constructed by solving the
orthogonality relations [13, 24, 29].

Contrary to Gauß-product rules, the abscissae of cubature rules are not
derived by quadratures of lower dimension. They are rather found individually
for a certain pair of degree of precision and dimension. Due to this freedom
of velocity discretization, off-lattice Boltzmann methods are required to
apply them in the LBM. There is only a very limited number of publications,
which display simulations using off-lattice velocity sets [30, 31, 32].

The present manuscript therefore derives velocity sets from cubature rules and
explores them using the semi-Lagrangian lattice Boltzmann method (SLLBM,
[33]), with applications to both weakly and fully compressible flows. In
contrast to Eulerian time integration schemes like finite difference [34, 35],
finite volume [36, 37], discontinuous Galerkin LBM schemes[38, 39], the SLLBM
inherits the Lagrangian time integration along characteristics from the LBM
and recovers the off-lattice distribution function values by interpolation.
Recent works [40, 41, 42, 43, 44, 45, 46] provided evidence that the SLLBM is
a promising off-lattice Boltzmann method for the simulation of both weakly and
compressible flows. However, these works used either on-lattice velocity sets,
e.g. D2Q9, or D3Q27, or velocity sets that had been derived by the Gauß-
product rule, e.g. D2Q25 from 1D Gauß-Hermite quadrature. This gap is closed
by the present article. This work shows that cubature-based velocity sets
significantly enhance off-lattice Boltzmann simulations in terms of both
efficiency and accuracy.

The remainder of this manuscript is structured as follows. Quadrature and
cubature in the LBM are briefly recapitulated in section 2 with a list of all
investigated velocity sets in this article. Section 3 details the semi-
Lagrangian lattice Boltzmann model for both weakly and fully compressible
flows. The results section 4 studies three test cases: a compressible two-
dimensional shock-vortex interaction and the three-dimensional Taylor-Green
vortex both in the weakly and in the fully compressible regime, each of them
requiring different equilibria and velocity sets. Sections 5 and 6 provide
discussion and conclusion.

##  2 Cubature in lattice Boltzmann methods

Following [17], we consider the approximation of an integral of function
$\mathcal{F}$

$I(\mathcal{F})=\int_{\Omega}\omega(\mathbf{x})\mathcal{F}(\mathbf{x})d\mathbf{x},$
(1)

with $\Omega\subset\mathbb{R}^{D}$, weight function $\omega(\mathbf{x})\geq 0$
and dimensions $D\geq 2$ by quadrature or cubature rules of the form

$\mathcal{C}(\mathcal{F})=\sum_{i=0}^{Q-1}w_{i}\mathcal{F}(\mathbf{x}_{i}),$
(2)

with number of abscissae $Q$ and discrete weights $w_{i}$.

The degree of a quadrature is defined as the largest integer $d$ that yields
$I(\mathcal{F})=\mathcal{C}(\mathcal{F})$ for all monomials

$\prod_{i=0}^{D-1}x_{i}^{j_{i}}\quad\mathrm{with}\quad\sum_{i=0}^{D-1}j_{i}\leq
d.$ (3)

of degree $\leq d$.

These cubature formulas are applied to the distribution function of the force-
free BGK-Boltzmann equation

$\frac{\partial f}{\partial t}+\bm{\xi}\cdot\nabla
f=-\frac{1}{\lambda}\left(f-f^{\mathrm{eq}}\right),$ (4)

with (unshifted) relaxation time $\lambda$, particle distribution function
$f$, equilibrium distribution function $f^{\mathrm{eq}}$, and particle
velocities $\bm{\xi}$. To that end, the Hermite moments $\bm{a}^{(n)}$ are
gained in the form [12, 13]

$\bm{a}^{(n)}=\int_{\mathbb{R}^{D}}f\bm{\mathcal{H}}^{(n)}d\bm{\xi},$ (5)

where $\bm{\mathcal{H}}^{(n)}$ is a nth-order Hermite polynomial. The
distribution function $f$ is expressed as a finite Hermite series

$f\approx
f^{N}=\omega(\bm{\xi})\sum_{n=0}^{N}\frac{1}{n!}\bm{a}^{(n)}\bm{\mathcal{H}}^{(n)},$
(6)

with truncation order $N$ [13] with the result that moments of order $M$ are
exactly represented, if $M\leq N$. The velocity space is discretized by
replacing the integral of Eq. (5) by a weighted quadrature leading to

$\displaystyle\bm{a}^{(n)}$
$\displaystyle=\int_{\mathbb{R}^{D}}\omega(\bm{\xi})\frac{f}{\omega(\bm{\xi})}\bm{\mathcal{H}}^{(n)}d\bm{\xi}$
(7)
$\displaystyle=\sum_{i=0}^{Q-1}\frac{w_{i}f(\bm{\xi})\bm{\mathcal{H}}^{(n)}(\bm{\xi}_{i})}{\omega(\bm{\xi}_{i})}=\sum_{i=0}^{Q-1}{f_{i}\bm{\mathcal{H}}^{(n)}_{i}},$
(8)

with the replacements $f_{i}=w_{i}f(\bm{\xi})/\omega(\bm{\xi}_{i})$ and
$\bm{\mathcal{H}}^{(n)}(\bm{\xi}_{i})=\bm{\mathcal{H}}^{(n)}_{i}$. From [47]
it is known that moments of order $M$ can be exactly determined by quadrature
or cubature rules of degree $d\geq N~{}+~{}M$. In particular, for weakly
compressible flows the order of the moments is usually limited to $N=2$,
although $N=3$ cancels the cubic error when deriving the Navier-Stokes
equations by a Chapman-Enskog analysis [13]. In contrast, fully compressible
flows require even expansion order $N=4$. By overcoming the restriction that
abscissae $\bm{\xi}_{i}$ need to match a regular grid, cubature rules become
applicable, provided they approximate the weight function

$\omega(\bm{\xi})=\frac{1}{(2\pi)^{D/2}}e^{-|\bm{\xi}|^{2}/2}.$ (9)

In the literature, cubature rules are usually listed for the weight function
$\tilde{\omega}(x)=\mathrm{exp}\left(-|\mathbf{x}|^{2}\right)$. To obtain a
velocity set to approximate the moments in Eq. (7), the cubature’s abscissae
have to be scaled by $\sqrt{2}$. The lattice speed of sound $c_{s}$ of these
velocity sets is—unless otherwise specified—unity $c_{s}=1$, which also
implies the reference temperature $T_{0}=1$.

The following velocity sets were identified to be subject of an in-depth
examination:

  * 1.

First, the newly introduced degree-nine D2Q19 in comparison to the D2Q25
derived from Gauß-product rule for fully two-dimensional compressible flows.
The shape of D2Q19 is shown in Fig. 1

  * 2.

Second, in three dimensions velocity sets for both weakly and fully
compressible flows were chosen. The D3Q13 based on a icosahedron was already
used in a weakly compressible finite difference LBM with a triangular mesh
[32], whereas the D3Q21 derived from a dodecahedron has not been examined so
far. Both platonic solids, shown in Fig. 2, possess a high geometric isotropy,
i.e. the flow information encoded into the distribution function values is
transported by regularly spaced abscissae ending on the surface of a sphere.
The weights and abscissae of the D3Q13 and D3Q21 velocity sets are listed in
Tables 2 and 3, respectively. Both velocity sets were compared to the
customary D3Q27 on-lattice velocity set and in addition to a degree-seven
velocity set, which we call D3V27 [15]. This velocity set is not affected by
the cubic error in the stress tensor that troubles degree-five velocity sets.
Although in the same spirit, it is not identical to the velocity set presented
by Yudistiawan et al. [31].

  * 3.

Last, for three-dimensional compressible flows, a recently introduced D3Q45
velocity set was applied [48], based on a cubature rule by Konyaev [49, 50].
This off-lattice velocity set was compared in an off-lattice Boltzmann
simulation to the state of the art on-lattice counterparts D3Q103 and D3V107
applied to compressible on-lattice Boltzmann simulations.

Figure 1: Shape of the D2Q19 velocity set.

(a)  (b)

Figure 2: Shape of icosahedron a) and dodecahedron b), which are the bases for
the velocity sets D3Q13 and D3Q21 used in this work.

Table 1 designates the velocity sets, which were used for the simulations in
this article. The additional notation $E_{D,d}^{Q}$ follows the literature on
multivariate quadratures integrating the weight function
$\tilde{\omega}(x)=\mathrm{exp}\left(-|\mathbf{x}|^{2}\right)$, listing
dimension $D$, degree $d$ and number of abscissae $Q$. To identify the
suitability of the introduced velocity sets beforehand, the symmetry
conditions in B were successfully tested [51]. For example, a degree-five
velocity set only reproduces the symmetry condition up to fifth order, whereas
degree-nine velocity sets will also correctly reproduce the symmetry
conditions up to ninth order.

Table 1: Quadrature and cubature rules utilized in the present article. Notation $E_{D,n}^{Q}$ according to Stroud [15] with number of abscissae $Q$, dimension $D$, and degree of precision $d$. Name |  $E_{D,d}^{Q}$ |  Remarks | Sources
---|---|---|---
D2Q19 |  $E_{2,9}^{19}$ |  Cubature rule by Haegemans and Piessens |  [52], C
D2Q25 |  $E_{2,9}^{25}$ |  Gauß-product rule from 1D degree-nine Gauß-Hermite quadrature | [28]
D3Q13 |  $E_{3,5}^{13}$ |  Derived by the vertices of an icosahedron |  [15, 32], Table 2
D3Q21 |  $E_{3,5}^{21}$ |  Derived by the vertices of a dodecahedron |  [15], Table 3
D3Q27 |  $E_{3,5}^{27}$ |  Doubly applied Gauß-product rule from 1D Gauß-Hermite quadrature | [53]
D3V27 |  $E_{3,7}^{27}$ |  Cubature rule by Stroud and Secrest. | [15]
D3Q45 |  $E_{3,9}^{45}$ |  Cubature rule by Konyaev. |  [49, 50], E
D3Q103 |  $E_{3,9}^{103}$ |  On-lattice velocity set |  [25]
D3V107 |  $E_{3,9}^{107}$ |  On-lattice velocity set | [30]
Table 2: Abscissae $\xi_{i}$ and weights $w_{i}$ for the D3Q13 velocity set, based on an icosahedron [15, 32, 4] with $r^{2}=(5+\sqrt{5})/2$ and $s^{2}=(5-\sqrt{5})/2$. The lattice speed of sound is $c_{s}=1$. $i$ | $w_{i}$ | $\xi_{i}$
---|---|---
$0$ | $2/5$ | $(0,0,0)$
$1,\ldots,4$ | $1/20$ | $(0,\pm r,\pm s)$
$5,\ldots,8$ | $1/20$ | $(\pm s,0,\pm r)$
$9,\ldots,12$ | $1/20$ | $(\pm r,\pm s,0)$
Table 3: Abscissae $\xi_{i}$ and weights $w_{i}$ for the D3Q21 velocity set, based on a dodecahedron [15] with $\phi=(1+\sqrt{5})/2$. The lattice speed of sound $c_{s}$ is $c_{s}=\sqrt{3/5}$. $i$ | $w_{i}$ | $\xi_{i}$
---|---|---
$0$ | $2/5$ | $(0,0,0)$
$1,\ldots,8$ | $3/100$ | $(\pm 1,\pm 1,\pm 1)$
$9,\ldots,12$ | $3/100$ | $(0,\pm\phi,\pm 1/\phi)$
$13,\ldots,16$ | $3/100$ | $(\pm 1/\phi,0,\pm\phi)$
$17,\ldots,20$ | $3/100$ | $(\pm\phi,\pm 1/\phi,0)$

##  3 Model description

The lattice Boltzmann equation reads

$f_{i}(\mathbf{x}+\bm{\xi}_{i}\delta_{t},t+\delta_{t})=f_{i}(\mathbf{x},t)-\frac{1}{\tau}\left(f_{i}(\mathbf{x},t)-f_{i}^{\mathrm{eq}}(\mathbf{x},t)\right)$
(10)

with relaxation time $\tau=\nu/(c_{s}^{2}\delta_{t})+0.5$ depending on
kinematic viscosity $\nu$, and time step size $\delta_{t}$.

This work applies the discrete equilibrium distribution function based on an
expansion in terms of Hermite polynomials

$f_{i}^{\mathrm{eq},N}(\mathbf{x},t)=w_{i}\sum_{n=0}^{N}\frac{1}{n!}\bm{a}^{(n)}_{eq}(\mathbf{x},t):\bm{\mathcal{H}}_{i}^{(n)},$
(11)

with $"\\!\\!:\\!\\!"$ denoting full contraction, $w_{i}$ the discrete weights
and $N$ the expansion order. Both the moments $\bm{a}^{(n)}_{eq}$ and the
Hermite tensors $\bm{\mathcal{H}}_{i}^{(n)}$ are listed in the Appendix. In
the weakly compressible case the local temperature $\theta=T/T_{0}$ is set to
$\theta=1$ with vanishing terms for second and higher order moments. With the
help of Eq. 7 the moments density and momentum are obtained

$\displaystyle\rho$ $\displaystyle=\sum_{i=0}^{Q-1}f_{i}$ (12)
$\displaystyle\rho\mathbf{u}$
$\displaystyle=\sum_{i=0}^{Q-1}\bm{\xi}_{i}f_{i}.$ (13)

Compared to the standard LBM, the SLLBM replaces the node-to-node streaming
step by a cell-wise interpolation procedure using interpolation polynomials
$\psi$ of interpolation order $p$ to obtain the departure points, whose
locations are dictated by the respective reversed particle velocities
$f_{i}(\mathbf{x}-\delta_{t}\bm{\xi}_{i})$ , i.e. for all $\mathcal{M}$ points
$\mathbf{x}$ in each cell $\Xi$:

$f_{i}(\mathbf{x},t)=\sum_{j=1}^{\mathcal{M}}\hat{f}_{i\Xi j}(t)\psi_{\Xi
j}(\mathbf{x}),$ (14)

with $\hat{f}_{i\Xi j}$ denoting the distribution function value at the
support points in each cell.

Gauß-Lobatto-Chebyshev points were used for the distribution of support points
in the cells [44, 33].

For the compressible test cases in Sections 4.1 and 4.3 we employed the
recently introduced compressible SLLBM [44] with the following extensions. A
second distribution function $g$ is introduced to enable the calculation of
flows with variable heat capacity ratio $\gamma$ [47], following the same
collide-and-stream algorithm as applied to the distribution function $f$.

$g_{i}(\mathbf{x}+\bm{\xi}_{i}\delta_{t},t+\delta_{t})=g_{i}(\mathbf{x},t)-\frac{1}{\tau}\left(g_{i}(\mathbf{x},t)-g_{i}^{\mathrm{eq}}(\mathbf{x},t)\right)$
(15)

The equilibrium of g is determined by the relation

$g^{\mathrm{eq}}_{i}=\theta(2C_{v}-D)f^{\mathrm{eq}}_{i},$ (16)

with heat capacity at constant volume $C_{v}$. To complement Eqs. (12) and
(13), the local temperature $\theta$ is determined by

$2\rho
C_{v}\theta=\sum_{i=0}^{Q-1}\left(|\bm{\zeta}_{i}|^{2}\frac{f_{i}}{T_{0}}\\!+g_{i}\right),$
(17)

with peculiar particle velocity $\bm{\zeta}_{i}=\bm{\xi}_{i}-\mathbf{u}$.
The local speed of sound $c_{s}^{*}$ consequently depends on the local
relative temperature $\theta$ and on the heat capacity ratio $\gamma$, via

$c_{s}^{*}=c_{s}\sqrt{\theta\gamma}.$ (18)

The relaxation time $\tau$ in the compressible case is dependent on the local
pressure $P$, i.e. $\tau=\mu/(c_{s}^{2}\delta_{t}P)+0.5$ with dynamic
viscosity $\mu$. In addition, a quasi-equilibrium approach was used to vary
the Prandtl number $\mathrm{Pr}$. More details about the compressible SLLBM
solver are listed in [44]. The NATriuM solver was used for the off-lattice
Boltzman simulations [45], which is based on the finite element library
deal.ii [54]. The on-lattice Boltzmann simulations were performed using the
lettuce software [55, 56], being based on the GPU-accelerated machine learning
toolkit PyTorch [57].

##  4 Results

The present work considers three test cases to evaluate the introduced
velocity sets: the two-dimensional shock-vortex interaction, comparing the
degree-nine velocity sets D2Q19 and D2Q25. This flow challenges the solver by
acoustic emissions by the vortex downstream the shock. Next, the weakly
compressible three-dimensional Taylor-Green vortex [58] was simulated by
various 3D velocity sets, e.g. D3Q13, D3Q21, D3V27 and D3Q27 to explore the
capability of the velocity sets to deal with underresolved flows. Finally,
simulations of fully compressible 3D Taylor-Green vortex flows were explored
with Mach numbers 0.5, 1.0, 1.5, and 2.0. This test case exhibits shocklets
and viscous effects and features a standardized initialization. The on-lattice
simulations were performed using the D3Q103 and D3V107 velocity sets, while
the off-lattice simulations were run by the D3Q45 velocity set.

###  4.1 Compressible 2D shock vortex interaction

The compressible two-dimensional shock-vortex interaction was tested as a
first test case to compare the two velocity sets D2Q19 and D2Q25. This
benchmark was intensively studied by Inoue and Hattori [59], with the
following setup. Two regions with $\mathrm{Ma}_{a}=1.2$ and $\mathrm{Ma}_{b}$
determined by Rankine-Hugoniot conditions are divided by a steady shock. The
Reynolds number is defined as $\mathrm{Re}=c_{s,\infty}^{*}R/\nu_{\infty}$
with subscript $\infty$ denoting the inflow conditions. The flow field of the
vortex is given by

$u_{\theta}(r)=\sqrt{\gamma T}\mathrm{Ma_{v}}\,r\,\mathrm{exp}((1-r^{2})/2),$
(19)

where $\mathrm{Ma}_{v}$ denotes the Mach number of the vortex. The initial
pressure and density field are

$P(r)=\frac{1}{\gamma}\left(1-\frac{\gamma-1}{2}\mathrm{Ma_{v}}\mathrm{exp}(1-r^{2})\right)^{\gamma/(\gamma-1)}$
(20)

and

$\rho(r)=\left(1-\frac{\gamma-1}{2}\mathrm{Ma_{v}}\mathrm{exp}(1-r^{2})\right)^{1/(\gamma-1)}.$
(21)

The resolution was $256\times 256$ cells with polynomial order $p=4$, ie.
$1024\times 1024$ grid points for the physical domain of $60R\times 24R$.
Initially the centre of the vortex was located at $x_{\mathrm{v}}=32$, while
the shock was located at $x_{\mathrm{shock}}=30$. Likewise as in [44] we used
a stretched grid to spatially resolve the shock region. The time step size was
$\delta_{t}=0.002$ with characteristic time $t^{\prime}=R/c_{s,\infty}^{*}$.
The Prandtl number was $\mathrm{Pr}=0.75$ in all cases. Table 4 designates the
most important parameters of the simulations. Fig. 3 depicts the density
contours at $t=8$ for test case a) with $\mathrm{Ma_{v}}=0.5$,
$\mathrm{Re}=400$ and provides evidence that there is no significant
difference between the simulations run by the D2Q19 and D2Q25 velocity sets.

Table 4: Simulation setup of the compressible 2D shock vortex interaction. Test case |  | a) | b)
---|---|---|---
Reynolds number | $\mathrm{Re}$ | 400 | 800
Vortex Mach number | $\mathrm{Ma_{v}}$ | 0.5 | 0.25
Advection Mach number | $\mathrm{Ma_{a}}$ | 1.2 | 1.2
Relaxation time at $\theta=1.0$ | $\tau$ | 3.82 | 2.16
Minimal spatial resolution | $\delta_{x}$ | 0.003 | 0.003
Temporal resolution | $\delta_{t}$ | 0.002 | 0.002
Interpolation polynomial order | $p$ | 4 | 4
Truncation order of equilibrium | $N$ | 4 | 4
Figure 3: Density contours of the shock-vortex interaction with $Ma_{v}=0.5$
and $Re=400$ in the range $\rho\in[0.92,1.55]$ in 119 steps with velocity sets
D2Q19 (top) and D2Q25 (below). No significant difference is visible despite
the D2Q19’s 24 percent reduction in computational cost.

The evaluation of the sound pressure for test case b) at
$\mathrm{Ma_{v}}=0.25$ and $\mathrm{Re}=800$ confirms this observation, being
shown in Fig. 4. The sound pressure $\Delta P=(P-P_{B})/P_{B}$, with subscript
$B$ denoting the pressure downstream, was measured at time $t=8$ from the
center of the vortex along the $135^{\circ}$ inclined radius with respect to
the x-axis. As expected, the sound pressure of both SLLBM simulations
coincided well with the reference solution by Inoue and Hattori [59]. These
findings approve the cubature-derived D2Q19 as a potent alternative for
compressible off-lattice simulations to the D2Q25 being based on the Gauß-
product rule used in previous works [41, 44].

Figure 4: Comparison of D2Q19 and D2Q25 in terms of the sound pressure $\Delta
P=(P-P_{B})/P_{B}$ of the shock-vortex interaction with $Re=800,Ma_{v}=0.25$
and at time $t^{\prime}=8^{\prime}$ measured from the center of the vortex
along the radius with $135^{\circ}$ with respect to the x-axis. The SLLBM
match the reference by Inoue and Hattori [59] well.

###  4.2 Weakly compressible 3D incompressible Taylor-Green vortex

The three-dimensional Taylor Green vortex served as a test case to apply the
degree-five velocity sets D3Q13, D3Q21, and degree-seven D3V27 in comparison
to the customary degree-five D3Q27 using an isothermal configuration
$\theta=1$ with truncation order of the equilibrium $N=2$ for D3Q13, D3Q21,
and D3Q27. The D3V27 simulation instead was equipped by a third truncation
order $N=3$. The Mach number was $Ma=0.1$. The computational domain is
$S=[0,2\pi]^{3}$ and the following initial conditions were used

$\displaystyle\bm{u}(\bm{x},0)=\begin{pmatrix}\mathrm{sin}(x)\mathrm{cos}(y)\mathrm{cos}(z)\\\
-\mathrm{cos}(x)\mathrm{sin}(y)\mathrm{cos}(z)\\\ 0\end{pmatrix},$ (22)
$P(\bm{x},t=0)=\frac{1}{16}(\mathrm{cos}(2x)+\mathrm{cos}(2y))\mathrm{cos}(2z+2),$
(23)

The kinetic energy of the flow is defined as

$k=\frac{1}{(2\pi)^{3}}\int_{S}\ \frac{1}{2}|\mathbf{u}|^{2}\ d^{3}\mathbf{x}$
(24)

and the enstrophy as a measure of the dissipation with respect to the resolved
scales [60]

$\mathcal{E}=\frac{1}{(2\pi)^{3}}\int_{S}\ (\nabla\times\mathbf{u})^{2}\
d^{3}\mathbf{x}.$ (25)

As explained in detail in [33] and references therein, the scaled enstrophy
$\nu\mathcal{E}$ measures the physical dissipation in terms of the local shear
stresses. In turn, its deviation from the phenomenological dissipation rate
-dk/dt provides a measure for the numerical dissipation due to discretization.

First, a well-resolved flow was tested by setting the Reynolds number to
$Re=400$. The domain was discretized by $32^{3}$ cells with order of finite
element $p=4$, resulting in $128^{3}$ grid points. The velocity sets were
scaled in such a way that the time step size was $\delta t=0.0015$ in all
cases. A summary of the flow variables is listed in Table 5.

Table 5: Simulation setup of the weakly compressible 3D Taylor-Green vortex. Test case |  | a) | b)
---|---|---|---
Reynolds number | $\mathrm{Re}$ | 400 | 1600
Mach number | $\mathrm{Ma}$ | 0.1 | 0.1
Relaxation time | $\tau$ | 0.517 | 0.508
Spatial resolution | $\delta_{x}$ | 0.05 | 0.05
Temporal resolution | $\delta_{t}$ | 0.0015 | 0.00075
Interpolation polynomial order | $p$ | 4 | 4
Truncation order of equilibrium | $N$ | 2* | 2*
*$N=3$ for D3V27 |  |  |

Figure 5 depicts the dissipation $dk/dt$ and the scaled enstrophy
$\nu\mathcal{E}$ of the four focused velocity sets in comparison to the
reference by Brachet et al. [58]. It is shown that all velocity sets match the
reference solution for both dissipation and scaled enstrophy. Thus in case of
well-resolved simulations, even the lean D3Q13 velocity set commended itself
as a good choice for three-dimensional flow simulations. In case of the
underresolved simulations with Reynolds number $\mathrm{Re}=1600$ with time
step size $\delta_{t}=0.00075$ the situation is slightly different though,
displayed in Fig. 6 showing different evolutions of dissipation and scaled
enstrophy for each of the velocity sets. Initially, the dissipation followed
the DNS reference simulation. As the flow entered the fully turbulent regime,
the dissipation became more volatile than in the resolved case due to under-
resolution of small vortices. In this phase, the numerical discretization
provided dissipation in the spirit of an implicit LES [61]. Despite the
absence of an explicit subgrid model, all stencils produced macroscopic
dissipation rates that roughly followed the reference solution. The D3Q13 and
D3Q27 stencils even captured the location of the dissipation peak within 0.5
time units, whereas D3Q21 and D3V27 produced a plateau and a premature peak,
respectively. That said, the broader implications of this difference are
debatable, since the dissipation in the under-resolved regime is mostly a
result of numerical errors interacting with the finest resolved scales. In
contrast, the scaled enstrophy is a more meaningful observable as it directly
reflects the local shear stresses. Fig. 7 compares the scaled enstrophies
revealing that the enstrophy determined by simulations with the common D3Q27
most drastically underestimated the enstrophy in this configuration. Table 6
confirms this observation by listing the Euclidean norm of the difference
between the scaled enstrophies obtained by simulation and reference, i.e.

$\|\nu\mathcal{E}^{\mathrm{ref}}-\nu\mathcal{E}^{\mathrm{sim}}\|=\sqrt{\frac{1}{\hat{N}}\sum_{j=1}^{\hat{N}}\left(\nu\mathcal{E}^{\mathrm{ref}}_{j}-\nu\mathcal{E}^{\mathrm{sim}}_{j}\right)^{2}},$
(26)

where $\hat{N}$ is the number of discrete measurements. The degree-seven D3V27
simulations as well as the degree-five D3Q21 simulations both more accurately
predicted the scaled enstrophy.

Figure 5: Dissipation and scaled enstrophy over time from simulations of the three-dimensional Taylor-Green vortex at Reynolds number $\mathrm{Re}=400$ using the velocity sets D3Q13 a), D3Q21 b), D3V27 c), and D3Q27 d). For this well-resolved simulation at $N=128\textthreesuperior$ both the dissipation $-dk/dt$ and the scaled enstrophy $\nu\mathcal{E}$ matched the reference solution well for all velocity sets considered. Figure 6: Dissipation $-kE/dt$ and scaled enstrophy over time from simulations of the three-dimensional Taylor-Green vortex at Reynolds number $\mathrm{Re}=1600$ using the velocity sets D3Q13 a), D3Q21 b), D3V27 c), and D3Q27 d). This under-resolved simulation at $N=128\textthreesuperior$ reveals differences for the dissipation and scaled enstrophy for the respective velocity sets.  Figure 7: Comparison of the scaled enstrophies $\nu\mathcal{E}$ shown in Fig. 6. Table 6: Root mean squared error of the scaled enstrophy with respect to the reference [58]. | D3Q13 | D3Q21 | D3V27 | D3Q27
---|---|---|---|---
$\|\nu\mathcal{E}^{\mathrm{ref}}-\nu\mathcal{E}^{\mathrm{sim}}\|$ | 0.00141 | 0.00116 | 0.00107 | 0.00158

###  4.3 Compressible 3D Taylor-Green vortex

As a final test case the compressible 3D Taylor-Green vortex was simulated.
Recently, this test case has been by extensively investigated by Peng and Yang
[62]. The Reynolds number was $Re=400$, the heat capacity ratio was
$\gamma=1.4$, and the Prandtl number was $\mathrm{Pr}=0.75$. The initial
velocity field corresponds to Eq. (22), whereas the pressure is obtained by
$P=\rho T_{0}$ with $T_{0}=1$, while the density is given by

$\rho(\bm{x},t=0)=1.0+\frac{\mathcal{C}^{2}}{16}(\mathrm{cos}(2x)+\mathrm{cos}(2y))\mathrm{cos}(2z+2).$
(27)

This relation of temperature and density corresponds to the constant
temperature initial condition (CTIC) detailed in Peng and Yang [62]. The
nominator $\mathcal{C}$ can be either set to $\mathcal{C}=\mathrm{Ma}^{2}$ or
to $\mathcal{C}=1$; this work made use of the latter case. Table 7 lists the
most important simulation variables and Table 8 separately lists the time step
sizes and relaxation times at $\theta=1.0$ of the respective simulations with
the velocity sets D3Q103, D3V107, and D3Q45. The time step sizes of the on-
lattice velocity sets were dictated by the configuration, whereas the time
step sizes of the D3Q45 SLLBM simulation were adjustable. The results of the
D3Q45 velocity set is discussed in more depth in a follow-up manuscript [48].

Table 7: Simulation setup of the fully compressible 3D Taylor-Green vortex. Test case |  | a) | b) | c) | d)
---|---|---|---|---|---
Reynolds number | $\mathrm{Re}$ | 400 | 400 | 400 | 400
Mach number | $\mathrm{Ma}$ | 0.5 | 1.0 | 1.5 | 2.0
Spatial resolution | $\delta_{x}$ | 0.025 | 0.025 | 0.025 | 0.025
Truncation order of equilibrium | $N$ | 4 | 4 | 4 | 4

Figure 8 depicts the mean kinetic energy over time for a resolution of
$256^{3}$ obtained by the velocity sets D3Q103, D3V107 and by the SLLBM D3Q45
simulation. It displays that the reference by Peng and Yang was well captured
for Mach numbers $\mathrm{Ma}=0.5$ and $\mathrm{Ma}=1.0$ by all three velocity
sets. However, the on-lattice simulations became unstable for the larger Mach
numbers $\mathrm{Ma}=1.5$ and $\mathrm{Ma}=2.0$.

These results confirm the assumptions of other works which predicted the
capability of the on-lattice velocity sets D3Q103 and D3V107 to simulate
compressible flows [25, 30, 29, 63]. However, simulations applying these
velocity sets lack stability, at least in the present configuration for Mach
numbers exceeding one. In contrast, the variable time step size of the SLLBM
enabled stable simulations even at higher Mach numbers using a substantially
smaller velocity set with 45 discrete velocities.

Figure 8: On-lattice simulations utilizing the D3Q103 a) and the D3V107 in
comparison to the SLLBM D3Q45 simulation c). All velocity sets accurately
reproduced the kinetic energies for the medium Mach numbers $\mathrm{Ma}=0.5$
and $\mathrm{Ma}=1.0$, while the on-lattice simulations failed for the higher
Mach numbers $\mathrm{Ma}=1.5$ and $\mathrm{Ma}=2.0$. By contrast, the
flexible time step size enabled successful SLLBM simulations even at the
highest Mach numbers.  Table 8: Time step sizes (top) and relaxation times at
$\theta=1.0$ (bottom) of the compressible Taylor-Green vortex simulations.
Unstable simulations marked by asterisks.

a) Time step size $\delta_{t}$ | D3Q103 | D3V107 | D3Q45
---|---|---|---
$\mathrm{Ma}=0.5$ | 0.0121 | 0.0136 | 0.0165
$\mathrm{Ma}=1.0$ | 0.0243 | 0.0271 | 0.0329
$\mathrm{Ma}=1.5$ | 0.0364* | 0.0407* | 0.0178
$\mathrm{Ma}=2.0$ | 0.0485* | 0.0542* | 0.0119

b) Relaxation time $\tau$ | D3Q103 | D3V107 | D3Q45
---|---|---|---
$\mathrm{Ma}=0.5$ | 0.572 | 0.564 | 0.553
$\mathrm{Ma}=1.0$ | 0.644 | 0.629 | 0.606
$\mathrm{Ma}=1.5$ | 0.716* | 0.694* | 0.943
$\mathrm{Ma}=2.0$ | 0.789* | 0.758* | 1.681

##  5 Discussion

The discretization of the velocity space by quadrature is a key asset of the
lattice Boltzmann method. However, the coupling of the velocity space
discretization with the spatial and temporal discretization in on-lattice
Boltzmann methods is obstructive to find sufficiently small velocity sets
especially for high quadrature order. Even a recent systematic study as done
by Spiller and Dünweg did not yield smaller velocity sets than the already
known D3Q103 [64], whose number of abscissae therefore appears to be the lower
limit at this quadrature degree with equidistant nodes. The enormous size of
these high-degree sets prevented their application in actual simulations, at
least up to the present study. The simulations of the fully compressible
Taylor-Green vortex in the present work showed on-lattice velocity sets are
generally capable for compressible simulations, but they lack stability at
Mach numbers beyond unity. The biggest issue, from our point of view, is the
fixed time step size of most compressible on-lattice Boltzmann methods. An
exception to this are hybrid lattice Boltzmann methods, which solve the energy
equation by finite volume or finite difference methods [65, 66]. These methods
are capable to adjust the time step size by changing the reference
temperature, but they suffer from restrictions in the stability regions [67].
Despite the successes of recent works [26, 42, 68, 69, 70], the applicability
of compressible on-lattice Boltzmann methods remains vague. Although shifted
stencils have proven to extend the Mach number range of on-lattice Boltzmann
methods [71, 72] in certain situations, the Taylor-Green vortex at high Mach
numbers would still not be stable with static reference frames. Solely dynamic
shifts of the reference frame as presented by Coreixas and Latt [73] might be
an alternative for on-lattice Boltzmann methods, although they cause
additional computational costs and require further investigation.

Compared to on-lattice Boltzmann methods, the discretization of the velocity
space is a largely unexplored field of research in off-lattice Boltzmann
methods. We were able, however, to identify a large potential in equipping
off-lattice Boltzmann methods by velocity sets with significantly different
shapes. Instead of deriving stencils case-by-case, the research on
multivariate quadrature rules yielded suitable cubature rules as templates for
velocity sets. Consequently, this work’s main purpose was to explore these
velocity sets by off-lattice Boltzmann simulations.

The results in Section 4 clearly indicate that the spatial freedom obtained in
the collocation of abscissae, as done in cubature rules, reduces the number of
discrete velocities significantly. Simultaneously, the computational costs
halved, narrowing the gap between on-lattice and off-lattice Boltzmann methods
in efficiency. For example, the D2Q19 degree-nine velocity set is
approximately half the size of its on-lattice counterpart D2V37 with equal
quadrature order [24]. The difference between the D3Q45 and the D3Q103
velocity set is even larger, but the size of the D3Q45 is also superior
compared to a recently utilized degree-nine D3Q77 off-lattice velocity set for
compressible finite volume LBMs [74]. The D3Q45 simulations of the Taylor-
Green vortex proved to be stable and accurate even for high Mach number flows.
This can be attributed to the temporal discretization of off-lattice Boltzmann
methods, which is independent of the discretization of space and velocity
space.

Differences were also identified in terms of accuracy. The D3Q13 and D3Q21
based on icodahedra and dodecahedra, shown in Fig. 2, proved to be both
efficient and accurate off-lattice velocity sets for weakly compressible
flows. Both velocity sets showed a better enstrophy agreement with the
reference solution in underresolved flow simulations. The reason for this is
not necessarily that the precision of the cubature is higher, since the degree
of precision of the D3Q13, D3Q21, and D3Q27 velocity sets is identical.
Instead, we attribute the better agreement of the enstrophy with the reference
to the fewer advection equations that have to be solved connected with fewer
interpolation steps, leading to less artificial diffusion. Moreover, due to
their derivation from platonic solids, these velocity sets exhibit a high
geometric isotropy, which is favourable for the simulation of turbulent flows.
This is probably also the reason why the D3Q21 outperforms the D3Q13 velocity
set in terms of enstrophy. However, simulations with all of these degree-five
velocity sets are affected by the cubic error term, being dependent on the
Mach number. Obviously, this error term also spoils the simulations, since the
Mach number of Ma=0.1 is not negligible. The degree-seven velocity set D3V27
with the aim to eliminate these cubic errors in the model proved to be more
accurate than the on-lattice D3Q27, but coming at a comparable price. Despite
the numerical diffusion introduced by the SLLBM, we showed in past works that
the spatial convergence order corresponds to the order of the finite element
[45]. By contrast, the spatial convergence order of on-lattice Boltzmann
methods is restricted to second order.

Numerical diffusion and dispersion are critical issues of off-lattice
Boltzmann methods compared to on-lattice Boltzmann methods, resulting from the
interpolation in case of the SLLBM. High-order interpolation polynomials as
used in the present article, however, significantly lower the dissipative
effects in the simulation [45]. Contrary to that, on-lattice Boltzmann methods
are valued for low dissipation [75] due to the exact streaming step, but they
suffer from instabilities at high Reynolds numbers when applying the usual BGK
collision operator. To encounter these stability issues, a number of collision
models have been proposed to stabilize simulations by introducing minimal
diffusion without sacrificing the locality of the collision operator. Examples
are multi-relaxation time (MRT) [76, 77, 78, 79], central MRT [80, 81],
entropic [82, 83], regularized [84, 85], or cumulant models [86, 87]. Due to
the non-negligible amount of dissipation evoked by the streaming step of off-
lattice Boltzmann schemes, it is likely that these established on-lattice
collision models fulfil a different role in off-lattice LBM. Nevertheless,
previous work has utilized the semi-Lagrangian streaming step in combination
with a stabilized collision model [40]. Consequently, future research is
necessary to clarify the impact of collision models on off-lattice Boltzmann
methods.

Finally, some technical aspects need to be discussed. In particular, the
streaming step’s implementations of on-lattice and off-lattice Boltzmann
methods are significantly different. On the one hand, the exact streaming step
of on-lattice Boltzmann methods can naively be implemented by on-board
routines of many software packages, e.g. by “numpy.roll” in Numpy [88]. On the
other hand, the interpolation routines of the SLLBM require a thoughtful
implementation, especially when following our recommendation to employ high-
order polynomials with non-equidistant support points and cell-wise
organization. For this task, a matured finite element library like deal.ii is
helpful [54]. The execution time of the semi-Lagrangian streaming step is one
magnitude larger [45] for the same configuration and velocity set.. Despite
the computational overhead, off-lattice Boltzmann methods have proven
beneficial for unstructured meshes or compressible simulations [40, 45, 44].
More appropriate velocity sets as presented in this article further reduce the
gap in computation cost, as the computational complexity of the SLLBM scales
linearly with respect to the number of discrete velocities.

To sum up, off-lattice Boltzmann methods relax the LBM in terms of the
velocity space discretization. As the development of sparse numerical
cubatures is still progressing, the present work lays the foundations to
directly utilize new cubature rules in CFD simulations with off-lattice
Boltzmann methods.

##  6 Conclusion

This paper studied the utilization of cubature rules in off-lattice Boltzmann
methods for weakly and fully compressible flows. The deduced off-lattice
velocity sets presented in this article reduce the number of discrete
velocities and increase the accuracy of the simulation. In addition, fully
compressible off-lattice simulations with degree-nine velocity sets feature
better stability compared to the corresponding on-lattice counterparts. Taken
together, cubature rule-derived off-lattice velocity sets are a viable
alternative to the customary on-lattice Boltzmann velocity sets and should
have priority in off-lattice Boltzmann simulations.

## Acknowledgements

We gratefully acknowledge support for D.W. by German Research Foundation (DFG)
project FO 674/17-1. The simulations were performed using the Platform for
Scientific Computing at Bonn-Rhein-Sieg University, which is funded by the
German Ministry of Education and Research and the Ministry for Culture and
Science North Rhine-Westfalia (research grant 13FH156IN6).

##  Appendix A Hermite polynomials and moments

The scaled Hermite polynomials with $\hat{\xi_{i}}=\xi_{i}/c_{s}$ up to fourth
order read

$\displaystyle{\mathcal{H}_{i}}^{(0)}$ $\displaystyle=1$
$\displaystyle{\mathcal{H}_{i\alpha}}^{(1)}$
$\displaystyle=\frac{{{\hat{\xi}_{i\alpha}}}}{c_{s}}$
$\displaystyle{\mathcal{H}_{i\alpha\beta}}^{(2)}$
$\displaystyle=\frac{{{{\hat{\xi}_{i\alpha}}{\hat{\xi}_{i\beta}}}}-\delta_{\alpha\beta}}{c_{s}^{2}}$
$\displaystyle{\mathcal{H}_{i\alpha\beta\gamma}}^{(3)}$
$\displaystyle=\frac{{\hat{\xi}_{i\alpha}}{\hat{\xi}_{i\beta}}{\hat{\xi}_{i\gamma}}-({\hat{\xi}_{i\alpha}\delta_{\beta\gamma}}+{\hat{\xi}_{i\beta}\delta_{\alpha\gamma}}+{\hat{\xi}_{i\gamma}\delta_{\alpha\beta}})}{c_{s}^{3}}$
$\displaystyle{\mathcal{H}_{i\alpha\beta\gamma\delta}}^{(4)}$
$\displaystyle=\frac{{{\hat{\xi}_{i\alpha}}{\hat{\xi}_{i\beta}}{\hat{\xi}_{i\gamma}}{\hat{\xi}_{i\delta}}}-\mathcal{T}_{i}+(\delta_{\alpha\beta}\delta_{\gamma\delta}+\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma})}{c_{s}^{4}},$

with

$\begin{split}\mathcal{T}_{i}=\hat{\xi}_{i\alpha}\hat{\xi}_{i\beta}\delta_{\gamma\delta}+\hat{\xi}_{i\alpha}\hat{\xi}_{i\gamma}\delta_{\beta\delta}+\hat{\xi}_{i\alpha}\hat{\xi}_{i\delta}\delta_{\beta\gamma}+\\\
\hat{\xi}_{i\beta}\hat{\xi}_{i\gamma}\delta_{\alpha\delta}+\hat{\xi}_{i\beta}\hat{\xi}_{i\delta}\delta_{\alpha\gamma}+\hat{\xi}_{i\gamma}\hat{\xi}_{i\delta}\delta_{\alpha\beta}.\end{split}$

The moments of the Boltzmann equation up to fourth order are

$\displaystyle a^{(0)}_{eq}$ $\displaystyle=$ $\displaystyle\rho$
$\displaystyle a^{(1)}_{\alpha,eq}$ $\displaystyle=$ $\displaystyle\rho
u_{\alpha}$ $\displaystyle a^{(2)}_{\alpha\beta,eq}$ $\displaystyle=$
$\displaystyle\rho(u_{\alpha}u_{\beta}+T_{0}(\theta-1)\delta_{\alpha\beta})$
$\displaystyle a^{(3)}_{\alpha\beta\gamma,eq}$ $\displaystyle=$
$\displaystyle\rho\left[u_{\alpha}u_{\beta}u_{\gamma}+T_{0}(\theta-1)(\delta_{\alpha\beta}u_{\gamma}\right.$
$\displaystyle\left.\hskip
34.14322pt+\delta_{\alpha\gamma}u_{\beta}+\delta_{\beta\gamma}u_{\alpha})\right]$
$\displaystyle a^{(4)}_{\alpha\beta\gamma\delta,eq}$ $\displaystyle=$
$\displaystyle\rho[u_{\alpha}u_{\beta}u_{\gamma}u_{\delta}+T_{0}(\theta-1)((\delta_{\alpha\beta}\delta_{\gamma\delta}$
$\displaystyle+\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma})T_{0}(\theta-1)$
$\displaystyle+\delta_{\alpha\beta}u_{\gamma}u_{\delta}+\delta_{\alpha\gamma}u_{\beta}u_{\delta}+\delta_{\alpha\delta}u_{\beta}u_{\gamma}+$
$\displaystyle\delta_{\beta\gamma}u_{\alpha}u_{\delta}+\delta_{\beta\delta}u_{\alpha}u_{\gamma}+\delta_{\gamma\delta}u_{\alpha}u_{\beta})].\vspace{.3cm}$

##  Appendix B Symmetry conditions

The symmetry conditions from zeroth to ninth order are

$\displaystyle\sum_{i}w_{i}$ $\displaystyle=1$
$\displaystyle\sum_{i}w_{i}\xi_{i\alpha}$ $\displaystyle=0$
$\displaystyle\sum_{i}w_{i}\xi_{i\alpha}\xi_{i\beta}$
$\displaystyle=c_{s}^{2}\delta_{\alpha\beta}$
$\displaystyle\sum_{i}w_{i}\xi_{i\alpha}\xi_{i\beta}\xi_{i\gamma}$
$\displaystyle=0$
$\displaystyle\sum_{i}w_{i}\xi_{i\alpha}\xi_{i\beta}\xi_{i\gamma}\xi_{i\delta}$
$\displaystyle=c_{s}^{4}\Delta_{\alpha\beta\gamma\delta}$
$\displaystyle\sum_{i}w_{i}\xi_{i\alpha}\xi_{i\beta}\xi_{i\gamma}\xi_{i\delta}\xi_{i\epsilon}$
$\displaystyle=0$
$\displaystyle\sum_{i}w_{i}\xi_{i\alpha}\xi_{i\beta}\xi_{i\gamma}\xi_{i\delta}\xi_{i\epsilon}\xi_{i\zeta}$
$\displaystyle=c_{s}^{6}\Delta_{\alpha\beta\gamma\delta\epsilon\zeta}$
$\displaystyle\sum_{i}w_{i}\xi_{i\alpha}\xi_{i\beta}\xi_{i\gamma}\xi_{i\delta}\xi_{i\epsilon}\xi_{i\zeta}\xi_{i\eta}$
$\displaystyle=0$
$\displaystyle\sum_{i}w_{i}\xi_{i\alpha}\xi_{i\beta}\xi_{i\gamma}\xi_{i\delta}\xi_{i\epsilon}\xi_{i\zeta}\xi_{i\eta}\xi_{i\theta}$
$\displaystyle=c_{s}^{8}\Delta_{\alpha\beta\gamma\delta\epsilon\zeta\eta\theta}$
$\displaystyle\sum_{i}w_{i}\xi_{i\alpha}\xi_{i\beta}\xi_{i\gamma}\xi_{i\delta}\xi_{i\epsilon}\xi_{i\zeta}\xi_{i\eta}\xi_{i\theta}\xi_{i\iota}$
$\displaystyle=0,$

with

$\Delta_{\alpha\beta\gamma\delta}=\delta_{\alpha\beta}\delta_{\gamma\delta}+\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma},$

as well as

$\Delta_{\alpha\beta\gamma\delta\epsilon\zeta}=\delta_{\alpha\beta}\Delta_{\gamma\delta\epsilon\zeta}\\\
+\delta_{\alpha\gamma}\Delta_{\beta\delta\epsilon\zeta}+\delta_{\alpha\delta}\Delta_{\beta\gamma\epsilon\zeta}+\delta_{\alpha\epsilon}\Delta_{\beta\gamma\delta\zeta}+\delta_{\alpha\zeta}\Delta_{\beta\gamma\delta\epsilon},$

and

$\Delta_{\alpha\beta\gamma\delta\epsilon\zeta\eta\theta}=\delta_{\alpha\beta}\Delta_{\gamma\delta\epsilon\zeta\eta\theta}+\delta_{\alpha\gamma}\Delta_{\beta\delta\epsilon\zeta\eta\theta}+\delta_{\alpha\delta}\Delta_{\beta\gamma\epsilon\zeta\eta\theta}+\delta_{\alpha\epsilon}\Delta_{\beta\gamma\delta\zeta\eta\theta}\\\
+\delta_{\alpha\zeta}\Delta_{\beta\gamma\delta\epsilon\eta\theta}+\delta_{\alpha\eta}\Delta_{\beta\gamma\delta\epsilon\zeta\theta}+\delta_{\alpha\theta}\Delta_{\beta\gamma\delta\epsilon\zeta\eta}.$

##  Appendix C D2Q19 velocity set

$i$ | $w$ | $\xi_{x}$ | $\xi_{y}$
---|---|---|---
$0$ | 0.3168437267921905 | 0.0 | 0.0
$1,2$ | 0.10558878375062891 | $\pm 1.4869982213169028$ | 0.0
$3,4,5,6$ | 0.1024247123210936 | $\pm 0.775196278121181$ | $\pm 1.367469636752619$
$7,8,9,10$ | 0.00953510698543825 | $\pm 2.5175897644357486$ | $\pm 1.105629214668943$
$11,12$ | 0.006865104210104631 | $0.0$ | $\pm 2.9213306655318734$
$13,14,15,16$ | 0.002405335328939458 | $\pm 1.8663975507141328$ | $\pm 2.6987507639352253$
$17,18$ | 0.0003939393722285871 | $\pm 3.8358342053914734$ | $0.0$
$c_{s}=1$ |  |  |

##  Appendix D D3V27 velocity set

$i$ | $w$ | $\xi_{x}$ | $\xi_{y}$ | $\xi_{z}$ |
---|---|---|---|---|---
$0$ | 0.31247897198654906 | 0.0 | 0.0 | 0.0 |
$1,\ldots,8$ | 0.06338446047675325 | $\pm 1.1198362860638005$ | $\pm 1.1198362860638005$ |  $\pm$ 1.1198362860638005 |
$9,\ldots,14$ | 0.029035130153906134 | $\pm 2.358709038202103$ | 0.0 | 0.0 | (cyc)
$15,\ldots,26$ | 0.0005195469396656799 | $\pm 3.142130383387586$ | $\pm 3.142130383387586$ | 0.0 | (cyc)
$c_{s}=1$ |  |  |  |  |

##  Appendix E D3Q45 velocity set

$w$ | $\xi_{x}$ | $\xi_{y}$ | $\xi_{z}$
---|---|---|---
0.20740740740740618 | 0.0 | 0.0 | 0.0
0.05787037037037047 | 0.06386083877343968 | -1.2239121278243665 | -1.2239121278243665
0.05787037037037047 | -0.06386083877343968 | 1.2239121278243665 | 1.2239121278243665
0.05787037037037047 | 1.2239121278243665 | -0.06386083877343968 | 1.2239121278243665
0.05787037037037047 | -1.2239121278243665 | 0.06386083877343968 | -1.2239121278243665
0.05787037037037047 | 1.2239121278243665 | 1.2239121278243665 | -0.06386083877343968
0.05787037037037047 | -1.2239121278243665 | -1.2239121278243665 | 0.06386083877343968
0.05787037037037047 | 1.5766994272507744 | -0.5069610024977665 | -0.5069610024977665
0.05787037037037047 | -1.5766994272507744 | 0.5069610024977665 | 0.5069610024977665
0.05787037037037047 | 0.5069610024977665 | 0.5069610024977665 | -1.5766994272507744
0.05787037037037047 | -0.5069610024977665 | -0.5069610024977665 | 1.5766994272507744
0.05787037037037047 | -0.5069610024977665 | 1.5766994272507744 | -0.5069610024977665
0.05787037037037047 | 0.5069610024977665 | -1.5766994272507744 | 0.5069610024977665
0.00462962962962958 | 2.403092127540177 | 0.8892242114059369 | -1.5602655313772367
0.00462962962962958 | -2.403092127540177 | -0.8892242114059369 | 1.5602655313772367
0.00462962962962958 | -2.403092127540177 | 1.5602655313772367 | -0.8892242114059369
0.00462962962962958 | 2.403092127540177 | -1.5602655313772367 | 0.8892242114059369
0.00462962962962958 | -0.8892242114059369 | -2.403092127540177 | 1.5602655313772367
0.00462962962962958 | 0.8892242114059369 | 2.403092127540177 | -1.5602655313772367
0.00462962962962958 | -0.8892242114059369 | 1.5602655313772367 | -2.403092127540177
0.00462962962962958 | 0.8892242114059369 | -1.5602655313772367 | 2.403092127540177
0.00462962962962958 | -1.5602655313772367 | 0.8892242114059369 | 2.403092127540177
0.00462962962962958 | -1.5602655313772367 | 2.403092127540177 | 0.8892242114059369
0.00462962962962958 | 1.5602655313772367 | -0.8892242114059369 | -2.403092127540177
0.00462962962962958 | 1.5602655313772367 | -2.403092127540177 | -0.8892242114059369
0.00462962962962958 | 0.4744978678080795 | 0.4744978678080795 | 2.9239876105912574
0.00462962962962958 | 0.4744978678080795 | 2.9239876105912574 | 0.4744978678080795
0.00462962962962958 | -0.4744978678080795 | -0.4744978678080795 | -2.9239876105912574
0.00462962962962958 | -0.4744978678080795 | -2.9239876105912574 | -0.4744978678080795
0.00462962962962958 | 2.9239876105912574 | 0.4744978678080795 | 0.4744978678080795
0.00462962962962958 | -2.9239876105912574 | -0.4744978678080795 | -0.4744978678080795
0.00462962962962958 | 1.7320508075688787 | 1.7320508075688787 | 1.7320508075688787
0.00462962962962958 | -1.7320508075688787 | -1.7320508075688787 | -1.7320508075688787
0.0004629629629629939 | -2.7367507163016924 | 0.14279717659756475 | -2.7367507163016924
0.0004629629629629939 | 2.7367507163016924 | 2.7367507163016924 | -0.14279717659756475
0.0004629629629629939 | 2.7367507163016924 | -0.14279717659756475 | 2.7367507163016924
0.0004629629629629939 | -2.7367507163016924 | -2.7367507163016924 | 0.14279717659756475
0.0004629629629629939 | 0.14279717659756475 | -2.7367507163016924 | -2.7367507163016924
0.0004629629629629939 | -0.14279717659756475 | 2.7367507163016924 | 2.7367507163016924
0.0004629629629629939 | -3.5256070994177073 | 1.1335992635264445 | 1.1335992635264445
0.0004629629629629939 | 3.5256070994177073 | -1.1335992635264445 | -1.1335992635264445
0.0004629629629629939 | 1.1335992635264445 | -3.5256070994177073 | 1.1335992635264445
0.0004629629629629939 | -1.1335992635264445 | 3.5256070994177073 | -1.1335992635264445
0.0004629629629629939 | 1.1335992635264445 | 1.1335992635264445 | -3.5256070994177073
0.0004629629629629939 | -1.1335992635264445 | -1.1335992635264445 | 3.5256070994177073
$c_{s}=1$ |  |  |

## References

  * [1] Guy R. McNamara and Gianluigi Zanetti.  Use of the boltzmann equation to simulate lattice-gas automata.  Phys. Rev. Lett., 1988.
  * [2] F. J. Higuera and J. Jiménez.  Boltzmann approach to lattice gas simulations.  EPL, 1989.
  * [3] Shiyi Chen, Hudong Chen, Daniel Martnez, and William Matthaeus.  Lattice Boltzmann model for simulation of magnetohydrodynamics.  Phys. Rev. Lett., 1991.
  * [4] Timm Krüger, Halim Kusumaatmaja, Alexander Kuzmin, Orest Shardt, Goncalo Silva, and Erlend Magnus Viggen.  The lattice Boltzmann method: Principles and practice.  Springer, 2017.
  * [5] Pierre Lallemand, Li-shi Luo, Manfred Krafczyk, and Wen-an Yong.  The Lattice Boltzmann Method for Nearly Incompressible Flows.  J. Comput. Phys., page 109713, 2020.
  * [6] Tadeusz Platkowski and Reinhard Illner.  Discrete velocity models of the Boltzmann equation: A survey on the mathematical aspects of the theory.  SIAM Review, 30(2):213–255, 1988.
  * [7] Li Shi Luo.  Some recent results on discrete velocity models and ramifications for lattice Boltzmann equation.  Computer Physics Communications, 129(1):63–74, 2000.
  * [8] Mingde Su, Kun Xu, and M. S. Ghidaoui.  Low-Speed Flow Simulation by the Gas-Kinetic Scheme.  Journal of Computational Physics, 1999.
  * [9] Kun Xu and Juan Chen Huang.  A unified gas-kinetic scheme for continuum and rarefied flows.  Journal of Computational Physics, 229(20):7747–7764, 2010.
  * [10] Xiaoyi He, Xiaowen Shan, and Gary D. Doolen.  Discrete Boltzmann equation model for nonideal gases.  Phys. Rev. E - Stat. Physics, Plasmas, Fluids, Relat. Interdiscip. Top., 57(1):R13–R16, 1998.
  * [11] Christophe Coreixas, Gauthier Wissocq, Guillaume Puigt, Jean François Boussuge, and Pierre Sagaut.  Recursive regularization step for high-order lattice Boltzmann methods.  Phys. Rev. E, 96(3), 2017.
  * [12] Xiaowen Shan and Xiaoyi He.  Discretization of the velocity space in the solution of the Boltzmann equation.  Phys. Rev. Lett., 80(1):65–68, jan 1998.
  * [13] Xiaowen Shan, Xue Feng Yuan, and Hudong Chen.  Kinetic theory representation of hydrodynamics: A way beyond the Navier-Stokes equation.  J. Fluid Mech., 550(-1):413–441, feb 2006.
  * [14] Ienkaran Arasaratnam and Simon Haykin.  Cubature kalman filters.  IEEE Trans. Automat. Contr., 54(6):1254–1269, 2009.
  * [15] A. H. Stroud.  Approximate Calculation of Multiple Integrals.  Math. Comput., 1973.
  * [16] Ronald Cools and Philip Rabinowitz.  Monomial cubature rules since ”Stroud”: a compilation.  J. Comput. Appl. Math., 48(3):309–326, 1993.
  * [17] Ronald Cools.  An encyclopaedia of cubature formulas.  Journal of Complexity, 19(3):445 – 453, 2003.  Oberwolfach Special Issue.
  * [18] Xiaoyi He and Li Shi Luo.  Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation.  Phys. Rev. E - Stat. Physics, Plasmas, Fluids, Relat. Interdiscip. Top., 55(6):6811–6820, 1997.
  * [19] Alexander Thomas White and Chuh Khiun Chong.  Rotational invariance in the three-dimensional lattice Boltzmann method is dependent on the choice of lattice.  J. Comput. Phys., 230(16):6367–6378, 2011.
  * [20] Goncalo Silva and Viriato Semiao.  Truncation errors and the rotational invariance of three-dimensional lattice models in the lattice Boltzmann method.  J. Comput. Phys., 2014.
  * [21] Martin Bauer, Goncalo Silva, and Ulrich Rüde.  Truncation errors of the d3q19 lattice model for the lattice boltzmann method.  Journal of Computational Physics, 405:109111, 2020.
  * [22] Manjusha Namburi, Siddharth Krithivasan, and Santosh Ansumali.  Crystallographic Lattice Boltzmann Method.  Sci. Rep., 6, 2016.
  * [23] Y. H. Qian and S. A. Orszag.  Lattice BGK Models for the Navier-Stokes Equation: Nonlinear Deviation in Compressible Regimes.  Europhysics Letters (EPL), 21(3):255–259, jan 1993.
  * [24] Paulo C. Philippi, Luiz A. Hegele, Luís O.E. Dos Santos, and Rodrigo Surmas.  From the continuous to the lattice Boltzmann equation: The discretization problem and thermal models.  Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., 73(5), 2006.
  * [25] Xiaowen Shan.  The mathematical structure of the lattices of the lattice Boltzmann method.  J. Comput. Sci., 17:475–481, 2016.
  * [26] N. Frapolli, S. S. Chikatamarla, and I. V. Karlin.  Entropic lattice Boltzmann model for compressible flows.  Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., 92(6):061301, dec 2015.
  * [27] Jianping Meng and Yonghao Zhang.  Accuracy analysis of high-order lattice Boltzmann models for rarefied gas flows.  J. Comput. Phys., 230(3):835–849, 2011.
  * [28] Shyam S. Chikatamarla and Iliya V. Karlin.  Lattices for the lattice Boltzmann method.  Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., 79(4), 2009.
  * [29] Xiaowen Shan.  General solution of lattices for Cartesian lattice Bhatanagar-Gross-Krook models.  Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., 81(3):1–7, 2010\.
  * [30] R. Surmas, Capico E. Pico Ortiz, and P. C. Philippi.  Simulating thermohydrodynamics by finite difference solutions of the Boltzmann equation.  European Physical Journal: Special Topics, 171(1):81–90, 2009.
  * [31] Wahyu Perdana Yudistiawan, Sang Kyu Kwak, D. V. Patil, and Santosh Ansumali.  Higher-order Galilean-invariant lattice Boltzmann model for microflows: Single-component gas.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 82(4), 2010.
  * [32] Akinori Tamura, Keita Okuyama, Shiro Takahashi, and Masaya Ohtsuka.  Three-dimensional discrete-velocity BGK model for the incompressible Navier-Stokes equations.  Comput. Fluids, 40(1):149–155, 2011.
  * [33] Andreas Krämer, Knut Küllmer, Dirk Reith, Wolfgang Joppich, and Holger Foysi.  Semi-Lagrangian off-lattice Boltzmann method for weakly compressible flows.  Phys. Rev. E, 95(2), 2017.
  * [34] Nianzheng Cao, Shiyi Chen, Shi Jin, and Daniel Martínez.  Physical symmetry and lattice symmetry in the lattice Boltzmann method.  Phys. Rev. E - Stat. Physics, Plasmas, Fluids, Relat. Interdiscip. Top., 55(1):R21–R24, 1997.
  * [35] Kazem Hejranfar, Mohammad Hossein Saadat, and Sina Taheri.  High-order weighted essentially nonoscillatory finite-difference formulation of the lattice Boltzmann method in generalized curvilinear coordinates.  Phys. Rev. E, 95(2):23314, 2017.
  * [36] Francesca Nannelli and Sauro Succi.  The lattice Boltzmann equation on irregular lattices.  J. Stat. Phys., 68(3-4):401–407, 1992.
  * [37] Weidong Li and Li-Shi Luo.  An implicit block lu-sgs finite-volume lattice-boltzmann scheme for steady flows on arbitrary unstructured meshes.  Journal of Computational Physics, 327:503 – 518, 2016.
  * [38] Xing Shi, Jianzhong Lin, and Zhaosheng Yu.  Discontinuous Galerkin spectral element lattice Boltzmann method on triangular element.  Int. J. Numer. Methods Fluids, 42(11):1249–1261, 2003.
  * [39] Misun Min and Taehun Lee.  A spectral-element discontinuous Galerkin lattice Boltzmann method for nearly incompressible flows.  J. Comput. Phys., 230(1):245–259, 2011.
  * [40] G Di Ilio, B Dorschner, G Bella, S Succi, and I V Karlin.  Simulation of turbulent flows with the entropic multirelaxation time lattice Boltzmann method on body-fitted meshes.  J. Fluid Mech., 849:35–56, aug 2018.
  * [41] B. Dorschner, F. Bösch, and I. V. Karlin.  Particles on Demand for Kinetic Theory.  Phys. Rev. Lett., 121(13), 2018.
  * [42] Mohammad Hossein Saadat, Fabian Bösch, and Ilya V. Karlin.  Lattice Boltzmann model for compressible flows on standard lattices: Variable Prandtl number and adiabatic exponent.  Phys. Rev. E, 99(1):013306, 2019.
  * [43] E. Reyhanian, B. Dorschner, and I. V. Karlin.  Thermokinetic lattice Boltzmann model of nonideal fluids.  Phys. Rev. E, 102(2):1–5, 2020.
  * [44] Dominik Wilde, Andreas Krämer, Dirk Reith, and Holger Foysi.  Semi-Lagrangian lattice Boltzmann method for compressible flows.  Phys. Rev. E, 101:053306, May 2020.
  * [45] Andreas Krämer, Dominik Wilde, Knut Küllmer, Dirk Reith, Holger Foysi, and Wolfgang Joppich.  Lattice Boltzmann simulations on irregular grids: Introduction of the NATriuM library.  Comput. Math. with Appl., 79(1):34–54, 2020.
  * [46] M. H. Saadat, F. Bösch, and I. V. Karlin.  Semi-Lagrangian lattice Boltzmann model for compressible flows on unstructured meshes.  Physical Review E, 101(2):23311, feb 2020.
  * [47] Xiaobo Nie, Xiaowen Shan, and Hudong Chen.  Thermal lattice Boltzmann model for gases with internal degrees of freedom.  Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., 77(3):1–4, 2008\.
  * [48] Dominik Wilde, Andreas Krämer, Dirk Reith, and Holger Foysi.  High-order semi-Lagrangian kinetic scheme for compressible turbulence.  Submitted to Phys. Rev. E, 2020.
  * [49] SI Konyaev.  Ninth-order quadrature formulas invariant with respect to the icosahedral group.  In Doklady Akademii Nauk, volume 233, pages 784–787. Russian Academy of Sciences, 1977.
  * [50] James R. Van Zandt.  Efficient cubature rules.  ETNA - Electron. Trans. Numer. Anal., 51:219–239, 2019.
  * [51] U. Frisch, B. Hasslacher, and Y. Pomeau.  Lattice-gas automata for the Navier-Stokes equation.  Phys. Rev. Lett., 56(14):1505–1508, 1986.
  * [52] Ann Haegemans and Robert Piessens.  Construction of cubature formulas of degree seven and nine symmetric planar regions, using orthogonal polynomials.  SIAM Journal on Numerical Analysis, 14(3):492–508, 1977.
  * [53] A. H. Stroud.  Some Seventh Degree Integration Formulas for Symmetric Regions.  SIAM J. Numer. Anal., 4(1):37–44, 1967.
  * [54] Daniel Arndt, Wolfgang Bangerth, Bruno Blais, Thomas C. Clevenger, Marc Fehling, Alexander V. Grayver, Timo Heister, Luca Heltai, Martin Kronbichler, Matthias Maier, Peter Munch, Jean-Paul Pelteret, Reza Rastak, Ignacio Thomas, Bruno Turcksin, Zhuoran Wang, and David Wells.  The deal.II library, version 9.2.  Journal of Numerical Mathematics, 28(3):131–146, 2020.
  * [55] Andreas Krämer, Dominik Wilde, and Mario Bedrunka.  Lettuce: PyTorch-based Lattice Boltzmann Solver, April 2020.
  * [56] Mario Bedrunka, Dominik Wilde, Martin Kliemank, Dirk Reith, Holger Foysi, and Andreas Krämer.  Neural Lattice Boltzmann Framework for Machine Learning-Enhanced Simulations.
  * [57] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  PyTorch: An Imperative Style, High-Performance Deep Learning Library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [58] Marc E. Brachet, Daniel I. Meiron, Steven A. Orszag, B. G. Nickel, Rudolf H. Morf, and Uriel Frisch.  Small-scale structure of the Taylor–Green vortex.  J. Fluid Mech., 130(-1):411, may 1983.
  * [59] Osamu Inoue and Yuji Hattori.  Sound generation by shock-vortex interactions.  Journal of Fluid Mechanics, 380:81–116, 1999.
  * [60] Andreas Krämer, Dominik Wilde, Knut Küllmer, Dirk Reith, and Holger Foysi.  Pseudoentropic derivation of the regularized lattice Boltzmann method.  Phys. Rev. E, 100(2):023302, aug 2019.
  * [61] Nikolaus A. Adams and Stefan Hickel.  Implicit large-eddy simulation: Theory and application.  In Advances in Turbulence XII, pages 743–750. Springer Berlin Heidelberg, Berlin, Heidelberg, 2009.
  * [62] Naifu Peng and Yue Yang.  Effects of the Mach number on the evolution of vortex-surface fields in compressible Taylor-Green flows.  Phys. Rev. Fluids, 3(1):1–21, 2018.
  * [63] Christophe Coreixas.  High-order extension of the recursive regularized lattice Boltzmann method.  PhD thesis, INP Toulouse, 2018.
  * [64] Dominic Spiller and Burkhard Dünweg.  Semiautomatic construction of lattice boltzmann models.  Phys. Rev. E, 101:043310, Apr 2020.
  * [65] Yongliang Feng, Pierre Sagaut, and Wen Quan Tao.  A compressible lattice Boltzmann finite volume model for high subsonic and transonic flows on regular lattices.  Comput. Fluids, 131:45–55, 2016.
  * [66] Yongliang Feng, Pierre Boivin, Jérôme Jacob, and Pierre Sagaut.  Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows.  Journal of Computational Physics, 394:82–99, 2019.
  * [67] Florian Renard, G. Wissocq, Jean franccois Boussuge, and P. Sagaut.  A linear stability analysis of compressible hybrid lattice boltzmann methods.  arXiv: Computational Physics, 2020.
  * [68] Nicolò Frapolli, S. S. Chikatamarla, and I. V. Karlin.  Entropic lattice Boltzmann model for gas dynamics: Theory, boundary conditions, and implementation.  Phys. Rev. E, 93(6):063302, jun 2016.
  * [69] Mohammad Atif, Manjusha Namburi, and Santosh Ansumali.  Higher-order lattice Boltzmann model for thermohydrodynamics.  Physical Review E, 98(5):053311, 2018.
  * [70] Jonas Latt, Christophe Coreixas, Joël Beny, and Andrea Parmigiani.  Efficient supersonic flow simulations using lattice Boltzmann methods based on numerical equilibria.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 378(2175):20190559, 2020.
  * [71] Nicolò Frapolli, S. S. Chikatamarla, and I. V. Karlin.  Lattice Kinetic Theory in a Comoving Galilean Reference Frame.  Physical Review Letters, 117(1):10604, jun 2016.
  * [72] S. A. Hosseini, C. Coreixas, N. Darabiha, and D. Thévenin.  Extensive analysis of the lattice Boltzmann method on shifted stencils.  Physical Review E, 100(6):1–15, 2019.
  * [73] Christophe Coreixas and Jonas Latt.  Compressible lattice Boltzmann methods with adaptive velocity stencils: An interpolation-free formulation.  Physics of Fluids, 2020.
  * [74] Zhaoli Guo and Jianchun Wang.  Simulation of three-dimensional compressible turbulence using the discrete unified gas kinetic scheme.  Phys. Fluids, 125104, 2020.
  * [75] Simon Marié, Denis Ricot, and Pierre Sagaut.  Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics.  J. Comput. Phys., 228(4):1056–1070, 2009.
  * [76] Y. H. Qian, D D’Humières, and P. Lallemand.  Lattice BGK Models for Navier-Stokes Equation.  Europhys. Lett., 17(6):479–484, feb 1992.
  * [77] Dominique D’Humières, Irina Ginzburg, Manfred Krafczyk, Pierre Lallemand, and Li Shi Luo.  Multiple-relaxation-time lattice Boltzmann models in three dimensions.  Philos. Trans. R. Soc. A Math. Phys. Eng. Sci., 360(1792):437–451, 2002.
  * [78] Paul J. Dellar.  Incompressible limits of lattice Boltzmann equations using multiple relaxation times.  J. Comput. Phys., 190(2):351–370, 2003.
  * [79] Xiaowen Shan and Hudong Chen.  A General Multiple-Relaxation-Time Boltzmann Collision Model.  Int. J. Mod. Phys. C, 18(04):635–643, apr 2007.
  * [80] Martin Geier, Andreas Greiner, and Jan G. Korvink.  Cascaded digital lattice Boltzmann automata for high Reynolds number flow.  Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., 73(6):1–10, 2006\.
  * [81] Xiaowen Shan.  Central-moment-based galilean-invariant multiple-relaxation-time collision model.  Phys. Rev. E, 100:043308, Oct 2019.
  * [82] I. V. Karlin, A. Ferrante, and H. C. Öttinger.  Perfect entropy functions of the Lattice Boltzmann method.  Europhys. Lett., 47(2):182–188, jul 1999.
  * [83] Mohammad Atif, Praveen Kumar Kolluru, Chakradhar Thantanapally, and Santosh Ansumali.  Essentially Entropic Lattice Boltzmann Model.  Phys. Rev. Lett., 119(24):1–5, 2017.
  * [84] Jonas Latt and Bastien Chopard.  Lattice Boltzmann method with regularized pre-collision distribution functions.  Math. Comput. Simul., 72(2-6):165–168, sep 2006.
  * [85] Orestis Malaspinas.  Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization.  pages 1–31, 2015.
  * [86] Martin Geier, Martin Schönherr, Andrea Pasquali, and Manfred Krafczyk.  The cumulant lattice Boltzmann equation in three dimensions: Theory and validation.  Comput. Math. with Appl., 70(4):507–547, aug 2015.
  * [87] Martin Geier, Stephan Lenz, Martin Schönherr, and Manfred Krafczyk.  Under-resolved and large eddy simulations of a decaying Taylor–Green vortex with the cumulant lattice Boltzmann method.  Theor. Comput. Fluid Dyn., 2020.
  * [88] Charles R. Harris, K. Jarrod Millman, St’efan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fern’andez del R’ıo, Mark Wiebe, Pearu Peterson, Pierre G’erard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.