# Simple Diffusion [solved, example code included]

Hi,

I’m still working my way into palabos, and to get a feeling how to use diffusion in more complex setups, I’d like to simulate a very simple diffusion setup: One drop of one fluid in the center, around it another fluid (later on I’ll need two fluids of differend viscosity, so I need multi-component which can completely mix - no fixed border between the lattices).

I’d now like to watch the random diffusion in this setup, but I’m not sure how to do that (rayleigh taylor looks far too complex for that simple task). The end-result with two fluids of equal viscosity should be complete mixing of both fluids, but I’m interested in the intermediate dynamics.

Hi,

The Rayleigh-Taylor code is certainly a good template for getting started with a multi-component simulation, even if your problem is different. In particular, you will probably want to re-use the strategy for creating the initial condition (defining the presence of fluid A in some areas and fluid B in some other areas).

I see two major differences between your problem and the simulation for the Rayleigh-Taylor instability. First, you don’t have any gravitation, so you must set the external force to zero. Second, you system is miscible, whereas the two fluids in the Rayleigh-Taylor problem are immiscible. You achieve miscible fluids by decreasing the amplitude G of the interaction potential. I don’t remember the critical value of G where you switch from the immiscible to the miscible regime, but you will easily find it by digging through the literature for the Shan/Chen model.

Hi jlatt,

Thank you for the pointers!

I’ll start off from Rayleigh-Taylor then and adapt it bit by bit (shutting down gravity shouldn’t be the problem, and G is in the initial conditions, so I now know where to start).

Are you interested in seeing the results here?

Hi,
Yes, I am absolutely interested in seeing the results! I hope that by exchanging our results here on the forum or on the public Wiki of LBMethod.org we can build a common know-how on the usage of Palabos and the various LB models.

The simple 3D diffusion works now. I have one more question though: Would I go massively wrong when I just used omega1 and omega2 for setting different viscosities?

Here’s the code I use:

``````
/* This file is part of the Palabos library.
* Copyright (C) 2009 Jonas Latt
* E-mail contact: jonas@lbmethod.org
* Copyright (C) 2010 Arne Babenhauserheide
* Contact: http://draketo.de/contact
* <http://www.lbmethod.org/palabos/>
*
* The library Palabos is free software: you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3 of the
*
* The library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

/** \file
* Simulation of simple diffusion.
* It is adapted from the RayleighTaylor simulation.
**/

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

using namespace plb;
using namespace std;

// Use double-precision arithmetics
typedef double T;
// Use a grid which additionally to the f's stores two variables for
//   the external force term.
#define DESCRIPTOR descriptors::ForcedShanChenD3Q19Descriptor

/// Initial condition: one fluid in center, the other one around it.
/** This functional is going to be used as an argument to the function "applyIndexed",
*  to setup the initial condition. For efficiency reasons, this approach should
*  always be preferred over explicit space loops in end-user codes.
*/
template<typename T, template<typename U> class Descriptor>
class TwoFluidInitializer : public OneCellIndexedFunctional3D<T,Descriptor> {
public:
TwoFluidInitializer(plint nx_, plint ny_, plint nz_, bool topLayer_, T force_)
: nx(nx_),
ny(ny_),
nz(nz_),
topLayer(topLayer_),
force(force_)
{ }
TwoFluidInitializer<T,Descriptor>* clone() const {
return new TwoFluidInitializer<T,Descriptor>(*this);
}
virtual void execute(plint iX, plint iY, plint iZ,
Cell<T,Descriptor>& cell) const {
T densityFluctuations = 1.e-2;
T almostNoFluid       = 1.e-1;
Array<T,2> zeroVelocity (0.,0.);

T rho = (T)1;

plint radiusSquared = pow(nx / 5, 2);

// Reduce the density where the fluid should not be.
T distanceFromCenter = pow(iX - nx/2, 2) + pow(iY - ny/2, 2) + pow(iZ - nz/2, 2);
if ( (topLayer && distanceFromCenter > radiusSquared) || (!topLayer && distanceFromCenter <= radiusSquared) ) {
rho = almostNoFluid;
}

iniCellAtEquilibrium(cell, rho, zeroVelocity);
}
private:
plint nx;
plint ny;
plint nz;
bool topLayer;
T force;
};

