# LBM BGK

Hi all!

I just started with the LBM. I correctly simulated 2D Poiseuille flow with a pressure gradient between the inlet and outlet (rho_in rho_out = 1.1 and = 1.0). For this, I used two distribution functions f [lx] [ly] [q] and f_old [lx] [ly] [q].

In the public void initialize (), I initialized. And the function void collision_streaming_bgk (), in the collision, I used f [x] [y] [k] éffectuer calculation, and the streaming part and the boundary conditions (pressure, bounce-back) I used f_old [x] [y] [k]. And prior to recalculate the macroscopic variables (density rho, and velocity_x velocity_y), I made a swapping between f_old [x] [y] [k] and f [x] [y] [k].

I would like to now work with only one distribution function f [x] [y] [k], without using f_old [x] [y] [k], do you have any suggestions for me as that should proceed.

Thank you for help!!!

typedef double T;

T *** f, ***f_old; // LBM populations (old and new)
T *feq; // Equilibrium pouplation
T **rho; // Fluid density
T **velocity_x; // Fluid velocity x-component
T **velocity_y; // Fluid velocity y-component

// This function computes the equilibrium function from the fluid density and velocity.

void equilibrium(T density, T vel_x, T vel_y) {

``````    T f1 = (T)3;
T f2 = (T)9/(T)2;
T f3 = (T)3/(T)2;

feq[0] = weights[0] * density * ((T)1                                                                          - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[1] = weights[1] * density * ((T)1 + f1 * vel_x                + f2 * SQ(vel_x)              - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[2] = weights[2] * density * ((T)1 + f1 * vel_y                + f2 * SQ(vel_y)              - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[3] = weights[3] * density * ((T)1 - f1 * vel_x                 + f2 * SQ(-vel_x)             - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[4] = weights[4] * density * ((T)1 - f1 * vel_y                 + f2 * SQ(-vel_y)             - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[5] = weights[5] * density * ((T)1 + f1 * ( vel_x + vel_y) + f2 * SQ( vel_x + vel_y) - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[6] = weights[6] * density * ((T)1 + f1 * (-vel_x + vel_y) + f2 * SQ(-vel_x + vel_y) - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[7] = weights[7] * density * ((T)1 - f1 * ( vel_x + vel_y) + f2 * SQ(-vel_x - vel_y)   - f3 * (SQ(vel_x) + SQ(vel_y)));
feq[8] = weights[8] * density * ((T)1 + f1 * ( vel_x - vel_y) + f2 * SQ( vel_x - vel_y)    - f3 * (SQ(vel_x) + SQ(vel_y)));
``````

}

// Allocate memory for the fluid density, velocity and Initialize the fluid density and velocity. Start with unit density and zero velocity.

void initialize() {

``````    rho          = new T*[lx];
velocity_x   = new T*[lx];
velocity_y   = new T*[lx];

for (int X = 0; X < lx; ++X) {

rho[X] = new T[ly];

velocity_x[X] = new T[ly];
velocity_y[X] = new T[ly];
}

for (int X = 0; X < lx; ++X) {

for (int Y = 0; Y < ly; ++Y) {

rho[X][Y] = 1.;

velocity_x[X][Y] = 0.;
velocity_y[X][Y] = 0.;
}
}

// Allocate memory for the populations and Initialise density distribution function with equilibrium fo zero density.

feq   = new T[Q];
f       = new T**[lx];
f_old = new T**[lx];

for (int X = 0; X < lx; ++X) {

f[X]        = new T*[ly];
f_old[X] = new T*[ly];

for (int Y = 0; Y < ly; ++Y) {

f[X][Y]     = new T[Q];
f_old[X][Y] = new T[Q];

for  (int iD = 0; iD < Q; ++iD) {

feq[iD] = 0.;

f[X][Y][iD]        = 0.;
f_old[X][Y][iD] = 0.;
}
}
}

for (int X = 0; X < lx; ++X) {

for (int Y = 0; Y < ly; ++Y) {

equilibrium(rho[X][Y], velocity_x[X][Y], velocity_y[X][Y]);

for  (int iD = 0; iD < Q; ++iD) {

f[X][Y][iD]        =  feq[iD];
f_old[X][Y][iD] =  feq[iD];
}
}
}
``````

}

/// Princal Lattice Boltzmann Method LBM BGK: Computed Collision, streaming and boundaries conditions steps.

void collision_streaming_bgk() {

``````    // Computation Collision, streaming and Boundaries conditions steps.

for (int X = 0; X < lx; ++X) {

int Xn = (X > 0)?(X-1):(lx-1);
int Xp = (X < lx-1)?(X+1):(0);

for (int Y = 0; Y < ly; ++Y) {

int Yn = (Y > 0)?(Y-1):(ly-1);
int Yp = (Y < ly-1)?(Y+1):(0);

// Compute equilibrium
// The pequilibrium population are computed.

equilibrium(rho[X][Y], velocity_x[X][Y], velocity_y[X][Y]);

// Compute new population
// This is the lattice boltzmann equation (combined collision and streaming) including external forcing.

f[X][Y][0] = f[X][Y][0] * ((T)1 - omega) + feq[0] * omega;
f[X][Y][1] = f[X][Y][1] * ((T)1 - omega) + feq[1] * omega;
f[X][Y][2] = f[X][Y][2] * ((T)1 - omega) + feq[2] * omega;
f[X][Y][3] = f[X][Y][3] * ((T)1 - omega) + feq[3] * omega;
f[X][Y][4] = f[X][Y][4] * ((T)1 - omega) + feq[4] * omega;
f[X][Y][5] = f[X][Y][5] * ((T)1 - omega) + feq[5] * omega;
f[X][Y][6] = f[X][Y][6] * ((T)1 - omega) + feq[6] * omega;
f[X][Y][7] = f[X][Y][7] * ((T)1 - omega) + feq[7] * omega;
f[X][Y][8] = f[X][Y][8] * ((T)1 - omega) + feq[8] * omega;

// Streaming steps

f_old[X][Y][0]   = f[X][Y][0];
f_old[Xp][Y][1]  = f[X][Y][1];
f_old[X][Yp][2]  = f[X][Y][2];
f_old[Xn][Y][3]  = f[X][Y][3];
f_old[X][Yn][4]  = f[X][Y][4];
f_old[Xp][Yp][5] = f[X][Y][5];
f_old[Xn][Yp][6] = f[X][Y][6];
f_old[Xn][Yn][7] = f[X][Y][7];
f_old[Xp][Yn][8] = f[X][Y][8];

}
}

// Dirichlet Pressure Boundaries conditions
// Zou and He computed Boundaries conditions

// Inlet (South side)

for (int X = 0; X < lx; ++X) {

int Y = 0;

rho[X][Y] = rho_in;

velocity_x[X][Y] = (T)0;
velocity_y[X][Y] = rho_in - (f_old[X][Y][0] + f_old[X][Y][1] + f_old[X][Y][3] + (T)2 * (f_old[X][Y][4] + f_old[X][Y][7] + f_old[X][Y][8]));

f_old[X][Y][2] = f_old[X][Y][4] + (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][5] = f_old[X][Y][7] + (T)1/(T)6 * velocity_y[X][Y] - (T)1/(T)2 * (f_old[X][Y][1]-f_old[X][Y][3]);
f_old[X][Y][6] = f_old[X][Y][8] + (T)1/(T)6 * velocity_y[X][Y] - (T)1/(T)2 * (f_old[X][Y][3]-f_old[X][Y][1]);
}

// Outlet (East side)

for (int X = 0; X < lx; ++X) {

int Y = ly-1;

rho[X][Y] = rho_out;

velocity_x[X][Y] = (T)0;
velocity_y[X][Y] = -rho_out + (f_old[X][Y][0] + f_old[X][Y][1] + f_old[X][Y][3] + (T)2 * (f_old[X][Y][2] + f_old[X][Y][5] + f_old[X][Y][6]));

f_old[X][Y][4] = f_old[X][Y][2] - (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][7] = f_old[X][Y][5] - (T)1/(T)6 * velocity_y[X][Y] + (T)1/(T)2 * (f_old[X][Y][1]-f_old[X][Y][3]);
f_old[X][Y][8] = f_old[X][Y][6] - (T)1/(T)6 * velocity_y[X][Y] + (T)1/(T)2 * (f_old[X][Y][3]-f_old[X][Y][1]);
}

// Bounce-back
// Due to the presence of the rigid walls at x = 0 and x = lx-1, the population have to be bounced back.

// left wall (West side)

for (int Y = 1; Y < ly-1; ++Y) {

int X = 0;

velocity_x[X][Y]   = (T)0;
velocity_y[X][Y]   = (T)0;

rho[X][Y]  = velocity_x[X][Y]+f_old[X][Y][0]+f_old[X][Y][2]+f_old[X][Y][4]+(T)2*(f_old[X][Y][3]+f_old[X][Y][6]+f_old[X][Y][7]);
f_old[X][Y][1] = f_old[X][Y][3]+(T)2/(T)3*velocity_x[X][Y];
f_old[X][Y][5] = f_old[X][Y][7]+(T)1/(T)6*velocity_x[X][Y]+(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])+(T)1/(T)2*velocity_y[X][Y];
f_old[X][Y][8] = f_old[X][Y][6]+(T)1/(T)6*velocity_x[X][Y]-(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])-(T)1/(T)2*velocity_y[X][Y];
}

// right wall (East side)

for (int Y = 1; Y < ly-1; ++Y) {

int X = lx-1;

velocity_x[X][Y]   = (T)0;
velocity_y[X][Y]   = (T)0;

rho[X][Y]  = -velocity_x[X][Y]+f_old[X][Y][0]+f_old[X][Y][2]+f_old[X][Y][4]+(T)2*(f_old[X][Y][1]+f_old[X][Y][5]+f_old[X][Y][8]);
f_old[X][Y][3] = f_old[X][Y][1]-(T)2/(T)3*velocity_x[X][Y];
f_old[X][Y][7] = f_old[X][Y][5]-(T)1/(T)6*velocity_x[X][Y]-(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])-(T)1/(T)2*velocity_y[X][Y];
f_old[X][Y][6] = f_old[X][Y][8]-(T)1/(T)6*velocity_x[X][Y]+(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])+(T)1/(T)2*velocity_y[X][Y];
}

// Corners nodes treatment

// Bottom left corner node
{
int X = 0;
int Y = 0;

rho[X][Y] = rho[X][Y+1];   //  Extrapolation 1st order
velocity_x[X][Y]   = (T)0;
velocity_y[X][Y]   = (T)0;

f_old[X][Y][1] = f_old[X][Y][3] + (T)2/(T)3 * velocity_x[X][Y];
f_old[X][Y][2] = f_old[X][Y][4] + (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][5] = f_old[X][Y][7] + (T)1/(T)6 * (velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][6] = (T)1/(T)12 * (-velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][8] = (T)1/(T)12 * (velocity_x[X][Y]-velocity_y[X][Y]);
f_old[X][Y][0] = (T)0;

f_old[X][Y][0] = rho[X][Y] - (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);
}

// Top left corner node
{
int X = 0;
int Y = ly-1;

rho[X][Y] = rho[X][Y-1];   //  Extrapolation 1st order
velocity_x[X][Y]   = (T)0;
velocity_y[X][Y]   = (T)0;

f_old[X][Y][1] = f_old[X][Y][3] + (T)2/(T)3 * velocity_x[X][Y];
f_old[X][Y][4] = f_old[X][Y][2] - (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][8] = f_old[X][Y][6] + (T)1/(T)6 * (velocity_x[X][Y]-velocity_y[X][Y]);
f_old[X][Y][5] = (T)1/(T)12 * (velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][7] = -(T)1/(T)12 * (velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][0] = (T)0;

f_old[X][Y][0] = rho[X][Y] - (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);
}

// Bottom right corner node
{
int X = lx-1;
int Y = 0;

rho[X][Y] = rho[X][Y+1];   //  Extrapolation 1st order
velocity_x[X][Y]   = (T)0;
velocity_y[X][Y]   = (T)0;

f_old[X][Y][3] = f_old[X][Y][1] - (T)2/(T)3 * velocity_x[X][Y];
f_old[X][Y][2] = f_old[X][Y][4] + (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][6] = f_old[X][Y][8] + (T)1/(T)6 * (-velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][5] = (T)1/(T)12 * (velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][7] = -(T)1/(T)12 * (velocity_x[X][Y]-velocity_y[X][Y]);
f_old[X][Y][0] = (T)0;

f_old[X][Y][0] = rho[X][Y] - (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);

}

// Top right corner node
{
int X = lx-1;
int Y = ly-1;

rho[X][Y] = rho[X][Y-1];   //  Extrapolation 1st order
velocity_x[X][Y]   = (T)0;
velocity_y[X][Y]   = (T)0;

f_old[X][Y][3] = f_old[X][Y][1] - (T)2/(T)3 * velocity_x[X][Y];
f_old[X][Y][4] = f_old[X][Y][2] - (T)2/(T)3 * velocity_y[X][Y];
f_old[X][Y][7] = f_old[X][Y][5] - (T)1/(T)6 * (velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][6] = (T)1/(T)12 * (-velocity_x[X][Y]+velocity_y[X][Y]);
f_old[X][Y][8] = (T)1/(T)12 * (velocity_x[X][Y]-velocity_y[X][Y]);
f_old[X][Y][0] = (T)0;

f_old[X][Y][0] = rho[X][Y] - (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);

}

// Swap popualations
// The present code used old and new populations which are swapped at the beginning of each time step.
// This way, the old populations are not overwritten during propagation.
// The resulting code is easier to write and debug.

T ***swap_tmp;

swap_tmp = f_old;
f_old    = f;
f        = swap_tmp;

// Compute fluid density and velocity
// This function computes the fluid density and velocity from the populations.

momenta();
``````

}

// Compute fluid density and velocity

void momenta() {

``````    for (int X = 0; X < lx; ++X) {

for (int Y = 0; Y < ly; ++Y) {

rho[X][Y]           =  f[X][Y][0] + f[X][Y][1] + f[X][Y][2] +  f[X][Y][3] + f[X][Y][4] + f[X][Y][5]  + f[X][Y][6] + f[X][Y][7] + f[X][Y][8];
velocity_x[X][Y] = (f[X][Y][1] + f[X][Y][5] + f[X][Y][8] - (f[X][Y][3] + f[X][Y][6] + f[X][Y][7])) /rho[X][Y];
velocity_y[X][Y] = (f[X][Y][2] + f[X][Y][5] + f[X][Y][6] - (f[X][Y][4] + f[X][Y][7] + f[X][Y][8])) /rho[X][Y];
}
}
``````

}

Hi and welcome to LBM!
If you are interested in implementing a “in place propagation” requiring only one array instead of two, you should have a look at http://arxiv.org/pdf/1111.0922.pdf, the section about the AA patern.
The idea is to use two functions: one compute a collision only and the other one compute a streaming-collision-streaming. This way there is no risk of overwriting some distribution functions, and you only need one set of distribution functions. This effectively reduces by a factor of 2 the memory used and no more swapping!
However this method is more complicated to implement. Thus you should ask yourself if you really need the memory save. After all, nowadays computers are packed with lot of ram. Also, handling boundaries with the AA pattern is difficult.
Whether you prefer to keep your code simple or to save memory is up to you.
Cheers
Nick

Thanks nick !!!