Hi Everyone,

I finally managed to solve (almost) the advection-diffusion problem. I also managed to couple the concentration field with a velocity field. I thought that maybe can be useful to somebody and I posted the code here. The code solve the concentration field within a poiseuille flow. The coupling is passive and uses a function from the palabos library.

[code=“cpp”]
#include
#include “palabos2D.h”
#include “palabos2D.hh”

using namespace std;
using namespace plb;

typedef double T;

#define NSDES descriptors::D2Q9Descriptor

// ************************************************************************
// Part for the NS part of the code. A classic poiseuille flow is simulated.
// Most of the code is taken from Palabos tutorial 1.5
// ************************************************************************

// Poiselle velocity analitical
// 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
template
class PoiseuilleVelocity {
public:
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;
}
private:
IncomprFlowParam parameters;
};

// Function to set up the channel
void channelSetup (
MultiBlockLattice2D<T,NSDES>& lattice,
IncomprFlowParam const& parameters,
OnLatticeBoundaryCondition2D<T,NSDES>& boundaryCondition )
{
// Create Velocity boundary conditions.
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);

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

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

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

}
// ****************************************************************************
// ****************************************************************************

// ************************************************************************
// ************************************************************************

void concentrationSetup(MultiBlockLattice2D<T,AVDES> &lattice,
IncomprFlowParam parameters) {

``````// Get lattice dimensions
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();

// Get velocity and rho for this lattice
Array<T,2> u0(0.,0.);

// Set-up the geometry
//Box2D topwall(0, nx-1, ny-1, ny-1);
//Box2D bottomwall(0, nx-1, 0,    0);
Box2D leftwall(0,1,   0,ny);
Box2D rightwall(nx,nx,   0,ny);

// Set the boundary-conditions

// Try to impose a different temperature
setBoundaryDensity(lattice, rightwall, 0.0);
setBoundaryDensity(lattice, leftwall, 1.0);

// Init lattice
initializeAtEquilibrium(lattice, lattice.getBoundingBox(), 0.0,u0);

// Command to init the lattice
lattice.periodicity().toggle(1,true);
lattice.initialize();
``````

}

// ************************************************************************
// ************************************************************************
// ************************************************************************

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

``````// Parameters for the dynamics
IncomprFlowParam<T> parameters(
(T) 1e-2,  // Reference velocity (the maximum velocity
//   in the Poiseuille profile) in lattice units.
(T) 100.,  // Reynolds number
100,       // Resolution of the reference length (channel height).
2.,        // Channel length in dimensionless variables
1.         // Channel height in dimensionless variables
);

Array<T,2> u2(2.0, 0.0);

// Istantiate the lattice
MultiBlockLattice2D<T, NSDES> nsLattice(parameters.getNx(), parameters.getNy(), new BGKdynamics<T, NSDES>(parameters.getOmega()));

// Istantiate the boundary conditions
OnLatticeBoundaryCondition2D<T,NSDES> *nsBC = createLocalBoundaryCondition2D<T, NSDES>();

// Set-up channel for NS equations
channelSetup(nsLattice, parameters, *nsBC);

//Set-up concentration field

// Main LBM cycle
ImageWriter<T> image("leeloo");
for (plint iT=0; iT < 10000; ++iT) {

if(iT%100==0) {
pcout << "Writing GIF... \n";
//image.writeScaledGif(createFileName("vel",iT, 2), *computeVelocityNorm(nsLattice), parameters.getNx(), parameters.getNy());
}
nsLattice.collideAndStream();
}//*/

// First attempt to couple the two physics

for (plint iT=0; iT < 10000; ++iT) {

if(iT%100==0) {
pcout << "Writing GIF... \n";
//image.writeScaledGif(createFileName("vel",iT, 2), *computeVelocityNorm(nsLattice), parameters.getNx(), parameters.getNy());
}
}//*/

pcout << "Hello palabos";
return 0;
``````

}

``````

The last step will be implement the advection diffusion with source. If any of you have some suggestions they will be appreciated.

Cheers,

Salvo``````
2 Likes

Hi!

I just want to thank you for sharing this example! I want to model a multicomponent advection-diffusion problem, and I just got started with Palabos. On the library examples directory the only problem that employs diffusion is the Boussinesq example, and it isn’t of much help. Your program is a good starting point for me!

Cheers,
Eduardo