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.

The code:

#include "palabos2D.h"
#ifndef PLB_PRECOMPILED // Unless precompiled version is used,
#include "palabos2D.hh"   // include full template code
#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 {
    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();
    IncomprFlowParam<T> parameters;

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

/// A functional, used to create an initial condition for with zero velocity,
///   and linearly decreasing pressure.
template<typename T>
class PoiseuilleDensityAndZeroVelocity {
    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();
    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.
                            lattice.getComponent(iLevel).getBoundingBox(), PoiseuilleVelocity<T>(parameters.getParameters(iLevel) ) );
		//Define the value of the imposed density at the outlet
						   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),
								 PoiseuilleDensityAndZeroVelocity<T>(parameters.getParameters(iLevel) ) );

/// 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::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);
    // adding refinements to management
    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;
    for (plint iBlock=0; iBlock<(plint)ranges.size()-1; ++iBlock){
	boundaryCondition = createLocalBoundaryCondition2D<T,DESCRIPTOR>();
	lattice( management, 
			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.
        // 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;


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.


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.

Thank you for your help.


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:

Velocity y-direction with grid refinement (only level 1, data range as without gridrefinment):

Velocity y-direction with grid refinement (only level 1, data range adapted data range):

Velocity y-direction with grid refinement (both levels of the grid):