Commit 80c258a6 authored by Christoph Rettinger's avatar Christoph Rettinger
Browse files

Fluid-particle coupling load balancing benchmark

parent 6d2c3d31
......@@ -4,6 +4,7 @@ add_subdirectory( DEM )
add_subdirectory( MeshDistance )
add_subdirectory( CouetteFlow )
add_subdirectory( FluidParticleCoupling )
add_subdirectory( FluidParticleCouplingWithLoadBalancing )
add_subdirectory( ForcesOnSphereNearPlaneInShearFlow )
add_subdirectory( GranularGas )
add_subdirectory( IntegratorAccuracy )
......
waLBerla_add_executable( NAME FluidParticleWorkloadEvaluation FILES FluidParticleWorkloadEvaluation.cpp DEPENDS blockforest boundary core field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
waLBerla_add_executable( NAME FluidParticleWorkloadDistribution FILES FluidParticleWorkloadDistribution.cpp DEPENDS blockforest boundary core field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
#pragma once
#include "lbm_mesapd_coupling/amr/BlockInfo.h"
namespace walberla {
namespace lbm_mesapd_coupling {
namespace amr {
/*
* Result from the workload evaluation as described in
* Rettinger, Ruede - "Dynamic Load Balancing Techniques for Particulate Flow Simulations", 2019, Computation
*/
real_t fittedLBMWeightEvaluationFunction(const BlockInfo& blockInfo)
{
uint_t Ce = blockInfo.numberOfCells;
uint_t F = blockInfo.numberOfFluidCells;
real_t weight = real_t(7.597476065046571e-06) * real_c(Ce) + real_t(8.95723566283202e-05) * real_c(F) + real_t(-0.1526111388616016);
return std::max(weight,real_t(0));
}
real_t fittedBHWeightEvaluationFunction(const BlockInfo& blockInfo)
{
uint_t Ce = blockInfo.numberOfCells;
uint_t NB = blockInfo.numberOfNearBoundaryCells;
real_t weight = real_t(1.3067711379655123e-07) * real_c(Ce) + real_t(0.0007289549127205142) * real_c(NB) + real_t(-0.1575698071795788);
return std::max(weight,real_t(0));
}
real_t fittedRPDWeightEvaluationFunction(const BlockInfo& blockInfo)
{
uint_t Pl = blockInfo.numberOfLocalParticles;
uint_t Pg = blockInfo.numberOfGhostParticles;
uint_t Sc = blockInfo.numberOfRPDSubCycles;
real_t cPlPg2 = real_t(2.402288635599054e-05);
real_t cPl = real_t(0.00040932622363097144);
real_t cPg = real_t(0.0007268941363125683);
real_t c = real_t(2.01883028312316e-19);
real_t weight = real_c(Sc) * ( cPlPg2 * real_c(Pl+Pg) * real_c(Pl+Pg) + cPl * real_c(Pl) + cPg * real_c(Pg) + c );
return std::max(weight,real_t(0));
}
real_t fittedCoup1WeightEvaluationFunction(const BlockInfo& blockInfo)
{
uint_t Ce = blockInfo.numberOfCells;
uint_t F = blockInfo.numberOfFluidCells;
uint_t Pl = blockInfo.numberOfLocalParticles;
uint_t Pg = blockInfo.numberOfGhostParticles;
real_t weight = real_t(5.610203409278647e-06) * real_c(Ce) + real_t(-7.257635845636656e-07) * real_c(F) + real_t(0.02049703546054693) * real_c(Pl) + real_t(0.04248208493809902) * real_c(Pg) + real_t(-0.26609470510074784);
return std::max(weight,real_t(0));
}
real_t fittedCoup2WeightEvaluationFunction(const BlockInfo& blockInfo)
{
uint_t Ce = blockInfo.numberOfCells;
uint_t F = blockInfo.numberOfFluidCells;
uint_t Pl = blockInfo.numberOfLocalParticles;
uint_t Pg = blockInfo.numberOfGhostParticles;
real_t weight = real_t(7.198479654682179e-06) * real_c(Ce) + real_t(1.178247475854302e-06) * real_c(F) + real_t(-0.0026401549115124628) * real_c(Pl) + real_t(0.008459646786179298) * real_c(Pg) + real_t(-0.001077320113275954);
return std::max(weight,real_t(0));
}
real_t fittedTotalWeightEvaluationFunction(const BlockInfo& blockInfo)
{
return fittedLBMWeightEvaluationFunction(blockInfo) + fittedBHWeightEvaluationFunction(blockInfo) +
fittedRPDWeightEvaluationFunction(blockInfo) + fittedCoup1WeightEvaluationFunction(blockInfo) +
fittedCoup2WeightEvaluationFunction(blockInfo);
}
} //namespace amr
} //namespace lbm_mesapd_coupling
} //namespace walberla
......@@ -28,28 +28,27 @@
namespace walberla {
namespace lbm_mesapd_coupling {
namespace amr {
struct BlockInfo {
// lbm quantities
uint_t numberOfCells;
uint_t numberOfFluidCells;
uint_t numberOfNearBoundaryCells;
// pe quantities
// rpd quantities
uint_t numberOfLocalParticles;
uint_t numberOfShadowParticles;
uint_t numberOfContacts;
uint_t numberOfGhostParticles;
// coupling quantities
uint_t numberOfMESAPDSubCycles;
uint_t numberOfRPDSubCycles;
BlockInfo()
: numberOfCells(0), numberOfFluidCells(0), numberOfNearBoundaryCells(0),
numberOfLocalParticles(0), numberOfShadowParticles(0), numberOfContacts(0), numberOfMESAPDSubCycles(0) {}
numberOfLocalParticles(0), numberOfGhostParticles(0), numberOfRPDSubCycles(0) {}
BlockInfo(const uint_t numCells, const uint_t numFluidCells, const uint_t numNearBoundaryCells,
const uint_t numLocalBodies, const uint_t numShadowParticles, const uint_t numContacts,
const uint_t numPeSubCycles)
const uint_t numLocalParticles, const uint_t numGhostParticles, const uint_t numRPDSubCycles)
: numberOfCells(numCells), numberOfFluidCells(numFluidCells), numberOfNearBoundaryCells(numNearBoundaryCells),
numberOfLocalParticles(numLocalBodies), numberOfShadowParticles(numShadowParticles), numberOfContacts(numContacts), numberOfMESAPDSubCycles(numPeSubCycles) {}
numberOfLocalParticles(numLocalParticles), numberOfGhostParticles(numGhostParticles), numberOfRPDSubCycles(numRPDSubCycles) {}
};
......@@ -57,19 +56,17 @@ inline
std::ostream& operator<<( std::ostream& os, const BlockInfo& bi )
{
os << bi.numberOfCells << " / " << bi.numberOfFluidCells << " / " << bi.numberOfNearBoundaryCells << " / "
<< bi.numberOfLocalParticles << " / "<< bi.numberOfShadowParticles << " / " << bi.numberOfContacts << " / "
<< bi.numberOfMESAPDSubCycles;
<< bi.numberOfLocalParticles << " / "<< bi.numberOfGhostParticles << " / " << bi.numberOfRPDSubCycles;
return os;
}
template< typename T, // Element type of SendBuffer
typename G> // Growth policy of SendBuffer
typename G> // Growth policy of SendBuffer
mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, const BlockInfo& info )
{
buf.addDebugMarker( "pca" );
buf << info.numberOfCells << info.numberOfFluidCells << info.numberOfNearBoundaryCells
<< info.numberOfLocalParticles << info.numberOfShadowParticles << info.numberOfContacts
<< info.numberOfMESAPDSubCycles;
<< info.numberOfLocalParticles << info.numberOfGhostParticles << info.numberOfRPDSubCycles;
return buf;
}
......@@ -78,10 +75,10 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, BlockInf
{
buf.readDebugMarker( "pca" );
buf >> info.numberOfCells >> info.numberOfFluidCells >> info.numberOfNearBoundaryCells
>> info.numberOfLocalParticles >> info.numberOfShadowParticles >> info.numberOfContacts
>> info.numberOfMESAPDSubCycles;
>> info.numberOfLocalParticles >> info.numberOfGhostParticles >> info.numberOfRPDSubCycles;
return buf;
}
} // namespace amr
} // namespace lbm_mesapd_coupling
} // namespace walberla
......@@ -26,20 +26,21 @@
#include "blockforest/BlockID.h"
#include "core/mpi/BufferSystem.h"
#include <mesa_pd/data/ParticleStorage.h>
#include "mesa_pd/data/Flags.h"
#include <map>
namespace walberla {
namespace lbm_mesapd_coupling {
namespace amr {
typedef std::map<blockforest::BlockID, BlockInfo> InfoCollection;
typedef std::pair<blockforest::BlockID, BlockInfo> InfoCollectionPair;
using InfoCollection = std::map<blockforest::BlockID, BlockInfo> ;
using InfoCollectionPair = std::pair<blockforest::BlockID, BlockInfo>;
template <typename BoundaryHandling_T, typename ParticleAccessor_T>
void createWithNeighborhood(BlockForest& bf, const BlockDataID boundaryHandlingID,
const ParticleAccessor_T& accessor, const uint_t numberOfMESAPDSubCycles,
InfoCollection& ic ) {
void updateAndSyncInfoCollection(BlockForest& bf, const BlockDataID boundaryHandlingID,
const ParticleAccessor_T& accessor, const uint_t numberOfRPDSubCycles,
InfoCollection& ic ) {
ic.clear();
mpi::BufferSystem bs( MPIManager::instance()->comm(), 856 );
......@@ -66,22 +67,21 @@ void createWithNeighborhood(BlockForest& bf, const BlockDataID boundaryHandlingI
}
// evaluate MESAPD quantities
// count block local and (possible) shadow particles here
uint_t numLocalParticles = 0, numShadowParticles = 0;
// count block local and (possible) ghost particles here
uint_t numLocalParticles = 0, numGhostParticles = 0;
for (size_t idx = 0; idx < accessor.size(); ++idx) {
if (block->getAABB().contains(accessor.getPosition(idx))) {
// a particle within a block on the current process cannot be a ghost particle
WALBERLA_ASSERT(!mesa_pd::data::particle_flags::isSet(accessor.getFlags(idx), mesa_pd::data::particle_flags::GHOST));
using namespace walberla::mesa_pd::data::particle_flags;
if( isSet(accessor.getFlags(idx), GLOBAL) ) continue;
if ( block->getAABB().contains(accessor.getPosition(idx)) ) {
numLocalParticles++;
} else if (block->getAABB().sqDistance(accessor.getPosition(idx)) < accessor.getInteractionRadius(idx)*accessor.getInteractionRadius(idx)) {
numShadowParticles++;
numGhostParticles++;
}
}
// count contacts here
const uint_t numContacts = 0;
BlockInfo blockInfo(numCells, numFluidCells, numNearBoundaryCells, numLocalParticles, numShadowParticles, numContacts, numberOfMESAPDSubCycles);
BlockInfo blockInfo(numCells, numFluidCells, numNearBoundaryCells, numLocalParticles, numGhostParticles, numberOfRPDSubCycles);
InfoCollectionPair infoCollectionEntry(block->getId(), blockInfo);
ic.insert( infoCollectionEntry );
......@@ -102,6 +102,8 @@ void createWithNeighborhood(BlockForest& bf, const BlockDataID boundaryHandlingI
bs.setReceiverInfoFromSendBufferState(false, true);
bs.sendAll();
// info collection has to be distirbuted to neighboring processes such that later on when coarsening was applied,
// the weight of the coarsened block can be computed
for( auto recvIt = bs.begin(); recvIt != bs.end(); ++recvIt )
{
while( !recvIt.buffer().isEmpty() )
......@@ -132,7 +134,6 @@ void getBlockInfoFromInfoCollection( const PhantomBlock * block, const shared_pt
// check the above mentioned assumptions
WALBERLA_ASSERT_EQUAL(infoIt->second.numberOfLocalParticles, uint_t(0));
WALBERLA_ASSERT_EQUAL(infoIt->second.numberOfContacts, uint_t(0));
blockInfo = infoIt->second;
}
......@@ -168,18 +169,40 @@ void getBlockInfoFromInfoCollection( const PhantomBlock * block, const shared_pt
// check above mentioned assumptions
WALBERLA_ASSERT_EQUAL(childIt->second.numberOfLocalParticles, uint_t(0));
WALBERLA_ASSERT_EQUAL(childIt->second.numberOfContacts, uint_t(0));
}
// total number of cells remains unchanged
combinedInfo.numberOfFluidCells = uint_c(numFluidCells / uint_t(8)); //average
combinedInfo.numberOfNearBoundaryCells = uint_c( numNearBoundaryCells / uint_t(8) ); //average
combinedInfo.numberOfLocalParticles = uint_t(0);
combinedInfo.numberOfContacts = uint_t(0); //sum
// number of pe sub cycles stays the same
// number of rpd sub cycles stays the same
blockInfo = combinedInfo;
}
}
/*
* Provides mapping from phantom block to weight evaluation via info collection
*/
class WeightEvaluationFunctor
{
public:
WeightEvaluationFunctor(const shared_ptr<InfoCollection>& ic,
const std::function<real_t(const BlockInfo&)> & weightEvaluationFct) :
ic_(ic), weightEvaluationFct_(weightEvaluationFct) {}
real_t operator()(const PhantomBlock * block)
{
BlockInfo blockInfo;
getBlockInfoFromInfoCollection(block, ic_, blockInfo);
return weightEvaluationFct_(blockInfo);
}
private:
shared_ptr<InfoCollection> ic_;
std::function<real_t(const BlockInfo&)> weightEvaluationFct_;
};
} // namepace amr
} // namespace lbm_mesapd_coupling
} // namespace walberla
//======================================================================================================================
//
// 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 MetisAssignmentFunctor.cpp
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#include "MetisAssignmentFunctor.h"
namespace walberla {
namespace lbm_mesapd_coupling {
namespace amr {
void MetisAssignmentFunctor::operator()( std::vector< std::pair< const PhantomBlock *, walberla::any > > & blockData,
const PhantomBlockForest & /*phantomBlockForest*/)
{
for (auto &it : blockData) {
const PhantomBlock * block = it.first;
//only change of one level is supported!
WALBERLA_ASSERT_LESS( std::abs(int_c(block->getLevel()) - int_c(block->getSourceLevel())), 2 );
std::vector<int64_t> metisVertexWeights(weightEvaluationFctVector_.size());
for( auto con = uint_t(0); con < weightEvaluationFctVector_.size(); ++con )
{
real_t vertexWeight = std::max(weightEvaluationFctVector_[con](block), blockBaseWeight_);
int64_t metisVertexWeight = int64_c( weightMultiplicator_ * vertexWeight );
WALBERLA_ASSERT_GREATER(metisVertexWeight, int64_t(0));
metisVertexWeights[con] = metisVertexWeight;
}
blockforest::DynamicParMetisBlockInfo info( metisVertexWeights );
info.setVertexCoords(it.first->getAABB().center() );
real_t blockVolume = it.first->getAABB().volume();
real_t approximateEdgeLength = std::cbrt( blockVolume );
int64_t faceNeighborWeight = int64_c(approximateEdgeLength * approximateEdgeLength ); //common face
int64_t edgeNeighborWeight = int64_c( approximateEdgeLength ); //common edge
int64_t cornerNeighborWeight = int64_c( 1 ); //common corner
int64_t vertexSize = int64_c(blockVolume);
WALBERLA_ASSERT_GREATER(vertexSize, int64_t(0));
info.setVertexSize( vertexSize );
for( const uint_t idx : blockforest::getFaceNeighborhoodSectionIndices() )
{
for( auto nb = uint_t(0); nb < it.first->getNeighborhoodSectionSize(idx); ++nb )
{
auto neighborBlockID = it.first->getNeighborId(idx,nb);
info.setEdgeWeight(neighborBlockID, faceNeighborWeight );
}
}
for( const uint_t idx : blockforest::getEdgeNeighborhoodSectionIndices() )
{
for( auto nb = uint_t(0); nb < it.first->getNeighborhoodSectionSize(idx); ++nb )
{
auto neighborBlockID = it.first->getNeighborId(idx,nb);
info.setEdgeWeight(neighborBlockID, edgeNeighborWeight );
}
}
for( const uint_t idx : blockforest::getCornerNeighborhoodSectionIndices() )
{
for( auto nb = uint_t(0); nb < it.first->getNeighborhoodSectionSize(idx); ++nb )
{
auto neighborBlockID = it.first->getNeighborId(idx,nb);
info.setEdgeWeight(neighborBlockID, cornerNeighborWeight );
}
}
it.second = info;
}
}
} // namespace amr
} // namespace lbm_mesapd_coupling
} // namespace walberla
//======================================================================================================================
//
// 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 MetisAssignmentFunctor.h
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "blockforest/PhantomBlockForest.h"
#include "blockforest/PhantomBlock.h"
#include "blockforest/loadbalancing/DynamicParMetis.h"
#include <functional>
#include <vector>
namespace walberla {
namespace lbm_mesapd_coupling {
namespace amr {
class MetisAssignmentFunctor
{
public:
using WeightEvaluationFct = std::function<real_t(const PhantomBlock *)>;
explicit MetisAssignmentFunctor( const WeightEvaluationFct& weightEvaluationFct)
: weightEvaluationFctVector_(1, weightEvaluationFct) {}
explicit MetisAssignmentFunctor( const std::vector< WeightEvaluationFct > & weightEvaluationFctVector )
: weightEvaluationFctVector_(weightEvaluationFctVector) {}
void operator()( std::vector< std::pair< const PhantomBlock *, walberla::any > > & blockData, const PhantomBlockForest & phantomBlockForest);
uint_t getNumberOfConstrains() const { return weightEvaluationFctVector_.size(); }
inline void setBlockBaseWeight( const real_t blockBaseWeight ){ blockBaseWeight_ = blockBaseWeight; }
inline real_t getBlockBaseWeight() const { return blockBaseWeight_; }
inline void setWeightMultiplicator( const real_t weightMultiplicator ){ weightMultiplicator_ = weightMultiplicator; }
private:
std::vector< WeightEvaluationFct > weightEvaluationFctVector_;
real_t blockBaseWeight_ = real_t(1);
real_t weightMultiplicator_ = real_t(1);
};
} // namespace amr
} // namespace lbm_mesapd_coupling
} // namespace walberla
//======================================================================================================================
//
// 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 WeightAssignmentFunctor.h
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "blockforest/loadbalancing/PODPhantomData.h"
#include "blockforest/PhantomBlockForest.h"
#include "blockforest/PhantomBlock.h"
#include <functional>
namespace walberla {
namespace lbm_mesapd_coupling {
namespace amr {
class WeightAssignmentFunctor
{
public:
using PhantomBlockWeight = walberla::blockforest::PODPhantomWeight<double>;
using WeightEvaluationFct = std::function<real_t(const PhantomBlock *)>;
explicit WeightAssignmentFunctor( const WeightEvaluationFct & weightEvaluationFct ) :
weightEvaluationFct_(weightEvaluationFct) {}
void operator()( std::vector< std::pair< const PhantomBlock *, walberla::any > > & blockData, const PhantomBlockForest & )
{
for (auto &it : blockData) {
const PhantomBlock * block = it.first;
//only change of one level is supported!
WALBERLA_ASSERT_LESS( std::abs(int_c(block->getLevel()) - int_c(block->getSourceLevel())), 2 );
real_t blockWeight = std::max(weightEvaluationFct_(block), blockBaseWeight_);
it.second = PhantomBlockWeight( double_c( blockWeight ) );
}
}
inline void setBlockBaseWeight( const real_t blockBaseWeight ) { blockBaseWeight_ = blockBaseWeight; }
inline real_t getBlockBaseWeight() const { return blockBaseWeight_; }
private:
WeightEvaluationFct weightEvaluationFct_;
real_t blockBaseWeight_ = real_t(1);
};
} // namespace amr
} // namespace lbm_mesapd_coupling
} // namespace walberla
//======================================================================================================================
//
// 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 WeightEvaluationFunctions.h
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "blockforest/PhantomBlock.h"
namespace walberla {
namespace lbm_mesapd_coupling {
namespace amr {
real_t defaultWeightEvaluationFunction(const PhantomBlock * /*block*/)
{
return real_t(1);
}
} // namespace amr
} // namespace lbm_mesapd_coupling
} // namespace walberla
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment