# Poisseuille current - grid refinement

Hello everybody,

I combined the grid refinement with the poiseuille showcase.
I can run a simulation with a poiseuille velocity profile and a grid refinement. When I look at the solution I get a velocity in y-direction in the area of the corners of the grid refinement.

As I cannot find a mistake in the code, I would be glad if you take a look at it.

Thanks a lot.
Cheers
Konstantin

The code:

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

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

typedef double T;
#define DESCRIPTOR D2Q9Descriptor

/// Velocity on the parabolic Poiseuille profile
// With this formula the width of the channel has to be 1
T poiseuilleVelocity(plint iY, IncomprFlowParam<T> const& parameters) {
T y = (T)iY / parameters.getResolution();
return 4.*parameters.getLatticeU() * (y-y*y);
}

/// Linearly decreasing pressure profile
T poiseuillePressure(plint iX, IncomprFlowParam<T> const& parameters) {
T Lx = parameters.getNx()-1;
T Ly = parameters.getNy()-1;
return 8.*parameters.getLatticeNu()*parameters.getLatticeU() / (Ly*Ly) * (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 velocity for the boundary conditions
template<typename T>
class PoiseuilleVelocity {
public:
PoiseuilleVelocity(IncomprFlowParam<T> parameters_)
: parameters(parameters_)
{ }
void operator()(plint iX, plint iY, Array<T,2>& u) const {
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
}
private:
IncomprFlowParam<T> parameters;
};

/// 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) 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, T& rho, Array<T,2>& u) const {
rho = poiseuilleDensity(iX,parameters);
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
}
private:
IncomprFlowParam<T> parameters;
};

// Initialisation of the domain.
/// A functional, used to instantiate bounce-back nodes at the locations of the cylinder

//Change to MultiGridLattice2D
//+Convective Refinement Parameters - data wrapper ConvectiveRefinementParameters instead of IncomprFlowParam<T>
void cylinderSetup( MultiGridLattice2D<T,DESCRIPTOR>& lattice,
ConvectiveRefinementParameters<T> const& parameters,
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )
{
//for each lattice inside the multigrid
for (plint iLevel=0; iLevel<lattice.getNumLevels(); ++iLevel){
const plint nx = parameters.getParameters(iLevel).getNx();
const plint ny = parameters.getParameters(iLevel).getNy();

Box2D inlet(0, 0, 2, ny-2);
Box2D outlet(nx-1, nx-1, 1, ny-2);

//Velocity boundary condition on the bottom wall
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice.getComponent(iLevel), Box2D(0, nx-1, 0, 0) );
//Velocity boundary condition on the top wall
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice.getComponent(iLevel), Box2D(0, nx-1, ny-1, ny-1) );

//Velocity boundary condition at the inlet
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice.getComponent(iLevel), inlet );
//Pressure boundary condition at the outlet
boundaryCondition.setPressureConditionOnBlockBoundaries (lattice.getComponent(iLevel), outlet );

//Define the value of the imposed velocity on all nodes which have previously been
//	defined to be velocity boundary nodes.
setBoundaryVelocity(lattice.getComponent(iLevel),
lattice.getComponent(iLevel).getBoundingBox(), PoiseuilleVelocity<T>(parameters.getParameters(iLevel) ) );

//Define the value of the imposed density at the outlet
setBoundaryDensity(lattice.getComponent(iLevel),
Box2D(nx-1, nx-1, 1, ny-2), PoiseuilleDensity<T>(parameters.getParameters(iLevel) ) );

// Initialize all cells at an equilibrium distribution, with a velocity and density
//   value of the analytical Poiseuille solution.
initializeAtEquilibrium (lattice.getComponent(iLevel),
lattice.getComponent(iLevel).getBoundingBox(),
PoiseuilleDensityAndZeroVelocity<T>(parameters.getParameters(iLevel) ) );

}

lattice.initialize();

}

/// Produce a VTK snapshot of the velocity-norm.
void writeVTKs(MultiBlockLattice2D<T,DESCRIPTOR> lattice, IncomprFlowParam<T> const& parameters, plint iter, string fName)
{
T dx = parameters.getDeltaX();
T dt = parameters.getDeltaT();

VtkImageOutput2D<T> vtkOut(createFileName(fName, iter, 6), dx);
vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", dx/dt);
vtkOut.writeData<2,float>(*computeVelocity(lattice), "velocity", dx/dt);
vtkOut.writeData<float>(*computeDensity(lattice), "density", dx/dt);

}

