Hello,
we have problems with running 3D simulations. For 2D case everything works OK, for 3D we get NaNs in some cells.
The 2D and 3D code looks almost similar, could someone look if we do some obvious errors please ?
Below there is the code for 3D case.
Comments only in Polish (sorry) but the code is quite clean.
file lbm3D_png.cpp:
[code=“cpp”]
#include “palabos3D.h”
#include “palabos3D.hh”
#include
#include
#include
#include <png++/png.hpp>
#include “command_options.hpp”
#include “tools.hpp”
using namespace plb;
using namespace std;
typedef double T;
#define DESCRIPTOR plb::descriptors::D3Q27Descriptor
/********************************************************************************/
/*
Klasa poniżej służy do tego, żeby na podstawie obrazka powiedzieć, które komorki
sceny zawierają płyn, a które są komórkami otoczenia.
Operator () zwraca true, jeśli komorka o zadanych współrzędnych NIE JEST komórką
kanału z płynem.
Dla przypadków trójwymiarowych kanał jest z góry i z dołu obudowany warstwami
komórek nie należących do kanału (takie “pokrywy”).
Klasa jest potrzebna dla funkcji defineDynamics()
Zerżnięta żywcem z palabosowego tutoriala 2
http://www.palabos.org/documentation/tutorial/multi-block.html
*/
template
class BounceBackNodes : public DomainFunctional3D
{
public:
BounceBackNodes(const png::image< png::rgb_pixel > & image, plint ZSIZE)
: _image(image), _zsize(ZSIZE)
{} ;
/// Return true for all cells outside the channel, on which bounce-back
/// dynamics must be instantiated.
virtual bool operator() (plint iX, plint iY, plint iZ) const
{
if ( 0 == iZ || (_zsize-1) == iZ ) return true ; // "pokrywy" z góry i z dołu
png::rgb_pixel pixel = _image.get_pixel(iX, _image.get_height() - iY - 1) ;
if ( (0 == pixel.red) && (0 == pixel.green) && (0 == pixel.blue) ) // black pixel
return false ;
else return true ;
} ;
virtual BounceBackNodes<T>* clone() const
{
return new BounceBackNodes<T>(*this);
} ;
private:
const png::image< png::rgb_pixel > & _image ;
const plint _zsize ;
};
/********************************************************************************/
/*
Klasa jest odwrotnością klasy BounceBackNodes - operator () zwraca true,
jeśli komórka o zadanych współrzędnych NALEŻY do kanału z płynem.
Klasa jest potrzebna do automatycznego podzielenia sceny na kwadratowe poddomeny.
Też zerżnąłem z palabosowego tutoriala 2.
*/
// TODO: Może ta klasa powinna dziedziczyć z BounceBackNodes ?
template
class FluidNodes
{
public:
FluidNodes(const png::image< png::rgb_pixel > & image, plint ZSIZE)
: _image(image), _zsize(ZSIZE)
{} ;
bool operator() (plint iX, plint iY, plint iZ) const
{
return ! BounceBackNodes<T>(_image, _zsize)(iX, iY, iZ) ;
} ;
private:
const png::image< png::rgb_pixel > & _image ;
const plint _zsize ;
};
/********************************************************************************/
/*
Funkcja ustawiająca początkowe parametry symulacji.
Dokładny odpowiednik funkcji o tej samej nazwie z przykladu poiseulle.
Jedyne różnice dotyczą położenia wlotów/wylotów i ustawianych wartości.
*/
void channelSetup (
MultiBlockLattice3D<T,DESCRIPTOR> &lattice,
const png::image< png::rgb_pixel > & image ,
IncomprFlowParam const& parameters,
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition,
double speed,
double rho )
{
Box3D main_box = lattice.getBoundingBox() ;
plint x0 = main_box.x0 ;
plint x1 = main_box.x1 ;
plint y0 = main_box.y0 ;
plint y1 = main_box.y1 ;
plint z0 = main_box.z0 ;
plint z1 = main_box.z1 ;
pcout << "main_box: x0 = " << x0 << ", x1 = " << x1 ;
pcout << ", y0 = " << y0 << ", y1 = " << y1 ;
pcout << ", z0 = " << z0 << ", z1 = " << z1 ;
pcout << "\n" ;
// To są cztery ściany odpowiadające bokom obrazka. Ściany rozciągają się wzdłuż osi Z, czyli obrazek jest położony
// na płaszczyźnie XY i wyciągnięty w kierunku osi Z. Poniższe ściany odpowiadają liniom o szerokości jednego piksela
// wzdłuż wszystkich czterech boków obrazka.
Box3D left_wall (x0, x0, y0, y1, z0, z1) ;
Box3D right_wall (x1, x1, y0, y1, z0, z1) ;
Box3D top_wall (x0, x1, y1, y1, z0, z1) ;
Box3D bottom_wall (x0, x1, y0, y0, z0, z1) ;
/*
mail od RS z poniedziałek, marzec 11, 2013 1:52
Jeśli chodzi o warunki brzegowe to dookoła kanału: bounce-back
na wlocie (wszystkie wloty kanałów po lewej stronie domeny) “inlet”: warunek stałej prędkości dirichleta (przyjmowany jako domyślny typ dla setVelocity…)
na wylocie (wszystkie wyloty kanałów na górze, na prawo i na dole): warunek dla prędkości neumana “outflow”.
W załączonym pliku to co zrobiliśmy z Piotrkiem - podobne ustawienia dla pojedynczego kanału 3D. Parametry do obliczeń na początku mogą być takie jak w załączonym przykładzie (może się nie wysypie 
).
/
/
Poniżej ustawiam TYPY warunków granicznych.
Ten w połączeniu z powyższymi definicjami left_wall itd. kod odpowiada liniom 128-155 z pliku pouseuille.cpp
Prościej nie umiem, bo będzie jeszcze gorzej.
*/
// Wloty tylko na lewej ścianie
boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, left_wall ) ;
// Reszta to wyloty
boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, right_wall , boundary::outflow ) ;
boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, top_wall , boundary::outflow ) ;
boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, bottom_wall, boundary::outflow ) ;
/*
To poniżej odpowiada liniom 157-171 z pliku pouseuille.cpp.
Różnica jest taka, że ja założyłem, że wartości prędkości są stałe, a w pouseuille.cpp wartość prędkości zależy
od odległości od ścianki kanału.
Poza tym w pouseuille.cpp oni automagicznie liczą wartości prędkości i gestości na podstawie obiektu typu
IncomprFlowParam, a mi sie wydawało, że łatwiej będzie te prędkość i gęstość policzyć ręcznie, bo mniej
parametrów wtedy trzeba znać.
*/
// zeroVelocity jest argumentem initializeAtEquilibrium(), czyli w jakiś sposób dotyczy chyba komórek wewnątrz
// kanału (bez brzegowych wlotów/wylotów)
Array<T,3> zeroVelocity ( 0.0 , 0.0, 0.0 ); // w nawiasach są kolejne współrzędne prędkości po x,y,z
// constantVelocity dotyczy wartości prędkości na wlotach z lewej strony obrazka
Array<T,3> constantVelocity( speed, 0.0, 0.0 ); // w nawiasach są kolejne współrzędne prędkości po x,y,z
// To poniżej ustawia na wszystkich wylotach prędkość równą 0, a na wlotach równą constantVelocity (przynajmniej
// tak mi się wydaje)
setBoundaryVelocity( lattice, lattice.getBoundingBox(), zeroVelocity ) ;
setBoundaryVelocity( lattice, left_wall, constantVelocity ) ;
T constantDensity = (T)rho;
// W pouseuille.cpp zamiast jawnego podawania gęstości i prędkości oni uzywają
// PoiseuilleDensityAndZeroVelocity<T>(IncomprFlowParam<T>), który te parametry jakoś wewnętrznie liczy.
initializeAtEquilibrium ( lattice, lattice.getBoundingBox(), constantDensity, zeroVelocity );
// Poniższej linii nie ma w pouseuille.cpp, bo oni całą scenę wypełnili kanałem i nie mają komórek typu bounceback
defineDynamics(lattice, lattice.getBoundingBox(),
new BounceBackNodes<T>(image, lattice.getNz()), // to mówi, które komórki są odbijającymi
new BounceBack<T,DESCRIPTOR>); // to definiuje ich typ
lattice.initialize(); // to jest linia 174 z poiseulle.cpp
}
/********************************************************************************/
/*
Funkcja zapisująca computeVelocityNorm() do obrazka gif.
Zerżnięta żywcem skądś…
*/
void writeGifs(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
{
const plint imSize = 1000; // jakiś współczynnik skalowania
const plint nx = lattice.getNx();
const plint ny = lattice.getNy();
const plint nz = lattice.getNz();
Box3D slice(0, nx-1, 0, ny-1, nz/2, nz/2);
ImageWriter imageWriter(“leeloo”);
imageWriter.writeScaledGif(createFileName(“u”, iter, 6),
*computeVelocityNorm(lattice, slice),
imSize, imSize );
imageWriter.writeScaledGif(createFileName(“rho”, iter, 6),
*computeDensity(lattice, slice),
imSize, imSize );
}
/********************************************************************************/
/*
Funkcja zapisująca do pliku vtk parametry wybrane z linii poleceń.
Też skądś zerżnąłem.
*/
void writeVTK(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter, const UserOptions & opts)
{
VtkImageOutput3D vtkOut(createFileName(“vtk”, iter, 6));
if ( opts.save_vtk_velocity() )
vtkOut.writeData<3,float>(*computeVelocity (lattice), "velocity" ) ;
if ( opts.save_vtk_density() )
vtkOut.writeData< float>(*computeDensity (lattice), "density" ) ;
if ( opts.save_vtk_rho_bar() )
vtkOut.writeData< float>(*computeRhoBar (lattice), "rho_bar" ) ;
if ( opts.save_vtk_kinetic_energy() )
vtkOut.writeData< float>(*computeKineticEnergy (lattice), "kinetic_energy" ) ;
if ( opts.save_vtk_velocity_norm() )
vtkOut.writeData< float>(*computeVelocityNorm (lattice), "velocity_norm" ) ;
if ( opts.save_vtk_strain_rate_from_stress() )
vtkOut.writeData<6,float>(*computeStrainRateFromStress(lattice), "strain_rate_from_stress" ) ;
// vtkOut.writeData< float>(computeDeviatoricStress (lattice), “deviatoric_stress” ) ; // TODO: nie ma tego
} ;
/
http://www.palabos.org/documentation/userguide/appendix-functions.html
do wyboru są:
computeDensity(lattice, domain)
Compute the density rho on each cell.
computeRhoBar(lattice, domain)
Compute the value rhoBar on each cell.
computeKineticEnergy(lattice, domain)
Compute the kinetic energy (half the squared velocity-norm) on each cell.
computeVelocityNorm(lattice, domain)
Compute the velocity-norm on each cell.
computeVelocityComponent(lattice, domain, iComponent)
Compute the component iComponent of the velocity (iComponent ranges from 0 to 1 in 2D, and from 0 to 2 in 3D) on each cell.
computeVelocity(lattice, domain)
Compute all components of the velocity on each cell.
computeDeviatoricStress(lattice, domain)
Compute the off-equilibrium part of the stress tensor Pi on each cell. The result is model-dependent and is determined by the dynamics object defined at the cell location. The result is stored in a 3-component (2D) or 6-component (3D) tensor-field (the number of components is reduced because the tensor is symmetric). To know the index of a specific component, use a notation like SymmetricTensorImpl<T,3>::xz for the the xy-component of a 3D symmetric tensor.
computeStrainRateFromStress(lattice, domain)
Compute the strain rate S=1/2(grad(u)+grad(u)^T), assuming that it is proportional to the deviatoric stress. The result is model-dependent and is determined by the dynamics object defined at the cell location. The result is stored in a 3-component (2D) or 6-component (3D) tensor-field (the number of components is reduced because the tensor is symmetric). To know the index of a specific component, use a notation like SymmetricTensorImpl<T,3>::xz for the the xy-component of a 3D symmetric tensor.
computePopulation(lattice, domain, iPop)
*/
/*********************************************************************************/
int main(int argc, char* argv[])
{
/*
Linie odtąd do początku kodu PP odpowiadają liniom 217-219 z pliku poiseulle.cpp
Musiałem dopisać parsowanie argumentów użytkownika (–block_size_x itd) i wczytywanie obrazka.
*/
plbInit(&argc, &argv);
// Linia poniżej jest moja, służy do parsowania wszystkich opcji linii poleceń, np. --block_size_x itd.
UserOptions opts(argc, argv) ; // parse user command line arguments
global::directories().setOutputDir("./tmp/");
/*
Odtąd jest wczytywanie obrazka.
*/
string file_path;
png::image< png::rgb_pixel > image(opts.file_path());
plint width = image.get_width () ; // Oba w pikselach
plint height = image.get_height() ;
plint ZSIZE = opts.channel_depth() +2 ; // +2, bo z góry i z dołu będzie przykryte warstwami odbijającymi (bounceback)
/*
Koniec wczytywania obrazka
*/
/*
Parametry brzegowe
*/
//IncomprFlowParam parameters(LatticeU, reynoldsNumber, resolution, Lx, Ly, Lz);
IncomprFlowParam parameters( opts.max_velocity(), opts.Re(), opts.resolution(), 0, 0, 0);
writeFullLogFile(parameters);
/*
To poniżej jest potrzebne do podzielenia sceny na kwadratowe poddomeny, żeby się
przypadek ładnie zrównoleglał. W pouseuille.cpp oni nie mają takich problemów, bo u nich
wszystkie komórki sceny są komórkami z płynem i trzeba na nich wykonywać obliczenia.
U nas komórek z płynem jest raptem kilka procent, dlatego trzeba wyrzucić wszystkie niepotrzebne
komórki typu bounceback, co pozwoli na znaczne (10x) zmniejszenie zużycia pamięci i przyspieszenie
obliczeń.
*/
MultiScalarField3D flagMatrix(width, height, ZSIZE);
setToFunction(flagMatrix, flagMatrix.getBoundingBox(), FluidNodes(image, ZSIZE));
plint envelopeWidth = 1;
MultiBlockManagement3D sparseBlockManagement =
computeSparseManagement (
*plb::reparallelize(flagMatrix, opts.block_size_x(), opts.block_size_y(), ZSIZE), // rozmiar bloku po Z rowny wysokości
envelopeWidth );
/*
To poniżej odpowiada liniom 238-240 z pliku poiseulle.cpp, ale znów różnice wynikają stąd,
że my mamy dużo nieużywanych komórek.
*/
// Instantiate the multi-block, based on the created block distribution and
// on default parameters.
MultiBlockLattice3D<T, DESCRIPTOR> lattice (
sparseBlockManagement,
defaultMultiBlockPolicy3D().getBlockCommunicator(),
defaultMultiBlockPolicy3D().getCombinedStatistics(),
defaultMultiBlockPolicy3D().getMultiCellAccess<T,DESCRIPTOR>(),
new BGKdynamics<T,DESCRIPTOR>( parameters.getOmega() )
);
pcout << getMultiBlockInfo(lattice) << std::endl; // To tylko dla mnie
/*
Cała reszta to samo, co linie 241-273 z pliku pouseuille.cpp, ale z wymienionymi poniżej różnicami.
/
// TODO: zupelnie nie rozumiem tych warunków brzegowych
// patrz http://www.palabos.org/documentation/userguide/boundary-conditions.html
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>
boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();
/*
Różnice w stosunku do pouseuille.cpp:
- nie ma obrazka (parametr image)
- ja w funkcji channelSetup() ustawiam na sztywno prędkości i gęstość, a w pliku pouseuille.cpp
prędkości i gęstość są liczone automagicznie obiektem klasy IncomprFlowParam
*/
channelSetup(lattice, image, parameters, *boundaryCondition, opts.inlet_velocity(), opts.rho());
/*
U nas liczbe kroków i momenty zapisów definiuje się w linii poleceń, w pouseuille.cpp są jakoś zaszyte w obiekcie klasy IncomprFlowParam
*/
plint maxT = opts.n_time_steps () ;
plint imageIter = opts.save_gif_steps() ;
plint vtkIter = opts.save_vtk_steps() ;
pcout << "Starting simulation: " << maxT << " time steps, saving gif interval " << imageIter << " time steps, saving VTK interval " << vtkIter << " time steps\n" ;
// Main loop over time iterations.
for (plint iT=0; iT<=maxT; ++iT)
{
if ( opts.save_gif(iT) ) {
pcout << "Saving Gif at time step " << iT << endl;
writeGifs(lattice, iT);
}
if ( opts.save_vtk(iT) ) {
pcout << "Saving VTK at time step " << iT << endl;
writeVTK(lattice, iT, opts);
}
lattice.collideAndStream();
}
delete boundaryCondition;
}
file command_options.hpp
[code="cpp"]
#ifndef __COMMAND_OPTIONS_HPP__
#define __COMMAND_OPTIONS_HPP__
#include <string>
class UserOptions
{
public:
UserOptions(int argc, char ** argv) ;
~UserOptions() ;
std::string file_path () const ;
int block_size_x () const ;
int block_size_y () const ;
int n_time_steps () const ;
int resolution () const ;
int save_gif_steps () const ;
int channel_depth () const ;
int save_vtk_steps () const ;
double Re () const ;
double inlet_velocity () const ;
double max_velocity () const ;
double rho () const ;
bool save_gif (int time_step) const ;
bool save_vtk (int time_step) const ;
bool save_vtk_velocity () const ;
bool save_vtk_density () const ;
bool save_vtk_rho_bar () const ;
bool save_vtk_kinetic_energy () const ;
bool save_vtk_velocity_norm () const ;
bool save_vtk_strain_rate_from_stress () const ;
private:
std::string _file_path ;
int _block_size_x ;
int _block_size_y ;
int _n_time_steps ;
int _resolution ;
int _save_gif_steps ;
int _channel_depth ;
int _save_vtk_steps ;
bool _save_vtk_option; // czy podano opcje --save_vtk_steps
double _Re ;
double _inlet_velocity ;
double _max_velocity ;
double _rho ;
bool _save_vtk_velocity ;
bool _save_vtk_density ;
bool _save_vtk_rho_bar ;
bool _save_vtk_kinetic_energy ;
bool _save_vtk_velocity_norm ;
bool _save_vtk_strain_rate_from_stress ;
} ;
#endif
file command_options.cpp
[code=“cpp”]
#include “command_options.hpp”
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include
#include “palabos2D.h”
#include “palabos2D.hh”
using namespace std ;
enum VTK_SAVE_SWITCHES {
SAVE_VTK_VELOCITY = 0 ,
SAVE_VTK_DENSITY ,
SAVE_VTK_RHO_BAR ,
SAVE_VTK_KINETIC_ENERGY ,
SAVE_VTK_VELOCITY_NORM ,
SAVE_VTK_STRAIN_RATE_FROM_STRESS ,
SAVE_VTK_ALL
} ;
UserOptions::
UserOptions (int argc, char ** argv)
{
_block_size_x = 15 ; // Wartości domyślne
_block_size_y = 15 ;
_n_time_steps = 1 ;
_resolution = 0 ; // to chyba spowoduje błąd - i dobrze.
_save_gif_steps = 0 ;
_channel_depth = 1 ;
_save_vtk_steps = 0 ;
_save_vtk_option = false ;
_Re = 0.0 ;
_inlet_velocity = 0.0 ;
_max_velocity = 0.0 ;
_rho = 0.0 ;
_save_vtk_velocity = false ;
_save_vtk_density = false ;
_save_vtk_rho_bar = false ;
_save_vtk_kinetic_energy = false ;
_save_vtk_velocity_norm = false ;
_save_vtk_strain_rate_from_stress = false ;
while (1)
{
static struct option long_options[] =
{
{"file" , required_argument, 0, 'f' },
{"block_size_x" , required_argument, 0, 'x' },
{"block_size_y" , required_argument, 0, 'y' },
{"n_time_steps" , required_argument, 0, 'n' },
{"resolution" , required_argument, 0, 'r' },
{"save_gif_steps", required_argument, 0, 's' },
{"channel_depth" , required_argument, 0, 'd' },
{"save_vtk_steps", required_argument, 0, 'v' },
{"Re" , required_argument, 0, 'R' },
{"inlet_velocity", required_argument, 0, 'i' },
{"max_velocity" , required_argument, 0, 'm' },
{"rho" , required_argument, 0, 'h' },
{"save_vtk_velocity" , optional_argument, 0, SAVE_VTK_VELOCITY },
{"save_vtk_density" , optional_argument, 0, SAVE_VTK_DENSITY },
{"save_vtk_rho_bar" , optional_argument, 0, SAVE_VTK_RHO_BAR },
{"save_vtk_kinetic_energy" , optional_argument, 0, SAVE_VTK_KINETIC_ENERGY },
{"save_vtk_velocity_norm" , optional_argument, 0, SAVE_VTK_VELOCITY_NORM },
{"save_vtk_strain_rate_from_stress", optional_argument, 0, SAVE_VTK_STRAIN_RATE_FROM_STRESS },
{"save_vtk_all" , optional_argument, 0, SAVE_VTK_ALL },
{0, 0, 0, 0}
};
char c = getopt_long (argc, argv, "f:", long_options, NULL);
/* Detect the end of the options. */
if (c == -1)
break;
switch (c)
{
case 'f':
//printf ("option -f with value `%s'\n", optarg);
_file_path = optarg ;
break;
case 'x':
//printf ("option -x with value `%d'\n", atoi(optarg));
_block_size_x = atoi(optarg) ;
break;
case 'y':
//printf ("option -y with value `%d'\n", atoi(optarg));
_block_size_y = atoi(optarg) ;
break;
case 'n':
_n_time_steps = atoi(optarg) ;
break ;
case 'r':
_resolution = atoi(optarg) ;
break ;
case 's':
_save_gif_steps = atoi(optarg) ;
break ;
case 'd':
_channel_depth = atoi(optarg) ;
break ;
case 'v':
_save_vtk_steps = atoi(optarg) ;
_save_vtk_option = true ;
break ;
case 'R':
_Re = atof(optarg) ; break ;
case 'i':
_inlet_velocity = atof(optarg) ; break ;
case 'm':
_max_velocity = atof(optarg) ; break ;
case 'h':
_rho = atof(optarg) ; break ;
#define PROCESS_OPTIONAL_ARGUMENT(option_name)
if (NULL == optarg || string(“no”) != optarg )
{
option_name = true ;
} else {
option_name = false ;
} ;
case SAVE_VTK_VELOCITY:
PROCESS_OPTIONAL_ARGUMENT( _save_vtk_velocity ) ;
break ;
case SAVE_VTK_DENSITY:
PROCESS_OPTIONAL_ARGUMENT( _save_vtk_density ) ;
break ;
case SAVE_VTK_RHO_BAR:
PROCESS_OPTIONAL_ARGUMENT( _save_vtk_rho_bar ) ;
break ;
case SAVE_VTK_KINETIC_ENERGY:
PROCESS_OPTIONAL_ARGUMENT( _save_vtk_kinetic_energy ) ;
break ;
case SAVE_VTK_VELOCITY_NORM:
PROCESS_OPTIONAL_ARGUMENT( _save_vtk_velocity_norm ) ;
break ;
case SAVE_VTK_STRAIN_RATE_FROM_STRESS:
PROCESS_OPTIONAL_ARGUMENT( _save_vtk_strain_rate_from_stress ) ;
break ;
#define SET_VTK_OPTIONS( val )
_save_vtk_velocity = (val) ;
_save_vtk_density = (val) ;
_save_vtk_rho_bar = (val) ;
_save_vtk_kinetic_energy = (val) ;
_save_vtk_velocity_norm = (val) ;
_save_vtk_strain_rate_from_stress = (val) ;
case SAVE_VTK_ALL:
if (NULL == optarg || string("no") != optarg )
{
SET_VTK_OPTIONS( true ) ;
} else {
SET_VTK_OPTIONS( false ) ;
} ;
break ;
default:
abort ();
}
}
if (false == _save_vtk_option) _save_vtk_steps = _n_time_steps ;
if ( "" == _file_path)
{
plb::pcout << "\nPodaj sciezke do pliku z obrazkiem.\n\n" ;
exit(1) ;
} ;
} ;
UserOptions::
~UserOptions()
{} ;
bool UserOptions::
save_gif( int time_step ) const
{
if (0 == _save_gif_steps ) return false ;
if (0 == time_step % _save_gif_steps ) return true ;
return false ;
} ;
bool UserOptions::
save_vtk( int time_step ) const
{
if (0 == _save_vtk_steps ) return false ;
if (0 == time_step % _save_vtk_steps ) return true ;
return false ;
} ;
#define RET_METHOD_CONST(type,name)
type UserOptions::
name() const
{
return _ ## name ;
} ;
RET_METHOD_CONST( string, file_path )
#define INT_RET_METHOD_CONST(name) RET_METHOD_CONST(int,name)
INT_RET_METHOD_CONST( block_size_x )
INT_RET_METHOD_CONST( block_size_y )
INT_RET_METHOD_CONST( n_time_steps )
INT_RET_METHOD_CONST( resolution )
INT_RET_METHOD_CONST( save_gif_steps )
INT_RET_METHOD_CONST( channel_depth )
INT_RET_METHOD_CONST( save_vtk_steps )
#define BOOL_RET_METHOD_CONST(name) RET_METHOD_CONST(bool,name)
BOOL_RET_METHOD_CONST( save_vtk_velocity )
BOOL_RET_METHOD_CONST( save_vtk_density )
BOOL_RET_METHOD_CONST( save_vtk_rho_bar )
BOOL_RET_METHOD_CONST( save_vtk_kinetic_energy )
BOOL_RET_METHOD_CONST( save_vtk_velocity_norm )
BOOL_RET_METHOD_CONST( save_vtk_strain_rate_from_stress )
#define DOUBLE_RET_METHOD_CONST(name) RET_METHOD_CONST(double,name)
DOUBLE_RET_METHOD_CONST( Re )
DOUBLE_RET_METHOD_CONST( inlet_velocity )
DOUBLE_RET_METHOD_CONST( max_velocity )
DOUBLE_RET_METHOD_CONST( rho )
file Makefile
[code="make"]
palabosRoot = ../../palabos/palabos-v1.3r0/
# Name of source files in current directory to compile and link with Palabos
projectFiles = lbm3D_png.cpp command_options.cpp
# Set optimization flags on/off
optimize = true
# Set debug mode and debug flags on/off
debug = true
# Set profiling flags on/off
profile = false
# Set MPI-parallel mode on/off (parallelism in cluster-like environment)
MPIparallel = true
# Set SMP-parallel mode on/off (shared-memory parallelism)
SMPparallel = false
# Decide whether to include calls to the POSIX API. On non-POSIX systems,
# including Windows, this flag must be false, unless a POSIX environment is
# emulated (such as with Cygwin).
usePOSIX = true
# Path to external libraries (other than Palabos)
libraryPaths = -L/usr/lib/x86_64-linux-gnu
# Path to inlude directories (other than Palabos)
includePaths = -I/usr/include/libpng12 -I../common
# Dynamic and static libraries (other than Palabos)
libraries = -lpng12
# Compiler to use without MPI parallelism
serialCXX = g++
# Compiler to use with MPI parallelism
parallelCXX = mpicxx
# General compiler flags (e.g. -Wall to turn on all warnings on g++)
compileFlags = -Wall -Wnon-virtual-dtor
# General linker flags (don't put library includes into this flag)
linkFlags =
# Compiler flags to use when optimization mode is on
optimFlags = -O3
# Compiler flags to use when debug mode is on
debugFlags = -g
# Compiler flags to use when profile mode is on
profileFlags = -pg
##########################################################################
# All code below this line is just about forwarding the options
# to SConstruct. It is recommended not to modify anything there.
##########################################################################
SCons = $(palabosRoot)/scons/scons.py -j 2 -f $(palabosRoot)/SConstruct
SConsArgs = palabosRoot=$(palabosRoot) \
projectFiles="$(projectFiles)" \
optimize=$(optimize) \
debug=$(debug) \
profile=$(profile) \
MPIparallel=$(MPIparallel) \
SMPparallel=$(SMPparallel) \
usePOSIX=$(usePOSIX) \
serialCXX=$(serialCXX) \
parallelCXX=$(parallelCXX) \
compileFlags="$(compileFlags)" \
linkFlags="$(linkFlags)" \
optimFlags="$(optimFlags)" \
debugFlags="$(debugFlags)" \
profileFlags="$(profileFlags)" \
libraryPaths="$(libraryPaths)" \
includePaths="$(includePaths)" \
libraries="$(libraries)"
compile:
python $(SCons) $(SConsArgs)
clean:
python $(SCons) -c $(SConsArgs)
/bin/rm -vf `find $(palabosRoot) -name '*~'`
bash script running simulation
#!/bin/bash
#
# Skrypt uruchamiający program lbm3D_png.
#
# Skrypt służy przede wszystkim do pamiętania wszystkich opcji programu lbm3D_png.
#
# Dodatkowo, jeśli do skryptu zostaną przekazane argumenty, to są przekazywane
# bezpośrednio do programu. Mogą więc być użyte do nadpisania wartości ustawionych
# w skrypcie.
#
# Ścieżka do skryptu jest potrzebna, jeśli byłby uruchamiany z innego katalogu niż bieżący.
script_path=${0%/*}
LBM_SIM="$script_path/lbm3D_png"
# UWAGA - nie jestem pewien, czy opcje są przetwarzane w kolejności
# ich pojawiania się.
#
# --save_gif_steps=0 wyłącza całkowicie zapis do plików gif. Jeśli ta opcja nie
# jest podana, zapis gif jest jeden na koncu symulacji. Zapis do plików gif można
# wymusić podając jako jeden z argumentów skryptu --save_gif_steps=10
RUN_COMMAND="mpirun -np 2 $LBM_SIM \
--file=rys6.png \
--channel_depth=5 \
--block_size_x=10 \
--block_size_y=10 \
--n_time_steps=5000 \
--resolution=15 \
--Re=0.56 \
--inlet_velocity=0.0125 \
--max_velocity=0.025 \
--rho=1 \
--save_gif_steps=1 \
--save_vtk_velocity=yes \
--save_vtk_density=no \
--save_vtk_rho_bar=no \
--save_vtk_kinetic_energy=no \
--save_vtk_velocity_norm=no \
--save_vtk_strain_rate_from_stress=no \
$*"
echo "running $RUN_COMMAND"
$RUN_COMMAND
image used as case definition