void rayleighTaylorSetup( MultiBlockLattice3D<T, DESCRIPTOR>& innerFluid,
MultiBlockLattice3D<T, DESCRIPTOR>& outerFluid,
T rho0, T rho1,
T force )
{
// The setup is: heavy fluid in the center, light around it.
int nx = innerFluid.getNx();
int ny = innerFluid.getNy();
int nz = innerFluid.getNz();

// Initialize center fluid.
applyIndexed(innerFluid, Box3D(0, nx-1, 0, ny-1, 0, nz-1),
new TwoFluidInitializer<T,DESCRIPTOR>(nx, ny, nz, true, force) );
// Initialize outer fluid.
applyIndexed(outerFluid, Box3D(0, nx-1, 0, ny-1, 0, nz-1),
new TwoFluidInitializer<T,DESCRIPTOR>(nx, ny, nz, false, T()) );

outerFluid.initialize();
innerFluid.initialize();
}

void writeGifs(MultiBlockLattice3D<T, DESCRIPTOR>& innerFluid,
MultiBlockLattice3D<T, DESCRIPTOR>& outerFluid, int iT)
{
const plint imSize = 600;
const plint nx = innerFluid.getNx();
const plint ny = innerFluid.getNy();
const plint nz = innerFluid.getNz();
Box3D slice(0, nx-1, 0, ny-1, nz/2, nz/2);

ImageWriter<T> imageWriter("leeloo.map");
imageWriter.writeScaledGif(createFileName("rho_heavy_", iT, 6),
*computeDensity(innerFluid, slice), imSize,imSize);
imageWriter.writeScaledGif(createFileName("rho_light_", iT, 6),
*computeDensity(outerFluid, slice), imSize,imSize);
}

void writeVTK(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
plint iter)
{
double dx = 0.002; // 1/N
double dt = 0.00004; // dx² -> http://www.lbmethod.org/_media/howtos:lbunits.pdf
VtkImageOutput3D<double> vtkOut(createFileName("vtk", iter, 6), dx);
vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", dx/dt);
vtkOut.writeData<float>(*computeDensity(lattice), "density", 1.);
//vtkOut.writeData<float>(*subtract(*computeDensity(lattice),(T)1.), "density", 1.);
vtkOut.writeData<3,float>(*computeVelocity(lattice), "velocity", dx/dt);
}

