Commit 281ce66f authored by Christoph Rettinger's avatar Christoph Rettinger
Browse files

Merge branch 'pecoupling_amr_lb_functionality' into 'master'

Functionality for Adaptive Mesh Refinement and Load Balancing with pe coupling

See merge request walberla/walberla!121
parents 0b3c8e30 e6b9bbe7
waLBerla_add_executable( NAME WorkloadEvaluation FILES WorkloadEvaluation.cpp DEPENDS blockforest boundary core field lbm pe pe_coupling postprocessing stencil timeloop vtk )
if( WALBERLA_BUILD_WITH_PARMETIS )
waLBerla_add_executable( NAME AMRSedimentSettling FILES AMRSedimentSettling.cpp DEPENDS blockforest boundary core field lbm pe pe_coupling postprocessing stencil timeloop vtk )
endif()
\ No newline at end of file
add_subdirectory( AdaptiveMeshRefinementFluidParticleCoupling )
add_subdirectory( ComplexGeometry )
add_subdirectory( DEM )
add_subdirectory( MeshDistance )
......
......@@ -84,6 +84,7 @@
#include "lbm/refinement/BoundarySetup.h"
#include "lbm/refinement/PdfFieldSyncPackInfo.h"
#include "lbm/refinement/TimeStep.h"
#include "lbm/refinement/VorticityBasedLevelDetermination.h"
#include "lbm/sweeps/CellwiseSweep.h"
#include "lbm/sweeps/SplitPureSweep.h"
#include "lbm/sweeps/SplitSweep.h"
......@@ -1064,109 +1065,6 @@ Set<SUID> Pseudo2DBlockStateDetermination::operator()( const std::vector< std::p
}
template< typename VectorField_T, typename Filter_T, bool Pseudo2D = false >
class VorticityRefinement // used as a 'BlockForest::RefreshMinTargetLevelDeterminationFunction'
{
public:
VorticityRefinement( const ConstBlockDataID & fieldId, const Filter_T & filter,
const real_t upperLimit, const real_t lowerLimit, const uint_t maxLevel ) :
fieldId_( fieldId ), filter_( filter ),
upperLimit_( upperLimit * upperLimit ), lowerLimit_( lowerLimit * lowerLimit ), maxLevel_( maxLevel )
{}
void operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
std::vector< const Block * > & blocksAlreadyMarkedForRefinement,
const BlockForest & forest );
private:
ConstBlockDataID fieldId_;
Filter_T filter_;
real_t upperLimit_;
real_t lowerLimit_;
uint_t maxLevel_;
}; // class VorticityRefinement
template< typename VectorField_T, typename Filter_T, bool Pseudo2D >
void VorticityRefinement< VectorField_T, Filter_T, Pseudo2D >::operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
std::vector< const Block * > &, const BlockForest & )
{
for( auto it = minTargetLevels.begin(); it != minTargetLevels.end(); ++it )
{
const Block * const block = it->first;
const VectorField_T * u = block->template getData< VectorField_T >( fieldId_ );
if( u == nullptr )
{
it->second = uint_t(0);
continue;
}
CellInterval interval = u->xyzSize();
Cell expand( cell_idx_c(-1), cell_idx_c(-1), Pseudo2D ? cell_idx_t(0) : cell_idx_c(-1) );
interval.expand( expand );
const cell_idx_t one( cell_idx_t(1) );
const real_t half( real_c(0.5) );
bool refine( false );
bool coarsen( true );
filter_( *block );
WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ( interval,
if( filter_(x,y,z) && filter_(x+one,y,z) && filter_(x-one,y,z) && filter_(x,y+one,z) && filter_(x,y-one,z) &&
( Pseudo2D || (filter_(x,y,z+one) && filter_(x,y,z-one)) ) )
{
const Vector3< real_t >& xa = u->get(x+one,y,z);
const Vector3< real_t >& xb = u->get(x-one,y,z);
const Vector3< real_t >& ya = u->get(x,y+one,z);
const Vector3< real_t >& yb = u->get(x,y-one,z);
const Vector3< real_t > za = Pseudo2D ? Vector3< real_t >(0) : u->get(x,y,z+one);
const Vector3< real_t > zb = Pseudo2D ? Vector3< real_t >(0) : u->get(x,y,z-one);
// ATTENTION: dx/y/z is assumed to be equal to '1'!
const real_t duzdy = half * ( ya[2] - yb[2] );
const real_t duydz = half * ( za[1] - zb[1] );
const real_t duxdz = half * ( za[0] - zb[0] );
const real_t duzdx = half * ( xa[2] - xb[2] );
const real_t duydx = half * ( xa[1] - xb[1] );
const real_t duxdy = half * ( ya[0] - yb[0] );
const Vector3< real_t > curl( duzdy - duydz, duxdz - duzdx, duydx - duxdy );
const auto curlSqr = curl.sqrLength();
if( curlSqr > lowerLimit_ )
{
coarsen = false;
if( curlSqr > upperLimit_ )
refine = true;
}
}
)
if( refine && block->getLevel() < maxLevel_ )
{
WALBERLA_ASSERT( !coarsen );
it->second = block->getLevel() + uint_t(1);
}
if( coarsen && block->getLevel() > uint_t(0) )
{
WALBERLA_ASSERT( !refine );
it->second = block->getLevel() - uint_t(1);
}
}
}
// used as a 'BlockForest::RefreshMinTargetLevelDeterminationFunction
void keepInflowOutflowAtTheSameLevel( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
std::vector< const Block * > & blocksAlreadyMarkedForRefinement,
......@@ -2578,8 +2476,8 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
{
const real_t lowerLimit = configBlock.getParameter< real_t >( "curlLowerLimit" );
const real_t upperLimit = configBlock.getParameter< real_t >( "curlUpperLimit" );
VorticityRefinement< VelocityField_T, field::FlagFieldEvaluationFilter<FlagField_T>, Is2D<LatticeModel_T>::value > vorticityRefinement(
lbm::refinement::VorticityBasedLevelDetermination< field::FlagFieldEvaluationFilter<FlagField_T>, Is2D<LatticeModel_T>::value > vorticityRefinement(
velocityFieldId, flagFieldFilter, upperLimit, lowerLimit, configBlock.getParameter< uint_t >( "maxLevel", uint_t(0) ) );
minTargetLevelDeterminationFunctions.add( vorticityRefinement );
......
//======================================================================================================================
//
// 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 VorticityBasedLevelDetermination.h
//! \ingroup lbm
//! \author Florian Schornbaum <florian.schornbaum@fau.de>
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "blockforest/BlockForest.h"
#include "core/math/Vector3.h"
#include "domain_decomposition/BlockDataID.h"
#include "field/GhostLayerField.h"
#include <vector>
namespace walberla {
namespace lbm {
namespace refinement {
/*!\brief Level determination for refinement check based on (scaled) vorticity values
*
* If (scaled) vorticity magnitude is below lowerLimit in all cells of a block, that block could be coarsened.
* If the (scaled) vorticity value is above the upperLimit for at least one cell, that block gets marked for refinement.
* Else, the block remains on the current level.
*
* The scaling originates from neglecting the actual mesh size on the block to obtain different vorticity values for
* different mesh sizes.
*/
template< typename Filter_T, bool Pseudo2D = false >
class VorticityBasedLevelDetermination // used as a 'BlockForest::RefreshMinTargetLevelDeterminationFunction'
{
public:
typedef GhostLayerField< Vector3<real_t>, 1 > VectorField_T;
VorticityBasedLevelDetermination( const ConstBlockDataID & fieldID, const Filter_T & filter,
const real_t upperLimit, const real_t lowerLimit, const uint_t maxLevel ) :
fieldID_( fieldID ), filter_( filter ),
upperLimit_( upperLimit * upperLimit ), lowerLimit_( lowerLimit * lowerLimit ), maxLevel_( maxLevel )
{}
void operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
std::vector< const Block * > & blocksAlreadyMarkedForRefinement,
const BlockForest & forest );
private:
ConstBlockDataID fieldID_;
Filter_T filter_;
real_t upperLimit_;
real_t lowerLimit_;
uint_t maxLevel_;
}; // class VorticityRefinement
template< typename Filter_T, bool Pseudo2D >
void VorticityBasedLevelDetermination< Filter_T, Pseudo2D >::operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
std::vector< const Block * > &, const BlockForest & )
{
for( auto it = minTargetLevels.begin(); it != minTargetLevels.end(); ++it )
{
const Block * const block = it->first;
const VectorField_T * u = block->template getData< VectorField_T >( fieldID_ );
if( u == nullptr )
{
it->second = uint_t(0);
continue;
}
WALBERLA_ASSERT_GREATER_EQUAL(u->nrOfGhostLayers(), uint_t(1));
CellInterval interval = u->xyzSize();
Cell expand( cell_idx_c(-1), cell_idx_c(-1), Pseudo2D ? cell_idx_t(0) : cell_idx_c(-1) );
interval.expand( expand );
const cell_idx_t one( cell_idx_t(1) );
const real_t half( real_c(0.5) );
bool refine( false );
bool coarsen( true );
filter_( *block );
const cell_idx_t xSize = cell_idx_c( interval.xSize() );
const cell_idx_t ySize = cell_idx_c( interval.ySize() );
const cell_idx_t zSize = cell_idx_c( interval.zSize() );
for( cell_idx_t z = cell_idx_t(0); z < zSize; ++z ) {
for (cell_idx_t y = cell_idx_t(0); y < ySize; ++y) {
for (cell_idx_t x = cell_idx_t(0); x < xSize; ++x) {
if( filter_(x,y,z) && filter_(x+one,y,z) && filter_(x-one,y,z) && filter_(x,y+one,z) && filter_(x,y-one,z) &&
( Pseudo2D || (filter_(x,y,z+one) && filter_(x,y,z-one)) ) )
{
const Vector3< real_t > xa = u->get(x+one,y,z);
const Vector3< real_t > xb = u->get(x-one,y,z);
const Vector3< real_t > ya = u->get(x,y+one,z);
const Vector3< real_t > yb = u->get(x,y-one,z);
const Vector3< real_t > za = Pseudo2D ? Vector3< real_t >(0) : u->get(x,y,z+one);
const Vector3< real_t > zb = Pseudo2D ? Vector3< real_t >(0) : u->get(x,y,z-one);
const real_t duzdy = half * ( ya[2] - yb[2] );
const real_t duydz = half * ( za[1] - zb[1] );
const real_t duxdz = half * ( za[0] - zb[0] );
const real_t duzdx = half * ( xa[2] - xb[2] );
const real_t duydx = half * ( xa[1] - xb[1] );
const real_t duxdy = half * ( ya[0] - yb[0] );
// since no dx was used in the gradient calculation, this value is actually curl * dx
// this is needed to have changing curl values on different grid levels
const Vector3< real_t > curl( duzdy - duydz, duxdz - duzdx, duydx - duxdy );
const auto curlSqr = curl.sqrLength();
if( curlSqr > lowerLimit_ )
{
coarsen = false;
if( curlSqr > upperLimit_ )
refine = true;
}
}
}
}
}
if( refine && block->getLevel() < maxLevel_ )
{
WALBERLA_ASSERT( !coarsen );
it->second = block->getLevel() + uint_t(1);
}
if( coarsen && block->getLevel() > uint_t(0) )
{
WALBERLA_ASSERT( !refine );
it->second = block->getLevel() - uint_t(1);
}
}
}
} // namespace refinement
} // namespace lbm
} // namespace walberla
......@@ -29,3 +29,4 @@
#include "TimeStep.h"
#include "TimeStepPdfPackInfo.h"
#include "TimeTracker.h"
#include "VorticityBasedLevelDetermination.h"
......@@ -24,6 +24,7 @@
#pragma once
#include "amr/all.h"
#include "discrete_particle_methods/all.h"
#include "mapping/all.h"
#include "geometry/all.h"
......
//======================================================================================================================
//
// 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 BlockInfo.h
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include <core/mpi/RecvBuffer.h>
#include <core/mpi/SendBuffer.h>
#include <ostream>
namespace walberla {
namespace pe_coupling {
struct BlockInfo {
// lbm quantities
uint_t numberOfCells;
uint_t numberOfFluidCells;
uint_t numberOfNearBoundaryCells;
// pe quantities
uint_t numberOfLocalBodies;
uint_t numberOfShadowBodies;
uint_t numberOfContacts;
// coupling quantities
uint_t numberOfPeSubCycles;
BlockInfo()
: numberOfCells(0), numberOfFluidCells(0), numberOfNearBoundaryCells(0),
numberOfLocalBodies(0), numberOfShadowBodies(0), numberOfContacts(0),
numberOfPeSubCycles(0) {}
BlockInfo(const uint_t numCells, const uint_t numFluidCells, const uint_t numNearBoundaryCells,
const uint_t numLocalBodies, const uint_t numShadowBodies, const uint_t numContacts,
const uint_t numPeSubCycles)
: numberOfCells(numCells), numberOfFluidCells(numFluidCells), numberOfNearBoundaryCells(numNearBoundaryCells),
numberOfLocalBodies(numLocalBodies), numberOfShadowBodies(numShadowBodies), numberOfContacts(numContacts),
numberOfPeSubCycles(numPeSubCycles) {}
};
inline
std::ostream& operator<<( std::ostream& os, const BlockInfo& bi )
{
os << bi.numberOfCells << " / " << bi.numberOfFluidCells << " / " << bi.numberOfNearBoundaryCells << " / "
<< bi.numberOfLocalBodies << " / " << bi.numberOfShadowBodies << " / " << bi.numberOfContacts << " / "
<< bi.numberOfPeSubCycles;
return os;
}
template< typename T, // Element type 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.numberOfLocalBodies << info.numberOfShadowBodies << info.numberOfContacts
<< info.numberOfPeSubCycles;
return buf;
}
template< typename T> // Element type of SendBuffer
mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, BlockInfo& info )
{
buf.readDebugMarker( "pca" );
buf >> info.numberOfCells >> info.numberOfFluidCells >> info.numberOfNearBoundaryCells
>> info.numberOfLocalBodies >> info.numberOfShadowBodies >> info.numberOfContacts
>> info.numberOfPeSubCycles;
return buf;
}
} // namespace pe_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 InfoCollection.h
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#include "InfoCollection.h"
namespace walberla {
namespace pe_coupling {
void getBlockInfoFromInfoCollection( const PhantomBlock * block, const shared_ptr<InfoCollection>& ic,
BlockInfo & blockInfo )
{
WALBERLA_ASSERT_NOT_NULLPTR(block);
if (block->sourceBlockIsLarger())
{
// block is a result of refinement -> BlockInfo object only available for the father block
// there should be no particles on the block (otherwise it would not have been refined)
// and refinement in LBM does not change the number of cells
// we assume that the number of fluid and near boundary cells also stays the same
// (ATTENTION: not true for blocks intersecting with a boundary!)
// -> we can use the information of the father block for weight assignment
auto infoIt = ic->find( block->getId().getFatherId() );
WALBERLA_CHECK_UNEQUAL( infoIt, ic->end(), "Father block with ID " << block->getId().getFatherId() << " not found in info collection!" );
// check the above mentioned assumptions
WALBERLA_ASSERT_EQUAL(infoIt->second.numberOfLocalBodies, uint_t(0));
WALBERLA_ASSERT_EQUAL(infoIt->second.numberOfShadowBodies, uint_t(0));
WALBERLA_ASSERT_EQUAL(infoIt->second.numberOfContacts, uint_t(0));
blockInfo = infoIt->second;
}
else if (block->sourceBlockHasTheSameSize())
{
auto infoIt = ic->find( block->getId() );
WALBERLA_CHECK_UNEQUAL( infoIt, ic->end(), "Block with ID " << block->getId() << " not found in info collection!" );
blockInfo = infoIt->second;
}
else
{
// source block of block is smaller
// block is a result of coarsening -> BlockInfo object is available on all 8 child blocks
// there should be no particles on the block (otherwise it would not have been coarsened)
// and refinement in LBM does not change the number of cells
// we assume that the number of fluid and near boundary cells will be the average of all 8 child blocks
// -> we can use the information of the child blocks for weight assignment
blockforest::BlockID childIdForInit(block->getId(), 0);
auto childForInitIt = ic->find( childIdForInit );
WALBERLA_CHECK_UNEQUAL( childForInitIt, ic->end(), "Child block with ID " << childIdForInit << " not found in info collection!" );
BlockInfo combinedInfo = childForInitIt->second;
uint_t numFluidCells(0), numNearBoundaryCells(0);
for (uint_t child = 0; child < 8; ++child)
{
blockforest::BlockID childId(block->getId(), child);
auto childIt = ic->find( childId );
WALBERLA_CHECK_UNEQUAL( childIt, ic->end(), "Child block with ID " << childId << " not found in info collection!" );
numFluidCells += childIt->second.numberOfFluidCells;
numNearBoundaryCells += childIt->second.numberOfNearBoundaryCells;
// check above mentioned assumptions
WALBERLA_ASSERT_EQUAL(childIt->second.numberOfLocalBodies, uint_t(0));
WALBERLA_ASSERT_EQUAL(childIt->second.numberOfShadowBodies, 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.numberOfLocalBodies = uint_t(0);
combinedInfo.numberOfShadowBodies = uint_t(0);
combinedInfo.numberOfContacts = uint_t(0); //sum
// number of pe sub cycles stays the same
blockInfo = combinedInfo;
}
}
} // namespace pe_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 InfoCollection.h
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "BlockInfo.h"
#include "blockforest/BlockForest.h"
#include "blockforest/BlockID.h"
#include "core/mpi/BufferSystem.h"
#include "pe/ccd/ICCD.h"
#include "pe/fcd/IFCD.h"
#include "pe/rigidbody/BodyStorage.h"
#include <map>
namespace walberla {
namespace pe_coupling {
typedef std::map<blockforest::BlockID, BlockInfo> InfoCollection;
typedef std::pair<blockforest::BlockID, BlockInfo> InfoCollectionPair;
template <typename BoundaryHandling_T>
void createWithNeighborhood(BlockForest& bf, const BlockDataID boundaryHandlingID,
const BlockDataID bodyStorageID, const BlockDataID ccdID, const BlockDataID fcdID,
const uint_t numberOfPeSubCycles,
InfoCollection& ic )
{
ic.clear();
mpi::BufferSystem bs( MPIManager::instance()->comm(), 856 );
for (auto blockIt = bf.begin(); blockIt != bf.end(); ++blockIt)
{
auto * block = static_cast<blockforest::Block*> (&(*blockIt));
// evaluate LBM quantities
BoundaryHandling_T * boundaryHandling = blockIt->getData< BoundaryHandling_T >( boundaryHandlingID );
auto xyzSize = boundaryHandling->getFlagField()->xyzSize();
const uint_t numCells = xyzSize.numCells();
uint_t numFluidCells(0), numNearBoundaryCells(0);
for( auto cellIt = xyzSize.begin(); cellIt != xyzSize.end(); ++cellIt)
{
if( boundaryHandling->isDomain(*cellIt) )
{
++numFluidCells;
}
if( boundaryHandling->isNearBoundary(*cellIt))
{
++numNearBoundaryCells;
}
}
// evaluate PE quantities