Simulating a canal with sliding gate


I’m trying to write a simple program in Palabos, that would simulate a canal filled with water, with a sliding gate in the middle.
It involves changing the collision and streaming functions for the part with the gate. I have a formula, that allow me to compute the distribution functions in that place (in all papers I read it was defined as f) and I could write it in matlab or simple C++, but I want to use the Palabos library for that purpose.
From what I read in the user’s guide, it should be possible with use of D2Q9 model and some data processor, but I have no idea which and how.

Could you help me with the task?


in what way do you want to change the collision and streaming. Unless you provide some additional details it will be difficult to help.

Thank you for your reply.

First some terms:
f - distribution function
fin - pre-collision distribution (what flows in in time t)
fout - post-collision distribution (what will flow out after time t)
feq - local equilibrium distribution function

The collision in LBGK model is computed using the formula: fout = fin + 1/tau * (feq-fin)

And then streaming is basically fin = fout with imposed proper boundary conditions.

If we assume fin(i,j,k) means pre-collision distribution at point i[/i] of the lattice, in direction k, water height h at point i[/i] is the sum of fin(i,j,k) for every k (in D2Q9 case it’d be k=0,1,2,…,8)
h(i,j) = sum(f(i,j,:))

It’s nice and easy equation for regular canals, but if we want to add an obstacle in the middle of it, some values of f change. What I’d like to know is if it’s possible to directly affect these values during the collision phase, so they’re computed using different formula

I’m not sure to understand in what way the collision-streaming steps are changed. But if you have a purely local operation to do, you can define a new collision (you can inspire yourself by the collisions available in the basicDynamics/isoThermalDynamics.h*) instead of using a dataProcessor. To then use the define dynamics using

defineDynamics(lattice, domain, new dynamics<T,DESCRIPTOR>);

I hope it helps,

Thank you, I’ll look into it. Yes, it’s purely local.

How can I directly access distribution function?

I want to apply these formulas on the left side of the gate and very similar ones on the right side:

f[sub]1[/sub] = f[sub]5[/sub] + 2hu[sub]x[/sub]/3v
f[sub]2[/sub] = hu[sub]x[/sub]/6v + f[sub]6[/sub] + (f[sub]7[/sub]-f[sub]3[/sub])/2
f[sub]8[/sub] = hu[sub]x[/sub]/6v + f[sub]4[/sub] + (f[sub]3[/sub]-f[sub]7[/sub])/2

where u is the flow in the axe x or y


you can use the code for the Zou-He boundary condition boudaryConditions/zouHeDynamics.h (if I’m not mistaken) as inspiration. Usually the f_i are accessed through the operator “[ ]” of the Cell objects.

Hmm, thanks. I’m not sure if I can just edit them, or should I recalculate anything after that?

BTW, what directions numbers Palabos uses for d2q9 model??
Something like that?

  f4  f3 f2
    \ | /
f5 - f0 - f1
    / | \
  f6  f7 f8

or is the (0,0) direction f8?

I was reading both dynamics and boundary conditions codes and they’re all pretty long. What way would you recommend for joining two lattices (canals) with the formulas mentioned above? It’s all purely local, depending on distribution functions values and two parameters I can calculate beforehand.
I’m looking for the simplest method possible, because I have little time and I’ll have to write very similar code for other obstacles.

Thank you for all your time, without your advices I’d be lost. The library is quite big, hard to get a grasp of.


the numbering in palabos for d2q9 is

f1 f8 f7
\ | /
f2 - f0 - f6
/ |
f3 f4 f5

If you want to write something similar to what you posted before, one thing that you can so is just copy paste the Zou-He BC that is located in the zouHeDynamics.h* and the zouHeDynamicsXd.h* files and replace ZouHe by another name (it should be done quickly). The only method that you have to replace is the collide method of the zouHeVelocityDynamics. And then you should be ok.

There is no collide method implemented in zouHeVelocityDynamics. It’s inherited from CompositeDynamics. Should I override it? Or just edit the completePopulations method?

Also, can I somehow access v (dx/dt) from inside of zouHeVelocityDynamics?

And is it possible to pass some kind of an argument to the class? Otherwise, I’d have to write 2 classes for each obstacle (left and right side). If I can pass an argument, I could write one that’d choose the appropriate function basing on it.

Oups… Sorry you are completely right… I meant CompletePopulations… I think that the best you can do is rename the ZouHeDynamics class to make it do what you want by rewriting the completePopulations method.

The dynamics do not contain global information of this kind. If you want to have access to the units conversion dx or dt in the dynamics you should have a member variable that is assigned to the desired value in the constructor for example.

Another possibility is to use global definitions. An example of the kind of classes you have to write can be found in carreauGlobalDefs.h and carreauGlobalDefs.cpp. Then you can define the members in your main for example by using the following example (here you want to set a global parameter named Nu0 at a value parameters.getNu0())


To access the value stored in Nu0 you then write something like


For the direction of the wall you can use templates (just like in the ZouHeDynamics where you have the direction and orientation template parameters).

I copied ZouHe* files, changed the names, classes and added necessary includes, but to avoid confusion, lets just assume I’m editing ZouHe classes.

How can I pass an argument to the ZouHeVelocityDynamics constructor? In the program, I create OnLatticeBoundaryCondition2D object using createZouHeBoundaryCondition2D function.
createZouHeBoundaryCondition2D is defined in zouHeBoundary2D.hh file as “return new BoundaryConditionInstantiator2D”
that’s where I get lost, because BoundaryConditionInstantiator2D class is defined in boundaryInstantiator2D.h and it doesn’t contain any constructor.
And the only predecessor is OnLatticeBoundaryCondition2D defined in file boundaryCondition2D.h which doesn’t have any constructor either.

Obviously I missed something, cause the program works and constructor body is evaluated. Could you help me find out what?

I’ll look into the global definitions for now, but I’ll need to pass some arguments to the constructor anyway.

I just encountered another problem. To calculate new f values, I need to know the height difference on both sides of the gate. The easiest way would be to store all these values in the VelocityDynamics class, but completePopulations function is constant, so the access is read-only.
Even if I used the global variable, I would need to know the exact coordinates of the cell. Can I access these informations?

Because both sides of the gate are on different lattices, the best way would be to store pointer to the other side of the gate in each of them, get information on the actual cell position and access appropriate cell in the other lattice to read the distribution values. Is it possible?

Or is it not possible and I should use data processor instead?

Ok, I learned that I can’t access coordinates from boundary condition dynamics, so I switched to data processors.

I wrote a simple OneCellIndexed functional that does nothing and am trying to run it with each iteration, but I keep getting segmentation fault. Is there any method that allows me to define an indexed functional that’d be run with each iteration? (something like integrateProcessingFunctional?) Or should I create a new functional object for each iteration?

Another thing I’d like to do is to disable streaming between 2 columns (in this case nx/2 and nx/2+1). Can it be done? Or should I define 2 separate lattices and apply functional to the end of one and beginning of another?

#include "palabos2D.h"
#include "palabos2D.hh"
#include <cstdlib>
#include <iostream>

using namespace plb;
using namespace std;

typedef double T;

#define DESCRIPTOR descriptors::D2Q9Descriptor

template<typename T, template<typename U> class Descriptor>
class Gate : public OneCellIndexedFunctional2D<T,Descriptor> {
    Gate(plint ny_, bool left_)
        : ny(ny_),
    { }
    Gate<T,Descriptor>* clone() const {
        return new Gate<T,Descriptor>(*this);
    virtual void execute(plint iX, plint iY, Cell<T,Descriptor>& cell) const {
//	T h = 0.;
//	for (plint i=0; i<9; i++)
//		h += cell[i];
//	cout << "Height in (" << iX << "," << iY << ") = " << h << endl;
    plint ny;
    bool left;

Gate<T, DESCRIPTOR> *gleft, *gright;

void canalSetup( MultiBlockLattice2D<T, DESCRIPTOR>& water,
                          T rho0,
                          T force )
    int nx = water.getNx();
    int ny = water.getNy();
    // Bounce-back on bottom wall
    defineDynamics(water, Box2D(0,nx-1, 0,0), new BounceBack<T, DESCRIPTOR>(rho0) );
    // Bounce-back on top wall
    defineDynamics(water, Box2D(0,nx-1, ny-1,ny-1), new BounceBack<T, DESCRIPTOR>(rho0) );

    // External force
//    setExternalVector(water, water.getBoundingBox(),DESCRIPTOR<T>::ExternalField::forceBeginsAt, Array<T,2>(0.,-force));


int main(int argc, char *argv[])
    plbInit(&argc, &argv);
    const T omega = 1.0;
    const int nx   = 1200;
    const int ny   = 400;
    T force        = 0.15/(T)ny;
    T rho0 = 0.;
    const int maxIter  = 16000;
    const int saveIter = 100;

    MultiBlockLattice2D<T, DESCRIPTOR> water (nx,ny, new BGKdynamics<T, DESCRIPTOR>(omega) );


    canalSetup(water, rho0, force);
    pcout << "Starting simulation" << endl;
    for (int iT=0; iT<maxIter; ++iT) {
        if (iT%saveIter==0) {
    		ImageWriter<T> imageWriter("leeloo");
			imageWriter.writeScaledGif(createFileName("u", iT, 6), *computeVelocityNorm(water));

        gleft = new Gate<T, DESCRIPTOR>(ny,true);
        gright = new Gate<T, DESCRIPTOR>(ny,false);
	applyIndexed(water, Box2D(nx/2, nx/2, 0, ny-1), gleft);
	applyIndexed(water, Box2D(nx/2+1, nx/2+1, 0, ny-1), gright);;


You never delete your gleft and gright pointers which may be a problem. A you sure it is not related with your segmentation fault?

Ah sorry, I didn’t update the information. I don’t get the segmentation fault anymore. I had it when I tried to use one gate may times, or tried to delete it. I think palabos is taking care of all the pointers that I pass as arguments, so it’s ok.
I already switched from dynamics to data processors though.

I’m trying to write a data processor that’d connect two separate lattices, but with no success. Don’t know how to turn off boundary condition completely, it always gets bounced back.
I decided this topic is way too messy and created another -,2921