diff options
| -rw-r--r-- | apps/adrian/cylinder2d/Makefile | 105 | ||||
| -rw-r--r-- | apps/adrian/cylinder2d/cylinder2d.cpp | 286 | ||||
| -rw-r--r-- | apps/adrian/cylinder2d/definitions.mk | 30 | ||||
| -rw-r--r-- | apps/adrian/cylinder2d/module.mk | 29 | ||||
| -rw-r--r-- | apps/adrian/poiseuille2d/poiseuille2d.cpp | 48 |
5 files changed, 458 insertions, 40 deletions
diff --git a/apps/adrian/cylinder2d/Makefile b/apps/adrian/cylinder2d/Makefile new file mode 100644 index 0000000..a953954 --- /dev/null +++ b/apps/adrian/cylinder2d/Makefile @@ -0,0 +1,105 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias Krause +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# <http://www.openlb.net/> +# +# 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. + +########################################################################### +## definitions + +include definitions.mk + +include $(ROOT)/global.mk + +OBJECTS := $(foreach file, $(SRC), $(PWD)/$(file:.cpp=.o)) +DEPS := $(foreach file, $(SRC), $(PWD)/$(file:.cpp=.d)) + +########################################################################### +## all + +all : depend compile link + + +########################################################################### +## dependencies + +depend : $(DEPS) + +$(PWD)/%.d : %.cpp + @echo Create dependencies for $< + @$(SHELL) -ec '$(CXX) -M $(CXXFLAGS) $(IDIR) $< \ + | sed -e "s!$*\.o!$(PWD)\/$*\.o!1" > .tmpfile; \ + cp -f .tmpfile $@;' + +########################################################################### +## compile + +compile : $(OBJECTS) + +$(PWD)/%.o: %.cpp + @echo Compile $< + $(CXX) $(CXXFLAGS) $(IDIR) -c $< -o $@ + +########################################################################### +## clean + +clean : cleanrub cleanobj cleandep + +cleanrub: + @echo Clean rubbish files + @rm -f *~ core .tmpfile tmp/*.* $(OUTPUT) + @rm -f tmp/vtkData/*.* tmp/vtkData/data/*.* tmp/imageData/*.* tmp/gnuplotData/*.* tmp/gnuplotData/data/*.* + +cleanobj: + @echo Clean object files + @rm -f $(OBJECTS) + +cleandep: + @echo Clean dependencies files + @rm -f $(DEPS) + +cleanbuild: + @cd $(ROOT); \ + $(MAKE) cleanlib; + +########################################################################### +## update lib + +$(ROOT)/$(LIBDIR)/lib$(LIB).a : + @cd $(ROOT); \ + $(MAKE) all + +########################################################################### +## link + +link: $(OUTPUT) + +$(OUTPUT): $(OBJECTS) $(ROOT)/$(LIBDIR)/lib$(LIB).a + @echo Link $@ + $(CXX) $(foreach file, $(SRC), $(file:.cpp=.o)) $(LDFLAGS) -L$(ROOT)/$(LIBDIR) -l$(LIB) -lz -o $@ + +########################################################################### +## include dependencies + +ifneq "$(strip $(wildcard *.d))" "" + include $(foreach file,$(DEPS),$(file)) +endif + +########################################################################### +########################################################################### 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 + * <http://www.openlb.net/> + * + * 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 <vector> + +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<T,DESCRIPTOR>& 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<T,2> wallExtend {lx+physSpacing, physSpacing/2}; + const Vector<T,2> wallOrigin {-physSpacing/4, -physSpacing/4}; + + IndicatorCuboid2D<T> lowerWall(wallExtend, wallOrigin); + sGeometry.rename(1,2,lowerWall); + + IndicatorCuboid2D<T> upperWall(wallExtend, wallOrigin + Vector<T,2> {0,ly}); + sGeometry.rename(1,2,upperWall); + } + + // Set material number for inflow and outflow + { + const Vector<T,2> extend { physSpacing/2, ly}; + const Vector<T,2> origin {-physSpacing/4, -physSpacing/4}; + + IndicatorCuboid2D<T> inflow(extend, origin); + sGeometry.rename(1,3,inflow); + + IndicatorCuboid2D<T> outflow(extend, origin + Vector<T,2> {lx,0}); + sGeometry.rename(1,4,outflow); + } + + // Set material number for vertically centered obstacle + { + const Vector<T,2> origin {1.25, ly/2}; + IndicatorCircle2D<T> 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<T,DESCRIPTOR>& grid, + Dynamics<T, DESCRIPTOR>& bulkDynamics, + sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& 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<T, DESCRIPTOR>()); + sLattice.defineDynamics(sGeometry, 1, &bulkDynamics); // bulk + sLattice.defineDynamics(sGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>()); + 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<T,T> rho(-p0/lx*DESCRIPTOR<T>::invCs2, 0, p0*DESCRIPTOR<T>::invCs2+1); + + const T maxVelocity = converter.getCharLatticeVelocity(); + const T radius = ly/2; + std::vector<T> axisPoint{0, ly/2}; + std::vector<T> axisDirection{1, 0}; + Poiseuille2D<T> 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<T,DESCRIPTOR>& grid, + int iT, Timer<T>& timer, bool hasConverged) +{ + OstreamManager clout(std::cout,"getResults"); + + auto& converter = grid.getConverter(); + auto& sLattice = grid.getSuperLattice(); + auto& sGeometry = grid.getSuperGeometry(); + + SuperVTMwriter2D<T> vtmWriter(prefix + "cylinder2d"); + SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity(sLattice, converter); + SuperLatticePhysPressure2D<T, DESCRIPTOR> pressure(sLattice, converter); + SuperLatticeGeometry2D<T, DESCRIPTOR> 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<T,2> coarseOrigin {0.0, 0.0}; + const Vector<T,2> coarseExtend {lx, ly}; + IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin); + + Grid2D<T,DESCRIPTOR> coarseGrid(coarseCuboid, N, baseTau, Re); + prepareGeometry(coarseGrid); + + const Vector<T,2> fineExtend {3.0, 1.5}; + const Vector<T,2> 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<T, DESCRIPTOR> coarseBulkDynamics( + coarseGrid.getConverter().getLatticeRelaxationFrequency(), + instances::getBulkMomenta<T, DESCRIPTOR>()); + + sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> coarseSBoundaryCondition(coarseGrid.getSuperLattice()); + createLocalBoundaryCondition2D<T, DESCRIPTOR>(coarseSBoundaryCondition); + + prepareLattice(coarseGrid, coarseBulkDynamics, coarseSBoundaryCondition); + + BGKdynamics<T, DESCRIPTOR> fineBulkDynamics( + fineGrid.getConverter().getLatticeRelaxationFrequency(), + instances::getBulkMomenta<T, DESCRIPTOR>()); + + sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition(fineGrid.getSuperLattice()); + createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition); + + const Vector<T,2> fineExtend2 {0.6, 0.4}; + const Vector<T,2> 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<T, DESCRIPTOR> fineBulkDynamics2( + fineGrid2.getConverter().getLatticeRelaxationFrequency(), + instances::getBulkMomenta<T, DESCRIPTOR>()); + + sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition2(fineGrid2.getSuperLattice()); + createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition2); + + prepareLattice(fineGrid2, fineBulkDynamics2, fineSBoundaryCondition2); + + clout << "starting simulation..." << endl; + Timer<T> timer( + coarseGrid.getConverter().getLatticeTime(maxPhysT), + coarseGrid.getSuperGeometry().getStatistics().getNvoxel()); + util::ValueTracer<T> 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(); +} diff --git a/apps/adrian/cylinder2d/definitions.mk b/apps/adrian/cylinder2d/definitions.mk new file mode 100644 index 0000000..6584906 --- /dev/null +++ b/apps/adrian/cylinder2d/definitions.mk @@ -0,0 +1,30 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias Krause +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# <http://www.openlb.net/> +# +# 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. + + +########################################################################### +########################################################################### +## DEFINITIONS TO BE CHANGED + +ROOT := ../../.. +SRC := cylinder2d.cpp +OUTPUT := cylinder2d diff --git a/apps/adrian/cylinder2d/module.mk b/apps/adrian/cylinder2d/module.mk new file mode 100644 index 0000000..1190482 --- /dev/null +++ b/apps/adrian/cylinder2d/module.mk @@ -0,0 +1,29 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2017 Markus Mohrhard +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# <http://www.openlb.net/> +# +# 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. + +current_dir := $(dir $(word $(words $(MAKEFILE_LIST)),$(MAKEFILE_LIST))) + +include global.mk +include rules.mk +include $(addsuffix definitions.mk, $(current_dir)) + +$(eval $(call sample,$(current_dir)$(OUTPUT),$(addprefix $(current_dir), $(SRC)))) diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp index 7112f58..b6c2f52 100644 --- a/apps/adrian/poiseuille2d/poiseuille2d.cpp +++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp @@ -37,10 +37,10 @@ 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 T lx = 4.0; // length of the channel +const T ly = 1.0; // height of the channel const int N = 50; // resolution of the model -const T Re = 200.; // Reynolds number +const T Re = 100.; // 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 @@ -83,13 +83,6 @@ void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid) sGeometry.rename(1,4,outflow); } - // Set material number for vertically centered obstacle - { - const Vector<T,2> origin {1.25, ly/2}; - IndicatorCircle2D<T> obstacle(origin, 0.1); - sGeometry.rename(1,2,obstacle); - } - sGeometry.clean(); sGeometry.innerClean(); sGeometry.checkForErrors(); @@ -170,7 +163,7 @@ void getResults(const std::string& prefix, vtmWriter.createMasterFile(); } - if (iT%20==0) { + if (iT%100==0) { vtmWriter.write(iT); } @@ -195,8 +188,8 @@ int main(int argc, char* argv[]) Grid2D<T,DESCRIPTOR> coarseGrid(coarseCuboid, N, baseTau, Re); prepareGeometry(coarseGrid); - const Vector<T,2> fineExtend {3.0, 1.5}; - const Vector<T,2> fineOrigin {0.8, (ly-fineExtend[1])/2}; + const Vector<T,2> fineExtend {2.0, 0.5}; + const Vector<T,2> fineOrigin {1.0, (ly-fineExtend[1])/2}; auto& fineGrid = coarseGrid.refine(fineOrigin, fineExtend); prepareGeometry(fineGrid); @@ -221,33 +214,14 @@ int main(int argc, char* argv[]) sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition(fineGrid.getSuperLattice()); createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition); - const Vector<T,2> fineExtend2 {0.6, 0.4}; - const Vector<T,2> 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<T, DESCRIPTOR> fineBulkDynamics2( - fineGrid2.getConverter().getLatticeRelaxationFrequency(), - instances::getBulkMomenta<T, DESCRIPTOR>()); - - sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition2(fineGrid2.getSuperLattice()); - createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition2); - - prepareLattice(fineGrid2, fineBulkDynamics2, fineSBoundaryCondition2); - clout << "starting simulation..." << endl; Timer<T> timer( coarseGrid.getConverter().getLatticeTime(maxPhysT), coarseGrid.getSuperGeometry().getStatistics().getNvoxel()); util::ValueTracer<T> converge( - fineGrid2.getConverter().getLatticeTime(physInterval), + fineGrid.getConverter().getLatticeTime(physInterval), residuum); timer.start(); @@ -271,14 +245,8 @@ int main(int argc, char* argv[]) iT, timer, converge.hasConverged()); - getResults( - "fine2_", - fineGrid2, - iT, - timer, - converge.hasConverged()); - converge.takeValue(fineGrid2.getSuperLattice().getStatistics().getAverageEnergy(), true); + converge.takeValue(fineGrid.getSuperLattice().getStatistics().getAverageEnergy(), true); } timer.stop(); |