Initialize at openlb-1-3

7 files changed, 791 insertions, 0 deletions

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
+#  <http://www.openlb.net/>
+#
+#  This program is free software; you can redistribute it and/or
+#  modify it under the terms of the GNU General Public License
+#  as published by the Free Software Foundation; either version 2
+#  of the License, or (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public
+#  License along with this program; if not, write to the Free
+#  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+#  Boston, MA  02110-1301, USA.
+
+###########################################################################
+## definitions
+
+include definitions.mk
+
+include $(ROOT)/global.mk
+
+OBJECTS    := $(foreach file, $(SRC), $(PWD)/$(file:.cpp=.o))
+DEPS       := $(foreach file, $(SRC), $(PWD)/$(file:.cpp=.d))
+
+###########################################################################
+## all
+
+all : depend compile link
+
+
+###########################################################################
+## dependencies
+
+depend : $(DEPS)
+
+$(PWD)/%.d : %.cpp
+	@echo Create dependencies for $<
+	@$(SHELL) -ec '$(CXX) -M $(CXXFLAGS) $(IDIR) $< \
+                       | sed -e "s!$*\.o!$(PWD)\/$*\.o!1" > .tmpfile; \
+                       cp -f .tmpfile $@;'
+
+###########################################################################
+## compile
+
+compile : $(OBJECTS)
+
+$(PWD)/%.o: %.cpp
+	@echo Compile $<
+	$(CXX) $(CXXFLAGS) $(IDIR) -c $< -o $@
+
+###########################################################################
+## clean
+
+clean : cleanrub cleanobj cleandep
+
+cleanrub:
+	@echo Clean rubbish files
+	@rm -f *~ core .tmpfile tmp/*.* $(OUTPUT)
+	@rm -f tmp/vtkData/*.* tmp/vtkData/data/*.* tmp/imageData/*.* tmp/gnuplotData/*.* tmp/gnuplotData/data/*.*
+
+cleanobj:
+	@echo Clean object files
+	@rm -f $(OBJECTS)
+
+cleandep:
+	@echo Clean dependencies files
+	@rm -f $(DEPS)
+
+cleanbuild:
+	@cd $(ROOT); \
+	 $(MAKE) cleanlib;
+
+###########################################################################
+## update lib
+
+$(ROOT)/$(LIBDIR)/lib$(LIB).a :
+	@cd $(ROOT); \
+	$(MAKE) all
+
+###########################################################################
+## link
+
+link: $(OUTPUT)
+
+$(OUTPUT): $(OBJECTS) $(ROOT)/$(LIBDIR)/lib$(LIB).a
+	@echo Link $@
+	$(CXX) $(foreach file, $(SRC), $(file:.cpp=.o)) $(LDFLAGS) -L$(ROOT)/$(LIBDIR) -l$(LIB) -lz -o $@
+
+###########################################################################
+## include dependencies
+
+ifneq "$(strip $(wildcard *.d))" ""
+   include $(foreach file,$(DEPS),$(file))
+endif
+
+###########################################################################
+###########################################################################
diff --git a/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
+#  <http://www.openlb.net/>
+#
+#  This program is free software; you can redistribute it and/or
+#  modify it under the terms of the GNU General Public License
+#  as published by the Free Software Foundation; either version 2
+#  of the License, or (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public
+#  License along with this program; if not, write to the Free
+#  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+#  Boston, MA  02110-1301, USA.
+
+
+###########################################################################
+###########################################################################
+## DEFINITIONS TO BE CHANGED
+
+ROOT            := ../../..
+SRC             := 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
+#  <http://www.openlb.net/>
+#
+#  This program is free software; you can redistribute it and/or
+#  modify it under the terms of the GNU General Public License
+#  as published by the Free Software Foundation; either version 2
+#  of the License, or (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public
+#  License along with this program; if not, write to the Free
+#  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+#  Boston, MA  02110-1301, USA.
+
+current_dir := $(dir $(word $(words $(MAKEFILE_LIST)),$(MAKEFILE_LIST)))
+
+include global.mk
+include rules.mk
+include $(addsuffix definitions.mk, $(current_dir))
+
+$(eval $(call sample,$(current_dir)$(OUTPUT),$(addprefix $(current_dir), $(SRC))))
diff --git a/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
+ *  <http://www.openlb.net/>
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; either version 2
+ *  of the License, or (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public
+ *  License along with this program; if not, write to the Free
+ *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ *  Boston, MA  02110-1301, USA.
+ */
+
+/* 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 <vector>
+#include <cmath>
+#include <iostream>
+#include <fstream>
+
+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<vector<T>> values_DNS;
+
+template <typename T, typename _DESCRIPTOR>
+class Tgv3D : public AnalyticalF3D<T,T> {
+
+protected:
+  T u0;
+
+// initial solution of the TGV
+public:
+  Tgv3D(UnitConverter<T,_DESCRIPTOR> const& converter, T frac) : AnalyticalF3D<T,S>(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<T>& 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<T,DESCRIPTOR>& sLattice,
+                    UnitConverter<T,DESCRIPTOR> const& converter,
+                    Dynamics<T, DESCRIPTOR>& bulkDynamics,
+                    SuperGeometry3D<T>& superGeometry)
+{
+
+  OstreamManager clout(std::cout,"prepareLattice");
+  clout << "Prepare Lattice ..." << std::endl;
+
+  /// Material=0 -->do nothing
+  sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
+  /// Material=1 -->bulk dynamics
+  sLattice.defineDynamics(superGeometry, 1, &bulkDynamics);
+
+  sLattice.initialize();
+  clout << "Prepare Lattice ... OK" << std::endl;
+}
+
+void setBoundaryValues(SuperLattice3D<T, DESCRIPTOR>& sLattice,
+                       UnitConverter<T,DESCRIPTOR> const& converter,
+                       SuperGeometry3D<T>& superGeometry)
+{
+  OstreamManager clout(std::cout,"setBoundaryValues");
+
+  AnalyticalConst3D<T,T> rho(1.);
+  Tgv3D<T,DESCRIPTOR> uSol(converter, 1);
+
+  sLattice.defineRhoU(superGeometry, 1, rho, uSol);
+  sLattice.iniEquilibrium(superGeometry, 1, rho, uSol);
+
+  sLattice.initialize();
+}
+
+// Interpolate the data points to the output interval
+void getDNSValues()
+{
+  string file_name;
+  //Brachet, Marc E., et al. "Small-scale structure of the Taylor–Green vortex." Journal of Fluid Mechanics 130 (1983): 411-452; Figure 7
+  if (abs(Re - 800.0) < numeric_limits<T>::epsilon() && maxPhysT <= 10.0 + numeric_limits<T>::epsilon()) {
+    file_name= "Re800_Brachet.inp";
+  }
+  else if (abs(Re - 1600.0) < numeric_limits<T>::epsilon() && maxPhysT <= 10.0 + numeric_limits<T>::epsilon()) {
+    file_name = "Re1600_Brachet.inp";
+  }
+  else if (abs(Re - 3000.0) < numeric_limits<T>::epsilon() && maxPhysT <= 10.0 + numeric_limits<T>::epsilon()) {
+    file_name = "Re3000_Brachet.inp";
+  }
+  else {
+    std::cout<<"Reynolds number not supported or maxPhysT>10: DNS plot will be disabled"<<std::endl;
+    plotDNS = false;
+    return;
+  }
+  std::ifstream data(file_name);
+  std::string line;
+  std::vector<vector<T>> parsedDat;
+  while (std::getline(data, line)) {
+    std::stringstream lineStream(line);
+    std::string cell;
+    std::vector<T> parsedRow;
+    while (std::getline(lineStream, cell, ' ')) {
+      parsedRow.push_back(atof(cell.c_str()));
+
+    }
+    if (parsedDat.size() > 0 && parsedRow.size() > 1) {
+      parsedRow.push_back((parsedRow[1] - parsedDat[parsedDat.size() - 1][1]) / (parsedRow[0] - parsedDat[parsedDat.size()-1][0]));
+      parsedRow.push_back(parsedDat[parsedDat.size()-1][1] - parsedRow[2] * parsedDat[parsedDat.size()-1][0]);
+    }
+
+    parsedDat.push_back(parsedRow);
+  }
+
+  int steps = maxPhysT / gnuplotSave + 1.5;
+  for (int i=0; i < steps; i++) {
+    std::vector<T> inValues_temp;
+    inValues_temp.push_back(i * gnuplotSave);
+    if (inValues_temp[0] < parsedDat[0][0]) {
+      inValues_temp.push_back(parsedDat[1][2] * inValues_temp[0] + parsedDat[1][3]);
+    }
+    else if (inValues_temp[0] > parsedDat[parsedDat.size()-1][0]) {
+      inValues_temp.push_back(parsedDat[parsedDat.size()-1][2] * inValues_temp[0] +
+                              parsedDat[parsedDat.size()-1][3]);
+    }
+    else {
+
+      for (size_t j=0; j < parsedDat.size()-1; j++)  {
+        if (inValues_temp[0] > parsedDat[j][0] && inValues_temp[0] < parsedDat[j+1][0]) {
+          inValues_temp.push_back(parsedDat[j+1][2] * inValues_temp[0] + parsedDat[j+1][3]);
+        }
+      }
+    }
+    values_DNS.push_back(inValues_temp);
+  }
+}
+
+
+
+void getResults(SuperLattice3D<T, DESCRIPTOR>& sLattice,
+                UnitConverter<T,DESCRIPTOR> const& converter, int iT,
+                SuperGeometry3D<T>& superGeometry, Timer<double>& timer,
+                Dynamics<T, DESCRIPTOR>* bulkDynamics)
+{
+  OstreamManager clout(std::cout,"getResults");
+
+  SuperVTMwriter3D<T> vtmWriter("tgv3d");
+
+  if (iT == 0) {
+    // Writes the geometry, cuboid no. and rank no. as vti file for visualization
+    SuperLatticeGeometry3D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
+    SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid(sLattice);
+    SuperLatticeRank3D<T, DESCRIPTOR> rank(sLattice);
+    superGeometry.rename(0,2);
+    vtmWriter.write(geometry);
+    vtmWriter.write(cuboid);
+    vtmWriter.write(rank);
+    vtmWriter.createMasterFile();
+    if (plotDNS==true) {
+      getDNSValues();
+    }
+  }
+  if (iT%converter.getLatticeTime(vtkSave) == 0) {
+    SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity(sLattice, converter);
+    SuperLatticePhysPressure3D<T, DESCRIPTOR> pressure(sLattice, converter);
+    vtmWriter.addFunctor( velocity );
+    vtmWriter.addFunctor( pressure );
+    vtmWriter.write(iT);
+
+    // write output of velocity as JPEG
+    SuperEuklidNorm3D<T, DESCRIPTOR> normVel( velocity );
+    BlockReduction3D2D<T> planeReduction( normVel, {0, 0, 1} );
+    heatmap::write(planeReduction, iT);
+
+    timer.update(iT);
+    timer.printStep(2);
+    sLattice.getStatistics().print(iT,converter.getPhysTime( iT ));
+  }
+
+  static Gnuplot<T> gplot("Turbulence_Dissipation_Rate");
+
+  if (iT%converter.getLatticeTime(gnuplotSave) == 0) {
+
+    int input[3];
+    T output[1];
+
+#if defined (finiteDiff)
+    std::list<int> matNumber;
+    matNumber.push_back(1);
+    SuperLatticePhysDissipationFD3D<T, DESCRIPTOR> diss(superGeometry, sLattice, matNumber, converter);
+#if !defined (DNS)
+    SuperLatticePhysEffectiveDissipationFD3D<T, DESCRIPTOR> effectiveDiss(superGeometry, sLattice, matNumber,
+                                                                          converter, *(dynamic_cast<LESDynamics<T,DESCRIPTOR>*>(bulkDynamics)));
+#endif
+#else
+    SuperLatticePhysDissipation3D<T, DESCRIPTOR> diss(sLattice, converter);
+#if !defined (DNS)
+    SuperLatticePhysEffevtiveDissipation3D<T, DESCRIPTOR> effectiveDiss(sLattice, converter, smagoConst, *(dynamic_cast<LESDynamics<T,DESCRIPTOR>*>(bulkDynamics)));
+#endif
+#endif
+    SuperIntegral3D<T> integralDiss(diss, superGeometry, 1);
+    integralDiss(output, input);
+    T diss_mol = output[0];
+    diss_mol /= volume;
+    T diss_eff = diss_mol;
+
+#if !defined (DNS)
+    SuperIntegral3D<T> integralEffectiveDiss(effectiveDiss, superGeometry, 1);
+    integralEffectiveDiss(output, input);
+    diss_eff = output[0];
+    diss_eff /= volume;
+#endif
+
+    T diss_eddy = diss_eff - diss_mol;
+    if(plotDNS==true) {
+     int step = converter.getPhysTime(iT) / gnuplotSave + 0.5;
+     gplot.setData(converter.getPhysTime(iT), {diss_mol, diss_eddy, diss_eff, values_DNS[step][1]}, {"molecular dissipation rate", "eddy dissipation rate", "effective dissipation rate" ,"Brachet et al."}, "bottom right");
+    } else {
+     gplot.setData(converter.getPhysTime(iT), {diss_mol, diss_eddy, diss_eff}, {"molecular dissipation rate", "eddy dissipation rate", "effective dissipation rate"}, "bottom right");
+    }
+    gplot.writePNG();
+  }
+
+  /// write pdf at last time step
+  if (iT == converter.getLatticeTime(maxPhysT)-1) {
+    gplot.writePDF();
+  }
+  return;
+}
+
+int main(int argc, char* argv[])
+{
+
+  /// === 1st Step: Initialization ===
+  olbInit(&argc, &argv);
+  singleton::directories().setOutputDir("./tmp/");
+  OstreamManager clout( std::cout,"main" );
+
+  UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR> converter(
+    int {int(std::nearbyint(N/(2*pi)))},        // resolution: number of voxels per charPhysL
+    (T)   0.507639, // 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)   1./Re,    // 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("tgv3d");
+
+#ifdef PARALLEL_MODE_MPI
+  const int noOfCuboids = 2 * singleton::mpi().getSize();
+#else
+  const int noOfCuboids = 1;
+#endif
+
+  CuboidGeometry3D<T> cuboidGeometry(0, 0, 0, converter.getConversionFactorLength(), N, N, N, noOfCuboids);
+
+  cuboidGeometry.setPeriodicity(true, true, true);
+
+  HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
+
+  // === 2nd Step: Prepare Geometry ===
+  SuperGeometry3D<T> superGeometry(cuboidGeometry, loadBalancer, 3);
+  prepareGeometry(superGeometry);
+
+  /// === 3rd Step: Prepare Lattice ===
+  SuperLattice3D<T, DESCRIPTOR> sLattice(superGeometry);
+
+  std::unique_ptr<Dynamics<T, DESCRIPTOR>> bulkDynamics;
+  const T omega = converter.getLatticeRelaxationFrequency();
+#if defined(RLB)
+  bulkDynamics.reset(new RLBdynamics<T, DESCRIPTOR>(omega, instances::getBulkMomenta<T, DESCRIPTOR>()));
+#elif defined(DNS)
+  bulkDynamics.reset(new BGKdynamics<T, DESCRIPTOR>(omega, instances::getBulkMomenta<T, DESCRIPTOR>()));
+#elif defined(WALE)
+  bulkDynamics.reset(new WALEBGKdynamics<T, DESCRIPTOR>(omega, instances::getBulkMomenta<T, DESCRIPTOR>(),
+      smagoConst));
+#elif defined(ShearSmagorinsky)
+  bulkDynamics.reset(new ShearSmagorinskyBGKdynamics<T, DESCRIPTOR>(omega, instances::getBulkMomenta<T, DESCRIPTOR>(),
+      smagoConst));
+#elif defined(Krause)
+  bulkDynamics.reset(new KrauseBGKdynamics<T, DESCRIPTOR>(omega, instances::getBulkMomenta<T, DESCRIPTOR>(),
+      smagoConst));
+#elif defined(ConsistentStrainSmagorinsky)
+  bulkDynamics.reset(new ConStrainSmagorinskyBGKdynamics<T, DESCRIPTOR>(omega, instances::getBulkMomenta<T, DESCRIPTOR>(),
+      smagoConst));
+#else //DNS Simulation
+
+  bulkDynamics.reset(new SmagorinskyBGKdynamics<T, DESCRIPTOR>(omega, instances::getBulkMomenta<T, DESCRIPTOR>(),
+      smagoConst));
+#endif
+
+  prepareLattice(sLattice, converter, *bulkDynamics, superGeometry);
+
+#if defined(WALE)
+  std::list<int> mat;
+  mat.push_back(1);
+  std::unique_ptr;SuperLatticeF3D<T, DESCRIPTOR>> functor(new SuperLatticeVelocityGradientFD3D<T, DESCRIPTOR>(superGeometry, sLattice, mat));
+#endif
+
+  /// === 4th Step: Main Loop with Timer ===
+  Timer<double> timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel() );
+  timer.start();
+
+  // === 5th Step: Definition of Initial and Boundary Conditions ===
+  setBoundaryValues(sLattice, converter, superGeometry);
+
+  for (int iT = 0; iT <= converter.getLatticeTime(maxPhysT); ++iT) {
+#if defined(WALE)
+    sLattice.defineField<descriptors::VELO_GRAD>(superGeometry, 1, *functor);
+#endif
+
+    /// === 6th Step: Computation and Output of the Results ===
+    getResults(sLattice, converter, iT, superGeometry, timer, bulkDynamics.get());
+
+    /// === 7th Step: Collide and Stream Execution ===
+    sLattice.collideAndStream();
+  }
+  timer.stop();
+  timer.printSummary();
+
+  return 0;
+}