/* Lattice Boltzmann sample, written in C++, using the OpenLB
* library
*
* Copyright (C) 2007, 2012 Jonas Latt, Mathias J. Krause
* Vojtech Cvrcek, Peter Weisbrod
* E-mail contact: info@openlb.net
* The most recent release of OpenLB can be downloaded at
*
*
* This program 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 2
* of the License, or (at your option) any later version.
*
* This program 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, write to the Free
* Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301, USA.
*/
#include "olb2D.h"
#include "olb2D.hh"
#include
#include
#include
#include
#include
using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace std;
typedef double T;
#define DESCRIPTOR D2Q9Descriptor
const T lx = 2.; // length of the channel
const T ly = 1.; // height of the channel
int N = 20; // resolution of the model
const T Re = 10.; // Reynolds number
const T maxPhysT = 20.; // max. simulation time in s, SI unit
const T physInterval = 0.25; // interval for the convergence check in s
const T residuum = 1e-5; // residuum for the convergence check
const T tuner = 0.99; // for partialSlip only: 0->bounceBack, 1->freeSlip
void prepareGeometry(UnitConverter const& converter,
SuperGeometry2D& superGeometry)
{
OstreamManager clout(std::cout,"prepareGeometry");
clout << "Prepare Geometry ..." << std::endl;
superGeometry.rename(0,2);
superGeometry.rename(2,1,0,1);
Vector extend;
Vector origin;
T physSpacing = converter.getPhysDeltaX();
// Set material number for inflow
extend[1] = ly;
extend[0] = physSpacing / 2;
origin[0] -= physSpacing / 4;
IndicatorCuboid2D inflow(extend, origin);
superGeometry.rename(1,3,1,inflow);
// Set material number for outflow
origin[0] = lx - physSpacing / 4;
IndicatorCuboid2D outflow(extend, origin);
superGeometry.rename(1,4,1,outflow);
// Removes all not needed boundary voxels outside the surface
superGeometry.clean();
// Removes all not needed boundary voxels inside the surface
superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();
clout << "Prepare Geometry ... OK" << std::endl;
}
void prepareCoarseLattice(UnitConverter const& converter,
SuperLattice2D& sLattice,
Dynamics& bulkDynamics,
sOnLatticeBoundaryCondition2D& sBoundaryCondition,
SuperGeometry2D& superGeometry)
{
OstreamManager clout(std::cout,"prepareLattice");
clout << "Prepare coarse lattice ..." << std::endl;
const T omega = converter.getLatticeRelaxationFrequency();
sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics());
sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); // bulk
sLattice.defineDynamics(superGeometry, 2, &instances::getBounceBack());
sLattice.defineDynamics(superGeometry, 3, &bulkDynamics); // inflow
sBoundaryCondition.addVelocityBoundary(superGeometry, 3, omega);
T Lx = converter.getLatticeLength(lx);
T Ly = converter.getLatticeLength(ly);
T p0 =8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly);
AnalyticalLinear2D rho(-p0/lx*DESCRIPTOR::invCs2, 0, p0*DESCRIPTOR::invCs2+1);
AnalyticalConst2D iniRho(1);
AnalyticalConst2D iniU(std::vector {0.,0.});
sLattice.defineRhoU(superGeometry, 1, iniRho, iniU);
sLattice.iniEquilibrium(superGeometry, 1, iniRho, iniU);
sLattice.defineRhoU(superGeometry, 2, iniRho, iniU);
sLattice.iniEquilibrium(superGeometry, 2, iniRho, iniU);
T maxVelocity = converter.getCharLatticeVelocity();
T distance2Wall = converter.getConversionFactorLength();
Poiseuille2D u(superGeometry, 3, maxVelocity, distance2Wall);
sLattice.defineRhoU(superGeometry, 3, rho, u);
sLattice.iniEquilibrium(superGeometry, 3, rho, u);
sLattice.initialize();
clout << "Prepare Lattice ... OK" << std::endl;
}
void prepareFineLattice(UnitConverter const& converter,
SuperLattice2D& sLattice,
Dynamics& bulkDynamics,
sOnLatticeBoundaryCondition2D& sBoundaryCondition,
SuperGeometry2D& superGeometry)
{
OstreamManager clout(std::cout,"prepareLattice");
clout << "Prepare Lattice ..." << std::endl;
const T omega = converter.getLatticeRelaxationFrequency();
sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics());
sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); // bulk
sLattice.defineDynamics(superGeometry, 2, &instances::getBounceBack());
sLattice.defineDynamics(superGeometry, 4, &bulkDynamics); // outflow
sBoundaryCondition.addPressureBoundary(superGeometry, 4, omega);
T Lx = converter.getLatticeLength(lx);
T Ly = converter.getLatticeLength(ly);
T p0 =8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly);
AnalyticalLinear2D rho(-p0/lx*DESCRIPTOR::invCs2, 0, p0*DESCRIPTOR::invCs2+1);
AnalyticalConst2D iniRho(1);
AnalyticalConst2D iniU(std::vector {0.,0.});
sLattice.defineRhoU(superGeometry, 1, iniRho, iniU);
sLattice.iniEquilibrium(superGeometry, 1, iniRho, iniU);
sLattice.defineRhoU(superGeometry, 2, iniRho, iniU);
sLattice.iniEquilibrium(superGeometry, 2, iniRho, iniU);
sLattice.defineRhoU(superGeometry, 4, iniRho, iniU);
sLattice.iniEquilibrium(superGeometry, 4, iniRho, iniU);
sLattice.initialize();
clout << "Prepare Lattice ... OK" << std::endl;
}
void getResults(const std::string& prefix,
SuperLattice2D& sLattice,
UnitConverter const& converter, int iT,
SuperGeometry2D& superGeometry, Timer& timer, bool hasConverged)
{
OstreamManager clout(std::cout,"getResults");
SuperVTMwriter2D vtmWriter(prefix + "poiseuille2d");
SuperLatticePhysVelocity2D velocity(sLattice, converter);
SuperLatticePhysPressure2D pressure(sLattice, converter);
SuperLatticeGeometry2D geometry( sLattice, superGeometry );
vtmWriter.addFunctor(geometry);
vtmWriter.addFunctor(velocity);
vtmWriter.addFunctor(pressure);
const int vtmIter = converter.getLatticeTime(maxPhysT/100.);
const int statIter = converter.getLatticeTime(maxPhysT/10.);
if ( iT==0 ) {
vtmWriter.createMasterFile();
}
if (iT%vtmIter==0 || hasConverged) {
vtmWriter.write(iT);
}
if (iT%statIter==0 || hasConverged) {
timer.update(iT);
timer.printStep();
sLattice.getStatistics().print(iT,converter.getPhysTime(iT));
}
}
template class DESCRIPTOR>
class Grid {
private:
std::unique_ptr> _converter;
std::unique_ptr> _cuboids;
std::unique_ptr> _balancer;
std::unique_ptr> _geometry;
std::unique_ptr> _lattice;
public:
Grid(IndicatorF2D& domainF, int resolution, T tau, int re):
_converter(new UnitConverterFromResolutionAndRelaxationTime(
resolution, // resolution: number of voxels per charPhysL
tau, // latticeRelaxationTime: relaxation time, have to be greater than 0.5!
T{1}, // charPhysLength: reference length of simulation geometry
T{1}, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
T{1./re}, // physViscosity: physical kinematic viscosity in __m^2 / s__
T{1})), // physDensity: physical density in __kg / m^3__
_cuboids(new CuboidGeometry2D(
domainF,
_converter->getConversionFactorLength(),
1)),
_balancer(new HeuristicLoadBalancer(
*_cuboids)),
_geometry(new SuperGeometry2D(
*_cuboids,
*_balancer,
2)),
_lattice(new SuperLattice2D(
*_geometry))
{
_converter->print();
}
UnitConverter& getConverter()
{
return *_converter;
}
CuboidGeometry2D& getCuboidGeometry()
{
return *_cuboids;
}
LoadBalancer& getLoadBalancer()
{
return *_balancer;
}
SuperGeometry2D& getSuperGeometry()
{
return *_geometry;
}
SuperLattice2D& getSuperLattice()
{
return *_lattice;
}
};
int main(int argc, char* argv[])
{
olbInit(&argc, &argv);
singleton::directories().setOutputDir("./tmp/");
OstreamManager clout(std::cout,"main");
const T coarseDeltaX = 1./N;
Vector coarseOrigin { 0.0 , 0.0 };
Vector coarseExtend { 0.5 * lx, ly };
IndicatorCuboid2D coarseCuboid(coarseExtend, coarseOrigin);
Vector fineOrigin { 0.5 * lx - coarseDeltaX, 0.0 };
Vector fineExtend { 0.5 * lx + coarseDeltaX, ly };
IndicatorCuboid2D fineCuboid(fineExtend, fineOrigin);
Grid coarseGrid(
coarseCuboid,
N,
0.8,
Re);
Grid fineGrid(
fineCuboid,
2*coarseGrid.getConverter().getResolution(),
2.0*coarseGrid.getConverter().getLatticeRelaxationTime() - 0.5,
Re);
prepareGeometry(coarseGrid.getConverter(), coarseGrid.getSuperGeometry());
prepareGeometry(fineGrid.getConverter(), fineGrid.getSuperGeometry());
Dynamics* coarseBulkDynamics;
coarseBulkDynamics = new BGKdynamics(
coarseGrid.getConverter().getLatticeRelaxationFrequency(),
instances::getBulkMomenta());
sOnLatticeBoundaryCondition2D coarseSBoundaryCondition(coarseGrid.getSuperLattice());
createLocalBoundaryCondition2D(coarseSBoundaryCondition);
prepareCoarseLattice(
coarseGrid.getConverter(),
coarseGrid.getSuperLattice(),
*coarseBulkDynamics,
coarseSBoundaryCondition,
coarseGrid.getSuperGeometry());
Dynamics* fineBulkDynamics;
fineBulkDynamics = new BGKdynamics(
fineGrid.getConverter().getLatticeRelaxationFrequency(),
instances::getBulkMomenta());
sOnLatticeBoundaryCondition2D fineSBoundaryCondition(fineGrid.getSuperLattice());
createLocalBoundaryCondition2D(fineSBoundaryCondition);
prepareFineLattice(
fineGrid.getConverter(),
fineGrid.getSuperLattice(),
*fineBulkDynamics,
fineSBoundaryCondition,
fineGrid.getSuperGeometry());
clout << "starting simulation..." << endl;
Timer timer(
coarseGrid.getConverter().getLatticeTime(maxPhysT),
coarseGrid.getSuperGeometry().getStatistics().getNvoxel());
util::ValueTracer converge(
coarseGrid.getConverter().getLatticeTime(physInterval),
residuum);
timer.start();
for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); iT += 2) {
if (converge.hasConverged()) {
clout << "Simulation converged." << endl;
getResults(
"coarse_",
coarseGrid.getSuperLattice(),
coarseGrid.getConverter(),
iT,
coarseGrid.getSuperGeometry(),
timer,
converge.hasConverged());
break;
}
coarseGrid.getSuperLattice().collideAndStream();
getResults(
"coarse_",
coarseGrid.getSuperLattice(),
coarseGrid.getConverter(),
iT / 2,
coarseGrid.getSuperGeometry(),
timer,
converge.hasConverged());
fineGrid.getSuperLattice().collideAndStream();
getResults(
"fine_",
fineGrid.getSuperLattice(),
fineGrid.getConverter(),
iT,
fineGrid.getSuperGeometry(),
timer,
converge.hasConverged());
fineGrid.getSuperLattice().collideAndStream();
getResults(
"fine_",
fineGrid.getSuperLattice(),
fineGrid.getConverter(),
iT+1,
fineGrid.getSuperGeometry(),
timer,
converge.hasConverged());
converge.takeValue(coarseGrid.getSuperLattice().getStatistics().getAverageEnergy(), true);
}
timer.stop();
timer.printSummary();
}