I hope this code can help people to get into working with palabos.
Please feel free to integrate it.
[code=cpp]
/** A Simulation of blood in a vessel.
*
- Copyright © 2010 Arne Babenhauserheide
- E-mail contact: arne_bab@web.de
- This file and the library Palabos are 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
- License, or (at your option) any later version.
- 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/.
*/
/** Main author: Arne Babenhauserheide
*/
/** Notes on the code
I comment this file extensively to make it easy to follow and understand. My expectition on myself is that
- You should be able to understand what it does from the comments alone, without having to look at the code.
- The comments should explain the reasons why things are done how they are done, and they should describe those internals of Palabos, which are useful for writing other simulations.
It uses the aneurysm STL (from readStl), because I don’t know the status of the file I use for my real simulation.
This file is massively inspired by readStl/sparseDomainFromSTL.cpp which I use as base for my simulation.
*/
//// Includes ////
// The first step is including the Palabos headers for a 3D simulation.
// They are found via scons (the custom makefile).
#include “palabos3D.h” // this file points to includes inside palabos.
#ifndef PLB_PRECOMPILED // if we don’t use a precompiled version
#include “palabos3D.hh” // we also need the header templates.
// They contain some of the implementations; if I understand it correctly the
// reason for that is optimization and template usage.
#endif
//// Constants which likely won’t change. ////
// variables are defined in the main function
// so they can optionally be passed to the function.
/* The next step is telling Palabos which kind of lattice we want to use.
For a 3D simulation the simplest full featured one has 19 degrees of freedom:
D3Q19
Using the preprocessor here allows us to easily change the lattice later on.
*/
#define DESCRIPTOR descriptors::D3Q19Descriptor
// also we want to use double for numbers
typedef double T;
// additionally we define the density the STL importer (voxellizer) should use for walls
// (bounce back walls) and the density for the inlet “boxes” (truncation).
const int voxellizer_bbDensity = 1;
const int truncation_bbDensity = 1;
// Now we enter the palabos and std namespaces, so we don’t have to type plb::
// or std:: all the time. => ensure that there never are clashes between palabor
// classes and the std.
using namespace std;
using namespace plb;
//// Helper functions ////
// Now we define the functions for outputting data.
// These aren’t part of the simulation itself, so I keep them seperate.
// Since they depend on functions from palabos, I can’t easily move them into
// another file.
/** write a gif.
Remember that DESCRIPTOR is the lattice type,
And T is the type of number we want to use.
incomprFlowParam are used to convert lattice units into real units.
The iter tells us at which iteration we are, so we can choose the filename accordingly.
*/
void writeGif(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters, plint iter)
{
// First set the desired image size
const plint imSize = 600;
// Then get the distance between lattice nodes.
const plint nx = parameters.getNx();
const plint nz = parameters.getNz();
// Set the y position of the slice we want to see.
plint yPos = parameters.getResolution()*0.25;
// And get a slice (3D rectangle) whose nodes we want to see
Box3D slice(0, nx-1, yPos, yPos, 0, nz-1);
// Get the image writer with the leeloo colormap.
// For alternative styles see
// src/io/colormaps.cpp - generateMap(std::string mapName)
ImageWriter<T> imageWriter("leeloo");
// finally write the gif (needs imagemagick, as it first writes a .pnm
// and then uses ‘convert’ to turn it into a gif
// -> http://www.imagemagick.org/ ).
imageWriter.writeScaledGif( createFileName("u", iter, 6),
*computeVelocityNorm (lattice, slice),
imSize, imSize );
}
/** Check if a given cell of the lattice is NOT a bounce back node. */
bool isNotBounceBack(Cell<T,DESCRIPTOR> const& cell)
{
return !( typeid(cell.getDynamics()) == typeid(BounceBack<T,DESCRIPTOR> const&) );
}
//// The Simulation ////
// Now we get to the really interesting part: the simulation itself.
/** Add borders, in- and outlets to a lattice (after we added the parsed STL)
Since the template setup looks a bit daunting:
T is the kind of number to use for calculation and Descriptor the type of lattice.
<typename U> class Descriptor allows using different lattices (see DESCRIPTOR).
This generates the code for all types we might need.
lattice is taken as reference, to enable us to directly edit it.
*/
template <typename T, template class Descriptor>
void addBordersToLattice( MultiBlockLattice3D<T,Descriptor>& lattice,
plint N)
{
// First we need an empty boundary condition
// see http://www.lbmethod.org/palabos/documentation.userguide/boundary-conditions.html#the-class-onlatticeboundaryconditionxd
OnLatticeBoundaryCondition3D<T, Descriptor>* boundaryCondition
= createLocalBoundaryCondition3D<T,Descriptor>();
/* Then prepare the pressure gradients at the different inlets.
These specify how strong the inflow will be. */
const T gradP1 = 1.0e-3;
// Now we add sheets for the inlets and outlets, including long boxes before them
// (why? - maybe so the flow is restricted to the STL or such? They get BounceBack nodes assigned)
Box3D inlet1( (plint)(0.1*N), (plint)(0.3*N), // x0, x1
(plint)(0.08*N), (plint)(0.3*N), // y0, y1
(plint)(0.05*N), (plint)(0.05*N) ); // z0, z1
// The long box before it has the same x and y dimensions,
// but goes through the whole box in the z direction.
Box3D beforeInlet1(inlet1.x0, inlet1.x1,
inlet1.y0, inlet1.y1,
0, inlet1.z1-1);
// … additional in- and outlets can de added.
// Now we add the boundary conditions on the inlets and their boxes to the empty boundary condition,
// so they become real boundaries (until here, they are only boxes without meaning)
// A Boundary2N is a boundary facing into the negative direction of the z axis
// (the axi are numbered beginning with 0).
boundaryCondition->addPressureBoundary2N(inlet1, lattice);
// Also we add the fixed density (which gets translated into pressure and so into inflow).
// It is higher than the fluid density (1.0) by the pressure gradient we want at the inflow.
// Remember again, that T is just the format used for numbers and Descriptor is the lattice type we use.
setBoundaryDensity(lattice, inlet1, (T)1.+gradP1);
// Then we add bounce back conditions at the box before the inlet (which gives us walls)
// Dynamics are the interaction of lattice elements. In this they define bounce back behaviour.
defineDynamics(lattice, beforeInlet1, new BounceBack<T,Descriptor>(truncation_bbDensity));
// additional outlets
Box3D outlet1( (plint)(0.2*N), (plint)(0.36*N),
(plint)(0.58*N), (plint)(0.74*N),
(plint)(0.5*N), (plint)(0.5*N) );
Box3D beforeOutlet1(outlet1.x0, outlet1.x1, outlet1.y0, outlet1.y1,
outlet1.z0-(plint)(0.1*N), outlet1.z1-1);
Box3D outlet2( (plint)(0.24*N), (plint)(0.24*N),
(plint)(0.47*N), (plint)(0.64*N),
(plint)(0.82*N), (plint)(0.96*N) );
Box3D beforeOutlet2((plint)(0.170*N), outlet2.x1-1,
outlet2.y0, outlet2.y1, outlet2.z0, outlet2.z1);
boundaryCondition->addPressureBoundary2N(outlet1, lattice);
setBoundaryDensity(lattice,outlet1,(T)1.);
defineDynamics(lattice, beforeOutlet1, new BounceBack<T,Descriptor>(truncation_bbDensity));
boundaryCondition->addPressureBoundary0N(outlet2, lattice);
setBoundaryDensity(lattice,outlet2,(T)1.);
defineDynamics(lattice, beforeOutlet2, new BounceBack<T,Descriptor>(truncation_bbDensity));
// finally we initalize the whole lattice at equlibrium (not counting the inlets)
// lattice.getBoundingBox() gives a Box3D around the whole lattice.
// Array<T,3>(x, y, z) is a flexible size 3D array filled with the given numbers.
initializeAtEquilibrium(lattice, lattice.getBoundingBox(), (T)1., Array<T,3>(0.,0.,0.));
// With this the lattice we got the reference to is prepared. We can now use it for simulation.
}
// Now we get to the controlling function.
/** Prepare and run the simulation
This function is adapted from the main function, so we can include and call it from the outside.
*/
int runSimulation(int argc, char* argv[], // commandline arguments
T uMax = (T) 1e-2, T Re = (T) 100, // flow parameters
plint N = 127, float lx = 1., float ly = 1., float lz = 1., // Simulation properties
plint maxIter = 6001, // How long should the simulation run?
plint imageIter = 30, plint statusIter = 10, // output paramaters
int limitDepth = 3, string stlFile = “aneurysm.stl”// STL import parameters
) {
// first we initialize Palabos, to process the commandline arguments
plbInit(&argc, &argv);
// Now we define the output directory
global::directories().setOutputDir("./tmp/");
// And process the function arguments.
// the first are the basic parameters for the simulation
IncomprFlowParam<T> parameters(uMax, Re, N, lx, ly, lz);
// Now we read in the STL file and turn it into a lattice
TreeReader3D<int>* treeReader =
createUniformTreeReader3D(stlFile, parameters.getNx(), limitDepth);
MultiBlockLattice3D<T,DESCRIPTOR>* lattice = generateSparseMultiBlock (
*treeReader,
new BounceBack<T,DESCRIPTOR>(voxellizer_bbDensity));
// now we get the fluid parts
bool fluidFlag = true;
dynamicsFromOctree( *lattice, *treeReader,
new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()), fluidFlag);
// add walls and inlets
addBordersToLattice(*lattice, parameters.getNx());
// and define the walls as bounce back
bool wallFlag = false;
dynamicsFromOctree( *lattice, *treeReader,
new BounceBack<T,DESCRIPTOR>(voxellizer_bbDensity), wallFlag);
// finally we cleanup no longer needed objects.
delete treeReader;
// now we give some info about the simulation before we start
pcout << "Structure of the multi-block, with an octree decomposed to depth "
<< limitDepth << ":" << endl;
pcout << getMultiBlockInfo(*lattice);
// isNotBounceBack is the function to decide which nodes to count.
pcout << "Number of fluid nodes=" << count(*lattice, isNotBounceBack) << endl;
// so now it's time to actually start calculating
// first we start the timer, so we can gather statistics how much time we need per step.
global::timer("main").start();
// and enter the loop over all iterations we want to do
for (int iT=0; iT<maxIter; ++iT) {
// The first part of the loop is statistics and image output
if (iT%statusIter==0) {
pcout << "Iteration " << iT << endl;
}
if (iT%imageIter==0) {
global::timer("main").stop();
pcout << "Writing GIF at time step " << iT << endl;
writeGif(*lattice, parameters, iT);
global::timer("main").start();
}
// Then we do the real calculation: one Lattice-Boltzmann iteration step
lattice->collideAndStream();
// Finally we get stats on the efficiency and exactness of the simulation.
if (iT%statusIter==0) {
pcout << "Rms velocity="
<< setprecision(10) << sqrt(2.*getStoredAverageEnergy<T>(*lattice))
<< "; av rho="
<< setprecision(10) << getStoredAverageDensity<T>(*lattice) << endl;
pcout << "Elapsed time: " << global::timer("main").getTime() / iT
<< " seconds per time step" << endl;
}
}
// After the loop we don't need the lattice anymore.
delete lattice;
}
/** Call the simulation run */
int main(int argc, char* argv[]) {
// The simplest way is just using the default simulation.
runSimulation(argc, argv);
}
