From 0a5b7370db31d9f39f3595933b717d1fdc457896 Mon Sep 17 00:00:00 2001
From: Michael Kuron <mkuron@icp.uni-stuttgart.de>
Date: Tue, 26 Mar 2019 09:58:39 +0100
Subject: [PATCH] Added functionality to LBM boundaries to optionally obtain
the hydrodynamic force on those boundaries. Also added tests.
---
src/lbm/boundary/DynamicUBB.h | 28 ++-
src/lbm/boundary/NoSlip.h | 37 ++-
src/lbm/boundary/ParserUBB.h | 94 +++++---
src/lbm/boundary/SimpleUBB.h | 45 +++-
src/lbm/boundary/UBB.h | 70 ++++--
tests/lbm/CMakeLists.txt | 7 +
tests/lbm/boundary/BoundaryForce.cpp | 228 +++++++++++++++++++
tests/lbm/boundary/BoundaryForceCouette.cpp | 237 ++++++++++++++++++++
8 files changed, 685 insertions(+), 61 deletions(-)
create mode 100644 tests/lbm/boundary/BoundaryForce.cpp
create mode 100644 tests/lbm/boundary/BoundaryForceCouette.cpp
diff --git a/src/lbm/boundary/DynamicUBB.h b/src/lbm/boundary/DynamicUBB.h
index 67425aa05..945cc36f3 100644
--- a/src/lbm/boundary/DynamicUBB.h
+++ b/src/lbm/boundary/DynamicUBB.h
@@ -51,12 +51,14 @@ namespace lbm {
// 2. A member function "Vector3< real_t > operator()( const Vector3< real_t > & x, const real_t t )" that is called for every
// boundary link treated by "treatDirection". The arguments are the position 'x' of the boudnary cell in the simulation space and the current time 't'.
// The functon is supposed to return the velocity used by the boundary treatment.
-template< typename LatticeModel_T, typename flag_t, typename VelocityFunctor_T, bool AdaptVelocityToExternalForce = false >
+template< typename LatticeModel_T, typename flag_t, typename VelocityFunctor_T, bool AdaptVelocityToExternalForce = false, bool StoreForce = false >
class DynamicUBB : public Boundary<flag_t>
{
typedef lbm::PdfField< LatticeModel_T > PDFField;
typedef typename LatticeModel_T::Stencil Stencil;
+ typedef GhostLayerField< Vector3<real_t>, 1 > ForceField;
+
public:
static const bool threadsafe = true;
@@ -83,6 +85,9 @@ public:
{
velocity_( time_ );
}
+
+ if (StoreForce)
+ force_ = make_shared<ForceField>( pdfField_->xSize(), pdfField_->ySize(), pdfField_->zSize(), pdfField_->nrOfGhostLayers(), field::zyxf );
}
DynamicUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField,
const uint_t level, const VelocityFunctor_T & velocity, const AABB & aabb ) :
@@ -100,6 +105,9 @@ public:
time_ = timeTracker_->getTime( level_ );
velocity_( time_ );
}
+
+ if (StoreForce)
+ force_->setWithGhostLayer( Vector3<real_t>() );
}
void afterBoundaryTreatment() const {}
@@ -131,6 +139,8 @@ public:
WALBERLA_ASSERT_EQUAL( mask & this->mask_, this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
// current implementation of this boundary condition (DynamicUBB)
+ const real_t pdf_old = pdfField_->get( x, y, z, Stencil::idx[dir] );
+
const Vector3< real_t > pos( origin_[0] + real_c(nx) * dx_[0],
origin_[1] + real_c(ny) * dx_[1],
origin_[2] + real_c(nz) * dx_[2] );
@@ -160,6 +170,21 @@ public:
real_c(stencil::cy[ dir ]) * vel[1] +
real_c(stencil::cz[ dir ]) * vel[2] ) );
}
+
+ if (StoreForce && pdfField_->isInInnerPart( Cell(x,y,z) ))
+ {
+ const real_t forceMEM = pdf_old + pdfField_->get( nx, ny, nz, Stencil::invDirIdx(dir) );
+ Vector3<real_t> force( real_c( stencil::cx[dir] ) * forceMEM,
+ real_c( stencil::cy[dir] ) * forceMEM,
+ real_c( stencil::cz[dir] ) * forceMEM );
+ force_->get( nx, ny, nz ) += force;
+ }
+ }
+
+ const typename ForceField::value_type & getForce( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const
+ {
+ static_assert(StoreForce, "this member function is only available if the fourth template argument on the class is true");
+ return force_->get(x,y,z);
}
private:
@@ -177,6 +202,7 @@ private:
uint_t level_;
VelocityFunctor_T velocity_;
+ shared_ptr<ForceField> force_;
}; // class DynamicUBB
diff --git a/src/lbm/boundary/NoSlip.h b/src/lbm/boundary/NoSlip.h
index 0b9d3e32a..dc769d345 100644
--- a/src/lbm/boundary/NoSlip.h
+++ b/src/lbm/boundary/NoSlip.h
@@ -42,7 +42,7 @@ namespace lbm {
-template< typename LatticeModel_T, typename flag_t >
+template< typename LatticeModel_T, typename flag_t, bool StoreForce = false >
class NoSlip : public Boundary<flag_t>
{
protected:
@@ -50,6 +50,8 @@ protected:
typedef PdfField< LatticeModel_T > PDFField;
typedef typename LatticeModel_T::Stencil Stencil;
+ typedef GhostLayerField< Vector3<real_t>, 1 > ForceField;
+
public:
static const bool threadsafe = true;
@@ -60,11 +62,21 @@ public:
NoSlip( const BoundaryUID& boundaryUID, const FlagUID& uid, PDFField* const pdfField ) :
- Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ) { WALBERLA_ASSERT_NOT_NULLPTR( pdfField_ ); }
+ Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField )
+ {
+ WALBERLA_ASSERT_NOT_NULLPTR( pdfField_ );
+
+ if (StoreForce)
+ force_ = make_shared<ForceField>( pdfField_->xSize(), pdfField_->ySize(), pdfField_->zSize(), pdfField_->nrOfGhostLayers(), field::zyxf );
+ }
void pushFlags( std::vector< FlagUID >& uids ) const { uids.push_back( uid_ ); }
- void beforeBoundaryTreatment() const {}
+ void beforeBoundaryTreatment() const
+ {
+ if (StoreForce)
+ force_->setWithGhostLayer( Vector3<real_t>() );
+ }
void afterBoundaryTreatment() const {}
template< typename Buffer_T >
@@ -95,7 +107,24 @@ public:
WALBERLA_ASSERT_EQUAL( mask & this->mask_, this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
// current implementation of this boundary condition (NoSlip)
+ const real_t pdf_old = pdfField_->get( x, y, z, Stencil::idx[dir] );
+
pdfField_->get( nx, ny, nz, Stencil::invDirIdx(dir) ) = pdfField_->get( x, y, z, Stencil::idx[dir] );
+
+ if (StoreForce && pdfField_->isInInnerPart( Cell(x,y,z) ))
+ {
+ const real_t forceMEM = pdf_old + pdfField_->get( nx, ny, nz, Stencil::invDirIdx(dir) );
+ Vector3<real_t> force( real_c( stencil::cx[dir] ) * forceMEM,
+ real_c( stencil::cy[dir] ) * forceMEM,
+ real_c( stencil::cz[dir] ) * forceMEM );
+ force_->get( nx, ny, nz ) += force;
+ }
+ }
+
+ const typename ForceField::value_type & getForce( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const
+ {
+ static_assert(StoreForce, "this member function is only available if the third template argument on the class is true");
+ return force_->get(x,y,z);
}
private:
@@ -104,6 +133,8 @@ private:
PDFField* const pdfField_;
+ shared_ptr<ForceField> force_;
+
}; // class NoSlip
diff --git a/src/lbm/boundary/ParserUBB.h b/src/lbm/boundary/ParserUBB.h
index 82bfff5e2..415f36bc4 100644
--- a/src/lbm/boundary/ParserUBB.h
+++ b/src/lbm/boundary/ParserUBB.h
@@ -47,12 +47,14 @@ namespace lbm {
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce = false >
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce = false, bool StoreForce = false >
class ParserUBB : public Boundary<flag_t>
{
typedef lbm::PdfField< LatticeModel_T > PDFField;
typedef typename LatticeModel_T::Stencil Stencil;
+ typedef GhostLayerField< Vector3<real_t>, 1 > ForceField;
+
public:
static const bool threadsafe = true;
@@ -118,6 +120,9 @@ public:
{
if( timeTracker_ )
time_ = timeTracker_->getTime( level_ );
+
+ if (StoreForce)
+ force_->setWithGhostLayer( Vector3<real_t>() );
}
void afterBoundaryTreatment() const {}
@@ -145,6 +150,12 @@ public:
inline Vector3<real_t> getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const;
inline Vector3<real_t> getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z, real_t t ) const;
+ const typename ForceField::value_type & getForce( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const
+ {
+ static_assert(StoreForce, "this member function is only available if the fourth template argument on the class is true");
+ return force_->get(x,y,z);
+ }
+
private:
const FlagUID uid_;
@@ -161,11 +172,12 @@ private:
shared_ptr<ParserField> parserField_;
shared_ptr<VelocityField> velocityField_;
+ shared_ptr<ForceField> force_;
}; // class ParserUBB
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
-inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::Parser( const Config::BlockHandle & config )
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
+inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce>::Parser::Parser( const Config::BlockHandle & config )
: parsers_(), equations_(), timeDependent_( false )
{
if( !config )
@@ -194,8 +206,8 @@ inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::
}
}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
-inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::Parser( std::array< std::string, 3 > & equations )
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
+inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce>::Parser::Parser( std::array< std::string, 3 > & equations )
: parsers_(), equations_( equations ), timeDependent_( false )
{
if( equations_[0].length() > 0 )
@@ -218,16 +230,16 @@ inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::
}
}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
-inline ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::Velocity::Velocity( const Config::BlockHandle & config )
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
+inline ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::Velocity::Velocity( const Config::BlockHandle & config )
{
velocity_[0] = ( config && config.isDefined( "x" ) ) ? config.getParameter<real_t>( "x" ) : real_c(0.0);
velocity_[1] = ( config && config.isDefined( "y" ) ) ? config.getParameter<real_t>( "y" ) : real_c(0.0);
velocity_[2] = ( config && config.isDefined( "z" ) ) ? config.getParameter<real_t>( "z" ) : real_c(0.0);
}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
-Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::operator()( const Vector3< real_t > & x, const real_t t ) const
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
+Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce>::Parser::operator()( const Vector3< real_t > & x, const real_t t ) const
{
std::map< std::string, double > symbols;
symbols["x"] = x[0];
@@ -242,8 +254,8 @@ Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce
return v;
}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
-Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::operator()( const Vector3< real_t > & x ) const
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
+Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce>::Parser::operator()( const Vector3< real_t > & x ) const
{
WALBERLA_ASSERT( !timeDependent_ );
@@ -261,8 +273,8 @@ Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
-inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::ParserUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField,
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
+inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce>::ParserUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField,
FlagField<flag_t> * const flagField, const shared_ptr< TimeTracker > & timeTracker,
const uint_t level, const AABB & aabb )
: Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ), timeTracker_( timeTracker ), time_( real_t(0) ), level_( level )
@@ -285,13 +297,16 @@ inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::ParserUB
parserField_ = make_shared<ParserField> ( pdfField->xSize(), pdfField->ySize(), pdfField->zSize(), pdfField->nrOfGhostLayers(), field::zyxf );
velocityField_ = make_shared<VelocityField>( pdfField->xSize(), pdfField->ySize(), pdfField->zSize(), pdfField->nrOfGhostLayers(), field::zyxf );
}
+
+ if (StoreForce)
+ force_ = make_shared<ForceField>( pdfField_->xSize(), pdfField_->ySize(), pdfField_->zSize(), pdfField_->nrOfGhostLayers(), field::zyxf );
}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
template< typename Buffer_T >
-inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::packCell( Buffer_T & buffer, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::packCell( Buffer_T & buffer, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
{
if( parserField_->get( x, y, z ) )
{
@@ -306,9 +321,9 @@ inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::p
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
template< typename Buffer_T >
-inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCell( Buffer_T & buffer, const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::registerCell( Buffer_T & buffer, const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
{
bool isparser;
buffer >> isparser;
@@ -329,16 +344,16 @@ inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::r
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
-inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::ParserUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField,
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
+inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce>::ParserUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField,
FlagField<flag_t> * const flagField, const uint_t level, const AABB & aabb )
: ParserUBB( boundaryUID, uid, pdfField, flagField, nullptr, level, aabb )
{}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
-inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
const BoundaryConfiguration & parser )
{
WALBERLA_ASSERT_EQUAL( dynamic_cast< const Parser * >( &parser ), &parser );
@@ -369,8 +384,8 @@ inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::r
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
-inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & parser )
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & parser )
{
WALBERLA_ASSERT_EQUAL( dynamic_cast< const Parser * >( &parser ), &parser );
WALBERLA_ASSERT_NOT_NULLPTR( parserField_ );
@@ -408,9 +423,9 @@ inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::r
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
template< typename CellIterator >
-inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCells( const flag_t, const CellIterator & begin, const CellIterator & end,
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::registerCells( const flag_t, const CellIterator & begin, const CellIterator & end,
const BoundaryConfiguration & parser )
{
WALBERLA_ASSERT_EQUAL( dynamic_cast< const Parser * >( &parser ), &parser );
@@ -451,20 +466,20 @@ inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::r
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
-inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::unregisterCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::unregisterCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
{
parserField_->get(x,y,z) = nullptr;
}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
#ifndef NDEBUG
-inline void ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::ParserUBB::treatDirection( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const stencil::Direction dir,
+inline void ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce>::ParserUBB::treatDirection( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const stencil::Direction dir,
const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t mask )
#else
-inline void ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::ParserUBB::treatDirection( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const stencil::Direction dir,
+inline void ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce>::ParserUBB::treatDirection( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const stencil::Direction dir,
const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t /*mask*/ )
#endif
{
@@ -475,6 +490,8 @@ inline void ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Par
WALBERLA_ASSERT_EQUAL( mask & this->mask_, this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
// current implementation of this boundary condition (ParserUBB)
+ const real_t pdf_old = pdfField_->get( x, y, z, Stencil::idx[dir] );
+
Vector3<real_t> velocity;
if( parserField_->get(nx,ny,nz) )
{
@@ -512,18 +529,27 @@ inline void ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Par
real_c(stencil::cy[ dir ]) * vel[1] +
real_c(stencil::cz[ dir ]) * vel[2] ) );
}
+
+ if (StoreForce && pdfField_->isInInnerPart( Cell(x,y,z) ))
+ {
+ const real_t forceMEM = pdf_old + pdfField_->get( nx, ny, nz, Stencil::invDirIdx(dir) );
+ Vector3<real_t> force( real_c( stencil::cx[dir] ) * forceMEM,
+ real_c( stencil::cy[dir] ) * forceMEM,
+ real_c( stencil::cz[dir] ) * forceMEM );
+ force_->get( nx, ny, nz ) += force;
+ }
}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
-inline Vector3<real_t> ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
+inline Vector3<real_t> ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const
{
return getValue( x, y, z, time_ );
}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
-inline Vector3<real_t> ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z, real_t t ) const
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
+inline Vector3<real_t> ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z, real_t t ) const
{
Vector3<real_t> velocity;
if( parserField_->get(x,y,z) )
diff --git a/src/lbm/boundary/SimpleUBB.h b/src/lbm/boundary/SimpleUBB.h
index da6e7839b..ba07eedc4 100644
--- a/src/lbm/boundary/SimpleUBB.h
+++ b/src/lbm/boundary/SimpleUBB.h
@@ -45,12 +45,14 @@ namespace lbm {
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce = false >
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce = false, bool StoreForce = false >
class SimpleUBB : public Boundary<flag_t>
{
typedef PdfField< LatticeModel_T > PDFField;
typedef typename LatticeModel_T::Stencil Stencil;
+ typedef GhostLayerField< Vector3<real_t>, 1 > ForceField;
+
public:
static const bool threadsafe = true;
@@ -61,14 +63,30 @@ public:
SimpleUBB( const BoundaryUID& boundaryUID, const FlagUID& uid, PDFField* const pdfField, const Vector3< real_t > & velocity ) :
- Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ), velocity_( velocity ) { WALBERLA_ASSERT_NOT_NULLPTR( pdfField_ ); }
+ Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ), velocity_( velocity )
+{
+ WALBERLA_ASSERT_NOT_NULLPTR( pdfField_ );
+
+ if (StoreForce)
+ force_ = make_shared<ForceField>( pdfField_->xSize(), pdfField_->ySize(), pdfField_->zSize(), pdfField_->nrOfGhostLayers(), field::zyxf );
+}
SimpleUBB( const BoundaryUID& boundaryUID, const FlagUID& uid, PDFField* const pdfField, const real_t x, const real_t y, const real_t z ) :
- Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ), velocity_( x, y, z ) { WALBERLA_ASSERT_NOT_NULLPTR( pdfField_ ); }
+ Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ), velocity_( x, y, z )
+{
+ WALBERLA_ASSERT_NOT_NULLPTR( pdfField_ );
+
+ if (StoreForce)
+ force_ = make_shared<ForceField>( pdfField_->xSize(), pdfField_->ySize(), pdfField_->zSize(), pdfField_->nrOfGhostLayers(), field::zyxf );
+}
void pushFlags( std::vector< FlagUID >& uids ) const { uids.push_back( uid_ ); }
- void beforeBoundaryTreatment() const {}
+ void beforeBoundaryTreatment() const
+ {
+ if (StoreForce)
+ force_->setWithGhostLayer( Vector3<real_t>() );
+ }
void afterBoundaryTreatment() const {}
template< typename Buffer_T >
@@ -99,6 +117,8 @@ public:
WALBERLA_ASSERT_EQUAL( mask & this->mask_, this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
// current implementation of this boundary condition (SimpleUBB)
+ const real_t pdf_old = pdfField_->get( x, y, z, Stencil::idx[dir] );
+
if( LatticeModel_T::compressible )
{
const auto density = pdfField_->getDensity(x,y,z);
@@ -122,6 +142,21 @@ public:
real_c(stencil::cy[ dir ]) * velocity[1] +
real_c(stencil::cz[ dir ]) * velocity[2] ) );
}
+
+ if (StoreForce && pdfField_->isInInnerPart( Cell(x,y,z) ))
+ {
+ const real_t forceMEM = pdf_old + pdfField_->get( nx, ny, nz, Stencil::invDirIdx(dir) );
+ Vector3<real_t> force( real_c( stencil::cx[dir] ) * forceMEM,
+ real_c( stencil::cy[dir] ) * forceMEM,
+ real_c( stencil::cz[dir] ) * forceMEM );
+ force_->get( nx, ny, nz ) += force;
+ }
+ }
+
+ const typename ForceField::value_type & getForce( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const
+ {
+ static_assert(StoreForce, "this member function is only available if the fourth template argument on the class is true");
+ return force_->get(x,y,z);
}
private:
@@ -132,6 +167,8 @@ private:
const Vector3< real_t > velocity_;
+ shared_ptr<ForceField> force_;
+
}; // class SimpleUBB
diff --git a/src/lbm/boundary/UBB.h b/src/lbm/boundary/UBB.h
index 361322414..5ce5fcc06 100644
--- a/src/lbm/boundary/UBB.h
+++ b/src/lbm/boundary/UBB.h
@@ -46,7 +46,7 @@ namespace lbm {
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce = false >
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce = false, bool StoreForce = false >
class UBB : public Boundary<flag_t>
{
typedef PdfField< LatticeModel_T > PDFField;
@@ -54,6 +54,8 @@ class UBB : public Boundary<flag_t>
typedef GhostLayerField< Vector3<real_t>, 1 > VelField;
+ typedef GhostLayerField< Vector3<real_t>, 1 > ForceField;
+
public:
static const bool threadsafe = true;
@@ -89,7 +91,7 @@ public:
void pushFlags( std::vector< FlagUID > & uids ) const { uids.push_back( uid_ ); }
- void beforeBoundaryTreatment() const {}
+ inline void beforeBoundaryTreatment() const;
void afterBoundaryTreatment() const {}
template< typename Buffer_T >
@@ -110,19 +112,26 @@ public:
inline const Vector3<real_t> & getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const { return vel_->get(x,y,z); }
+ const typename ForceField::value_type & getForce( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const
+ {
+ static_assert(StoreForce, "this member function is only available if the fourth template argument on the class is true");
+ return force_->get(x,y,z);
+ }
+
private:
const FlagUID uid_;
PDFField* const pdfField_;
shared_ptr<VelField> vel_;
+ shared_ptr<ForceField> force_;
}; // class UBB
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
-inline UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::Velocity::Velocity( const Config::BlockHandle & config )
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
+inline UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::Velocity::Velocity( const Config::BlockHandle & config )
{
velocity_[0] = ( config && config.isDefined( "x" ) ) ? config.getParameter<real_t>( "x" ) : real_c(0.0);
velocity_[1] = ( config && config.isDefined( "y" ) ) ? config.getParameter<real_t>( "y" ) : real_c(0.0);
@@ -131,8 +140,8 @@ inline UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::Velocity::Ve
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
-inline UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::UBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField* const pdfField, FlagField<flag_t> * const flagField ) :
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
+inline UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::UBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField* const pdfField, FlagField<flag_t> * const flagField ) :
Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField )
{
@@ -141,30 +150,42 @@ inline UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::UBB( const B
vel_ = make_shared<VelField>( pdfField_->xSize(), pdfField_->ySize(), pdfField_->zSize(), flagField->nrOfGhostLayers(), field::zyxf );
else
vel_ = make_shared<VelField>( pdfField_->xSize(), pdfField_->ySize(), pdfField_->zSize(), pdfField_->nrOfGhostLayers(), field::zyxf );
+
+ if (StoreForce)
+ force_ = make_shared<ForceField>( pdfField_->xSize(), pdfField_->ySize(), pdfField_->zSize(), pdfField_->nrOfGhostLayers(), field::zyxf );
+}
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
+inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::beforeBoundaryTreatment() const
+{
+ if (StoreForce)
+ force_->setWithGhostLayer( Vector3<real_t>() );
}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
template< typename Buffer_T >
-inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::packCell( Buffer_T & buffer, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::packCell( Buffer_T & buffer, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
{
buffer << vel_->get( x, y, z );
}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
template< typename Buffer_T >
-inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCell( Buffer_T & buffer, const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::registerCell( Buffer_T & buffer, const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
{
buffer >> vel_->get( x, y, z );
}
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
-inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
+inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
const BoundaryConfiguration & velocity )
{
WALBERLA_ASSERT_EQUAL( dynamic_cast< const Velocity * >( &velocity ), &velocity );
@@ -177,8 +198,8 @@ inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registe
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
-inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & velocity )
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
+inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & velocity )
{
WALBERLA_ASSERT_EQUAL( dynamic_cast< const Velocity * >( &velocity ), &velocity );
WALBERLA_ASSERT_NOT_NULLPTR( vel_ );
@@ -191,9 +212,9 @@ inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registe
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
template< typename CellIterator >
-inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCells( const flag_t, const CellIterator & begin, const CellIterator & end,
+inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::registerCells( const flag_t, const CellIterator & begin, const CellIterator & end,
const BoundaryConfiguration & velocity )
{
WALBERLA_ASSERT_EQUAL( dynamic_cast< const Velocity * >( &velocity ), &velocity );
@@ -207,12 +228,12 @@ inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registe
-template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
#ifndef NDEBUG
-inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::treatDirection( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const stencil::Direction dir,
+inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::treatDirection( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const stencil::Direction dir,
const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t mask )
#else
-inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::treatDirection( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const stencil::Direction dir,
+inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce >::treatDirection( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const stencil::Direction dir,
const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t /*mask*/ )
#endif
{
@@ -223,6 +244,8 @@ inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::treatDi
WALBERLA_ASSERT_EQUAL( mask & this->mask_, this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
// current implementation of this boundary condition (UBB)
+ const real_t pdf_old = pdfField_->get( x, y, z, Stencil::idx[dir] );
+
if( LatticeModel_T::compressible )
{
const auto density = pdfField_->getDensity(x,y,z);
@@ -246,6 +269,15 @@ inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::treatDi
real_c(stencil::cy[ dir ]) * velocity[1] +
real_c(stencil::cz[ dir ]) * velocity[2] ) );
}
+
+ if (StoreForce && pdfField_->isInInnerPart( Cell(x,y,z) ))
+ {
+ const real_t forceMEM = pdf_old + pdfField_->get( nx, ny, nz, Stencil::invDirIdx(dir) );
+ Vector3<real_t> force( real_c( stencil::cx[dir] ) * forceMEM,
+ real_c( stencil::cy[dir] ) * forceMEM,
+ real_c( stencil::cz[dir] ) * forceMEM );
+ force_->get( nx, ny, nz ) += force;
+ }
}
diff --git a/tests/lbm/CMakeLists.txt b/tests/lbm/CMakeLists.txt
index 7fef30af0..dd9496fce 100644
--- a/tests/lbm/CMakeLists.txt
+++ b/tests/lbm/CMakeLists.txt
@@ -23,6 +23,13 @@ waLBerla_execute_test( NAME SimpleDiffusionDirichletDef COMMAND $<TARGET_FILE:Si
waLBerla_execute_test( NAME SimpleDiffusionDirichletRef COMMAND $<TARGET_FILE:SimpleDiffusionDirichlet> -l 4 -w 4 -o 1.2 -d 3 -t 3 -c 0 -r 1 )
waLBerla_execute_test( NAME SimpleDiffusionDirichletCav COMMAND $<TARGET_FILE:SimpleDiffusionDirichlet> -l 5 -w 3 -o 1.2 -d 3 -t 333 -c 1 )
+waLBerla_compile_test( FILES boundary/BoundaryForce.cpp DEPENDS blockforest timeloop )
+waLBerla_execute_test( NAME BoundaryForceShort COMMAND $<TARGET_FILE:BoundaryForce> 10 PROCESSES 2 CONFIGURATIONS Debug DebugOptimized )
+waLBerla_execute_test( NAME BoundaryForce COMMAND COMMAND $<TARGET_FILE:BoundaryForce> PROCESSES 2 CONFIGURATIONS Release RelWithDbgInfo )
+
+waLBerla_compile_test( FILES boundary/BoundaryForceCouette.cpp DEPENDS blockforest timeloop )
+waLBerla_execute_test( NAME BoundaryForceCouetteShort COMMAND $<TARGET_FILE:BoundaryForceCouette> 10 PROCESSES 2 CONFIGURATIONS Debug DebugOptimized )
+waLBerla_execute_test( NAME BoundaryForceCouette COMMAND COMMAND $<TARGET_FILE:BoundaryForceCouette> PROCESSES 2 CONFIGURATIONS Release RelWithDbgInfo )
waLBerla_compile_test( FILES boundary/DiffusionDirichlet.cpp DEPENDS field blockforest geometry timeloop vtk )
waLBerla_execute_test( NAME DiffusionDirichletTest1 COMMAND $<TARGET_FILE:DiffusionDirichlet> -l 16 -w 1 -o 1.2 -v 0.1 -e 0.097711 -t 50 )
diff --git a/tests/lbm/boundary/BoundaryForce.cpp b/tests/lbm/boundary/BoundaryForce.cpp
new file mode 100644
index 000000000..08b41401c
--- /dev/null
+++ b/tests/lbm/boundary/BoundaryForce.cpp
@@ -0,0 +1,228 @@
+//======================================================================================================================
+//
+// 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 BoundaryForce.cpp
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//! \brief Calculates the drag force on a spherical no-slip boundary in a flow
+//! due to a velocity boundary condition and compares with Stokes' law.
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "boundary/BoundaryHandling.h"
+
+#include "core/mpi/Environment.h"
+#include "core/mpi/MPIManager.h"
+#include "core/mpi/Reduce.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "domain_decomposition/SharedSweep.h"
+
+#include "field/AddToStorage.h"
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+
+#include "geometry/initializer/BoundaryFromBody.h"
+#include "geometry/initializer/BoundaryFromDomainBorder.h"
+#include "geometry/InitBoundaryHandling.h"
+
+#include "lbm/boundary/SimpleUBB.h"
+#include "lbm/communication/PdfFieldPackInfo.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/D3Q19.h"
+#include "lbm/sweeps/CellwiseSweep.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+using namespace walberla;
+
+
+
+//////////////
+// TYPEDEFS //
+//////////////
+
+using LatticeModel_T = lbm::D3Q19<lbm::collision_model::TRT>;
+using Stencil_T = LatticeModel_T::Stencil;
+using CommunicationStencil_T = LatticeModel_T::CommunicationStencil;
+using PdfField_T = lbm::PdfField<LatticeModel_T>;
+
+using flag_t = walberla::uint8_t;
+using FlagField_T = FlagField<flag_t>;
+
+typedef lbm::SimpleUBB< LatticeModel_T, flag_t, false, true > UBB_Sphere_T;
+typedef lbm::SimpleUBB< LatticeModel_T, flag_t > UBB_Wall_T;
+typedef boost::tuples::tuple< UBB_Sphere_T, UBB_Wall_T > BoundaryConditions_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+
+
+
+///////////
+// FLAGS //
+///////////
+
+const FlagUID Fluid_Flag( "fluid" );
+const FlagUID UBB_Sphere_Flag( "sphere");
+const FlagUID UBB_Wall_Flag( "wall" );
+
+
+
+///////////////////////
+// BOUNDARY HANDLING //
+///////////////////////
+
+class MyBoundaryHandling
+{
+public:
+
+ MyBoundaryHandling( const BlockDataID & flagFieldId, const BlockDataID & pdfFieldId, const Vector3<real_t> & velocity ) :
+ flagFieldId_( flagFieldId ), pdfFieldId_( pdfFieldId ), velocity_( velocity ) {}
+
+ BoundaryHandling_T * operator()( IBlock* const block ) const;
+
+private:
+
+ const BlockDataID flagFieldId_;
+ const BlockDataID pdfFieldId_;
+
+ const Vector3<real_t> velocity_;
+
+}; // class MyBoundaryHandling
+
+BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block ) const
+{
+ FlagField_T * flagField = block->getData< FlagField_T >( flagFieldId_ );
+ PdfField_T * pdfField = block->getData< PdfField_T > ( pdfFieldId_ );
+
+ const flag_t fluid = flagField->registerFlag( Fluid_Flag ); // register the fluid flag at the flag field
+
+ return new BoundaryHandling_T( "boundary handling", flagField, fluid,
+ boost::tuples::make_tuple( UBB_Sphere_T( "sphere", UBB_Sphere_Flag, pdfField, Vector3<real_t>() ),
+ UBB_Wall_T( "wall", UBB_Wall_Flag, pdfField, velocity_ ) ) );
+}
+
+
+
+//////////
+// MAIN //
+//////////
+
+int main( int argc, char ** argv )
+{
+ mpi::Environment env( argc, argv );
+
+ int processes( MPIManager::instance()->numProcesses() );
+ if( processes != 1 && processes != 2 && processes != 4 && processes != 8 )
+ WALBERLA_ABORT( "The number of processes must be equal to either 1, 2, 4, or 8!" );
+
+ const uint_t xBlocks = ( processes == 8 ) ? uint_c(2) : uint_c(1);
+ const uint_t yBlocks = ( processes >= 4 ) ? uint_c(2) : uint_c(1);
+ const uint_t zBlocks = ( processes >= 2 ) ? uint_c(2) : uint_c(1);
+
+ const uint_t L = uint_t(64);
+
+ auto blocks = blockforest::createUniformBlockGrid( xBlocks, yBlocks, zBlocks,
+ L / xBlocks,
+ L / yBlocks,
+ L / zBlocks,
+ real_c(1),
+ true,
+ false, false, true);
+
+ const real_t omega = real_c(0.06);
+ LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ) );
+
+ const Vector3<real_t> velocity(real_t(0), real_t(0), real_t(0.01));
+ BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel );
+ BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" );
+ BlockDataID boundaryHandlingId = blocks->addBlockData< BoundaryHandling_T >( MyBoundaryHandling( flagFieldId, pdfFieldId, velocity ),
+ "boundary handling" );
+
+ // set boundary conditions
+ geometry::initializer::BoundaryFromBody< BoundaryHandling_T > bodyInitializer( *blocks, boundaryHandlingId );
+ const real_t R = real_t(5.4);
+ geometry::Sphere sphere( Vector3<real_t>( real_t(L/2), real_t(L/2), real_t(L/2) ), R );
+ bodyInitializer.template init< geometry::Sphere >( sphere, UBB_Sphere_Flag );
+ geometry::initializer::BoundaryFromDomainBorder< BoundaryHandling_T > borderInitializer( *blocks, boundaryHandlingId );
+ borderInitializer.init( UBB_Wall_Flag, stencil::N, cell_idx_t(-1) );
+ borderInitializer.init( UBB_Wall_Flag, stencil::S, cell_idx_t(-1) );
+ borderInitializer.init( UBB_Wall_Flag, stencil::E, cell_idx_t(-1) );
+ borderInitializer.init( UBB_Wall_Flag, stencil::W, cell_idx_t(-1) );
+ geometry::setNonBoundaryCellsToDomain< BoundaryHandling_T >( *blocks, boundaryHandlingId );
+
+ uint_t timeSteps = uint_c(200);
+ if( argc > 1 )
+ timeSteps = boost::lexical_cast<uint_t>( argv[1] );
+ SweepTimeloop timeloop( blocks->getBlockStorage(), timeSteps );
+
+ blockforest::communication::UniformBufferedScheme< CommunicationStencil_T > communication( blocks );
+ communication.addPackInfo( make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >( pdfFieldId ) );
+
+ timeloop.add() << BeforeFunction( communication, "communication" )
+ << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingId ), "boundary handling" );
+
+ timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldId, flagFieldId, Fluid_Flag ) ), "LB stream & collide" );
+
+ timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
+
+ timeloop.run();
+
+ // get the force acting on the sphere
+ real_t fx = real_t(0), fy = real_t(0), fz = real_t(0);
+ for( auto block = blocks->begin(); block != blocks->end(); ++block )
+ {
+ FlagField_T * flagField = block->getData< FlagField_T >( flagFieldId );
+ BoundaryHandling_T * boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingId );
+ auto buid = boundaryHandling->getBoundaryUID( UBB_Sphere_Flag );
+ auto & bc = boundaryHandling->getBoundaryCondition<UBB_Sphere_T>( buid );
+ auto flag = flagField->getFlag(UBB_Sphere_Flag);
+
+ // need to iterate over ghost cells too because they may contain part of the force
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ_OMP(flagField, omp parallel for schedule(static) reduction(+:fx) reduction(+:fy) reduction(+:fz),
+ {
+ if( flagField->isFlagSet(x, y, z, flag))
+ {
+ const Vector3<real_t> & cellForce = bc.getForce(x, y, z);
+ fx += cellForce[0];
+ fy += cellForce[1];
+ fz += cellForce[2];
+ }
+ });
+ }
+ Vector3<real_t> force(fx, fy, fz);
+
+ mpi::allReduceInplace( force[0], mpi::SUM );
+ mpi::allReduceInplace( force[1], mpi::SUM );
+ mpi::allReduceInplace( force[2], mpi::SUM );
+
+ real_t visc = lbm::collision_model::viscosityFromOmega( omega );
+ Vector3<real_t> stokes = 6 * math::PI * visc * R * velocity;
+
+ WALBERLA_LOG_RESULT_ON_ROOT("Expected force: " << stokes);
+ WALBERLA_LOG_RESULT_ON_ROOT("Actual force: " << force);
+ real_t err = (force-stokes).length() / stokes.length();
+ WALBERLA_LOG_RESULT_ON_ROOT("Relative error: " << (real_t(100)*err) << " %");
+
+ if( timeSteps > 100 )
+ {
+ WALBERLA_CHECK_LESS(err, 0.1);
+ WALBERLA_LOG_PROGRESS_ON_ROOT("Test succeeded");
+ }
+
+ return EXIT_SUCCESS;
+}
diff --git a/tests/lbm/boundary/BoundaryForceCouette.cpp b/tests/lbm/boundary/BoundaryForceCouette.cpp
new file mode 100644
index 000000000..1acf8f8c7
--- /dev/null
+++ b/tests/lbm/boundary/BoundaryForceCouette.cpp
@@ -0,0 +1,237 @@
+//======================================================================================================================
+//
+// 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 BoundaryForceCouette.cpp
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \brief Calculates the force on top and bottom boundary in plane Couette flow
+//! and checks correctness against analytical force
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "boundary/BoundaryHandling.h"
+
+#include "core/mpi/Environment.h"
+#include "core/mpi/MPIManager.h"
+#include "core/mpi/Reduce.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "domain_decomposition/SharedSweep.h"
+
+#include "field/AddToStorage.h"
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+
+#include "geometry/initializer/BoundaryFromDomainBorder.h"
+#include "geometry/InitBoundaryHandling.h"
+
+#include "lbm/boundary/NoSlip.h"
+#include "lbm/boundary/SimpleUBB.h"
+#include "lbm/communication/PdfFieldPackInfo.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/D3Q19.h"
+#include "lbm/sweeps/CellwiseSweep.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+using namespace walberla;
+
+
+
+//////////////
+// TYPEDEFS //
+//////////////
+
+using LatticeModel_T = lbm::D3Q19<lbm::collision_model::TRT>;
+using Stencil_T = LatticeModel_T::Stencil;
+using CommunicationStencil_T = LatticeModel_T::CommunicationStencil;
+using PdfField_T = lbm::PdfField<LatticeModel_T>;
+
+using flag_t = walberla::uint8_t;
+using FlagField_T = FlagField<flag_t>;
+
+typedef lbm::NoSlip< LatticeModel_T, flag_t, true > BottomWall_T;
+typedef lbm::SimpleUBB< LatticeModel_T, flag_t, false, true > TopWall_T;
+typedef boost::tuples::tuple< BottomWall_T, TopWall_T > BoundaryConditions_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+
+
+
+///////////
+// FLAGS //
+///////////
+
+const FlagUID Fluid_Flag( "fluid" );
+const FlagUID BottomWall_Flag( "bottom wall" );
+const FlagUID TopWall_Flag( "top wall" );
+
+
+
+///////////////////////
+// BOUNDARY HANDLING //
+///////////////////////
+
+class MyBoundaryHandling
+{
+public:
+
+ MyBoundaryHandling( const BlockDataID & flagFieldId, const BlockDataID & pdfFieldId, const Vector3<real_t> & velocity ) :
+ flagFieldId_( flagFieldId ), pdfFieldId_( pdfFieldId ), velocity_( velocity ) {}
+
+ BoundaryHandling_T * operator()( IBlock* const block ) const;
+
+private:
+
+ const BlockDataID flagFieldId_;
+ const BlockDataID pdfFieldId_;
+
+ const Vector3<real_t> velocity_;
+
+}; // class MyBoundaryHandling
+
+BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block ) const
+{
+ FlagField_T * flagField = block->getData< FlagField_T >( flagFieldId_ );
+ PdfField_T * pdfField = block->getData< PdfField_T > ( pdfFieldId_ );
+
+ const flag_t fluid = flagField->registerFlag( Fluid_Flag ); // register the fluid flag at the flag field
+
+ return new BoundaryHandling_T( "boundary handling", flagField, fluid,
+ boost::tuples::make_tuple( BottomWall_T( "bottom", BottomWall_Flag, pdfField ),
+ TopWall_T( "top", TopWall_Flag, pdfField, velocity_ ) ) );
+}
+
+
+
+//////////
+// MAIN //
+//////////
+
+int main( int argc, char ** argv )
+{
+ mpi::Environment env( argc, argv );
+
+ int processes( MPIManager::instance()->numProcesses() );
+ if( processes != 1 && processes != 2 && processes != 4 && processes != 8 )
+ WALBERLA_ABORT( "The number of processes must be equal to either 1, 2, 4, or 8!" );
+
+ const uint_t xBlocks = ( processes == 8 ) ? uint_t(2) : uint_t(1);
+ const uint_t yBlocks = ( processes >= 4 ) ? uint_t(2) : uint_t(1);
+ const uint_t zBlocks = ( processes >= 2 ) ? uint_t(2) : uint_t(1);
+
+ const uint_t L = uint_t(16);
+
+ auto blocks = blockforest::createUniformBlockGrid( xBlocks, yBlocks, zBlocks,
+ L / xBlocks,
+ L / yBlocks,
+ L / zBlocks,
+ real_t(1),
+ true,
+ true, true, false);
+
+ const real_t omega = real_c(1);
+ LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ) );
+
+ const Vector3<real_t> velocity(real_t(0.01), real_t(0), real_t(0));
+ BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel );
+ BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" );
+ BlockDataID boundaryHandlingId = blocks->addBlockData< BoundaryHandling_T >( MyBoundaryHandling( flagFieldId, pdfFieldId, velocity ),
+ "boundary handling" );
+
+ // set boundary conditions
+ geometry::initializer::BoundaryFromBody< BoundaryHandling_T > bodyInitializer( *blocks, boundaryHandlingId );
+ geometry::initializer::BoundaryFromDomainBorder< BoundaryHandling_T > borderInitializer( *blocks, boundaryHandlingId );
+ borderInitializer.init( TopWall_Flag, stencil::T, cell_idx_t(-1) );
+ borderInitializer.init( BottomWall_Flag, stencil::B, cell_idx_t(-1) );
+ geometry::setNonBoundaryCellsToDomain< BoundaryHandling_T >( *blocks, boundaryHandlingId );
+
+ uint_t timeSteps = uint_c(2000);
+ if( argc > 1 )
+ timeSteps = boost::lexical_cast<uint_t>( argv[1] );
+ SweepTimeloop timeloop( blocks->getBlockStorage(), timeSteps );
+
+ blockforest::communication::UniformBufferedScheme< CommunicationStencil_T > communication( blocks );
+ communication.addPackInfo( make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >( pdfFieldId ) );
+
+ timeloop.add() << BeforeFunction( communication, "communication" )
+ << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingId ), "boundary handling" );
+
+ timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldId, flagFieldId, Fluid_Flag ) ), "LB stream & collide" );
+
+ timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
+
+ timeloop.run();
+
+ // get the force acting on the walls
+ real_t fXBottom = real_t(0), fXTop = real_t(0);
+ for( auto block = blocks->begin(); block != blocks->end(); ++block )
+ {
+ FlagField_T * flagField = block->getData< FlagField_T >( flagFieldId );
+ BoundaryHandling_T * boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingId );
+ auto topWallUid = boundaryHandling->getBoundaryUID( TopWall_Flag );
+ auto bottomWallUid = boundaryHandling->getBoundaryUID( BottomWall_Flag );
+ auto & topBc = boundaryHandling->getBoundaryCondition<TopWall_T>( topWallUid );
+ auto & bottomBc = boundaryHandling->getBoundaryCondition<BottomWall_T>( bottomWallUid );
+ auto topFlag = flagField->getFlag(TopWall_Flag);
+ auto bottomFlag = flagField->getFlag(BottomWall_Flag);
+
+ // need to iterate over ghost cells too because they may contain part of the force
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ_OMP(flagField, omp parallel for schedule(static) reduction(+:fXTop) reduction(+:fXBottom),
+ {
+ if( flagField->isFlagSet(x, y, z, topFlag))
+ {
+ const Vector3<real_t> & cellForce = topBc.getForce(x, y, z);
+ fXTop += cellForce[0];
+ }
+
+ if( flagField->isFlagSet(x, y, z, bottomFlag))
+ {
+ const Vector3<real_t> & cellForce = bottomBc.getForce(x, y, z);
+ fXBottom += cellForce[0];
+ }
+ });
+ }
+
+ mpi::allReduceInplace( fXTop, mpi::SUM );
+ mpi::allReduceInplace( fXBottom, mpi::SUM );
+
+ real_t viscosity = lbm::collision_model::viscosityFromOmega( omega );
+ real_t shearRate = velocity[0] / real_c(L);
+ real_t wallArea = real_c(L) * real_c(L);
+ real_t analyticalForce = shearRate * viscosity * wallArea;
+
+ WALBERLA_LOG_RESULT_ON_ROOT("Expected force: " << analyticalForce);
+ WALBERLA_LOG_RESULT_ON_ROOT("Actual force top: " << fXTop);
+ WALBERLA_LOG_RESULT_ON_ROOT("Actual force bottom: " << fXBottom);
+
+ real_t errorTop = std::fabs( ( std::fabs(fXTop) - analyticalForce ) / analyticalForce );
+ real_t errorBottom = std::fabs( ( std::fabs(fXBottom) - analyticalForce ) / analyticalForce );
+
+ WALBERLA_LOG_RESULT_ON_ROOT("Relative error top: " << (real_t(100)*errorTop) << " %");
+ WALBERLA_LOG_RESULT_ON_ROOT("Relative error bottom: " << (real_t(100)*errorBottom) << " %");
+
+
+ if( timeSteps > 100 )
+ {
+ WALBERLA_CHECK_LESS(errorTop, 0.01, "Error in force at top wall is too high! value = " << fXTop << " vs analytical = " << analyticalForce);
+ WALBERLA_CHECK_LESS(errorBottom, 0.01, "Error in force at bottom wall is too high! value = " << fXBottom << " vs analytical = " << analyticalForce);
+ WALBERLA_LOG_PROGRESS_ON_ROOT("Test succeeded");
+ }
+
+ return EXIT_SUCCESS;
+}
--
GitLab