Particles on Demand method: theoretical analysis, simplification techniques and model extensions

The Particles on Demand method [B. Dorschner, F. B\"{o}sch and I. V. Karlin, {\it Phys. Rev. Lett.} {\bf 121}, 130602 (2018)] was recently formulated with a conservative finite volume discretization and validated against challenging benchmarks. In this work, we rigorously analyze the properties of the reference frame transformation and its implications on the accuracy of the model. Based on these considerations, we propose strategies to boost the efficiency of the scheme and to reduce the computational cost. Additionally, we generalize the model such that it includes a tunable Prandlt number via quasi-equilibrium relaxation. Finally, we adapt concepts from the multi-scale semi-Lagrangian lattice Boltzmann formulation to the proposed framework, further improving the potential and the operating range of the kinetic model. Numerical simulations of high Mach compressible flows demonstrate excellent accuracy and stability of the model over a wide range of conditions.


I. INTRODUCTION
The understanding of the nature of high-speed compressible flows has been a long sought goal in the scientific and engineering community.An accurate prediction of complex hydrodynamic features is crucial in modern research, as well as in technology, with examples such as the interpretation of astrophysical jets, captured in the images of deep space telescopes [1,2] and the design of air-frames and propulsion systems of high-Mach low-altitude flying vehicles [3].Throughout the history of computational fluid dynamics (CFD), a number of numerical approaches has been suggested for the simulation of high-speed flows, including artificial viscosity methods [4], total variation diminishing (TVD) [5], essentially non-oscillatory (ENO) schemes [6,7] and weighted ENO (WENO) schemes [8,9].The challenging nature of these flows renders the field an active research area [10,11], with developments such as positivity preserving limiters and targeted ENO (TENO) schemes [12], extending the domain of CFD towards even more exotic hydrodynamics [13][14][15].
In contrast to conventional CFD, the lattice Boltzmann method (LBM) addresses the evolution of hydrodynamic fields through the dynamics of a fully discrete kinetic system of designed particles associated with the discrete velocities c i , i = 0, . . ., Q − 1.The state is described in terms of the populations f i (x, t), which evolve in time and space by a simple algorithm "stream along links c i and collide at the nodes x in discrete time t".LBM has evolved into a versatile tool for the simulation of complex flows including transitional flows [16], flows in complex moving geometries [17], thermal and convective flows [18][19][20], multiphase and multicomponent flows [21][22][23][24], reactive flows [25] and rarefied gas [26], to mention a few recent instances; see [27][28][29] for a discussion of LBM and its application areas.However, despite the high efficiency and low numerical dissipation of LBM for nearly incompressible flows, the domain of high-speed compressible flows presents a number of severe challenges [18,[30][31][32][33].The main directions to extend conventional LBM towards the compressible realm includes standard lattices LBM augmented with correction terms [34][35][36][37], multi-speed lattices [38][39][40][41][42] and hybrid approaches [43][44][45][46].
A common feature of the conventional LBM is the propagation of the populations with fixed discrete velocities, which translates as fixing the reference frame "at rest".It is well known that, when the fluid velocity significantly deviates from the frame velocity, errors and numerical instabilities corrupt the solution, impeding the applicability of LBM to high-Mach flows [30,34,47].A remedy was the introduction of uniformly shifted lattices, which amounts to a constant shift of the reference frame, at every grid point of the numerical domain [48].The concept demonstrated excellent performance for predominately unidirectional compressible flows, shifting the operational domain of the method in par with the chosen reference velocity shift [35,42,48].While the concept of the uniform frame shift maintains key advantages of the scheme, such as simplicity and exact propagation, its potential diminishes for flows exhibiting large variations in flow velocity and temperature, due to the inevitable presence of strong deviations between the velocity of the actual flow and the imposed reference frame.
In contrast with the conventional LBM, the recently proposed Particles on Demand (PonD) method reformulates the kinetic equations in a space-time adaptive reference frame, dictated by the actual local fluid velocity and temperature [49].Two key elements were introduced with the PonD method: Firstly, PonD uses a consistent representation of populations in different reference frames, an operation termed as reference frame transformation.Secondly, a predictor-corrector iteration loop was applied, which ensured the realization of the propagation and collision step in the local co-moving reference frame, thereby optimizing accuracy and stability.Early realizations of PonD employed a semi-Lagrangian discretization, providing off-lattice flexibility to accommodate a varying reference frame, and validated the central concepts with a series of benchmarks, including multiphase and rarefied flows [50][51][52][53][54][55][56].However, the semi-Lagrangian method is prone to errors in conservation of mass, momentum and energy, deteriorating the accuracy of the solution in the presence of discontinuities (shock waves) [57].As a remedy to these shortcomings, a finite volume formulation of PonD was proposed in [57], following the discretization of the discrete unified gas kinetic scheme (DUGKS) [58][59][60].The resulting conservative scheme, combined with a reference frame transformation based on Grad's projection of particles populations, demonstrated excellent performance in an array of challenging hypersonic compressible benchmark flows, including extreme hydrodynamic features such as the formation of near-vacuum regions.
In this paper, we aim at a further development of the finite-volume formulation of PonD, targeting strategies that simplify the scheme and enhance efficiency.A detailed analysis of the solution methodology is presented, along with the requirements to be met by the reference frame transformation.The scheme is extended to include a forcing term, as well as a variable Prandtl number.Finally, we combine the idea of the multiscale framework suggested in [54] , with the Grad's projection frame transformation.The theoretical findings are validated in a series of numerical experiments along with extensive benchmarking of the scheme in challenging hydrodynamic flows.
The paper is organized as follows.The formulation of the kinetic equations in an adaptive reference frame is laid out in detail in Sec.II.Sec.III presents the kinetic model, which allows for variable adiabatic exponent and Prandlt number.Subsequently, Sec.IV describes the numerical discretization of the model.The model is extensively benchmarked in Sec.V, along with demonstration of important notions of the reference frame transformation.Finally, concluding remarks are provided in Sec.VI.

II. ADAPTIVE REFERENCE FRAME FORMULATION A. Discrete velocities
Without a loss of generality, we consider discrete speeds in two dimensions formed by tensor products of roots of Hermite polynomials c iα , The model is characterized by the lattice temperature T L and the weights W i associated with the vectors (1), where w iα are weights of the Gauss-Hermite quadrature.
The discrete velocities and the associated weights are shown in Table I.With the discrete speeds (1), the particles' velocities v λ ref i are defined relative to a reference frame λ ref , specified by the frame velocity u ref and the frame temperature T ref , The optimal reference frame is the comoving reference frame, which is specified by the local temperature T ref = T (x, t) and the local flow velocity u ref = u(x, t).

B. Reference frame transformation
A critical element of our construction is the transformation of the populations f λ i , defined with respect to a λ reference frame, to a different reference frame λ , λ = {u , T }. ( In this work, we follow the strategy of [57].Let us denote M λ k a moment tensor of order k, The reference frame transformation is then defined by the condition of invariance of the moments of orders k = 0, 1, . . ., K, where K denotes the maximal moment order which is required to be frame invariant.The transformed populations are then sought as a Grad's projection, where H (n) (c i ) correspond to the Hermite polynomials of the lattice velocities and the expansion coefficients α (n) (m; λ ) are calculated such that the moment invariant system (7) is satisfied (detailed in Appendix A).The latter depend on the vector of frame invariant moments m = {M 0 , . . ., M K } and the target reference frame λ .As a shorthand notation for the reference frame transformation, we use the following formula,

C. Solution methodology
We consider a simple kinetic model, which recovers compressible hydrodynamics under the restriction of fixed adiabatic exponent and Prandtl number.The kinetic evolution can be formulated in an arbitrary constant, four-parametric, reference frame λ, where Ω λ f,i is a collision kernel of the populations.The populations f λ(x,t) i , are described with respect to a local reference frame λ(x, t), which generally differs from the monitoring frame λ.The reference frame transformation connects the populations between these two reference frames, To simplify the notation, we drop the space-time dependence of the reference frame and λ is reserved for the local reference frame, λ = λ(x, t).The overbar and subscripts shall be used to denote monitoring reference frames, uniform throughout the domain (e.g.λ).Inserting the transformed populations (11) into the evolution (10) recovers the final equation of PonD, We emphasize that the kinetic equation ( 12) is formulated with respect to a uniform frame λ.Therefore, Eq. ( 12) constitutes a typical kinetic equation, with constant characteristics, amenable to usual numerical realizations in the context of LBM, such as integration along characteristics.The necessary element is the introduction of the frame transformation operator.The effect of the varying reference frame is evident in the operation G λ i,λ f λ , inside the non-local gradient operations.The above observation is crucial and determines the requirements that must be satisfied by the reference frame transformation.
Comments are in order: • The kinetic equations can be formulated in principle with respect to any arbitrary reference frame.As such, different monitoring points can employ different reference frames.The consistency of the evolution in different reference frames is established by proper reference frame transformations.In the limit of infinite discrete velocities, the solution in every frame would be identical, i.e. no reason to do that.But in discrete systems, the accuracy of the solution depends on the proximity of the imposed frame with the actual local frame, dictated by the local flow conditions.Thus, with this procedure, we maximize the accuracy of a given model across the domain.
• The direct formulation of a kinetic equation with adaptive velocities leads to additional "forcing" terms, containing derivatives with respect to the particles velocities.A thorough discussion in this direction can be found in [61].In the context of PonD, the solution methodology consists of a set of equations, each one at its own, spatially uniform, reference frame.This strategy avoids the explicit requirement for the computation of "forcing" terms.The price to be paid instead, amounts to the operation of reference frame transformations.
• Which domain needs to be transformed around a "frame-generating" monitoring point?Thanks to the hyperbolicity of the system, only the numerical domain of dependence needs to be transformed.Of course, this procedure for elliptic type of equations would be prohibitively computationally demanding.

D. Hydrodynamic limit analysis
We analyze the governing kinetic equation ( 12) with the Chapman-Enskog method and investigate the consistency requirements for the moment invariant system.We rewrite Eq. ( 12) in terms of a Bhatnagar-Gross-Krook (BGK) collision operator and a small parameter for the relaxation time τ 1 , Following the conventional notation, we introduce the following multiscale expansion, We inject the expansions into the governing equations and separate the dynamics according to different orders of , O( 1 ) : t (G λ i,λ f λ, (1) ) (2) . ( At the O( 0 ) order we obtain the equilibrium populations, which implies the following solvability constraints,

Equilibrium moments
The functional form of the equilibrium moments is the basic element of the analysis and determines the recovered hydrodynamic equations.We underline that all the lattices discussed in this work and listed in Table I, reproduce the pertinent equilibrium moments as their Maxwell-Boltzmann (MB) continuous counterparts in the comoving reference frame.For example, even for the standard D2Q9 lattice, the evaluation with the comoving reference frame frame retrieves the following moments, where the MB moments are, Here, overline denotes symmetrization.While all pertinent equilibrium moments in the comoving reference frame are accurate, the same conclusion does not necessarily hold when a different, non comoving, arbitrary reference frame λ = {u, T } is used for the evaluation.Indeed, the crucial difference between various lattices rests with the frame invariance of the equilibrium moments.
In the non comoving reference frame λ, the equilibrium populations f eq,λ are no longer given by the simple expression ρW i , but must be computed.This can easily be accomplished by the reference frame transformation, operating on the vector of the lattice weights W , from the comoving frame λ to the λ frame.By construction, the equilibrium moments in the λ frame match the MB moments, if they are frame invariant.For example, the following relation holds for a third-order Grad's projection, sustained by the D2Q16 lattice, However, deviations occur for higher order moments, not included in the set of frame invariant moments.Continuing with the same example, the fourth-order equilibrium moment in the λ frame becomes, The explicit form of the deviation can be computed via algebraic manipulations of Eqs. ( 33)- (34), where, As indicated by the previous expressions, the deviations vanish when the monitoring reference frame λ approach the comoving reference frame λ.

Full invariant moment system
Let us consider the case where the isotropy of the lattice supports the frame invariance of moments, up to fourth-order.This case corresponds, in particular, to the D2Q25 lattice mentioned in Table I.The zeroth order moment evaluation of Eq. ( 17) leads to the following equation, ∂ t ρ eq,λ + ∇ • J eq,λ = 0, (39) where, Both moments belong to the frame invariant system of the reference frame transformation.As such, they can be evaluated equally well at the comoving reference frame, J eq,λ = J eq,λ = ρu.
We substitute Eqs. ( 42), (43) into Eq.( 39) and recover the continuity equation, ∂ The same reasoning applies to the rest of the conserved moments of Eq. ( 17), since all pertinent moments are frame invariant.The momentum and energy conservation laws at the Euler level are as follows, ∂

∂
(2) where The compressible NSF equations are recovered from the summation of the O( 1 ), O( 2 ) contributions, where π is the pressure tensor, S the strain rate tensor, q is the heat flux, and 3. Third-order invariant moment system We continue with the analysis of a third order moment invariant system, which corresponds to the D2Q16 lattice.The Euler level dynamics and the NSF density and momentum contributions include moments which are frame invariant.Thus, Eqs. ( 45)-( 46) and the NSF density and momentum contributions ( 47)-( 48) are obtained accurately.However, the NSF energy contribution includes the flux of energy flux tensor, which is not frame invariant.The evaluation of this moment at the monitoring frame λ, instead of the comoving frame λ, induces an error, as seen from Eq. ( 37).This deviation gives rise to a diffusive error term at the NSF energy equation (49), As expected from this analysis, the D2Q16 lattice demonstrates excellent performance in inviscid Euler gas dynamic systems, even at the presence of very strong discontinuities.In the presence of important viscous effects, the error term (last term in Eq.( 60)) will affect the accuracy of the solution.The magnitude of the error scales with the spatial variation of the reference frame, or in other words the gradients of the velocity and temperature field.As shown in subsequent numerical simulations, benchmark viscous hydrodynamic flows can be accurately captured with D2Q16 lattice, suggesting that the magnitude of the error term is rather weak.However, as the velocity and temperature gradients grow, the error terms manifests in the solution.Sec.V B provides further discussion on this topic, with the aid of numerical simulations.

Second-order invariant moment system
Finally, we examine the hydrodynamic properties of a second order frame invariant moment system, with a representative example being the D2Q9 lattice.Such a model cannot support the full energy flux tensor and higher order tensors (Q, R) as frame invariant moments.Following the same reasoning with the previous cases, we observe that error terms are introduced in the energy equation of the Euler-level dynamics (46) and the momentum and energy equations of the NSF-level dynamics ( 48), (49).For relatively smooth flows without shocks, numerical evidence suggest that the effect of the error terms of the D2Q9 model is rather small.A prominent example is the case of an advected vortex, which has been shown to be captured accurately even for vortex and advection Mach numbers at the range of Ma v = 0.8 and Ma a = 100 [49].However, for hydrodynamic flows with shocks, the errors due to the frame variation are nonnegligible.Importantly, the shock dynamics at the Euler level are not described accurately, which translate into errors in the shock propagation speed.

Summary of observations
We summarize the domain of validity of the kinetic model, according to the order of the frame invariant moment system: • Second-order (D2Q9): Appropriate for smooth regions of compressible flows.
We take advantage of the above hierarchy to achieve the best efficiency with our framework.In particular, the low-order model can be applied to the smooth regions of the flow and the high-order model in the regions of steep hydrodynamic gradients.In the following section, we discuss the coupling of the different models, in the spirit of [54].

E. Multiscale frame transformation
The Grad's projection approach for the reference frame transformation is advantageous, in terms of stability and efficiency.In the following, we demonstrate an additional benefit, which is the deployment of different lattices throughout the domain, with minimal change in the framework and limited computational overhead.In essence, we combine the core idea of the multiscale concept [54] along with the Grad projection frame transformation.
Let us consider two velocity sets of different order, where q < Q.We distinguish two different operations coupled with the frame transformation from λ to λ : • Lifting: The lifting operation switches from the lower-order q-model to the higher-order Q-model, requiring thus a map, • Projection: The projection operation switches from the higher-order Q-model to the lower-order qmodel, requiring thus a map, The construction of both operations amounts to identifying the proper expansion coefficients of the Grad's expansion.We recall that the expansion coefficients are function of the moments and the reference frame λ .

Lifting
For the lifting operation, we can identify the list of moments m q→Q , required for the reference frame transformation , as a composition, where m q is operationally available from f λ q and m Q−q constitutes the remaining unknown higher order moments.The lifting operation consists in specifying the respective contributions as, where M q,λ is the q × q matrix of the populations to moments map.With the required moments identified, the construction of the lifted populations proceeds similarly with Sec.II B. The expansion coefficients are computed from the moments m q→Q and the reference frame λ , The lifted populations f λ Q can then be found from Grad's expansion, where and K are the weights, the Hermite polynomials and the order of expansion of the high-order model respectively.

Projection
In the projection step, the high-order population f λ Q contains the subset of the q linearly independent moments, which is required for the construction of the loworder population f λ q .Hence, in contrast to the lifting procedure, there is no missing information and the loworder moment vector m Q→q is operationally available Similarly to the lifting operation, the projected populations f λ q are given by the Grad's expansion, where W q i , H (n) (c q i ) and k are the weights, the Hermite polynomials and the order of expansion of the low-order model respectively.

III. VARIABLE ADIABATIC EXPONENT AND PRANDTL NUMBER A. Kinetic model
The kinetic model can be extended towards a variable adiabatic exponent via the two-population approach [42].The second set of populations (g-populations) is designed to carry the internal energy associated with non-translational degrees of freedom, and thus enable an adjustable adiabatic exponent γ = C p /C v , where C p = C v + 1 is the specific heat of ideal gas at constant pressure and C v is the specific heat at constant volume [62,63].The governing kinetic equations can be written as follows, Additionally, the collision operators can accommodate an intermediate relaxation to quasi-equilibrium states, thus enabling a variable Prandtl number [42,64], where f * i , g * i are the quasi-equilibria of the f − and g− populations and the relaxation time τ 2 determines the Prandtl number.The local conservation laws for the density ρ, momentum ρu and the total energy ρE are, where the total energy of ideal gas is, The equilibrium populations in the comoving reference frame are as follows, The expressions for the dynamic viscosity, bulk viscosity and thermal conductivity are [42], The Prandtl number is therefore, For Pr < 1, the quasi-equilibria are designed to conserve the centered heat flux, resulting in the following expressions, where e i = v i − u, Q is the non-equilibrium third-order flux tensor and ς is the energy flux associated with the internal degrees of freedom, e i (g i − g eq i ). (87)

B. Comments on g− populations
We note that the concepts presented so far apply equally well for the g− populations, with the sole difference being the required frame invariant moments, which have to be supported by the corresponding g− lattice.A Chapman-Enskog analysis [42] shows that g− equilibrium moments up to second order are enough to recover the NSF equations.Therefore, the D2Q9 is safely employed in this work for the g− populations.The reference frame transformation (Sec.II B) and its multiscale realization (Sec.II E) apply equally well for the g− populations, taking into account that the maximal frame invariant moment is second order.

IV. NUMERICAL IMPLEMENTATION A. Finite volume discretization
We proceed with the finite-volume discretization, in the spirit of PonD-DUGKS framework [57,58].In accord with the notions above, the kinetic equation can be formulated in an arbitrary reference frame λ.We first present the discretization for Pr = 1.The extension for variable Prandtl number is explained in the following section.

Updating rule
The evolution of the kinetic model ( 70)-( 71) can be discretized as follows, The update equations are derived from the integration of the continuous equations ( 70)-( 71), formulated in the reference frame λ, in a control volume centered at x j , with volume V j , from time t n to t n+1 = t n + δt, using the midpoint rule for the convection term and the trapezoidal rule for the collision term [58].To remove the implicitness, DUGKS scheme adopts the variable transformation from the standard LBM practice, [65,66] where φ stands for the f -and g-populations and Ω φ,i are the collision BGK kernels defined in Eqs. ( 70)- (71).The fluxes of the populations F λ φ,i (x j , t n+1/2 ) across the surface of the control volume are defined as, where n is the outward unit vector normal to the surface.Finally, we remark that within the finite volume context, the populations and the collision terms are cell-averaged quantities, The reference frame which is used for the evolution of the populations at (x j , t n ), is set to the comoving frame, from the known flow velocity and temperature,

Flux evaluation
The key element of the update equations ( 88),( 89) is the evaluation of the flux term, F λ φ,i (x j , t n+1/2 ), which contains the unknown populations φ λ i (x b , t n+1/2 ) at the cell interface x b and time t n+1/2 .The frame which shall be used for the flux evaluation is λ F = {u F , T F }, with the frame velocity u F and temperature T F constructed by the average frame of the adjacent cell centers to the interface x 1 , x 2 , The integration of Eqs. ( 70)-( 71) along the characteristics for half-time step shows that the required populations φ λ F i (x b , t n+1/2 ), are connected with the known populations at time t n through the following equation [58], where, Eq. ( 96) is essentially a half-time semi-Lagrangian step, with the final point located at the interface x b , at t n+1/2 .The populations φ+,λ F i and the spatial gradients σ λ F i = ∇ φ+,λ F i are subsequently evaluated in the neighbouring cells of the interface, at time t n .In this work, Van Leer and minmod slope limiters were used for the computation of the spatial derivatives [67,68].We also note that the reference frame transformation is applied, to express the required populations from their original reference frame to the target reference frame λ F .The populations are reconstructed at the departure point x = x b − v λ F i δt/2, with the MUSCL scheme [69], According to Eq. ( 96), we obtain the φλ F i populations at the interface x b and time The density, momentum and temperature at (x b , t n+1/2 ) are finally computed by With the calculated macroscopic fields (ρ, u, T ) at (x b , t n+1/2 ), the equilibrium populations φ eq,λF i (ρ, u, T ) can be computed and subsequently also the populations φ λ F i (x b , t n+1/2 ), after inversion of Eq. ( 97).We remind that the equilibrium populations can be obtained through the reference frame transformation (34), Finally, the fluxes which are required to update the cell centers populations can be found from summation over the faces of the cell and proper reference frame transformation, where x b,c designates the center of the c-th face of the cell, n c is the outwards normal vector and λ is the reference frame of the evolution of the cell.

Summary of the algorithm
Based on the previous steps, we summarize the evolution procedure from time t n to t n+1 : 1. Initial data (cell centers x j ) • Given (ρ, u, T ) , comoving reference frame λ = {u(x j , t n ), T (x j , t n )} and populations.φ λ i (x j , t n ) • Calculation of the φ+,λ i (x j , t n ) populations, according to eq. (98).

Calculation of the fluxes (Loop over cell faces x b )
• Set reference frame λ F at interface and time (x b , t n+1/2 ) • Calculation of the populations φ λ F i (x b , t n+1/2 ) according to procedure in Sec.IV A 2.

Population update (Loop over cell centers x j )
• Computation of the fluxes to the local reference frame of the cell, Eq. ( 105) and update the populations through Eqs.(88), (89).
We stress the crucial difference between the proposed realization and the scheme suggested in [57], which is the absence of iterations within the flux evaluation step.
We remind that a semi-Lagrangian step is executed to retrieve the populations at the cell faces, according to Eq. ( 96).In this work, the reference frame for the above step is set from the average reference frames of the neighbouring cell centers (eq.( 94)) and the flux evaluation is performed explicitly.With this approach, it is necessary to obtain non-comoving equilibrium populations, φ eq,λF (ρ, u, T ), to finalize the flux evaluation.The scheme in [57] suggested an iterative predictor-corrector procedure, such that the flux calculation is realized in the comoving reference frame.While computationally demanding, the iteration procedure operates only with the simple comoving equilibrium populations Eqs. ( 78), (79).A further analysis of this aspect via numerical simulations is provided in Sec.V D.

B. Imlementation of variable Prandtl number
The quasi-equilibrium relaxation can be implemented as a forcing term in the kinetic equations.We follow a typical approach in the context of DUGKS [60] and realize the quasi-equilibrium relaxation via the Strangsplitting method [70]: 1. Quasi-equilibrium relaxation of the populations in the cell centers (half-time step), 2. Update step without quasi-equilibrium relaxation, Eqs. ( 88), (89).
3. Quasi-equilibrium relaxation of the populations in the cell centers (half-time step), as in step 1.
We note that half-time relaxations steps occur in each cell center and are local operations.By construction, the quasi-equilibrium relaxation conserves the flow velocity and temperature, and the populations remain in their comoving reference frame.

C. Multiscale implementation
The presented framework can be implemented in a multiscale setting with minimal changes in the algorithm.The different lattices are deployed adaptively in the simulation domain following a switching criterion.According to Sec.II D, the switching criterion is a function of the hydrodynamic gradients, with the high-order lattice being activated in the regions of steep gradients.In this work, the switching function consists of threshold criteria on the numerically computed flow velocity and temperature gradients.The different lattices are updated normally as presented in the previous section, with the difference being that the reference frame transformations in the vicinity of the interface regions are replaced by the multiscale frame transformations (presented in Sec.II E).

D. Boundary conditions
The boundary conditions (BCs) are enforced in the current work via the ghost node approach [71].First, the density, flow velocity and temperature are determined at the ghost cell C G (see Fig. 1).For fixed values at the wall, e.g.no slip velocity u w , the value u G at the ghost cell is, where u B is the corresponding value at the boundary cell.To impose zero normal gradient condition, e.g. for density computation, we enforce With the macroscopic values (ρ G , u G , T G ) defined, the reference frame of the ghost cell λ G is set to the comoving reference frame, λ G = {u G , T G }.The equilibrium populations are then, The approximation of non-equilibrium contributions follows the implementation of [42].In particular, the firstorder non-equilibrium moments are estimated from the Chapman-Enskog solution, and they depend on the local hydrodynamic gradients.The pertinent non-equilibrium moments of the f − populations are [42], Q (1) = −τ 2 ρ G T G (∇T I) + uP (1) . (112) The zeroth up to second order non-equilibrium moments of the g− populations, M g,0 , M g,1 , M g,2 are estimated as [42], M (1) The hydrodynamic gradients are evaluated with a second-order centered scheme, based on previous time step quantities.The non-equilibrium populations are computed from their non-equilibrium moments, according to the Grad's projection procedure, g,0 , M g,1 , M g,2 ; λ G )H (n) (c i ). (117)

V. RESULTS AND DISCUSSION
In this section, we validate the model with 1D/2D Euler gas dynamics benchmarks, and viscous flows to assess the Prandtl number as well as the accuracy of the wall BCs.Subsequently, we focus on the shock structure problem and demonstrate numerically the implications of the moment analysis of Sec.II D. The framework is then implemented with the multiscale setting (Sec.II E), via the deployment of different lattices across the simulation domain.We conclude this section with a summary of our observations and discussion of the model capabilities.We remind that the g− populations evolve with the D2Q9 lattice.Unless stated otherwise, the numerical parameters of the simulations are the following.The time step δt is such that the Courant-Friedrichs-Lewy (CFL) number is CFL = max |v iα |(δt/δx) = 0.2, where δx is the grid resolution.The adiabatic exponent is γ = 1.4.Additionally, the viscosity for the Euler flows is low enough such that the results remain invariant (typically µ ∼ O(10 −3 − 10 −2 )).Finally, we note that the formulation of the initial and boundary conditions are based on non-dimensional variables, scaled with appropriate reference density, velocity and pressure.

A. Euler gas dynamics
We validate the model using the D2Q16 lattice and a third-order Grad's projection for the moment transformation.According to the moment analysis in Sec.II D, the hydrodynamics at the Euler level should be captured accurately.Indeed, the model performs very well against a series of 1D Riemann problems, involving low density-near vacuum regions and very strong discontinuities.While all benchmarks of the previous work [57] were tested, we present here two representative 1D examples.The 2D cases include a high Mach Riemann problem, a Mach 3 flow over a step obstacle and a shock diffraction over a corner.

Strong shock tube
We consider the case of a strong shock tube [72], where the ratio between the temperature of the left and right side is 10 5 .The initial conditions for this problem are, (ρ, u x , p) = (1, 0, 1000), 0 ≤ x < 0.5, (1, 0, 0.01), 0.5 This problem, characterized by the strong temperature discontinuity, probes the robustness and accuracy of the numerical methods.The results of the simulation, at t = 0.012 and L = 800, are shown in Fig. 2. Overall, a very good agreement with the exact solution is noted.
The results of the density field, as well as density contours near the center of the domain, are depicted in Fig. 4. The initial conditions of the Riemann problem lead to shock wave interaction and the formation of complex patterns.The results show a very good agreement with the reference solutions in [74,75].

Mach 3 flow over step
In this problem, a uniform Mach 3 flow is imposed on a wind tunnel containing a step [76].A transient shock wave develops from the step, reflects at the walls and forms a complicated flow pattern.The computational domain is bounded by a [0, 3] × [0, 1] rectangle, while the step is located at (0.6, 0) and has a height of ∆y = 0.2.Initially a gas with γ = 1.4 is spatially uniform, with the following hydrodynamic conditions, (ρ, u x , u y , p) = (1.4,3, 0, 1). (121) The same conditions are imposed as inflow BCs at the left boundary x = 0 and outflow BCs at the right bound-

Shock diffraction over corner
Here we investigate the shock diffraction problem, in which a shock wave flows over a backward facing corner [15].The hydrodynamic patterns of this problem have been studied theoretically, experimentally and via simulations.From the numerical standpoint however, this problem has been challenging due to the development of negative pressure and/or density around the corner.We follow the conventional setup of the problem: the computational domain consists of the union of [0, 1] × [6,11] and [1,13] × [0, 11] rectangles.Initially, a Ma = 5.09 right-moving shock wave is located at x = 0.5 and 6 ≤ y ≤ 11 and propagates into undisturbed air, with density 1.4 and pressure 1.For the BCs, we use inflow with the initial conditions at x = 0, 0 ≤ y ≤ 11, outflow at x = 13, 0 ≤ y ≤ 11, 1 ≤ x ≤ 13, y = 0 and 0 ≤ x ≤ 13, y = 11.Reflective BCs are applied at the walls of the domain 0 ≤ x ≤ 1, y = 6 and x = 1, 0 ≤ y ≤ 6.The results, for resolution [390,330] and t = 2.3, are shown in Fig. 6 and compare very well with the reference results from [15].

B. Viscous flows
In this section, we focus on hydrodynamic flows with important viscous effects.The discussion pivots around the accuracy and the limitations of the third-order moment invariant system, sustained by the D2Q16 lattice.

Channel flow
We begin with an isothermal channel flow with Reynolds number of Re = 100, to assess the wall BCs.The results for a simulation with 80 grid points, shown in Fig. 7, demonstrate an excellent agreement with the analytical solution.Additionally, a convergence order study with respect to the L 2 error, verifies a second order spatial convergence of the scheme.

Thermal Couette flow
The thermal Couette flow is a benchmark test case to probe the viscous heat dissipation and the Prandtl number.The upper wall with the higher temperature T H is in motion with a constant speed u 0 , while the lower wall is at rest and at a temperature T C .The analytical solution for the temperature is, where Ec = u 2 0 /(C p ∆T ) is the Eckert number and ∆T = T H − T C .No slip and constant temperature BCs are applied at the top and bottom walls, while periodic BCs are enforced in the horizontal direction.The parameters for the simulations are Ma = u 0 / √ γT C = 0.5, L = 150, Re = ρu 0 L/µ = 100, T C = 1.Fig. 8 shows the temperature profiles for three different Prandtl numbers (Pr = 0.5, 0.7, 1.0) and different Eckert numbers (Ec = 4, 20, 40), which are in very good agreement with the analytical solution.We note that the simulations have been performed with the D2Q16 lattice and thus a third-order Grad's projection frame transformation.The accuracy of the results suggest that the error term in the energy equation ( 60) is negligibly small.

Viscous shock tube
In this problem we probe the performance of our model with the viscous shock tube test, proposed by Daru and 0.0 0.2 0.4 0.6 0.8 1.0 x/L Tenaud [77].A 2D shock tube [0, 1] × [0, 1] is initialized with the following conditions, (ρ, u x , u y , p) = (120, 0, 0, 120/γ), 0 ≤ x < 0.5, (1.2, 0, 0, 1.2/γ), 0.5 where γ = 1.4, the Prandtl number is set to Pr = 0.73 and the viscosity is set such that Reynolds number is Re = 200.No slip and adiabatic BCs are applied at  [78].x-TP (y-TP) correspond to the x-coordinate (ycoordinate) of the triple point.x-PV corresponds to the horizontal axis intersection of the line passing through the primary vortex.y-PV corresponds to the height of the primary vortex.
x-TP y-TP x-PV y-PV Reference 0.58 0.137 0.78 0.166 D2Q16 0.58 0.133 0.774 0.168 D2Q25 0.58 0.134 0.775 0.168 the walls of the shock tube.Due to the symmetric configuration of the problem, the actual simulated geometry consists of the [0, 1] × [0, 0.5] domain, with symmetric conditions applied on the top boundary.The initial flow conditions create a right propagating shock wave of Ma = 2.37, a contact discontinuity and an expanding rarefaction wave towards both directions.It is noted that the motion of the shock wave induces a non-negligible boundary layer along the horizontal wall of the tube.The boundary layer interacts with the incident and reflected shock, forming a complicated flow pattern.
The results of the density contours, the pressure and temperature fields for resolution of [500, 250] at t = 1, are shown in Fig. 9. Additionally, we repeat the simulation with the D2Q25 and a fourth-order frame invariant moment system and compare the results.For the comparison, we report metrics suggested from [78].In particular, Table II summarizes the coordinates associated with the triple point and the primary vortex.The comparison with the reference data show a very good match of both D2Q16 and D2Q25 simulations.The same conclusion is drawn from Fig. 10, which plots the density distribution along the solid wall.

Shock structure problem
The problems so far have demonstrated very good accuracy of the third-order frame invariant moment system and the associated D2Q16 lattice.The following benchmark involves steep hydrodynamic gradients and clearly demonstrates the limitations of the D2Q16.At the same time, the expansion of the frame invariant moment system from third-order to fourth-order (D2Q25) restores the accuracy.
The shock structure problem is a classical problem in kinetic theory of gases, in which non-equilibrium effects dominate the flow [79].We consider a quasi onedimensional plane shock wave, with an initial step of density, velocity and temperature at the center of the computational domain.The upstream and downstream flow values are connected through the Rankine-Hugoniot conditions [80].The upstream mean free path for hard sphere molecules is defined as, where p 1 , α 1 , µ 1 are the pressure, the speed of sound and the viscosity of the gas upstream of the shock, respectively.The viscosity varies with the temperature as, where for the case of hard spheres s = 0.5.The steady-state non-dimensional density, temperature, normal stress and heat flux are defined as follows, where the subscripts 1 and 2 indicate the upstream and downstream values, respectively.The Prandtl number is set to Pr = 2/3 and the adiabatic exponent of monoatomic ideal gas to γ = 5/3.The results reported for this case are the steady-state solutions and compared with the results of Ohwada [81].The origin of the coordinate system is the point with ρ n = 0.5 and x n = x/0.5 √ πλ 1 is used as the reduced coordinate.We consider first the shock structure profiles for a Mach number Ma = 1.2.Two simulations are performed with different reference frame transformation orders.In particular, we compare the performance of the D2Q16 and D2Q25 lattices, using third-and fourth-order Grad's projection respectively.The results for the density, temperature, normal stress and heat flux profiles are shown in Fig. 11.It is evident that both models perform very accurately, compared with the reference data.
We continue with the shock structure at a higher Mach of Ma = 1.6 and repeat the numerical experiments with the different frame transformation orders.The results are summarized in Fig. 12.Here, the third-order model clearly shows deviations in all the profiles, with the errors being prominent in the temperature and heat flux profiles, at the upstream part of the shock.In this case, the deviation terms due to the frame variant R eq moment are sustained, due to the steep gradients of velocity and temperature within the shock profile.Including the R eq moment list into the frame invariant list, i.e. the fourth-order model, recovers the accuracy of the model and achieves very good agreement with the reference results.

C. Multiscale framework
The final topic of interest is the multiscale extension of the scheme, with the deployment of different lattices across the domain.The switching criterion is a threshold on the local flow velocity and temperature gradients.The high-order lattice is activated at the portion of the domain with high gradients, while the low-order lattice everywhere else.

Lax tube
We demonstrate a D2Q9/D2Q16 model, with the simulation of the Lax problem [82].The initial conditions are the following, (ρ, u x , p) = (0.445, 0.698, 3.528), 0 ≤ x < 0.5, (0.5, 0, 0.571), 0.5 The simulation is performed with L = 600, until t = 0.14.Fig. 13 shows the solution obtained by the D2Q9 and D2Q16 lattices independently.While the D2Q16 model is in excellent agreement with the analytical solution, the D2Q9 model develops deviations, which manifest as overestimated density between the shock wave and the contact discontinuity.The discrepancy in the Euler level is expected for the case of D2Q9 and therefore a secondorder moment invariant system.D2Q25.Symbols: reference data from [78] and the contact discontinuity.The results of the multiscale model match again very well with the analytical solution.
(128) The results for the density profile are presented at t = 1.8 and L = 800.Fig. 15 shows the solutions of the D2Q9 and D2Q16 models and the comparison with a reference solution, obtained with characteristic-based 5th order WENO, RK4 temporal integration and resolution of 5000 points [84].Apart from a small underestimation of the post-shock waves amplitudes, it is evident that the D2Q16 captures very well the shock location and the high frequency waves.In contrast, the D2Q9 model clearly deviates from the reference solution.Fig. 16 captures the evolution of the D2Q9/D2Q16 solution and compares it with the pure D2Q16 solution.The multiscale model is almost indistinguishable from the D2Q16 model and thus with the reference solution also.It is also interesting to observe that D2Q16 is activated only in narrow regions of the domain, as shown by the spikes in Fig. 16.

High Mach Astrophysical jet
As a final test case, we consider an astrophysical jet of Mach 30, without radiative cooling [13].This case is an example of actual gas flows revealed from images of the Hubble Space Telescope and therefore is of high scientific interest.Following the configuration in [13], we initialize the computational domain [0, 2] × [−0.5, 0.5] with the following conditions, (ρ, u x , u y , p)
Outflow BCs are used around the domain, except the left boundary, where the prescribed fixed conditions are imposed.The simulation was performed with resolution [1200, 600].We compare the results between the D2Q9, D2Q16 and the multi-scale D2Q9/D2Q16 models.Fig. 17 shows a comparison of the pressure, density and temperature fields between the D2Q16 and D2Q9/D2Q16 solution.The propagation of the bow shock into the surrounding medium, as well as the developed Rayleigh Taylor instabilities within the jet cocoon, are captured in very good agreement between the two simulations.Fig. 18 depicts the distribution of the D2Q16 lattice in the computational domain.For a quantitative comparison, Fig. 19 plots the density field across three horizontal cuts of the domain.The multi-scale D2Q9/D2Q16 model and the pure D2Q16 are in excellent agreement.On the contrary, the pure D2Q9 model evolves with clear deviations, as shown in Fig. 19.

D. Discussion
We summarize the main strategies that we adopted to increase the efficiency of the PonD method, with minimal sacrifice of accuracy.The pivotal point is the identification of the frame invariant moment system for the f -and g-populations, according to the target hydrodynamic system.According to the analysis of Sec.II D, the frame invariant moment system for the g-populations should include up to second order moments.Hence, irrespective of the f -lattice, we used in all simulations in this work the D2Q9 lattice for the g-populations, decreasing the computational cost for both the g− populations update and the g-reference frame transformations.Numerical experiments with different f -and g-lattices did not reveal any appreciable effect on the stability and the accuracy of the scheme.
The multiscale formulation enables the deployment of a low-order lattice for the f -populations, in regions with smooth flow velocity and temperature variations.In accordance with observations in [54], the stability and accuracy of the solutions are well-maintained.The efficiency gains from this approach are naturally case dependent.We note that the different lattices communicate solely through the Grad's reference frame transformation, which renders the transition from a single lattice to a multiscale model easy to program and highly efficient.The last element which differs from the PonD formulations in [57] is the absence of iterations within the flux calculation, as discussed in IV A 3. We demonstrate a comparison between the iterative and the current formulation through the Shu-Osher problem [84].Fig. 20 shows the results from the two schemes, for different CFL numbers and resolutions.One observes that for high CFL and coarse domains, the iterative scheme is marginally more accurate than its explicit counterpart.For moderate CFL and resolved domains the two solutions are almost indistinguishable.Additional numerical experiments confirm the above observations.We can conclude that for resolved simulations (spatially and temporally), the non-iterative flux calculation can be safely employed.

VI. CONCLUSIONS
In this work, we presented the PonD formulation with an emphasis on the requirements of the reference frame transformation.According to the target hydrodynamic equations, conventional LBM models on a static reference frame require a set of equilibrium moment constraints.In constrast, PonD utilizes an adaptive comoving reference frame with the pertinent equilibrium moment constraints being automatically satisfied by exact equilibrium populations.However, the target hydrodynamic equations introduce requirements on the frame invariant moment system of the reference frame transformation.The framework presented on this work is a finite volume discretization of the governing kinetic equations in an adaptive reference frame.In comparison with conventional finite volume LBMs (such as conventional DUGKS), the cost to be paid for the adaptive formulation amounts to the reference frame transformations.The benefit of this approach is enhanced accuracy, stability and an increased operating window in terms of Mach number and temperature.Additionally, a multiscale extension can easily be incorporated and results in further efficiency gains.Further high Mach 3D simulations with the presence of curved boundaries shall be the focus of future work.Appendix A: Hermite polynomials The Hermite polynomials, up to fourth order and with discrete velocities scaled such that T L = 1, are the following, We insert Eq. (B3) into Eq.(B1) and obtain the following equation, where m λ Ω,j denotes the moments from the collision operator.Next, we invoke the reference frame invariance of the moments, where λ(x, t) denotes the local reference frame.We insert the moments evaluated from the local reference frame (B5) into Eq.(B4), , and retrieve the following equation, , is the definition of the reference frame transformation, from the local frame λ(x, t) to the monitoring frame λ.Thus, Eq. (B7) is the final PonD equation, where for convenience, the summation over repeated indices is not explicit, G λ i,λ(x,t) f λ(x,t) = [G λ λ(x,t) ] i,k f λ(x,t) k (B10)

Appendix C: Conservation properties
Without loss of generality, we consider a face at x I , with a unit normal vector pointing at the x-direction.According to the presented scheme, the fluxes have been calculated with a reference frame λ and are calculated as, F λ g,i (x I ) = v λ ix g λ i (x I ).(C2) The f − populations at the left and right neighbouring cells x L , x R are updated due to the fluxes F λ f,i (x I ) as, and accordingly the g− populations, The finite volume is strictly conservative with respect to mass, momentum and total energy when, δρ R + δρ L = 0, (C13) δ(ρu) R + δ(ρu) L = 0, (C14) δ(ρE) R + δ(ρE) L = 0.
(C15) Substituting from the above expressions we arrive at the following constraints, The constraints are satisfied if the following moments of the f − populations are invariant upon reference frame transformation , and the following for the g− populations, Appendix D: Forcing scheme We consider the continuous kinetic equation, with a forcing term, The body force can be expressed as The force can be incorporated by the Strang-Splitting approach, ∂ t f i = F i , (for δt/2) (D3) The intermediate step is the kinetic update without body force.In the two half-time forcing steps, the distribution function and the macroscopic velocity are updated as,

FIG. 1 .
FIG. 1. Schematic for the implementation of boundary conditions.

FIG. 5 .
FIG. 5. Density profiles for the Mach 3 flow over forward step.Six snapshots are shown at equal time intervals from t = 0 to t = 3.

FIG. 6 .
FIG. 6. Pressure (top) and density (bottom) profiles for the Mach 5.09 shock diffraction problem over a corner, at t = 2.3.30 equidistant contours are superimposed on the fields.

2 FIG. 7 .
FIG. 7. Top: Force driven flow over channel at Re = 100.Symbols correspond to simulation results and solid line to analytical solution.Bottom: convergence order analysis.Solid line corresponds to L2 error norm.Dashed line indicates a second order L2 error norm slope.

FIG. 11 .
FIG.11.The shock structure problem with Ma 1.2.Density, temperature profiles (top), normal stress and heat flux profiles (bottom).Third-order transformation at the left column, fourth-order at the right column.

FIG. 12 .FIG. 13 .
FIG.12.The shock structure problem with Ma 1.6.Density, temperature profiles (top), normal stress and heat flux profiles (bottom).3rd order transformation at the left column, 4th order at the right column.

FIG. 14 .
FIG. 14. Density profile for the Lax tube problem, at t = 0.14.The red line corresponds to the D2Q9/D2Q16 model and the black dashed line to the analytical solution.The grey dashed line indicates the occupancy regions of the D2Q16 lattice (φ = 1) and of the D2Q9 lattice (φ = 0).

FIG. 16 .
FIG. 16.Density profile for the Shu-Osher problem, at different times.The red line corresponds to the D2Q9/D2Q16 model and the black dashed line to the D2Q16 solution.The grey dashed line indicates the occupancy regions of the D2Q16 lattice (φ = 1) and of the D2Q9 lattice (φ = 0).

FIG. 20 .
FIG.20.Density profile for the Shu-Osher problem, comparing the scheme with (red dashed line) and without (blue solid line) iterations.Black solid line corresponds to the reference results[84].Results are shown for different CFL numbers and resolution.

TABLE I .
Lattice temperature TL, roots of Hermite polynomials ciα and weights wiα of the D = 1 Gauss-Hermite quadrature, and nomenclature.

TABLE II .
Accuracy criteria for viscous shock tube, according to