# LBM Examples

## 2D flow around a cylinder

This example introduces the LBM module in 2D for a flow around a cylinder in 2d. The example code is as follows:
```                                    struct UserData
{
Array<Cell *> Left;
Array<Cell *> Right;
Array<double> Vel;
double        vmax;
double        rho;
};

void Setup (LBM::Domain & dom, void * UD)
{
UserData & dat = (*static_cast<UserData *>(UD));
// Cells with prescribed velocity
#ifdef USE_OMP
#pragma omp parallel for schedule(static) num_threads(dom.Nproc)
#endif
for (size_t i=0; i<dat.Left.Size(); ++i)
{
Cell * c = dat.Left[i];
if (c->IsSolid) continue;
double rho = (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[3]+c->F[6]+c->F[7]))/(1.0-dat.Vel[i]);
c->F[1] = c->F[3] + (2.0/3.0)*rho*dat.Vel[i];
c->F[5] = c->F[7] + (1.0/6.0)*rho*dat.Vel[i] - 0.5*(c->F[2]-c->F[4]);
c->F[8] = c->F[6] + (1.0/6.0)*rho*dat.Vel[i] + 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
}

// Cells with prescribed density
#ifdef USE_OMP
#pragma omp parallel for schedule(static) num_threads(dom.Nproc)
#endif
for (size_t i=0; i<dat.Right.Size(); ++i)
{
Cell * c = dat.Right[i];
if (c->IsSolid) continue;
double vx = -1.0 + (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[1]+c->F[5]+c->F[8]))/c->RhoBC;
c->F[3] = c->F[1] - (2.0/3.0)*c->RhoBC*vx;
c->F[7] = c->F[5] - (1.0/6.0)*c->RhoBC*vx + 0.5*(c->F[2]-c->F[4]);
c->F[6] = c->F[8] - (1.0/6.0)*c->RhoBC*vx - 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
}
}

int main(int argc, char **argv) try
{
size_t nproc = 1;
if (argc==2) nproc=atoi(argv[1]);
double u_max  = 0.1;                // Poiseuille's maximum velocity
double Re     = 100000.0;                  // Reynold's number
size_t nx = 800;
size_t ny = 100;
int radius = ny/10 + 1;           // radius of inner circle (obstacle)
double nu     = u_max*(2*radius)/Re; // viscocity
LBM::Domain Dom(D2Q9, nu, iVec3_t(nx,ny,1), 1.0, 1.0);
UserData dat;
Dom.UserData = &dat;

dat.vmax = u_max;
//Assigning the left and right cells
for (size_t i=0;i<ny;i++)
{
dat.Left .Push(Dom.Lat[0].GetCell(iVec3_t(0   ,i,0)));
dat.Right.Push(Dom.Lat[0].GetCell(iVec3_t(nx-1,i,0)));

// set parabolic profile
double L  = ny - 2;                       // channel width in cell units
double yp = i - 1.5;                      // ordinate of cell
double vx = dat.vmax*4/(L*L)*(L*yp - yp*yp); // horizontal velocity
dat.Vel.Push(vx);
dat.Left [i]->RhoBC = 1.0;
//dat.Right[i]->VelBC = 0.08,0.0,0.0;
dat.Right[i]->RhoBC = 1.0;

}
dat.rho  = 1.0;

// set inner obstacle
int obsX   = ny/2;   // x position
int obsY   = ny/2+3; // y position
for (size_t i=0;i<nx;i++)
for (size_t j=0;j<ny;j++)
{
{
Dom.Lat[0].GetCell(iVec3_t(i,j,0))->IsSolid = true;
}
}

//Assigning solid boundaries at top and bottom
for (size_t i=0;i<nx;i++)
{
Dom.Lat[0].GetCell(iVec3_t(i,0   ,0))->IsSolid = true;
Dom.Lat[0].GetCell(iVec3_t(i,ny-1,0))->IsSolid = true;
}

double rho0 = 1.0;
Vec3_t v0(0.08,0.0,0.0);

//Initializing values
for (size_t i=0;i<Dom.Lat[0].Ncells;i++)
{
Dom.Lat[0].Cells[i]->Initialize(rho0, v0);
}

//Solving
Dom.Time = 0.0;
Dom.Solve(40000.0,80.0,Setup,NULL,"tlbm01",true,nproc);

}
MECHSYS_CATCH
```

### Boundary Conditions

