From 3d447eb5e493c112682e69e5fcd5c91807fcc7fc Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Wed, 16 Jan 2019 15:06:20 +0100 Subject: Rename example to match its content, extract poiseuille2d --- apps/adrian/cylinder2d/cylinder2d.cpp | 286 ++++++++++++++++++++++++++++++++++ 1 file changed, 286 insertions(+) create mode 100644 apps/adrian/cylinder2d/cylinder2d.cpp (limited to 'apps/adrian/cylinder2d/cylinder2d.cpp') diff --git a/apps/adrian/cylinder2d/cylinder2d.cpp b/apps/adrian/cylinder2d/cylinder2d.cpp new file mode 100644 index 0000000..0cdfa7b --- /dev/null +++ b/apps/adrian/cylinder2d/cylinder2d.cpp @@ -0,0 +1,286 @@ +/* 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" +#ifndef OLB_PRECOMPILED +#include "olb2D.hh" +#endif + +#include + +using namespace olb; +using namespace olb::descriptors; + +typedef double T; + +#define DESCRIPTOR D2Q9Descriptor + +const T lx = 8.0; // length of the channel +const T ly = 2.0; // height of the channel +const int N = 50; // resolution of the model +const T Re = 200.; // Reynolds number +const T baseTau = 0.8; // Relaxation time of coarsest grid +const T maxPhysT = 60.; // 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(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/2}; + const Vector wallOrigin {-physSpacing/4, -physSpacing/4}; + + IndicatorCuboid2D lowerWall(wallExtend, wallOrigin); + sGeometry.rename(1,2,lowerWall); + + IndicatorCuboid2D upperWall(wallExtend, wallOrigin + Vector {0,ly}); + sGeometry.rename(1,2,upperWall); + } + + // Set material number for inflow and outflow + { + const Vector extend { physSpacing/2, ly}; + 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 obstacle + { + const Vector origin {1.25, ly/2}; + IndicatorCircle2D obstacle(origin, 0.1); + sGeometry.rename(1,2,obstacle); + } + + sGeometry.clean(); + sGeometry.innerClean(); + sGeometry.checkForErrors(); + sGeometry.print(); + + clout << "Prepare Geometry ... OK" << std::endl; +} + +void prepareLattice(Grid2D& grid, + Dynamics& bulkDynamics, + sOnLatticeBoundaryCondition2D& sBoundaryCondition) +{ + OstreamManager clout(std::cout,"prepareLattice"); + clout << "Prepare lattice ..." << std::endl; + + auto& converter = grid.getConverter(); + auto& sGeometry = grid.getSuperGeometry(); + auto& sLattice = grid.getSuperLattice(); + + const T omega = converter.getLatticeRelaxationFrequency(); + + sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics()); + sLattice.defineDynamics(sGeometry, 1, &bulkDynamics); // bulk + sLattice.defineDynamics(sGeometry, 2, &instances::getBounceBack()); + sLattice.defineDynamics(sGeometry, 3, &bulkDynamics); // inflow + sLattice.defineDynamics(sGeometry, 4, &bulkDynamics); // outflow + + sBoundaryCondition.addVelocityBoundary(sGeometry, 3, omega); + sBoundaryCondition.addVelocityBoundary(sGeometry, 4, omega); + + const T Lx = converter.getLatticeLength(lx); + const T Ly = converter.getLatticeLength(ly); + const T p0 = 8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly); + AnalyticalLinear2D rho(-p0/lx*DESCRIPTOR::invCs2, 0, p0*DESCRIPTOR::invCs2+1); + + const T maxVelocity = converter.getCharLatticeVelocity(); + const T radius = ly/2; + std::vector axisPoint{0, ly/2}; + std::vector axisDirection{1, 0}; + Poiseuille2D u(axisPoint, axisDirection, maxVelocity, radius); + + sLattice.defineRhoU(sGeometry, 1, rho, u); + sLattice.iniEquilibrium(sGeometry, 1, rho, u); + sLattice.defineRhoU(sGeometry, 2, rho, u); + sLattice.iniEquilibrium(sGeometry, 2, rho, u); + + sLattice.defineRhoU(sGeometry, 3, rho, u); + sLattice.iniEquilibrium(sGeometry, 3, rho, u); + sLattice.defineRhoU(sGeometry, 4, rho, u); + sLattice.iniEquilibrium(sGeometry, 4, rho, u); + + sLattice.initialize(); + + clout << "Prepare lattice ... OK" << std::endl; +} + +void getResults(const std::string& prefix, + Grid2D& grid, + int iT, Timer& timer, bool hasConverged) +{ + 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); + vtmWriter.addFunctor(geometry); + vtmWriter.addFunctor(velocity); + vtmWriter.addFunctor(pressure); + + const int statIter = converter.getLatticeTime(maxPhysT/10.); + + if (iT==0) { + vtmWriter.createMasterFile(); + } + + if (iT%20==0) { + vtmWriter.write(iT); + } + + if (iT%statIter==0 || hasConverged) { + timer.update(iT); + timer.printStep(); + + sLattice.getStatistics().print(iT,converter.getPhysTime(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, baseTau, Re); + prepareGeometry(coarseGrid); + + const Vector fineExtend {3.0, 1.5}; + const Vector fineOrigin {0.8, (ly-fineExtend[1])/2}; + + auto& fineGrid = coarseGrid.refine(fineOrigin, fineExtend); + prepareGeometry(fineGrid); + + auto refinedOverlap = fineGrid.getRefinedOverlap(); + coarseGrid.getSuperGeometry().rename(1,0,*refinedOverlap); + coarseGrid.getSuperGeometry().rename(2,0,*refinedOverlap); + + BGKdynamics coarseBulkDynamics( + coarseGrid.getConverter().getLatticeRelaxationFrequency(), + instances::getBulkMomenta()); + + sOnLatticeBoundaryCondition2D coarseSBoundaryCondition(coarseGrid.getSuperLattice()); + createLocalBoundaryCondition2D(coarseSBoundaryCondition); + + prepareLattice(coarseGrid, coarseBulkDynamics, coarseSBoundaryCondition); + + BGKdynamics fineBulkDynamics( + fineGrid.getConverter().getLatticeRelaxationFrequency(), + instances::getBulkMomenta()); + + sOnLatticeBoundaryCondition2D fineSBoundaryCondition(fineGrid.getSuperLattice()); + createLocalBoundaryCondition2D(fineSBoundaryCondition); + + const Vector fineExtend2 {0.6, 0.4}; + const Vector fineOrigin2 {1.05, (ly-fineExtend2[1])/2}; + + auto& fineGrid2 = fineGrid.refine(fineOrigin2, fineExtend2); + prepareGeometry(fineGrid2); + + auto refinedOverlap2 = fineGrid2.getRefinedOverlap(); + fineGrid.getSuperGeometry().rename(1,0,*refinedOverlap2); + fineGrid.getSuperGeometry().rename(2,0,*refinedOverlap2); + + prepareLattice(fineGrid, fineBulkDynamics, fineSBoundaryCondition); + + BGKdynamics fineBulkDynamics2( + fineGrid2.getConverter().getLatticeRelaxationFrequency(), + instances::getBulkMomenta()); + + sOnLatticeBoundaryCondition2D fineSBoundaryCondition2(fineGrid2.getSuperLattice()); + createLocalBoundaryCondition2D(fineSBoundaryCondition2); + + prepareLattice(fineGrid2, fineBulkDynamics2, fineSBoundaryCondition2); + + clout << "starting simulation..." << endl; + Timer timer( + coarseGrid.getConverter().getLatticeTime(maxPhysT), + coarseGrid.getSuperGeometry().getStatistics().getNvoxel()); + util::ValueTracer converge( + fineGrid2.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.collideAndStream(); + + getResults( + "coarse_", + coarseGrid, + iT, + timer, + converge.hasConverged()); + getResults( + "fine_", + fineGrid, + iT, + timer, + converge.hasConverged()); + getResults( + "fine2_", + fineGrid2, + iT, + timer, + converge.hasConverged()); + + converge.takeValue(fineGrid2.getSuperLattice().getStatistics().getAverageEnergy(), true); + } + + timer.stop(); + timer.printSummary(); +} -- cgit v1.2.3