int main(int argc, char *argv[])
{
plbInit(&argc, &argv);
global::directories().setOutputDir("./tmp/");

const T omega1 = 1.0;
const T omega2 = 1.2;
const int nx   = 50;
const int ny   = 50;
const int nz   = 50;
const T G      = 0.32;
T force        = 0.;//1.e-3;
const int maxIter  = 10000;
const int saveIter = 5;
const int vtkIter = 100;
const int statIter = 50;

// Use regularized BGK dynamics to improve numerical stability (but note that
//   BGK dynamics works well too).
MultiBlockLattice3D<T, DESCRIPTOR> innerFluid (
nx,ny,nz, new ExternalMomentRegularizedBGKdynamics<T, DESCRIPTOR>(omega1) );
MultiBlockLattice3D<T, DESCRIPTOR> outerFluid (
nx,ny,nz, new ExternalMomentRegularizedBGKdynamics<T, DESCRIPTOR>(omega2) );

// Make all directions periodic.
innerFluid.periodicity().toggle(0, true);
innerFluid.periodicity().toggle(1, true);
innerFluid.periodicity().toggle(2, true);
outerFluid.periodicity().toggle(0, true);
outerFluid.periodicity().toggle(1, true);
outerFluid.periodicity().toggle(2, true);

T rho1 = 1.; // Fictitious density experienced by the partner fluid on a Bounce-Back node.
T rho0 = 0.; // Fictitious density experienced by the partner fluid on a Bounce-Back node.

// Store a pointer to all lattices (two in the present application) in a vector to
//   create the Shan/Chen coupling therm. The heavy fluid being at the first place
//   in the vector, the coupling term is going to be executed at the end of the call
//   to collideAndStream() or stream() for the heavy fluid.
vector<MultiBlockLattice3D<T, DESCRIPTOR>* > blockLattices;
blockLattices.push_back(&innerFluid);
blockLattices.push_back(&outerFluid);

// The argument "constOmegaValues" to the Shan/Chen processor is optional,
//   and is used for efficiency reasons only. It tells the data processor
//   that the relaxation times are constant, and that their inverse must be
//   computed only once.
std::vector<T> constOmegaValues;
constOmegaValues.push_back(omega1);
constOmegaValues.push_back(omega2);
plint processorLevel = 1;
integrateProcessingFunctional (
new ShanChenMultiComponentProcessor3D<T,DESCRIPTOR>(G),//,constOmegaValues),
Box3D(0,nx-1,1,ny-2,0,nz-1),
blockLattices, processorLevel );

rayleighTaylorSetup(innerFluid, outerFluid, rho0, rho1, force);

pcout << "Starting simulation" << endl;
// Main loop over time iterations.
for (int iT=0; iT<maxIter; ++iT) {
if (iT%saveIter==0) {
writeGifs(innerFluid, outerFluid, iT);
}
if (iT%vtkIter==0) {
writeVTK(innerFluid, iT);
}

// Time iteration for the light fluid.
outerFluid.collideAndStream();
// Time iteration for the heavy fluid must come after the light fluid,
//   because the coupling is executed here. You should understand this as follows.
//   The effect of the coupling is to compute the interaction force between
//   species, and to precompute density and momentum for each species. This must
//   be executed *before* collide-and-streaming the fluids, because the collision
//   step needs to access all these values. In the present case, it is done after
//   both collide-and-stream step, which means, before the collide-and-stream of
//   the next iteration (it's the same if you are before or after; the important
//   point is not to be between the two collide-and-streams of the light and heavy
//   fluid. As for the initial condition, the coupling is initially performed once
//   during the function call to innerFluid.initialize().
innerFluid.collideAndStream();

if (iT%statIter==0) {
pcout << "Average energy fluid one = "
<< getStoredAverageEnergy<T>(innerFluid);
pcout << ", average energy fluid two = "
<< getStoredAverageEnergy<T>(outerFluid) << endl;
}
}
}

``````

Hi,

Nice to hear that it works. As far as i know, multi-component fluids behave the same as single-component ones with respect to viscosity. So yes, omega1 and omega2 are the parameters you need to change in order to adjust the viscosities. However, although I don’t have much experience in the field, I think that you quickly run into stability issues when the two viscosities are too different. The simplest thing to do is to try.

Hi,

It indeed gets unstable at some point, but it seems to work when the differences in viscosity aren’t too big. Many thanks!

Besides: I changed the subject to include the information, that there’s example code to find here.

Hi ArneBab,

For some weird reason, I can’t compile your sample code. I get the following error message:

``````
python ../../../scons/scons.py -j 6 -f ../../../SConstruct palabosRoot=../../.. projectFiles="Sim2_g.cpp" optimize=true debug=false profile=false MPIparallel=true SMPparallel=false usePOSIX=true serialCXX=g++ parallelCXX=mpicxx compileFlags="-Wall -Wnon-virtual-dtor -DPLB_MAC_OS_X" linkFlags="" optimFlags="-O3" debugFlags="-g" profileFlags="-pg" libraryPaths="" includePaths="" libraries=""
scons: Building targets ...
mpicxx -o Sim2_g.o -c -Wall -Wnon-virtual-dtor -DPLB_MAC_OS_X -O3 -DPLB_MPI_PARALLEL -DPLB_USE_POSIX -I/home/Tiago/palabos-v1.5r1/src -I/home/Tiago/palabos-v1.5r1/externalLibraries Sim2_g.cpp
Sim2_g.cpp: In instantiation of ‘void TwoFluidInitializer<T, Descriptor>::execute(plb::plint, plb::plint, plb::plint, plb::Cell<T, Descriptor>&) const [with T = double; Descriptor = plb::descriptors::ForcedShanChenD3Q19Descriptor; plb::plint = long int]’:
Sim2_g.cpp:226:1:   required from here
Sim2_g.cpp:76:53: error: no matching function for call to ‘iniCellAtEquilibrium(plb::Cell<double, plb::descriptors::ForcedShanChenD3Q19Descriptor>&, double&, plb::Array<double, 2ul>&)’
iniCellAtEquilibrium(cell, rho, zeroVelocity);
^
Sim2_g.cpp:76:53: note: candidates are:
from /home/Tiago/palabos-v1.5r1/src/palabos3D.hh:27,
from Sim2_g.cpp:29:
/home/Tiago/palabos-v1.5r1/src/core/cell.hh:124:6: note: void plb::iniCellAtEquilibrium(plb::Cell<T, Descriptor>&, T, const plb::Array<T, Descriptor<T>::d>&) [with T = double; Descriptor = plb::descriptors::ForcedShanChenD3Q19Descriptor]
void iniCellAtEquilibrium(Cell<T,Descriptor>& cell, T density, Array<T,Descriptor<T>::d> const& velocity) {
^
/home/Tiago/palabos-v1.5r1/src/core/cell.hh:124:6: note:   no known conversion for argument 3 from ‘plb::Array<double, 2ul>’ to ‘const plb::Array<double, 3ul>&’
/home/Tiago/palabos-v1.5r1/src/core/cell.hh:135:6: note: template<class T, template<class U> class Descriptor> void plb::iniCellAtEquilibrium(plb::Cell<T, Descriptor>&, T, const plb::Array<T, Descriptor<T>::d>&, T)
void iniCellAtEquilibrium(Cell<T,Descriptor>& cell, T density, Array<T,Descriptor<T>::d> const& velocity, T temperature) {
^
/home/Tiago/palabos-v1.5r1/src/core/cell.hh:135:6: note:   template argument deduction/substitution failed:
Sim2_g.cpp:76:53: note:   candidate expects 4 arguments, 3 provided
iniCellAtEquilibrium(cell, rho, zeroVelocity);
^
Sim2_g.cpp:62:11: warning: unused variable ‘densityFluctuations’ [-Wunused-variable]
T densityFluctuations = 1.e-2;
^
scons: *** [Sim2_g.o] Error 1
scons: building terminated because of errors.
Makefile:87: recipe for target 'compile' failed
make: *** [compile] Error 2

``````

Does someone know what might be happening?

@asauser

I found the same problem, under Debian Linux 8 64-bit. However, since I am interested in modelling a 2D diffusion problem, I adapted ArneBab’s code to two dimensions and it compiled flawlessly. First of all, hat tip to ArneBab for sharing the code! With the due permission, I post the 2D variant of it:

[code=“cpp”]
/* This file is part of the Palabos library.

• E-mail contact: jonas@lbmethod.org
• Contact: [draketo.de]
• [www.lbmethod.org]
• The library Palabos is free software: you can redistribute it and/or
• modify it under the terms of the GNU General Public License as
• published by the Free Software Foundation, either version 3 of the
• The library is distributed in the hope that it will be useful,
• but WITHOUT ANY WARRANTY; without even the implied warranty of
• MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
• GNU General Public License for more details.
• You should have received a copy of the GNU General Public License
• along with this program. If not, see [www.gnu.org].
*/

/** \file

• Simulation of simple diffusion.
• It is adapted from the RayleighTaylor simulation.
**/

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

using namespace plb;
using namespace std;

// Use double-precision arithmetics
typedef double T;
// Use a grid which additionally to the f’s stores two variables for
// the external force term.
#define DESCRIPTOR descriptors::ForcedShanChenD2Q9Descriptor

/// Initial condition: one fluid in center, the other one around it.
/** This functional is going to be used as an argument to the function “applyIndexed”,

• to setup the initial condition. For efficiency reasons, this approach should

• always be preferred over explicit space loops in end-user codes.
/
template<typename T, template class Descriptor>
class TwoFluidInitializer : public OneCellIndexedFunctional2D<T,Descriptor> {
public:
TwoFluidInitializer(plint nx_, plint ny_, bool topLayer_, T force_)
: nx(nx_),
ny(ny_),
topLayer(topLayer_),
force(force_)
{ }
TwoFluidInitializer<T,Descriptor>
clone() const {
return new TwoFluidInitializer<T,Descriptor>(*this);
}
virtual void execute(plint iX, plint iY,
Cell<T,Descriptor>& cell) const {
T densityFluctuations = 1.e-2;
T almostNoFluid = 1.e-1;
Array<T,2> zeroVelocity (0.,0.);

``````T rho = (T)1;
``````

plint radiusSquared = pow(nx / 5, 2);

``````// Reduce the density where the fluid should not be.
``````

T distanceFromCenter = pow(iX - nx/2, 2) + pow(iY - ny/2, 2);
if ( (topLayer && distanceFromCenter > radiusSquared) || (!topLayer && distanceFromCenter <= radiusSquared) ) {
rho = almostNoFluid;
}

``````iniCellAtEquilibrium(cell, rho, zeroVelocity);
``````

}
private:
plint nx;
plint ny;
bool topLayer;
T force;
};

void rayleighTaylorSetup( MultiBlockLattice2D<T, DESCRIPTOR>& innerFluid,
MultiBlockLattice2D<T, DESCRIPTOR>& outerFluid,
T rho0, T rho1,
T force )
{
// The setup is: heavy fluid in the center, light around it.
int nx = innerFluid.getNx();
int ny = innerFluid.getNy();

``````// Initialize center fluid.
applyIndexed(innerFluid, Box2D(0, nx-1, 0, ny-1),
new TwoFluidInitializer<T,DESCRIPTOR>(nx, ny, true, force) );
// Initialize outer fluid.
applyIndexed(outerFluid, Box2D(0, nx-1, 0, ny-1),
new TwoFluidInitializer<T,DESCRIPTOR>(nx, ny, false, T()) );

outerFluid.initialize();
innerFluid.initialize();
``````

}

void writeGifs(MultiBlockLattice2D<T, DESCRIPTOR>& innerFluid,
MultiBlockLattice2D<T, DESCRIPTOR>& outerFluid, int iT)
{
const plint imSize = 600;
const plint nx = innerFluid.getNx();
const plint ny = innerFluid.getNy();
Box2D slice(0, nx-1, 0, ny-1);

``````ImageWriter<T> imageWriter("leeloo.map");
imageWriter.writeScaledGif(createFileName("rho_heavy_", iT, 6),
*computeDensity(innerFluid, slice), imSize,imSize);
imageWriter.writeScaledGif(createFileName("rho_light_", iT, 6),
*computeDensity(outerFluid, slice), imSize,imSize);
``````

}

int main(int argc, char *argv[])
{
plbInit(&argc, &argv);
global::directories().setOutputDir("./tmp/");

``````const T omega1 = 1.0;
const T omega2 = 1.2;
const int nx   = 50;
const int ny   = 50;
const T G      = 0.32;
T force        = 0.;//1.e-3;
const int maxIter  = 5000;
const int saveIter = 5;
//const int vtkIter = 100;
const int statIter = 50;

// Use regularized BGK dynamics to improve numerical stability (but note that
//   BGK dynamics works well too).
MultiBlockLattice2D<T, DESCRIPTOR> innerFluid (
nx,ny, new ExternalMomentRegularizedBGKdynamics<T, DESCRIPTOR>(omega1) );
MultiBlockLattice2D<T, DESCRIPTOR> outerFluid (
nx,ny, new ExternalMomentRegularizedBGKdynamics<T, DESCRIPTOR>(omega2) );

// Make all directions periodic.
innerFluid.periodicity().toggle(0, true);
innerFluid.periodicity().toggle(1, true);
outerFluid.periodicity().toggle(0, true);
outerFluid.periodicity().toggle(1, true);

T rho1 = 1.; // Fictitious density experienced by the partner fluid on a Bounce-Back node.
T rho0 = 0.; // Fictitious density experienced by the partner fluid on a Bounce-Back node.

// Store a pointer to all lattices (two in the present application) in a vector to
//   create the Shan/Chen coupling therm. The heavy fluid being at the first place
//   in the vector, the coupling term is going to be executed at the end of the call
//   to collideAndStream() or stream() for the heavy fluid.
vector<MultiBlockLattice2D<T, DESCRIPTOR>* > blockLattices;
blockLattices.push_back(&innerFluid);
blockLattices.push_back(&outerFluid);

// The argument "constOmegaValues" to the Shan/Chen processor is optional,
//   and is used for efficiency reasons only. It tells the data processor
//   that the relaxation times are constant, and that their inverse must be
//   computed only once.
std::vector<T> constOmegaValues;
constOmegaValues.push_back(omega1);
constOmegaValues.push_back(omega2);
plint processorLevel = 1;
integrateProcessingFunctional (
new ShanChenMultiComponentProcessor2D<T,DESCRIPTOR>(G),//,constOmegaValues),
Box2D(0,nx-1,1,ny-2),
blockLattices, processorLevel );

rayleighTaylorSetup(innerFluid, outerFluid, rho0, rho1, force);

pcout << "Starting simulation" << endl;
// Main loop over time iterations.
for (int iT=0; iT<maxIter; ++iT) {
if (iT%saveIter==0) {
writeGifs(innerFluid, outerFluid, iT);
}

// Time iteration for the light fluid.
outerFluid.collideAndStream();
// Time iteration for the heavy fluid must come after the light fluid,
//   because the coupling is executed here. You should understand this as follows.
//   The effect of the coupling is to compute the interaction force between
//   species, and to precompute density and momentum for each species. This must
//   be executed *before* collide-and-streaming the fluids, because the collision
//   step needs to access all these values. In the present case, it is done after
//   both collide-and-stream step, which means, before the collide-and-stream of
//   the next iteration (it's the same if you are before or after; the important
//   point is not to be between the two collide-and-streams of the light and heavy
//   fluid. As for the initial condition, the coupling is initially performed once
//   during the function call to innerFluid.initialize().
innerFluid.collideAndStream();

if (iT%statIter==0) {
pcout << "Average energy fluid one = "
<< getStoredAverageEnergy<T>(innerFluid);
pcout << ", average energy fluid two = "
<< getStoredAverageEnergy<T>(outerFluid) << endl;
}
}
``````

}

Hi

How do you find out the G and force value? Is there any formula or relationship? Also, how do we find delta t in the current code?

I am trying to write a piece of code which can apply a periodic external force on the liquids. Now, I have written the function to do it but the value is dependent on time but in the code given in the showCases (MultiComponent2d), there is no mention of delta t. I can understand it is because of some lattice unit conversion but can someone explain it to me?

Vishwesh

The issue is that the variable `zeroVelocity` is defined as `Array<T,2>` and is not compatible with the function definition of `iniCellAtEquilibrium` as given in `src/core/cell.h`

``````void iniCellAtEquilibrium(Cell<T,Descriptor>& cell, T density, Array<T,Descriptor<T>::d> const& velocity);
``````

The solution is to replace that line with

``````Array<T,3> zeroVelocity (0.,0.,0.);
``````