It can be noticed that both an structure `UserData` and one function `Setup` are defined. These two items allow the user to interact with the internal for loop representing the time evolution inside the function `Solve`. The `UserData` contains the information the user may want to control or report, in this case it contains five elements. It is worth noting that the boundary conditions of this problem involves a velocity boundary condition for the fluid entering trough the left limit and a pressure (density) boundary condition in the right hand size.

```                                    struct UserData
{
Array Left;
Array Right;
Array Vel;
double        vmax;
double        rho;
};
```

The first array `Left` will store the pointers to the cells in the left boundary of the domain, and `Right` does the same with the right hand size. The array of doubles `Vel` will store the velocity profile which in the case of a pipe flow is the well known parabolic profile. Finally `vmax` and `rho` are both the maximum velocity of the profile and the density for the right hand size boundary.

Boundary conditions will be imposed with the function `Setup`:

```                                    void Setup (LBM::Domain & dom, void * UD)
{
UserData & dat = (*static_cast<UserData *>(UD));
// Cells with prescribed velocity
#ifdef USE_OMP
#pragma omp parallel for schedule(static) num_threads(dom.Nproc)
#endif
for (size_t i=0; i<dat.Left.Size(); ++i)
{
Cell * c = dat.Left[i];
if (c->IsSolid) continue;
double rho = (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[3]+c->F[6]+c->F[7]))/(1.0-dat.Vel[i]);
c->F[1] = c->F[3] + (2.0/3.0)*rho*dat.Vel[i];
c->F[5] = c->F[7] + (1.0/6.0)*rho*dat.Vel[i] - 0.5*(c->F[2]-c->F[4]);
c->F[8] = c->F[6] + (1.0/6.0)*rho*dat.Vel[i] + 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
}

// Cells with prescribed density
#ifdef USE_OMP
#pragma omp parallel for schedule(static) num_threads(dom.Nproc)
#endif
for (size_t i=0; i<dat.Right.Size(); ++i)
{
Cell * c = dat.Right[i];
if (c->IsSolid) continue;
double vx = -1.0 + (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[1]+c->F[5]+c->F[8]))/c->RhoBC;
c->F[3] = c->F[1] - (2.0/3.0)*c->RhoBC*vx;
c->F[7] = c->F[5] - (1.0/6.0)*c->RhoBC*vx + 0.5*(c->F[2]-c->F[4]);
c->F[6] = c->F[8] - (1.0/6.0)*c->RhoBC*vx - 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
}
}
```

The first line of the function will cast the void pointer `UD` into one of the clase `UserData`. The next stage is a for loop going trough all the cells in the left hand size which is the one with imposed velocity. That for is parallelized by the `pragma` precompiler command which takes the variable `Nproc` of the `LBM::Domain` class which will be set up later. Inside the loop some elements are important:

```                                    		if (c->IsSolid) continue;
double rho = (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[3]+c->F[6]+c->F[7]))/(1.0-dat.Vel[i]);
c->F[1] = c->F[3] + (2.0/3.0)*rho*dat.Vel[i];
c->F[5] = c->F[7] + (1.0/6.0)*rho*dat.Vel[i] - 0.5*(c->F[2]-c->F[4]);
c->F[8] = c->F[6] + (1.0/6.0)*rho*dat.Vel[i] + 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
```

The flag `IsSolid` of the `Cell` class is a boolean variable determining if the particle is a solid node or not. If it is solid the boundary condition does not apply and the bounce-back condition will be imposed. The next lines are the formulas for the distribution functions `F` going in the directions 1,5 and 8 which are the only ones with components going to the right direction. They will ensure that the velocity on those nodes is the same given by the array `dat.Vel` which contains the parabolic velocity profile as mentioned before. These formulas as well as the bounce back condition can be found in the LBM literature with ease.

Similarly the next for loop goes through the right boundary. The formulas are similar but now they ensure that the density on those nodes is the ones given by the variable `RhoBC`.

```                                    		Cell * c = dat.Right[i];
if (c->IsSolid) continue;
double vx = -1.0 + (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[1]+c->F[5]+c->F[8]))/c->RhoBC;
c->F[3] = c->F[1] - (2.0/3.0)*c->RhoBC*vx;
c->F[7] = c->F[5] - (1.0/6.0)*c->RhoBC*vx + 0.5*(c->F[2]-c->F[4]);
c->F[6] = c->F[8] - (1.0/6.0)*c->RhoBC*vx - 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
```

### The main scope

The first line of the main scope declares an integer variable `nproc` which will be imported from the command line when the program is called.

```                                        size_t nproc = 1;
if (argc==2) nproc=atoi(argv[1]);
```

The next steps assign the maximun velocity of the parabolic velocity profile in the left hand size, the Reynolds number, the dimensions of the array `nx,ny` the radius of the obstacle and the viscosity `nu` obtained from the Reynolds number and the velocity. Then the LBM domain is created. The first argument from the constructor is the type of LBM method which in this case is the `D2Q9` scheme whose definition can be easily found in the LBM literature. The second argument is the viscosity of the fluid. The third is a vector of integers containing the lattice dimensions. The last two arguments are the grid size and time step which are taken as one. Finally an object of the `UserData` class is declared and assigned to the variable `Dom.UserData`. This step is necessary for the application of the boundary conditions explained above.

```                                        double u_max  = 0.1;                // Poiseuille's maximum velocity
double Re     = 100000.0;                  // Reynold's number
size_t nx = 800;
size_t ny = 100;
int radius = ny/10 + 1;           // radius of inner circle (obstacle)
double nu     = u_max*(2*radius)/Re; // viscocity
LBM::Domain Dom(D2Q9, nu, iVec3_t(nx,ny,1), 1.0, 1.0);
UserData dat;
Dom.UserData = &dat;
```

Next, the arrays `Left, Right` and `Vel` from the `UserData`. For it the command `GetCell` is useful to get a pointer to the cell in the coordinates given by the integer vector serving as argument. Also the density for the cells at the right hand size are imposed.

```                                        dat.vmax = u_max;
//Assigning the left and right cells
for (size_t i=0;i<ny;i++)
{
dat.Left .Push(Dom.Lat[0].GetCell(iVec3_t(0   ,i,0)));
dat.Right.Push(Dom.Lat[0].GetCell(iVec3_t(nx-1,i,0)));

// set parabolic profile
double L  = ny - 2;                       // channel width in cell units
double yp = i - 1.5;                      // ordinate of cell
double vx = dat.vmax*4/(L*L)*(L*yp - yp*yp); // horizontal velocity
dat.Vel.Push(vx);
dat.Left [i]->RhoBC = 1.0;
//dat.Right[i]->VelBC = 0.08,0.0,0.0;
dat.Right[i]->RhoBC = 1.0;

}
dat.rho  = 1.0;
```

The obstacle is obtained by going through the whole domain and setting some cells as solid if they fulfill the circle equation. The obstacle is positioned above the middle line to ensure the appearance of turbulent eddies. Also the upper and bottom boundaries are set as solid as well to model pipe flow.

```                                    	int obsX   = ny/2;   // x position
int obsY   = ny/2+3; // y position
for (size_t i=0;i<nx;i++)
for (size_t j=0;j<ny;j++)
{
{
Dom.Lat[0].GetCell(iVec3_t(i,j,0))->IsSolid = true;
}
}

//Assigning solid boundaries at top and bottom
for (size_t i=0;i<nx;i++)
{
Dom.Lat[0].GetCell(iVec3_t(i,0   ,0))->IsSolid = true;
Dom.Lat[0].GetCell(iVec3_t(i,ny-1,0))->IsSolid = true;
}

```

Finally an initial velocity and density are assigned for the fluid for all the cells. Then the function `Solve` is called. This function is very similar as the one called in the example of two polyhedra colliding in the DEM domain. the one difference is that is does not require a time step since it was declared in the constructor. Again the first argument is the final time and the second is the timespan for outputs. The third is the name of the function `Setup` containing the boundary condition formulas. This function is going to be called at each time step. The fourth argument is a `NULL` pointer which can be filled with another function similar to `Setup` with the only difference being that this one is going to be called only during each output step. The fifth argument is the file key. The fourth is true if you want to visualize the results and the last one is the number of processors for parallelazation.

```                                        double rho0 = 1.0;
Vec3_t v0(0.08,0.0,0.0);

//Initializing values
for (size_t i=0;i<Dom.Lat[0].Ncells;i++)
{
Dom.Lat[0].Cells[i]->Initialize(rho0, v0);
}

//Solving
Dom.Time = 0.0;
Dom.Solve(40000.0,80.0,Setup,NULL,"tlbm01",true,nproc);
```

### Running the simulation and visualizing the results

After following the compiling instructions the example should be run by typing:

```                                user@mach: ~/msys_tutorials \$ ./test_cylinder 8
```

Where the number 8 is the number of processors. Please check your hardware specifications before running.

Visualizing requires again the Visit software as in the colliding example. The only difference is that now the pseudocolor layer should be chosen for the velocity field. Likewise the velocity as vectors can also be plotted. The velocity field should look like this video: