Added new LBM tutorial for backward facing step setup

apps/tutorials/lbm/05_BackwardFacingStep.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla 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 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file 05_BackwardFacingStep.cpp
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Amin Nabikhani
+//
+//======================================================================================================================
+
+
+#include "blockforest/all.h"
+#include "core/all.h"
+#include "domain_decomposition/all.h"
+#include "field/all.h"
+#include "geometry/all.h"
+#include "lbm/all.h"
+#include "timeloop/all.h"
+
+namespace walberla {
+
+//! [typedefs]
+using LatticeModel_T = lbm::D2Q9<lbm::collision_model::SRT>;
+using Stencil_T = LatticeModel_T::Stencil;
+using CommunicationStencil_T = LatticeModel_T::CommunicationStencil;
+//! [typedefs]
+
+using PdfField_T = lbm::PdfField<LatticeModel_T>;
+
+using flag_t = walberla::uint8_t;
+using FlagField_T = FlagField<flag_t>;
+
+//**********************************************************************************************************************
+/*!
+*   Class for determining (and logging) the normalized reattachment length after the backward-facing step
+*/
+//**********************************************************************************************************************
+class ReattachmentLengthFinder
+{
+public:
+
+   ReattachmentLengthFinder( const shared_ptr< StructuredBlockStorage > & blocks,
+                             const BlockDataID & pdfFieldID, const BlockDataID & flagFieldID, const FlagUID & fluid,
+                             const std::string & filename, const uint_t checkFrequency, const Vector3< real_t > stepSize) :
+      blocks_( blocks ),
+      pdfFieldID_( pdfFieldID ), flagFieldID_( flagFieldID ), fluid_( fluid ),
+      filename_( filename ), checkFrequency_( checkFrequency ), stepSize_( stepSize ),
+      executionCounter_( uint_c(0) )
+   {
+      // open and close file on root process - purpose: If the file already exists, its content will be erased.
+      WALBERLA_ROOT_SECTION()
+      {
+         std::ofstream fileLocBottom( filename_.c_str());
+         fileLocBottom << "Time,[Locations on the BOTTOM Wall that Reattachment Occures (Normalized with Step Height) ]" << std::endl;
+         fileLocBottom.close();
+      }
+   }
+
+   void operator()()
+   {
+
+      ++executionCounter_;
+      if( ( executionCounter_ - uint_c(1) ) % checkFrequency_ != 0 )
+         return;
+
+      // variables for storing the process local reattachment location
+      std::vector<real_t> reattachmentLocations;
+
+      // iterate all blocks stored locally on this process
+      for( auto block = blocks_->begin(); block != blocks_->end(); ++block )
+      {
+         // retrieve the pdf field and the flag field from the block
+         PdfField_T  *  pdfField = block->getData< PdfField_T >( pdfFieldID_ );
+         FlagField_T * flagField = block->getData< FlagField_T >( flagFieldID_ );
+
+         // retrieve the lattice model from the pdf field
+         const auto & latticeModel = pdfField->latticeModel();
+
+         // get the flag that marks a cell as being fluid
+         auto fluid = flagField->getFlag( fluid_ );
+
+         // iterate all cells of the current block
+         for( auto cell = pdfField->beginXYZ(); cell != pdfField->end(); ++cell )
+         {
+            const cell_idx_t x = cell.x();
+            const cell_idx_t y = cell.y();
+            const cell_idx_t z = cell.z();
+
+            // If the current cell is marked as being fluid ...
+            if( flagField->isFlagSet( x, y, z, fluid ) )
+            {
+               Cell currentPosition( x, y, z );
+               blocks_->transformBlockLocalToGlobalCell( currentPosition, *block );
+
+               // only consider the bottom row
+               if(currentPosition[1] == uint_t(0))
+               {
+                  Vector3< real_t > velocity;
+                  Vector3< real_t > vel_left;
+                  Vector3< real_t > vel_right;
+                  getVelocity( velocity, latticeModel, cell );
+                  pdfField->getVelocity( vel_left, x-int_c(1), y, z );
+                  pdfField->getVelocity( vel_right, x+int_c(1), y, z );
+                  if( (vel_left[0] * vel_right[0]) < real_c(0) && velocity[0] > real_c(0))
+                  {
+                     real_t xPosition = blocks_->getBlockLocalCellCenter(*block, currentPosition)[0];
+                     reattachmentLocations.push_back( (xPosition - stepSize_[0]) / stepSize_[1] );
+                  }
+               }
+            }
+         }
+      }
+
+      WALBERLA_ROOT_SECTION() {
+         std::ofstream fileLocBottom(filename_.c_str(), std::ios_base::app);
+         fileLocBottom << executionCounter_ << " ";
+         for (auto i = reattachmentLocations.begin(); i != reattachmentLocations.end(); i++) {
+            fileLocBottom << *i << " ";
+         }
+         fileLocBottom << std::endl;
+         fileLocBottom.close();
+      }
+   }
+
+private:
+
+   shared_ptr< StructuredBlockStorage > blocks_;
+   BlockDataID  pdfFieldID_;
+   BlockDataID flagFieldID_;
+   FlagUID fluid_;
+   std::string filename_;
+   uint_t checkFrequency_;
+   Vector3< real_t > stepSize_;
+
+   uint_t executionCounter_;
+};
+
+
+int main( int argc, char ** argv )
+{
+   walberla::Environment walberlaEnv( argc, argv );
+   
+   auto blocks = blockforest::createUniformBlockGridFromConfig( walberlaEnv.config() );
+
+   // read parameters
+   auto parameters = walberlaEnv.config()->getOneBlock( "Parameters" );
+   auto domain     = walberlaEnv.config()->getOneBlock( "DomainSetup" );
+   auto boundariesConfig = walberlaEnv.config()->getOneBlock( "Boundaries" );
+
+   //! [Params]
+   const real_t Re = parameters.getParameter< real_t >( "Re", real_c( 1000 ) );
+   const real_t uLB = parameters.getParameter< real_t >( "uLB", real_c( 0.01 )  );
+   //! [Params]
+   //! [height]
+   const Vector3<real_t> domainsize = domain.getParameter< Vector3<real_t> >( "cellsPerBlock", Vector3<real_t>() );
+  //! [height]
+   const uint_t timesteps = parameters.getParameter< uint_t > ( "timesteps", uint_c( 100000 ) );   
+
+   const Vector3<real_t> inletVelocity(uLB, real_c(0), real_c(0));
+   //!  [Omega]
+   const real_t viscosity = (domainsize[1] * uLB) / Re;
+   const real_t omega = lbm::collision_model::omegaFromViscosity(viscosity);
+   //!  [Omega] 
+   WALBERLA_LOG_INFO_ON_ROOT( "Re =  " << Re );
+   WALBERLA_LOG_INFO_ON_ROOT( "uLB =  " << uLB );
+   WALBERLA_LOG_INFO_ON_ROOT( "timesteps =  " << timesteps );
+   WALBERLA_LOG_INFO_ON_ROOT( "viscosity =  " << viscosity );
+   WALBERLA_LOG_INFO_ON_ROOT( "omega =  " << omega );
+				
+   const real_t remainingTimeLoggerFrequency = parameters.getParameter< real_t >( "remainingTimeLoggerFrequency", real_c(3) ); // in seconds
+
+   // create fields
+   LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::SRT( omega ) );
+   BlockDataID pdfFieldID = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel, inletVelocity, real_t(1) );
+   BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" );
+
+   // create and initialize boundary handling
+   const FlagUID fluidFlagUID( "Fluid" );
+
+
+   using BHFactory = lbm::DefaultBoundaryHandlingFactory< LatticeModel_T, FlagField_T > ;
+
+   BlockDataID boundaryHandlingID = BHFactory::addBoundaryHandlingToStorage( blocks, "boundary handling", flagFieldID, pdfFieldID, fluidFlagUID,
+                                                                             inletVelocity,
+                                                                             boundariesConfig.getParameter< Vector3<real_t> >( "velocity1", Vector3<real_t>() ),
+                                                                             boundariesConfig.getParameter< real_t > ( "pressure0", real_c( 1.0 ) ),
+                                                                             boundariesConfig.getParameter< real_t > ( "pressure1", real_c( 1.0 ) ) );
+  //! [geomboundary]
+   geometry::initBoundaryHandling<BHFactory::BoundaryHandling>( *blocks, boundaryHandlingID, boundariesConfig );
+  //! [geomboundary]
+   geometry::setNonBoundaryCellsToDomain<BHFactory::BoundaryHandling> ( *blocks, boundaryHandlingID );
+
+   // create time loop
+   SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps );
+
+   // create communication for PdfField
+   blockforest::communication::UniformBufferedScheme< CommunicationStencil_T > communication( blocks );
+   communication.addPackInfo( make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >( pdfFieldID ) );
+
+   // add LBM sweep and communication to time loop
+   timeloop.add() << BeforeFunction( communication, "communication" )
+                  << Sweep( BHFactory::BoundaryHandling::getBlockSweep( boundaryHandlingID ), "boundary handling" );
+   timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, fluidFlagUID ) ), "LB stream & collide" );
+
+   //add ReattachmentLengthFinder logger
+
+   //!  [Logger]
+   std::string loggingFileName = "ReattachmentLengthLogging_Re" + std::to_string(uint_c(Re)) + ".txt";
+   uint_t loggingFrequency = parameters.getParameter< uint_t >( "loggingFrequency", uint_c(1) );
+   Vector3<real_t> stepSize = boundariesConfig.getOneBlock("Body").getParameter< Vector3<real_t> >("max", Vector3<real_t>() );
+
+   timeloop.addFuncAfterTimeStep( ReattachmentLengthFinder( blocks, pdfFieldID, flagFieldID, fluidFlagUID,
+                                                            loggingFileName, loggingFrequency, stepSize ),
+                                  "reattachment length finder" );
+   //!  [Logger]
+
+   // log remaining time
+   timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger( timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency ), "remaining time logger" );
+
+   // add VTK output to time loop
+   lbm::VTKOutput< LatticeModel_T, FlagField_T >::addToTimeloop( timeloop, blocks, walberlaEnv.config(), pdfFieldID, flagFieldID, fluidFlagUID );
+
+   timeloop.run();
+
+   return EXIT_SUCCESS;
+}
+}
+
+int main( int argc, char ** argv )
+{
+   walberla::main(argc, argv);
+}

