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/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 +++++++++++++++++++++++++++
7 files changed, 791 insertions(+)
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
(limited to 'examples/turbulence/tgv3d')
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;
+}
--
cgit v1.2.3