my program with many of errors!
#include “xflows3D.h”
#ifndef XFL_PRECOMPILED // Unless precompiled version is used,
#include “xflows3D.hh” // include full template code
#endif
#include
#include
#include
#include
#include
using namespace xfl;
using namespace xfl::descriptors;
using namespace std;
typedef double T;
#define DESCRIPTOR D2Q9Descriptor
/// Velocity on the parabolic Poiseuille profile
T poiseuilleVelocity(xint iY, xint iZ, IncomprFlowParam const& parameters) {
T y = (T)iY / parameters.getResolution();
T z = (T)iZ / parameters.getResolution();
return 4.parameters.getLatticeU() * (y-yy) * (z-z*z);
}
/// Linearly decreasing pressure profile
T poiseuillePressure(xint iX, xint iZ,IncomprFlowParam const& parameters) {
T Lx = parameters.getNx()-1;
T Ly = parameters.getNy()-1;
T Lz = parameters.getNz()-1;
return 8.*parameters.getLatticeNu()*parameters.getLatticeU() / (Ly*Ly) * (Lz*Lz) * (Lx/(T)2-(T)iX);
}
/// Convert pressure to density according to ideal gas law
T poiseuilleDensity(xint iX, xint iZ, IncomprFlowParam const& parameters) {
return poiseuillePressure(iX,parameters)*DESCRIPTOR::invCs2 + (T)1;
}
/// A functional, used to initialize the density for the boundary conditions
template
class PoiseuilleDensity {
public:
PoiseuilleDensity(IncomprFlowParam parameters_)
: parameters(parameters_)
{ }
T operator()(xint iX, xint iY, xint iZ) const {
return poiseuilleDensity(iX,parameters);
}
private:
IncomprFlowParam parameters;
};
/// A functional, used to create an initial condition for with zero velocity,
/// and linearly decreasing pressure.
template
class PoiseuilleDensityAndZeroVelocity {
public:
PoiseuilleDensityAndZeroVelocity(IncomprFlowParam parameters_)
: parameters(parameters_)
{ }
void operator()(xint iX, xint iY,xint iZ, T& rho, Array<T,3>& u) const {
rho = poiseuilleDensity(iX,parameters);
u[0] = T();
u[1] = T();
u[2] = T();
}
private:
IncomprFlowParam parameters;
};
void channelSetup( BlockLatticeBase3D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters,
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition )
{
const xint nx = parameters.getNx();
const xint ny = parameters.getNy();
const xint nz = parameters.getNz();
// Velocity boundary condition on bottom wall.
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box3D(0, nx-1, 0, 0,0,0) );
// Velocity boundary condition on top wall.
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box3D(0, nx-1, ny-1, ny-1,nz-1,nz-1) );
// Pressure resp. velocity boundary condition on the inlet and outlet.
if (inletOutlet == pressure) {
// Note: pressure boundary conditions are currently implemented
// only for edges of the boundary, and not for corner nodes.
boundaryCondition.setPressureConditionOnBlockBoundaries (
lattice, Box3D(0,0, 1,ny-2, 0,0) );
boundaryCondition.setPressureConditionOnBlockBoundaries (
lattice, Box3D(nx-1,nx-1, 1,ny-2,nz-1,nz-1) );
}
else {
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box3D(0,0, 1,ny-2, 0, 0) );
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box3D(nx-1,nx-1, 1,ny-2,nz-1,nz-1) );
}
setBoundaryDensity (
lattice, lattice.getBoundingBox(),
PoiseuilleDensity<T>(parameters) );
setBoundaryVelocity (
lattice, lattice.getBoundingBox(),
PoiseuilleVelocity<T>(parameters) );
initializeAtEquilibrium (
lattice, lattice.getBoundingBox(),
PoiseuilleDensityAndZeroVelocity<T>(parameters) );
// Call initialize to get the lattice ready for the simulation.
lattice.initialize();
}
/// Produce a GIF snapshot of the velocity-norm.
template
void writeGifs(BlockLatticeT& lattice,
IncomprFlowParam const& parameters, xint iter)
{
const xint imSize = 600;
const xint nx = parameters.getNx();
const xint ny = parameters.getNy();
const xint nz = parameters.getNz();
Box3D slice(0, nx-1, 0, ny-1, nz/2, nz/2);
ImageWriter<T> imageWriter("leeloo");
imageWriter.writeScaledGif( createFileName("uz", iter, 6),
*computeVelocityComponent (lattice, slice, zComponent),
imSize, imSize );
imageWriter.writeScaledGif( createFileName("uNorm", iter, 6),
*computeVelocityNorm (lattice, slice),
imSize, imSize );
imageWriter.writeScaledGif( createFileName("omega", iter, 6),
*computeNorm(*computeVorticity (
*computeVelocity(lattice, slice)) ),
imSize, imSize );
}
template
void writeVTK(BlockLatticeT& lattice,
IncomprFlowParam const& parameters, xint iter)
{
T dx = parameters.getDeltaX();
T dt = parameters.getDeltaT();
VtkImageOutput3D vtkOut(createFileName(“vtk”, iter, 6), dx);
vtkOut.writeData<3,float>(*computeVelocity(lattice), “velocity”, dx/dt);
vtkOut.writeData<3,float>(*computeVorticity(*computeVelocity(lattice)), “vorticity”, 1./dt);
}
int main(int argc, char* argv[]) {
xflInit(&argc, &argv);
global::directories().setOutputDir("./tmp/");
IncomprFlowParam<T> parameters(
(T) 2e-2, // uMax
(T) 5., // Re
60, // N
3., // lx
1. // ly
);
const T logT = (T)0.1;
const T imSave = (T)0.5;
const T vtkSave = (T)2.;
const T maxT = (T)5.1;
writeLogFile(parameters, "3D Poiseuille flow");
MultiBlockLattice3D<T, DESCRIPTOR> lattice (
parameters.getNx(), parameters.getNy(), parameters.getNz(),
new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>*
//boundaryCondition = createInterpBoundaryCondition3D<T,DESCRIPTOR>();
boundaryCondition = createLocalBoundaryCondition3D<T,DESCRIPTOR>();
channelSetup(lattice, parameters, *boundaryCondition);
// Main loop over time iterations.
for (xint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
if (iT%parameters.nStep(imSave)==0) {
xcout << "Saving Gif ..." << endl;
writeGif(lattice, iT);
}
if (iT%parameters.nStep(vtkSave)==0 && iT>0) {
xcout << "Saving VTK file ..." << endl;
writeVTK(lattice, parameters, iT);
}
// Lattice Boltzmann iteration step.
lattice.collideAndStream();
}
delete boundaryCondition;
}