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/tgv3d.cpp | 421 ++++++++++++++++++++++++++++++++++++ 1 file changed, 421 insertions(+) create mode 100644 examples/turbulence/tgv3d/tgv3d.cpp (limited to 'examples/turbulence/tgv3d/tgv3d.cpp') diff --git a/examples/turbulence/tgv3d/tgv3d.cpp b/examples/turbulence/tgv3d/tgv3d.cpp new file mode 100644 index 0000000..5a1fda6 --- /dev/null +++ b/examples/turbulence/tgv3d/tgv3d.cpp @@ -0,0 +1,421 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2017 Mathias J. Krause, Patrick Nathan, Alejandro C. Barreto, Marc Haußmann + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program; if not, write to the Free + * Software Foundation, Inc., 51 Franklin Street, Fifth Floor, + * Boston, MA 02110-1301, USA. + */ + +/* tgv3d.cpp: + * The Taylor-Green-Vortex (TGV) is one of the simplest configuration, + * where you can investigate the generation of small structures and the resulting turbulence. + * The 2pi periodic box domain and the single mode initial conditions contribute to the simplicity. + * In consequence, the TGV is a common benchmark case for + * Direct Numerical Simulations (DNS) and Large Eddy Simulations (LES). + * + * This example shows the usage and the effects of different subgrid scale turbulence models. + * The molecular dissipation rate, the eddy dissipation rate and + * the effective dissipation rate are calculated and plotted over the simulation time. + * This results can be compared with a published DNS solution, e.g. + * Brachet, Marc E., et al. "Small-scale structure of the Taylor–Green vortex." + * Journal of Fluid Mechanics 130 (1983): 411-452. + */ + + +#include "olb3D.h" +#ifndef OLB_PRECOMPILED // Unless precompiled version is used, +#include "olb3D.hh" // include full template code +#endif +#include +#include +#include +#include + +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace olb::util; +using namespace std; + +typedef double T; + +// Choose your turbulent model of choice +//#define RLB +#define Smagorinsky +//#define WALE +//#define ConsistentStrainSmagorinsky +//#define ShearSmagorinsky +//#define Krause +//#define DNS + +#define finiteDiff //for N<256 + +#ifdef ShearSmagorinsky +#define DESCRIPTOR ShearSmagorinskyD3Q19Descriptor +#elif defined (WALE) +#define DESCRIPTOR WALED3Q19Descriptor +#else +#define DESCRIPTOR D3Q19<> +#endif + +// Global constants +const T pi = 4.0 * std::atan(1.0); +const T volume = pow(2. * pi, 3.); // volume of the 2pi periodic box + + +// Parameters for the simulation setup +const T maxPhysT = 10; // max. simulation time in s, SI unit +int N = 128; // resolution of the model +T Re = 800; // defined as 1/kinematic viscosity +T smagoConst = 0.1; // Smagorisky Constant, for ConsistentStrainSmagorinsky smagoConst = 0.033 +T vtkSave = 0.25; // time interval in s for vtk output +T gnuplotSave = 0.1; // time interval in s for gnuplot output + +bool plotDNS = true; //available for Re=800, Re=1600, Re=3000 (maxPhysT<=10) +vector> values_DNS; + +template +class Tgv3D : public AnalyticalF3D { + +protected: + T u0; + +// initial solution of the TGV +public: + Tgv3D(UnitConverter const& converter, T frac) : AnalyticalF3D(3) + { + u0 = converter.getCharLatticeVelocity(); + }; + + bool operator()(T output[], const S input[]) override + { + T x = input[0]; + T y = input[1]; + T z = input[2]; + + output[0] = u0 * sin(x) * cos(y) * cos(z); + output[1] = -u0 * cos(x) * sin(y) * cos(z); + output[2] = 0; + + return true; + }; +}; + +void prepareGeometry(SuperGeometry3D& superGeometry) +{ + OstreamManager clout(std::cout,"prepareGeometry"); + clout << "Prepare Geometry ..." << std::endl; + + superGeometry.rename(0,1); + superGeometry.communicate(); + + /// Removes all not needed boundary voxels outside the surface + superGeometry.clean(); + /// Removes all not needed boundary voxels inside the surface + superGeometry.innerClean(); + superGeometry.checkForErrors(); + + superGeometry.print(); + clout << "Prepare Geometry ... OK" << std::endl; + return; +} + +// Set up the geometry of the simulation +void prepareLattice(SuperLattice3D& sLattice, + UnitConverter const& converter, + Dynamics& bulkDynamics, + SuperGeometry3D& superGeometry) +{ + + OstreamManager clout(std::cout,"prepareLattice"); + clout << "Prepare Lattice ..." << std::endl; + + /// Material=0 -->do nothing + sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics()); + /// Material=1 -->bulk dynamics + sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); + + sLattice.initialize(); + clout << "Prepare Lattice ... OK" << std::endl; +} + +void setBoundaryValues(SuperLattice3D& sLattice, + UnitConverter const& converter, + SuperGeometry3D& superGeometry) +{ + OstreamManager clout(std::cout,"setBoundaryValues"); + + AnalyticalConst3D rho(1.); + Tgv3D 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::epsilon() && maxPhysT <= 10.0 + numeric_limits::epsilon()) { + file_name= "Re800_Brachet.inp"; + } + else if (abs(Re - 1600.0) < numeric_limits::epsilon() && maxPhysT <= 10.0 + numeric_limits::epsilon()) { + file_name = "Re1600_Brachet.inp"; + } + else if (abs(Re - 3000.0) < numeric_limits::epsilon() && maxPhysT <= 10.0 + numeric_limits::epsilon()) { + file_name = "Re3000_Brachet.inp"; + } + else { + std::cout<<"Reynolds number not supported or maxPhysT>10: DNS plot will be disabled"<> parsedDat; + while (std::getline(data, line)) { + std::stringstream lineStream(line); + std::string cell; + std::vector 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 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& sLattice, + UnitConverter const& converter, int iT, + SuperGeometry3D& superGeometry, Timer& timer, + Dynamics* bulkDynamics) +{ + OstreamManager clout(std::cout,"getResults"); + + SuperVTMwriter3D vtmWriter("tgv3d"); + + if (iT == 0) { + // Writes the geometry, cuboid no. and rank no. as vti file for visualization + SuperLatticeGeometry3D geometry(sLattice, superGeometry); + SuperLatticeCuboid3D cuboid(sLattice); + SuperLatticeRank3D rank(sLattice); + 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 velocity(sLattice, converter); + SuperLatticePhysPressure3D pressure(sLattice, converter); + vtmWriter.addFunctor( velocity ); + vtmWriter.addFunctor( pressure ); + vtmWriter.write(iT); + + // write output of velocity as JPEG + SuperEuklidNorm3D normVel( velocity ); + BlockReduction3D2D planeReduction( normVel, {0, 0, 1} ); + heatmap::write(planeReduction, iT); + + timer.update(iT); + timer.printStep(2); + sLattice.getStatistics().print(iT,converter.getPhysTime( iT )); + } + + static Gnuplot gplot("Turbulence_Dissipation_Rate"); + + if (iT%converter.getLatticeTime(gnuplotSave) == 0) { + + int input[3]; + T output[1]; + +#if defined (finiteDiff) + std::list matNumber; + matNumber.push_back(1); + SuperLatticePhysDissipationFD3D diss(superGeometry, sLattice, matNumber, converter); +#if !defined (DNS) + SuperLatticePhysEffectiveDissipationFD3D effectiveDiss(superGeometry, sLattice, matNumber, + converter, *(dynamic_cast*>(bulkDynamics))); +#endif +#else + SuperLatticePhysDissipation3D diss(sLattice, converter); +#if !defined (DNS) + SuperLatticePhysEffevtiveDissipation3D effectiveDiss(sLattice, converter, smagoConst, *(dynamic_cast*>(bulkDynamics))); +#endif +#endif + SuperIntegral3D integralDiss(diss, superGeometry, 1); + integralDiss(output, input); + T diss_mol = output[0]; + diss_mol /= volume; + T diss_eff = diss_mol; + +#if !defined (DNS) + SuperIntegral3D 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 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 cuboidGeometry(0, 0, 0, converter.getConversionFactorLength(), N, N, N, noOfCuboids); + + cuboidGeometry.setPeriodicity(true, true, true); + + HeuristicLoadBalancer loadBalancer(cuboidGeometry); + + // === 2nd Step: Prepare Geometry === + SuperGeometry3D superGeometry(cuboidGeometry, loadBalancer, 3); + prepareGeometry(superGeometry); + + /// === 3rd Step: Prepare Lattice === + SuperLattice3D sLattice(superGeometry); + + std::unique_ptr> bulkDynamics; + const T omega = converter.getLatticeRelaxationFrequency(); +#if defined(RLB) + bulkDynamics.reset(new RLBdynamics(omega, instances::getBulkMomenta())); +#elif defined(DNS) + bulkDynamics.reset(new BGKdynamics(omega, instances::getBulkMomenta())); +#elif defined(WALE) + bulkDynamics.reset(new WALEBGKdynamics(omega, instances::getBulkMomenta(), + smagoConst)); +#elif defined(ShearSmagorinsky) + bulkDynamics.reset(new ShearSmagorinskyBGKdynamics(omega, instances::getBulkMomenta(), + smagoConst)); +#elif defined(Krause) + bulkDynamics.reset(new KrauseBGKdynamics(omega, instances::getBulkMomenta(), + smagoConst)); +#elif defined(ConsistentStrainSmagorinsky) + bulkDynamics.reset(new ConStrainSmagorinskyBGKdynamics(omega, instances::getBulkMomenta(), + smagoConst)); +#else //DNS Simulation + + bulkDynamics.reset(new SmagorinskyBGKdynamics(omega, instances::getBulkMomenta(), + smagoConst)); +#endif + + prepareLattice(sLattice, converter, *bulkDynamics, superGeometry); + +#if defined(WALE) + std::list mat; + mat.push_back(1); + std::unique_ptr;SuperLatticeF3D> functor(new SuperLatticeVelocityGradientFD3D(superGeometry, sLattice, mat)); +#endif + + /// === 4th Step: Main Loop with Timer === + Timer 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(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