# Boundary condition on edges

Dear all,
I’m trying to simulate a turbulent three-dimensional flux around a
cylinder. I need to impose a given velocity at the inlet while a fixed density
condition is enforced at the outlet and periodic boundary
conditions should be set on the remaining boundaries (i. e. the sides
along y and z directions).

In order to configure in this way the boundary conditions, as soon as
the lattice is instantiated by using the following instruction

MultiBlockLattice3D<T, DESCRIPTOR> lattice(90, 30, 30,new
BGKdynamics<T,DESCRIPTOR>(1.96464) );

I activate periodic boundaries along y and z as follows

lattice.periodicity().toggle(1, true);
lattice.periodicity().toggle(2, true);

and finally I set inlet and outlet conditions with the following code

Box3D right(nx-1,nx-1, 1, ny-2, 1, nz-2);
Box3D left(0,0, 1, ny-2, 1, nz-2);

boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice,left);
boundaryCondition.setPressureConditionOnBlockBoundaries(lattice,right);

setBoundaryVelocity ( lattice, left, Array<T,3>(0.001,0.0,0.0) );
setBoundaryDensity (lattice, right, 1. );

initializeAtEquilibrium (lattice, lattice.getBoundingBox(),1.,Array<T,3>(0.0,0.0,0.0));

lattice.initialize();

My problem is that in this way I can’t control the boundary condition
enforced on the edges placed on left and right side of the box and in
these edges the velocity remains constant at 0 all along the
simulation.
I tried to extend the definition of the left and side Box3D so that
the edges are included but this seems to be in contrast with the
periodicity imposed on the other sides and the program crashes.
How should I set up the boundary conditions so that I can impose a
specific condition on the edges?

Hi,

You are right that the definition of the box “right” and “left” must be extended to include the edges. You should probably use a local boundary condition in order to combine periodicity with inlet. Which boundary condition do you use? If your program crashes even with a local boundary condition, can you post your code, or at least a simplified version of your code?

Thanks!

Thank you very much for your answer, I’ve just started using Palabos, I don’t know how to impose a local boundary condition, I’ll have a look at that.
My code is very simple, I’ve only extended to 3D the example cylinder2d, anyway I’m posting it, this is the version of the code which gives the problem on the edges:

#include “palabos3D.h”
#ifndef PLB_PRECOMPILED // Unless precompiled version is used,
#include “palabos3D.hh” // include full template code
#endif
#include
#include
#include
#include
#include

using namespace plb;
using namespace plb::descriptors;
using namespace std;

typedef double T;
#define DESCRIPTOR D3Q19Descriptor

/// A functional, used to initialize a pressure boundary to constant density
template
class ConstantDensity {
public:
ConstantDensity(T density_)
: density(density_)
{ }
T operator()(plint iX, plint iY, plint iZ) const {
return density;
}
private:
T density;
};

template
class CylinderShapeDomain3D : public plb::DomainFunctional3D {
public:
CylinderShapeDomain3D(plb::plint cx_, plb::plint cy_, plb::plint cz_, plb::plint radius)
: cx(cx_),
cy(cy_),
cz(cz_),
{ }
virtual bool operator() (plb::plint iX, plb::plint iY, plb::plint iZ) const {
return plb::util::sqr(iX-cx) + plb::util::sqr(iZ-cz) <= radiusSqr;
}
virtual CylinderShapeDomain3D* clone() const {
return new CylinderShapeDomain3D(*this);
}
private:
plb::plint cx;
plb::plint cy;
plb::plint cz;
};

/// A functional, used to instantiate bounce-back nodes at the locations of the cylinder
void cylinderSetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters,
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition )
{
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();
Box3D right(nx-1,nx-1, 1, ny-2, 1, nz-2);
Box3D left(0,0, 1, ny-2, 1, nz-2);

``````boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice,left);
boundaryCondition.setPressureConditionOnBlockBoundaries(lattice,right);

setBoundaryVelocity (lattice, left, Array<T,3>(0.001,0.0,0.0) );

setBoundaryDensity (lattice, right, ConstantDensity<T>(1.) );

initializeAtEquilibrium (lattice, lattice.getBoundingBox(),1.,Array<T,3>(0.0,0.0,0.0));

plint cx     = 30.;
plint cy     = 30.;
plint cz     = 30.;

defineDynamics(lattice, lattice.getBoundingBox(),
new plb::BounceBack<T,DESCRIPTOR>);

lattice.initialize();
``````

}

