Skip to content
Snippets Groups Projects
Commit 0a5b7370 authored by Michael Kuron's avatar Michael Kuron :mortar_board: Committed by Christoph Rettinger
Browse files

Added functionality to LBM boundaries to optionally obtain the hydrodynamic... 

Added functionality to LBM boundaries to optionally obtain the hydrodynamic force on those boundaries. Also added tests.
parent 9863e539
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 685 additions and 61 deletions
src/lbm/boundary/DynamicUBB.h
+ 27
− 1
View file @ 0a5b7370
...@@ -51,12 +51,14 @@ namespace lbm { ...@@ -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 // 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'. // 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. // 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> class DynamicUBB : public Boundary<flag_t>
{ {
typedef lbm::PdfField< LatticeModel_T > PDFField; typedef lbm::PdfField< LatticeModel_T > PDFField;
typedef typename LatticeModel_T::Stencil Stencil; typedef typename LatticeModel_T::Stencil Stencil;
typedef GhostLayerField< Vector3<real_t>, 1 > ForceField;
public: public:
static const bool threadsafe = true; static const bool threadsafe = true;
...@@ -83,6 +85,9 @@ public: ...@@ -83,6 +85,9 @@ public:
{ {
velocity_( time_ ); 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, DynamicUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField,
const uint_t level, const VelocityFunctor_T & velocity, const AABB & aabb ) : const uint_t level, const VelocityFunctor_T & velocity, const AABB & aabb ) :
...@@ -100,6 +105,9 @@ public: ...@@ -100,6 +105,9 @@ public:
time_ = timeTracker_->getTime( level_ ); time_ = timeTracker_->getTime( level_ );
velocity_( time_ ); velocity_( time_ );
} }
if (StoreForce)
force_->setWithGhostLayer( Vector3<real_t>() );
} }
void afterBoundaryTreatment() const {} void afterBoundaryTreatment() const {}
...@@ -131,6 +139,8 @@ public: ...@@ -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 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) // 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], const Vector3< real_t > pos( origin_[0] + real_c(nx) * dx_[0],
origin_[1] + real_c(ny) * dx_[1], origin_[1] + real_c(ny) * dx_[1],
origin_[2] + real_c(nz) * dx_[2] ); origin_[2] + real_c(nz) * dx_[2] );
...@@ -160,6 +170,21 @@ public: ...@@ -160,6 +170,21 @@ public:
real_c(stencil::cy[ dir ]) * vel[1] + real_c(stencil::cy[ dir ]) * vel[1] +
real_c(stencil::cz[ dir ]) * vel[2] ) ); 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: private:
...@@ -177,6 +202,7 @@ private: ...@@ -177,6 +202,7 @@ private:
uint_t level_; uint_t level_;
VelocityFunctor_T velocity_; VelocityFunctor_T velocity_;
shared_ptr<ForceField> force_;
}; // class DynamicUBB }; // class DynamicUBB
......
src/lbm/boundary/NoSlip.h
+ 34
− 3
View file @ 0a5b7370
...@@ -42,7 +42,7 @@ namespace lbm { ...@@ -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> class NoSlip : public Boundary<flag_t>
{ {
protected: protected:
...@@ -50,6 +50,8 @@ protected: ...@@ -50,6 +50,8 @@ protected:
typedef PdfField< LatticeModel_T > PDFField; typedef PdfField< LatticeModel_T > PDFField;
typedef typename LatticeModel_T::Stencil Stencil; typedef typename LatticeModel_T::Stencil Stencil;
typedef GhostLayerField< Vector3<real_t>, 1 > ForceField;
public: public:
static const bool threadsafe = true; static const bool threadsafe = true;
...@@ -60,11 +62,21 @@ public: ...@@ -60,11 +62,21 @@ public:
NoSlip( const BoundaryUID& boundaryUID, const FlagUID& uid, PDFField* const pdfField ) : 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 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 {} void afterBoundaryTreatment() const {}
template< typename Buffer_T > template< typename Buffer_T >
...@@ -95,7 +107,24 @@ public: ...@@ -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 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) // 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] ); 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: private:
...@@ -104,6 +133,8 @@ private: ...@@ -104,6 +133,8 @@ private:
PDFField* const pdfField_; PDFField* const pdfField_;
shared_ptr<ForceField> force_;
}; // class NoSlip }; // class NoSlip
......
src/lbm/boundary/ParserUBB.h
+ 60
− 34
View file @ 0a5b7370
...@@ -47,12 +47,14 @@ namespace lbm { ...@@ -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> class ParserUBB : public Boundary<flag_t>
{ {
typedef lbm::PdfField< LatticeModel_T > PDFField; typedef lbm::PdfField< LatticeModel_T > PDFField;
typedef typename LatticeModel_T::Stencil Stencil; typedef typename LatticeModel_T::Stencil Stencil;
typedef GhostLayerField< Vector3<real_t>, 1 > ForceField;
public: public:
static const bool threadsafe = true; static const bool threadsafe = true;
...@@ -118,6 +120,9 @@ public: ...@@ -118,6 +120,9 @@ public:
{ {
if( timeTracker_ ) if( timeTracker_ )
time_ = timeTracker_->getTime( level_ ); time_ = timeTracker_->getTime( level_ );
if (StoreForce)
force_->setWithGhostLayer( Vector3<real_t>() );
} }
void afterBoundaryTreatment() const {} void afterBoundaryTreatment() const {}
...@@ -145,6 +150,12 @@ public: ...@@ -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 ) const;
inline Vector3<real_t> getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z, real_t t ) 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: private:
const FlagUID uid_; const FlagUID uid_;
...@@ -161,11 +172,12 @@ private: ...@@ -161,11 +172,12 @@ private:
shared_ptr<ParserField> parserField_; shared_ptr<ParserField> parserField_;
shared_ptr<VelocityField> velocityField_; shared_ptr<VelocityField> velocityField_;
shared_ptr<ForceField> force_;
}; // class ParserUBB }; // class ParserUBB
template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce> template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::Parser( const Config::BlockHandle & config ) inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce>::Parser::Parser( const Config::BlockHandle & config )
: parsers_(), equations_(), timeDependent_( false ) : parsers_(), equations_(), timeDependent_( false )
{ {
if( !config ) if( !config )
...@@ -194,8 +206,8 @@ inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser:: ...@@ -194,8 +206,8 @@ inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::
} }
} }
template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce> template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::Parser( std::array< std::string, 3 > & equations ) inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce>::Parser::Parser( std::array< std::string, 3 > & equations )
: parsers_(), equations_( equations ), timeDependent_( false ) : parsers_(), equations_( equations ), timeDependent_( false )
{ {
if( equations_[0].length() > 0 ) if( equations_[0].length() > 0 )
...@@ -218,16 +230,16 @@ inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser:: ...@@ -218,16 +230,16 @@ inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::
} }
} }
template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce > template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
inline ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::Velocity::Velocity( const Config::BlockHandle & config ) 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_[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_[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); velocity_[2] = ( config && config.isDefined( "z" ) ) ? config.getParameter<real_t>( "z" ) : real_c(0.0);
} }
template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce> template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::operator()( const Vector3< real_t > & x, const real_t t ) const 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; std::map< std::string, double > symbols;
symbols["x"] = x[0]; symbols["x"] = x[0];
...@@ -242,8 +254,8 @@ Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce ...@@ -242,8 +254,8 @@ Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce
return v; return v;
} }
template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce> template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::operator()( const Vector3< real_t > & x ) const Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce, StoreForce>::Parser::operator()( const Vector3< real_t > & x ) const
{ {
WALBERLA_ASSERT( !timeDependent_ ); WALBERLA_ASSERT( !timeDependent_ );
...@@ -261,8 +273,8 @@ Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce ...@@ -261,8 +273,8 @@ Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce
template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce> template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce>
inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::ParserUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField, 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, FlagField<flag_t> * const flagField, const shared_ptr< TimeTracker > & timeTracker,
const uint_t level, const AABB & aabb ) const uint_t level, const AABB & aabb )
: Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ), timeTracker_( timeTracker ), time_( real_t(0) ), level_( level ) : 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 ...@@ -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 ); 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 ); 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 > 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 ) ) if( parserField_->get( x, y, z ) )
{ {
...@@ -306,9 +321,9 @@ inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::p ...@@ -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 > 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; bool isparser;
buffer >> isparser; buffer >> isparser;
...@@ -329,16 +344,16 @@ inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::r ...@@ -329,16 +344,16 @@ 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>
inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::ParserUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField, 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 ) FlagField<flag_t> * const flagField, const uint_t level, const AABB & aabb )
: ParserUBB( boundaryUID, uid, pdfField, flagField, nullptr, level, aabb ) : ParserUBB( boundaryUID, uid, pdfField, flagField, nullptr, level, aabb )
{} {}
template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce > template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
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, 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 ) const BoundaryConfiguration & parser )
{ {
WALBERLA_ASSERT_EQUAL( dynamic_cast< const Parser * >( &parser ), &parser ); WALBERLA_ASSERT_EQUAL( dynamic_cast< const Parser * >( &parser ), &parser );
...@@ -369,8 +384,8 @@ inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::r ...@@ -369,8 +384,8 @@ 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 >
inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & parser ) 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_EQUAL( dynamic_cast< const Parser * >( &parser ), &parser );
WALBERLA_ASSERT_NOT_NULLPTR( parserField_ ); WALBERLA_ASSERT_NOT_NULLPTR( parserField_ );
...@@ -408,9 +423,9 @@ inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::r ...@@ -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 > 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 ) const BoundaryConfiguration & parser )
{ {
WALBERLA_ASSERT_EQUAL( dynamic_cast< const Parser * >( &parser ), &parser ); WALBERLA_ASSERT_EQUAL( dynamic_cast< const Parser * >( &parser ), &parser );
...@@ -451,20 +466,20 @@ inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::r ...@@ -451,20 +466,20 @@ 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 >
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 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; 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 #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 ) const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t mask )
#else #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*/ ) const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t /*mask*/ )
#endif #endif
{ {
...@@ -475,6 +490,8 @@ inline void ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Par ...@@ -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 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) // current implementation of this boundary condition (ParserUBB)
const real_t pdf_old = pdfField_->get( x, y, z, Stencil::idx[dir] );
Vector3<real_t> velocity; Vector3<real_t> velocity;
if( parserField_->get(nx,ny,nz) ) if( parserField_->get(nx,ny,nz) )
{ {
...@@ -512,18 +529,27 @@ inline void ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Par ...@@ -512,18 +529,27 @@ inline void ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Par
real_c(stencil::cy[ dir ]) * vel[1] + real_c(stencil::cy[ dir ]) * vel[1] +
real_c(stencil::cz[ dir ]) * vel[2] ) ); 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 > template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
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 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_ ); return getValue( x, y, z, time_ );
} }
template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce > template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
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 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; Vector3<real_t> velocity;
if( parserField_->get(x,y,z) ) if( parserField_->get(x,y,z) )
......
src/lbm/boundary/SimpleUBB.h
+ 41
− 4
View file @ 0a5b7370
...@@ -45,12 +45,14 @@ namespace lbm { ...@@ -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> class SimpleUBB : public Boundary<flag_t>
{ {
typedef PdfField< LatticeModel_T > PDFField; typedef PdfField< LatticeModel_T > PDFField;
typedef typename LatticeModel_T::Stencil Stencil; typedef typename LatticeModel_T::Stencil Stencil;
typedef GhostLayerField< Vector3<real_t>, 1 > ForceField;
public: public:
static const bool threadsafe = true; static const bool threadsafe = true;
...@@ -61,14 +63,30 @@ public: ...@@ -61,14 +63,30 @@ public:
SimpleUBB( const BoundaryUID& boundaryUID, const FlagUID& uid, PDFField* const pdfField, const Vector3< real_t > & velocity ) : 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 ) : 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 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 {} void afterBoundaryTreatment() const {}
template< typename Buffer_T > template< typename Buffer_T >
...@@ -99,6 +117,8 @@ public: ...@@ -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 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) // current implementation of this boundary condition (SimpleUBB)
const real_t pdf_old = pdfField_->get( x, y, z, Stencil::idx[dir] );
if( LatticeModel_T::compressible ) if( LatticeModel_T::compressible )
{ {
const auto density = pdfField_->getDensity(x,y,z); const auto density = pdfField_->getDensity(x,y,z);
...@@ -122,6 +142,21 @@ public: ...@@ -122,6 +142,21 @@ public:
real_c(stencil::cy[ dir ]) * velocity[1] + real_c(stencil::cy[ dir ]) * velocity[1] +
real_c(stencil::cz[ dir ]) * velocity[2] ) ); 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: private:
...@@ -132,6 +167,8 @@ private: ...@@ -132,6 +167,8 @@ private:
const Vector3< real_t > velocity_; const Vector3< real_t > velocity_;
shared_ptr<ForceField> force_;
}; // class SimpleUBB }; // class SimpleUBB
......
src/lbm/boundary/UBB.h
+ 51
− 19
View file @ 0a5b7370
...@@ -46,7 +46,7 @@ namespace lbm { ...@@ -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> class UBB : public Boundary<flag_t>
{ {
typedef PdfField< LatticeModel_T > PDFField; typedef PdfField< LatticeModel_T > PDFField;
...@@ -54,6 +54,8 @@ class UBB : public Boundary<flag_t> ...@@ -54,6 +54,8 @@ class UBB : public Boundary<flag_t>
typedef GhostLayerField< Vector3<real_t>, 1 > VelField; typedef GhostLayerField< Vector3<real_t>, 1 > VelField;
typedef GhostLayerField< Vector3<real_t>, 1 > ForceField;
public: public:
static const bool threadsafe = true; static const bool threadsafe = true;
...@@ -89,7 +91,7 @@ public: ...@@ -89,7 +91,7 @@ public:
void pushFlags( std::vector< FlagUID > & uids ) const { uids.push_back( uid_ ); } void pushFlags( std::vector< FlagUID > & uids ) const { uids.push_back( uid_ ); }
void beforeBoundaryTreatment() const {} inline void beforeBoundaryTreatment() const;
void afterBoundaryTreatment() const {} void afterBoundaryTreatment() const {}
template< typename Buffer_T > template< typename Buffer_T >
...@@ -110,19 +112,26 @@ public: ...@@ -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); } 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: private:
const FlagUID uid_; const FlagUID uid_;
PDFField* const pdfField_; PDFField* const pdfField_;
shared_ptr<VelField> vel_; shared_ptr<VelField> vel_;
shared_ptr<ForceField> force_;
}; // class UBB }; // class UBB
template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce > template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
inline UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::Velocity::Velocity( const Config::BlockHandle & config ) 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_[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_[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 ...@@ -131,8 +140,8 @@ inline UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::Velocity::Ve
template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce > template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce, bool StoreForce >
inline UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::UBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField* const pdfField, FlagField<flag_t> * const flagField ) : 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 ) Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField )
{ {
...@@ -141,30 +150,42 @@ inline UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::UBB( const B ...@@ -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 ); vel_ = make_shared<VelField>( pdfField_->xSize(), pdfField_->ySize(), pdfField_->zSize(), flagField->nrOfGhostLayers(), field::zyxf );
else else
vel_ = make_shared<VelField>( pdfField_->xSize(), pdfField_->ySize(), pdfField_->zSize(), pdfField_->nrOfGhostLayers(), field::zyxf ); 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 > 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 ); 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 > 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 ); 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 >
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, 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 ) const BoundaryConfiguration & velocity )
{ {
WALBERLA_ASSERT_EQUAL( dynamic_cast< const Velocity * >( &velocity ), &velocity ); WALBERLA_ASSERT_EQUAL( dynamic_cast< const Velocity * >( &velocity ), &velocity );
...@@ -177,8 +198,8 @@ inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registe ...@@ -177,8 +198,8 @@ 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 >
inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & velocity ) 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_EQUAL( dynamic_cast< const Velocity * >( &velocity ), &velocity );
WALBERLA_ASSERT_NOT_NULLPTR( vel_ ); WALBERLA_ASSERT_NOT_NULLPTR( vel_ );
...@@ -191,9 +212,9 @@ inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registe ...@@ -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 > 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 ) const BoundaryConfiguration & velocity )
{ {
WALBERLA_ASSERT_EQUAL( dynamic_cast< const Velocity * >( &velocity ), &velocity ); WALBERLA_ASSERT_EQUAL( dynamic_cast< const Velocity * >( &velocity ), &velocity );
...@@ -207,12 +228,12 @@ inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registe ...@@ -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 #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 ) const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t mask )
#else #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*/ ) const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz, const flag_t /*mask*/ )
#endif #endif
{ {
...@@ -223,6 +244,8 @@ inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::treatDi ...@@ -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 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) // current implementation of this boundary condition (UBB)
const real_t pdf_old = pdfField_->get( x, y, z, Stencil::idx[dir] );
if( LatticeModel_T::compressible ) if( LatticeModel_T::compressible )
{ {
const auto density = pdfField_->getDensity(x,y,z); const auto density = pdfField_->getDensity(x,y,z);
...@@ -246,6 +269,15 @@ inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::treatDi ...@@ -246,6 +269,15 @@ inline void UBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::treatDi
real_c(stencil::cy[ dir ]) * velocity[1] + real_c(stencil::cy[ dir ]) * velocity[1] +
real_c(stencil::cz[ dir ]) * velocity[2] ) ); 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;
}
} }
......
tests/lbm/CMakeLists.txt
+ 7
− 0
View file @ 0a5b7370
...@@ -23,6 +23,13 @@ waLBerla_execute_test( NAME SimpleDiffusionDirichletDef COMMAND $<TARGET_FILE:Si ...@@ -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 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_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_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 ) waLBerla_execute_test( NAME DiffusionDirichletTest1 COMMAND $<TARGET_FILE:DiffusionDirichlet> -l 16 -w 1 -o 1.2 -v 0.1 -e 0.097711 -t 50 )
......
tests/lbm/boundary/BoundaryForce.cpp 0 → 100644
+ 228
− 0
View file @ 0a5b7370
//======================================================================================================================
//
// 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;
}
tests/lbm/boundary/BoundaryForceCouette.cpp 0 → 100644
+ 237
− 0
View file @ 0a5b7370
//======================================================================================================================
//
// 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;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment