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_),
radiusSqr(plb::util::sqr(radius))
{ }
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;
plb::plint radiusSqr;
};
/// 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.;
plint radius = 10.;
defineDynamics(lattice, lattice.getBoundingBox(),
new CylinderShapeDomain3D<T>(cx,cy,cz,radius),
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;
}