diff --git a/apps/tutorials/lbm/05_BackwardFacingStep.cpp b/apps/tutorials/lbm/05_BackwardFacingStep.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e939f2ee2202d663d158b2ef1b2202da37f2c8ae
--- /dev/null
+++ b/apps/tutorials/lbm/05_BackwardFacingStep.cpp
@@ -0,0 +1,243 @@
+//======================================================================================================================
+//
+// 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);
+}
diff --git a/apps/tutorials/lbm/05_BackwardFacingStep.dox b/apps/tutorials/lbm/05_BackwardFacingStep.dox
new file mode 100644
index 0000000000000000000000000000000000000000..52f7160d5bc6d80f128c81f3f5762151625e5dab
--- /dev/null
+++ b/apps/tutorials/lbm/05_BackwardFacingStep.dox
@@ -0,0 +1,115 @@
+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
diff --git a/apps/tutorials/lbm/05_BackwardFacingStep.prm b/apps/tutorials/lbm/05_BackwardFacingStep.prm
new file mode 100644
index 0000000000000000000000000000000000000000..100099d60d2eb6782096979ad4f07814836190ee
--- /dev/null
+++ b/apps/tutorials/lbm/05_BackwardFacingStep.prm
@@ -0,0 +1,67 @@
+
+//! [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;
+ }
+
+}
diff --git a/apps/tutorials/lbm/CMakeLists.txt b/apps/tutorials/lbm/CMakeLists.txt
index 67d4df6877e3a103ca54280632c543849a8811e1..3f1183f2a59a40c16ccd41ca7ab7aba33f04ece9 100644
--- a/apps/tutorials/lbm/CMakeLists.txt
+++ b/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()
\ No newline at end of file
+endif()
+
+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
diff --git a/doc/Mainpage.dox b/doc/Mainpage.dox
index 0f35c3422b1e59f367c7339373bc4c57931c6a97..6402c2aae6a05ed234a36393dba85851a512216b 100644
--- a/doc/Mainpage.dox
+++ b/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
diff --git a/doc/pics/tutorial_lbm05_BackwardFacingStep_Re400Result.png b/doc/pics/tutorial_lbm05_BackwardFacingStep_Re400Result.png
new file mode 100644
index 0000000000000000000000000000000000000000..5bec5bdf30f65a559c87c36ca143be7619374045
Binary files /dev/null and b/doc/pics/tutorial_lbm05_BackwardFacingStep_Re400Result.png differ