void writeGifs(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters, int iter)
{
const plint imSize = 600;
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();

``````Box3D zSlice(0, nx-1, 0, ny-1, nz/2, nz/2);
Box3D ySlice(0, nx-1, ny/2, ny/2, 0, nz-1);
Box3D extendedZslice(0, nx-1, 0, ny-1, nz/2-1, nz/2+1);
ImageWriter<T> imageWriter("leeloo");

imageWriter.writeScaledGif( createFileName("zSlice_u", iter, 6),
*computeVelocityNorm (lattice, zSlice),
imSize, imSize );
imageWriter.writeScaledGif( createFileName("ySlice_u", iter, 6),
*computeVelocityNorm (lattice, ySlice),
imSize, imSize );
imageWriter.writeScaledGif( createFileName("zSlice_omega", iter, 6),
*computeNorm (
*computeVorticity(*computeVelocity(lattice, extendedZslice)),
zSlice ),
imSize, imSize );
``````

}

void writeVTK(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters, plint iter)
{
T dx = parameters.getDeltaX();
T dt = parameters.getDeltaT();
VtkImageOutput3D vtkOut(createFileName(“vtk”, iter, 6), dx);
vtkOut.writeData(*computeVelocityNorm(lattice), “velocityNorm”, 1.);
vtkOut.writeData<3,float>(*computeVelocity(lattice), “velocity”, 1.);

``````vtkOut.writeData<float>(*computeDensity(lattice), "density", 1.);

vtkOut.writeData<3,float>(*computeVorticity(*computeVelocity(lattice)), "vorticity", 1./dt);
``````

}

int main(int argc, char* argv[]) {
plbInit(&argc, &argv);

``````global::directories().setOutputDir("./tmp/");

IncomprFlowParam<T> parameters(
(T) 1e-2,  // uMax
(T) 100.,  // Re
100,       // N resolution, a lattice of length 1 has N points
1.4,        // lx
0.05,         // ly
0.6         //lz
);
const T logT     = (T)0.02;
const T imSave   = (T)0.06;
const T vtkSave  = (T)0.06;
const T maxT     = (T)2.;

writeLogFile(parameters, "Poiseuille flow");

T cSmago = 0.32;
MultiBlockLattice3D<T, DESCRIPTOR> lattice (
parameters.getNx(), parameters.getNy(), parameters.getNz(),
new SmagorinskyBGKdynamics<T,DESCRIPTOR>(parameters.getOmega(), cSmago) );

//MultiBlockLattice3D<T, DESCRIPTOR> lattice (
//	parameters.getNx(), parameters.getNy(), parameters.getNz(),
//    new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

//Periodicity in y and z
lattice.periodicity().toggle(1, true);
lattice.periodicity().toggle(2, true);

OnLatticeBoundaryCondition3D<T,DESCRIPTOR>*
boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();

cylinderSetup(lattice, parameters, *boundaryCondition);

// Main loop over time iterations.
for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
// At this point, the state of the lattice corresponds to the
//   discrete time iT. However, the stored averages (getStoredAverageEnergy
//   and getStoredAverageDensity) correspond to the previous time iT-1.

if (iT%parameters.nStep(imSave)==0) {
pcout << "Saving Gif ..." << endl;
writeGifs(lattice, parameters, iT);
}

if (iT%parameters.nStep(vtkSave)==0 && iT>0) {
pcout << "Saving VTK file ..." << endl;
writeVTK(lattice, parameters, iT);
}

if (iT%parameters.nStep(logT)==0) {
pcout << "step " << iT
<< "; t=" << iT*parameters.getDeltaT();
}

// Lattice Boltzmann iteration step.
lattice.collideAndStream();

// At this point, the state of the lattice corresponds to the
//   discrete time iT+1, and the stored averages are upgraded to time iT.
if (iT%parameters.nStep(logT)==0) {
pcout << "; av energy ="
<< setprecision(10) << getStoredAverageEnergy<T>(lattice)
<< "; av rho ="
<< getStoredAverageDensity<T>(lattice) << endl;
}
}

delete boundaryCondition;
``````

}

My understanding is that with boundaryCondition.setVelocityConditionOnBlockBoundaries you also define values on edges and periodic boundary condition does not like it : it crashes. However if you use a BC without edges like

``````

``````

then it runs with periodic BC

``````
lattice.periodicity().toggle(0, true);
lattice.periodicity().toggle(1, true);

``````

which makes sense since threre is no edges if I have periodic BC.

Now the trouble is that even with this trick, results have trouble. If someone knows why, please let me know.

Here is my 3D channel with no-slip BC on the bottom, slip on the top and periodic in x and y (z is vertical direction) with density BC for in and out.

``````
#include "palabos3D.h"
#ifndef PLB_PRECOMPILED // Unless precompiled version is used,
#include "palabos3D.hh"   // include full template code
#endif
#include <vector>
#include <iostream>
#include <iomanip>

using namespace plb;
using namespace std;

typedef double T;
#define DESCRIPTOR plb::descriptors::D3Q19Descriptor
#define DYNAMICS BGKdynamics

/// Linearly decreasing pressure profile
T poiseuillePressure(plint iX, IncomprFlowParam<T> const& parameters) {
T Lx = parameters.getNx()-1;
T Ly = parameters.getNy()-1;
T Lz = parameters.getNz()-1;
return 8.*parameters.getLatticeNu()*parameters.getLatticeU() / (Lz*Lz) * (Lx/(T)2-(T)iX);
}

/// Convert pressure to density according to ideal gas law
T poiseuilleDensity(plint iX, IncomprFlowParam<T> const& parameters) {
return poiseuillePressure(iX,parameters)*DESCRIPTOR<T>::invCs2 + (T)1;
}

/// A functional, used to initialize the density for the boundary conditions
template<typename T>
class PoiseuilleDensity {
public:
PoiseuilleDensity(IncomprFlowParam<T> parameters_)
: parameters(parameters_)
{ }
T operator()(plint iX, plint iY,plint iZ) const {
return poiseuilleDensity(iX,parameters);
}
private:
IncomprFlowParam<T> parameters;
};

/// A functional, used to create an initial condition for with zero velocity,
///   and linearly decreasing pressure.
template<typename T>
class PoiseuilleDensityAndZeroVelocity {
public:
PoiseuilleDensityAndZeroVelocity(IncomprFlowParam<T> parameters_)
: parameters(parameters_)
{ }
void operator()(plint iX, plint iY, plint iZ, T& rho, Array<T,3>& u) const {
rho = poiseuilleDensity(iX,parameters);
u[0] = T();
u[1] = T();
u[2] = T();
}
private:
IncomprFlowParam<T> parameters;
};

void channelSetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
IncomprFlowParam<T> const& parameters,
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition )
{
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();
const plint N = parameters.getResolution();

Box3D top(0, nx-1, 0, ny-1, nz-1, nz-1);
Box3D bottom(0, nx-1, 0, ny-1, 0, 0);
Box3D in(0,0,0,ny-1,0,nz-1);
Box3D out(nx-1,nx-1,0,ny-1,0,nz-1);
Box3D in_edge(0, 0, 0, ny-2, nz-1, nz-1);
Box3D out_edge(nx-1, nx-1, 0, ny-1, nz-1, nz-1);

// Create pressure boundary conditions
boundaryCondition.setPressureConditionOnBlockBoundaries(lattice,in);
boundaryCondition.setPressureConditionOnBlockBoundaries(lattice,out);

setBoundaryDensity (
lattice, lattice.getBoundingBox(),
PoiseuilleDensity<T>(parameters) );

// velocity

initializeAtEquilibrium (
lattice, lattice.getBoundingBox(),
PoiseuilleDensityAndZeroVelocity<T>(parameters) );

lattice.initialize();
}

void writeGifs(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
{
const plint imSize = 600;
const plint nx = lattice.getNx();
const plint ny = lattice.getNy();
const plint nz = lattice.getNz();
Box3D slice(0, nx-1, ny/2, ny/2, 0, nz-1);
ImageWriter<T> imageWriter("leeloo");
imageWriter.writeScaledGif(createFileName("u", iter, 6),
*computeVelocityNorm(lattice, slice),
imSize, imSize );
imageWriter.writeScaledGif(createFileName("d", iter, 6),
*computeDensity(lattice, slice),
imSize, imSize );
}
int main(int argc, char* argv[]) {
plbInit(&argc, &argv);

global::directories().setOutputDir("./tmp/");

// Use the class IncomprFlowParam to convert from
//   dimensionless variables to lattice units, in the
//   context of incompressible flows.
IncomprFlowParam<T> parameters(
(T) 1e-2,  // Reference velocity (the maximum velocity
//   in the Poiseuille profile) in lattice units.
(T) 100.,  // Reynolds number
50,        // Resolution of the reference length (channel height).
3.,        // Channel length in dimensionless variables
1.,        // Channel height in dimensionless variables
1.         // Channel depth in dimensionless variables
);
const T imSave   = (T)0.02; // Time intervals at which to save GIF
//   images, in dimensionless time units.
const T maxT     = (T)2.2;  // Total simulation time, in dimensionless
//   time units.
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();

writeLogFile(parameters, "3D Poiseuille flow");

MultiBlockLattice3D<T, DESCRIPTOR> lattice (
parameters.getNx(), parameters.getNy(), parameters.getNz(),
new DYNAMICS<T,DESCRIPTOR>(parameters.getOmega()) );

lattice.periodicity().toggle(0, true);
lattice.periodicity().toggle(1, true);

OnLatticeBoundaryCondition3D<T,DESCRIPTOR>*
boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();

channelSetup(lattice, parameters, *boundaryCondition);

// Main loop over time iterations.
for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
if (iT%parameters.nStep(imSave)==0) {
pcout << "Saving Gif at time step " << iT << endl;
writeGifs(lattice, iT);
}
// Execute lattice Boltzmann iteration.
lattice.collideAndStream();
}

delete boundaryCondition;
}

``````

Thank you Olivier for your explanation, the validity of which I can confirm. Indeed, in the current version of Palabos, if you use the “setXXConditionOnBlockBoundaries” function to create the boundary condition, then the code is incompatible with periodic boundaries. The reason is that Palabos uses extrapolations on edges and corners (with nearest-neighbor access); under the assumption of a periodic domain, the code gets confused and tries to access non-existing neighbors. This is an issue, and we’ll try to fix it as soon as possible.

In the meantime, Olivier’s workaround offers a neat way of getting a working code. I’m just getting a bit confused with Olivier’s code, as it starts by using the “setPressureConditionOnBlockBoundaries” function (which as we said should, for now, not be used in this context), before invoking “addVelocityBoundary2P”. Anyway, here’s an adaptation of Olivier’s code which seems to work just fine (I created a little obstacle in order to get an interesting flow pattern).

``````
#include "palabos3D.h"
#ifndef PLB_PRECOMPILED // Unless precompiled version is used,
#include "palabos3D.hh"   // include full template code
#endif
#include <vector>
#include <iostream>
#include <iomanip>

using namespace plb;
using namespace std;

typedef double T;
#define DESCRIPTOR plb::descriptors::D3Q19Descriptor
#define DYNAMICS BGKdynamics

/// Linearly decreasing pressure profile
T poiseuillePressure(plint iX, IncomprFlowParam<T> const& parameters) {
T Lx = parameters.getNx()-1;
T Ly = parameters.getNy()-1;
T Lz = parameters.getNz()-1;
return 8.*parameters.getLatticeNu()*parameters.getLatticeU() / (Lz*Lz) * (Lx/(T)2-(T)iX);
}

/// Convert pressure to density according to ideal gas law
T poiseuilleDensity(plint iX, IncomprFlowParam<T> const& parameters) {
return poiseuillePressure(iX,parameters)*DESCRIPTOR<T>::invCs2 + (T)1;
}

/// A functional, used to initialize the density for the boundary conditions
template<typename T>
class PoiseuilleDensity {
public:
PoiseuilleDensity(IncomprFlowParam<T> parameters_)
: parameters(parameters_)
{ }
T operator()(plint iX, plint iY,plint iZ) const {
return poiseuilleDensity(iX,parameters);
}
private:
IncomprFlowParam<T> parameters;
};

/// A functional, used to create an initial condition for with zero velocity,
///   and linearly decreasing pressure.
template<typename T>
class PoiseuilleDensityAndZeroVelocity {
public:
PoiseuilleDensityAndZeroVelocity(IncomprFlowParam<T> parameters_)
: parameters(parameters_)
{ }
void operator()(plint iX, plint iY, plint iZ, T& rho, Array<T,3>& u) const {
rho = poiseuilleDensity(iX,parameters);
u[0] = T();
u[1] = T();
u[2] = T();
}
private:
IncomprFlowParam<T> parameters;
};

void channelSetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
IncomprFlowParam<T> const& parameters,
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition )
{
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();
const plint N = parameters.getResolution();

Box3D in(0,0,0,ny-1,0,nz-1);
Box3D out(nx-1,nx-1,0,ny-1,0,nz-1);

defineDynamics(lattice, Box3D(nx/4,nx/4, ny/2-2, ny/2+5, nz/2-2, nz/2+5),
new BounceBack<T,DESCRIPTOR>());

// Create pressure boundary conditions

setBoundaryDensity (
lattice, lattice.getBoundingBox(),
PoiseuilleDensity<T>(parameters) );

initializeAtEquilibrium (
lattice, lattice.getBoundingBox(),
PoiseuilleDensityAndZeroVelocity<T>(parameters) );

lattice.initialize();
}

void writeGifs(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
{
const plint imSize = 600;
const plint nx = lattice.getNx();
const plint ny = lattice.getNy();
const plint nz = lattice.getNz();
Box3D slice(0, nx-1, ny/2, ny/2, 0, nz-1);
ImageWriter<T> imageWriter("leeloo");
imageWriter.writeScaledGif(createFileName("u", iter, 6),
*computeVelocityNorm(lattice, slice),
imSize, imSize );
imageWriter.writeScaledGif(createFileName("d", iter, 6),
*computeDensity(lattice, slice),
imSize, imSize );
}
int main(int argc, char* argv[]) {
plbInit(&argc, &argv);

global::directories().setOutputDir("./tmp/");

// Use the class IncomprFlowParam to convert from
//   dimensionless variables to lattice units, in the
//   context of incompressible flows.
IncomprFlowParam<T> parameters(
(T) 1e-2,  // Reference velocity (the maximum velocity
//   in the Poiseuille profile) in lattice units.
(T) 200.,  // Reynolds number
30,        // Resolution of the reference length (channel height).
4.,        // Channel length in dimensionless variables
1.,        // Channel height in dimensionless variables
1.         // Channel depth in dimensionless variables
);
const T imSave   = (T)0.2; // Time intervals at which to save GIF
//   images, in dimensionless time units.
const T maxT     = (T)10.2; // Total simulation time, in dimensionless
//   time units.
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();

writeLogFile(parameters, "3D Poiseuille flow");

MultiBlockLattice3D<T, DESCRIPTOR> lattice (
parameters.getNx(), parameters.getNy(), parameters.getNz(),
new DYNAMICS<T,DESCRIPTOR>(parameters.getOmega()) );

lattice.periodicity().toggle(1, true);
lattice.periodicity().toggle(2, true);

OnLatticeBoundaryCondition3D<T,DESCRIPTOR>*
boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();

channelSetup(lattice, parameters, *boundaryCondition);

// Main loop over time iterations.
for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
if (iT%parameters.nStep(imSave)==0) {
pcout << "Saving Gif at time step " << iT << endl;
writeGifs(lattice, iT);
}
// Execute lattice Boltzmann iteration.
lattice.collideAndStream();
}

delete boundaryCondition;
}

``````

Thanks Jonas for your example. Could you make it with a no slip BC on the top and periodic BC in x and y. Here what I tried but it has problems at corners. They reduce slowly with time steps but I guess we should be able to avoid them from the beginning.

``````
Box3D top(0,nx-1,0,ny-1,nz-1,nz-1);
Box3D in(0,0,0,ny-1,0,nz-2);               // with nz-1 it is worst
Box3D out(nx-1,nx-1,0,ny-1,0,nz-2);

``````

and

``````
lattice.periodicity().toggle(0, true);
lattice.periodicity().toggle(1, true);

``````

TIA

Olivier.

Hi Olivier,

I don’t precisely understand the setup you want to get. Is it freeslip on top and noslip on bottom? And periodicity in x- and y-direction? But then, why do you define the regions “in” and “out”, if the x-direction is simply periodic?

Hi Johan,

I need something to move the fuild. I find a difference of pression more suitable than imposed velocity. A body force would be the best, but up to now I didn’t succeded to make it work.

Olivier.

hie there, i have come up with a code for lamina flow in a pipe and i need to validate my results. my code is giving me very inflated results . please i need help in understanding where i went wrong with my dimensions and boundary conditions.

here is my code:

#include “palabos3D.h”
#include “palabos3D.hh”
#include <math.h>
#include
#include

using namespace plb;
typedef double T;

#define DESCRIPTOR descriptors::D3Q19Descriptor

static std::string outDir("./tmp/");

void terminalOutput()
{
pcout << “Flow in a lateral pipe simulation…” << std::endl << “\n”;

pcout << "Creating the geometry. " << std::endl << “\n”;

}

template
class BuildCylinder {
public:
BuildCylinder(plint ny_, plint nz_)
: ny(ny_), nz(nz_)

``````{ }

int operator()(plint iX, plint iY, plint iZ)
{
if ((iY-(ny-1)/2)*(iY-(ny-1)/2)+(iZ-(nz-1)/2)*(iZ-(nz-1)/2) >= (ny-1)*(ny-1)/4){return 1;}
return 0;
}
``````

private:
plint ny;
plint nz;
};

int main(int argc, char* argv[]) {

``````	plbInit(&argc, &argv);

terminalOutput();

T rho_LB = 1.;
T lx = 10e-3;
T ly = 1e-3;
T lz = 1e-3;               // Size of computational domain.

T nu = 1e-6;
T pTau = 1;
T nu_LB = (pTau - 0.5)*DESCRIPTOR<T>::cs2;
T rho = 1000;
T resolution = 20;
T dx = ly/(resolution-1);
T dt = nu_LB/nu * (dx * dx);

plint nx = util::roundToInt(lx / dx) + 1;
plint ny = util::roundToInt(ly / dx) + 1;
plint nz = util::roundToInt(lz / dx) + 1;

//Array<T,3> inletVelocity = Array<T,3>(0.2, 0., 0.);

MultiScalarField3D<int> *solid = new MultiScalarField3D<int>(nx, ny, nz, 0);
setToFunction(*solid, solid->getBoundingBox(), BuildCylinder<T>(ny, nz));

pcout << "Creating pipe image output" << std::endl << "\n";

VtkImageOutput3D<T> vtkOut("./tmp/pipe", 1);
vtkOut.writeData<float>(*copyConvert<int,double>(*solid, solid->getBoundingBox()), "tag", 1.0);

// Create Object that holds distribution function + dynamics

T cSmago = 0.; // 0 = laminar, 0.14 = turbulent
Dynamics<T,DESCRIPTOR> *Dynamics = 0;
T omega = 1;

if (util::isZero(cSmago)){Dynamics = new RegularizedBGKdynamics<T,DESCRIPTOR>(omega);}
else {Dynamics = new SmagorinskyBGKdynamics<T,DESCRIPTOR>(omega, cSmago);}

MultiBlockLattice3D<T,DESCRIPTOR> *Lattice = new MultiBlockLattice3D<T,DESCRIPTOR>((MultiBlock3D&) *solid);
defineDynamics(*Lattice, Lattice->getBoundingBox(), Dynamics->clone());
delete Dynamics;

pcout << "Computing boundary conditions" << std::endl << "\n";
// Boundary Conditions

OnLatticeBoundaryCondition3D<T,DESCRIPTOR> *bound = createLocalBoundaryCondition3D<T,DESCRIPTOR>();
bound->setVelocityConditionOnBlockBoundaries(*Lattice, Box3D(0., 0., 0., ny-1, 0., nz-1), boundary::dirichlet);
T inletVelocity = 2e-2;
Array<T,3> inletVelocity_LB = Array<T,3>(inletVelocity*(dt/dx), 0., 0.); // in lattice units
setBoundaryVelocity(*Lattice, Box3D(0., 0. , 0., ny-1, 0., nz-1), inletVelocity_LB);
bound->setVelocityConditionOnBlockBoundaries(*Lattice, Box3D(nx-1, nx-1, 0., ny-1, 0., nz-1), boundary::neumann);
defineDynamics(*Lattice, *solid, Lattice->getBoundingBox(), new BounceBack<T,DESCRIPTOR>(rho_LB), 1);

pcout << "Running the simulation" << std::endl << "\n";
// Initialize (i.e. perform one streaming step)
initializeAtEquilibrium(*Lattice, Lattice->getBoundingBox(), rho_LB, inletVelocity_LB);
Lattice->initialize();

//writeVTK(*Lattice, *bound, i, outDir);

for (plint i = 0; i <= 10000; i++)
{
Lattice->collideAndStream();
if (i % 100 == 0){
std::string fname = createFileName("./tmp/Domain_", i, 8);
VtkImageOutput3D<T> vtkOut(fname, 1);
std::unique_ptr<MultiTensorField3D<T,3> > v = computeVelocity(*Lattice, Lattice->getBoundingBox());
vtkOut.writeData<3,float>(*v, "velocity", (dx/dt));
std::unique_ptr<MultiScalarField3D<T> > p = computeDensity(*Lattice, Lattice->getBoundingBox());
vtkOut.writeData<float>(*p, "pressure", (dx*dx/dt*dt)*rho);}
}
}``````