/*
* Lattice Boltzmann grid refinement sample, written in C++,
* using the OpenLB library
*
* Copyright (C) 2019 Adrian Kummerländer
* 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"
#ifndef OLB_PRECOMPILED
#include "olb2D.hh"
#endif
#include
using namespace olb;
typedef double T;
#define DESCRIPTOR descriptors::D2Q9Descriptor
/// Setup geometry relative to cylinder diameter as defined by [SchaeferTurek96]
const T cylinderD = 0.1;
const int N = 10; // resolution of the cylinder
const T deltaR = cylinderD / N; // coarse lattice spacing
const T lx = 22*cylinderD + deltaR; // length of the channel
const T ly = 4.1*cylinderD + deltaR; // height of the channel
const T cylinderX = 2*cylinderD;
const T cylinderY = 2*cylinderD + deltaR/2;
const T Re = 100.; // Reynolds number
const T tau = 0.51; // relaxation time
const T maxPhysT = 16.; // max. simulation time in s, SI unit
const Characteristics PhysCharacteristics(
0.1, // char. phys. length
1.0, // char. phys. velocity
0.1/Re, // phsy. kinematic viscosity
1.0); // char. phys. density
void prepareGeometry(Grid2D& grid, Vector origin, Vector extend)
{
OstreamManager clout(std::cout,"prepareGeometry");
clout << "Prepare Geometry ..." << std::endl;
auto& converter = grid.getConverter();
auto& sGeometry = grid.getSuperGeometry();
sGeometry.rename(0,1);
const T physSpacing = converter.getPhysDeltaX();
// Set material number for channel walls
{
const Vector wallExtend { extend[0]+physSpacing/2, physSpacing/2 };
const Vector wallOrigin = origin - physSpacing/4;
IndicatorCuboid2D lowerWall(wallExtend, wallOrigin);
sGeometry.rename(1,2,lowerWall);
}
{
const Vector wallExtend { extend[0]+physSpacing/2, physSpacing/2 };
const Vector wallOrigin { origin[0]-physSpacing/4, extend[1]-physSpacing/4 };
IndicatorCuboid2D upperWall(wallExtend, wallOrigin);
sGeometry.rename(1,2,upperWall);
}
// Set material number for inflow and outflow
{
const Vector inflowExtend { physSpacing/2, extend[1]+physSpacing/4 };
const Vector inflowOrigin = origin - physSpacing/4;
IndicatorCuboid2D inflow(inflowExtend, inflowOrigin);
sGeometry.rename(1,3,inflow);
}
{
const Vector outflowExtend { physSpacing/2, extend[1]+physSpacing/4 };
const Vector outflowOrigin { extend[0]-physSpacing/4, origin[0]-physSpacing/4 };
IndicatorCuboid2D outflow(outflowExtend, outflowOrigin);
sGeometry.rename(1,4,outflow);
}
// Set material number for vertically centered cylinder
{
const Vector cylinderOrigin = origin + Vector {cylinderX, cylinderY};
IndicatorCircle2D obstacle(cylinderOrigin, cylinderD/2);
sGeometry.rename(1,5,obstacle);
}
sGeometry.clean();
sGeometry.innerClean();
sGeometry.checkForErrors();
clout << "Prepare Geometry ... OK" << std::endl;
}
void disableRefinedArea(Grid2D& coarseGrid,
RefiningGrid2D& fineGrid)
{
auto& sGeometry = coarseGrid.getSuperGeometry();
auto refinedOverlap = fineGrid.getRefinedOverlap();
sGeometry.reset(*refinedOverlap);
}
void prepareLattice(Grid2D& grid)
{
OstreamManager clout(std::cout,"prepareLattice");
clout << "Prepare lattice ..." << std::endl;
auto& converter = grid.getConverter();
auto& sGeometry = grid.getSuperGeometry();
auto& sLattice = grid.getSuperLattice();
Dynamics& bulkDynamics = grid.addDynamics(
std::unique_ptr>(
new BGKdynamics(
grid.getConverter().getLatticeRelaxationFrequency(),
instances::getBulkMomenta())));
sOnLatticeBoundaryCondition2D& sBoundaryCondition = grid.getOnLatticeBoundaryCondition();
createInterpBoundaryCondition2D(sBoundaryCondition);
const T omega = converter.getLatticeRelaxationFrequency();
sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics());
sLattice.defineDynamics(sGeometry, 1, &bulkDynamics); // bulk
sLattice.defineDynamics(sGeometry, 2, &bulkDynamics); // walls
sLattice.defineDynamics(sGeometry, 3, &bulkDynamics); // inflow
sLattice.defineDynamics(sGeometry, 4, &bulkDynamics); // outflow
sLattice.defineDynamics(sGeometry, 5, &instances::getBounceBack()); // cylinder
sBoundaryCondition.addVelocityBoundary(sGeometry, 2, omega);
sBoundaryCondition.addVelocityBoundary(sGeometry, 3, omega);
sBoundaryCondition.addPressureBoundary(sGeometry, 4, omega);
AnalyticalConst2D rho0(1.0);
AnalyticalConst2D u0(0.0, 0.0);
auto materials = sGeometry.getMaterialIndicator({1, 2, 3, 4});
sLattice.defineRhoU(materials, rho0, u0);
sLattice.iniEquilibrium(materials, rho0, u0);
sLattice.initialize();
clout << "Prepare lattice ... OK" << std::endl;
sGeometry.print();
}
void setBoundaryValues(Grid2D& grid, int iT)
{
auto& converter = grid.getConverter();
auto& sGeometry = grid.getSuperGeometry();
auto& sLattice = grid.getSuperLattice();
const int iTmaxStart = converter.getLatticeTime(0.4*16);
const int iTupdate = 5;
if ( iT % iTupdate == 0 && iT <= iTmaxStart ) {
PolynomialStartScale StartScale(iTmaxStart, 1);
T iTvec[1] { T(iT) };
T frac[1] { };
StartScale(frac, iTvec);
const T maxVelocity = converter.getCharLatticeVelocity() * 3./2. * frac[0];
Poiseuille2D u(sGeometry, 3, maxVelocity, deltaR/2);
sLattice.defineU(sGeometry, 3, u);
}
}
void getResults(Grid2D& grid,
const std::string& prefix,
int iT)
{
auto& converter = grid.getConverter();
auto& sLattice = grid.getSuperLattice();
auto& sGeometry = grid.getSuperGeometry();
SuperVTMwriter2D vtmWriter(prefix);
SuperLatticePhysVelocity2D velocity(sLattice, converter);
SuperLatticePhysPressure2D pressure(sLattice, converter);
SuperLatticeGeometry2D geometry(sLattice, sGeometry);
SuperLatticeKnudsen2D knudsen(sLattice);
vtmWriter.addFunctor(geometry);
vtmWriter.addFunctor(velocity);
vtmWriter.addFunctor(pressure);
vtmWriter.addFunctor(knudsen);
if (iT==0) {
vtmWriter.createMasterFile();
}
vtmWriter.write(iT);
}
void takeMeasurements(Grid2D& grid)
{
static T maxDrag = 0.0;
OstreamManager clout(std::cout,"measurement");
auto& sLattice = grid.getSuperLattice();
auto& sGeometry = grid.getSuperGeometry();
auto& converter = grid.getConverter();
SuperLatticePhysPressure2D pressure(sLattice, converter);
AnalyticalFfromSuperF2D intpolatePressure(pressure, true);
SuperLatticePhysDrag2D dragF(sLattice, sGeometry, 5, converter);
const T radiusCylinder = cylinderD/2;
const T point1[2] { cylinderX - radiusCylinder, cylinderY };
const T point2[2] { cylinderX + radiusCylinder, cylinderY };
T pressureInFrontOfCylinder, pressureBehindCylinder;
intpolatePressure(&pressureInFrontOfCylinder, point1);
intpolatePressure(&pressureBehindCylinder, point2);
T pressureDrop = pressureInFrontOfCylinder - pressureBehindCylinder;
clout << "pressureDrop=" << pressureDrop;
const int input[3] {};
T drag[dragF.getTargetDim()] {};
dragF(drag, input);
if (drag[0] > maxDrag) {
maxDrag = drag[0];
};
clout << "; drag=" << drag[0] << "; maxDrag: " << maxDrag << "; lift=" << drag[1] << endl;
}
int main(int argc, char* argv[])
{
olbInit(&argc, &argv);
singleton::directories().setOutputDir("./tmp/");
OstreamManager clout(std::cout,"main");
const Vector coarseOrigin {0.0, 0.0};
const Vector coarseExtend {lx, ly};
IndicatorCuboid2D coarseCuboid(coarseExtend, coarseOrigin);
Grid2D coarseGrid(
coarseCuboid,
RelaxationTime(tau),
N,
PhysCharacteristics);
const Vector domainOrigin = coarseGrid.getSuperGeometry().getStatistics().getMinPhysR(0);
const Vector domainExtend = coarseGrid.getSuperGeometry().getStatistics().getPhysExtend(0);
prepareGeometry(coarseGrid, domainOrigin, domainExtend);
const auto coarseDeltaX = coarseGrid.getConverter().getPhysDeltaX();
const Vector fineExtend {10*cylinderD, domainExtend[1]-4*coarseDeltaX};
const Vector fineOrigin {0.5*cylinderD, (domainExtend[1]-fineExtend[1])/2};
auto& fineGrid = coarseGrid.refine(fineOrigin, fineExtend);
prepareGeometry(fineGrid, domainOrigin, domainExtend);
disableRefinedArea(coarseGrid, fineGrid);
const Vector fineExtend2 {6*cylinderD, fineGrid.getExtend()[1]-4*coarseDeltaX};
const Vector fineOrigin2 {0.75*cylinderD, (domainExtend[1]-fineExtend2[1])/2};
auto& fineGrid2 = fineGrid.refine(fineOrigin2, fineExtend2);
prepareGeometry(fineGrid2, domainOrigin, domainExtend);
disableRefinedArea(fineGrid, fineGrid2);
const Vector fineExtend3 {4*cylinderD, 2*cylinderD};
const Vector fineOrigin3 {1*cylinderD, (domainExtend[1]-fineExtend3[1])/2};
auto& fineGrid3 = fineGrid2.refine(fineOrigin3, fineExtend3);
prepareGeometry(fineGrid3, domainOrigin, domainExtend);
disableRefinedArea(fineGrid2, fineGrid3);
const Vector fineExtend4 {1.25*cylinderD, 1.25*cylinderD};
const Vector fineOrigin4 {cylinderX - fineExtend4[0]/2, cylinderY - fineExtend4[1]/2};
auto& fineGrid4 = fineGrid3.refine(fineOrigin4, fineExtend4);
prepareGeometry(fineGrid4, domainOrigin, domainExtend);
disableRefinedArea(fineGrid3, fineGrid4);
coarseGrid.forEachGrid(prepareLattice);
clout << "Total number of active cells: " << coarseGrid.getActiveVoxelN() << endl;
clout << "Starting simulation..." << endl;
const int statIter = coarseGrid.getConverter().getLatticeTime(0.05);
Timer timer(
coarseGrid.getConverter().getLatticeTime(maxPhysT),
coarseGrid.getSuperGeometry().getStatistics().getNvoxel());
timer.start();
for (int iT = 0; iT <= coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) {
setBoundaryValues(coarseGrid, iT);
coarseGrid.collideAndStream();
if (iT%statIter == 0) {
timer.update(iT);
timer.printStep();
coarseGrid.forEachGrid("cylinder2d", [&](Grid2D& grid, const std::string& id) {
getResults(grid, id, iT);
});
takeMeasurements(fineGrid4);
}
}
timer.stop();
timer.printSummary();
}