Initialize at openlb-1-3

1 files changed, 421 insertions, 0 deletions

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;
+}