Does the multiresolution lattice Boltzmann method allow to deal with waves passing through mesh jumps?

We consider an adaptive multiresolution-based lattice Boltzmann scheme, which we have recently introduced and studied from the perspective of the error control and the theory of the equivalent equations. This numerical strategy leads to high compression rates, error control and its high accuracy has been explained on uniform and dynamically adaptive grids. However, one key issue with non-uniform meshes within the framework of lattice Boltzmann schemes is to properly handle acoustic waves passing through a level jump of the grid. It usually yields spurious effects, in particular reflected waves. In this paper, we propose a simple mono-dimensional test-case for the linear wave equation with a fixed adapted mesh characterized by a potentially large level jump. We investigate this configuration with our original strategy and prove that we can handle and control the amplitude of the reflected wave, which is of fourth order in the space step of the finest mesh. Numerical illustrations show that the proposed strategy outperforms the existing methods in the literature and allow to assess the ability of the method to handle the mesh jump properly.


Introduction
In [3,4] we have proposed a novel approach to perform time-dynamic grid adaptation and to devise lattice Boltzmann schemes to be deployed on such adaptive meshes using multiresolution analysis. This framework guarantees reduced memory footprints for problems with steep solution and potential gains in terms of computational time. This is obtained by storing the solution solely at the local level of refinement, which is essentially coarse in the case of problems involving shocks, see [6]. Moreover, the preeminent feature of this strategy is that it yields a precise bound on the additional numerical error caused by the compression of the solution on a dynamically adapted grid coupled with the numerical scheme, thanks to the information on the local smoothness of the solution provided by multiresolution. These peculiarities have been thoroughly investigated in these two works. The next question was to clarify -besides the previously elucidated error control -the potential alteration of the physics of the simulated system by the adaptive method, which has been fully studied in [2], concluding that our method of choice does not alter the reference lattice Boltzmann scheme up to local contributions of order four in the space step. This performance seems unprecedented in the literature. This analysis has been carried out on uniform or locally uniform grids for smooth solutions using the equivalent equations [7]. When non-smooth solutions are obtained, one must utilize the time-adaptive refinement to allow for the representation of singularities at the finest level of resolution.
The question which stimulated the present work has been raised by Pierre Lallemand and François Dubois during a seminary at the IHP in Paris and is the following: how does the scheme introduced in our previous works [2][3][4] behave once an acoustic wave is forced to pass through grids of different resolutions, like the one on Fig. 1? We emphasize that in the context of dynamic mesh adaptation, the only case in which waves could penetrate through a level jump is when they are locally very smooth. For any numerical solver dealing with adapted meshes, beyond lattice Boltzmann schemes, it is known that a wave passing through a grid transition normally splits into two parts: the first one propagating through the interface and the second one which is reflected back. The phenomenological reason for this is the different acoustic impedance of two media made up of grids at different levels of refinement. Generally speaking, the aim is to devise numerical approaches which minimize the amplitude of the reflected waves. This is particularly important in applications such as aeroacoustics [9], where spurious currents on the density field can have important consequences on the solution of the problem.
Many strategies, both for the representation of an adaptive grid and for the construction of adaptive lattice Boltzmann methods have been explored. For a review on the different techniques and the issues linked with grid transitions, the reader can consult [11] and [10]. There is no widespread consensus on a standard test case to analyze the issue of reflected waves: we therefore consider a basic one-dimensional configuration which already introduces the main difficulties of more realistic multidimensional systems (e.g. the D2Q9 scheme) to be found in the available literature, arising from the applications. We also tested 1 our approach against the two-dimensional acoustic pulse test case from [1,9] yielding results -not presented in this note -fully compatible with the simpler analysis introduced here.

Numerical method
Let us briefly sketch the structure of our numerical method in order to make the paper selfcontained. The interested reader can consult [2][3][4] for more precision.

