Sudden Expansion flow

I am a new user of palabos.
I try to solve sudden Expansion flow problem.
In order to solve it I took the tutorial code of Poiseuille flow and changed the geometry.
However when running the simulation it seems that the walls are not solid and I still recieve a tube tunnel in the result.
The code is attached and I will be greatfull for any asistance.
Thank you.

#include “palabos2D.h”
#include “palabos2D.hh”

/* Code 1.5 in the Palabos tutorial

using namespace plb;
using namespace std;

typedef double T;
#define DESCRIPTOR plb::descriptors::D2Q9Descriptor

const int nx = 1200;
const int ny = 400;

/// Velocity on the parabolic Poiseuille profile
T poiseuilleVelocity(plint iY, IncomprFlowParam const& parameters) {
T y = (T)iY / parameters.getResolution();
return 4.parameters.getLatticeU() * (y-yy);

/// A functional, used to initialize the velocity for the boundary conditions
class PoiseuilleVelocity {
PoiseuilleVelocity(IncomprFlowParam parameters_)
: parameters(parameters_)
{ }
/// This version of the operator returns the velocity only,
/// to instantiate the boundary condition.
void operator()(plint iX, plint iY, Array<T,2>& u) const {
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
/// This version of the operator returns also a constant value for
/// the density, to create the initial condition.
void operator()(plint iX, plint iY, T& rho, Array<T,2>& u) const {
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
rho = (T)1;
IncomprFlowParam parameters;

void channelSetup (
MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters,
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )

Box2D topLow(0, nx/2-1, ny*5/8-1, ny*5/8-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, topLow, boundary::freeslip );
Box2D topSide(nx/2-1, nx/2-1, ny*5/8-1 , ny-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, topSide, boundary::freeslip );
Box2D topWall(nx/2-1, nx-1,  ny-1, ny-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, topWall, boundary::freeslip );
Box2D bottomLow(0, nx/2-1, ny*3/8-1, ny*3/8-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, bottomLow, boundary::freeslip );
Box2D bottomSide(nx/2-1, nx/2-1, ny*3/8-1, 0);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, bottomSide, boundary::freeslip );
Box2D bottomWall(nx/2-1, nx-1, 0, 0);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, bottomWall, boundary::freeslip );
Box2D inlet(0, 0, ny*3/8+2, ny*5/8+2);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, inlet );
Box2D outlet(nx-1, nx-1, 2, ny-2);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, outlet); //, boundary::outflow );
boundaryCondition.addInternalVelocityCornerPN( nx/2, ny*5/8, lattice, boundary::freeslip );
boundaryCondition.addInternalVelocityCornerPP( nx/2, ny*3/8, lattice, boundary::freeslip );
boundaryCondition.addExternalVelocityCornerPN( nx/2, ny, lattice, boundary::freeslip );
boundaryCondition.addExternalVelocityCornerPP( nx/2, 0, lattice, boundary::freeslip );

// Specify the boundary velocity.
setBoundaryVelocity (
        lattice, lattice.getBoundingBox(),
        PoiseuilleVelocity<T>(parameters) );

// Create the initial condition.
initializeAtEquilibrium (
       lattice, lattice.getBoundingBox(), PoiseuilleVelocity<T>(parameters) );



void writeGifs(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter)
const plint imSize = 600;
ImageWriter imageWriter(“leeloo”);
imageWriter.writeScaledGif(createFileName(“u”, iter, 6),
imSize, imSize );

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


// Use the class IncomprFlowParam to convert from
//   dimensionless variables to lattice units, in the
//   context of incompressible flows.
IncomprFlowParam<T> parameters(
        (T) 1e-1,  // Reference velocity (the maximum velocity
                   //   in the Poiseuille profile) in lattice units.
        (T) 100.,  // Reynolds number
        300,       // Resolution of the reference length (channel height).
        3.,        // Channel length in dimensionless variables
        1.         // Channel height in dimensionless variables
const T imSave   = (T)0.1;  // Time intervals at which to save GIF
                            //   images, in dimensionless time units.
const T maxT     = (T)3.1;  // Total simulation time, in dimensionless
                            //   time units.

writeLogFile(parameters, "Poiseuille flow");

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

    boundaryCondition = createLocalBoundaryCondition2D<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 && iT>0) {
        pcout << "Saving Gif at time step " << iT << endl;
        writeGifs(lattice, iT);
    // Execute lattice Boltzmann iteration.

delete boundaryCondition;


please… someone…
I just want to change the geometry from box to those shown in the picture.

And I don’t understand why I dont get “solid walls” with the code that I wrote.

Thank you in advence…

// Velocity boundary condition on bottom wall.
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(1, nx/8-1, ny3/8, ny3/8) );
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(nx/8-1, nx/8-1, ny3/8, 1) );
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(nx/8-1, nx-1, 1, 1) );
// Velocity boundary condition on top wall.
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(1, nx/8-1, ny
5/8, ny5/8) );
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(nx/8-1, nx/8-1, ny
5/8, ny-1) );
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(nx/8-1, nx-1, ny-1, ny-1) );

boundaryCondition.addInternalVelocityCornerPN( nx/8-1, ny*3/8, lattice );
boundaryCondition.addInternalVelocityCornerPP( nx/8-1, ny*5/8, lattice );
boundaryCondition.addExternalVelocityCornerPP( nx/8-1, 1, lattice, boundary::freeslip );
boundaryCondition.addExternalVelocityCornerPN( nx/8-1, ny-1, lattice, boundary::freeslip );