Dear all,
I am pretty new to Lattice Boltzmann and recently I upgraded my 2D code to 3D. However, it is not clear to me how to handle corners and edges, preferably in a generic way.
Currently, I am using the Finite-difference velocity gradient method (BC4) from Latt & Chopard’s Straight velocity boundaries in the lattice Boltzmann method. More precisely, I impose my desired velocity using the following code:
#define DELTA(a,b) (a == b)
{
// Mass fluxes with centered FD
T danub[dim][dim];
danub[0][0] = 0.5*(lab(ix+1,iy,iz).rho*lab(ix+1,iy,iz).u[0] - lab(ix-1,iy,iz).rho*lab(ix-1,iy,iz).u[0]);
danub[0][1] = 0.5*(lab(ix+1,iy,iz).rho*lab(ix+1,iy,iz).u[1] - lab(ix-1,iy,iz).rho*lab(ix-1,iy,iz).u[1]);
danub[1][0] = 0.5*(lab(ix,iy+1,iz).rho*lab(ix,iy+1,iz).u[0] - lab(ix,iy-1,iz).rho*lab(ix,iy-1,iz).u[0]);
danub[1][1] = 0.5*(lab(ix,iy+1,iz).rho*lab(ix,iy+1,iz).u[1] - lab(ix,iy-1,iz).rho*lab(ix,iy-1,iz).u[1]);
if (dim == 3) {
danub[0][2] = 0.5*(lab(ix+1,iy,iz).rho*lab(ix+1,iy,iz).u[2] - lab(ix-1,iy,iz).rho*lab(ix-1,iy,iz).u[2]);
danub[1][2] = 0.5*(lab(ix,iy+1,iz).rho*lab(ix,iy+1,iz).u[2] - lab(ix,iy-1,iz).rho*lab(ix,iy-1,iz).u[2]);
danub[2][0] = 0.5*(lab(ix,iy,iz+1).rho*lab(ix,iy,iz+1).u[0] - lab(ix,iy,iz-1).rho*lab(ix,iy,iz-1).u[0]);
danub[2][1] = 0.5*(lab(ix,iy,iz+1).rho*lab(ix,iy,iz+1).u[1] - lab(ix,iy,iz-1).rho*lab(ix,iy,iz-1).u[1]);
danub[2][2] = 0.5*(lab(ix,iy,iz+1).rho*lab(ix,iy,iz+1).u[2] - lab(ix,iy,iz-1).rho*lab(ix,iy,iz-1).u[2]);
}
// Special treatment for boundaries, where we must do one-sided FD
if (lab(ix,iy,iz).boundary & X0) {
danub[0][0] = 0.5*(4*lab(ix+1,iy,iz).rho*lab(ix+1,iy,iz).u[0]-lab(ix+2,iy,iz).rho*lab(ix+2,iy,iz).u[0]-3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[0]);
danub[0][1] = 0.5*(4*lab(ix+1,iy,iz).rho*lab(ix+1,iy,iz).u[1]-lab(ix+2,iy,iz).rho*lab(ix+2,iy,iz).u[1]-3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[1]);
if (dim == 3)
danub[0][2] = 0.5*(4*lab(ix+1,iy,iz).rho*lab(ix+1,iy,iz).u[2]-lab(ix+2,iy,iz).rho*lab(ix+2,iy,iz).u[2]-3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[2]);
}
if (lab(ix,iy,iz).boundary & XN) {
danub[0][0] = 0.5*(-4*lab(ix-1,iy,iz).rho*lab(ix-1,iy,iz).u[0]+lab(ix-2,iy,iz).rho*lab(ix-2,iy,iz).u[0]+3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[0]);
danub[0][1] = 0.5*(-4*lab(ix-1,iy,iz).rho*lab(ix-1,iy,iz).u[1]+lab(ix-2,iy,iz).rho*lab(ix-2,iy,iz).u[1]+3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[1]);
if (dim == 3)
danub[0][2] = 0.5*(-4*lab(ix-1,iy,iz).rho*lab(ix-1,iy,iz).u[2]+lab(ix-2,iy,iz).rho*lab(ix-2,iy,iz).u[2]+3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[2]);
}
if (lab(ix,iy,iz).boundary & Y0) {
danub[1][0] = 0.5*(4*lab(ix,iy+1,iz).rho*lab(ix,iy+1,iz).u[0]-lab(ix,iy+2,iz).rho*lab(ix,iy+2,iz).u[0]-3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[0]);
danub[1][1] = 0.5*(4*lab(ix,iy+1,iz).rho*lab(ix,iy+1,iz).u[1]-lab(ix,iy+2,iz).rho*lab(ix,iy+2,iz).u[1]-3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[1]);
if (dim == 3) {
danub[1][2] = 0.5*(4*lab(ix,iy+1,iz).rho*lab(ix,iy+1,iz).u[2]-lab(ix,iy+2,iz).rho*lab(ix,iy+2,iz).u[2]-3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[2]);
}
}
if (lab(ix,iy,iz).boundary & YN) {
danub[1][0]=0.5*(-4*lab(ix,iy-1,iz).rho*lab(ix,iy-1,iz).u[0]+lab(ix,iy-2,iz).rho*lab(ix,iy-2,iz).u[0]+3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[0]);
danub[1][1]=0.5*(-4*lab(ix,iy-1,iz).rho*lab(ix,iy-1,iz).u[1]+lab(ix,iy-2,iz).rho*lab(ix,iy-2,iz).u[1]+3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[1]);
if (dim == 3)
danub[1][2] = 0.5*(-4*lab(ix,iy-1,iz).rho*lab(ix,iy-1,iz).u[2]+lab(ix,iy-2,iz).rho*lab(ix,iy-2,iz).u[2]+3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[2]);
}
if (lab(ix,iy,iz).boundary & Z0) {
danub[2][0] = 0.5*(4*lab(ix,iy,iz+1).rho*lab(ix,iy,iz+1).u[0]-lab(ix,iy,iz+2).rho*lab(ix,iy,iz+2).u[0]-3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[0]);
danub[2][1] = 0.5*(4*lab(ix,iy,iz+1).rho*lab(ix,iy,iz+1).u[1]-lab(ix,iy,iz+2).rho*lab(ix,iy,iz+2).u[1]-3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[1]);
danub[2][2] = 0.5*(4*lab(ix,iy,iz+1).rho*lab(ix,iy,iz+1).u[2]-lab(ix,iy,iz+2).rho*lab(ix,iy,iz+2).u[2]-3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[2]);
}
if (lab(ix,iy,iz).boundary & ZN) {
danub[2][0] = 0.5*(-4*lab(ix,iy,iz-1).rho*lab(ix,iy,iz-1).u[0]+lab(ix,iy,iz-2).rho*lab(ix,iy,iz-2).u[0]+3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[0]);
danub[2][1] = 0.5*(-4*lab(ix,iy,iz-1).rho*lab(ix,iy,iz-1).u[1]+lab(ix,iy,iz-2).rho*lab(ix,iy,iz-2).u[1]+3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[1]);
danub[2][2] = 0.5*(-4*lab(ix,iy,iz-1).rho*lab(ix,iy,iz-1).u[2]+lab(ix,iy,iz-2).rho*lab(ix,iy,iz-2).u[2]+3*lab(ix,iy,iz).rho*lab(ix,iy,iz).u[2]);
}
// Compute equilibria
T feq[pop];
for (int k = 0; k < pop; k++)
feq[k] = o(ix,iy,iz).equilibrium(k);
// Compute the population
for (int k = 0; k < pop; k++)
{
T sum = 0;
for (int a = 0; a < dim; a++)
for (int b = 0; b < dim; b++)
sum += (c[a][k]*c[b][k] - cs2*DELTA(a,b)) * danub[a][b];
o(ix,iy,iz).f[k] = feq[k] - (t[k]*tau / cs2) * sum;
}
}
where ix,iy,iz are the coordinates of the grid cell, lab(ix,iy,iz) is a object that gives me access to the data (including ghosts for coordinates < 0 or >= n). The code works for both 2D and 3D and, of course, if can be made much more efficient by expanding the last for loop analytically.
Before calling the function above to set the populations, I use the following code to find my density (this is for D3Q19):
static void setBoundaryCellPressure(int wallPosition, Real rho, Real u[], Real f[])
{
if (wallPosition == X0) {
u[0] = 1 - 1.0 / rho * (f[3] + f[4] + f[5] + f[6] + f[15] + f[16] + f[17] + f[18] + f[0] + 2*(f[2] + f[11] + f[12] + f[13] + f[14]));
u[1] = u[2] = 0;
}
if (wallPosition == XN) {
u[0] = -1 + 1.0 / rho * (f[3] + f[4] + f[5] + f[6] + f[15] + f[16] + f[17] + f[18] + f[0] + 2*(f[1] + f[7] + f[8] + f[9] + f[10]));
u[1] = u[2] = 0;
}
if (wallPosition == Y0){
u[0] = u[2] = 0;
u[1] = 1 - 1.0 / rho * (f[1] + f[2] + f[5] + f[6] + f[9] + f[10] + f[13] + f[14] + f[0] + 2*(f[4] + f[8] + f[12] + f[17] + f[18]));
}
if (wallPosition == YN){
u[0] = u[2] = 0;
u[1] = -1 + 1.0 / rho * (f[1] + f[2] + f[5] + f[6] + f[9] + f[10] + f[13] + f[14] + f[0] + 2*(f[3] + f[7] + f[11] + f[15] + f[16]));
}
if (wallPosition == Z0){
u[0] = u[1] = 0;
u[2] = 1 - 1.0 / rho * (f[1] + f[2] + f[3] + f[4] + f[7] + f[11] + f[12] + f[8] + f[0] + 2*(f[6] + f[10] + f[14] + f[16] + f[18]));
}
if (wallPosition == ZN){
u[0] = u[1] = 0;
u[2] = -1 + 1.0 / rho * (f[1] + f[2] + f[3] + f[4] + f[7] + f[11] + f[12] + f[8] + f[0] + 2*(f[5] + f[9] + f[13] + f[15] + f[17]));
}
}
My question is how to extrapolate the values for the densities for corners and edges. I have read in the “Straight velocity boundaries in the lattice Boltzmann method” that “a second-order accurate extrapolation is used to evaluate the pressure in the 4 corners in 2D problems, and in the 8 corners and the 12 edges in 3D problems”, but I don’t know how to do that … I would be very grateful if anybody could point me some paper, formula or piece of code which explains how this can be done.
Thanks in advance,
Serban