Hello everybody,
I combined the grid refinement with the poiseuille showcase.
I can run a simulation with a poiseuille velocity profile and a grid refinement. When I look at the solution I get a velocity in y-direction in the area of the corners of the grid refinement.
As I cannot find a mistake in the code, I would be glad if you take a look at it.
Thanks a lot.
Cheers
Konstantin
The code:
#include "palabos2D.h"
#ifndef PLB_PRECOMPILED // Unless precompiled version is used,
#include "palabos2D.hh" // include full template code
#endif
#include <vector>
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
using namespace plb;
using namespace plb::descriptors;
using namespace std;
typedef double T;
#define DESCRIPTOR D2Q9Descriptor
/// Velocity on the parabolic Poiseuille profile
// With this formula the width of the channel has to be 1
T poiseuilleVelocity(plint iY, IncomprFlowParam<T> const& parameters) {
T y = (T)iY / parameters.getResolution();
return 4.*parameters.getLatticeU() * (y-y*y);
}
/// Linearly decreasing pressure profile
T poiseuillePressure(plint iX, IncomprFlowParam<T> const& parameters) {
T Lx = parameters.getNx()-1;
T Ly = parameters.getNy()-1;
return 8.*parameters.getLatticeNu()*parameters.getLatticeU() / (Ly*Ly) * (Lx/(T)2-(T)iX);
}
/// Convert pressure to density according to ideal gas law
T poiseuilleDensity(plint iX, IncomprFlowParam<T> const& parameters) {
return poiseuillePressure(iX,parameters)*DESCRIPTOR<T>::invCs2 + (T)1;
}
/// A functional, used to initialize the velocity for the boundary conditions
template<typename T>
class PoiseuilleVelocity {
public:
PoiseuilleVelocity(IncomprFlowParam<T> parameters_)
: parameters(parameters_)
{ }
void operator()(plint iX, plint iY, Array<T,2>& u) const {
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
}
private:
IncomprFlowParam<T> parameters;
};
/// A functional, used to initialize the density for the boundary conditions
template<typename T>
class PoiseuilleDensity {
public:
PoiseuilleDensity(IncomprFlowParam<T> parameters_)
: parameters(parameters_)
{ }
T operator()(plint iX, plint iY) const {
return poiseuilleDensity(iX,parameters);
}
private:
IncomprFlowParam<T> parameters;
};
/// A functional, used to create an initial condition for with zero velocity,
/// and linearly decreasing pressure.
template<typename T>
class PoiseuilleDensityAndZeroVelocity {
public:
PoiseuilleDensityAndZeroVelocity(IncomprFlowParam<T> parameters_)
: parameters(parameters_)
{ }
void operator()(plint iX, plint iY, T& rho, Array<T,2>& u) const {
rho = poiseuilleDensity(iX,parameters);
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
}
private:
IncomprFlowParam<T> parameters;
};
// Initialisation of the domain.
/// A functional, used to instantiate bounce-back nodes at the locations of the cylinder
//Change to MultiGridLattice2D
//+Convective Refinement Parameters - data wrapper ConvectiveRefinementParameters instead of IncomprFlowParam<T>
void cylinderSetup( MultiGridLattice2D<T,DESCRIPTOR>& lattice,
ConvectiveRefinementParameters<T> const& parameters,
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )
{
//for each lattice inside the multigrid
for (plint iLevel=0; iLevel<lattice.getNumLevels(); ++iLevel){
const plint nx = parameters.getParameters(iLevel).getNx();
const plint ny = parameters.getParameters(iLevel).getNy();
Box2D inlet(0, 0, 2, ny-2);
Box2D outlet(nx-1, nx-1, 1, ny-2);
//Velocity boundary condition on the bottom wall
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice.getComponent(iLevel), Box2D(0, nx-1, 0, 0) );
//Velocity boundary condition on the top wall
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice.getComponent(iLevel), Box2D(0, nx-1, ny-1, ny-1) );
//Velocity boundary condition at the inlet
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice.getComponent(iLevel), inlet );
//Pressure boundary condition at the outlet
boundaryCondition.setPressureConditionOnBlockBoundaries (lattice.getComponent(iLevel), outlet );
//Define the value of the imposed velocity on all nodes which have previously been
// defined to be velocity boundary nodes.
setBoundaryVelocity(lattice.getComponent(iLevel),
lattice.getComponent(iLevel).getBoundingBox(), PoiseuilleVelocity<T>(parameters.getParameters(iLevel) ) );
//Define the value of the imposed density at the outlet
setBoundaryDensity(lattice.getComponent(iLevel),
Box2D(nx-1, nx-1, 1, ny-2), PoiseuilleDensity<T>(parameters.getParameters(iLevel) ) );
// Initialize all cells at an equilibrium distribution, with a velocity and density
// value of the analytical Poiseuille solution.
initializeAtEquilibrium (lattice.getComponent(iLevel),
lattice.getComponent(iLevel).getBoundingBox(),
PoiseuilleDensityAndZeroVelocity<T>(parameters.getParameters(iLevel) ) );
}
lattice.initialize();
}
/// Produce a VTK snapshot of the velocity-norm.
void writeVTKs(MultiBlockLattice2D<T,DESCRIPTOR> lattice, IncomprFlowParam<T> const& parameters, plint iter, string fName)
{
T dx = parameters.getDeltaX();
T dt = parameters.getDeltaT();
VtkImageOutput2D<T> vtkOut(createFileName(fName, iter, 6), dx);
vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", dx/dt);
vtkOut.writeData<2,float>(*computeVelocity(lattice), "velocity", dx/dt);
vtkOut.writeData<float>(*computeDensity(lattice), "density", dx/dt);
}
template<typename T, template<typename U> class Descriptor>
void writeVTKs(MultiGridLattice2D<T,Descriptor> const& lattice, ConvectiveRefinementParameters<T> const& parameters, plint iter) {
for (plint iLevel=0; iLevel<lattice.getNumLevels(); ++iLevel) {
std::stringstream fname("level_");
fname << "level_" << iLevel << "_";
writeVTKs(lattice.getComponent(iLevel), parameters.getParameters(iLevel), iter, fname.str());
}
}
int main(int argc, char* argv[]) {
plbInit(&argc, &argv);
global::directories().setOutputDir("./tmp/");
//NEW
global::setDefaultMultiScaleManager(new ConvectiveMultiScaleManager());
plint referenceLevel = 0;
plint behaviorLevel = referenceLevel;
plint numLevel = 2;
plint N=100;
IncomprFlowParam<T> baseParameters(
(T) 0.001, // uMax
(T) 100., // Re
N, // N
6., // lx
1. // ly
);
const T logT = (T)0.02;
const T imSave = (T)1000.1;
const T vtkSave = (T)0.05;
const T maxT = (T)5.001;
writeLogFile(baseParameters, "GridRefinement and Initialization with the poiseuille profil");
// we choose a convective refinement dx=dt
ConvectiveRefinementParameters<T> convectiveParameters (numLevel,referenceLevel,baseParameters );
Dynamics<T,DESCRIPTOR>* backgroundDynamics = new BGKdynamics<T,DESCRIPTOR>(baseParameters.getOmega());
// main object for the grid refinement
MultiGridManagement2D management(baseParameters.getNx(),baseParameters.getNy(), numLevel, referenceLevel);
// adding refinements to management
management.refine(0,Box2D(2*N,4*N,N/3,2*N/3)); // Block in the middle (lenghth 2N - height N/3)
// Parallelization. Here, the domain is parallelized along the x-direction. The position
// of the different slices is explicitly created, and provided in units of the
// finest lattice.
std::vector<plint> parallelizeIn;
plint procNumber = global::mpi().getSize();
std::vector<std::pair<plint,plint> > ranges;
plint finestNy = convectiveParameters[numLevel-1].getNy()-1;
util::linearRepartition(0,finestNy,procNumber,ranges);
for (plint iBlock=0; iBlock<(plint)ranges.size()-1; ++iBlock){
parallelizeIn.push_back(ranges[iBlock].second);
}
management.parallelizeX(parallelizeIn);
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>*
boundaryCondition = createLocalBoundaryCondition2D<T,DESCRIPTOR>();
MultiGridLattice2D<T,DESCRIPTOR>
lattice( management,
defaultMultiGridPolicy2D().getBlockCommunicator<T>(management.getNumLevels()),
defaultMultiGridPolicy2D().getCombinedStatistics(management.getNumLevels()),
backgroundDynamics, behaviorLevel );
cylinderSetup(lattice, convectiveParameters, *boundaryCondition);
// output file
plb_ofstream outfile("out.dat");
// Main loop over time iterations.
for (plint iT=0; iT*baseParameters.getDeltaT()<maxT; ++iT) {
if (iT%baseParameters.nStep(vtkSave)==0) {
pcout << "Saving VTK ..." << endl;
writeVTKs(lattice, convectiveParameters, iT);
pcout << endl;
}
if (iT%baseParameters.nStep(logT)==0) {
pcout << "step " << iT
<< "; t=" << iT*baseParameters.getDeltaT();
}
// Lattice Boltzmann iteration step.
lattice.collideAndStream();
// At this point, the state of the lattice corresponds to the
// discrete time iT+1, and the stored averages are upgraded to time iT.
if (iT%baseParameters.nStep(logT)==0) {
pcout << "; av energy ="
<< setprecision(10) << getStoredAverageEnergy<T>(lattice)
<< "; av rho ="
<< getStoredAverageDensity<T>(lattice) << endl;
pcout << "; Omega =" << baseParameters.getOmega() << endl;
}
}
delete boundaryCondition;
}