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/thermal/porousPlate2d/Makefile | 113 ++++ examples/thermal/porousPlate2d/definitions.mk | 30 + examples/thermal/porousPlate2d/module.mk | 29 + .../thermal/porousPlate2d/thermalPorousPlate2d.cpp | 524 +++++++++++++++ examples/thermal/porousPlate3d/Makefile | 113 ++++ examples/thermal/porousPlate3d/definitions.mk | 30 + examples/thermal/porousPlate3d/module.mk | 29 + .../thermal/porousPlate3d/thermalPorousPlate3d.cpp | 494 ++++++++++++++ examples/thermal/rayleighBenard2d/Makefile | 113 ++++ examples/thermal/rayleighBenard2d/definitions.mk | 30 + examples/thermal/rayleighBenard2d/module.mk | 29 + .../thermal/rayleighBenard2d/rayleighBenard2d.cpp | 346 ++++++++++ examples/thermal/rayleighBenard3d/Makefile | 113 ++++ examples/thermal/rayleighBenard3d/definitions.mk | 30 + examples/thermal/rayleighBenard3d/module.mk | 29 + .../thermal/rayleighBenard3d/rayleighBenard3d.cpp | 336 ++++++++++ examples/thermal/squareCavity2d/Makefile | 113 ++++ examples/thermal/squareCavity2d/definitions.mk | 30 + examples/thermal/squareCavity2d/module.mk | 29 + examples/thermal/squareCavity2d/squareCavity2d.cpp | 708 +++++++++++++++++++++ examples/thermal/squareCavity3d/Makefile | 113 ++++ examples/thermal/squareCavity3d/definitions.mk | 30 + examples/thermal/squareCavity3d/module.mk | 29 + examples/thermal/squareCavity3d/squareCavity3d.cpp | 518 +++++++++++++++ 24 files changed, 3958 insertions(+) create mode 100644 examples/thermal/porousPlate2d/Makefile create mode 100644 examples/thermal/porousPlate2d/definitions.mk create mode 100644 examples/thermal/porousPlate2d/module.mk create mode 100644 examples/thermal/porousPlate2d/thermalPorousPlate2d.cpp create mode 100644 examples/thermal/porousPlate3d/Makefile create mode 100644 examples/thermal/porousPlate3d/definitions.mk create mode 100644 examples/thermal/porousPlate3d/module.mk create mode 100644 examples/thermal/porousPlate3d/thermalPorousPlate3d.cpp create mode 100644 examples/thermal/rayleighBenard2d/Makefile create mode 100644 examples/thermal/rayleighBenard2d/definitions.mk create mode 100644 examples/thermal/rayleighBenard2d/module.mk create mode 100644 examples/thermal/rayleighBenard2d/rayleighBenard2d.cpp create mode 100644 examples/thermal/rayleighBenard3d/Makefile create mode 100644 examples/thermal/rayleighBenard3d/definitions.mk create mode 100644 examples/thermal/rayleighBenard3d/module.mk create mode 100644 examples/thermal/rayleighBenard3d/rayleighBenard3d.cpp create mode 100644 examples/thermal/squareCavity2d/Makefile create mode 100644 examples/thermal/squareCavity2d/definitions.mk create mode 100644 examples/thermal/squareCavity2d/module.mk create mode 100644 examples/thermal/squareCavity2d/squareCavity2d.cpp create mode 100644 examples/thermal/squareCavity3d/Makefile create mode 100644 examples/thermal/squareCavity3d/definitions.mk create mode 100644 examples/thermal/squareCavity3d/module.mk create mode 100644 examples/thermal/squareCavity3d/squareCavity3d.cpp (limited to 'examples/thermal') diff --git a/examples/thermal/porousPlate2d/Makefile b/examples/thermal/porousPlate2d/Makefile new file mode 100644 index 0000000..345217e --- /dev/null +++ b/examples/thermal/porousPlate2d/Makefile @@ -0,0 +1,113 @@ +# 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 := +OUTPUT := thermalPorousPlate2d + +########################################################################### +## definitions + +include $(ROOT)/global.mk + +OBJECTS := $(foreach file, $(SRC) $(OUTPUT), $(PWD)/$(file).o) +DEPS := $(foreach file, $(SRC) $(OUTPUT), $(PWD)/$(file).d) + +########################################################################### +## all + +all : depend compile updatelib 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: + @echo Clean olb main + @cd $(ROOT); \ + $(MAKE) cleanbuild; + +########################################################################### +## update lib + +updatelib : + @cd $(ROOT); \ + $(MAKE) all; + +########################################################################### +## link + +link: $(OUTPUT) + +$(OUTPUT): $(OBJECTS) $(ROOT)/$(LIBDIR)/lib$(LIB).a + @echo Link $@ + $(CXX) $(foreach file, $(SRC), $(file).o) $@.o $(LDFLAGS) -L$(ROOT)/$(LIBDIR) $(LIBS) -o $@ + +########################################################################### +## include dependencies + +ifneq "$(strip $(wildcard *.d))" "" + include $(foreach file,$(DEPS),$(file)) +endif + +########################################################################### +########################################################################### diff --git a/examples/thermal/porousPlate2d/definitions.mk b/examples/thermal/porousPlate2d/definitions.mk new file mode 100644 index 0000000..d5858ad --- /dev/null +++ b/examples/thermal/porousPlate2d/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 := thermalPorousPlate2d.cpp +OUTPUT := thermalPorousPlate2d diff --git a/examples/thermal/porousPlate2d/module.mk b/examples/thermal/porousPlate2d/module.mk new file mode 100644 index 0000000..1190482 --- /dev/null +++ b/examples/thermal/porousPlate2d/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/thermal/porousPlate2d/thermalPorousPlate2d.cpp b/examples/thermal/porousPlate2d/thermalPorousPlate2d.cpp new file mode 100644 index 0000000..e7d859e --- /dev/null +++ b/examples/thermal/porousPlate2d/thermalPorousPlate2d.cpp @@ -0,0 +1,524 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2008 Orestis Malaspinas + * 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. + */ + +/* rayleighBenard2d.cpp: + * Rayleigh-Benard convection rolls in 2D, simulated with + * the thermal LB model by Z. Guo e.a., between a hot plate at + * the bottom and a cold plate at the top. + */ + + +#include "olb2D.h" +#include "olb2D.hh" // use only generic version! +#include +#include +#include +#include +#include + + +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace std; + +typedef double T; + +#define NSDESCRIPTOR D2Q9 +#define TDESCRIPTOR D2Q5 + +// #define TemperatureBoundary +// #define RegularizedTemperatureBoundary +#define RegularizedHeatFluxBoundary + +// const int maxIter = 1000000; +// const int saveIter = 5000; + +// Parameters for the simulation setup +const T lx = 1.; // length of the channel +const T ly = 1.; // height of the channel +int N = 20; // resolution of the model +T tau = 1.; // relaxation time +const T Re = 20.; // Reynolds number +const T Ra = 100.; // Rayleigh number +const T Pr = 0.71; // Prandtl number +const T maxPhysT = 1e4; // max. simulation time in s, SI unit +const T epsilon = 1.e-7; // precision of the convergence (residuum) + +const T Tcold = 273.15; +const T Thot = 274.15; + +int iT; + +// analytical solution from point light source in infinte domain +// appliacation from R3 to R1. +// effective for x in R3, only the distance to (0,0) is needed. +// documentation e.g. Biomedical Optics, Lihong V. Wang Hsin-I Wu +template +class AnalyticalVelocityPorousPlate2D : public AnalyticalF2D { +private: + T _Re; + T _u0; + T _v0; + T _ly; +public: + AnalyticalVelocityPorousPlate2D(T Re, T u0, T v0, T ly) : AnalyticalF2D(2), + _Re(Re), _u0(u0), _v0(v0), _ly(ly) + { + this->getName() = "AnalyticalVelocityPorousPlate2D"; + }; + + bool operator()(T output[2], const S x[2]) + { + output[0] = _u0*((exp(_Re* x[1] / _ly) - 1) / (exp(_Re) - 1)); + output[1] = _v0; + return true; + }; +}; + +template +class AnalyticalTemperaturePorousPlate2D : public AnalyticalF2D { +private: + T _Re; + T _Pr; + T _ly; + T _T0; + T _deltaT; +public: + AnalyticalTemperaturePorousPlate2D(T Re, T Pr, T ly, T T0, T deltaT) : AnalyticalF2D(1), + _Re(Re), _Pr(Pr), _ly(ly), _T0(T0), _deltaT(deltaT) + { + this->getName() = "AnalyticalTemperaturePorousPlate2D"; + }; + + bool operator()(T output[1], const S x[2]) + { + output[0] = _T0 + _deltaT*((exp(_Pr*_Re*x[1] / _ly) - 1) / (exp(_Pr*_Re) - 1)); + return true; + }; +}; + +template +class AnalyticalHeatFluxPorousPlate2D : public AnalyticalF2D { +private: + T _Re; + T _Pr; + T _deltaT; + T _ly; + T _lambda; +public: + AnalyticalHeatFluxPorousPlate2D(T Re, T Pr, T deltaT, T ly,T lambda) : AnalyticalF2D(2), + _Re(Re), _Pr(Pr), _deltaT(deltaT), _ly(ly), _lambda(lambda) + { + this->getName() = "AnalyticalHeatFluxPorousPlate2D"; + }; + + bool operator()(T output[2], const S x[2]) + { + output[0] = 0; + output[1] = - _lambda * _Re * _Pr * _deltaT / _ly * (exp(_Pr * _Re * x[1] / _ly))/(exp(_Pr * _Re) - 1); + return true; + }; +}; + +void error( SuperGeometry2D& superGeometry, + SuperLattice2D& NSlattice, + SuperLattice2D& ADlattice, + ThermalUnitConverter const& converter, + T Re ) +{ + OstreamManager clout( std::cout, "error" ); + + int input[1] = { }; + T result[1] = { }; + + auto indicatorF = superGeometry.getMaterialIndicator({1, 2, 3}); + + T u_Re = Re * converter.getPhysViscosity() / converter.getCharPhysLength(); + AnalyticalVelocityPorousPlate2D uSol(Re, converter.getCharPhysVelocity(), u_Re, converter.getCharPhysLength()); + SuperLatticePhysVelocity2D u(NSlattice,converter); + + SuperAbsoluteErrorL2Norm2D absVelocityErrorNormL2(u, uSol, indicatorF); + absVelocityErrorNormL2(result, input); + clout << "velocity-L2-error(abs)=" << result[0]; + SuperRelativeErrorL2Norm2D relVelocityErrorNormL2(u, uSol, indicatorF); + relVelocityErrorNormL2(result, input); + clout << "; velocity-L2-error(rel)=" << result[0] << std::endl; + + AnalyticalTemperaturePorousPlate2D tSol(Re, Pr, converter.getCharPhysLength(), converter.getCharPhysLowTemperature(), converter.getCharPhysTemperatureDifference()); + SuperLatticePhysTemperature2D t(ADlattice,converter); + + SuperAbsoluteErrorL2Norm2D absTemperatureErrorNormL2(t, tSol, indicatorF); + absTemperatureErrorNormL2(result, input); + clout << "temperature-L2-error(abs)=" << result[0]; + SuperRelativeErrorL2Norm2D relTemperatureErrorNormL2(t, tSol, indicatorF); + relTemperatureErrorNormL2(result, input); + clout << "; temperature-L2-error(rel)=" << result[0] << std::endl; + + AnalyticalHeatFluxPorousPlate2D HeatFluxSol(Re, Pr, converter.getCharPhysTemperatureDifference(), converter.getCharPhysLength(), converter.getThermalConductivity()); + SuperLatticePhysHeatFlux2D HeatFlux(ADlattice,converter); + + SuperAbsoluteErrorL2Norm2D absHeatFluxErrorNormL2(HeatFlux, HeatFluxSol, indicatorF); + absHeatFluxErrorNormL2(result, input); + clout << "heatFlux-L2-error(abs)=" << result[0]; + SuperRelativeErrorL2Norm2D relHeatFluxErrorNormL2(HeatFlux, HeatFluxSol, indicatorF); + relHeatFluxErrorNormL2(result, input); + clout << "; heatFlux-L2-error(rel)=" << result[0] << std::endl; +} + +/// Stores geometry information in form of material numbers +void prepareGeometry(SuperGeometry2D& superGeometry, + ThermalUnitConverter const& converter) +{ + + OstreamManager clout(std::cout,"prepareGeometry"); + clout << "Prepare Geometry ..." << std::endl; + + superGeometry.rename(0,2); + superGeometry.rename(2,1,0,1); + + std::vector extend( 2, T(0) ); + extend[0] = lx+2*converter.getPhysLength(1); + extend[1] = converter.getPhysLength(1); + std::vector origin( 2, T(0) ); + origin[0] = -converter.getPhysLength(1); + IndicatorCuboid2D bottom(extend, origin); + /// Set material number for bottom + superGeometry.rename(2,3,1,bottom); + + /// 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( ThermalUnitConverter const& converter, + SuperLattice2D& NSlattice, + SuperLattice2D& ADlattice, + Dynamics &bulkDynamics, + Dynamics& advectionDiffusionBulkDynamics, + sOnLatticeBoundaryCondition2D& NSboundaryCondition, + sOnLatticeBoundaryCondition2D& TboundaryCondition, + SuperGeometry2D& superGeometry ) +{ + + OstreamManager clout(std::cout,"prepareLattice"); + + T Tomega = converter.getLatticeThermalRelaxationFrequency(); + T NSomega = converter.getLatticeRelaxationFrequency(); + + /// define lattice Dynamics + clout << "defining dynamics" << endl; + + ADlattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics()); + NSlattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics()); + + ADlattice.defineDynamics(superGeometry, 1, &advectionDiffusionBulkDynamics); + ADlattice.defineDynamics(superGeometry, 2, &advectionDiffusionBulkDynamics); + ADlattice.defineDynamics(superGeometry, 3, &advectionDiffusionBulkDynamics); + NSlattice.defineDynamics(superGeometry, 1, &bulkDynamics); + NSlattice.defineDynamics(superGeometry, 2, &bulkDynamics); + NSlattice.defineDynamics(superGeometry, 3, &bulkDynamics); + + + /// sets boundary + NSboundaryCondition.addVelocityBoundary(superGeometry, 2, NSomega); + NSboundaryCondition.addVelocityBoundary(superGeometry, 3, NSomega); + + #ifdef TemperatureBoundary + TboundaryCondition.addTemperatureBoundary(superGeometry, 2, Tomega); + TboundaryCondition.addTemperatureBoundary(superGeometry, 3, Tomega); + #endif + #ifdef RegularizedTemperatureBoundary + TboundaryCondition.addRegularizedTemperatureBoundary(superGeometry, 2, Tomega); + TboundaryCondition.addRegularizedTemperatureBoundary(superGeometry, 3, Tomega); + #endif + #ifdef RegularizedHeatFluxBoundary + T heatFlux[2]; + T input[2] = {0.,1.}; + AnalyticalHeatFluxPorousPlate2D HeatFluxSol(Re, Pr, converter.getCharPhysTemperatureDifference(), converter.getCharPhysLength(), converter.getThermalConductivity()); + HeatFluxSol(heatFlux, input); + T temp = converter.getLatticeSpecificHeatCapacity(converter.getPhysSpecificHeatCapacity())*(converter.getLatticeThermalRelaxationTime() - 0.5) / converter.getLatticeThermalRelaxationTime(); + heatFlux[0] = converter.getLatticeHeatFlux(heatFlux[0]) / temp; + heatFlux[1] = converter.getLatticeHeatFlux(heatFlux[1]) / temp; + TboundaryCondition.addRegularizedHeatFluxBoundary(superGeometry, 2, Tomega, heatFlux); + TboundaryCondition.addRegularizedTemperatureBoundary(superGeometry, 3, Tomega); + #endif +} + +void setBoundaryValues(ThermalUnitConverter const& converter, + SuperLattice2D& NSlattice, + SuperLattice2D& ADlattice, + int iT, SuperGeometry2D& superGeometry) +{ + + if (iT == 0) { + + // typedef advectionDiffusionLbHelpers TlbH; + + /// for each material set the defineRhoU and the Equilibrium + std::vector zero(2,T()); + AnalyticalConst2D u(zero); + AnalyticalConst2D rho(1.); + AnalyticalConst2D force(zero); + + T u_Re = converter.getLatticeVelocity( Re * converter.getPhysViscosity() / converter.getCharPhysLength() ); + + AnalyticalConst2D u_top(converter.getCharLatticeVelocity(), u_Re); + AnalyticalConst2D u_bot(0.0, u_Re); + + NSlattice.defineRhoU(superGeometry, 1, rho, u); + NSlattice.iniEquilibrium(superGeometry, 1, rho, u); + NSlattice.defineField(superGeometry, 1, force); + NSlattice.defineRhoU(superGeometry, 2, rho, u_top); + NSlattice.iniEquilibrium(superGeometry, 2, rho, u_top); + NSlattice.defineField(superGeometry, 2, force); + NSlattice.defineRhoU(superGeometry, 3, rho, u_bot); + NSlattice.iniEquilibrium(superGeometry, 3, rho, u_bot); + NSlattice.defineField(superGeometry, 3, force); + + AnalyticalConst2D Cold(converter.getLatticeTemperature(Tcold)); + AnalyticalConst2D Hot(converter.getLatticeTemperature(Thot)); + + ADlattice.defineRho(superGeometry, 1, Cold); + ADlattice.iniEquilibrium(superGeometry, 1, Cold, u); + ADlattice.defineField(superGeometry, 1, u); + ADlattice.defineRhoU(superGeometry, 2, Hot, u_top); + ADlattice.iniEquilibrium(superGeometry, 2, Hot, u_top); + ADlattice.defineField(superGeometry, 2, u); + ADlattice.defineRhoU(superGeometry, 3, Cold, u_bot); + ADlattice.iniEquilibrium(superGeometry, 3, Cold, u_bot); + ADlattice.defineField(superGeometry, 3, u); + + /// Make the lattice ready for simulation + NSlattice.initialize(); + ADlattice.initialize(); + } +} + +void getResults(ThermalUnitConverter const& converter, + SuperLattice2D& NSlattice, + SuperLattice2D& ADlattice, int iT, + SuperGeometry2D& superGeometry, + Timer& timer, + bool converged) +{ + + OstreamManager clout(std::cout,"getResults"); + + SuperVTMwriter2D vtkWriter("thermalPorousPlate2d"); + SuperLatticeGeometry2D geometry2(NSlattice, superGeometry); + SuperLatticePhysVelocity2D velocity(NSlattice, converter); + SuperLatticePhysPressure2D pressure(NSlattice, converter); + SuperLatticePhysTemperature2D temperature(ADlattice, converter); + SuperLatticePhysHeatFlux2D heatflux(ADlattice, converter); + + T u_Re = Re * converter.getPhysViscosity() / converter.getCharPhysLength(); + AnalyticalVelocityPorousPlate2D uSol(Re, converter.getCharPhysVelocity(), u_Re, converter.getCharPhysLength()); + SuperLatticeFfromAnalyticalF2D uSolLattice(uSol,NSlattice); + AnalyticalTemperaturePorousPlate2D TSol(Re, Pr, converter.getCharPhysLength(), converter.getCharPhysLowTemperature(), converter.getCharPhysTemperatureDifference()); + SuperLatticeFfromAnalyticalF2D TSolLattice(TSol,ADlattice); + AnalyticalHeatFluxPorousPlate2D HeatFluxSol(Re, Pr, converter.getCharPhysTemperatureDifference(), converter.getCharPhysLength(), converter.getThermalConductivity()); + SuperLatticeFfromAnalyticalF2D HeatFluxSolLattice(HeatFluxSol,ADlattice); + + vtkWriter.addFunctor( geometry2 ); + vtkWriter.addFunctor( pressure ); + vtkWriter.addFunctor( velocity ); + vtkWriter.addFunctor( temperature ); + vtkWriter.addFunctor( heatflux ); + vtkWriter.addFunctor( uSolLattice ); + vtkWriter.addFunctor( TSolLattice ); + vtkWriter.addFunctor( HeatFluxSolLattice ); + + const int vtkIter = converter.getLatticeTime(10.); + + if (iT == 0) { + /// Writes the converter log file + // writeLogFile(converter,"thermalPorousPlate2d"); + T tmpIn[2] = {0.,1.}; + T tmpOut[2]; + HeatFluxSol(tmpOut,tmpIn); + clout << converter.getLatticeHeatFlux(tmpOut[0]) << " " << converter.getLatticeHeatFlux(tmpOut[1]) << endl; + clout << tmpOut[0] << " " << tmpOut[1] << endl; + + /// Writes the geometry, cuboid no. and rank no. as vti file for visualization + SuperLatticeGeometry2D geometry(NSlattice, superGeometry); + SuperLatticeCuboid2D cuboid(NSlattice); + SuperLatticeRank2D rank(NSlattice); + vtkWriter.write(geometry); + vtkWriter.write(cuboid); + vtkWriter.write(rank); + + vtkWriter.createMasterFile(); + } + + /// Writes the VTK files + if (iT%vtkIter == 0 || converged) { + NSlattice.getStatistics().print(iT,converter.getPhysTime(iT)); + timer.print(iT); + error(superGeometry, NSlattice, ADlattice, converter, Re); + + cout << endl << endl; + vtkWriter.write(iT); + + ///writes Jpeg + //SuperEuklidNorm2D normVel(velocity); + BlockReduction2D2D planeReduction( temperature, 600, BlockDataSyncMode::ReduceOnly ); + // write output of velocity as JPEG + heatmap::plotParam jpeg_Param; + jpeg_Param.maxValue = Thot; + jpeg_Param.minValue = Tcold; + heatmap::write(planeReduction, iT, jpeg_Param); + } + +} + +int main(int argc, char *argv[]) +{ + + /// === 1st Step: Initialization === + OstreamManager clout(std::cout,"main"); + olbInit(&argc, &argv); + singleton::directories().setOutputDir("./tmp/"); + + if (argc >= 2) { + N = atoi(argv[1]); + } + if (argc == 3) { + tau = atof(argv[2]); + } + + ThermalUnitConverter const converter( + (T) 1.0 / N, // physDeltaX + (T) 1.0 / N * 1.0 / 1e-3 * (tau - 0.5) / 3 / N, // physDeltaT + (T) 1.0, // charPhysLength + (T) sqrt( 9.81 * Ra * 1e-3 * 1e-3 / Pr / 9.81 / (Thot - Tcold) / pow(1.0, 3) * (Thot - Tcold) * 1.0 ), // charPhysVelocity + (T) 1e-3, // physViscosity + (T) 1.0, // physDensity + (T) 0.03, // physThermalConductivity + (T) Pr * 0.03 / 1e-3 / 1.0, // physSpecificHeatCapacity + (T) Ra * 1e-3 * 1e-3 / Pr / 9.81 / (Thot - Tcold) / pow(1.0, 3), // physThermalExpansionCoefficient + (T) Tcold, // charPhysLowTemperature + (T) Thot // charPhysHighTemperature + ); + converter.print(); + + /// === 2nd Step: Prepare Geometry === + std::vector extend(2,T()); + extend[0] = lx; + extend[1] = ly; + std::vector origin(2,T()); + IndicatorCuboid2D cuboid(extend, origin); + + /// Instantiation of a cuboidGeometry with weights +#ifdef PARALLEL_MODE_MPI + const int noOfCuboids = singleton::mpi().getSize(); +#else + const int noOfCuboids = 1; +#endif + CuboidGeometry2D cuboidGeometry(cuboid, converter.getPhysDeltaX(), noOfCuboids); + cuboidGeometry.setPeriodicity(true, false); + + /// Instantiation of a loadBalancer + HeuristicLoadBalancer loadBalancer(cuboidGeometry); + + /// Instantiation of a superGeometry + SuperGeometry2D superGeometry(cuboidGeometry, loadBalancer, 2); + + prepareGeometry(superGeometry, converter); + + /// === 3rd Step: Prepare Lattice === + + SuperLattice2D ADlattice(superGeometry); + SuperLattice2D NSlattice(superGeometry); + + sOnLatticeBoundaryCondition2D NSboundaryCondition(NSlattice); + createLocalBoundaryCondition2D(NSboundaryCondition); + + sOnLatticeBoundaryCondition2D TboundaryCondition(ADlattice); + createAdvectionDiffusionBoundaryCondition2D(TboundaryCondition); + + ForcedBGKdynamics NSbulkDynamics( + converter.getLatticeRelaxationFrequency(), + instances::getBulkMomenta()); + + BGKdynamics TbulkDynamics ( + converter.getLatticeThermalRelaxationFrequency(), + instances::getAdvectionDiffusionBulkMomenta() + ); + + // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!// + // This coupling must be necessarily be put on the Navier-Stokes lattice!! + // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!// + + std::vector dir{0.0, 1.0}; + + T boussinesqForcePrefactor = 9.81 / converter.getConversionFactorVelocity() * converter.getConversionFactorTime() * + converter.getCharPhysTemperatureDifference() * converter.getPhysThermalExpansionCoefficient(); + + NavierStokesAdvectionDiffusionCouplingGenerator2D + coupling(0, converter.getLatticeLength(lx), 0, converter.getLatticeLength(ly), + boussinesqForcePrefactor, converter.getLatticeTemperature(Tcold), 1., dir); + + NSlattice.addLatticeCoupling(coupling, ADlattice); + + prepareLattice(converter, + NSlattice, ADlattice, + NSbulkDynamics, TbulkDynamics, + NSboundaryCondition, TboundaryCondition, superGeometry ); + + /// === 4th Step: Main Loop with Timer === + Timer timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel() ); + timer.start(); + + util::ValueTracer converge(converter.getLatticeTime(1.0),epsilon); + for (iT = 0; iT < converter.getLatticeTime(maxPhysT); ++iT) { + + if (converge.hasConverged()) { + clout << "Simulation converged." << endl; + getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged()); + break; + } + + /// === 5th Step: Definition of Initial and Boundary Conditions === + setBoundaryValues(converter, NSlattice, ADlattice, iT, superGeometry); + + /// === 6th Step: Collide and Stream Execution === + NSlattice.collideAndStream(); + NSlattice.executeCoupling(); + ADlattice.collideAndStream(); + + /// === 7th Step: Computation and Output of the Results === + getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged()); + converge.takeValue(NSlattice.getStatistics().getAverageEnergy()); + } + + timer.stop(); + timer.printSummary(); +} diff --git a/examples/thermal/porousPlate3d/Makefile b/examples/thermal/porousPlate3d/Makefile new file mode 100644 index 0000000..e2c7ac8 --- /dev/null +++ b/examples/thermal/porousPlate3d/Makefile @@ -0,0 +1,113 @@ +# 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 := +OUTPUT := thermalPorousPlate3d + +########################################################################### +## definitions + +include $(ROOT)/global.mk + +OBJECTS := $(foreach file, $(SRC) $(OUTPUT), $(PWD)/$(file).o) +DEPS := $(foreach file, $(SRC) $(OUTPUT), $(PWD)/$(file).d) + +########################################################################### +## all + +all : depend compile updatelib 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/*.* line.ghia.txt + +cleanobj: + @echo Clean object files + @rm -f $(OBJECTS) + +cleandep: + @echo Clean dependencies files + @rm -f $(DEPS) + +cleanbuild: + @echo Clean olb main + @cd $(ROOT); \ + $(MAKE) cleanbuild; + +########################################################################### +## update lib + +updatelib : + @cd $(ROOT); \ + $(MAKE) all; + +########################################################################### +## link + +link: $(OUTPUT) + +$(OUTPUT): $(OBJECTS) $(ROOT)/$(LIBDIR)/lib$(LIB).a + @echo Link $@ + $(CXX) $(foreach file, $(SRC), $(file).o) $@.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/thermal/porousPlate3d/definitions.mk b/examples/thermal/porousPlate3d/definitions.mk new file mode 100644 index 0000000..bceea4d --- /dev/null +++ b/examples/thermal/porousPlate3d/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 := thermalPorousPlate3d.cpp +OUTPUT := thermalPorousPlate3d diff --git a/examples/thermal/porousPlate3d/module.mk b/examples/thermal/porousPlate3d/module.mk new file mode 100644 index 0000000..1190482 --- /dev/null +++ b/examples/thermal/porousPlate3d/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/thermal/porousPlate3d/thermalPorousPlate3d.cpp b/examples/thermal/porousPlate3d/thermalPorousPlate3d.cpp new file mode 100644 index 0000000..c197228 --- /dev/null +++ b/examples/thermal/porousPlate3d/thermalPorousPlate3d.cpp @@ -0,0 +1,494 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2008 Orestis Malaspinas + * 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. + */ + +/* rayleighBenard2d.cpp: + * Rayleigh-Benard convection rolls in 2D, simulated with + * the thermal LB model by Z. Guo e.a., between a hot plate at + * the bottom and a cold plate at the top. + */ + +#include "olb3D.h" +#include "olb3D.hh" // use only generic version! +#include +#include +#include +#include +#include + + +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace std; + +typedef double T; + + +#define NSDESCRIPTOR D3Q19 +#define TDESCRIPTOR D3Q7 + + + + +// const int maxIter = 1000000; +// const int saveIter = 5000; + +// Parameters for the simulation setup +const T lx = 1.0; // length of the channel +const T ly = 1.0; // height of the channel +T lz = 1.0; +int N = 20; // resolution of the model +T tau = 1.; // relaxation time +const T Re = 5.; // Reynolds number +const T Ra = 100.; // Rayleigh number +const T Pr = 0.71; // Prandtl number +const T maxPhysT = 1e4; // max. simulation time in s, SI unit +const T epsilon = 1.e-7; // precision of the convergence (residuum) + +const T Tcold = 273.15; +const T Thot = 274.15; + + +template +class AnalyticalVelocityPorousPlate3D : public AnalyticalF3D { +private: + T _Re; + T _u0; + T _v0; + T _ly; +public: + AnalyticalVelocityPorousPlate3D(T Re, T u0, T v0, T ly) : AnalyticalF3D(3), + _Re(Re), _u0(u0), _v0(v0), _ly(ly) + { + this->getName() = "AnalyticalVelocityPorousPlate3D"; + }; + + bool operator()(T output[3], const S x[3]) + { + output[0] = _u0*((exp(_Re* x[1] / _ly) - 1) / (exp(_Re) - 1)); + output[1] = _v0; + output[2] = 0.0; + return true; + }; +}; + +template +class AnalyticalTemperaturePorousPlate3D : public AnalyticalF3D { +private: + T _Re; + T _Pr; + T _ly; + T _T0; + T _deltaT; +public: + AnalyticalTemperaturePorousPlate3D(T Re, T u0, T ly, T T0, T deltaT) : AnalyticalF3D(1), + _Re(Re), _Pr(Pr), _ly(ly), _T0(T0), _deltaT(deltaT) + { + this->getName() = "AnalyticalTemperaturePorousPlate3D"; + }; + + bool operator()(T output[1], const S x[3]) + { + output[0] = _T0 + _deltaT*((exp(_Pr*_Re*x[1] / _ly) - 1) / (exp(_Pr*_Re) - 1)); + return true; + }; +}; + +template +class AnalyticalHeatFluxPorousPlate3D : public AnalyticalF3D { +private: + T _Re; + T _Pr; + T _deltaT; + T _ly; + T _lambda; +public: + AnalyticalHeatFluxPorousPlate3D(T Re, T Pr, T deltaT, T ly,T lambda) : AnalyticalF3D(3), + _Re(Re), _Pr(Pr), _deltaT(deltaT), _ly(ly), _lambda(lambda) + { + this->getName() = "AnalyticalHeatFluxPorousPlate3D"; + }; + + bool operator()(T output[3], const S x[3]) + { + output[0] = 0; + output[1] = - _lambda * _Re * _Pr * _deltaT / _ly * (exp(_Pr * _Re * x[1] / _ly))/(exp(_Pr * _Re) - 1); + output[2] = 0; + return true; + }; +}; + +void error( SuperGeometry3D& superGeometry, + SuperLattice3D& NSlattice, + SuperLattice3D& ADlattice, + ThermalUnitConverter const& converter, + T Re ) +{ + OstreamManager clout( std::cout, "error" ); + + int input[1] = { }; + T result[1] = { }; + + auto indicatorF = superGeometry.getMaterialIndicator(1); + + T u_Re = Re * converter.getPhysViscosity() / converter.getCharPhysLength(); + AnalyticalVelocityPorousPlate3D uSol(Re,converter.getCharPhysVelocity(), u_Re, converter.getCharPhysLength()); + SuperLatticePhysVelocity3D u(NSlattice,converter); + + SuperAbsoluteErrorL2Norm3D absVelocityErrorNormL2(u, uSol, indicatorF); + absVelocityErrorNormL2(result, input); + clout << "velocity-L2-error(abs)=" << result[0]; + SuperRelativeErrorL2Norm3D relVelocityErrorNormL2(u, uSol, indicatorF); + relVelocityErrorNormL2(result, input); + clout << "; velocity-L2-error(rel)=" << result[0] << std::endl; + + AnalyticalTemperaturePorousPlate3D tSol(Re, Pr, converter.getCharPhysLength(), converter.getCharPhysLowTemperature(), converter.getCharPhysTemperatureDifference()); + SuperLatticePhysTemperature3D t(ADlattice, converter); + + SuperAbsoluteErrorL2Norm3D absTemperatureErrorNormL2(t, tSol, indicatorF); + absTemperatureErrorNormL2(result, input); + clout << "temperature-L2-error(abs)=" << result[0]; + SuperRelativeErrorL2Norm3D relTemperatureErrorNormL2(t, tSol, indicatorF); + relTemperatureErrorNormL2(result, input); + clout << "; temperature-L2-error(rel)=" << result[0] << std::endl; + + AnalyticalHeatFluxPorousPlate3D HeatFluxSol(Re, Pr, converter.getCharPhysTemperatureDifference(), converter.getCharPhysLength(), converter.getThermalConductivity()); + SuperLatticePhysHeatFlux3D HeatFlux(ADlattice,converter); + + SuperAbsoluteErrorL2Norm3D absHeatFluxErrorNormL2(HeatFlux, HeatFluxSol, indicatorF); + absHeatFluxErrorNormL2(result, input); + clout << "heatFlux-L2-error(abs)=" << result[0]; + SuperRelativeErrorL2Norm3D relHeatFluxErrorNormL2(HeatFlux, HeatFluxSol, indicatorF); + relHeatFluxErrorNormL2(result, input); + clout << "; heatFlux-L2-error(rel)=" << result[0] << std::endl; +} + +/// Stores geometry information in form of material numbers +void prepareGeometry(SuperGeometry3D& superGeometry, + ThermalUnitConverter const&converter) +{ + + OstreamManager clout(std::cout,"prepareGeometry"); + clout << "Prepare Geometry ..." << std::endl; + + superGeometry.rename(0,2); + superGeometry.rename(2,1,0,1,0); + //superGeometry.clean(); + + std::vector extend( 3, T(0) ); + extend[0] = lx; + extend[1] = converter.getPhysLength(1); + extend[2] = lz; + std::vector origin( 3, T(0) ); + IndicatorCuboid3D bottom(extend, origin); + /// Set material number for bottom + superGeometry.rename(2,3,1,bottom); + + /// 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( ThermalUnitConverter const& converter, + SuperLattice3D& NSlattice, + SuperLattice3D& ADlattice, + Dynamics &bulkDynamics, + Dynamics& advectionDiffusionBulkDynamics, + sOnLatticeBoundaryCondition3D& NSboundaryCondition, + sOnLatticeBoundaryCondition3D& TboundaryCondition, + SuperGeometry3D& superGeometry ) +{ + + OstreamManager clout(std::cout,"prepareLattice"); + + T Tomega = converter.getLatticeThermalRelaxationFrequency(); + T NSomega = converter.getLatticeRelaxationFrequency(); + + /// define lattice Dynamics + clout << "defining dynamics" << endl; + + ADlattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics()); + NSlattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics()); + + ADlattice.defineDynamics(superGeometry.getMaterialIndicator({1, 2, 3}), &advectionDiffusionBulkDynamics); + NSlattice.defineDynamics(superGeometry.getMaterialIndicator({1, 2, 3}), &bulkDynamics); + + /// sets boundary + NSboundaryCondition.addVelocityBoundary(superGeometry.getMaterialIndicator({2, 3}), NSomega); + TboundaryCondition.addTemperatureBoundary(superGeometry.getMaterialIndicator({2, 3}), Tomega); +} + +void setBoundaryValues(ThermalUnitConverter const& converter, + SuperLattice3D& NSlattice, + SuperLattice3D& ADlattice, + int iT, SuperGeometry3D& superGeometry) +{ + + if (iT == 0) { + +// typedef advectionDiffusionLbHelpers TlbH; + + /// for each material set the defineRhoU and the Equilibrium + std::vector zero(3,T()); + AnalyticalConst3D u(zero); + AnalyticalConst3D rho(1.); + AnalyticalConst3D force(zero); + + T u_Re = converter.getLatticeVelocity( Re * converter.getPhysViscosity() / converter.getCharPhysLength() ); + AnalyticalConst3D u_top(converter.getCharLatticeVelocity(), u_Re, 0.0); + AnalyticalConst3D u_bot(0.0, u_Re, 0.0); + + NSlattice.defineRhoU(superGeometry, 1, rho, u); + NSlattice.iniEquilibrium(superGeometry, 1, rho, u); + NSlattice.defineField(superGeometry, 1, force); + NSlattice.defineRhoU(superGeometry, 2, rho, u_top); + NSlattice.iniEquilibrium(superGeometry, 2, rho, u_top); + NSlattice.defineField(superGeometry, 2, force); + NSlattice.defineRhoU(superGeometry, 3, rho, u_bot); + NSlattice.iniEquilibrium(superGeometry, 3, rho, u_bot); + NSlattice.defineField(superGeometry, 3, force); + + AnalyticalConst3D Cold(converter.getLatticeTemperature(Tcold)); + AnalyticalConst3D Hot(converter.getLatticeTemperature(Thot)); + + ADlattice.defineRho(superGeometry, 1, Cold); + ADlattice.iniEquilibrium(superGeometry, 1, Cold, u); + ADlattice.defineField(superGeometry, 1, u); + ADlattice.defineRho(superGeometry, 2, Hot); + ADlattice.iniEquilibrium(superGeometry, 2, Hot, u); + ADlattice.defineField(superGeometry, 2, u); + ADlattice.defineRho(superGeometry, 3, Cold); + ADlattice.iniEquilibrium(superGeometry, 3, Cold, u); + ADlattice.defineField(superGeometry, 3, u); + + /// Make the lattice ready for simulation + NSlattice.initialize(); + ADlattice.initialize(); + } +} + +void getResults(ThermalUnitConverter const& converter, + SuperLattice3D& NSlattice, + SuperLattice3D& ADlattice, int iT, + SuperGeometry3D& superGeometry, + Timer& timer, + bool converged) +{ + + OstreamManager clout(std::cout,"getResults"); + + SuperVTMwriter3D vtkWriter("thermalPorousPlate3d"); + + SuperLatticePhysVelocity3D velocity(NSlattice, converter); + SuperLatticePhysPressure3D pressure(NSlattice, converter); + SuperLatticePhysTemperature3D temperature(ADlattice, converter); + + AnalyticalHeatFluxPorousPlate3D HeatFluxSol(Re, Pr, converter.getCharPhysTemperatureDifference(), converter.getCharPhysLength(), converter.getThermalConductivity()); + SuperLatticePhysHeatFlux3D HeatFlux1(ADlattice,converter); + SuperLatticeFfromAnalyticalF3D HeatFluxSolLattice(HeatFluxSol,ADlattice); + vtkWriter.addFunctor( pressure ); + vtkWriter.addFunctor( velocity ); + vtkWriter.addFunctor( temperature ); + vtkWriter.addFunctor( HeatFlux1 ); + vtkWriter.addFunctor( HeatFluxSolLattice ); + + const int vtkIter = converter.getLatticeTime(100.); + + if (iT == 0) { + /// Writes the converter log file + //writeLogFile(converter,"thermalPorousPlate3d"); + + /// Writes the geometry, cuboid no. and rank no. as vti file for visualization + SuperLatticeGeometry3D geometry(NSlattice, superGeometry); + SuperLatticeCuboid3D cuboid(NSlattice); + SuperLatticeRank3D rank(NSlattice); + vtkWriter.write(geometry); + vtkWriter.write(cuboid); + vtkWriter.write(rank); + + vtkWriter.createMasterFile(); + } + + /// Writes the VTK files + if (iT%vtkIter == 0 || converged) { + NSlattice.getStatistics().print(iT,converter.getPhysTime(iT)); + timer.print(iT); + error(superGeometry, NSlattice, ADlattice, converter, Re); + + vtkWriter.write(iT); + + BlockReduction3D2D planeReduction( temperature, {0,0,1}, 600, BlockDataSyncMode::ReduceOnly ); + // write output as JPEG + heatmap::plotParam jpeg_Param; + jpeg_Param.maxValue = Thot; + jpeg_Param.minValue = Tcold; + heatmap::write(planeReduction, iT, jpeg_Param); + + /* BlockLatticeReduction3D planeReduction(temperature); + BlockGifWriter gifWriter; + gifWriter.write(planeReduction, iT, "temperature");*/ + } + +} + + +int main(int argc, char *argv[]) +{ + + /// === 1st Step: Initialization === + OstreamManager clout(std::cout,"main"); + olbInit(&argc, &argv); + singleton::directories().setOutputDir("./tmp/"); + + fstream f; + if (argc >= 2) { + N = atoi(argv[1]); + } + if (argc == 3) { + tau = atof(argv[2]); + } + + ThermalUnitConverter const converter( + (T) 1.0 / N, // physDeltaX + (T) 1.0 / N * 1.0 / 1e-3 * (tau - 0.5) / 3 / N, // physDeltaT + (T) 1.0, // charPhysLength + (T) sqrt( 9.81 * Ra * 1e-3 * 1e-3 / Pr / 9.81 / (Thot - Tcold) / pow(1.0, 3) * (Thot - Tcold) * 1.0 ), // charPhysVelocity + (T) 1e-3, // physViscosity + (T) 1.0, // physDensity + (T) 0.03, // physThermalConductivity + (T) Pr * 0.03 / 1e-3 / 1.0, // physSpecificHeatCapacity + (T) Ra * 1e-3 * 1e-3 / Pr / 9.81 / (Thot - Tcold) / pow(1.0, 3), // physThermalExpansionCoefficient + (T) Tcold, // charPhysLowTemperature + (T) Thot // charPhysHighTemperature + ); + converter.print(); + + /// === 2nd Step: Prepare Geometry === + lz = converter.getPhysDeltaX() * 3.; // depth of the channel + std::vector extend(3,T()); + extend[0] = lx; + extend[1] = ly; + extend[2] = lz; + std::vector origin(3,T()); + IndicatorCuboid3D cuboid(extend, origin); + + /// Instantiation of a cuboidGeometry with weights +#ifdef PARALLEL_MODE_MPI + const int noOfCuboids = singleton::mpi().getSize(); +#else + const int noOfCuboids = 7; +#endif + CuboidGeometry3D cuboidGeometry(cuboid, converter.getPhysDeltaX(), noOfCuboids); + cuboidGeometry.setPeriodicity(true,false, true); + + /// Instantiation of a loadBalancer + HeuristicLoadBalancer loadBalancer(cuboidGeometry); + + /// Instantiation of a superGeometry + SuperGeometry3D superGeometry(cuboidGeometry, loadBalancer, 2); + + prepareGeometry(superGeometry, converter); + + /// === 3rd Step: Prepare Lattice === + + SuperLattice3D ADlattice(superGeometry); + SuperLattice3D NSlattice(superGeometry); + + sOnLatticeBoundaryCondition3D NSboundaryCondition(NSlattice); + createLocalBoundaryCondition3D(NSboundaryCondition); + + sOnLatticeBoundaryCondition3D TboundaryCondition(ADlattice); + createAdvectionDiffusionBoundaryCondition3D(TboundaryCondition); + + ForcedBGKdynamics NSbulkDynamics( + converter.getLatticeRelaxationFrequency(), + instances::getBulkMomenta()); + + AdvectionDiffusionBGKdynamics TbulkDynamics ( + converter.getLatticeThermalRelaxationFrequency(), + instances::getAdvectionDiffusionBulkMomenta() + ); + + // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!// + // This coupling must be necessarily be put on the Navier-Stokes lattice!! + // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!// + + /* int nx = converter.numCells(lx) + 1; + int ny = converter.numCells(ly) + 1; + int nz = converter.numCells(lz) + 1; + */ + std::vector dir{0.0, 1.0, 0.0}; + + T boussinesqForcePrefactor = 9.81 / converter.getConversionFactorVelocity() * converter.getConversionFactorTime() * + converter.getCharPhysTemperatureDifference() * converter.getPhysThermalExpansionCoefficient(); + + NavierStokesAdvectionDiffusionCouplingGenerator3D + coupling(0, converter.getLatticeLength(lx), 0, converter.getLatticeLength(ly), 0, converter.getLatticeLength(lz), + boussinesqForcePrefactor, converter.getLatticeTemperature(Tcold), 1., dir); + + NSlattice.addLatticeCoupling(coupling, ADlattice); + + prepareLattice(converter, + NSlattice, ADlattice, + NSbulkDynamics, TbulkDynamics, + NSboundaryCondition, TboundaryCondition, superGeometry ); + + + /// === 4th Step: Main Loop with Timer === + Timer timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel() ); + timer.start(); + + util::ValueTracer converge(converter.getLatticeTime(1.0),epsilon); + for (int iT = 0; iT < converter.getLatticeTime(maxPhysT); ++iT) { + + if (converge.hasConverged()) { + clout << "Simulation converged." << endl; + getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged()); + break; + } + + /// === 5th Step: Definition of Initial and Boundary Conditions === + setBoundaryValues(converter, NSlattice, ADlattice, iT, superGeometry); + + /// === 6th Step: Collide and Stream Execution === + + NSlattice.collideAndStream(); + + NSlattice.executeCoupling(); + ADlattice.collideAndStream(); + + /// === 7th Step: Computation and Output of the Results === + getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged()); + converge.takeValue(ADlattice.getStatistics().getAverageEnergy()); + } + + timer.stop(); + timer.printSummary(); +} diff --git a/examples/thermal/rayleighBenard2d/Makefile b/examples/thermal/rayleighBenard2d/Makefile new file mode 100644 index 0000000..d4d2464 --- /dev/null +++ b/examples/thermal/rayleighBenard2d/Makefile @@ -0,0 +1,113 @@ +# 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 := +OUTPUT := rayleighBenard2d + +########################################################################### +## definitions + +include $(ROOT)/global.mk + +OBJECTS := $(foreach file, $(SRC) $(OUTPUT), $(PWD)/$(file).o) +DEPS := $(foreach file, $(SRC) $(OUTPUT), $(PWD)/$(file).d) + +########################################################################### +## all + +all : depend compile updatelib 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/*.* + +cleanobj: + @echo Clean object files + @rm -f $(OBJECTS) + +cleandep: + @echo Clean dependencies files + @rm -f $(DEPS) + +cleanbuild: + @echo Clean olb main + @cd $(ROOT); \ + $(MAKE) cleanbuild; + +########################################################################### +## update lib + +updatelib : + @cd $(ROOT); \ + $(MAKE) all; + +########################################################################### +## link + +link: $(OUTPUT) + +$(OUTPUT): $(OBJECTS) $(ROOT)/$(LIBDIR)/lib$(LIB).a + @echo Link $@ + $(CXX) $(foreach file, $(SRC), $(file).o) $@.o $(LDFLAGS) -L$(ROOT)/$(LIBDIR) -l$(LIB) -o $@ + +########################################################################### +## include dependencies + +ifneq "$(strip $(wildcard *.d))" "" + include $(foreach file,$(DEPS),$(file)) +endif + +########################################################################### +########################################################################### diff --git a/examples/thermal/rayleighBenard2d/definitions.mk b/examples/thermal/rayleighBenard2d/definitions.mk new file mode 100644 index 0000000..049cf73 --- /dev/null +++ b/examples/thermal/rayleighBenard2d/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 := rayleighBenard2d.cpp +OUTPUT := rayleighBenard2d diff --git a/examples/thermal/rayleighBenard2d/module.mk b/examples/thermal/rayleighBenard2d/module.mk new file mode 100644 index 0000000..1190482 --- /dev/null +++ b/examples/thermal/rayleighBenard2d/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/thermal/rayleighBenard2d/rayleighBenard2d.cpp b/examples/thermal/rayleighBenard2d/rayleighBenard2d.cpp new file mode 100644 index 0000000..1fe1385 --- /dev/null +++ b/examples/thermal/rayleighBenard2d/rayleighBenard2d.cpp @@ -0,0 +1,346 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2008 Orestis Malaspinas + * 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. + */ + +/* rayleighBenard2d.cpp: + * Rayleigh-Benard convection rolls in 2D, simulated with + * the thermal LB model by Z. Guo e.a., between a hot plate at + * the bottom and a cold plate at the top. + */ + + +#include "olb2D.h" +#include "olb2D.hh" // use only generic version! + +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace std; + +typedef double T; + +#define NSDESCRIPTOR D2Q9 +#define TDESCRIPTOR D2Q5 + +// Parameters for the simulation setup +const T lx = 2.; // length of the channel +const T ly = 1.; // height of the channel +const int N = 10; // resolution of the model +const T Ra = 1e4; // Rayleigh number +const T Pr = 0.71; // Prandtl numbe