apps/tutorials/lbm/05_BackwardFacingStep.dox

+namespace walberla {
+
+/**
+\page tutorial_lbm05 Tutorial - LBM 5:  Backward-facing step
+
+\brief A configurable application for simulation of backward-facing step problem
+
+\section tutorial05_overview Overview
+
+The aim of this tutorial is to show how to build and solve the backward-facing step model using lattice Boltzman method in waLBerla.
+The "01_BasicLBM" case is used as the foundation of the current work. Therefore, most of the functionalities have already been introduced and discussed in LBM 1 tutorial.
+Here the main focus is on the following areas:
+
+- parameterization with the help of the Reynolds number
+- setting up the step geometry with the help of the input file
+- evaluation of the recirculation length
+
+\image html tutorial_lbm05_BackwardFacingStep_Re400Result.png "Backward-facing step (implemented by a box shape in configuration file) in a 2D channel, Re=400"  width=1200px
+
+\section tutorial05_lbmdgeneralsetup General Setup
+
+A 2D channel flow with a backward-facing step is simulated.
+The application is configured with the Reynolds number and lattice velocity as input and reattachment location as the output.
+The input parameters are specified in a configuration file while the calculation and writing of reattachment locations are performed in a functor `ReattachmentLengthFinder` which is implemented inside the program.
+
+\section tutorial05_lbmdatastructures Lattice Boltzmann Model
+
+Since the simulations are carried out in 2D, the standard D2Q9 stencil with SRT collision model is used which is already implemented in the 'lbm' module.
+
+\snippet 05_BackwardFacingStep.cpp typedefs
+
+
+\section tutorial05_parameterization Parameterization with Reynolds number
+
+The **main()** function must be modified so that the program could accept the Reynolds number and lattice velocity as input parameters.
+These parameters are specified in the 'parameters' section of the config file.
+
+\snippet 05_BackwardFacingStep.prm Parameters
+
+The following lines inside the main() function reads these two values from the parameters section of the configuration file:
+
+\snippet 05_BackwardFacingStep.cpp Params
+
+The height of the channel is considered as the reference length in this simulation (it may differ in other applications).
+This value is read from the DomainSetup section in the configuration file where the dimensions of the channel are specified.
+
+\snippet 05_BackwardFacingStep.cpp height
+
+Finally, viscosity consequently **omega** are calculated with the Reynolds number, lattice velocity and reference length.
+
+\snippet 05_BackwardFacingStep.cpp Omega
+
+
+\section tutorial05_geometry Geometry
+
+Since the step geometry is a plain rectangular area, the simplest approach is to create it by geometry module in walberla.
+This module offers capability to read boundaries of a random geometry from images, voxel files, coordinates of verticies, etc.
+Using this module, obstacles of basic shapes could be conveniently positioned inside the domain.
+It is also easier to have the program to read the geometry specifications from the Boundaries section of the configuration file.
+This is implemented by reading and storing the Boundaries block of the configuration file in 'boundariesConfig' object and passing it to a convenience function provided in the geometry class to initialize the boundaries.
+
+\snippet 05_BackwardFacingStep.cpp geomboundary
+Here a subblock 'Body' is created inside 'Boundaries' section in the configuration file in order to create a box (rectangle in 2D) using two diagonal verticies.
+
+\snippet 05_BackwardFacingStep.prm geometry
+
+For more details about specifying boundaries using configuration file, please refer to the documentation of walberla::geometry::initBoundaryHandling().
+
+
+\section tutorial05_evaluation Evaluation of the recirculation length
+
+In order to find the reattachment location precisely, the velocity on each lattice cell on the bottom boundary line of the domain following the step is examined.
+The idea is to locate the exact position on the bottom boundary in which the streamwise velocity component switches sign.
+This mechanism is implemented by a functor named `ReattachmentLengthFinder` and is passed to a method of timeloop 'addFuncBeforeTimeStep()'' to create the logger.
+
+\snippet 05_BackwardFacingStep.cpp Logger
+
+After running the program, the locations of reattachment against timestep are written to 'ReattachmentLengthLogging_Re_[Re].txt' in the working directory.
+Note that there might be more than one reattachment location before the flow fully develops along the channel, and all are given in the file.
+This simply means that it is expected to have multiple occurances of seperation and reattachment at the same time along the bottom boundary of the channel following the step in the early stages.
+However, most of them are smeared later as the flow starts to develop.
+The logging frequency can also be adjusted by 'checkFrequency' which is passed to the `ReattachmentLengthFinder` functor.
+
+\section tutorial05_results Results
+
+The following results of the normalized reattachment length were obtained by a numerical study, using the here developed application.
+Note that it might take a long time until the simulation has reached a steady state and the final values for this length can be obtained.
+Also, the step has to be long enough to feature a converged, parabolic flow profile before the expansion, and the domain itself has to be long enough to not perturb the results by the vicinity of the outflow boundary condition.
+Reference results were taken from Ercan Erturk, “Numerical solutions of 2-D steady incompressible flow over a backward-facing step, Part I: High Reynolds number solutions”, Computers & Fluids, Vol. 37, 2008.
+
+<table>
+<caption id="multi_row">Normalized recirculation length for different Re</caption>
+<tr><th>Re     <th>simulation  <th>literature <th>% of difference
+<tr><td>   100 <td>  2.880 <td>  2.922 <td> 1.44
+<tr><td>   200 <td>  4.940 <td>  4.982 <td> 0.84
+<tr><td>   300 <td>  6.700 <td>  6.751 <td> 0.76
+<tr><td>   400 <td>  8.180 <td>  8.237 <td> 0.69
+<tr><td>   500 <td>  9.360 <td>  9.420 <td> 0.64
+<tr><td>   600 <td> 10.300 <td> 10.349 <td> 0.47
+<tr><td>   700 <td> 11.080 <td> 11.129 <td> 0.44
+<tr><td>   800 <td> 11.780 <td> 11.834 <td> 0.46
+<tr><td>   900 <td> 12.420 <td> 12.494 <td> 0.59
+<tr><td>  1000 <td> 13.040 <td> 13.121 <td> 0.62
+</table>
+
+\section tutorial05_outlook Outlook
+
+To improve this application, parallelization with a proper synchronization of the reattachment lengths could be added.
+Also, a stopping criterion based on the convergence of some relevant quantity could be added to avoid a manual termination of the simulation.
+
+\tableofcontents
+
+*/
+
+}// namespace walberla
\ No newline at end of file

