/* 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
const int N = 16; // resolution of the model
const T Re = 10.; // Reynolds number
const T maxPhysT = 40.; // 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
void prepareGeometry(UnitConverter const& converter,
SuperGeometry2D& superGeometry)
{
OstreamManager clout(std::cout,"prepareGeometry");
clout << "Prepare Geometry ..." << std::endl;
superGeometry.rename(0,2);
// Set material number for bounce back boundaries
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 statIter = converter.getLatticeTime(maxPhysT/10.);
if (iT==0) {
vtmWriter.createMasterFile();
}
if (iT%100==0) {
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, has 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;
}
T getScalingFactor()
{
return 4. - _converter->getLatticeRelaxationFrequency();
}
T getInvScalingFactor()
{
return 1./getScalingFactor();
}
};
T order4interpolation(T fm3, T fm1, T f1, T f3)
{
return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3);
}
T order3interpolation(T fm1, T f1, T f3)
{
return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3;
}
template class DESCRIPTOR>
void coupleC2F(Grid& coarse, Grid& fine)
{
auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
auto& fineLattice = fine.getSuperLattice().getBlockLattice(0);
const int x = coarseLattice.getNx()-2;
for (int y=0; y < coarseLattice.getNy(); ++y) {
for (int i=0; i < DESCRIPTOR::q; ++i) {
fineLattice.get(0, 2*y)[i] = coarseLattice.get(x,y)[i];
}
}
for (int y=2; y < coarseLattice.getNy()-2; ++y) {
for (int i=0; i < DESCRIPTOR::q; ++i) {
fineLattice.get(0,1+2*y)[i] = order4interpolation(
coarseLattice.get(x,y-1)[i],
coarseLattice.get(x,y)[i],
coarseLattice.get(x,y+1)[i],
coarseLattice.get(x,y+2)[i]
);
}
}
for (int i=0; i < DESCRIPTOR::q; ++i) {
fineLattice.get(0,1)[i] = order3interpolation(
coarseLattice.get(x,0)[i],
coarseLattice.get(x,1)[i],
coarseLattice.get(x,2)[i]
);
fineLattice.get(0,fineLattice.getNy()-2)[i] = order3interpolation(
coarseLattice.get(x,coarseLattice.getNy()-1)[i],
coarseLattice.get(x,coarseLattice.getNy()-2)[i],
coarseLattice.get(x,coarseLattice.getNy()-3)[i]
);
}
}
template class DESCRIPTOR>
void computeFneq(const Cell& cell, T fNeq[DESCRIPTOR::q])
{
T rho{};
T u[2]{};
cell.computeRhoU(rho, u);
lbHelpers::computeFneq(cell, fNeq, rho, u);
}
template class DESCRIPTOR>
void computeFeq(const Cell& cell, T fEq[DESCRIPTOR::q])
{
T rho{};
T u[2]{};
cell.computeRhoU(rho, u);
const T uSqr = u[0]*u[0] + u[1]*u[1];
for (int i=0; i < DESCRIPTOR::q; ++i) {
fEq[i] = lbHelpers::equilibrium(i, rho, u, uSqr);
}
}
template class DESCRIPTOR>
T computeRestrictedFneq(const BlockLatticeView2D& lattice, int x, int y)
{
T sum = T{0};
for (int i=0; i < DESCRIPTOR::q; ++i) {
T fNeq[DESCRIPTOR::q]{};
computeFneq(lattice.get(x + DESCRIPTOR::c[i][0], y + DESCRIPTOR::c[i][1]), fNeq);
sum += fNeq[i];
}
return sum / DESCRIPTOR::q;
}
template class DESCRIPTOR>
void coupleF2C(Grid& coarse, Grid& fine)
{
auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
auto& fineLattice = fine.getSuperLattice().getBlockLattice(0);
const int x = coarseLattice.getNx()-1;
for (int y=0; y < coarseLattice.getNy(); ++y) {
T fEq[DESCRIPTOR::q]{};
computeFeq(fineLattice.get(2,2*y), fEq);
const T restrictedFneq = computeRestrictedFneq(fineLattice, 2, 2*y);
for (int i=0; i < DESCRIPTOR::q; ++i) {
coarseLattice.get(x,y)[i] = fEq[i] + coarse.getInvScalingFactor() * restrictedFneq;
}
}
}
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(
fineGrid.getConverter().getLatticeTime(physInterval),
residuum);
timer.start();
for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) {
if (converge.hasConverged()) {
clout << "Simulation converged." << endl;
break;
}
coarseGrid.getSuperLattice().collideAndStream();
fineGrid.getSuperLattice().collideAndStream();
coupleC2F(coarseGrid, fineGrid);
fineGrid.getSuperLattice().collideAndStream();
coupleC2F(coarseGrid, fineGrid);
coupleF2C(coarseGrid, fineGrid);
getResults(
"coarse_",
coarseGrid.getSuperLattice(),
coarseGrid.getConverter(),
iT,
coarseGrid.getSuperGeometry(),
timer,
converge.hasConverged());
getResults(
"fine_",
fineGrid.getSuperLattice(),
fineGrid.getConverter(),
iT,
fineGrid.getSuperGeometry(),
timer,
converge.hasConverged());
converge.takeValue(fineGrid.getSuperLattice().getStatistics().getAverageEnergy(), true);
}
timer.stop();
timer.printSummary();
}