From 50412516091e95b5c6fdb89fde0ce96532c80452 Mon Sep 17 00:00:00 2001
From: Christoph Rettinger <christoph.rettinger@fau.de>
Date: Thu, 16 Nov 2017 15:32:45 +0100
Subject: [PATCH] [API]: changed pe coupling mapping functions interfaces,
added new mapping functions, adapted test cases and benchmarks
---
.../MotionSingleHeavySphere.cpp | 10 +-
src/pe_coupling/mapping/BodyBBMapping.cpp | 10 +-
src/pe_coupling/mapping/BodyBBMapping.h | 2 +-
src/pe_coupling/mapping/BodyMapping.h | 113 ++++--
.../momentum_exchange_method/BodyMapping.h | 32 +-
.../restoration/PDFReconstruction.h | 4 +-
.../restoration/Reconstructor.h | 2 +-
.../BodyAndVolumeFractionMapping.h | 16 +-
.../PSMUtility.h | 10 +-
.../utility/LubricationCorrection.h | 4 +-
.../BodyAtBlockBoarderCheck.cpp | 2 +-
.../DragForceSphereMEM.cpp | 14 +-
.../DragForceSphereMEMRefinement.cpp | 12 +-
.../LubricationCorrectionMEM.cpp | 18 +-
.../PeriodicParticleChannelMEM.cpp | 5 +-
.../SegreSilberbergMEM.cpp | 14 +-
.../momentum_exchange_method/SquirmerTest.cpp | 379 +++++++++---------
.../TaylorCouetteFlowMEM.cpp | 2 +-
.../TorqueSphereMEM.cpp | 14 +-
.../DragForceSpherePSM.cpp | 43 +-
.../DragForceSpherePSMRefinement.cpp | 43 +-
.../SegreSilberbergPSM.cpp | 91 +----
.../TorqueSpherePSM.cpp | 42 +-
23 files changed, 411 insertions(+), 471 deletions(-)
diff --git a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
index b037a5cbf..0774932f3 100644
--- a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
+++ b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
@@ -900,9 +900,9 @@ int main( int argc, char **argv )
//initialization of the PDFs inside the particles, important for PSM M3
if( psmVariant == PSMVariant::SC1W1 || psmVariant == PSMVariant::SC2W1 || psmVariant == PSMVariant::SC3W1 )
- pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
+ pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
else
- pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
+ pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
}
else
@@ -912,11 +912,11 @@ int main( int argc, char **argv )
if( memVariant == MEMVariant::CLI )
{
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_CLI_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_CLI_Flag );
}else if ( memVariant == MEMVariant::MR ){
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_MR_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_MR_Flag );
}else{
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_BB_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_BB_Flag );
}
}
diff --git a/src/pe_coupling/mapping/BodyBBMapping.cpp b/src/pe_coupling/mapping/BodyBBMapping.cpp
index 7cec1f1fd..393f87781 100644
--- a/src/pe_coupling/mapping/BodyBBMapping.cpp
+++ b/src/pe_coupling/mapping/BodyBBMapping.cpp
@@ -28,7 +28,7 @@
namespace walberla {
namespace pe_coupling {
-CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, const shared_ptr<StructuredBlockStorage> & blockStorage,
+CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, StructuredBlockStorage & blockStorage,
const uint_t numberOfGhostLayersToInclude )
{
@@ -38,16 +38,16 @@ CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, const
if( body->isFinite() )
{
- blockStorage->getCellBBFromAABB( cellBB, body->getAABB(), blockStorage->getLevel(block) );
+ blockStorage.getCellBBFromAABB( cellBB, body->getAABB(), blockStorage.getLevel(block) );
} else
{
- blockStorage->getCellBBFromAABB( cellBB, body->getAABB().getIntersection( blockStorage->getDomain() ), blockStorage->getLevel(block) );
+ blockStorage.getCellBBFromAABB( cellBB, body->getAABB().getIntersection( blockStorage.getDomain() ), blockStorage.getLevel(block) );
}
cellBB.xMin() -= cell_idx_t(1); cellBB.yMin() -= cell_idx_t(1); cellBB.zMin() -= cell_idx_t(1);
cellBB.xMax() += cell_idx_t(1); cellBB.yMax() += cell_idx_t(1); cellBB.zMax() += cell_idx_t(1);
- CellInterval blockBB = blockStorage->getBlockCellBB( block );
+ CellInterval blockBB = blockStorage.getBlockCellBB( block );
cell_idx_t layers = cell_idx_c( numberOfGhostLayersToInclude );
@@ -56,7 +56,7 @@ CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, const
cellBB.intersect( blockBB );
- blockStorage->transformGlobalToBlockLocalCellInterval( cellBB, block );
+ blockStorage.transformGlobalToBlockLocalCellInterval( cellBB, block );
return cellBB;
}
diff --git a/src/pe_coupling/mapping/BodyBBMapping.h b/src/pe_coupling/mapping/BodyBBMapping.h
index a015ac972..59b58680f 100644
--- a/src/pe_coupling/mapping/BodyBBMapping.h
+++ b/src/pe_coupling/mapping/BodyBBMapping.h
@@ -32,7 +32,7 @@ namespace walberla {
namespace pe_coupling {
-CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, const shared_ptr<StructuredBlockStorage> & blockStorage,
+CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, StructuredBlockStorage & blockStorage,
const uint_t numberOfGhostLayersToInclude );
diff --git a/src/pe_coupling/mapping/BodyMapping.h b/src/pe_coupling/mapping/BodyMapping.h
index 6df7a1341..3c689ebb2 100644
--- a/src/pe_coupling/mapping/BodyMapping.h
+++ b/src/pe_coupling/mapping/BodyMapping.h
@@ -37,24 +37,16 @@
namespace walberla {
namespace pe_coupling {
-
-
+// general mapping functions for a given single body on a given single block
template< typename BoundaryHandling_T >
-void mapBody( const pe::BodyID & body, IBlock & block, const shared_ptr<StructuredBlockStorage> & blockStorage,
- const BlockDataID & boundaryHandlingID, const FlagUID & obstacle,
- const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
+void mapBody( const pe::BodyID & body, IBlock & block, StructuredBlockStorage & blockStorage,
+ BoundaryHandling_T * boundaryHandlingPtr, const FlagUID & obstacle )
{
- WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage->getBlockStorage()) );
-
- WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
+ WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage.getBlockStorage()) );
+ WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandlingPtr );
- if( ( fixedBodiesOnly && !body->isFixed() ) /*|| ( moBodiesOnly && !isMOBody( body ) )*/ )
- return;
-
- BoundaryHandling_T * boundaryHandling = block.getData< BoundaryHandling_T >( boundaryHandlingID );
- auto * flagField = boundaryHandling->getFlagField();
+ auto * flagField = boundaryHandlingPtr->getFlagField();
- WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandling );
WALBERLA_ASSERT_NOT_NULLPTR( flagField );
WALBERLA_ASSERT( flagField->flagExists( obstacle ) );
@@ -62,10 +54,10 @@ void mapBody( const pe::BodyID & body, IBlock & block, const shared_ptr<Structur
CellInterval cellBB = getCellBB( body, block, blockStorage, flagField->nrOfGhostLayers() );
- Vector3<real_t> startCellCenter = blockStorage->getBlockLocalCellCenter( block, cellBB.min() );
- const real_t dx = blockStorage->dx( blockStorage->getLevel(block) );
- const real_t dy = blockStorage->dy( blockStorage->getLevel(block) );
- const real_t dz = blockStorage->dz( blockStorage->getLevel(block) );
+ Vector3<real_t> startCellCenter = blockStorage.getBlockLocalCellCenter( block, cellBB.min() );
+ const real_t dx = blockStorage.dx( blockStorage.getLevel(block) );
+ const real_t dy = blockStorage.dy( blockStorage.getLevel(block) );
+ const real_t dz = blockStorage.dz( blockStorage.getLevel(block) );
real_t cz = startCellCenter[2];
for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
@@ -77,58 +69,119 @@ void mapBody( const pe::BodyID & body, IBlock & block, const shared_ptr<Structur
for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
{
if( body->containsPoint(cx,cy,cz) )
- boundaryHandling->forceBoundary( obstacleFlag, x, y, z );
+ boundaryHandlingPtr->forceBoundary( obstacleFlag, x, y, z );
cx += dx;
}
cy += dy;
}
cz += dz;
}
+}
+template< typename BoundaryHandling_T >
+void mapBody( const pe::BodyID & body, IBlock & block, StructuredBlockStorage & blockStorage,
+ const BlockDataID & boundaryHandlingID, const FlagUID & obstacle )
+{
+ WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage.getBlockStorage()) );
+ BoundaryHandling_T * boundaryHandling = block.getData< BoundaryHandling_T >( boundaryHandlingID );
+
+ mapBody(body, block, blockStorage, boundaryHandling, obstacle );
}
+// mapping function to map all bodies from the body storage - with certain properties - onto all blocks
template< typename BoundaryHandling_T >
-void mapBodies( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, const BlockDataID & bodyStorageID, const FlagUID & obstacle,
+void mapBodies( StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
+ const BlockDataID & bodyStorageID, const FlagUID & obstacle,
const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
{
- for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
+ for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle, fixedBodiesOnly, moBodiesOnly );
+ {
+ WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
+
+ if( fixedBodiesOnly && bodyIt->isFixed() )
+ continue;
+
+ mapBody<BoundaryHandling_T>( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle );
+ }
}
}
+
+// mapping function to map all global bodies - with certain properties - onto all blocks
template< typename BoundaryHandling_T >
-void mapGlobalBodies( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle,
+void mapGlobalBodies( StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
+ pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle,
const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
{
- for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
+ for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
{
- for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt)
+ for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt )
{
- mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle, fixedBodiesOnly, moBodiesOnly );
+ WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
+
+ if( fixedBodiesOnly && bodyIt->isFixed() )
+ continue;
+
+ mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle );
}
}
}
+// mapping function to map a given single global body onto all blocks
template< typename BoundaryHandling_T >
void mapGlobalBody( const id_t globalBodySystemID,
- const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle,
- const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
+ StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
+ pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle )
{
- for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
+ for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
{
auto bodyIt = globalBodyStorage.find( globalBodySystemID );
if( bodyIt != globalBodyStorage.end() )
{
- mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle, fixedBodiesOnly, moBodiesOnly );
+ mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle );
}
}
}
+// mapping function to map all global bodies - with certain properties - onto a given single block
+template< typename BoundaryHandling_T >
+void mapGlobalBodiesOnBlock( IBlock & block,
+ StructuredBlockStorage & blockStorage, BoundaryHandling_T * boundaryHandlingPtr,
+ pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle,
+ const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
+{
+ for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt)
+ {
+ WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
+
+ if( fixedBodiesOnly && bodyIt->isFixed() )
+ continue;
+
+ mapBody< BoundaryHandling_T >( *bodyIt, block, blockStorage, boundaryHandlingPtr, obstacle );
+ }
+}
+
+
+// mapping function to map a given single global body onto a given single block
+template< typename BoundaryHandling_T >
+void mapGlobalBodyOnBlock( const id_t globalBodySystemID, IBlock & block,
+ StructuredBlockStorage & blockStorage, BoundaryHandling_T * boundaryHandlingPtr,
+ pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle )
+{
+ auto bodyIt = globalBodyStorage.find( globalBodySystemID );
+ if( bodyIt != globalBodyStorage.end() )
+ {
+ mapBody< BoundaryHandling_T >( *bodyIt, block, blockStorage, boundaryHandlingPtr, obstacle );
+ }
+
+}
+
+
} // namespace pe_coupling
} // namespace walberla
diff --git a/src/pe_coupling/momentum_exchange_method/BodyMapping.h b/src/pe_coupling/momentum_exchange_method/BodyMapping.h
index 732b31fe6..6149ec635 100644
--- a/src/pe_coupling/momentum_exchange_method/BodyMapping.h
+++ b/src/pe_coupling/momentum_exchange_method/BodyMapping.h
@@ -108,7 +108,7 @@ void BodyMapping< BoundaryHandling_T >::operator()( IBlock * const block ) const
// policy: every body manages only its own flags
- CellInterval cellBB = getCellBB( *bodyIt, *block, blockStorage_, flagField->nrOfGhostLayers() );
+ CellInterval cellBB = getCellBB( *bodyIt, *block, *blockStorage_, flagField->nrOfGhostLayers() );
Vector3<real_t> startCellCenter = blockStorage_->getBlockLocalCellCenter( *block, cellBB.min() );
const real_t dx = blockStorage_->dx( blockStorage_->getLevel(*block) );
@@ -174,12 +174,12 @@ void BodyMapping< BoundaryHandling_T >::operator()( IBlock * const block ) const
////////////////////////////
template< typename BoundaryHandling_T >
-void mapMovingBody( const pe::BodyID body, IBlock & block, const shared_ptr<StructuredBlockStorage> & blockStorage,
+void mapMovingBody( const pe::BodyID body, IBlock & block, StructuredBlockStorage & blockStorage,
const BlockDataID & boundaryHandlingID, const BlockDataID & bodyFieldID, const FlagUID & obstacle )
{
typedef Field< pe::BodyID, 1 > BodyField_T;
- WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage->getBlockStorage()) );
+ WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage.getBlockStorage()) );
if( body->isFixed() /*|| !body->isFinite()*/ )
return;
@@ -198,10 +198,10 @@ void mapMovingBody( const pe::BodyID body, IBlock & block, const shared_ptr<Stru
CellInterval cellBB = getCellBB( body, block, blockStorage, flagField->nrOfGhostLayers() );
- Vector3<real_t> startCellCenter = blockStorage->getBlockLocalCellCenter( block, cellBB.min() );
- const real_t dx = blockStorage->dx( blockStorage->getLevel(block) );
- const real_t dy = blockStorage->dy( blockStorage->getLevel(block) );
- const real_t dz = blockStorage->dz( blockStorage->getLevel(block) );
+ Vector3<real_t> startCellCenter = blockStorage.getBlockLocalCellCenter( block, cellBB.min() );
+ const real_t dx = blockStorage.dx( blockStorage.getLevel(block) );
+ const real_t dy = blockStorage.dy( blockStorage.getLevel(block) );
+ const real_t dz = blockStorage.dz( blockStorage.getLevel(block) );
real_t cz = startCellCenter[2];
for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
@@ -230,10 +230,10 @@ void mapMovingBody( const pe::BodyID body, IBlock & block, const shared_ptr<Stru
template< typename BoundaryHandling_T >
-void mapMovingBodies( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, const BlockDataID & bodyStorageID,
+void mapMovingBodies( StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID, const BlockDataID & bodyStorageID,
const BlockDataID & bodyFieldID, const FlagUID & obstacle )
{
- for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
+ for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
{
for (auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt)
mapMovingBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, bodyFieldID, obstacle );
@@ -241,10 +241,10 @@ void mapMovingBodies( const shared_ptr<StructuredBlockStorage> & blockStorage, c
}
template< typename BoundaryHandling_T >
-void mapMovingGlobalBodies( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, pe::BodyStorage & globalBodyStorage,
- const BlockDataID & bodyFieldID, const FlagUID & obstacle )
+void mapMovingGlobalBodies( StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
+ pe::BodyStorage & globalBodyStorage, const BlockDataID & bodyFieldID, const FlagUID & obstacle )
{
- for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
+ for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
{
for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt)
{
@@ -255,13 +255,13 @@ void mapMovingGlobalBodies( const shared_ptr<StructuredBlockStorage> & blockStor
template< typename BoundaryHandling_T >
void mapMovingGlobalBody( const id_t globalBodySystemID,
- const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, pe::BodyStorage & globalBodyStorage,
- const BlockDataID & bodyFieldID, const FlagUID & obstacle )
+ StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
+ pe::BodyStorage & globalBodyStorage, const BlockDataID & bodyFieldID, const FlagUID & obstacle )
{
- for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
+ for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
{
auto bodyIt = globalBodyStorage.find( globalBodySystemID );
- if( bodyIt != globalBodyStorage.end() )
+ if( bodyIt != globalBodyStorage.end() )
{
mapMovingBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, bodyFieldID, obstacle );
}
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h b/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
index 9d751bf33..f7ca0d80a 100644
--- a/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
+++ b/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
@@ -125,7 +125,7 @@ void PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructer_T >
if( bodyIt->isFixed() )
continue;
- CellInterval cellBB = getCellBB( *bodyIt, *block, blockStorage_, numberOfGhostLayersToInclude );
+ CellInterval cellBB = getCellBB( *bodyIt, *block, *blockStorage_, numberOfGhostLayersToInclude );
for( auto cellIt = cellBB.begin(); cellIt != cellBB.end(); ++cellIt )
{
@@ -160,7 +160,7 @@ void PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructer_T >
if( bodyIt->isFixed() )
continue;
- CellInterval cellBB = getCellBB( *bodyIt, *block, blockStorage_, numberOfGhostLayersToInclude );
+ CellInterval cellBB = getCellBB( *bodyIt, *block, *blockStorage_, numberOfGhostLayersToInclude );
for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
{
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h b/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h
index 8acb17a1e..7f9ce1017 100644
--- a/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h
+++ b/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h
@@ -263,7 +263,7 @@ public:
typedef lbm::PdfField< LatticeModel_T > PdfField_T;
typedef Field< pe::BodyID, 1 > BodyField_T;
- ExtrapolationReconstructor( const shared_ptr<StructuredBlockStorage> blockStorage, const BlockDataID & boundaryHandlingID,
+ ExtrapolationReconstructor( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID,
const BlockDataID & pdfFieldID, const BlockDataID & bodyFieldID,
const ExtrapolationDirectionFinder_T & extrapolationDirectionFinder,
const bool & enforceNoSlipConstraintAfterExtrapolation = false )
diff --git a/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
index 88c1092f1..1e461ebbd 100644
--- a/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
+++ b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
@@ -45,7 +45,7 @@ namespace pe_coupling {
* As several bodies could intersect with one cell, the pairs are stored in a vector with the size of the amount of intersecting bodies.
*
*/
-void mapPSMBodyAndVolumeFraction( const pe::BodyID body, IBlock & block, const shared_ptr<StructuredBlockStorage> & blockStorage,
+void mapPSMBodyAndVolumeFraction( const pe::BodyID body, IBlock & block, StructuredBlockStorage & blockStorage,
const BlockDataID bodyAndVolumeFractionFieldID )
{
typedef std::pair< pe::BodyID, real_t > BodyAndVolumeFraction_T;
@@ -63,9 +63,9 @@ void mapPSMBodyAndVolumeFraction( const pe::BodyID body, IBlock & block, const s
// get the cell's center
Vector3<real_t> cellCenter;
- cellCenter = blockStorage->getBlockLocalCellCenter( block, cell );
+ cellCenter = blockStorage.getBlockLocalCellCenter( block, cell );
- const real_t fraction = overlapFractionPe( *body, cellCenter, blockStorage->dx( blockStorage->getLevel( block ) ) );
+ const real_t fraction = overlapFractionPe( *body, cellCenter, blockStorage.dx( blockStorage.getLevel( block ) ) );
// if the cell intersected with the body, store a pointer to that body and the corresponding volume fraction in the field
if( fraction > real_t(0) )
@@ -100,7 +100,7 @@ public:
typedef GhostLayerField< std::vector< BodyAndVolumeFraction_T >, 1 > BodyAndVolumeFractionField_T;
BodyAndVolumeFractionMapping( const shared_ptr<StructuredBlockStorage> & blockStorage,
- const shared_ptr<pe::BodyStorage> globalBodyStorage,
+ const shared_ptr<pe::BodyStorage> & globalBodyStorage,
const BlockDataID & bodyStorageID,
const BlockDataID & bodyAndVolumeFractionFieldID,
const bool mapFixedBodies = false,
@@ -180,7 +180,7 @@ void BodyAndVolumeFractionMapping::initialize()
// only PSM and finite bodies (no planes, etc.) are mapped
if( /*!isPSMBody( *bodyIt ) ||*/ ( bodyIt->isFixed() && !mapFixed_ ) )
continue;
- mapPSMBodyAndVolumeFraction( *bodyIt, *blockIt, blockStorage_, bodyAndVolumeFractionFieldID_ );
+ mapPSMBodyAndVolumeFraction( *bodyIt, *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_ );
lastUpdatedPositionMap_.insert( std::pair< walberla::id_t, Vector3< real_t > >( bodyIt->getSystemID(), bodyIt->getPosition() ) );
}
@@ -188,7 +188,7 @@ void BodyAndVolumeFractionMapping::initialize()
{
for( auto globalBodyIt = globalBodyStorage_->begin(); globalBodyIt != globalBodyStorage_->end(); ++globalBodyIt )
{
- mapPSMBodyAndVolumeFraction( *globalBodyIt, *blockIt, blockStorage_, bodyAndVolumeFractionFieldID_ );
+ mapPSMBodyAndVolumeFraction( *globalBodyIt, *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_ );
}
}
}
@@ -252,7 +252,7 @@ void BodyAndVolumeFractionMapping::updatePSMBodyAndVolumeFraction( const pe::Bod
}
// get bounding box of body
- CellInterval cellBB = getCellBB( body, block, blockStorage_, oldBodyAndVolumeFractionField->nrOfGhostLayers() );
+ CellInterval cellBB = getCellBB( body, block, *blockStorage_, oldBodyAndVolumeFractionField->nrOfGhostLayers() );
// if body has not moved (specified by some epsilon), just reuse old fraction values
if( body->getLinearVel().sqrLength() < velocityUpdatingEpsilonSquared_ &&
@@ -315,7 +315,7 @@ void BodyAndVolumeFractionMapping::updateGlobalPSMBodyAndVolumeFraction( const p
WALBERLA_ASSERT_NOT_NULLPTR( oldBodyAndVolumeFractionField );
// get bounding box of body
- CellInterval cellBB = getCellBB( body, block, blockStorage_, oldBodyAndVolumeFractionField->nrOfGhostLayers() );
+ CellInterval cellBB = getCellBB( body, block, *blockStorage_, oldBodyAndVolumeFractionField->nrOfGhostLayers() );
// copy values of global body to new field
for( auto cellIt = cellBB.begin(); cellIt != cellBB.end(); ++cellIt )
diff --git a/src/pe_coupling/partially_saturated_cells_method/PSMUtility.h b/src/pe_coupling/partially_saturated_cells_method/PSMUtility.h
index 526c0f994..e39c3b7b2 100644
--- a/src/pe_coupling/partially_saturated_cells_method/PSMUtility.h
+++ b/src/pe_coupling/partially_saturated_cells_method/PSMUtility.h
@@ -46,7 +46,7 @@ template < typename LatticeModel_T, int Weighting_T >
Vector3<real_t> getPSMMacroscopicVelocity( const IBlock & block,
lbm::PdfField< LatticeModel_T > * pdfField,
GhostLayerField< std::vector< std::pair< pe::BodyID, real_t > >, 1 > * bodyAndVolumeFractionField,
- const shared_ptr<StructuredBlockStorage> & blockStorage,
+ StructuredBlockStorage & blockStorage,
const Cell & cell )
{
static_assert( LatticeModel_T::compressible == false, "Only works with incompressible models!" );
@@ -60,7 +60,7 @@ Vector3<real_t> getPSMMacroscopicVelocity( const IBlock & block,
for( auto bodyFracIt = bodyAndVolumeFractionField->get( cell ).begin(); bodyFracIt != bodyAndVolumeFractionField->get( cell ).end(); ++bodyFracIt )
{
- const Vector3< real_t > coordsCellCenter = blockStorage->getBlockLocalCellCenter( block, cell );
+ const Vector3< real_t > coordsCellCenter = blockStorage.getBlockLocalCellCenter( block, cell );
const real_t eps = (*bodyFracIt).second;
@@ -89,7 +89,7 @@ Vector3<real_t> getPSMMacroscopicVelocity( const IBlock & block,
* Only the velocity of cells intersecting with bodies is set, pure fluid cells remain unchanged.
*/
template < typename LatticeModel_T, int Weighting_T >
-void initializeDomainForPSM( const shared_ptr<StructuredBlockStorage> & blockStorage,
+void initializeDomainForPSM( StructuredBlockStorage & blockStorage,
const BlockDataID & pdfFieldID, const BlockDataID & bodyAndVolumeFractionFieldID )
{
typedef lbm::PdfField< LatticeModel_T > PdfField_T;
@@ -97,7 +97,7 @@ void initializeDomainForPSM( const shared_ptr<StructuredBlockStorage> & blockSto
typedef GhostLayerField< std::vector< BodyAndVolumeFraction_T >, 1 > BodyAndVolumeFractionField_T;
// iterate all blocks with an iterator 'block'
- for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
+ for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
{
// get the field data out of the block
PdfField_T* pdfField = blockIt->getData< PdfField_T > ( pdfFieldID );
@@ -115,7 +115,7 @@ void initializeDomainForPSM( const shared_ptr<StructuredBlockStorage> & blockSto
for( auto bodyFracIt = bodyAndVolumeFractionField->get(cell).begin(); bodyFracIt != bodyAndVolumeFractionField->get(cell).end(); ++bodyFracIt )
{
- const Vector3< real_t > coordsCellCenter = blockStorage->getBlockLocalCellCenter( *blockIt, cell );
+ const Vector3< real_t > coordsCellCenter = blockStorage.getBlockLocalCellCenter( *blockIt, cell );
const real_t eps = (*bodyFracIt).second;
diff --git a/src/pe_coupling/utility/LubricationCorrection.h b/src/pe_coupling/utility/LubricationCorrection.h
index 0903c3ffd..e64c74030 100644
--- a/src/pe_coupling/utility/LubricationCorrection.h
+++ b/src/pe_coupling/utility/LubricationCorrection.h
@@ -42,8 +42,8 @@ class LubricationCorrection
public:
// constructor
- LubricationCorrection ( const shared_ptr<StructuredBlockStorage> & blockStorage, shared_ptr<pe::BodyStorage> globalBodyStorage, const BlockDataID & bodyStorageID,
- real_t dynamicViscosity, real_t cutOffDistance = real_t(2) / real_t(3) )
+ LubricationCorrection ( const shared_ptr<StructuredBlockStorage> & blockStorage, const shared_ptr<pe::BodyStorage> & globalBodyStorage,
+ const BlockDataID & bodyStorageID, real_t dynamicViscosity, real_t cutOffDistance = real_t(2) / real_t(3) )
: blockStorage_ ( blockStorage )
, globalBodyStorage_( globalBodyStorage )
, bodyStorageID_( bodyStorageID )
diff --git a/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp b/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
index 03cd4c395..16e0afc65 100644
--- a/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
@@ -354,7 +354,7 @@ int main( int argc, char **argv )
MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
// initially map pe bodies into the LBM simulation
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
///////////////
// TIME LOOP //
diff --git a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
index 095fed5f1..581f51deb 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
@@ -13,7 +13,7 @@
// 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 DragForceSphereMEMPe.cpp
+//! \file DragForceSphereMEM.cpp
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
@@ -64,7 +64,7 @@
#include <iomanip>
#include <iostream>
-namespace drag_force_sphere_mem_pe
+namespace drag_force_sphere_mem
{
///////////
@@ -527,13 +527,13 @@ int main( int argc, char **argv )
if( method == MEMVariant::CLI )
{
// uses a higher order boundary condition (CLI)
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
}else if ( method == MEMVariant::MR ){
// uses a higher order boundary condition (MR)
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_MR_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_MR_Flag );
}else{
// uses standard bounce back boundary conditions
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag );
}
// since external forcing is applied, the evaluation of the velocity has to be carried out directly after the streaming step
@@ -608,8 +608,8 @@ int main( int argc, char **argv )
}
-} //namespace drag_force_sphere_mem_pe
+} //namespace drag_force_sphere_mem
int main( int argc, char **argv ){
- drag_force_sphere_mem_pe::main(argc, argv);
+ drag_force_sphere_mem::main(argc, argv);
}
diff --git a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
index 40fb413f6..49ac5f453 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
@@ -13,7 +13,7 @@
// 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 DragForceSphereMEMPeRefinement.cpp
+//! \file DragForceSphereMEMRefinement.cpp
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
@@ -71,7 +71,7 @@
#include <iomanip>
#include <iostream>
-namespace drag_force_sphere_mem_pe_refinement
+namespace drag_force_sphere_mem_refinement
{
///////////
@@ -549,10 +549,10 @@ int main( int argc, char **argv )
if( MO_CLI )
{
// uses a higher order boundary condition (CLI)
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
}else{
// uses standard bounce back boundary conditions
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag );
}
@@ -647,8 +647,8 @@ int main( int argc, char **argv )
}
-} //namespace drag_force_sphere_mem_pe_refinement
+} //namespace drag_force_sphere_mem_refinement
int main( int argc, char **argv ){
- drag_force_sphere_mem_pe_refinement::main(argc, argv);
+ drag_force_sphere_mem_refinement::main(argc, argv);
}
diff --git a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
index 89068030d..62a3cfd2c 100644
--- a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
@@ -13,7 +13,7 @@
// 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 LubricationCorrectionMEMPe.cpp
+//! \file LubricationCorrectionMEM.cpp
//! \ingroup pe_coupling
//! \author Kristina Pickl <kristina.pickl@fau.de>
//! \author Dominik Bartuschat
@@ -58,7 +58,7 @@
#include "pe_coupling/utility/all.h"
-namespace lubrication_correction_mem_pe
+namespace lubrication_correction_mem
{
///////////
@@ -127,7 +127,7 @@ class EvaluateLubricationForce
public:
EvaluateLubricationForce( const shared_ptr< StructuredBlockStorage > & blocks,
const BlockDataID & bodyStorageID,
- pe::BodyStorage & globalBodyStorage,
+ const shared_ptr<pe::BodyStorage> & globalBodyStorage,
uint_t id1, uint_t id2, pe::Vec3 vel, real_t nu_L, real_t radius,
SweepTimeloop* timeloop, bool print, bool sphSphTest, bool sphWallTest )
: blocks_( blocks ), bodyStorageID_( bodyStorageID ), globalBodyStorage_( globalBodyStorage ),
@@ -286,7 +286,7 @@ private:
pe::SphereID sphereI = ( *curSphereIt );
if ( sphereI->getID() == id1_ )
{
- for( auto globalBodyIt = globalBodyStorage_.begin(); globalBodyIt != globalBodyStorage_.end(); ++globalBodyIt)
+ for( auto globalBodyIt = globalBodyStorage_->begin(); globalBodyIt != globalBodyStorage_->end(); ++globalBodyIt)
{
if( globalBodyIt->getID() == id2_ )
{
@@ -368,7 +368,7 @@ private:
shared_ptr< StructuredBlockStorage > blocks_;
const BlockDataID bodyStorageID_;
- pe::BodyStorage & globalBodyStorage_;
+ shared_ptr<pe::BodyStorage> globalBodyStorage_;
uint_t id1_;
uint_t id2_;
@@ -922,7 +922,7 @@ int main( int argc, char **argv )
}
// map pe bodies into the LBM simulation
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
///////////////
// TIME LOOP //
@@ -970,7 +970,7 @@ int main( int argc, char **argv )
timeloop.addFuncAfterTimeStep( pe_coupling::LubricationCorrection( blocks, globalBodyStorage, bodyStorageID, nu_L ), "Lubrication Force" );
// perform lubrication evaluation
- timeloop.addFuncAfterTimeStep( EvaluateLubricationForce( blocks, bodyStorageID, *globalBodyStorage, id1, id2, velocity,
+ timeloop.addFuncAfterTimeStep( EvaluateLubricationForce( blocks, bodyStorageID, globalBodyStorage, id1, id2, velocity,
nu_L, radius, &timeloop, fileIO, sphSphTest, sphWallTest ), "Evaluate Lubrication Force" );
// reset the forces and apply a constant velocity
@@ -992,10 +992,10 @@ int main( int argc, char **argv )
return 0;
}
-} //namespace lubrication_correction_mem_pe
+} //namespace lubrication_correction_mem
int main( int argc, char **argv ){
- lubrication_correction_mem_pe::main(argc, argv);
+ lubrication_correction_mem::main(argc, argv);
}
diff --git a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
index ba33e553f..2e351fa79 100644
--- a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
@@ -502,10 +502,11 @@ int main( int argc, char **argv )
// special care has to be taken here that only the global spheres (not the planes) are mapped since the planes would overwrite the already set boundary flags
for( auto globalBodyIt = globalBodiesToBeMapped.begin(); globalBodyIt != globalBodiesToBeMapped.end(); ++globalBodyIt )
{
- pe_coupling::mapGlobalBody< BoundaryHandling_T >( *globalBodyIt, blocks, boundaryHandlingID, *globalBodyStorage, NoSlip_Flag, false );
+ pe_coupling::mapGlobalBody< BoundaryHandling_T >( *globalBodyIt, *blocks, boundaryHandlingID, *globalBodyStorage, NoSlip_Flag );
}
+
// moving bodies are handled by the momentum exchange method
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
///////////////
// TIME LOOP //
diff --git a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
index 3bff481c7..055c68b0d 100644
--- a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
@@ -13,7 +13,7 @@
// 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 SegreSilberbergMEMPe.cpp
+//! \file SegreSilberbergMEM.cpp
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
@@ -65,7 +65,7 @@
#include "field/vtk/all.h"
#include "lbm/vtk/all.h"
-namespace segre_silberberg_mem_pe
+namespace segre_silberberg_mem
{
///////////
@@ -609,13 +609,13 @@ int main( int argc, char **argv )
if( MO_CLI )
{
// uses a higher order boundary condition (CLI)
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
}else if ( MO_MR ){
// uses a higher order boundary condition (MR)
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_MR_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_MR_Flag );
}else{
// uses standard bounce back boundary conditions
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag );
}
///////////////
@@ -820,8 +820,8 @@ int main( int argc, char **argv )
return EXIT_SUCCESS;
}
-} // namespace segre_silberberg_mem_pe
+} // namespace segre_silberberg_mem
int main( int argc, char **argv ){
- segre_silberberg_mem_pe::main(argc, argv);
+ segre_silberberg_mem::main(argc, argv);
}
diff --git a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
index f0c248c14..5680cb5ab 100644
--- a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
@@ -66,6 +66,11 @@
#include "lbm/vtk/Velocity.h"
#include "vtk/VTKOutput.h"
+
+namespace squirmer
+{
+
+
///////////
// USING //
///////////
@@ -75,45 +80,44 @@ using walberla::uint_t;
// PDF field, flag field & body field
typedef lbm::force_model::None ForceModel_T;
-typedef lbm::D3Q19< lbm::collision_model::TRT, false, ForceModel_T > LatticeModel_T;
+typedef lbm::D3Q19<lbm::collision_model::TRT, false, ForceModel_T> LatticeModel_T;
-typedef LatticeModel_T::Stencil Stencil_T;
-typedef lbm::PdfField< LatticeModel_T > PdfField_T;
+typedef LatticeModel_T::Stencil Stencil_T;
+typedef lbm::PdfField<LatticeModel_T> PdfField_T;
-typedef walberla::uint8_t flag_t;
-typedef FlagField< flag_t > FlagField_T;
-typedef GhostLayerField< pe::BodyID, 1 > PeBodyField_T;
+typedef walberla::uint8_t flag_t;
+typedef FlagField<flag_t> FlagField_T;
+typedef GhostLayerField<pe::BodyID, 1> PeBodyField_T;
const uint_t FieldGhostLayers = 1;
// boundary handling
-typedef pe_coupling::SimpleBB < LatticeModel_T, FlagField_T > MO_BB_T;
+typedef pe_coupling::SimpleBB<LatticeModel_T, FlagField_T> MO_BB_T;
-typedef boost::tuples::tuple< MO_BB_T > BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef boost::tuples::tuple<MO_BB_T> BoundaryConditions_T;
+typedef BoundaryHandling<FlagField_T, Stencil_T, BoundaryConditions_T> BoundaryHandling_T;
-typedef boost::tuple<pe::Squirmer> BodyTypeTuple ;
+typedef boost::tuple<pe::Squirmer> BodyTypeTuple;
///////////
// FLAGS //
///////////
-const FlagUID Fluid_Flag ( "fluid" );
-const FlagUID MO_BB_Flag ( "moving obstacle BB" );
-const FlagUID FormerMO_BB_Flag ( "former moving obstacle BB" );
+const FlagUID Fluid_Flag("fluid");
+const FlagUID MO_BB_Flag("moving obstacle BB");
+const FlagUID FormerMO_BB_Flag("former moving obstacle BB");
/////////////////////////////////////
// BOUNDARY HANDLING CUSTOMIZATION //
/////////////////////////////////////
-class MyBoundaryHandling
-{
+class MyBoundaryHandling {
public:
- MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID, const BlockDataID & bodyFieldID ) :
- flagFieldID_( flagFieldID ), pdfFieldID_( pdfFieldID ), bodyFieldID_ ( bodyFieldID ) {}
+ MyBoundaryHandling(const BlockDataID &flagFieldID, const BlockDataID &pdfFieldID, const BlockDataID &bodyFieldID) :
+ flagFieldID_(flagFieldID), pdfFieldID_(pdfFieldID), bodyFieldID_(bodyFieldID) {}
- BoundaryHandling_T * operator()( IBlock* const block, const StructuredBlockStorage* const storage ) const;
+ BoundaryHandling_T *operator()(IBlock *const block, const StructuredBlockStorage *const storage) const;
private:
@@ -127,141 +131,139 @@ private:
// VTK OUTPUT HELPER FUNCTION //
////////////////////////////////
-BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const StructuredBlockStorage * const storage ) const
-{
- WALBERLA_ASSERT_NOT_NULLPTR( block );
- WALBERLA_ASSERT_NOT_NULLPTR( storage );
+BoundaryHandling_T *
+MyBoundaryHandling::operator()(IBlock *const block, const StructuredBlockStorage *const storage) const {
+ WALBERLA_ASSERT_NOT_NULLPTR(block);
+ WALBERLA_ASSERT_NOT_NULLPTR(storage);
- FlagField_T * flagField = block->getData< FlagField_T >( flagFieldID_ );
- PdfField_T * pdfField = block->getData< PdfField_T > ( pdfFieldID_ );
- PeBodyField_T * bodyField = block->getData< PeBodyField_T >( bodyFieldID_ );
+ FlagField_T *flagField = block->getData<FlagField_T>(flagFieldID_);
+ PdfField_T *pdfField = block->getData<PdfField_T>(pdfFieldID_);
+ PeBodyField_T *bodyField = block->getData<PeBodyField_T>(bodyFieldID_);
- const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
+ const auto fluid = flagField->flagExists(Fluid_Flag) ? flagField->getFlag(Fluid_Flag) : flagField->registerFlag(
+ Fluid_Flag);
- BoundaryHandling_T * handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
- boost::tuples::make_tuple( MO_BB_T ( "MO_BB", MO_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+ BoundaryHandling_T *handling = new BoundaryHandling_T("fixed obstacle boundary handling", flagField, fluid,
+ boost::tuples::make_tuple(
+ MO_BB_T("MO_BB", MO_BB_Flag, pdfField, flagField,
+ bodyField, fluid, *storage, *block)));
- handling->fillWithDomain( FieldGhostLayers );
+ handling->fillWithDomain(FieldGhostLayers);
return handling;
}
-shared_ptr< vtk::VTKOutput> createFluidFieldVTKWriter( shared_ptr< StructuredBlockForest > & blocks,
- const BlockDataID & pdfFieldId, const BlockDataID & flagFieldId,
- std::string fname )
-{
- auto pdfFieldVTKWriter = vtk::createVTKOutput_BlockData( blocks, fname, uint_c(1), uint_c(0) );
+shared_ptr<vtk::VTKOutput> createFluidFieldVTKWriter(shared_ptr<StructuredBlockForest> &blocks,
+ const BlockDataID &pdfFieldId, const BlockDataID &flagFieldId,
+ std::string fname) {
+ auto pdfFieldVTKWriter = vtk::createVTKOutput_BlockData(blocks, fname, uint_c(1), uint_c(0));
- field::FlagFieldCellFilter< FlagField_T > fluidFilter( flagFieldId );
- fluidFilter.addFlag( Fluid_Flag );
- pdfFieldVTKWriter->addCellInclusionFilter( fluidFilter );
+ field::FlagFieldCellFilter<FlagField_T> fluidFilter(flagFieldId);
+ fluidFilter.addFlag(Fluid_Flag);
+ pdfFieldVTKWriter->addCellInclusionFilter(fluidFilter);
- auto velocityWriter = make_shared< lbm::VelocityVTKWriter< LatticeModel_T, float > >( pdfFieldId, "VelocityFromPDF" );
- pdfFieldVTKWriter->addCellDataWriter( velocityWriter );
+ auto velocityWriter = make_shared<lbm::VelocityVTKWriter<LatticeModel_T, float> >(pdfFieldId, "VelocityFromPDF");
+ pdfFieldVTKWriter->addCellDataWriter(velocityWriter);
return pdfFieldVTKWriter;
}
-class ResetPosition
-{
+class ResetPosition {
public:
-
- explicit ResetPosition( const shared_ptr<StructuredBlockStorage> & blockStorage,
- const BlockDataID & bodyStorageID,
- real_t dt )
- : blockStorage_( blockStorage )
- , bodyStorageID_( bodyStorageID )
- , dt_( dt )
- {}
-
- void operator()()
- {
- for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
- {
- for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- bodyIt->setPosition( bodyIt->getPosition() - bodyIt->getLinearVel() * dt_ );
+
+ explicit ResetPosition(const shared_ptr<StructuredBlockStorage> &blockStorage,
+ const BlockDataID &bodyStorageID,
+ real_t dt)
+ : blockStorage_(blockStorage), bodyStorageID_(bodyStorageID), dt_(dt) {}
+
+ void operator()() {
+ for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt) {
+ for (auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_);
+ bodyIt != pe::BodyIterator::end(); ++bodyIt) {
+ bodyIt->setPosition(bodyIt->getPosition() - bodyIt->getLinearVel() * dt_);
}
}
}
-
+
private:
-
+
shared_ptr<StructuredBlockStorage> blockStorage_;
- const BlockDataID & bodyStorageID_;
-
+ const BlockDataID &bodyStorageID_;
+
real_t dt_;
-
+
}; // class ResetPosition
-Vector3<real_t> analytical_velocity(Vector3<real_t> r, real_t b1, real_t b2, Vector3<real_t> e, real_t radius)
-{
+Vector3<real_t> analytical_velocity(Vector3<real_t> r, real_t b1, real_t b2, Vector3<real_t> e, real_t radius) {
auto r_len = r.length();
auto rs = r.getNormalized();
- auto frac = radius/r_len;
- auto frac2 = frac*frac;
- auto frac3 = frac2*frac;
- auto frac4 = frac3*frac;
+ auto frac = radius / r_len;
+ auto frac2 = frac * frac;
+ auto frac3 = frac2 * frac;
+ auto frac4 = frac3 * frac;
- return b1 * frac3 * ((e*rs) * rs - e / real_c(3.0)) + ( frac4 - frac2 ) * b2 * 0.5 * ( 3 * (e * rs) * (e * rs) - 1 ) * rs + frac4 * b2 * (e*rs) * ( (e*rs) * rs - e);
+ return b1 * frac3 * ((e * rs) * rs - e / real_c(3.0)) +
+ (frac4 - frac2) * b2 * 0.5 * (3 * (e * rs) * (e * rs) - 1) * rs + frac4 * b2 * (e * rs) * ((e * rs) * rs - e);
}
-int main( int argc, char **argv )
-{
+int main(int argc, char **argv) {
debug::enterTestMode();
- mpi::Environment env( argc, argv );
+ mpi::Environment env(argc, argv);
auto processes = MPIManager::instance()->numProcesses();
-
- if( processes != 1 && processes != 2 && processes != 4 && processes != 8)
- {
+
+ if (processes != 1 && processes != 2 && processes != 4 && processes != 8) {
std::cerr << "Number of processes must be equal to either 1, 2, 4, or 8!" << std::endl;
return EXIT_FAILURE;
}
-
+
bool shortrun = false;
bool writevtk = false;
- for( int i = 1; i < argc; ++i )
- {
- if( std::strcmp( argv[i], "--shortrun" ) == 0 ) { shortrun = true; continue; }
- if( std::strcmp( argv[i], "--vtk" ) == 0 ) { writevtk = true; continue; }
+ for (int i = 1; i < argc; ++i) {
+ if (std::strcmp(argv[i], "--shortrun") == 0) {
+ shortrun = true;
+ continue;
+ }
+ if (std::strcmp(argv[i], "--vtk") == 0) {
+ writevtk = true;
+ continue;
+ }
}
-
+
///////////////////////////
// SIMULATION PROPERTIES //
///////////////////////////
-
+
real_t R = 3;
uint_t L = 24;
uint_t T = 1000;
- if (shortrun)
- {
- L = 10;
- T = 10;
+ if (shortrun) {
+ L = 10;
+ T = 10;
}
const real_t dx = real_c(1.0);
const real_t dt = real_c(1.0);
- const real_t omega = lbm::collision_model::omegaFromViscosity( real_c(1.0/6.0) );
+ const real_t omega = lbm::collision_model::omegaFromViscosity(real_c(1.0 / 6.0));
const real_t v0 = real_c(0.01);
const real_t beta = real_c(5.0);
-
+
///////////////////////////
// BLOCK STRUCTURE SETUP //
///////////////////////////
- const uint_t XBlocks = (processes >= 2) ? uint_t( 2 ) : uint_t( 1 );
- const uint_t YBlocks = (processes >= 4) ? uint_t( 2 ) : uint_t( 1 );
- const uint_t ZBlocks = (processes == 8) ? uint_t( 2 ) : uint_t( 1 );
+ const uint_t XBlocks = (processes >= 2) ? uint_t(2) : uint_t(1);
+ const uint_t YBlocks = (processes >= 4) ? uint_t(2) : uint_t(1);
+ const uint_t ZBlocks = (processes == 8) ? uint_t(2) : uint_t(1);
const uint_t XCells = L / XBlocks;
const uint_t YCells = L / YBlocks;
const uint_t ZCells = L / ZBlocks;
// create fully periodic domain
- auto blocks = blockforest::createUniformBlockGrid( XBlocks, YBlocks, ZBlocks, XCells, YCells, ZCells, dx, true,
- true, true, true );
-
+ auto blocks = blockforest::createUniformBlockGrid(XBlocks, YBlocks, ZBlocks, XCells, YCells, ZCells, dx, true,
+ true, true, true);
+
////////
// PE //
////////
@@ -269,35 +271,39 @@ int main( int argc, char **argv )
auto globalBodyStorage = make_shared<pe::BodyStorage>();
pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
- auto ccdID = blocks->addBlockData(pe::ccd::createHashGridsDataHandling( globalBodyStorage, bodyStorageID ), "CCD");
- auto fcdID = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
-
+ auto ccdID = blocks->addBlockData(pe::ccd::createHashGridsDataHandling(globalBodyStorage, bodyStorageID), "CCD");
+ auto fcdID = blocks->addBlockData(
+ pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
+
pe::cr::HCSITS cr(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID);
-
+
/////////////////
// PE COUPLING //
/////////////////
// connect to pe
- const real_t overlap = real_c( 1.5 ) * dx;
+ const real_t overlap = real_c(1.5) * dx;
- if( R > real_c( L ) * real_c(0.5) - overlap )
- {
+ if (R > real_c(L) * real_c(0.5) - overlap) {
std::cerr << "Periodic sphere is too large and would lead to invalid PE state!" << std::endl;
return EXIT_FAILURE;
}
-
- boost::function<void(void)> syncCall = boost::bind( pe::syncShadowOwners<BodyTypeTuple>, boost::ref(blocks->getBlockForest()), bodyStorageID, static_cast<WcTimingTree*>(NULL), overlap, false );
-
- const auto myMat = pe::createMaterial( "myMat", real_c(1), real_t(0), real_t(1), real_t(1), real_t(0), real_t(1), real_t(1), real_t(0), real_t(0) );
+
+ boost::function<void(void)> syncCall = boost::bind(pe::syncShadowOwners<BodyTypeTuple>,
+ boost::ref(blocks->getBlockForest()), bodyStorageID,
+ static_cast<WcTimingTree *>(NULL), overlap, false);
+
+ const auto myMat = pe::createMaterial("myMat", real_c(1), real_t(0), real_t(1), real_t(1), real_t(0), real_t(1),
+ real_t(1), real_t(0), real_t(0));
// create the squirmer in the middle of the domain
- const Vector3<real_t> position (real_c(L)*real_c(0.5), real_c(L)*real_c(0.5), real_c(L)*real_c(0.5));
+ const Vector3<real_t> position(real_c(L) * real_c(0.5), real_c(L) * real_c(0.5), real_c(L) * real_c(0.5));
const Vector3<real_t> up(0.0, 0.0, 1.0);
const Vector3<real_t> orientation(1.0, 0.0, 0.0);
- const auto w = std::acos(up*orientation);
+ const auto w = std::acos(up * orientation);
math::Quaternion<real_t> q(up % orientation, w);
- auto squirmer = pe::createSquirmer(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, position, R, v0, beta, myMat);
+ auto squirmer = pe::createSquirmer(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, position, R, v0,
+ beta, myMat);
if (squirmer)
squirmer->setOrientation(q);
@@ -308,88 +314,96 @@ int main( int argc, char **argv )
////////////////////////
// create the lattice model
- LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ) );
+ LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::TRT::constructWithMagicNumber(omega));
// add PDF field ( uInit = <0,0,0>, rhoInit = 1 )
- BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel, FieldGhostLayers );
+ BlockDataID pdfFieldID = lbm::addPdfFieldToStorage<LatticeModel_T>(blocks, "pdf field (zyxf)", latticeModel,
+ FieldGhostLayers);
// add flag field
- BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" );
+ BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>(blocks, "flag field");
// add body field
- BlockDataID bodyFieldID = field::addToStorage< PeBodyField_T >( blocks, "body field", NULL, field::zyxf );
+ BlockDataID bodyFieldID = field::addToStorage<PeBodyField_T>(blocks, "body field", NULL, field::zyxf);
// add boundary handling & initialize outer domain boundaries
- BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >(
- MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
+ BlockDataID boundaryHandlingID = blocks->addStructuredBlockData<BoundaryHandling_T>(
+ MyBoundaryHandling(flagFieldID, pdfFieldID, bodyFieldID), "boundary handling");
///////////////
// TIME LOOP //
///////////////
// create the timeloop
- SweepTimeloop timeloop( blocks->getBlockStorage(), T );
-
+ SweepTimeloop timeloop(blocks->getBlockStorage(), T);
+
// sweep for updating the pe body mapping into the LBM simulation
timeloop.add()
- << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag, FormerMO_BB_Flag ),
- "Body Mapping" );
-
+ << Sweep(pe_coupling::BodyMapping<BoundaryHandling_T>(blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,
+ MO_BB_Flag, FormerMO_BB_Flag),
+ "Body Mapping");
+
// sweep for restoring PDFs in cells previously occupied by pe bodies
- typedef pe_coupling::EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T > Reconstructor_T;
- Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID );
+ typedef pe_coupling::EquilibriumReconstructor<LatticeModel_T, BoundaryHandling_T> Reconstructor_T;
+ Reconstructor_T reconstructor(blocks, boundaryHandlingID, pdfFieldID, bodyFieldID);
timeloop.add()
- << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,
- reconstructor, FormerMO_BB_Flag, Fluid_Flag ), "PDF Restore" );
+ << Sweep(pe_coupling::PDFReconstruction<LatticeModel_T, BoundaryHandling_T, Reconstructor_T>(blocks,
+ boundaryHandlingID,
+ bodyStorageID,
+ bodyFieldID,
+ reconstructor,
+ FormerMO_BB_Flag,
+ Fluid_Flag),
+ "PDF Restore");
// setup of the LBM communication for synchronizing the pdf field between neighboring blocks
- boost::function< void () > commFunction;
+ boost::function<void()> commFunction;
- blockforest::communication::UniformBufferedScheme< stencil::D3Q27 > scheme( blocks );
- scheme.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldID ) );
+ blockforest::communication::UniformBufferedScheme<stencil::D3Q27> scheme(blocks);
+ scheme.addPackInfo(make_shared<field::communication::PackInfo<PdfField_T> >(pdfFieldID));
commFunction = scheme;
// uses standard bounce back boundary conditions
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag );
-
- auto sweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
-
+ pe_coupling::mapMovingBodies<BoundaryHandling_T>(*blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,
+ MO_BB_Flag);
+
+ auto sweep = lbm::makeCellwiseSweep<LatticeModel_T, FlagField_T>(pdfFieldID, flagFieldID, Fluid_Flag);
+
// collision sweep
- timeloop.add() << Sweep( lbm::makeCollideSweep( sweep ), "cell-wise LB sweep (collide)" );
+ timeloop.add() << Sweep(lbm::makeCollideSweep(sweep), "cell-wise LB sweep (collide)");
// add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
- timeloop.add() << BeforeFunction( commFunction, "LBM Communication" )
- << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
+ timeloop.add() << BeforeFunction(commFunction, "LBM Communication")
+ << Sweep(BoundaryHandling_T::getBlockSweep(boundaryHandlingID), "Boundary Handling");
// streaming
- timeloop.add() << Sweep( lbm::makeStreamSweep( sweep ), "cell-wise LB sweep (stream)" );
-
+ timeloop.add() << Sweep(lbm::makeStreamSweep(sweep), "cell-wise LB sweep (stream)");
+
// add pe timesteps
- timeloop.addFuncAfterTimeStep( pe_coupling::TimeStep( blocks, bodyStorageID, cr, syncCall, dt ), "pe Time Step" );
- timeloop.addFuncAfterTimeStep( ResetPosition( blocks, bodyStorageID, dt ), "Reset pe positions" );
+ timeloop.addFuncAfterTimeStep(pe_coupling::TimeStep(blocks, bodyStorageID, cr, syncCall, dt), "pe Time Step");
+ timeloop.addFuncAfterTimeStep(ResetPosition(blocks, bodyStorageID, dt), "Reset pe positions");
- timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
+ timeloop.addFuncAfterTimeStep(RemainingTimeLogger(timeloop.getNrOfTimeSteps()), "Remaining Time Logger");
////////////////////////
// EXECUTE SIMULATION //
////////////////////////
-
+
WALBERLA_LOG_INFO_ON_ROOT("Starting test...");
WcTimingPool timeloopTiming;
- timeloop.run( timeloopTiming );
+ timeloop.run(timeloopTiming);
timeloopTiming.logResultOnRoot();
- if (writevtk)
- {
+ if (writevtk) {
WALBERLA_LOG_INFO_ON_ROOT("Writing out VTK file...");
- auto pdfFieldVTKWriter = createFluidFieldVTKWriter( blocks, pdfFieldID, flagFieldID, "fluid_field" );
- auto wf = vtk::writeFiles( pdfFieldVTKWriter );
+ auto pdfFieldVTKWriter = createFluidFieldVTKWriter(blocks, pdfFieldID, flagFieldID, "fluid_field");
+ auto wf = vtk::writeFiles(pdfFieldVTKWriter);
wf();
- auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", 1, 0 );
- flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID, "FlagField" ) );
- vtk::writeFiles( flagFieldVTK )();
+ auto flagFieldVTK = vtk::createVTKOutput_BlockData(blocks, "flag_field", 1, 0);
+ flagFieldVTK->addCellDataWriter(make_shared<field::VTKWriter<FlagField_T> >(flagFieldID, "FlagField"));
+ vtk::writeFiles(flagFieldVTK)();
}
-
+
if (squirmer) // this object only exists on the node that contains position
{
WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getSquirmerVelocity(), v0);
@@ -398,7 +412,7 @@ int main( int argc, char **argv )
WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getPosition(), position);
WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getQuaternion(), q);
}
-
+
WALBERLA_LOG_INFO_ON_ROOT("Checking result...");
auto b1 = real_c(1.5) * v0;
auto b2 = beta * b1;
@@ -409,66 +423,71 @@ int main( int argc, char **argv )
real_t abs_tolerance = real_c(0.0026);
real_t rel_tolerance = real_c(0.10);
- for ( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
- {
- IBlock & block = *blockIt;
- auto src = block.getData< PdfField_T >(pdfFieldID);
+ for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) {
+ IBlock &block = *blockIt;
+ auto src = block.getData<PdfField_T>(pdfFieldID);
- auto flagField = block.getData<FlagField_T> (flagFieldID);
+ auto flagField = block.getData<FlagField_T>(flagFieldID);
auto flag = flagField->getFlag(Fluid_Flag);
WALBERLA_FOR_ALL_CELLS_XYZ(src, {
Vector3<real_t> pos;
- blocks->getBlockLocalCellCenter( block, Cell(x,y,z), pos[0], pos[1], pos[2] );
+ blocks->getBlockLocalCellCenter(block, Cell(x, y, z), pos[0], pos[1], pos[2]);
- if( flagField->isFlagSet(x, y, z, flag) ) //check for fluid/non-boundary
+ if (flagField->isFlagSet(x, y, z, flag)) //check for fluid/non-boundary
{
auto r = pos - squirmer_pos;
-
+
auto v_ana = analytical_velocity(r, b1, b2, e, radius);
-
+
// Superposition with the nearest neighbours (faces)
- for( auto s = stencil::D3Q7::beginNoCenter(); s != stencil::D3Q7::end(); ++s )
- {
+ for (auto s = stencil::D3Q7::beginNoCenter(); s != stencil::D3Q7::end(); ++s) {
Vector3<real_t> imidpoint(squirmer_pos);
imidpoint[0] += real_c(s.cx() * int_c(L));
imidpoint[1] += real_c(s.cy() * int_c(L));
imidpoint[2] += real_c(s.cz() * int_c(L));
- v_ana += analytical_velocity( pos - imidpoint, b1, b2, e, radius );
+ v_ana += analytical_velocity(pos - imidpoint, b1, b2, e, radius);
}
-
- if( !shortrun && r.length() > 2*radius ) // Exclude the volume around the squirmer, as the analytical solution is a far field solution only
+
+ if (!shortrun && r.length() > 2 *
+ radius) // Exclude the volume around the squirmer, as the analytical solution is a far field solution only
{
auto v_data = src->getVelocity(x, y, z);
auto diff = v_data - v_ana;
-
- for( uint_t d = 0; d < 3; ++d )
- {
- if( std::abs( diff[d] / v_ana[d] ) > rel_tolerance && std::abs( diff[d] ) > abs_tolerance )
- {
- WALBERLA_LOG_DEVEL("Difference too large in " << Cell(x,y,z) << " (" << r.length() << "). Expected " << v_ana << ", got " << v_data << ", relative " << std::abs(diff[d] / v_ana[d]) << ", absolute " << std::abs(diff[d]));
+
+ for (uint_t d = 0; d < 3; ++d) {
+ if (std::abs(diff[d] / v_ana[d]) > rel_tolerance && std::abs(diff[d]) > abs_tolerance) {
+ WALBERLA_LOG_DEVEL(
+ "Difference too large in " << Cell(x, y, z) << " (" << r.length() << "). Expected " << v_ana
+ << ", got " << v_data << ", relative "
+ << std::abs(diff[d] / v_ana[d]) << ", absolute "
+ << std::abs(diff[d]));
}
}
}
-
- if( writevtk )
- {
+
+ if (writevtk) {
// store the value into the PDF field so we can write it to VTK later
- src->setToEquilibrium( x, y, z, v_ana );
+ src->setToEquilibrium(x, y, z, v_ana);
}
}
});
}
-
- if (writevtk)
- {
+
+ if (writevtk) {
WALBERLA_LOG_INFO_ON_ROOT("Writing out reference VTK file...");
- auto pdfFieldVTKWriter = createFluidFieldVTKWriter( blocks, pdfFieldID, flagFieldID, "fluid_field_ref" );
- auto wf = vtk::writeFiles( pdfFieldVTKWriter );
+ auto pdfFieldVTKWriter = createFluidFieldVTKWriter(blocks, pdfFieldID, flagFieldID, "fluid_field_ref");
+ auto wf = vtk::writeFiles(pdfFieldVTKWriter);
wf();
}
-
+
WALBERLA_LOG_INFO_ON_ROOT("Completed test.");
-
+
return 0;
}
+
+} // namespace squirmer
+
+int main( int argc, char **argv ){
+ squirmer::main(argc, argv);
+}
\ No newline at end of file
diff --git a/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp b/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
index b7d4fc43c..f4c2df525 100644
--- a/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
@@ -339,7 +339,7 @@ int main( int argc, char **argv )
// map pe bodies into the LBM simulation
// moving bodies are handled by the momentum exchange method, thus act here as velocity boundary condition
- pe_coupling::mapMovingGlobalBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, *globalBodyStorage, bodyFieldID, MO_Flag );
+ pe_coupling::mapMovingGlobalBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, *globalBodyStorage, bodyFieldID, MO_Flag );
///////////////
// TIME LOOP //
diff --git a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
index 4e41079db..96bb5a31e 100644
--- a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
@@ -13,7 +13,7 @@
// 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 TorqueSphereMEMPe.cpp
+//! \file TorqueSphereMEM.cpp
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
@@ -64,7 +64,7 @@
#include <iomanip>
#include <iostream>
-namespace torque_sphere_mem_pe
+namespace torque_sphere_mem
{
@@ -476,13 +476,13 @@ int main( int argc, char **argv )
if( MO_CLI )
{
// uses a higher order boundary condition (CLI)
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
}else if ( MO_MR ){
// uses a higher order boundary condition (MR)
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_MR_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_MR_Flag );
}else{
// uses standard bounce back boundary conditions
- pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag );
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag );
}
shared_ptr< TorqueEval > torqueEval = make_shared< TorqueEval >( &timeloop, &setup, blocks, bodyStorageID, fileIO );
@@ -563,8 +563,8 @@ int main( int argc, char **argv )
}
-} //namespace torque_sphere_mem_pe
+} //namespace torque_sphere_mem
int main( int argc, char **argv ){
- torque_sphere_mem_pe::main(argc, argv);
+ torque_sphere_mem::main(argc, argv);
}
diff --git a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp
index 8444242fe..79e3c2f49 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSM.cpp
@@ -13,7 +13,7 @@
// 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 DragForceSpherePSMPe.cpp
+//! \file DragForceSpherePSM.cpp
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
@@ -58,12 +58,13 @@
#include "pe/basic.h"
#include "pe_coupling/partially_saturated_cells_method/all.h"
+#include "pe_coupling/utility/all.h"
#include <vector>
#include <iomanip>
#include <iostream>
-namespace drag_force_sphere_psm_pe
+namespace drag_force_sphere_psm
{
///////////
@@ -252,12 +253,12 @@ private:
{
for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
{
- velocity_sum += pe_coupling::getPSMMacroscopicVelocity<LatticeModel_T, 1>( *blockIt, pdfField, bodyAndVolumeFractionField, blocks_, *fieldIt )[0];
+ velocity_sum += pe_coupling::getPSMMacroscopicVelocity<LatticeModel_T, 1>( *blockIt, pdfField, bodyAndVolumeFractionField, *blocks_, *fieldIt )[0];
}
}else{
for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
{
- velocity_sum += pe_coupling::getPSMMacroscopicVelocity<LatticeModel_T, 2>( *blockIt, pdfField, bodyAndVolumeFractionField, blocks_, *fieldIt )[0];
+ velocity_sum += pe_coupling::getPSMMacroscopicVelocity<LatticeModel_T, 2>( *blockIt, pdfField, bodyAndVolumeFractionField, *blocks_, *fieldIt )[0];
}
}
}
@@ -287,30 +288,6 @@ private:
};
-class ResetForce
-{
-public:
- ResetForce( const shared_ptr< StructuredBlockStorage > & blocks,
- const BlockDataID & bodyStorageID )
- : blocks_( blocks ), bodyStorageID_( bodyStorageID )
- { }
-
- void operator()()
- {
- for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
- {
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- bodyIt->resetForceAndTorque();
- }
- }
- }
-private:
- shared_ptr< StructuredBlockStorage > blocks_;
- const BlockDataID bodyStorageID_;
-};
-
-
//////////
// MAIN //
//////////
@@ -463,11 +440,11 @@ int main( int argc, char **argv )
// initialize the PDF field for PSM
if( method == PSMVariant::SC1W1 || method == PSMVariant::SC2W1 || method == PSMVariant::SC3W1 )
{
- pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
+ pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
}
else
{
- pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
+ pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
}
///////////////
@@ -553,7 +530,7 @@ int main( int argc, char **argv )
timeloop.addFuncAfterTimeStep( SharedFunctor< DragForceEvaluator >(forceEval), "drag force evaluation" );
// resetting force
- timeloop.addFuncAfterTimeStep( ResetForce( blocks, bodyStorageID ), "reset force on sphere");
+ timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesResetter( blocks, bodyStorageID ), "reset force on sphere");
timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
@@ -600,8 +577,8 @@ int main( int argc, char **argv )
return 0;
}
-} //namespace drag_force_sphere_psm_pe
+} //namespace drag_force_sphere_psm
int main( int argc, char **argv ){
- drag_force_sphere_psm_pe::main(argc, argv);
+ drag_force_sphere_psm::main(argc, argv);
}
diff --git a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
index 329c56767..f29f7d650 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
@@ -13,7 +13,7 @@
// 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 DragForceSpherePSMRefinementPe.cpp
+//! \file DragForceSpherePSMRefinement.cpp
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
@@ -61,6 +61,7 @@
#include "pe/vtk/SphereVtkOutput.h"
#include "pe_coupling/partially_saturated_cells_method/all.h"
+#include "pe_coupling/utility/all.h"
#include "vtk/all.h"
#include "field/vtk/all.h"
@@ -70,7 +71,7 @@
#include <iomanip>
#include <iostream>
-namespace drag_force_sphere_psm_pe_refinement
+namespace drag_force_sphere_psm_refinement
{
///////////
@@ -371,12 +372,12 @@ private:
{
for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
{
- velocity_sum += cellVolume * pe_coupling::getPSMMacroscopicVelocity<LatticeModel_T, 1>( *blockIt, pdfField, bodyAndVolumeFractionField, blocks_, *fieldIt )[0];
+ velocity_sum += cellVolume * pe_coupling::getPSMMacroscopicVelocity<LatticeModel_T, 1>( *blockIt, pdfField, bodyAndVolumeFractionField, *blocks_, *fieldIt )[0];
}
}else{
for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
{
- velocity_sum += cellVolume * pe_coupling::getPSMMacroscopicVelocity<LatticeModel_T, 2>( *blockIt, pdfField, bodyAndVolumeFractionField, blocks_, *fieldIt )[0];
+ velocity_sum += cellVolume * pe_coupling::getPSMMacroscopicVelocity<LatticeModel_T, 2>( *blockIt, pdfField, bodyAndVolumeFractionField, *blocks_, *fieldIt )[0];
}
}
}
@@ -408,30 +409,6 @@ private:
};
-class ResetForce
-{
- public:
- ResetForce( const shared_ptr< StructuredBlockStorage > & blocks,
- const BlockDataID & bodyStorageID )
- : blocks_( blocks ), bodyStorageID_( bodyStorageID )
- { }
-
- void operator()()
- {
- for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
- {
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- bodyIt->resetForceAndTorque();
- }
- }
- }
-private:
- shared_ptr< StructuredBlockStorage > blocks_;
- const BlockDataID bodyStorageID_;
-};
-
-
//////////
// MAIN //
//////////
@@ -581,11 +558,11 @@ int main( int argc, char **argv )
// initialize the PDF field for PSM
if( method == PSMVariant::SC1W1 || method == PSMVariant::SC2W1 || method == PSMVariant::SC3W1 )
{
- pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
+ pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
}
else
{
- pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
+ pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
}
///////////////
// TIME LOOP //
@@ -632,7 +609,7 @@ int main( int argc, char **argv )
timeloop.addFuncAfterTimeStep( SharedFunctor< DragForceEvaluator >(forceEval), "drag force evaluation" );
// resetting force
- timeloop.addFuncAfterTimeStep( ResetForce( blocks, bodyStorageID ), "reset force on sphere");
+ timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesResetter( blocks, bodyStorageID ), "reset force on sphere");
if( vtkIO )
{
@@ -702,8 +679,8 @@ int main( int argc, char **argv )
}
-} //namespace drag_force_sphere_psm_pe_refinement
+} //namespace drag_force_sphere_psm_refinement
int main( int argc, char **argv ){
- drag_force_sphere_psm_pe_refinement::main(argc, argv);
+ drag_force_sphere_psm_refinement::main(argc, argv);
}
diff --git a/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
index ef3378bca..f05d9b43d 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
@@ -13,7 +13,7 @@
// 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 SegreSilberbergPSMPe.cpp
+//! \file SegreSilberbergPSM.cpp
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
@@ -64,7 +64,7 @@
#include "field/vtk/all.h"
#include "lbm/vtk/all.h"
-namespace segre_silberberg_psm_pe
+namespace segre_silberberg_psm
{
///////////
@@ -353,71 +353,6 @@ class SteadyStateCheck
real_t posNew_;
};
-// When using N LBM steps and one (larger) pe step in a single simulation step, the force applied on the pe bodies in each LBM sep has to be scaled
-// by a factor of 1/N before running the pe simulation. This corresponds to an averaging of the force and torque over the N LBM steps and is said to
-// damp oscillations. Usually, N = 2.
-// See Ladd - " Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation", 1994, p. 302
-class AverageForce
-{
- public:
- AverageForce( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & bodyStorageID,
- const uint_t lbmSubCycles )
- : blocks_( blocks ), bodyStorageID_( bodyStorageID ), invLbmSubCycles_( real_c(1) / real_c( lbmSubCycles ) ) { }
-
- void operator()()
- {
- pe::Vec3 force = pe::Vec3(0.0);
- pe::Vec3 torque = pe::Vec3(0.0);
-
- for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
- {
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- force = invLbmSubCycles_ * bodyIt->getForce();
- torque = invLbmSubCycles_ * bodyIt->getTorque();
-
- bodyIt->resetForceAndTorque();
-
- bodyIt->setForce ( force );
- bodyIt->setTorque( torque );
- }
- }
- }
- private:
- shared_ptr< StructuredBlockStorage > blocks_;
- const BlockDataID bodyStorageID_;
- real_t invLbmSubCycles_;
-};
-
-
-// The flow in the channel is here driven by a body force and is periodic
-// Thus, the pressure gradient, usually encountered in a channel flow, is not present here
-// The buoyancy force on the body due to this pressure gradient has to be added 'artificially'
-// F_{buoyancy} = - V_{body} * grad ( p ) = V_{body} * \rho_{fluid} * a
-// ( V_{body} := volume of body, a := acceleration driving the flow )
-class BuoyancyForce
-{
- public:
- BuoyancyForce( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & bodyStorageID,
- const Vector3<real_t> pressureForce )
- : blocks_( blocks ), bodyStorageID_( bodyStorageID ), pressureForce_( pressureForce ) { }
-
- void operator()()
- {
-
- for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
- {
- for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
- {
- bodyIt->addForce ( pressureForce_ );
- }
- }
- }
- private:
- shared_ptr< StructuredBlockStorage > blocks_;
- const BlockDataID bodyStorageID_;
- const Vector3<real_t> pressureForce_;
-};
//////////
// MAIN //
@@ -621,11 +556,11 @@ int main( int argc, char **argv )
// initialize the PDF field for PSM
if( SC1W1 || SC2W1 || SC3W1 )
{
- pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
+ pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
}
else
{
- pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
+ pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
}
///////////////
@@ -679,13 +614,17 @@ int main( int argc, char **argv )
}
// average the forces acting on the bodies over the number of LBM steps
- timeloop.addFuncAfterTimeStep( AverageForce( blocks, bodyStorageID, numLbmSubCycles ), "Force averaging for several LBM steps" );
+ timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesScaler( blocks, bodyStorageID, real_t(1)/real_c(numLbmSubCycles) ), "Force averaging for several LBM steps" );
// add pressure force contribution
- timeloop.addFuncAfterTimeStep( BuoyancyForce( blocks, bodyStorageID,
- Vector3<real_t>( math::M_PI / real_c(6) * setup.forcing * setup.particleDiam * setup.particleDiam * setup.particleDiam ,
- real_c(0), real_c(0) ) ),
- "Buoyancy force" );
+ // The flow in the channel is here driven by a body force and is periodic
+ // Thus, the pressure gradient, usually encountered in a channel flow, is not present here
+ // The buoyancy force on the body due to this pressure gradient has to be added 'artificially'
+ // F_{buoyancy} = - V_{body} * grad ( p ) = V_{body} * \rho_{fluid} * a
+ // ( V_{body} := volume of body, a := acceleration driving the flow )
+ Vector3<real_t> buoyancyForce(math::M_PI / real_t(6) * setup.forcing * setup.particleDiam * setup.particleDiam * setup.particleDiam ,
+ real_t(0), real_t(0));
+ timeloop.addFuncAfterTimeStep( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, buoyancyForce ), "Buoyancy force" );
// add pe timesteps
timeloop.addFuncAfterTimeStep( pe_coupling::TimeStep( blocks, bodyStorageID, cr, syncCall, dt_pe, pe_interval ), "pe Time Step" );
@@ -775,8 +714,8 @@ int main( int argc, char **argv )
return EXIT_SUCCESS;
}
-} // namespace segre_silberberg_psm_pe
+} // namespace segre_silberberg_psm
int main( int argc, char **argv ){
- segre_silberberg_psm_pe::main(argc, argv);
+ segre_silberberg_psm::main(argc, argv);
}
diff --git a/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp
index 75fd18cba..31921f2f0 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/TorqueSpherePSM.cpp
@@ -13,7 +13,7 @@
// 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 TorqueSpherePSMPe.cpp
+//! \file TorqueSpherePSM.cpp
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
@@ -57,12 +57,13 @@
#include "pe/basic.h"
#include "pe_coupling/partially_saturated_cells_method/all.h"
+#include "pe_coupling/utility/all.h"
#include <vector>
#include <iomanip>
#include <iostream>
-namespace torque_sphere_psm_pe
+namespace torque_sphere_psm
{
@@ -219,36 +220,10 @@ class TorqueEval
};
-class ResetForce
-{
- public:
- ResetForce( const shared_ptr< StructuredBlockStorage > & blocks,
- const BlockDataID & bodyStorageID )
- : blocks_( blocks ), bodyStorageID_( bodyStorageID )
- { }
-
- void operator()()
- {
- for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
- {
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- bodyIt->resetForceAndTorque();
- }
- }
- }
-
- private:
- shared_ptr< StructuredBlockStorage > blocks_;
- const BlockDataID bodyStorageID_;
-};
-
-
//////////
// MAIN //
//////////
-
//*******************************************************************************************************************
/*!\brief Testcase that checks the torque acting on a constantly rotating sphere in the center of a cubic domain
*
@@ -273,7 +248,6 @@ class ResetForce
*/
//*******************************************************************************************************************
-
int main( int argc, char **argv )
{
debug::enterTestMode();
@@ -413,11 +387,11 @@ int main( int argc, char **argv )
// initialize the PDF field for PSM
if( SC1W1 || SC2W1 || SC3W1 )
{
- pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
+ pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
}
else
{
- pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
+ pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
}
///////////////
@@ -481,7 +455,7 @@ int main( int argc, char **argv )
}
// resetting force
- timeloop.addFuncAfterTimeStep( ResetForce( blocks, bodyStorageID ), "reset force on sphere");
+ timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesResetter( blocks, bodyStorageID ), "reset force on sphere");
timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
@@ -528,8 +502,8 @@ int main( int argc, char **argv )
}
-} //namespace torque_sphere_psm_pe
+} //namespace torque_sphere_psm
int main( int argc, char **argv ){
- torque_sphere_psm_pe::main(argc, argv);
+ torque_sphere_psm::main(argc, argv);
}
--
GitLab