From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Mon, 24 Jun 2019 14:43:36 +0200 Subject: Initialize at openlb-1-3 --- examples/turbulence/aorta3d/Makefile | 105 +++++++ examples/turbulence/aorta3d/aorta3d.cpp | 356 ++++++++++++++++++++++ examples/turbulence/aorta3d/aorta3d.stl | Bin 0 -> 132784 bytes examples/turbulence/aorta3d/definitions.mk | 30 ++ examples/turbulence/aorta3d/module.mk | 29 ++ examples/turbulence/nozzle3d/Makefile | 105 +++++++ examples/turbulence/nozzle3d/definitions.mk | 30 ++ examples/turbulence/nozzle3d/module.mk | 29 ++ examples/turbulence/nozzle3d/nozzle3d.cpp | 382 ++++++++++++++++++++++++ examples/turbulence/tgv3d/Makefile | 105 +++++++ examples/turbulence/tgv3d/Re1600_Brachet.inp | 73 +++++ examples/turbulence/tgv3d/Re3000_Brachet.inp | 69 +++++ examples/turbulence/tgv3d/Re800_Brachet.inp | 64 ++++ examples/turbulence/tgv3d/definitions.mk | 30 ++ examples/turbulence/tgv3d/module.mk | 29 ++ examples/turbulence/tgv3d/tgv3d.cpp | 421 +++++++++++++++++++++++++++ examples/turbulence/venturi3d/Makefile | 105 +++++++ examples/turbulence/venturi3d/definitions.mk | 30 ++ examples/turbulence/venturi3d/module.mk | 29 ++ examples/turbulence/venturi3d/venturi3d.cpp | 284 ++++++++++++++++++ examples/turbulence/venturi3d/venturi3d.xml | 41 +++ 21 files changed, 2346 insertions(+) create mode 100644 examples/turbulence/aorta3d/Makefile create mode 100644 examples/turbulence/aorta3d/aorta3d.cpp create mode 100644 examples/turbulence/aorta3d/aorta3d.stl create mode 100644 examples/turbulence/aorta3d/definitions.mk create mode 100644 examples/turbulence/aorta3d/module.mk create mode 100644 examples/turbulence/nozzle3d/Makefile create mode 100644 examples/turbulence/nozzle3d/definitions.mk create mode 100644 examples/turbulence/nozzle3d/module.mk create mode 100644 examples/turbulence/nozzle3d/nozzle3d.cpp create mode 100644 examples/turbulence/tgv3d/Makefile create mode 100644 examples/turbulence/tgv3d/Re1600_Brachet.inp create mode 100644 examples/turbulence/tgv3d/Re3000_Brachet.inp create mode 100644 examples/turbulence/tgv3d/Re800_Brachet.inp create mode 100644 examples/turbulence/tgv3d/definitions.mk create mode 100644 examples/turbulence/tgv3d/module.mk create mode 100644 examples/turbulence/tgv3d/tgv3d.cpp create mode 100644 examples/turbulence/venturi3d/Makefile create mode 100644 examples/turbulence/venturi3d/definitions.mk create mode 100644 examples/turbulence/venturi3d/module.mk create mode 100644 examples/turbulence/venturi3d/venturi3d.cpp create mode 100644 examples/turbulence/venturi3d/venturi3d.xml (limited to 'examples/turbulence') diff --git a/examples/turbulence/aorta3d/Makefile b/examples/turbulence/aorta3d/Makefile new file mode 100644 index 0000000..a953954 --- /dev/null +++ b/examples/turbulence/aorta3d/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 +# +# +# 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/examples/turbulence/aorta3d/aorta3d.cpp b/examples/turbulence/aorta3d/aorta3d.cpp new file mode 100644 index 0000000..356561b --- /dev/null +++ b/examples/turbulence/aorta3d/aorta3d.cpp @@ -0,0 +1,356 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2011-2014 Mathias J. Krause + * 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. + */ + +/* aorta3d.cpp: + * In this example the fluid flow through a bifurcation is + * simulated. The geometry is obtained from a mesh in stl-format. + * With Bouzidi boundary conditions the curved boundary is + * adequately mapped and initialized fully automatically. As + * dynamics a Smagorinsky turbulent BGK model is used to stabilize + * the simulation for low resolutions. As output the flux at the + * inflow and outflow region is computed. The results has been + * validated by comparison with other results obtained with FEM + * and FVM. + */ + + +#include "olb3D.h" +#ifndef OLB_PRECOMPILED // Unless precompiled version is used, +#include "olb3D.hh" // include full template code; +#endif +#include +#include +#include +#include + +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace olb::util; + +typedef double T; +#define DESCRIPTOR D3Q19<> + + +//simulation parameters +const int N = 40; // resolution of the model +const int M = 20; // time discretization refinement +const bool bouzidiOn = true; // choice of boundary condition +const T maxPhysT = 2.; // max. simulation time in s, SI unit + + +// Stores data from stl file in geometry in form of material numbers +void prepareGeometry( UnitConverter const& converter, IndicatorF3D& indicator, + STLreader& stlReader, SuperGeometry3D& superGeometry ) { + + OstreamManager clout( std::cout,"prepareGeometry" ); + clout << "Prepare Geometry ..." << std::endl; + + superGeometry.rename( 0,2,indicator ); + superGeometry.rename( 2,1,stlReader ); + + superGeometry.clean(); + + // Set material number for inflow + IndicatorCircle3D inflow( 0.218125 ,0.249987 ,0.0234818, 0., 1.,0., 0.0112342 ); + IndicatorCylinder3D layerInflow( inflow, 2.*converter.getConversionFactorLength() ); + superGeometry.rename( 2,3,1,layerInflow ); + + // Set material number for outflow0 + //IndicatorCircle3D outflow0(0.2053696,0.0900099,0.0346537, 2.5522,5.0294,-1.5237, 0.0054686 ); + IndicatorCircle3D outflow0( 0.2053696,0.0900099,0.0346537, 0.,-1.,0., 0.0054686 ); + IndicatorCylinder3D layerOutflow0( outflow0, 2.*converter.getConversionFactorLength() ); + superGeometry.rename( 2,4,1,layerOutflow0 ); + + // Set material number for outflow1 + //IndicatorCircle3D outflow1(0.2388403,0.0900099,0.0343228, -1.5129,5.1039,-2.8431, 0.0058006 ); + IndicatorCircle3D outflow1( 0.2388403,0.0900099,0.0343228, 0.,-1.,0., 0.0058006 ); + IndicatorCylinder3D layerOutflow1( outflow1, 2.*converter.getConversionFactorLength() ); + superGeometry.rename( 2,5,1,layerOutflow1 ); + + // Removes all not needed boundary voxels outside the surface + superGeometry.clean(); + // Removes all not needed boundary voxels inside the surface + superGeometry.innerClean( 3 ); + superGeometry.checkForErrors(); + + superGeometry.print(); + clout << "Prepare Geometry ... OK" << std::endl; +} + +// Set up the geometry of the simulation +void prepareLattice( SuperLattice3D& lattice, + UnitConverter const& converter, Dynamics& bulkDynamics, + sOnLatticeBoundaryCondition3D& bc, + sOffLatticeBoundaryCondition3D& offBc, + STLreader& stlReader, SuperGeometry3D& superGeometry ) { + + OstreamManager clout( std::cout,"prepareLattice" ); + clout << "Prepare Lattice ..." << std::endl; + + const T omega = converter.getLatticeRelaxationFrequency(); + + // material=0 --> do nothing + lattice.defineDynamics( superGeometry,0,&instances::getNoDynamics() ); + + // material=1 --> bulk dynamics + lattice.defineDynamics( superGeometry,1,&bulkDynamics ); + + if ( bouzidiOn ) { + // material=2 --> no dynamics + bouzidi zero velocity + lattice.defineDynamics( superGeometry,2,&instances::getNoDynamics() ); + offBc.addZeroVelocityBoundary( superGeometry,2,stlReader ); + // material=3 --> no dynamics + bouzidi velocity (inflow) + lattice.defineDynamics( superGeometry,3,&instances::getNoDynamics() ); + offBc.addVelocityBoundary( superGeometry,3,stlReader ); + } else { + // material=2 --> bounceBack dynamics + lattice.defineDynamics( superGeometry, 2, &instances::getBounceBack() ); + // material=3 --> bulk dynamics + velocity (inflow) + lattice.defineDynamics( superGeometry,3,&bulkDynamics ); + bc.addVelocityBoundary( superGeometry,3,omega ); + } + + // material=4,5 --> bulk dynamics + pressure (outflow) + lattice.defineDynamics( superGeometry.getMaterialIndicator({4, 5}),&bulkDynamics ); + bc.addPressureBoundary( superGeometry.getMaterialIndicator({4, 5}),omega ); + + // Initial conditions + AnalyticalConst3D rhoF( 1 ); + std::vector velocity( 3,T() ); + AnalyticalConst3D uF( velocity ); + + // Initialize all values of distribution functions to their local equilibrium + lattice.defineRhoU( superGeometry.getMaterialIndicator({1, 3, 4, 5}),rhoF,uF ); + lattice.iniEquilibrium( superGeometry.getMaterialIndicator({1, 3, 4, 5}),rhoF,uF ); + + // Lattice initialize + lattice.initialize(); + + clout << "Prepare Lattice ... OK" << std::endl; +} + +// Generates a slowly increasing sinuidal inflow +void setBoundaryValues( SuperLattice3D& sLattice, + sOffLatticeBoundaryCondition3D& offBc, + UnitConverter const& converter, int iT, + SuperGeometry3D& superGeometry ) { + + // No of time steps for smooth start-up + int iTperiod = converter.getLatticeTime( 0.5 ); + int iTupdate = 50; + + if ( iT%iTupdate == 0 ) { + // Smooth start curve, sinus + SinusStartScale nSinusStartScale( iTperiod,converter.getCharLatticeVelocity() ); + + // Creates and sets the Poiseuille inflow profile using functors + int iTvec[1]= {iT}; + T maxVelocity[1]= {T()}; + nSinusStartScale( maxVelocity,iTvec ); + CirclePoiseuille3D velocity( superGeometry,3,maxVelocity[0] ); + + if ( bouzidiOn ) { + offBc.defineU( superGeometry,3,velocity ); + } else { + sLattice.defineU( superGeometry,3,velocity ); + } + } +} + +// Computes flux at inflow and outflow +void getResults( SuperLattice3D& sLattice, + UnitConverter& converter, int iT, + Dynamics& bulkDynamics, + SuperGeometry3D& superGeometry, Timer& timer, STLreader& stlReader ) { + + OstreamManager clout( std::cout,"getResults" ); + + SuperVTMwriter3D vtmWriter( "aorta3d" ); + SuperLatticePhysVelocity3D velocity( sLattice, converter ); + SuperLatticePhysPressure3D pressure( sLattice, converter ); + vtmWriter.addFunctor( velocity ); + vtmWriter.addFunctor( pressure ); + + const int vtkIter = converter.getLatticeTime( .1 ); + const int statIter = converter.getLatticeTime( .1 ); + + if ( iT==0 ) { + // Writes the geometry, cuboid no. and rank no. as vti file for visualization + SuperLatticeGeometry3D geometry( sLattice, superGeometry ); + SuperLatticeCuboid3D cuboid( sLattice ); + SuperLatticeRank3D rank( sLattice ); + vtmWriter.write( geometry ); + vtmWriter.write( cuboid ); + vtmWriter.write( rank ); + + vtmWriter.createMasterFile(); + } + + // Writes the vtk files + if ( iT%vtkIter==0 ) { + vtmWriter.write( iT ); + + SuperEuklidNorm3D normVel( velocity ); + BlockReduction3D2D planeReduction( normVel, {0,0,1}, 600, BlockDataSyncMode::ReduceOnly ); + // write output as JPEG + heatmap::write(planeReduction, iT); + } + + // Writes output on the console + if ( iT%statIter==0 ) { + // Timer console output + timer.update( iT ); + timer.printStep(); + + // Lattice statistics console output + sLattice.getStatistics().print( iT,converter.getPhysTime( iT ) ); + + // Flux at the inflow and outflow region + std::vector materials = { 1, 3, 4, 5 }; + + IndicatorCircle3D inflow( 0.218125 ,0.249987-2.*converter.getConversionFactorLength() ,0.0234818, 0., -1.,0., 0.0112342+2*converter.getConversionFactorLength() ); + SuperPlaneIntegralFluxVelocity3D vFluxInflow( sLattice, converter, superGeometry, inflow, materials, BlockDataReductionMode::Discrete ); + vFluxInflow.print( "inflow","ml/s" ); + SuperPlaneIntegralFluxPressure3D pFluxInflow( sLattice, converter, superGeometry, inflow, materials, BlockDataReductionMode::Discrete ); + pFluxInflow.print( "inflow","N","mmHg" ); + + IndicatorCircle3D outflow0( 0.2053696,0.0900099+2.*converter.getConversionFactorLength(),0.0346537, 0.,1.,0., 0.0054686+2*converter.getConversionFactorLength() ); + SuperPlaneIntegralFluxVelocity3D vFluxOutflow0( sLattice, converter, superGeometry, outflow0, materials, BlockDataReductionMode::Discrete ); + vFluxOutflow0.print( "outflow0","ml/s" ); + SuperPlaneIntegralFluxPressure3D pFluxOutflow0( sLattice, converter, superGeometry, outflow0, materials, BlockDataReductionMode::Discrete ); + pFluxOutflow0.print( "outflow0","N","mmHg" ); + + IndicatorCircle3D outflow1( 0.2388403,0.0900099+2.*converter.getConversionFactorLength(),0.0343228, 0.,1.,0., 0.0058006+2*converter.getConversionFactorLength() ); + SuperPlaneIntegralFluxVelocity3D vFluxOutflow1( sLattice, converter, superGeometry, outflow1, materials, BlockDataReductionMode::Discrete ); + vFluxOutflow1.print( "outflow1","ml/s" ); + SuperPlaneIntegralFluxPressure3D pFluxOutflow1( sLattice, converter, superGeometry, outflow1, materials, BlockDataReductionMode::Discrete ); + pFluxOutflow1.print( "outflow1","N","mmHg" ); + + if ( bouzidiOn ) { + SuperLatticeYplus3D yPlus( sLattice, converter, superGeometry, stlReader, 3 ); + SuperMax3D yPlusMaxF( yPlus, superGeometry, 1 ); + int input[4]= {}; + T yPlusMax[1]; + yPlusMaxF( yPlusMax,input ); + clout << "yPlusMax=" << yPlusMax[0] << std::endl; + } + } + + if ( sLattice.getStatistics().getMaxU() > 0.3 ) { + clout << "PROBLEM uMax=" << sLattice.getStatistics().getMaxU() << std::endl; + vtmWriter.write( iT ); + std::exit( 0 ); + } +} + +int main( int argc, char* argv[] ) { + + // === 1st Step: Initialization === + olbInit( &argc, &argv ); + singleton::directories().setOutputDir( "./tmp/" ); + OstreamManager clout( std::cout,"main" ); + // display messages from every single mpi process + //clout.setMultiOutput(true); + + UnitConverter converter( + (T) 0.02246/N, // physDeltaX: spacing between two lattice cells in __m__ + (T) 0.02246/(M*N), // physDeltaT: time step in __s__ + (T) 0.02246, // charPhysLength: reference length of simulation geometry + (T) 0.45, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__ + (T) 0.003/1055., // physViscosity: physical kinematic viscosity in __m^2 / s__ + (T) 1055 // physDensity: physical density in __kg / m^3__ + ); + // Prints the converter log as console output + converter.print(); + // Writes the converter log in a file + converter.write("aorta3d"); + + // === 2nd Step: Prepare Geometry === + + // Instantiation of the STLreader class + // file name, voxel size in meter, stl unit in meter, outer voxel no., inner voxel no. + STLreader stlReader( "aorta3d.stl", converter.getConversionFactorLength(), 0.001, 0, true ); + IndicatorLayer3D extendedDomain( stlReader, converter.getConversionFactorLength() ); + + // Instantiation of a cuboidGeometry with weights +#ifdef PARALLEL_MODE_MPI + const int noOfCuboids = std::min( 16*N,2*singleton::mpi().getSize() ); +#else + const int noOfCuboids = 2; +#endif + CuboidGeometry3D cuboidGeometry( extendedDomain, converter.getConversionFactorLength(), noOfCuboids ); + + // Instantiation of a loadBalancer + HeuristicLoadBalancer loadBalancer( cuboidGeometry ); + + // Instantiation of a superGeometry + SuperGeometry3D superGeometry( cuboidGeometry, loadBalancer, 2 ); + + prepareGeometry( converter, extendedDomain, stlReader, superGeometry ); + + // === 3rd Step: Prepare Lattice === + SuperLattice3D sLattice( superGeometry ); + + SmagorinskyBGKdynamics bulkDynamics( converter.getLatticeRelaxationFrequency(), + instances::getBulkMomenta(), 0.1 ); + + // choose between local and non-local boundary condition + sOnLatticeBoundaryCondition3D sBoundaryCondition( sLattice ); + createInterpBoundaryCondition3D( sBoundaryCondition ); + // createLocalBoundaryCondition3D(sBoundaryCondition); + + sOffLatticeBoundaryCondition3D sOffBoundaryCondition( sLattice ); + createBouzidiBoundaryCondition3D ( sOffBoundaryCondition ); + + Timer timer1( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() ); + timer1.start(); + + prepareLattice( sLattice, converter, bulkDynamics, + sBoundaryCondition, sOffBoundaryCondition, + stlReader, superGeometry ); + + timer1.stop(); + timer1.printSummary(); + + // === 4th Step: Main Loop with Timer === + clout << "starting simulation..." << std::endl; + Timer timer( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() ); + timer.start(); + + for ( int iT = 0; iT <= converter.getLatticeTime( maxPhysT ); iT++ ) { + + // === 5th Step: Definition of Initial and Boundary Conditions === + setBoundaryValues( sLattice, sOffBoundaryCondition, converter, iT, superGeometry ); + + // === 6th Step: Collide and Stream Execution === + sLattice.collideAndStream(); + + // === 7th Step: Computation and Output of the Results === + getResults( sLattice, converter, iT, bulkDynamics, superGeometry, timer, stlReader ); + } + + timer.stop(); + timer.printSummary(); +} diff --git a/examples/turbulence/aorta3d/aorta3d.stl b/examples/turbulence/aorta3d/aorta3d.stl new file mode 100644 index 0000000..340fece Binary files /dev/null and b/examples/turbulence/aorta3d/aorta3d.stl differ diff --git a/examples/turbulence/aorta3d/definitions.mk b/examples/turbulence/aorta3d/definitions.mk new file mode 100644 index 0000000..d15305e --- /dev/null +++ b/examples/turbulence/aorta3d/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 +# +# +# 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 := aorta3d.cpp +OUTPUT := aorta3d diff --git a/examples/turbulence/aorta3d/module.mk b/examples/turbulence/aorta3d/module.mk new file mode 100644 index 0000000..1190482 --- /dev/null +++ b/examples/turbulence/aorta3d/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 +# +# +# 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/examples/turbulence/nozzle3d/Makefile b/examples/turbulence/nozzle3d/Makefile new file mode 100644 index 0000000..a953954 --- /dev/null +++ b/examples/turbulence/nozzle3d/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 +# +# +# 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/examples/turbulence/nozzle3d/definitions.mk b/examples/turbulence/nozzle3d/definitions.mk new file mode 100644 index 0000000..840f62d --- /dev/null +++ b/examples/turbulence/nozzle3d/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 +# +# +# 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 := nozzle3d.cpp +OUTPUT := nozzle3d diff --git a/examples/turbulence/nozzle3d/module.mk b/examples/turbulence/nozzle3d/module.mk new file mode 100644 index 0000000..1190482 --- /dev/null +++ b/examples/turbulence/nozzle3d/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 +# +# +# 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/examples/turbulence/nozzle3d/nozzle3d.cpp b/examples/turbulence/nozzle3d/nozzle3d.cpp new file mode 100644 index 0000000..c736f78 --- /dev/null +++ b/examples/turbulence/nozzle3d/nozzle3d.cpp @@ -0,0 +1,382 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2015 Mathias J. Krause, Patrick Nathan + * 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. + */ + +/* nozzle3d.cpp: + * This example examines a turbulent flow in a nozzle injection tube. At the + * main inlet, either a block profile or a power 1/7 profile is imposed as a + * Dirchlet velocity boundary condition, whereas at the outlet a + * Dirichlet pressure condition is set by p=0 (i.e. rho=1). + * + * The example shows the usage of turbulent models. + * + * The results of a simular simulation setup are publish in TODO + */ + +#include "olb3D.h" +#ifndef OLB_PRECOMPILED // Unless precompiled version is used +#include "olb3D.hh" // Include full template code +#endif +#include +#include +#include +#include + +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace olb::util; +using namespace std; + +typedef double T; + +// Choose your turbulent model of choice +//#define RLB +#define Smagorinsky //default +//#define ConsitentStrainSmagorinsky +//#define ShearSmagorinsky +//#define Krause + +#ifdef ShearSmagorinsky +#define DESCRIPTOR ShearSmagorinskyD3Q19Descriptor +#else +#define DESCRIPTOR D3Q19<> +#endif + +// Parameters for the simulation setup +const int N = 3; // resolution of the model, for RLB N>=5, others N>=2, but N>=5 recommended +const int M = 1; // time discretization refinement +const int inflowProfileMode = 0; // block profile (mode=0), power profile (mode=1) +const T maxPhysT = 200.; // max. simulation time in s, SI unit + +template +class TurbulentVelocity3D : public AnalyticalF3D { + +protected: + // block profile (mode=0), power profile (mode=1) + int _mode; + T rho; + T nu; + T u0; + T p0; + T charL; + T dx; + +public: + TurbulentVelocity3D( UnitConverter const& converter, int mode=0 ) : AnalyticalF3D( 3 ) { + _mode = mode; + u0 = converter.getCharLatticeVelocity(); + rho = converter.getPhysDensity(); + nu = converter.getPhysViscosity(); + charL = converter.getCharPhysLength(); + p0 = converter.getCharPhysPressure(); + dx = converter.getConversionFactorLength(); + + this->getName() = "turbulentVelocity3d"; + }; + + bool operator()( T output[], const S input[] ) override { + T y = input[1]; + T z = input[2]; + // block profile inititalization + T u_calc = u0; + // power profile inititalization + if ( _mode==1 ) { + T obst_y = 5.5+dx; + T obst_z = 5.5+dx; + T obst_r = 0.5; + + T B = 5.5; + T kappa = 0.4; + T ReTau = 183.6; + + u_calc = u0/7.*( 2.*nu*ReTau/( charL*kappa )*log( fabs( 2.*ReTau/charL*( obst_r - sqrt( pow( y - obst_y, 2. ) + + pow( z - obst_z, 2. ) ) )*1.5*( 1 + sqrt( pow( y - obst_y, 2. ) + + pow( z - obst_z, 2. ) )/obst_r )/( 1 + 2.*pow( sqrt( pow( y - obst_y, 2. ) + + pow( z - obst_z, 2. ) )/obst_r, 2. ) ) ) + B ) ); + } + T a = -1., b = 1.; + T nRandom = rand()/( T )RAND_MAX*( b-a ) + a; + + output[0] = u_calc+ 0.15*u0*nRandom; + output[1] = 0.15*u0*nRandom; + output[2] = 0.15*u0*nRandom; + return true; + }; +}; + + +void prepareGeometry( UnitConverter const& converter, IndicatorF3D& indicator, SuperGeometry3D& superGeometry ) { + + OstreamManager clout( std::cout,"prepareGeometry" ); + clout << "Prepare Geometry ..." << std::endl; + + // Sets material number for fluid and boundary + superGeometry.rename( 0,2,indicator ); + + Vector origin( T(), + 5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(), + 5.5*converter.getCharPhysLength()+converter.getConversionFactorLength() ); + + Vector extend( 4.*converter.getCharPhysLength()+5*converter.getConversionFactorLength(), + 5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(), + 5.5*converter.getCharPhysLength()+converter.getConversionFactorLength() ); + + IndicatorCylinder3D inletCylinder( extend, origin, converter.getCharPhysLength() ); + superGeometry.rename( 2,1,inletCylinder ); + + + origin[0]=4.*converter.getCharPhysLength(); + origin[1]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + origin[2]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + + extend[0]=54.*converter.getCharPhysLength(); + extend[1]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + extend[2]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + + IndicatorCylinder3D injectionTube( extend, origin, 5.5*converter.getCharPhysLength() ); + superGeometry.rename( 2,1,injectionTube ); + + origin[0]=converter.getConversionFactorLength(); + origin[1]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + origin[2]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + + extend[0]=T(); + extend[1]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + extend[2]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + + IndicatorCylinder3D cylinderIN( extend, origin, converter.getCharPhysLength() ); + superGeometry.rename( 1,3,cylinderIN ); + + + origin[0]=54.*converter.getCharPhysLength()-converter.getConversionFactorLength(); + origin[1]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + origin[2]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + + extend[0]=54.*converter.getCharPhysLength(); + extend[1]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + extend[2]=5.5*converter.getCharPhysLength()+converter.getConversionFactorLength(); + + IndicatorCylinder3D cylinderOUT( extend, origin, 5.5*converter.getCharPhysLength() ); + superGeometry.rename( 1,4,cylinderOUT ); + + // 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 prepareLattice( SuperLattice3D& sLattice, + UnitConverter const& converter, + Dynamics& bulkDynamics, + sOnLatticeBoundaryCondition3D& bc, + sOffLatticeBoundaryCondition3D& offBc, + SuperGeometry3D& superGeometry ) { + + + + OstreamManager clout( std::cout,"prepareLattice" ); + clout << "Prepare Lattice ..." << std::endl; + + const T omega = converter.getLatticeRelaxationFrequency(); + + // Material=0 -->do nothing + sLattice.defineDynamics( superGeometry, 0, &instances::getNoDynamics() ); + + // Material=1 -->bulk dynamics + // Material=3 -->bulk dynamics (inflow) + // Material=4 -->bulk dynamics (outflow) + sLattice.defineDynamics( superGeometry.getMaterialIndicator({1, 3, 4}), &bulkDynamics ); + + // Material=2 -->bounce back + sLattice.defineDynamics( superGeometry, 2, &instances::getBounceBack() ); + + bc.addVelocityBoundary( superGeometry, 3, omega ); + bc.addPressureBoundary( superGeometry, 4, omega ); + + clout << "Prepare Lattice ... OK" << std::endl; +} + +void setBoundaryValues( UnitConverter const&converter, + SuperLattice3D& lattice, SuperGeometry3D& superGeometry, int iT ) { + + OstreamManager clout( std::cout,"setBoundaryValues" ); + + if ( iT==0 ) { + AnalyticalConst3D rhoF( 1 ); + std::vector velocity( 3,T() ); + AnalyticalConst3D uF( velocity ); + + // Seeding of random fluctuations and definition of the velocity field + srand( time( nullptr ) ); + TurbulentVelocity3D uSol( converter, inflowProfileMode ); + + lattice.iniEquilibrium( superGeometry.getMaterialIndicator({1, 2, 4}), rhoF, uF ); + lattice.iniEquilibrium( superGeometry, 3, rhoF, uSol ); + + lattice.defineU( superGeometry, 3, uSol ); + lattice.defineRho( superGeometry, 4, rhoF ); + + // Make the lattice ready for simulation + lattice.initialize(); + } +} + +void getResults( SuperLattice3D& sLattice, + UnitConverter const& converter, int iT, + SuperGeometry3D& superGeometry, Timer& timer ) { + + OstreamManager clout( std::cout,"getResults" ); + SuperVTMwriter3D vtmWriter( "nozzle3d" ); + + if ( iT==0 ) { + // Writes the geometry, cuboid no. and rank no. as vti file for visualization + SuperLatticeGeometry3D geometry( sLattice, superGeometry ); + SuperLatticeCuboid3D cuboid( sLattice ); + SuperLatticeRank3D rank( sLattice ); + vtmWriter.write( geometry ); + vtmWriter.write( cuboid ); + vtmWriter.write( rank ); + vtmWriter.createMasterFile(); + } + + // Writes the vtk files + if ( iT%converter.getLatticeTime( maxPhysT/100. )==0 ) { + // Create the data-reading functors... + SuperLatticePhysVelocity3D velocity( sLattice, converter ); + SuperLatticePhysPressure3D pressure( sLattice, converter ); + vtmWriter.addFunctor( velocity ); + vtmWriter.addFunctor( pressure ); + vtmWriter.write( iT ); + + SuperEuklidNorm3D normVel( velocity ); + BlockReduction3D2D planeReduction( normVel, {0, 1, 0} ); + // write output as JPEG + heatmap::write(planeReduction, iT); + } + + // Writes output on the console + if ( iT%converter.getLatticeTime( maxPhysT/200. )==0 ) { + timer.update( iT ); + timer.printStep(); + sLattice.getStatistics().print( iT, converter.getPhysTime( iT ) ); + } +} + + + +int main( int argc, char* argv[] ) { + + // === 1st Step: Initialization === + + olbInit( &argc, &argv ); + singleton::directories().setOutputDir( "./tmp/" ); + OstreamManager clout( std::cout,"main" ); + // display messages from every single mpi process + // clout.setMultiOutput(true); + + UnitConverterFromResolutionAndRelaxationTime const converter( + int {N}, // resolution: number of voxels per charPhysL + (T) 0.500018, // latticeRelaxationTime: relaxation time, have 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) 0.0002, // physViscosity: physical kinematic viscosity in __m^2 / s__ + (T) 1.0 // physDensity: physical density in __kg / m^3__ + ); + // Prints the converter log as console output + converter.print(); + // Writes the converter log in a file + converter.write("nozzle3d"); + + Vector origin; + Vector extend( 54.*converter.getCharPhysLength(), 11.*converter.getCharPhysLength()+2.*converter.getConversionFactorLength(), 11.*converter.getCharPhysLength()+2.*converter.getConversionFactorLength() ); + + IndicatorCuboid3D cuboid( extend,origin ); + + CuboidGeometry3D cuboidGeometry( cuboid, converter.getConversionFactorLength(), singleton::mpi().getSize() ); + HeuristicLoadBalancer loadBalancer( cuboidGeometry ); + + // === 2nd Step: Prepare Geometry === + + SuperGeometry3D superGeometry( cuboidGeometry, loadBalancer, 2 ); + prepareGeometry( converter, cuboid, superGeometry ); + + // === 3rd Step: Prepare Lattice === + + SuperLattice3D sLattice( superGeometry ); + + Dynamics* bulkDynamics; + const T omega = converter.getLatticeRelaxationFrequency(); +#if defined(RLB) + bulkDynamics = new RLBdynamics( omega, instances::getBulkMomenta() ); +#elif defined(Smagorinsky) + bulkDynamics = new SmagorinskyBGKdynamics( omega, instances::getBulkMomenta(), + 0.15); +#elif defined(ShearSmagorinsky) + bulkDynamics = new ShearSmagorinskyBGKdynamics( omega, instances::getBulkMomenta(), + 0.15); +#elif defined(Krause) + bulkDynamics = new KrauseBGKdynamics( omega, instances::getBulkMomenta(), + 0.15); +#else //ConsitentStrainSmagorinsky + bulkDynamics = new ConStrainSmagorinskyBGKdynamics( omega, instances::getBulkMomenta(), + 0.05); +#endif + + sOnLatticeBoundaryCondition3D sBoundaryCondition( sLattice ); + createInterpBoundaryCondition3D ( sBoundaryCondition ); + + sOffLatticeBoundaryCondition3D sOffBoundaryCondition( sLattice ); + createBouzidiBoundaryCondition3D ( sOffBoundaryCondition ); + + prepareLattice( sLattice, converter, *bulkDynamics, sBoundaryCondition, sOffBoundaryCondition, superGeometry ); + + // === 4th Step: Main Loop with Timer === + + Timer timer( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() ); + timer.start(); + + for ( int iT = 0; iT <= converter.getLatticeTime( maxPhysT ); ++iT ) { + // === 5ath Step: Apply filter +#ifdef ADM + SuperLatticeADM3D admF( sLattice, 0.01, 2 ); + admF.execute( superGeometry, 1 ); +#endif + // === 5bth Step: Definition of Initial and Boundary Conditions === + setBoundaryValues( converter, sLattice, superGeometry, iT ); + + // === 6th Step: Collide and Stream Execution === + sLattice.collideAndStream(); + + // === 7th Step: Computation and Output of the Results === + getResults( sLattice, converter, iT, superGeometry, timer ); + } + timer.stop(); + timer.printSummary(); + delete bulkDynamics; +} diff --git a/examples/turbulence/tgv3d/Makefile b/examples/turbulence/tgv3d/Makefile new file mode 100644 index 0000000..a953954 --- /dev/null +++ b/examples/turbulence/tgv3d/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 +# +# +# 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/examples/turbulence/tgv3d/Re1600_Brachet.inp b/examples/turbulence/tgv3d/Re1600_Brachet.inp new file mode 100644 index 0000000..4d3eae4 --- /dev/null +++ b/examples/turbulence/tgv3d/Re1600_Brachet.inp @@ -0,0 +1,73 @@ +0.0140725544777 0.000476232358793 +0.285634821966 0.000477040579827 +0.557197089453 0.000477848800861 +0.796810854884 0.000478561937068 +1.02033151325 0.000508989081885 +1.24385217162 0.000539416226701 +1.46737282999 0.000569843371518 +1.69089348835 0.000600270516334 +1.93026954171 0.000660506754588 +2.15367134405 0.000720695450429 +2.34512464432 0.000780789061441 +2.53645908856 0.000870644223478 +2.77559742986 0.000990403563779 +2.9510764791 0.00105044963238 +3.12631781628 0.00117001880302 +3.26972950744 0.00125973133782 +3.44485198858 0.00140906205949 +3.63582986472 0.0015882018746 +3.81071463379 0.00179705569831 +3.96962515184 0.00200586197962 +4.12841681385 0.00224442981195 +4.28685190776 0.00257228229734 +4.44528700167 0.00290013478274 +4.55615591059 0.00313855998783 +4.66690596348 0.00340674674394 +4.82534105739 0.00373459922934 +4.93597225425 0.00403254753647 +5.09464506022 0.00430087691982 +5.28514751222 0.00459906293902 +5.42820263528 0.00477806012689 +5.60296854832 0.00501667550163 +5.7619979224 0.00519572023191 +5.90493418941 0.0054044789708 +6.0477516004 0.00564299926072 +6.19045015534 0.00591128110165 +6.30120020824 0.00617946785776 +6.41195026113 0.00644765461387 +6.52258145798 0.00674560292101 +6.61759497191 0.00695421903265 +6.74467584394 0.00713316867811 +6.87199442803 0.00725259522151 +6.99919415609 0.00740178331594 +7.11018192105 0.00761044697 +7.18922118395 0.00781901553924 +7.26826044685 0.00802758410847 +7.36291739268 0.00832548487319 +7.50490281142 0.00877233602027 +7.59955975725 0.00907023678499 +7.69397899101 0.00942766065176 +7.78875479287 0.00969579986545 +7.88329288266 0.0100234621812 +7.97794982849 0.0103213629459 +8.05627595519 0.0107085008213 +8.15069518895 0.0110659246881 +8.24535213477 0.0113638254528 +8.34012793664 0.0116319646665 +8.43478488246 0.0119298654312 +8.51334872123 0.0122574802045 +8.59203141603 0.0125553334268 +8.68692607392 0.0127937110895 +8.76608419286 0.0129725181077 +8.84571773594 0.0130322789218 +8.94144438607 0.0130623257273 +9.05409613551 0.0128543276671 +9.10273202481 0.0126759009882 +9.16746102116 0.0124677603007 +9.23207116148 0.0122893811642 +9.3128932649 0.0120515264681 +9.40945190728 0.0118732424164 +9.52186594466 0.0117247674582 +9.66599077202 0.0116359106868 +9.79402249232 0.0115767679241 +9.95364614657 0.0116070048992 diff --git a/examples/turbulence/tgv3d/Re3000_Brachet.inp b/examples/turbulence/tgv3d/Re3000_Brachet.inp new file mode 100644 index 0000000..7d90c96 --- /dev/null +++ b/examples/turbulence/tgv3d/Re3000_Brachet.inp @@ -0,0 +1,69 @@ +0.000000000000 0.000268206645697 +0.512243645318 0.000239976881208 +0.880830747993 0.000241075501484 +1.377478762 0.000302168342063 +1.72996868932 0.000333025241995 +2.14656072146 0.000364073206323 +2.62703991096 0.000484730372313 +3.07539639891 0.00063509805186 +3.39562032639 0.00075527755686 +3.69974707851 0.000905215341516 +3.98763335698 0.00114452349734 +4.19524870614 0.00144320491418 +4.43484345855 0.00177178790897 +4.62643328151 0.00207042155971 +4.80192592911 0.0023988134901 +4.97770517331 0.00260798123748 +5.15334111921 0.00287676107636 +5.3291203634 0.00308592882374 +5.44101245023 0.00320548736945 +5.60076616823 0.00341460735072 +5.76059153537 0.00359392128625 +5.87241197305 0.00374328587771 +6.00011463864 0.00395231032679 +6.09576625182 0.00416123924367 +6.17546398794 0.0043403143487 +6.25501842577 0.00457900154523 +6.41398400313 0.00511598802982 +6.54125677383 0.00550384875342 +6.62081121166 0.00574253594996 +6.74837057895 0.00601117249055 +6.87607324454 0.00622019693963 +6.95569933151 0.00642907809041 +7.08304375136 0.00678713276826 +7.14657266298 0.0070257721987 +7.22612710081 0.00726445939523 +7.30560988949 0.00753295263752 +7.48067264221 0.00804018084245 +7.6079454129 0.00842804156606 +7.734931587 0.00893512647269 +7.79831720033 0.00923337799464 +7.87751339242 0.00962109541995 +7.95670958451 0.0100088128453 +8.08340916202 0.0106351219349 +8.1465798279 0.0110227915941 +8.22584766913 0.0113807029737 +8.30461396633 0.0119472566735 +8.36749803562 0.0124541505158 +8.41435657871 0.0129609965919 +8.49348112164 0.013378520063 +8.57217576969 0.0139748798086 +8.66689594394 0.0145712873203 +8.73028155727 0.0148695388422 +8.77778494269 0.0151081305065 +8.80947774936 0.0152572562675 +8.87329325758 0.0153766715149 +8.93739536239 0.0153768625793 +9.00171241465 0.0152876355065 +9.0985104142 0.0150196676912 +9.1796411333 0.0146026218812 +9.2444597296 0.014304752488 +9.29332444884 0.013977029283 +9.35814304514 0.0136791598899 +9.43891551849 0.0134111443085 +9.51975964099 0.0131133226815 +9.61670093884 0.0127857427748 +9.7134989384 0.0125177749596 +9.81022528881 0.0122796131901 +9.92283386712 0.0121011112783 +10.0033197439 0.01195231988 diff --git a/examples/turbulence/tgv3d/Re800_Brachet.inp b/examples/turbulence/tgv3d/Re800_Brachet.inp new file mode 100644 index 0000000..d663d97 --- /dev/null +++ b/examples/turbulence/tgv3d/Re800_Brachet.inp @@ -0,0 +1,64 @@ +0.0130684040331 0.000922657941679 +0.285501702708 0.000923468755068 +0.557839611573 0.00095404118932 +0.862323886562 0.00095494739252 +1.13447101581 0.0010450430685 +1.43866912136 0.00113523413429 +1.71081625061 0.00122532981027 +2.0308490647 0.00137509181269 +2.22277336335 0.00149471063501 +2.49463432316 0.00167409117358 +2.7184188185 0.00185332862744 +2.95803802238 0.00209213701792 +3.18172712791 0.00230113609264 +3.32538418246 0.00248013507197 +3.54878711856 0.00277841900928 +3.70827888165 0.00301698892524 +3.83562427861 0.00328522507226 +3.96306506539 0.00352369959841 +4.09041046235 0.00379193574542 +4.21746968989 0.00414945675503 +4.34462430723 0.00447721614377 +4.45575343642 0.00480492783761 +4.53502236891 0.00507302089991 +4.69384640333 0.00551992216191 +4.80497553251 0.00584763385575 +4.93203476005 0.00620515486536 +5.05918937739 0.0065329142541 +5.2343250694 0.00689057834842 +5.32981026967 0.0070991958639 +5.52106683964 0.00742714603226 +5.66472389419 0.00760614501159 +5.80828555894 0.00781490561178 +6.14338996308 0.00826233151775 +6.38272299753 0.00859042477083 +6.59019583528 0.00885889939237 +6.76561769672 0.0091272786241 +6.94132572759 0.00930637299324 +7.10091288049 0.00951518128833 +7.24476071466 0.00963465702594 +7.40444325737 0.00981370370017 +7.50011923726 0.00996279797392 +7.59569982734 0.0101416538685 +7.70711512596 0.0103800806998 +7.80240954661 0.010648221457 +7.86603455019 0.0107972203409 +7.96132897084 0.0110653610981 +8.05681417111 0.0112739786136 +8.16822946973 0.0115124054449 +8.31179113447 0.011721166045 +8.43951809068 0.0118703557086 +8.58346131466 0.0119600698253 +8.67942346399 0.0120198792365 +8.80753197943 0.0120500224166 +8.9676914712 0.0120802609865 +9.09618154589 0.0119913576832 +9.19262064426 0.01190235899 +9.27312964429 0.0117835509811 +9.33780393578 0.0116051720355 +9.40238283747 0.0114565547108 +9.4829872273 0.011307985081 +9.56378239676 0.0110998922095 +9.7090610781 0.0107729436342 +9.82200261368 0.0105351845316 +9.96689973577 0.0103272824397 diff --git a/examples/turbulence/tgv3d/definitions.mk b/examples/turbulence/tgv3d/definitions.mk new file mode 100644 index 0000000..0117b56 --- /dev/null +++ b/examples/turbulence/tgv3d/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 +# +# +# 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 := tgv3d.cpp +OUTPUT := tgv3d diff --git a/examples/turbulence/tgv3d/module.mk b/examples/turbulence/tgv3d/module.mk new file mode 100644 index 0000000..1190482 --- /dev/null +++ b/examples/turbulence/tgv3d/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 +# +# +# 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/examples/turbulence/tgv3d/tgv3d.cpp b/examples/turbulence/tgv3d/tgv3d.cpp new file mode 100644 index 0000000..5a1fda6 --- /dev/null +++ b/examples/turbulence/tgv3d/tgv3d.cpp @@ -0,0 +1,421 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2017 Mathias J. Krause, Patrick Nathan, Alejandro C. Barreto, Marc Haußmann + * 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. + */ + +/* tgv3d.cpp: + * The Taylor-Green-Vortex (TGV) is one of the simplest configuration, + * where you can investigate the generation of small structures and the resulting turbulence. + * The 2pi periodic box domain and the single mode initial conditions contribute to the simplicity. + * In consequence, the TGV is a common benchmark case for + * Direct Numerical Simulations (DNS) and Large Eddy Simulations (LES). + * + * This example shows the usage and the effects of different subgrid scale turbulence models. + * The molecular dissipation rate, the eddy dissipation rate and + * the effective dissipation rate are calculated and plotted over the simulation time. + * This results can be compared with a published DNS solution, e.g. + * Brachet, Marc E., et al. "Small-scale structure of the Taylor–Green vortex." + * Journal of Fluid Mechanics 130 (1983): 411-452. + */ + + +#include "olb3D.h" +#ifndef OLB_PRECOMPILED // Unless precompiled version is used, +#include "olb3D.hh" // include full template code +#endif +#include +#include +#include +#include + +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace olb::util; +using namespace std; + +typedef double T; + +// Choose your turbulent model of choice +//#define RLB +#define Smagorinsky +//#define WALE +//#define ConsistentStrainSmagorinsky +//#define ShearSmagorinsky +//#define Krause +//#define DNS + +#define finiteDiff //for N<256 + +#ifdef ShearSmagorinsky +#define DESCRIPTOR ShearSmagorinskyD3Q19Descriptor +#elif defined (WALE) +#define DESCRIPTOR WALED3Q19Descriptor +#else +#define DESCRIPTOR D3Q19<> +#endif + +// Global constants +const T pi = 4.0 * std::atan(1.0); +const T volume = pow(2. * pi, 3.); // volume of the 2pi periodic box + + +// Parameters for the simulation setup +const T maxPhysT = 10; // max. simulation time in s, SI unit +int N = 128; // resolution of the model +T Re = 800; // defined as 1/kinematic viscosity +T smagoConst = 0.1; // Smagorisky Constant, for ConsistentStrainSmagorinsky smagoConst = 0.033 +T vtkSave = 0.25; // time interval in s for vtk output +T gnuplotSave = 0.1; // time interval in s for gnuplot output + +bool plotDNS = true; //available for Re=800, Re=1600, Re=3000 (maxPhysT<=10) +vector> values_DNS; + +template +class Tgv3D : public AnalyticalF3D { + +protected: + T u0; + +// initial solution of the TGV +public: + Tgv3D(UnitConverter const& converter, T frac) : AnalyticalF3D(3) + { + u0 = converter.getCharLatticeVelocity(); + }; + + bool operator()(T output[], const S input[]) override + { + T x = input[0]; + T y = input[1]; + T z = input[2]; + + output[0] = u0 * sin(x) * cos(y) * cos(z); + output[1] = -u0 * cos(x) * sin(y) * cos(z); + output[2] = 0; + + return true; + }; +}; + +void prepareGeometry(SuperGeometry3D& superGeometry) +{ + OstreamManager clout(std::cout,"prepareGeometry"); + clout << "Prepare Geometry ..." << std::endl; + + superGeometry.rename(0,1); + superGeometry.communicate(); + + /// 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; + return; +} + +// Set up the geometry of the simulation +void prepareLattice(SuperLattice3D& sLattice, + UnitConverter const& converter, + Dynamics& bulkDynamics, + SuperGeometry3D& superGeometry) +{ + + OstreamManager clout(std::cout,"prepareLattice"); + clout << "Prepare Lattice ..." << std::endl; + + /// Material=0 -->do nothing + sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics()); + /// Material=1 -->bulk dynamics + sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); + + sLattice.initialize(); + clout << "Prepare Lattice ... OK" << std::endl; +} + +void setBoundaryValues(SuperLattice3D& sLattice, + UnitConverter const& converter, + SuperGeometry3D& superGeometry) +{ + OstreamManager clout(std::cout,"setBoundaryValues"); + + Analytic