apps/tutorials/lbm/05_BackwardFacingStep.prm

+
+//! [Parameters]
+Parameters 
+{
+    Re              400;
+    timesteps       10000000;
+    uLB             0.05;
+    remainingTimeLoggerFrequency 3; // in seconds
+    loggingFrequency 50; // in time steps
+}
+//! [Parameters]
+
+DomainSetup
+{
+   blocks        <  1,  1,  1 >;
+   cellsPerBlock <  6000, 100, 1 >;
+   periodic      <  0,   0,  1 >;  
+}
+
+Boundaries 
+{
+	velocity0 < 0.01, 0, 0 >; // velocity of cells where Velocity0 boundary is set
+	velocity1 <   0, 0, 0 >; // velocity of cells where Velocity1 boundary is set
+	pressure0 1.1;		     // pressure of cells where Pressure0 boundary is set
+	pressure1 1.0;			 // pressure of cells where Pressure1 boundary is set
+    
+    
+	Border { direction W;    walldistance -1;  Velocity0 {} }
+	Border { direction E;    walldistance -1;  Pressure1 {} }
+	Border { direction S,N;  walldistance -1;  NoSlip    {} }
+//! [geometry]
+    Body
+    {
+       shape box;
+       min <0,0,0>;
+       max <2000,50,1>;
+       NoSlip {} 
+    }
+//! [geometry]
+}
+
+VTK 
+{
+   // for parameter documentation see src/vtk/Initialization.cpp
+   fluid_field
+   {
+      writeFrequency   20000;
+      
+      writers {
+         Velocity;
+         Density;
+      }
+      
+      baseFolder	vtk_05_BackwardFacingStep;
+   }
+   
+   flag_field
+   {
+      writeFrequency 10000000; // write only once
+      
+      writers {
+          FlagField;
+      }
+      baseFolder	vtk_05_BackwardFacingStep;
+   }
+
+}

apps/tutorials/lbm/CMakeLists.txt

@@ -21,4 +21,8 @@ waLBerla_add_executable ( NAME 04_LBComplexGeometry
                           FILES 04_LBComplexGeometry.cpp
                           DEPENDS blockforest boundary core field lbm mesh stencil timeloop vtk )
                           
-endif()
+endif()
\ No newline at end of file
+
+waLBerla_add_executable ( NAME 05_BackwardFacingStep
+                          FILES 05_BackwardFacingStep.cpp
+                          DEPENDS blockforest boundary core field lbm stencil timeloop vtk )
\ No newline at end of file

doc/Mainpage.dox

@@ -42,6 +42,8 @@ all the basic data strcutures and concepts of the framework.
   A full LBM simulation is built.
 - \ref tutorial_lbm04 \n
   A LBM simulation with complex geometries is built.
+- \ref tutorial_lbm05 \n
+  A LBM simulation with a backward-facing step is set up.
 
 \subsection codegen Code Generation