Investigation of an entropic stabilizer for the lattice-Boltzmann method

.


I. INTRODUCTION
The classical numerical methods for the simulation of hydrodynamic fluid flows rely explicitly on the macroscopic target equations. An alternative, implicit or indirect way of computing approximate solutions for these target equations is provided by the lattice-Boltzmann (LB) method [1,2]. It is based on the kinetic theory and relies on mesoscopic, Boltzmann model equations. In the standard schemes, the dynamic variables, single-particle distribution functions, evolve via Lagrangian dynamics: evolution of the distribution functions mimic the alternating collision and propagation events experienced by the underlying particles. It can then be shown, using, e.g., a multiple-scale analysis, that the collective evolution of the distributions gives rise to macroscopic balance equations which conform with the original target equations.
Due to its kinetic origins, the standard LB schemes carry more degrees of freedom than strictly needed, e.g., for the approximation of solutions to the Navier-stokes equation. To put it simply, these LB schemes involve more distributions than there are independent, hydrodynamic variables present in the target equations. Interestingly enough, many important method developments have emerged when this fact is embraced instead of considering the spare degrees of freedom just as an extra burden on computational resources. Notable examples include LB models for fluctuating hydrodynamics [3][4][5], recoloring steps of the color-gradient based LB multicomponent models [6][7][8], and tuning of the relaxation parameters for the high-order * Corresponding author: keijo.mattila@jyu.fi moments in multiple-relaxation-time (MRT) as well as in two-relaxation-time (TRT) schemes [9][10][11].
Tuning the independent relaxation parameters often aims at improving stability of LB schemes. Other stabilization techniques for LB schemes include utilization of Ehrenfest's coarse-graining idea [12], addition of artificial dissipation [13], nonequilibrium entropy limiters [14], and application of selective viscosity filters [15]. A review of some stabilization techniques is provided in Ref. [16]. In particular, the entropic LB-BGK (Bhatnagar-Gross-Krook) model [17] maximizes the entropy by locally tuning the single relaxation time. It can be perceived as a kind of subgrid-scale model altering the local viscosity wherever the computational grid is underresolved [18,19].
Recently a new extension to LB schemes was proposed: namely, the entropic stabilizer [20]. Unlike the entropic LB-BGK model mentioned above, this extension does not locally alter the viscosity, but rather relies on modifying the relaxation time for the higher-order moments (i.e., the moments beyond the stress tensor) which do not contribute to the viscosity. In this respect, the extension is akin to the already discussed relaxation parameter tuning for MRT and TRT schemes. For example, in TRT schemes the relaxation time for the odd moments is, in principle, and for isothermal flows, a free parameter [21]. The specific aspect of the entropic stabilizer is that the common relaxation time for the higher-order moments is not fixed, but locally computed so as to maximize the entropy in the post-collisional state. That is, this relaxation time varies both in space and time; it self-adapts according to the maximum entropy condition which, at first glance, seems convenient from a simulation point of view.
In fact, the entropic stabilizer is very liberal: due to the proposed self-adaption mechanism, the relaxation time for the higher-order moments is not bounded nor explicitly controllable. Our aim here is to numerically investigate the ramifications of the extreme freedom granted for this relaxation time. We proceed by first reviewing the lattice-Boltzmann method and the entropic stabilizer in Sec. II. There we also present how the well-known regularized LB scheme [3,[22][23][24][25] can be extended with the entropic stabilizer. The numerical experiments based on the perturbed double periodic shear layer flow, and the related observations, are then explained in Sec. III. Finally, some conclusions are presented.

II. LATTICE-BOLTZMANN METHOD
In accordance with its kinetic roots, the LB method operates by evolving a single-particle distribution function f i (r,t) ≡ f (r,c i ,t). The distribution function describes the probability of finding a particle, traveling with a velocity c i , from a position r at a given instant t. The standard LB-BGK scheme, for example, follows the governing equation where f eq i is an equilibrium function, f For isothermal flows, the standard second-order equilibrium function is with the kinetic projectors w i denote the lattice-dependent weight coefficients, δ αβ is the Kronecker delta, and the Einstein summation convention is implied with the repeated indexes. The so-called thermal velocity c T = c r /a s , also the speed of sound for isothermal flows, depends on the reference velocity, c r = r/ t, where r is the lattice spacing, and a s is the lattice dependent scaling factor [26,27].
The local mass and momentum density, ρ and ρu, are obtained as the first moments of the distribution function: This simple scheme provides approximate solutions for the Navier-Stokes equation with the kinematic viscosity ν = c 2 T (τ − t/2) [26] and pressure p = c 2 T ρ (an ideal gas equation of state).
Finally, it is common to decompose, at least conceptually, the evolution equation into streaming and collision steps. The collision step is completely local and involves only the relaxation procedure: Then, the streaming step is merely a shift of the distribution values and does not involve any computation:

A. Regularized LB scheme
A simple modification of the LB-BGK scheme operates according to the equation where is computed at the pre-collisional state, i.e., from f in i . This update procedure has been proposed by several authors [3,[22][23][24][25], and is here referred to as the (original) regularized LB scheme.
The zeroth, first, and third moments of s neq i are identically zero while The equilibrium function in Eq. (2) is only a second-order Hermite polynomial expansion of the true Maxwell-Boltzmann distribution f MB . In other words, it is a projection of only the first three moments of f MB onto the kinetic space. Equations (3) and (4) specify an equivalent projection of the nonequilibrium moments (the zeroth and first moments are zero).
To summarize, Eq. (3) describes propagation of filtered information along the ith characteristic line, i.e., a projection of the full distribution onto the velocity space generated by the first three Hermite polynomials when the discrete velocity sets are of second order (like the D2Q9 and D3Q19 velocity sets). In particular, all the third-and high-order ghost moments are filtered out. This does not happen, e.g., with the original BGK model.

B. Entropic stabilizer
In order to improve stability of the standard LB-BGK scheme, Karlin et al. [20] proposed a modification of Eq. (1): the essence is to relax the nonequilibrium part related to the stress tensor using τ while the remaining part, related to higher-order nonequilibrium moments, is relaxed using a separate relaxation time.
To begin with, a suitable moment space is constructed by defining q linearly independent combinations of f i . Using the inverse transformation, the distribution functions can now be represented as f i = k i + s i + h i , where k i , s i , and h i depend only on the kinematic or conserved variables, the stress tensor, and the higher-order moments, respectively. Moment spaces for the D2Q9 LB model, based on both natural and central moments computed with the microscopic c i and peculiar (c i − u) velocities, respectively, are presented in the Supplemental Material of Ref. [20]. Below the corresponding schemes are referred to as the natural and central moment space LB schemes.
Then the LB equation is reformulated in a general way: where β = t/2τ and the dimensionless relaxation parameter γ is associated with the higher-order moments. This formulation reduces to the LB-BGK scheme when γ = 2.
A regularized LB scheme is obtained when γ = 1/β, not with γ = 1 as suggested previously [20]. However, in the low-viscosity case, where β ≈ 1, the latter choice is a good approximation of the correct value. Now comes the pivotal aspect of the entropic stabilizer: γ is computed locally and at every time step so as to maximize the entropy in the post-collisional state. An approximation for the optimal relaxation parameter γ is given in Ref. [20]: where is an entropic scalar product. We will utilize Eq. (6) in all our numerical experiments and, accordingly, consider γ ≡ γ . Karlin et al. compared the above approximation with the numerically solved optimal γ and found excellent agreement [20]. Hence, the use of this approximation has no influence on our conclusions. By substituting Eq. (6) into Eq. (5), we obtain meaning that the higher-order nonequilibrium part is first carried over to the mirror state, −h neq i , and then relaxed towards zero while φ scales the resulting value. Stability is typically compromised in low-viscosity cases, i.e. β ≈ 1, where the post-collisional value of the higher-order nonequilibrium part is well approximated with φ h neq i . Note that φ is not bounded and it cannot be controlled in the proposed setting.
Therefore, in contrast with the regularized scheme, the high-order part of the distribution, h i , responsible for the third-and higher-order ghost moments, is not filtered out but instead relaxed towards equilibrium with a parameter γ . This high-order contribution can be further decomposed into two parts, h i = g i + r i , where g i is related to the thirdand fourth-order moments, as presented in the Supplemental Material of Ref. [20], and r i to the higher-order moments that do not fit into the moment space of a second-order discrete velocity set. Since second-order sets are not able to retrieve the macroscopic equations without third-order errors (u 3 ) [27][28][29], both the regularized and entropic stabilizer schemes require simulations to be performed at low Mach number regime. The inclusion of third-and fourth-order moments into the moment space of a second-order velocity set can be beneficial to stability [30]. Nevertheless, these third-and fourth-order moments do not affect the macroscopic behavior of a LB scheme, up to the usual Navier-Stokes level of description, and hence they must be considered as ghost moments.

C. Regularized LB scheme with entropic stabilizer
The above presented entropic stabilizer can be considered as an add-on or extension to LB schemes. It is composed of three main ingredients: (1) contributions from high-order moments, the so-called nonhydrodynamic or ghost modes, are collected into h i , (2) the high-order part h i is relaxed with a specific relaxation coefficient, and (3) this relaxation coefficient is dynamically computed so as to maximize the entropy in the post-collisional state.
As an add-on, the entropic stabilizer can be plugged into various basic LB schemes. Here we extend the regularized LB scheme of Sec. II A with the stabilizer. In the following numerical investigation, results computed with this alternative LB scheme are compared with those obtained with the schemes of Sec. II B in an effort to better understand the properties of the entropic stabilizer.
First, the construction of a suitable moment space is an inconvenience of the schemes presented in Sec. II B. It has to be done separately for each discrete-velocity set; particular labor is required with sets including large number of velocities [27,28,31]. Instead, the Hermite polynomial based kinetic projectors could be utilized. Specifically, the equilibrium part in Eq.  (5) is then applied with these definitions.
The first three moments of the newly defined h neq i are zero. In addition, these definitions are not dependent on a particular discrete-velocity set. The presented procedure for defining h neq i can also be extended to thermal LB models. Recall that this schemes reduces to the original regularized LB scheme presented in Sec. II A when γ = 1/β.

III. NUMERICAL EXPERIMENTS
In order to numerically investigate the properties of the entropic stabilizer, we will set up a perturbed double periodic shear layer flow as a benchmark case. There a small velocity perturbation, perpendicular to the shear flow direction, initiates a Kelvin-Helmholtz instability, and causes roll up of the antiparallel shear layers. Periodic boundary conditions are applied. The initial velocity field is given by U 0 defines amplitude of the inital velocities, x and y refer to node index coordinates, while L is the number of nodes in both x and y directions. The parameters λ and control width of the shear layer and the initial velocity perturbation, respectively [32]. Here we consider only the thin layer case: we set λ = 80 and = 0.05. The Reynolds number is defined as Re = U 0 L/ν and the characteristic time T = L/U 0 . The D2Q9 LB model is used in all numerical experiments presented below, and time is given in dimensionless units, i.e., t * = t/T .

A. Initialization procedure
In order to fully specify f i at t = 0, the initial density ρ 0 and the nonequilibrium part remain to be defined. For the density initialization, the simplest choice is a constant value, e.g., ρ 0 = 1. However, as was pointed out by Skordos [33], a consistent initial density field is obtained by solving the related pressure-Poisson equation (for an incompressible fluid). Mei et al. [34] proposed that both the consistent initial density field and the nonequilibrium part can be obtained simultaneously within the lattice-Boltzmann framework via a specific initial iteration procedure. Here this iteration procedure is not practical because the relatively low viscosity in the benchmark simulations leads to very slow convergence. An additional difficulty with the entropic stabilizers is that during the initial iteration procedure, the complicated interplay between the evolving density field and the adaptive γ may lead to divergence. A possible solution is to fix γ in the initialization procedure, but here this option is not investigated further.
Instead, we choose to solve the related pressure-Poisson equation with a non-LB, iterative scheme based on the pseudounsteady method or the method of false transients [35,36]; the spatial derivatives are approximated with second-orderaccurate, isotropic finite differences (see, e.g., Ref. [37] and references therein). From the code structure point of view, a straightforward implementation of this iterative scheme is identical to a typical implementation of a LB scheme. Hence, incorporating it as an initialization routine into an existing LB code requires a relatively small amount of work.
Compared to the Mei et al. approach, using a non-LB, iterative scheme benefits from the fact that the pressure-Poisson equation itself does not depend on the viscosity: here the equation is solved only for each combination of L and U 0 . Numerical solution of the pressure-Poisson equation for L = 512 and U 0 = 1/32 is presented in Fig. 1. The nonequilibrium part is then defined separately. The simplest choice is to set f neq i = 0, i.e., distributions are initialized with the equilibrium part only. However, a more refined approach for initializing f neq i exists [33,38]. First, based on the Chapman-Enskog analysis, the viscous stress tensor is measured, using second-order-accurate, isotropic finite-differences [37].
For the LB-BGK and regularized LB schemes we then set respectively. For both schemes, initialization with a constant density, ρ 0 = 1, leads to high-amplitude, long-wavelength oscillations while the consistent initial density field results in smooth variations. The same observation was made for the Taylor-Green vortex flow in Ref. [34]. In addition, the central moment space LB scheme seems to be sensitive to the initialization of s  density field and nonequilibrium part, according to Eq. (9), results in smooth evolution of the system state as expected. Hence, such an initialization procedure is adopted for all simulations presented below.

B. Stability
We continue with a simple stability experiment and measure the decay of the average kinetic energy for L = 128 and Re = 3×10 4 . The results are presented in Fig. 6: the LB-BGK becomes unstable very early, around t * = 0.3, while the central moment space LB scheme (Sec. II B) with γ = 1/β experiences instability close to t * = 2.5. The natural and central moment space LB schemes (Sec. II B, γ = 1/β) remain stable. These observations agree with the previously reported results (see Fig. 2 in Ref. [20]).
Moreover, both regularized LB schemes are also free from stability problems. This in turn conforms with the published results, where the original regularized LB scheme (Sec. II A) is shown to significantly improve stability of the LB-BGK scheme [24,39]. The evolution of the average kinetic energy appears almost identical for the four stable schemes.
In fact, in this benchmark case, all three LB schemes utilizing the entropic stabilizer seem virtually indestructible in the sense that the simulated flow fields do not appear to diverge, i.e., to blow up, which is usually considered the tell-tale sign of instability. In support of this observation, a staggering improvement of four orders of magnitude in stability was reported (till at least Re ∼ 10 7 with L = 128) when the central moment space LB scheme (Sec. II B) was utilized instead of the LB-BGK scheme [20]. Such an improvement is of course astounding and warrants further investigation.
As was pointed out by Minion et al. [32], a global measure like the average kinetic energy above, does not capture well the differences between local flow features produced by the schemes under inspection. Here a hint of the differences is given in Fig. 7 where the simulated vorticity field, at t * = 1, is presented for L = 256,Re = 8×10 4 , and the six schemes considered; the LB-BKG case is simulated with Re = 3. hand, the fields appear undisturbed when the computation is carried out using the two LB schemes which filter out the high-order nonequilibrium moments, i.e., the original regularized LB scheme (Sec. II A) and the central moment space LB scheme (Sec. II B) with γ = 1/β.
The appearance of additional, spurious vortices has been previously investigated in Ref. [32]: the three LB schemes utilizing the entropic stabilizer would there be categorized as robust schemes implying that, while remaining stable, they can produce unphysical results whenever the computational grid is under-resolved.
Here the spurious vortices seem to be dependent on the initial velocity U 0 in particular. This would suggest that, with LB schemes, the vortices are related to compressibility issues as was proposed by Dellar [40]. Furthermore, the central moment space LB scheme (Sec. II B) is especially sensitive to the initial conditions: the scheme has a distinct dependence on the simulation parameters with each of the four initialization procedures discussed in Sec. III A and the observed spurious vortices are markedly different in all four cases.

C. Magnitude of γ
The peculiar variable related to the entropic stabilizer is the relaxation parameter γ for the high-order moments. Figure 8 presents its relative magnitude for the LB scheme based on the central moments space (Sec. II B) and for various Reynolds numbers. The result, with a maximum in |γ |/L ∼ 1, implies that the maximum |γ | is inversely proportional to the lattice spacing r. If the leading-order error terms of schemes utilizing the entropic stabilizer are proportional to γ r 2 , i.e., in the same way as the error is dependent on τ in the LB-BGK scheme [41], the result would immediately suggest first-order spatial accuracy.
However, the isolated maximum values do not necessarily determine the global behavior of a simulated system. But it is clear from Fig. 8 that the maximum values of γ are far beyond the linear stability range [0,2] even in well-resolved simulation cases. This is not so surprising. First of all, it is reasonable to assume that lowering the Reynolds number or, in particular, improving the grid resolution will reduce the nonequilibrium part of the distribution. Hence the value h neq | h neq in Eq. (7) can become small. At the same time, in the hydrodynamic regime with well-resolved gradients, the low-order moments, including the stress tensor, should be much greater than the higher-order, nonhydrodynamic moments [42,43]. This allows, in particular, s neq | h neq h neq | h neq , i.e., γ (φ) can become large.

D. Convergence of error
In order to measure the effect of very large values of γ on the numerical behavior of a scheme at the system scale, we will determine the convergence rates of error with lattice-refinement studies. First, separately for each scheme and Reynolds number, a reference solution for the velocity field, u R = √ u R · u R , is computed using L = 2048. Then a relative L 2 -error norm is computed: where u S is the velocity field simulated with a lower resolution, and the summations involve only those nodes which correspond to the nodes of the coarsest lattice (L = 64). Figure 9 plots the numerical errors for LB-BGK, the natural and central moment space LB schemes (Sec. II B), and for the regularized LB schemes when Re = 2×10 4 . The schemes share a common convergence rate in the high-resolution cases, i.e., L = 512 and L = 1024. The numerical errors for the original regularized scheme (Sec. II A) and LB-BGK continue  to converge with the same rate in the low-resolution cases; the LB-BGK becomes unstable at L = 128.
The natural and central moment space LB schemes, on the other hand, experience a steep increase in the numerical error when going from L = 512 to L = 256. Furthermore, at L = 128 the numerical error with natural moments is greater than with Central moments. The simulations on the coarsest lattice, L = 64, are under-resolved to such an extent that all schemes struggle to capture the proper flow dynamics. The measured base or reference convergence rate is ∼ r 2.36 .
The same story is repeated in Fig. 10 where the numerical errors are plotted for Re = 4×10 4 . However, now the error for the central moment space LB scheme (Sec. II B) deviates from the common value already at L = 512. Moreover, the natural moment space LB scheme (Sec. II B) experiences an even greater increase in the error when going from L = 512 to L = 128. In all cases, the errors of the natural moment space LB scheme (Sec. II B) and the regularized LB scheme (Sec. II C) are very close to each other. All these observations clearly indicate that the value of the relaxation parameter γ has an effect on the leading-order error term and, hence, influences the order of accuracy of a scheme in a complex way. At the same time, the original regularized LB scheme (Sec. II A) appears to behave consistently over the full range of simulation parameters experimented.

E. Discussion
In an attempt to find a correlation between the values of γ and the numerical behavior of a scheme at the system scale, we measured the average magnitude of the relaxation parameter for the higher-order moments. The measured values for the central moment space LB scheme (Sec. II B) are presented in Fig. 11. The first observation is that in well-resolved cases the average value of γ ≈ 1/β. That is, in these cases the entropic stabilizer filters out, on average, the high-order nonequilibrium moments similarly to the original regularized LB scheme. An analogous observation was already made above from the numerical error measurements.
From Fig. 11 it is clear that the average value grows when the resolution is decreased or the Reynolds number is increased. When further comparing the average relaxation parameter and numerical error measurements for the central moment space LB scheme (Sec. II B), it is observed that when the magnitude of the average value exceeds 0.02 or 0.03, anomalies in the numerical error begin to emerge in this particular simulation setup. Moreover, the proportion of distribution updates where the higher-order moments are amplified, i.e., |1 − βγ | > 1, is plotted in Fig. 12. After a comparison with the error measurements, it is observed that about 1% of such updates is tolerated; a larger number of amplification occurrences will be manifested adversely in the numerical error.
Numerical analysis of LB schemes utilizing the entropic stabilizer, by theoretical means, is intricate. The main challenge is to estimate the order of magnitude of γ . For example, the standard Chapman-Enskog (CE) multiple scale analysis requires an estimate for the magnitude of γ with respect to the small expansion parameter . With the assumption γ ∼ 0 , the usual Navier-Stokes level CE analysis does not reveal any contribution from the entropic stabilizer to the hydrodynamic approximation. However, the maximum values of γ observed in Sec. III C clearly compromise such an assumption. Based on the above presented results, estimating the magnitude of γ using observed maximum values is perhaps overly restrictive while the relevance of the average value as an order of magnitude estimate remains unclear. Hence, more theoretical work on the numerical analysis is required.

IV. CONCLUSIONS
We have reviewed the entropic stabilizer recently proposed for improving the stability of LB schemes. The entropic stabilizer can be considered as an add-on or extension to basic LB schemes. For example, here we extended the well-known Regularized LB scheme with the entropic stabilizer. The main advantage offered by this extended regularization scheme is that it does not require the construction of a moment space. This is convenient especially with LB models involving a large number of discrete velocity vectors.
Using the perturbed double periodic shear layer flow as a benchmark case, we numerically investigated the behavior of the entropic stabilizer. The investigation was carried out by comparing numerical results from total of six distinct LB schemes: LB-BGK, the original regularized LB scheme (Sec. II A), the natural and central moment space LB schemes (Sec. II B), the central moment space LB scheme (Sec. II B) with γ = 1/β, and the extended regularized LB scheme (Sec. II C).
It is verified that all three LB schemes which utilize the entropic stabilizer seem unconditionally stable in the sense that the simulated flow fields do not appear to diverge. Our main objective, however, was to inspect how the unbounded, and not explicitly controllable, relaxation parameter for the higherorder moments influences simulated local flow features and the numerical accuracy of a scheme. It was observed that the freely evolving relaxation parameter can have an undesirable effect on both, especially in the under-resolved simulation cases for which the stabilizer is targeted.
In particular, the entropic stabilizer can degrade the order of accuracy of a LB scheme, and promote unphysical flow features like spurious vortexes in perturbed double periodic shear layer flows. Furthermore, it was observed that the natural and central moment space LB schemes (Sec. II B) differ in their behavior: the central moment space LB scheme appears as the most sensitive to initialization, the natural moment space LB scheme (Sec. II B) and the extended regularized LB scheme (Sec. II C) behave alike, and the central moment space LB scheme (Sec. II B) with γ = 1/β is less stable than the original regularized LB scheme (Sec. II A). Moreover, the original regularized LB scheme (Sec. II A) performed consistently over the whole range of simulation parameters used.
In summary, the extreme stability provided by the entropic stabilizer can potentially deceive users into simulating fluid flows with improper configuration parameters. To avoid mixing features that are genuine with those that are numerical artifacts, careful convergence studies must be carried out per case. This can be inconvenient and expensive. Hence, in addition to systematic numerical validation, more detailed theoretical analysis of the entropic stabilizer is still required in order to clearly understand its effect on the numerical properties of a scheme.
On the other hand, the adaptive relaxation parameter for higher-order moments could perhaps be used as an indicator. For example, it could indicate the need for grid refinement or it could simply act as a gauge for the validity of the simulations. This is a prospect for future work. Another option would be to explicitly enforce the adaptive relaxation parameter into a limited range, say, the linear stability range. The behavior of the entropic stabilizer in this case could also be a topic for future work.