template<typename T, template<typename U> class Descriptor>
void writeVTKs(MultiGridLattice2D<T,Descriptor> const& lattice, ConvectiveRefinementParameters<T> const& parameters, plint iter) {
for (plint iLevel=0; iLevel<lattice.getNumLevels(); ++iLevel) {
std::stringstream fname("level_");
fname << "level_" << iLevel << "_";
writeVTKs(lattice.getComponent(iLevel), parameters.getParameters(iLevel), iter, fname.str());
}
}

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

global::directories().setOutputDir("./tmp/");
//NEW
global::setDefaultMultiScaleManager(new ConvectiveMultiScaleManager());

plint referenceLevel = 0;
plint behaviorLevel = referenceLevel;
plint numLevel = 2;

plint N=100;

IncomprFlowParam<T> baseParameters(
(T) 0.001,  // uMax
(T) 100.,  // Re
N,       // N
6.,        // lx
1.         // ly
);

const T logT     = (T)0.02;
const T imSave   = (T)1000.1;
const T vtkSave  = (T)0.05;
const T maxT     = (T)5.001;

writeLogFile(baseParameters, "GridRefinement and Initialization with the poiseuille profil");

// we choose a convective refinement dx=dt
ConvectiveRefinementParameters<T> convectiveParameters (numLevel,referenceLevel,baseParameters );

Dynamics<T,DESCRIPTOR>* backgroundDynamics = new BGKdynamics<T,DESCRIPTOR>(baseParameters.getOmega());

// main object for the grid refinement
MultiGridManagement2D management(baseParameters.getNx(),baseParameters.getNy(), numLevel, referenceLevel);

management.refine(0,Box2D(2*N,4*N,N/3,2*N/3)); // Block in the middle (lenghth 2N - height N/3)

// Parallelization. Here, the domain is parallelized along the x-direction. The position
//   of the different slices is explicitly created, and provided in units of the
//   finest lattice.
std::vector<plint> parallelizeIn;
plint procNumber = global::mpi().getSize();
std::vector<std::pair<plint,plint> > ranges;
plint finestNy = convectiveParameters[numLevel-1].getNy()-1;
util::linearRepartition(0,finestNy,procNumber,ranges);
for (plint iBlock=0; iBlock<(plint)ranges.size()-1; ++iBlock){
parallelizeIn.push_back(ranges[iBlock].second);
}

management.parallelizeX(parallelizeIn);

OnLatticeBoundaryCondition2D<T,DESCRIPTOR>*
boundaryCondition = createLocalBoundaryCondition2D<T,DESCRIPTOR>();

MultiGridLattice2D<T,DESCRIPTOR>
lattice( management,
defaultMultiGridPolicy2D().getBlockCommunicator<T>(management.getNumLevels()),
defaultMultiGridPolicy2D().getCombinedStatistics(management.getNumLevels()),
backgroundDynamics, behaviorLevel );

cylinderSetup(lattice, convectiveParameters, *boundaryCondition);

// output file
plb_ofstream outfile("out.dat");

// Main loop over time iterations.
for (plint iT=0; iT*baseParameters.getDeltaT()<maxT; ++iT) {

if (iT%baseParameters.nStep(vtkSave)==0) {
pcout << "Saving VTK ..." << endl;
writeVTKs(lattice, convectiveParameters, iT);
pcout << endl;
}

if (iT%baseParameters.nStep(logT)==0) {
pcout << "step " << iT
<< "; t=" << iT*baseParameters.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%baseParameters.nStep(logT)==0) {
pcout << "; av energy ="
<< setprecision(10) << getStoredAverageEnergy<T>(lattice)
<< "; av rho ="
<< getStoredAverageDensity<T>(lattice) << endl;
pcout << "; Omega =" << baseParameters.getOmega() << endl;
}

}

delete boundaryCondition;
}

``````

Hello,

I signal that you have an error while declaring the inlet for the BC:

``````Box2D inlet(0, 0, 2, ny-2);
``````

should be

``````Box2D inlet(0, 0, 1, ny-2);
``````

I am verifying if there is no bug with the grid refinement, but it seems to me that even without grid refinement you are going to have a velocity in the y direction.

Cheers,
Daniel

Hi Daniel,

thank you very much for your answer. We are aware of the velocity in y-direction in the area of the inlet and the outlet. We observe this also without grid refinement.
But there is also a velocity in the corners of the grid refinement. I added some images with and without grid refinement.
This result doesn’t change with more iterations.

Cheers,
Konstantin

Edit: As you can’t see the data range on the images here. I add the links directly to the images:
Velocity y-direction without grid refinement: