/*
* 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 = 1.0;
const int N = 5; // resolution of the model
const T lx = 22 * cylinderD; // length of the channel
const T ly = 4.1 * cylinderD; // height of the channel
const T Re = 100.; // Reynolds number
const T uLat = 0.05; // lattice velocity
const T maxPhysT = 60.; // max. simulation time in s, SI unit
void prepareGeometry(Grid2D& grid)
{
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 bounce back boundaries
{
const Vector wallExtend {lx+physSpacing, physSpacing};
const Vector wallOrigin {-physSpacing/2, -physSpacing/2};
IndicatorCuboid2D lowerWall(wallExtend, wallOrigin);
sGeometry.rename(1,2,lowerWall);
IndicatorCuboid2D upperWall(wallExtend, wallOrigin + Vector {0,ly-physSpacing/2});
sGeometry.rename(1,2,upperWall);
}
// Set material number for inflow and outflow
{
const Vector extend { physSpacing/2, ly-physSpacing/2};
const Vector origin {-physSpacing/4, -physSpacing/4};
IndicatorCuboid2D inflow(extend, origin);
sGeometry.rename(1,3,inflow);
IndicatorCuboid2D outflow(extend, origin + Vector {lx,0});
sGeometry.rename(1,4,outflow);
}
// Set material number for vertically centered cylinder
{
const Vector origin {2*cylinderD, 2*cylinderD};
IndicatorCircle2D obstacle(origin, cylinderD/2);
sGeometry.rename(1,5,obstacle);
}
sGeometry.clean();
sGeometry.innerClean();
sGeometry.checkForErrors();
sGeometry.print();
clout << "Prepare Geometry ... OK" << std::endl;
}
void disableRefinedArea(Grid2D& coarseGrid,
RefiningGrid2D& fineGrid)
{
auto& sGeometry = coarseGrid.getSuperGeometry();
auto refinedOverlap = fineGrid.getRefinedOverlap();
sGeometry.rename(1,0,*refinedOverlap);
sGeometry.rename(2,0,*refinedOverlap);
sGeometry.rename(5,0,*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();
createLocalBoundaryCondition2D(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, 5});
sLattice.defineRhoU(materials, rho0, u0);
sLattice.iniEquilibrium(materials, rho0, u0);
sLattice.initialize();
clout << "Prepare lattice ... OK" << std::endl;
}
void setBoundaryValues(Grid2D& grid, int iT)
{
auto& converter = grid.getConverter();
auto& sGeometry = grid.getSuperGeometry();
auto& sLattice = grid.getSuperLattice();
const int iTmaxStart = converter.getLatticeTime(0.2*maxPhysT);
const int iTupdate = 10;
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() * frac[0];
const T radius = ly/2;
std::vector axisPoint{0, ly/2};
std::vector axisDirection{1, 0};
Poiseuille2D u(axisPoint, axisDirection, maxVelocity, radius);
sLattice.defineU(sGeometry, 3, u);
}
}
void getResults(const std::string& prefix,
Grid2D& grid,
int iT)
{
OstreamManager clout(std::cout,"getResults");
auto& converter = grid.getConverter();
auto& sLattice = grid.getSuperLattice();
auto& sGeometry = grid.getSuperGeometry();
SuperVTMwriter2D vtmWriter(prefix + "cylinder2d");
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();
}
if (iT%100==0) {
vtmWriter.write(iT);
}
}
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, N, LatticeVelocity(uLat), Re);
prepareGeometry(coarseGrid);
const auto coarseDeltaX = coarseGrid.getConverter().getPhysDeltaX();
const Vector fineExtend {8*cylinderD, ly-2*coarseDeltaX};
const Vector fineOrigin {0.75*cylinderD, coarseDeltaX};
auto& fineGrid = coarseGrid.refine(fineOrigin, fineExtend);
prepareGeometry(fineGrid);
const Vector fineExtendB {10*coarseDeltaX, ly};
const Vector fineOriginB {lx-10*coarseDeltaX, 0};
auto& fineOutflowGrid = coarseGrid.refine(fineOriginB, fineExtendB, false);
prepareGeometry(fineOutflowGrid);
{
const Vector origin = fineOutflowGrid.getOrigin();
const Vector extend = fineOutflowGrid.getExtend();
const Vector extendY = {0,extend[1]};
coarseGrid.addFineCoupling(fineOutflowGrid, origin, extendY);
coarseGrid.addCoarseCoupling(fineOutflowGrid, origin + Vector {coarseDeltaX,0}, extendY);
IndicatorCuboid2D refined(extend, origin + Vector {2*coarseDeltaX,0});
coarseGrid.getSuperGeometry().rename(1,0,refined);
coarseGrid.getSuperGeometry().rename(2,0,refined);
coarseGrid.getSuperGeometry().rename(4,0,refined);
}
disableRefinedArea(coarseGrid, fineGrid);
const Vector fineExtend2 {4*cylinderD, 2*cylinderD};
const Vector fineOrigin2 {1*cylinderD, 2*cylinderD-fineExtend2[1]/2};
auto& fineGrid2 = fineGrid.refine(fineOrigin2, fineExtend2);
prepareGeometry(fineGrid2);
disableRefinedArea(fineGrid, fineGrid2);
prepareLattice(coarseGrid);
prepareLattice(fineGrid);
prepareLattice(fineOutflowGrid);
prepareLattice(fineGrid2);
clout << "Starting simulation..." << endl;
Timer timer(
coarseGrid.getConverter().getLatticeTime(maxPhysT),
coarseGrid.getSuperGeometry().getStatistics().getNvoxel());
timer.start();
const int statIter = coarseGrid.getConverter().getLatticeTime(0.5);
for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) {
setBoundaryValues(coarseGrid, iT);
coarseGrid.collideAndStream();
getResults("level0_", coarseGrid, iT);
getResults("level1_", fineGrid, iT);
getResults("level1_outflow_", fineOutflowGrid, iT);
getResults("level2_", fineGrid2, iT);
if (iT%statIter == 0) {
timer.update(iT);
timer.printStep();
coarseGrid.getSuperLattice().getStatistics().print(iT, coarseGrid.getConverter().getPhysTime(iT));
}
}
timer.stop();
timer.printSummary();
}