I’m working on compressible LBM and as you know the standard BGK model can’t simulate flows with high Mach number.

Most of compressible LB models use discrete velocities which are different from D2Q9 (*) like D2Q7, D2Q16, D2Q12 and other types. Usually there are more than 9 directions, sometimes there are also some directions overlapping each other with different discrete velocities. I don’t know how I should apply boundary conditions (e.g. velocity inlet or pressure outlet) in these cases. I hope someone could help me with this and it would be a great help if someone would send me a sample code in compressible LBM.
Thanks

Simple lattice Boltzmann model for simulating flows with shock wave, Yan Guangwu et.al, Physical Review E, Vol. 59, No. 1, 1999

A multi-energy-level lattice Boltzmann model for the compressible Navier–Stokes equations, Int. J. Numer. Meth. Fluids 2007; 55:41–56

Two-dimensional lattice Boltzmann model for compressible flows with high Mach number, Yanbiao Gan et. al, Physica A 387 (2008) 1721–1732

u[i][j][0] is the X coponent of the velocity
u[i][j][1] is the Y coponent of the velocity
c[n] [0] is the X coponent of the lattice velocity of the nth distribution
c[n] [1] is the Y coponent of the lattice velocity of the nth distribution
D is the number of dimensions = 2 in 2D
Q is the number of distributions = 9 in 2DQ9

now we can say

//compute the density and the velocity
Rho[i][j] = 0;
u[i][j][Dimension = 0;
for (int iF = 0; iF < Q; ++iF)
{
Rho[i][j] = Rho[i][j]+f[iF];
u[i][j][Dimension] = u[i][j][Dimension]+f[iF]*c[iF][Dimension];
}
u[i][j][Dimension]=u[i][j][Dimension]/Rho[i][j];
//compute the new density.
Rho_new[i][j] = (Rho[i][j] + Direction * Rho[i][j] * u[i][j][Dimension];
//store the new density and velocity to calculate the nuknown distributions
Rho[i][j] = Rho_new[i][j];
for (int iD = 0; iD < D; ++iD)
{
u[i][j][iD]=u_boundary[i][j][iD];
}
//compute the unknown distribution functions
for (int iF = 0; iF < Q; ++iF)
{
if(c[iF][Dimension]==-Direction * abs(c[iF][Dimension]) )
f[iF][i][j] = f[opposite[iF]]
+fEq( iF , Rho[i][j] , u[i][j] )
-fEq( opposite[iF] , Rho[i][j] , u[i][j] );
}