From 25cc89e0d2f469243cb9b8d52774d1f14cdbd98a Mon Sep 17 00:00:00 2001
From: Michael Kuron <mkuron@icp.uni-stuttgart.de>
Date: Thu, 18 Jul 2019 22:26:03 +0200
Subject: [PATCH] [API] Add momentum conservation to the Momentum Exchange
Method and hook for PDF destruction
---
.../AMRSedimentSettling.cpp | 6 +-
.../WorkloadEvaluation.cpp | 6 +-
.../MotionSingleHeavySphere.cpp | 10 +--
.../momentum_exchange_method/BodyMapping.h | 48 ++++++++--
.../momentum_exchange_method/all.h | 1 +
.../destruction/Destroyer.h | 45 ++++++++++
.../destruction/all.h | 25 ++++++
.../restoration/PDFReconstruction.h | 38 +++++---
.../restoration/Reconstructor.h | 87 ++++++++-----------
.../BodyAtBlockBoarderCheck.cpp | 6 +-
.../BodyMappingTest.cpp | 4 +-
.../LubricationCorrectionMEM.cpp | 6 +-
.../PeriodicParticleChannelMEM.cpp | 6 +-
.../SegreSilberbergMEM.cpp | 18 ++--
.../SettlingSphereMEM.cpp | 68 ++++++++++-----
.../SettlingSphereMEMDynamicRefinement.cpp | 6 +-
.../SettlingSphereMEMStaticRefinement.cpp | 6 +-
.../momentum_exchange_method/SquirmerTest.cpp | 7 +-
18 files changed, 267 insertions(+), 126 deletions(-)
create mode 100644 src/pe_coupling/momentum_exchange_method/destruction/Destroyer.h
create mode 100644 src/pe_coupling/momentum_exchange_method/destruction/all.h
diff --git a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp
index 750a8e964..e7a4a0fea 100644
--- a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp
+++ b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp
@@ -1840,13 +1840,13 @@ int main( int argc, char **argv )
"pe Time Step", finestLevel );
// add sweep for updating the pe body mapping into the LBM simulation
- refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), blocks ),
+ refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), blocks ),
"Body Mapping", finestLevel );
// add 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 );
- refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks,
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID );
+ refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks, pdfFieldID,
boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), blocks ),
"PDF Restore", finestLevel );
diff --git a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp
index 588007edf..7b87d670a 100644
--- a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp
+++ b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp
@@ -751,14 +751,14 @@ int main( int argc, char **argv )
// sweep for updating the pe body mapping into the LBM simulation
timeloop.add()
- << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_CLI_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
+ << Sweep( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_CLI_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "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);
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID);
timeloop.add()
<< Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >
- ( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag, pe_coupling::selectRegularBodies, optimizeForSmallObstacleFraction ), "PDF Restore" );
+ ( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag, pe_coupling::selectRegularBodies, optimizeForSmallObstacleFraction ), "PDF Restore" );
shared_ptr<pe_coupling::BodiesForceTorqueContainer> bodiesFTContainer1 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
std::function<void(void)> storeForceTorqueInCont1 = std::bind(&pe_coupling::BodiesForceTorqueContainer::store, bodiesFTContainer1);
diff --git a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
index de0647508..725ff77b5 100644
--- a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
+++ b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
@@ -1195,20 +1195,20 @@ int main( int argc, char **argv )
// sweep for updating the pe body mapping into the LBM simulation
if( memVariant == MEMVariant::CLI )
- timeloop.add() << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MEM_CLI_Flag, FormerMEM_Flag ), "Body Mapping" );
+ timeloop.add() << Sweep( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MEM_CLI_Flag, FormerMEM_Flag ), "Body Mapping" );
else if ( memVariant == MEMVariant::MR )
- timeloop.add() << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MEM_MR_Flag, FormerMEM_Flag ), "Body Mapping" );
+ timeloop.add() << Sweep( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MEM_MR_Flag, FormerMEM_Flag ), "Body Mapping" );
else
- timeloop.add() << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MEM_BB_Flag, FormerMEM_Flag ), "Body Mapping" );
+ timeloop.add() << Sweep( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MEM_BB_Flag, FormerMEM_Flag ), "Body Mapping" );
// reconstruct missing PDFs
using ExtrapolationFinder_T = pe_coupling::SphereNormalExtrapolationDirectionFinder;
ExtrapolationFinder_T extrapolationFinder( blocks, bodyFieldID );
typedef pe_coupling::ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationFinder_T > Reconstructor_T;
- Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder, true );
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID, extrapolationFinder, true );
timeloop.add() << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >
- ( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMEM_Flag, Fluid_Flag ), "PDF Restore" );
+ ( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMEM_Flag, Fluid_Flag ), "PDF Restore" );
for( uint_t lbmcycle = 0; lbmcycle < numLBMSubCycles; ++lbmcycle ){
diff --git a/src/pe_coupling/momentum_exchange_method/BodyMapping.h b/src/pe_coupling/momentum_exchange_method/BodyMapping.h
index 045eeb5d1..4c46bc729 100644
--- a/src/pe_coupling/momentum_exchange_method/BodyMapping.h
+++ b/src/pe_coupling/momentum_exchange_method/BodyMapping.h
@@ -55,35 +55,55 @@ namespace pe_coupling {
*
* The 'mappingBodySelectorFct' can be used to decide which bodies should be mapped or not.
*/
-template< typename BoundaryHandling_T >
+template< typename LatticeModel_T, typename BoundaryHandling_T, typename Destroyer_T = NaNDestroyer<LatticeModel_T>, bool conserveMomentum = false >
class BodyMapping
{
public:
+ typedef lbm::PdfField< LatticeModel_T > PdfField_T;
typedef typename BoundaryHandling_T::FlagField FlagField_T;
typedef typename BoundaryHandling_T::flag_t flag_t;
typedef Field< pe::BodyID, 1 > BodyField_T;
BodyMapping( const shared_ptr<StructuredBlockStorage> & blockStorage,
+ const BlockDataID & pdfFieldID,
const BlockDataID & boundaryHandlingID,
const BlockDataID & bodyStorageID,
const shared_ptr<pe::BodyStorage> & globalBodyStorage,
const BlockDataID & bodyFieldID,
+ const Destroyer_T & destroyer,
const FlagUID & obstacle, const FlagUID & formerObstacle,
const std::function<bool(pe::BodyID)> & mappingBodySelectorFct = selectRegularBodies )
- : blockStorage_( blockStorage ), boundaryHandlingID_( boundaryHandlingID ),
+ : blockStorage_( blockStorage ), pdfFieldID_( pdfFieldID ), boundaryHandlingID_( boundaryHandlingID ),
bodyStorageID_(bodyStorageID), globalBodyStorage_( globalBodyStorage ), bodyFieldID_( bodyFieldID ),
- obstacle_( obstacle ), formerObstacle_( formerObstacle ), mappingBodySelectorFct_( mappingBodySelectorFct )
+ destroyer_( destroyer ), obstacle_( obstacle ), formerObstacle_( formerObstacle ),
+ mappingBodySelectorFct_( mappingBodySelectorFct )
+ {}
+
+ BodyMapping( const shared_ptr<StructuredBlockStorage> & blockStorage,
+ const BlockDataID & pdfFieldID,
+ const BlockDataID & boundaryHandlingID,
+ const BlockDataID & bodyStorageID,
+ const shared_ptr<pe::BodyStorage> & globalBodyStorage,
+ const BlockDataID & bodyFieldID,
+ const FlagUID & obstacle, const FlagUID & formerObstacle,
+ const std::function<bool(pe::BodyID)> & mappingBodySelectorFct = selectRegularBodies )
+ : blockStorage_( blockStorage ), pdfFieldID_( pdfFieldID ), boundaryHandlingID_( boundaryHandlingID ),
+ bodyStorageID_(bodyStorageID), globalBodyStorage_( globalBodyStorage ), bodyFieldID_( bodyFieldID ),
+ destroyer_( Destroyer_T() ), obstacle_( obstacle ), formerObstacle_( formerObstacle ),
+ mappingBodySelectorFct_( mappingBodySelectorFct )
{}
void operator()( IBlock * const block )
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
+ PdfField_T * pdfField = block->getData< PdfField_T >( pdfFieldID_ );
BoundaryHandling_T * boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingID_ );
FlagField_T * flagField = boundaryHandling->getFlagField();
BodyField_T * bodyField = block->getData< BodyField_T >( bodyFieldID_ );
+ WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandling );
WALBERLA_ASSERT_NOT_NULLPTR( flagField );
WALBERLA_ASSERT_NOT_NULLPTR( bodyField );
@@ -103,19 +123,19 @@ public:
for( auto bodyIt = pe::BodyIterator::begin(*block, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
- mapBodyAndUpdateMapping(bodyIt.getBodyID(), block, boundaryHandling, flagField , bodyField, obstacle, formerObstacle, dx, dy, dz);
+ mapBodyAndUpdateMapping(bodyIt.getBodyID(), block, pdfField, boundaryHandling, flagField , bodyField, obstacle, formerObstacle, dx, dy, dz);
}
for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt)
{
if( mappingBodySelectorFct_(bodyIt.getBodyID()))
- mapBodyAndUpdateMapping(bodyIt.getBodyID(), block, boundaryHandling, flagField , bodyField, obstacle, formerObstacle, dx, dy, dz);
+ mapBodyAndUpdateMapping(bodyIt.getBodyID(), block, pdfField, boundaryHandling, flagField , bodyField, obstacle, formerObstacle, dx, dy, dz);
}
}
private:
void mapBodyAndUpdateMapping(pe::BodyID body, IBlock * const block,
- BoundaryHandling_T * boundaryHandling, FlagField_T * flagField, BodyField_T * bodyField,
+ PdfField_T * pdfField, BoundaryHandling_T * boundaryHandling, FlagField_T * flagField, BodyField_T * bodyField,
const flag_t & obstacle, const flag_t & formerObstacle,
real_t dx, real_t dy, real_t dz)
{
@@ -156,6 +176,19 @@ private:
{
// set obstacle flag
boundaryHandling->forceBoundary( obstacle, x, y, z );
+
+ if( conserveMomentum && pdfField->isInInnerPart(Cell(x,y,z)) ) {
+ // this removes the fluid information (PDFs) from the simulation, together with the containing momentum
+ // to ensure momentum conservation, the momentum of this fluid cell has to be added to the particle
+ // see Aidun, C. K., Lu, Y., & Ding, E. J. (1998). Direct analysis of particulate suspensions with inertia using the discrete Boltzmann equation. Journal of Fluid Mechanics, 373, 287-311.
+
+ // force = momentum / dt, with dt = 1
+ Vector3<real_t> momentum = pdfField->getMomentumDensity(x, y, z);
+ body->addForceAtPos(momentum[0], momentum[1], momentum[2], cx, cy, cz);
+ }
+
+ // invalidate PDF values
+ destroyer_( x, y, z, block, pdfField );
}
}
// let pointer from body field point to this body
@@ -190,11 +223,14 @@ private:
shared_ptr<StructuredBlockStorage> blockStorage_;
+ const BlockDataID pdfFieldID_;
const BlockDataID boundaryHandlingID_;
const BlockDataID bodyStorageID_;
shared_ptr<pe::BodyStorage> globalBodyStorage_;
const BlockDataID bodyFieldID_;
+ Destroyer_T destroyer_;
+
const FlagUID obstacle_;
const FlagUID formerObstacle_;
diff --git a/src/pe_coupling/momentum_exchange_method/all.h b/src/pe_coupling/momentum_exchange_method/all.h
index a482e213d..5b0bfaf23 100644
--- a/src/pe_coupling/momentum_exchange_method/all.h
+++ b/src/pe_coupling/momentum_exchange_method/all.h
@@ -23,5 +23,6 @@
#pragma once
#include "boundary/all.h"
+#include "destruction/all.h"
#include "restoration/all.h"
#include "BodyMapping.h"
diff --git a/src/pe_coupling/momentum_exchange_method/destruction/Destroyer.h b/src/pe_coupling/momentum_exchange_method/destruction/Destroyer.h
new file mode 100644
index 000000000..68d9ae2a3
--- /dev/null
+++ b/src/pe_coupling/momentum_exchange_method/destruction/Destroyer.h
@@ -0,0 +1,45 @@
+//======================================================================================================================
+//
+// 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 Destroyer.h
+//! \ingroup pe_coupling
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//
+//======================================================================================================================
+
+
+#pragma once
+
+#include <limits>
+
+
+namespace walberla {
+namespace pe_coupling {
+
+template< typename LatticeModel_T >
+class NaNDestroyer
+{
+public:
+
+ typedef lbm::PdfField< LatticeModel_T > PdfField_T;
+
+ void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const /*block*/, PdfField_T * const pdfField ) {
+ for (auto d = uint_t(0); d < LatticeModel_T::Stencil::Size; ++d)
+ pdfField->get(x, y, z, d) = std::numeric_limits<typename PdfField_T::value_type>::quiet_NaN();
+ }
+};
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/momentum_exchange_method/destruction/all.h b/src/pe_coupling/momentum_exchange_method/destruction/all.h
new file mode 100644
index 000000000..9b2149a6a
--- /dev/null
+++ b/src/pe_coupling/momentum_exchange_method/destruction/all.h
@@ -0,0 +1,25 @@
+//======================================================================================================================
+//
+// 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 all.h
+//! \ingroup pe_coupling
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//! \brief Collective header file for module pe_coupling moving obstacle destruction
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "Destroyer.h"
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h b/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
index e65d2cc3c..0d6dcfe12 100644
--- a/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
+++ b/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
@@ -55,16 +55,18 @@ namespace pe_coupling {
*/
//**************************************************************************************************************************************
-template< typename LatticeModel_T, typename BoundaryHandling_T, typename Reconstructer_T >
+template< typename LatticeModel_T, typename BoundaryHandling_T, typename Reconstructer_T, bool conserveMomentum = false, bool includeGhostLayers = false >
class PDFReconstruction
{
public:
+ typedef lbm::PdfField< LatticeModel_T > PdfField_T;
typedef typename BoundaryHandling_T::FlagField FlagField_T;
typedef typename BoundaryHandling_T::flag_t flag_t;
typedef Field< pe::BodyID, 1 > BodyField_T;
inline PDFReconstruction( const shared_ptr<StructuredBlockStorage> & blockStorage,
+ const BlockDataID & pdfFieldID,
const BlockDataID & boundaryHandlingID,
const BlockDataID & bodyStorageID,
const shared_ptr<pe::BodyStorage> & globalBodyStorage,
@@ -73,7 +75,8 @@ public:
const FlagUID & formerObstacle, const FlagUID & fluid,
const std::function<bool(pe::BodyID)> & movingBodySelectorFct = selectRegularBodies,
const bool optimizeForSmallObstacleFraction = false ) :
- blockStorage_( blockStorage ), boundaryHandlingID_( boundaryHandlingID ), bodyStorageID_(bodyStorageID),
+ blockStorage_( blockStorage ), pdfFieldID_( pdfFieldID ),
+ boundaryHandlingID_( boundaryHandlingID ), bodyStorageID_(bodyStorageID),
globalBodyStorage_( globalBodyStorage ), bodyFieldID_( bodyFieldID ),
reconstructor_ ( reconstructor ), formerObstacle_( formerObstacle ), fluid_( fluid ),
movingBodySelectorFct_( movingBodySelectorFct ),
@@ -85,13 +88,23 @@ public:
private:
void reconstructPDFsInCells( const CellInterval & cells, IBlock * const block,
- FlagField_T * flagField, const flag_t & formerObstacle )
+ PdfField_T * pdfField, FlagField_T * flagField, BodyField_T * bodyField,
+ const flag_t & formerObstacle )
{
for( cell_idx_t z = cells.zMin(); z <= cells.zMax(); ++z ) {
for (cell_idx_t y = cells.yMin(); y <= cells.yMax(); ++y) {
for (cell_idx_t x = cells.xMin(); x <= cells.xMax(); ++x) {
if (isFlagSet(flagField->get(x,y,z), formerObstacle)) {
- reconstructor_(x,y,z,block);
+ reconstructor_(x,y,z,block,pdfField);
+
+ if( conserveMomentum && pdfField->isInInnerPart(Cell(x,y,z)) ) {
+ // the (artificially) added momentum in the restored fluid cell has to be subtracted from the former particle to ensure momentum conservation
+ Vector3<real_t> momentum = pdfField->getMomentumDensity(x,y,z);
+
+ // force = momentum / dt, with dt = 1
+ Vector3< real_t > cellCenter = blockStorage_->getBlockLocalCellCenter(*block, Cell(x,y,z) );
+ (*bodyField)(x,y,z)->addForceAtPos(-momentum, cellCenter);
+ }
}
}
}
@@ -117,6 +130,7 @@ private:
shared_ptr<StructuredBlockStorage> blockStorage_;
+ const BlockDataID pdfFieldID_;
const BlockDataID boundaryHandlingID_;
const BlockDataID bodyStorageID_;
shared_ptr<pe::BodyStorage> globalBodyStorage_;
@@ -134,16 +148,18 @@ private:
};
-template< typename LatticeModel_T, typename BoundaryHandling_T, typename Reconstructer_T >
-void PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructer_T >
+template< typename LatticeModel_T, typename BoundaryHandling_T, typename Reconstructer_T, bool conserveMomentum, bool includeGhostLayers >
+void PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructer_T, conserveMomentum, includeGhostLayers >
::operator()( IBlock * const block )
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
+ PdfField_T * pdfField = block->getData< PdfField_T >( pdfFieldID_ );
BoundaryHandling_T * boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingID_ );
FlagField_T * flagField = boundaryHandling->getFlagField();
BodyField_T * bodyField = block->getData< BodyField_T >( bodyFieldID_ );
+ WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandling );
WALBERLA_ASSERT_NOT_NULLPTR( flagField );
WALBERLA_ASSERT_NOT_NULLPTR( bodyField );
@@ -159,7 +175,7 @@ void PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructer_T >
// reconstruct all missing PDFs (only inside the domain, ghost layer values get communicated)
if( optimizeForSmallObstacleFraction_ )
{
- const uint_t numberOfGhostLayersToInclude = uint_t(0);
+ const uint_t numberOfGhostLayersToInclude = includeGhostLayers ? flagField->nrOfGhostLayers() : uint_t(0);
for( auto bodyIt = pe::BodyIterator::begin(*block, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
@@ -167,7 +183,7 @@ void PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructer_T >
continue;
CellInterval cellBB = getCellBB( bodyIt.getBodyID(), *block, *blockStorage_, numberOfGhostLayersToInclude );
- reconstructPDFsInCells(cellBB, block, flagField, formerObstacle );
+ reconstructPDFsInCells(cellBB, block, pdfField, flagField, bodyField, formerObstacle );
}
for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
{
@@ -175,13 +191,13 @@ void PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructer_T >
continue;
CellInterval cellBB = getCellBB( bodyIt.getBodyID(), *block, *blockStorage_, numberOfGhostLayersToInclude );
- reconstructPDFsInCells(cellBB, block, flagField, formerObstacle );
+ reconstructPDFsInCells(cellBB, block, pdfField, flagField, bodyField, formerObstacle );
}
}
else
{
- CellInterval cells = flagField->xyzSize();
- reconstructPDFsInCells(cells, block, flagField, formerObstacle );
+ CellInterval cells = includeGhostLayers ? flagField->xyzSizeWithGhostLayer() : flagField->xyzSize();
+ reconstructPDFsInCells(cells, block, pdfField, flagField, bodyField, formerObstacle );
}
// update the flags from formerObstacle to fluid (inside domain & in ghost layers)
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h b/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h
index fdb0804d1..ad2b8adfa 100644
--- a/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h
+++ b/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h
@@ -43,7 +43,7 @@ namespace pe_coupling {
* \brief Classes to be used together with the PDFReconstruction class to reconstruct PDFs
*
* Each reconstructor must exactly implement the member function
-* void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block );
+* void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField );
* that reconstructs all PDFs in a specific cell with indices x,y,z on a given block.
*
* Different variants are available:
@@ -77,48 +77,46 @@ public:
typedef Field< pe::BodyID, 1 > BodyField_T;
EquilibriumReconstructor( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID,
- const BlockDataID & pdfFieldID, const BlockDataID & bodyFieldID )
- : blockStorage_( blockStorage ), boundaryHandlingID_( boundaryHandlingID ), pdfFieldID_( pdfFieldID ), bodyFieldID_( bodyFieldID )
+ const BlockDataID & bodyFieldID )
+ : blockStorage_( blockStorage ), boundaryHandlingID_( boundaryHandlingID ), bodyFieldID_( bodyFieldID )
{}
- void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block );
+ void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField );
private:
- void setLocalEquilibrium( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block );
+ void setLocalEquilibrium( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField );
- real_t getLocalAverageDensity( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block ) const;
+ real_t getLocalAverageDensity( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField ) const;
shared_ptr<StructuredBlockStorage> blockStorage_;
const BlockDataID boundaryHandlingID_;
- const BlockDataID pdfFieldID_;
const BlockDataID bodyFieldID_;
};
template< typename LatticeModel_T, typename BoundaryHandling_T >
void EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T >
-::operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block )
+::operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField )
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
- setLocalEquilibrium( x, y, z, block);
+ setLocalEquilibrium( x, y, z, block, pdfField);
}
template< typename LatticeModel_T, typename BoundaryHandling_T >
void EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T >
-::setLocalEquilibrium( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block )
+::setLocalEquilibrium( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField )
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
- PdfField_T * pdfField = block->getData< PdfField_T >( pdfFieldID_ );
BodyField_T * bodyField = block->getData< BodyField_T >( bodyFieldID_ );
WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
WALBERLA_ASSERT_NOT_NULLPTR( bodyField );
WALBERLA_ASSERT_EQUAL( pdfField->xyzSize(), bodyField->xyzSize() );
- const real_t averageDensity = getLocalAverageDensity( x, y, z, block );
+ const real_t averageDensity = getLocalAverageDensity( x, y, z, block, pdfField );
real_t cx, cy, cz;
blockStorage_->getBlockLocalCellCenter( *block, Cell(x,y,z), cx, cy, cz );
@@ -131,11 +129,10 @@ void EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T >
template< typename LatticeModel_T, typename BoundaryHandling_T >
real_t EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T >
-::getLocalAverageDensity( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block ) const
+::getLocalAverageDensity( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField ) const
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
- PdfField_T * pdfField = block->getData< PdfField_T >( pdfFieldID_ );
BoundaryHandling_T * boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingID_ );
WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
@@ -175,25 +172,24 @@ public:
typedef Field< pe::BodyID, 1 > BodyField_T;
EquilibriumAndNonEquilibriumReconstructor( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID,
- const BlockDataID & pdfFieldID, const BlockDataID & bodyFieldID,
+ const BlockDataID & bodyFieldID,
const ExtrapolationDirectionFinder_T & extrapolationDirectionFinder )
: blockStorage_( blockStorage ), boundaryHandlingID_( boundaryHandlingID ),
- pdfFieldID_( pdfFieldID ), bodyFieldID_( bodyFieldID ),
+ bodyFieldID_( bodyFieldID ),
extrapolationDirectionFinder_( extrapolationDirectionFinder ),
- equilibriumReconstructor_( EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T >( blockStorage, boundaryHandlingID, pdfFieldID, bodyFieldID ) )
+ equilibriumReconstructor_( EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T >( blockStorage, boundaryHandlingID, bodyFieldID ) )
{ }
- void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block );
- void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, const Vector3<cell_idx_t> & extrapolationDirection );
+ void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField );
+ void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField, const Vector3<cell_idx_t> & extrapolationDirection );
private:
- void setNonEquilibrium( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, const Vector3<cell_idx_t> & extrapolationDirection );
+ void setNonEquilibrium( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField, const Vector3<cell_idx_t> & extrapolationDirection );
shared_ptr<StructuredBlockStorage> blockStorage_;
const BlockDataID boundaryHandlingID_;
- const BlockDataID pdfFieldID_;
const BlockDataID bodyFieldID_;
ExtrapolationDirectionFinder_T extrapolationDirectionFinder_;
@@ -203,37 +199,36 @@ private:
template< typename LatticeModel_T, typename BoundaryHandling_T, typename ExtrapolationDirectionFinder_T >
void EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T >
-::operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block )
+::operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField )
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
Vector3<cell_idx_t> extrapolationDirection(0);
extrapolationDirectionFinder_.getDirection( x, y, z, block, extrapolationDirection );
- (*this)(x, y, z, block, extrapolationDirection);
+ (*this)(x, y, z, block, pdfField, extrapolationDirection);
}
template< typename LatticeModel_T, typename BoundaryHandling_T, typename ExtrapolationDirectionFinder_T >
void EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T >
-::operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, const Vector3<cell_idx_t> & extrapolationDirection )
+::operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField, const Vector3<cell_idx_t> & extrapolationDirection )
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
- equilibriumReconstructor_( x, y, z, block );
+ equilibriumReconstructor_( x, y, z, block, pdfField );
if( extrapolationDirection != cell_idx_t(0) )
{
- setNonEquilibrium( x, y, z, block, extrapolationDirection);
+ setNonEquilibrium( x, y, z, block, pdfField, extrapolationDirection);
}
}
template< typename LatticeModel_T, typename BoundaryHandling_T, typename ExtrapolationDirectionFinder_T >
void EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T >
-::setNonEquilibrium( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, const Vector3<cell_idx_t> & extrapolationDirection )
+::setNonEquilibrium( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField, const Vector3<cell_idx_t> & extrapolationDirection )
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
- PdfField_T * pdfField = block->getData< PdfField_T >( pdfFieldID_ );
BoundaryHandling_T * boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingID_ );
WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
@@ -268,15 +263,15 @@ public:
typedef Field< pe::BodyID, 1 > BodyField_T;
ExtrapolationReconstructor( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID,
- const BlockDataID & pdfFieldID, const BlockDataID & bodyFieldID,
+ const BlockDataID & bodyFieldID,
const ExtrapolationDirectionFinder_T & extrapolationDirectionFinder,
const bool & enforceNoSlipConstraintAfterExtrapolation = false )
: blockStorage_( blockStorage ), boundaryHandlingID_( boundaryHandlingID ),
- pdfFieldID_( pdfFieldID ), bodyFieldID_( bodyFieldID ),
+ bodyFieldID_( bodyFieldID ),
enforceNoSlipConstraintAfterExtrapolation_( enforceNoSlipConstraintAfterExtrapolation ),
extrapolationDirectionFinder_( extrapolationDirectionFinder ),
alternativeReconstructor_( EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T >
- ( blockStorage, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationDirectionFinder ) )
+ ( blockStorage, boundaryHandlingID, bodyFieldID, extrapolationDirectionFinder ) )
{
if( enforceNoSlipConstraintAfterExtrapolation_ ) {
WALBERLA_CHECK((std::is_same<typename LatticeModel_T::Stencil, stencil::D3Q19>::value),
@@ -284,22 +279,21 @@ public:
}
}
- void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block );
+ void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField );
private:
- uint_t getNumberOfExtrapolationCells( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block,
+ uint_t getNumberOfExtrapolationCells( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField,
const Vector3<cell_idx_t> & extrapolationDirection ) const;
- void extrapolatePDFs( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block,
+ void extrapolatePDFs( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField,
const Vector3<cell_idx_t> & extrapolationDirection, const uint_t & numberOfCellsForExtrapolation);
- void enforceNoSlipConstraint( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block );
+ void enforceNoSlipConstraint( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField );
shared_ptr<StructuredBlockStorage> blockStorage_;
const BlockDataID boundaryHandlingID_;
- const BlockDataID pdfFieldID_;
const BlockDataID bodyFieldID_;
const bool enforceNoSlipConstraintAfterExtrapolation_;
@@ -312,37 +306,36 @@ private:
template< typename LatticeModel_T, typename BoundaryHandling_T, typename ExtrapolationDirectionFinder_T >
void ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T >
-::operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block )
+::operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField )
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
Vector3<cell_idx_t> extrapolationDirection(0);
extrapolationDirectionFinder_.getDirection( x, y, z, block, extrapolationDirection );
- uint_t numberOfCellsForExtrapolation = getNumberOfExtrapolationCells( x, y, z, block, extrapolationDirection );
+ uint_t numberOfCellsForExtrapolation = getNumberOfExtrapolationCells( x, y, z, block, pdfField, extrapolationDirection );
if( numberOfCellsForExtrapolation < uint_t(2) )
{
- alternativeReconstructor_( x, y, z, block, extrapolationDirection );
+ alternativeReconstructor_( x, y, z, block, pdfField, extrapolationDirection );
}
else
{
- extrapolatePDFs( x, y, z, block, extrapolationDirection, numberOfCellsForExtrapolation );
+ extrapolatePDFs( x, y, z, block, pdfField, extrapolationDirection, numberOfCellsForExtrapolation );
if( enforceNoSlipConstraintAfterExtrapolation_ )
{
- enforceNoSlipConstraint( x, y, z, block );
+ enforceNoSlipConstraint( x, y, z, block, pdfField );
}
}
}
template< typename LatticeModel_T, typename BoundaryHandling_T, typename ExtrapolationDirectionFinder_T >
uint_t ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T >
-::getNumberOfExtrapolationCells( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block,
+::getNumberOfExtrapolationCells( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField,
const Vector3<cell_idx_t> & extrapolationDirection ) const
{
WALBERLA_ASSERT_NOT_NULLPTR( block );
- PdfField_T * pdfField = block->getData< PdfField_T >( pdfFieldID_ );
BoundaryHandling_T * boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingID_ );
WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
@@ -371,12 +364,9 @@ uint_t ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, Extrapola
template< typename LatticeModel_T, typename BoundaryHandling_T, typename ExtrapolationDirectionFinder_T >
void ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T >
-::extrapolatePDFs( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block,
+::extrapolatePDFs( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const /*block*/, PdfField_T * const pdfField,
const Vector3<cell_idx_t> & extrapolationDirection, const uint_t & numberOfCellsForExtrapolation)
{
- WALBERLA_ASSERT_NOT_NULLPTR( block );
-
- PdfField_T * pdfField = block->getData< PdfField_T >( pdfFieldID_ );
WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
if( numberOfCellsForExtrapolation == uint_t(3) )
@@ -400,13 +390,12 @@ void ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, Extrapolati
template< typename LatticeModel_T, typename BoundaryHandling_T, typename ExtrapolationDirectionFinder_T >
void ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T >
-::enforceNoSlipConstraint( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block )
+::enforceNoSlipConstraint( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, PdfField_T * const pdfField )
{
//NOTE: this currently works only for D3Q19 stencils!
WALBERLA_ASSERT_NOT_NULLPTR( block );
- PdfField_T * pdfField = block->getData< PdfField_T >( pdfFieldID_ );
BodyField_T * bodyField = block->getData< BodyField_T >( bodyFieldID_ );
WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
diff --git a/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp b/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
index cfacb7cd3..5b07af01c 100644
--- a/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
@@ -350,14 +350,14 @@ int main( int argc, char **argv )
"pe Time Step", finestLevel );
// add sweep for updating the pe body mapping into the LBM simulation
- refinementTimestep->addPostStreamVoidFunction( lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID,
+ refinementTimestep->addPostStreamVoidFunction( lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID,
bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag ), blocks ),
"Body Mapping", finestLevel );
// add 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 );
- refinementTimestep->addPostStreamVoidFunction( lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks,
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID );
+ refinementTimestep->addPostStreamVoidFunction( lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks, pdfFieldID,
boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), blocks ),
"PDF Restore", finestLevel );
diff --git a/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp b/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp
index 0edda5942..3a12ec962 100644
--- a/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp
@@ -477,8 +477,8 @@ int main( int argc, char **argv )
// testing utilities
MappingChecker mappingChecker(blocks, bodyStorageID, globalBodyStorage, boundaryHandlingID, bodyFieldID, radius);
MappingResetter mappingResetter(blocks, boundaryHandlingID, bodyFieldID);
- pe_coupling::BodyMapping<BoundaryHandling_T> regularBodyMapper(blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies );
- pe_coupling::BodyMapping<BoundaryHandling_T> globalBodyMapper(blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectGlobalBodies );
+ pe_coupling::BodyMapping<LatticeModel_T, BoundaryHandling_T> regularBodyMapper(blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies );
+ pe_coupling::BodyMapping<LatticeModel_T, BoundaryHandling_T> globalBodyMapper(blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectGlobalBodies );
// sphere positions for test scenarios
diff --git a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
index 502cab175..f3b768189 100644
--- a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
@@ -931,13 +931,13 @@ int main( int argc, char **argv )
// sweep for updating the pe body mapping into the LBM simulation
timeloop.add()
- << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
+ << Sweep( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "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 );
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID );
timeloop.add()
- << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,
+ << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,
reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
// setup of the LBM communication for synchronizing the pdf field between neighboring blocks
std::function< void () > commFunction;
diff --git a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
index 380a10f79..41ebf9d9d 100644
--- a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
@@ -511,14 +511,14 @@ int main( int argc, char **argv )
// sweep for updating the pe body mapping into the LBM simulation
timeloop.add()
- << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ),
+ << Sweep( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ),
"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 );
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID );
timeloop.add()
- << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,
+ << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,
reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
// setup of the LBM communication for synchronizing the pdf field between neighboring blocks
diff --git a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
index 73a9d47cf..0042c361a 100644
--- a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
@@ -636,13 +636,13 @@ int main( int argc, char **argv )
if( MO_CLI )
{
timeloop.add()
- << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_CLI_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
+ << Sweep( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_CLI_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
}else if ( MO_MR ){
timeloop.add()
- << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_MR_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
+ << Sweep( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_MR_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
}else{
timeloop.add()
- << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_BB_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
+ << Sweep( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_BB_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
}
// sweep for restoring PDFs in cells previously occupied by pe bodies
@@ -650,24 +650,24 @@ int main( int argc, char **argv )
{
pe_coupling::FlagFieldNormalExtrapolationDirectionFinder< BoundaryHandling_T > extrapolationFinder( blocks, boundaryHandlingID );
typedef pe_coupling::ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, pe_coupling::FlagFieldNormalExtrapolationDirectionFinder<BoundaryHandling_T> > Reconstructor_T;
- Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder, true );
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID, extrapolationFinder, true );
timeloop.add()
<< Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >
- ( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
+ ( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
} else if ( eanReconstructor )
{
pe_coupling::SphereNormalExtrapolationDirectionFinder extrapolationFinder( blocks, bodyFieldID );
typedef pe_coupling::EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, pe_coupling::SphereNormalExtrapolationDirectionFinder > Reconstructor_T;
- Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder );
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID, extrapolationFinder );
timeloop.add()
<< Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >
- ( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
+ ( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
} else {
typedef pe_coupling::EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T > Reconstructor_T;
- Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID );
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID );
timeloop.add()
<< Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >
- ( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
+ ( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
}
shared_ptr<pe_coupling::BodiesForceTorqueContainer> bodiesFTContainer1 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
index b8c0b6a0b..f8a4e3c73 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
@@ -162,10 +162,12 @@ class SpherePropertyLogger
{
public:
SpherePropertyLogger( SweepTimeloop* timeloop, const shared_ptr< StructuredBlockStorage > & blocks,
+ const BlockDataID & pdfFieldID, const BlockDataID & boundaryHandlingID,
const BlockDataID & bodyStorageID, const std::string & fileName, bool fileIO,
- real_t dx_SI, real_t dt_SI, real_t diameter) :
- timeloop_( timeloop ), blocks_( blocks ), bodyStorageID_( bodyStorageID ), fileName_( fileName ), fileIO_(fileIO),
- dx_SI_( dx_SI ), dt_SI_( dt_SI ), diameter_( diameter ),
+ real_t dx_SI, real_t dt_SI, real_t diameter, real_t mass) :
+ timeloop_( timeloop ), blocks_( blocks ), pdfFieldID_( pdfFieldID ),
+ boundaryHandlingID_( boundaryHandlingID ), bodyStorageID_( bodyStorageID ), fileName_( fileName ), fileIO_(fileIO),
+ dx_SI_( dx_SI ), dt_SI_( dt_SI ), diameter_( diameter ), mass_( mass ),
position_( real_t(0) ), maxVelocity_( real_t(0) )
{
if ( fileIO_ )
@@ -186,14 +188,25 @@ public:
Vector3<real_t> pos(real_t(0));
Vector3<real_t> transVel(real_t(0));
+ Vector3<real_t> force(real_t(0));
+ Vector3<real_t> fluidMomentum(real_t(0));
for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
{
+ PdfField_T * pdfField = blockIt->getData<PdfField_T>(pdfFieldID_);
+ BoundaryHandling_T * bh = blockIt->getData<BoundaryHandling_T>(boundaryHandlingID_);
+ WALBERLA_FOR_ALL_CELLS_XYZ(pdfField,
+ if( bh->isDomain(Cell(x,y,z)) ) fluidMomentum += pdfField->getMomentumDensity(x,y,z);)
+
for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
{
pos = bodyIt->getPosition();
transVel = bodyIt->getLinearVel();
}
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ force += bodyIt->getForce();
+ }
}
WALBERLA_MPI_SECTION()
@@ -205,13 +218,21 @@ public:
mpi::allReduceInplace( transVel[0], mpi::SUM );
mpi::allReduceInplace( transVel[1], mpi::SUM );
mpi::allReduceInplace( transVel[2], mpi::SUM );
+
+ mpi::allReduceInplace( force[0], mpi::SUM );
+ mpi::allReduceInplace( force[1], mpi::SUM );
+ mpi::allReduceInplace( force[2], mpi::SUM );
+
+ mpi::allReduceInplace( fluidMomentum[0], mpi::SUM );
+ mpi::allReduceInplace( fluidMomentum[1], mpi::SUM );
+ mpi::allReduceInplace( fluidMomentum[2], mpi::SUM );
}
position_ = pos[2];
maxVelocity_ = std::max(maxVelocity_, -transVel[2]);
if( fileIO_ )
- writeToFile( timestep, pos, transVel);
+ writeToFile( timestep, pos, transVel, force, fluidMomentum);
}
real_t getPosition() const
@@ -225,7 +246,7 @@ public:
}
private:
- void writeToFile( uint_t timestep, const Vector3<real_t> & position, const Vector3<real_t> & velocity )
+ void writeToFile( uint_t timestep, const Vector3<real_t> & position, const Vector3<real_t> & velocity, const Vector3<real_t> & force, const Vector3<real_t> & fluidMomentum )
{
WALBERLA_ROOT_SECTION()
{
@@ -239,6 +260,9 @@ private:
file << timestep << "\t" << real_c(timestep) * dt_SI_ << "\t"
<< "\t" << scaledPosition[0] << "\t" << scaledPosition[1] << "\t" << scaledPosition[2] - real_t(0.5)
<< "\t" << velocity_SI[0] << "\t" << velocity_SI[1] << "\t" << velocity_SI[2]
+ << "\t" << force[0] << "\t" << force[1] << "\t" << force[2]
+ << "\t" << fluidMomentum[0] << "\t" << fluidMomentum[1] << "\t" << fluidMomentum[2]
+ << "\t" << mass_ * velocity[0] << "\t" << mass_ * velocity[1] << "\t" << mass_ * velocity[2]
<< "\n";
file.close();
}
@@ -246,10 +270,12 @@ private:
SweepTimeloop* timeloop_;
shared_ptr< StructuredBlockStorage > blocks_;
+ const BlockDataID pdfFieldID_;
+ const BlockDataID boundaryHandlingID_;
const BlockDataID bodyStorageID_;
std::string fileName_;
bool fileIO_;
- real_t dx_SI_, dt_SI_, diameter_;
+ real_t dx_SI_, dt_SI_, diameter_, mass_;
real_t position_;
real_t maxVelocity_;
@@ -428,7 +454,7 @@ int main( int argc, char **argv )
auto blocks = blockforest::createUniformBlockGrid( numberOfBlocksPerDirection[0], numberOfBlocksPerDirection[1], numberOfBlocksPerDirection[2],
cellsPerBlockPerDirection[0], cellsPerBlockPerDirection[1], cellsPerBlockPerDirection[2], dx,
0, false, false,
- false, false, false, //periodicity
+ true, true, true, //periodicity
false );
WALBERLA_LOG_INFO_ON_ROOT("Domain decomposition:");
@@ -464,6 +490,7 @@ int main( int argc, char **argv )
// create pe bodies
// bounding planes (global)
+ /*
const auto planeMaterial = pe::createMaterial( "myPlaneMat", real_t(8920), real_t(0), real_t(1), real_t(1), real_t(0), real_t(1), real_t(1), real_t(0), real_t(0) );
pe::createPlane( *globalBodyStorage, 0, Vector3<real_t>(1,0,0), Vector3<real_t>(0,0,0), planeMaterial );
pe::createPlane( *globalBodyStorage, 0, Vector3<real_t>(-1,0,0), Vector3<real_t>(real_c(domainSize[0]),0,0), planeMaterial );
@@ -471,6 +498,7 @@ int main( int argc, char **argv )
pe::createPlane( *globalBodyStorage, 0, Vector3<real_t>(0,-1,0), Vector3<real_t>(0,real_c(domainSize[1]),0), planeMaterial );
pe::createPlane( *globalBodyStorage, 0, Vector3<real_t>(0,0,1), Vector3<real_t>(0,0,0), planeMaterial );
pe::createPlane( *globalBodyStorage, 0, Vector3<real_t>(0,0,-1), Vector3<real_t>(0,0,real_c(domainSize[2])), planeMaterial );
+ */
// add the sphere
const auto sphereMaterial = pe::createMaterial( "mySphereMat", densitySphere , real_t(0.5), real_t(0.1), real_t(0.1), real_t(0.24), real_t(200), real_t(200), real_t(0), real_t(0) );
@@ -500,7 +528,7 @@ int main( int argc, char **argv )
BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >( MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
// map planes into the LBM simulation -> act as no-slip boundaries
- pe_coupling::mapBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectGlobalBodies );
+ //pe_coupling::mapBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectGlobalBodies );
// map pe bodies into the LBM simulation
pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectRegularBodies );
@@ -520,15 +548,15 @@ int main( int argc, char **argv )
// sweep for updating the pe body mapping into the LBM simulation
timeloop.add()
- << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
+ << Sweep( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
// sweep for restoring PDFs in cells previously occupied by pe bodies
pe_coupling::SphereNormalExtrapolationDirectionFinder extrapolationFinder( blocks, bodyFieldID );
typedef pe_coupling::EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, pe_coupling::SphereNormalExtrapolationDirectionFinder > Reconstructor_T;
- Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder );
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID, extrapolationFinder );
timeloop.add()
<< Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >
- ( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
+ ( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
shared_ptr<pe_coupling::BodiesForceTorqueContainer> bodiesFTContainer1 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
std::function<void(void)> storeForceTorqueInCont1 = std::bind(&pe_coupling::BodiesForceTorqueContainer::store, bodiesFTContainer1);
@@ -568,12 +596,6 @@ int main( int argc, char **argv )
}
- Vector3<real_t> gravitationalForce( real_t(0), real_t(0), -(densitySphere - densityFluid) * gravitationalAcceleration * sphereVolume );
- timeloop.addFuncAfterTimeStep( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, gravitationalForce ), "Gravitational force" );
-
- // add pe timesteps
- timeloop.addFuncAfterTimeStep( pe_coupling::TimeStep( blocks, bodyStorageID, cr, syncCall, real_t(1), numPeSubCycles ), "pe Time Step" );
-
// check for convergence of the particle position
std::string loggingFileName( baseFolder + "/LoggingSettlingSphere_");
loggingFileName += std::to_string(fluidType);
@@ -582,10 +604,16 @@ int main( int argc, char **argv )
{
WALBERLA_LOG_INFO_ON_ROOT(" - writing logging output to file \"" << loggingFileName << "\"");
}
- shared_ptr< SpherePropertyLogger > logger = walberla::make_shared< SpherePropertyLogger >( &timeloop, blocks, bodyStorageID,
- loggingFileName, fileIO, dx_SI, dt_SI, diameter );
+ shared_ptr< SpherePropertyLogger > logger = walberla::make_shared< SpherePropertyLogger >( &timeloop, blocks, pdfFieldID, boundaryHandlingID, bodyStorageID,
+ loggingFileName, fileIO, dx_SI, dt_SI, diameter, sphereVolume * densitySphere );
timeloop.addFuncAfterTimeStep( SharedFunctor< SpherePropertyLogger >( logger ), "Sphere property logger" );
+ Vector3<real_t> gravitationalForce( real_t(0), real_t(0), -(densitySphere - densityFluid) * gravitationalAcceleration * sphereVolume );
+ timeloop.addFuncAfterTimeStep( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, gravitationalForce ), "Gravitational force" );
+
+ // add pe timesteps
+ timeloop.addFuncAfterTimeStep( pe_coupling::TimeStep( blocks, bodyStorageID, cr, syncCall, real_t(1), numPeSubCycles ), "pe Time Step" );
+
if( vtkIOFreq != uint_t(0) )
{
// spheres
@@ -624,7 +652,7 @@ int main( int argc, char **argv )
WcTimingPool timeloopTiming;
- real_t terminationPosition = real_t(0.51) * diameter; // right before sphere touches the bottom wall
+ real_t terminationPosition = diameter; // right before sphere touches the bottom wall
// time loop
for (uint_t i = 0; i < timesteps; ++i )
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
index 9cecbe8b7..19fc95b56 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
@@ -714,14 +714,14 @@ int main( int argc, char **argv )
"pe Time Step", finestLevel );
// add sweep for updating the pe body mapping into the LBM simulation
- refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), blocks ),
+ refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), blocks ),
"Body Mapping", finestLevel );
// add sweep for restoring PDFs in cells previously occupied by pe bodies
pe_coupling::SphereNormalExtrapolationDirectionFinder extrapolationFinder( blocks, bodyFieldID );
typedef pe_coupling::EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, pe_coupling::SphereNormalExtrapolationDirectionFinder > Reconstructor_T;
- Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder );
- refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks,
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID, extrapolationFinder );
+ refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks, pdfFieldID,
boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), blocks ),
"PDF Restore", finestLevel );
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
index 6ba596359..e7bc2538a 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
@@ -651,15 +651,15 @@ int main( int argc, char **argv )
"pe Time Step", finestLevel );
// add sweep for updating the pe body mapping into the LBM simulation
- refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID,
+ refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< LatticeModel_T, BoundaryHandling_T >( blocks, pdfFieldID, boundaryHandlingID, bodyStorageID,
globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), blocks ),
"Body Mapping", finestLevel );
// add sweep for restoring PDFs in cells previously occupied by pe bodies
pe_coupling::SphereNormalExtrapolationDirectionFinder extrapolationFinder( blocks, bodyFieldID );
typedef pe_coupling::EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, pe_coupling::SphereNormalExtrapolationDirectionFinder > Reconstructor_T;
- Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder );
- refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks,
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, bodyFieldID, extrapolationFinder );
+ refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks, pdfFieldID,
boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), blocks ),
"PDF Restore", finestLevel );
diff --git a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
index c13ca3550..1fca0831a 100644
--- a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
@@ -336,15 +336,16 @@ int main(int argc, char **argv) {
// sweep for updating the pe body mapping into the LBM simulation
timeloop.add()
- << Sweep(pe_coupling::BodyMapping<BoundaryHandling_T>(blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,
+ << Sweep(pe_coupling::BodyMapping<LatticeModel_T, BoundaryHandling_T, pe_coupling::NaNDestroyer<LatticeModel_T>, true>(blocks, pdfFieldID, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,
MO_BB_Flag, FormerMO_BB_Flag, pe_coupling::selectRegularBodies),
"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);
+ Reconstructor_T reconstructor(blocks, boundaryHandlingID, bodyFieldID);
timeloop.add()
- << Sweep(pe_coupling::PDFReconstruction<LatticeModel_T, BoundaryHandling_T, Reconstructor_T>(blocks,
+ << Sweep(pe_coupling::PDFReconstruction<LatticeModel_T, BoundaryHandling_T, Reconstructor_T, true>(blocks,
+ pdfFieldID,
boundaryHandlingID,
bodyStorageID,
globalBodyStorage,
--
GitLab