Space and time discretization
For the sake of presentation, we consider a one-dimensional setting on the bounded domain Ω = [0, 3]. As in our previous work [2] and references therein, we take cells C ,k = 2 − k, 2 − (k + 1) for levels of resolution ∈ { min , . . . , max } and k ∈ {0, . . . , 3 × 2 ∆ − 1}, where ∆ = max − . The cell center is denoted by x ,k . Therefore, the space step of the grid at finest level of resolution is ∆x = 2 − max . Cells with different levels can be employed to cover the domain Ω of interest as in Fig. 1, which also illustrates the experimental configuration used in the sequel. The time discretization is uniform with step ∆t = ∆x/λ, where λ > 0 is the so-called "lattice velocity". The time-step is the same for all the levels ∈ { min , . . . , max } and determined by max . In many approaches available in the literature, level-wise time steps are employed. However, local time stepping would require a lot of attention and computational effort to preserve error control, thus preventing from a significant cost reduction (see [4], [5] and [12]).

Lattice Boltzmann method
The lattice Boltzmann method is a numerical algorithm based on q ∈ N velocities (e α ) α=q−1 α=0 ⊂ λZ compatible with the lattice, with dimensionless counterparts c α := e α /λ ∈ Z, α = 0, . . . , q − 1. The distribution function of the population moving with velocity e α shall be denoted f α . The algorithm is made up of two phases at each time step.

Collision phase
For every cell C ,k belonging to the mesh, the collision phase is performed locally and is a diagonal relaxation in the space of the moments where M ∈ GL(q, R), S ∈ M q (R) diagonal with rank(S) = q −q c where the non-zero 2 entries belong to ]0, 2]. Here, q c is the number of conserved moments and the vector of moments at equilibrium m eq is a function of these conserved moments.

Transport phase
For every cell C ,k belonging to the mesh, the transport phase does not mix different populations. For every α ∈ {0, . . . , q − 1} it comes under the form with weights recursively defined by The transport phase (1) needs information stored on ghost cells (dashed in Fig. 1), which thus must be updated. This is done by averaging for the ghost cells underneath more refined cells followed by the use of interpolations (from which (1) is recovered, see [2] for the details) for ghost cells above coarser cells. It should be observed that the proposed method exactly conserves the moments which are conserved by the original lattice Boltzmann scheme, as remarked in [3,4]. Moreover, we stress that the aim of our numerical strategy is not to change the features of the original (on the uniform mesh) lattice Boltzmann scheme, with its strengths and weaknesses, but rather to minimize its perturbation once introducing non-uniform meshes.

Target problem, results and discussion
We introduce a simple system which sustains travelling acoustic waves which are eventually sent against a level jump in the computational grid.

Equations and numerical scheme
The target equation is the linear wave equation with velocity c > 0 (taken c = 1/2 in the experiments, but any other value c < 1 would be fine) over the whole real line R which has been recast under the form of first order system for simulating it. For this test we simulate on the bounded domain Ω = Ω left ∪ Ω right with Ω left = [0, 2] and Ω right = [2, 3] with the initial datum u 0 (x) = exp(−100(x − 3/2) 2 ). The most bare-bone lattice Boltzmann scheme to handle such equation -yet yielding the difficulties of more sophisticated ones -is a three velocity scheme q = 3 with two conserved moments (u = m 0 and v = m 1 ) We consider λ = 1, p = 1.7 and a final time T = 1.5625 for each simulation. The initial data are taken at equilibrium.

Results
To quantify the amplitude of the reflected wave, we simulate on different configurations: • u jump is the solution obtained using the spatial discretization illustrated on Fig. 1, the left sub-domain Ω left is finely meshed with level left = max , whereas the right subdomain Ω right is coarsely meshed with level right = min ≤ left . We vary the jump gap jump := | left − right |. This configuration is the monodimensional equivalent of that presented by [11].
• u ref is the reference solution obtained on the uniform mesh at the finest level max .
• u coarse is the solution obtained on the uniform mesh at the coarsest level min .
We introduce the L 1 -norm of the exact solution u for the sake of normalizing: Here the double hat operator represents the application of the same interpolation used to update the ghost cells. This step is needed to compare solution on the same mesh. Using the analysis of [2], we can prove an estimate on the amplitude of the reflected wave D jump-refl (T ): Thus by computing the moment u and taking the difference, we deduce that for every time The CFL condition on the finest resolution imposes that information, thus the errors, propagate of one cell for one time step. The constant C > 0 carries the normalization The results of the numerical simulations are provided in Table 1. We have also run comparisons with the method from Fakhari and Lee [8] and that of Rohde et al. [13] presented in Fig. 2.

Discussion
Commenting on Table 1, we see that the reference scheme converges linearly E ref (T ) = O (∆x) once refining as expected from the analysis by the equivalent equations [7]. The error E coarse (T ) ≤  [2]. The same behavior is observed for the mesh with jump, namely for E jump (T ) and D jump (T ). Very interestingly D jump (T ) ≤ D coarse (T ): counter-intuitively this is a priori not granted due to the possible formation of waves reflected at the jump, even though only a part of the domain is coarsened. This gives a first indication about the fact that the reflected waves are perfectly mastered. The second indication comes from D jump-refl (T ) = O (∆x 4 ). This means that with our method, we are able to decrease the amplitude of the reflected waves with fourth-order convergence in the space step, in accordance with Prop. 1. The supra-convergence compared to D jump (T ) comes from the fact that at each time step, the reflected wave is generated only on the cell of Ω right next to the interface, so that it eventually propagates to the left inside the fine medium without additional amplification of the error. Observe that the convergence rates worsen for large ∆ and for small max due to the fact that we are no longer allowed to perform the Taylor expansions needed by Proposition 1, which are done at the current level of resolution . Indeed, in this case, one can no longer claim that 2 ∆ ∆x is O (∆x).
Compared to other methods, we can show with the same proof path than Prop. 1 that the Lax-Wendroff strategy by [8] yields D jump-refl (T ) = O (∆x 3 ), which is one order less than our method. This can also be qualitatively seen on Fig. 2. Concerning the approach by [13] where local timestepping is used, we observe that it yields quite large reflected waves. This waves are one order of magnitude larger than for the method by [8] and two orders of magnitude larger than our approach. However, the local time-stepping prevents us from applying the same theoretical study to this scheme. Mastering reflected waves at a high order of accuracy is important when our technique is extended to typical multidimensional applications. When simulating the incompressible Navier-Stokes equations via a quasi-incompressible D2Q9 scheme, spurious acoustic waves are of order O (∆x 2 ), thus controlling their reflection at order O (∆x 4 ) is a highly desirable feature of the scheme.   Figure 2. Results of the simulation on the mesh with jump (whole domain on the left, magnification on [1.6, 1.9] and on the y axis on the right). Initial solution in pale orange and solution at t = T in green (left subdomain) and blue (right subdomain). On the first row, we use our multiresolution scheme [2][3][4]. On the second row, we use the Lax-Wendroff scheme by Fakhari and Lee [8]. On the third row, the scheme with local time-stepping by Rohde et al. [13]. The simulation uses left = 10 and jump = 3.

Conclusions
In this contribution, we have briefly presented our adaptive lattice Boltzmann method, which is studied in Prop. 1 using the technique introduced in [2] to conclude that in case of a fixed mesh jump, the amplitude of the spuriously reflected waves is of order O (∆x 4 ). This fact is numerically verified and compared to the performance of other approaches available in the literature [8] and [13], showing that our method outperforms these traditional approaches. It is worthwhile observing that the original method [3,4] was conceived to be used with dynamically adapted meshes which automatically follow waves and fronts with finer discretizations once their lack of regularity justifies the depart from a coarse uniform mesh. Thus, in this case, we even do not expect the O (∆x 4 ) perturbation because fronts never cross level jumps but are precisely and successfully "chased" by the fine discretization. We observe that multiresolution is not relevant for systems developing homogeneous isotropic turbulence. However, dealing with spatially large problems where turbulent flows at high Reynolds number are present in a small portion of the domain could be interesting and advantageous. The analysis of this framework could be tackled in publications to come with a tailored data structure.
For the sake of a quick and effective presentation, we have restrained the study to the onedimensional setting. However, the generalization to higher spatial dimensions is straightforward and follows the indications of [2], since the operators involved in the multiresolution analysis are extended [3] by tensor product in the remaining directions of the space. We have proved in [2] that the multidimensional extension of the method retains the same accuracy levels in every spatial direction.