Commit c917016e authored by Lukas Werner's avatar Lukas Werner
Browse files

Adaptive grid refinement for LBM and MESA-PD

parent 25008165
//======================================================================================================================
//
// 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>
//! \author Lukas Werner <lks.werner@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 local curl
*
* 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.
*
* Parametes upperLimit corresponds to sigma_c, coarsenFactor to c, lowerLimit to c*sigma_c, lengthScaleWeight to r.
*/
template< typename Filter_T >
class CurlBasedLevelDetermination // used as a 'BlockForest::RefreshMinTargetLevelDeterminationFunction'
{
public:
typedef GhostLayerField< Vector3<real_t>, 1 > VectorField_T;
CurlBasedLevelDetermination(const ConstBlockDataID & fieldID, const StructuredBlockForest & structuredBlockForest,
const Filter_T & filter, const uint_t maxLevel,
const real_t upperLimit, const real_t lowerLimit, const real_t lengthScaleWeight = real_t(2)) :
fieldID_(fieldID), structuredBlockForest_(structuredBlockForest), filter_(filter), maxLevel_(maxLevel),
upperLimitSqr_(upperLimit*upperLimit), lowerLimitSqr_(lowerLimit*lowerLimit),
lengthScaleWeight_(lengthScaleWeight)
{
WALBERLA_CHECK_FLOAT_UNEQUAL(lengthScaleWeight_, real_t(0)); // else std::pow(x, y/0) is calculated further below
}
void operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
std::vector< const Block * > & blocksAlreadyMarkedForRefinement,
const BlockForest & forest );
private:
ConstBlockDataID fieldID_;
const StructuredBlockForest & structuredBlockForest_;
Filter_T filter_;
uint_t maxLevel_;
real_t upperLimitSqr_;
real_t lowerLimitSqr_;
real_t lengthScaleWeight_;
};
template< typename Filter_T >
void CurlBasedLevelDetermination< Filter_T >::operator()( std::vector< std::pair< const Block *,
uint_t > > & minTargetLevels, std::vector< const Block * > &, const BlockForest & forest) {
for(auto & minTargetLevel : minTargetLevels) {
const Block * const block = minTargetLevel.first;
const VectorField_T * u = block->template getData< VectorField_T >(fieldID_);
if(u == nullptr) {
minTargetLevel.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), cell_idx_c(-1));
interval.expand(expand);
const auto one = cell_idx_t(1);
const real_t dx = structuredBlockForest_.dx(forest.getLevel(*block));
const real_t dy = structuredBlockForest_.dy(forest.getLevel(*block));
const real_t dz = structuredBlockForest_.dz(forest.getLevel(*block));
const auto halfInvDx = real_t(0.5) * real_t(1) / dx;
const auto halfInvDy = real_t(0.5) * real_t(1) / dy;
const auto halfInvDz = real_t(0.5) * real_t(1) / dz;
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());
const real_t lengthScale = std::cbrt(dx*dy*dz);
const real_t weightedLengthScale = std::pow(lengthScale, (lengthScaleWeight_+1)/lengthScaleWeight_);
const real_t weightedLengthScaleSqr = weightedLengthScale*weightedLengthScale;
for (auto z = cell_idx_t(0); z < zSize; ++z) {
for (auto y = cell_idx_t(0); y < ySize; ++y) {
for (auto 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)
&& 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 = u->get(x,y,z+one);
const Vector3< real_t > zb = u->get(x,y,z-one);
const real_t duzdy = halfInvDy * (ya[2] - yb[2]);
const real_t duydz = halfInvDz * (za[1] - zb[1]);
const real_t duxdz = halfInvDz * (za[0] - zb[0]);
const real_t duzdx = halfInvDx * (xa[2] - xb[2]);
const real_t duydx = halfInvDx * (xa[1] - xb[1]);
const real_t duxdy = halfInvDy * (ya[0] - yb[0]);
const Vector3< real_t > curl( duzdy - duydz, duxdz - duzdx, duydx - duxdy );
const auto curlSqr = curl.sqrLength();
const auto curlSensorSqr = curlSqr * weightedLengthScaleSqr;
if (curlSensorSqr > lowerLimitSqr_) {
// curl is not small enough to coarsen, i.e. stay at least on the current level
coarsen = false;
if (curlSensorSqr > upperLimitSqr_) {
// curl is big enough for refinement
refine = true;
}
}
}
}
}
}
if (refine && block->getLevel() < maxLevel_) {
WALBERLA_ASSERT(!coarsen);
minTargetLevel.second = block->getLevel() + uint_t(1);
}
if (coarsen && block->getLevel() > uint_t(0)) {
WALBERLA_ASSERT(!refine);
minTargetLevel.second = block->getLevel() - uint_t(1);
}
}
}
} // namespace refinement
} // namespace lbm
} // namespace walberla
......@@ -30,3 +30,4 @@
#include "TimeStepPdfPackInfo.h"
#include "TimeTracker.h"
#include "VorticityBasedLevelDetermination.h"
#include "CurlBasedLevelDetermination.h"
\ No newline at end of file
//======================================================================================================================
//
// 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>
//! \author Lukas Werner <lks.werner@fau.de>
//! //
//======================================================================================================================
#pragma once
#include <core/mpi/RecvBuffer.h>
#include <core/mpi/SendBuffer.h>
#include <ostream>
namespace walberla {
namespace lbm_mesapd_coupling {
struct BlockInfo {
// lbm quantities
uint_t numberOfCells;
uint_t numberOfFluidCells;
uint_t numberOfNearBoundaryCells;
// pe quantities
uint_t numberOfLocalParticles;
uint_t numberOfShadowParticles;
uint_t numberOfContacts;
// coupling quantities
uint_t numberOfMESAPDSubCycles;
BlockInfo()
: numberOfCells(0), numberOfFluidCells(0), numberOfNearBoundaryCells(0),
numberOfLocalParticles(0), numberOfShadowParticles(0), numberOfContacts(0), numberOfMESAPDSubCycles(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)
: numberOfCells(numCells), numberOfFluidCells(numFluidCells), numberOfNearBoundaryCells(numNearBoundaryCells),
numberOfLocalParticles(numLocalBodies), numberOfShadowParticles(numShadowParticles), numberOfContacts(numContacts), numberOfMESAPDSubCycles(numPeSubCycles) {}
};
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;
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.numberOfLocalParticles << info.numberOfShadowParticles << info.numberOfContacts
<< info.numberOfMESAPDSubCycles;
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.numberOfLocalParticles >> info.numberOfShadowParticles >> info.numberOfContacts
>> info.numberOfMESAPDSubCycles;
return buf;
}
} // 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 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 <mesa_pd/data/ParticleStorage.h>
#include <map>
namespace walberla {
namespace lbm_mesapd_coupling {
typedef std::map<blockforest::BlockID, BlockInfo> InfoCollection;
typedef std::pair<blockforest::BlockID, BlockInfo> InfoCollectionPair;
template <typename BoundaryHandling_T, typename ParticleAccessor_T>
void createWithNeighborhood(BlockForest& bf, const BlockDataID boundaryHandlingID,
const ParticleAccessor_T& accessor, const uint_t numberOfMESAPDSubCycles,
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 MESAPD quantities
// count block local and (possible) shadow particles here
uint_t numLocalParticles = 0, numShadowParticles = 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));
numLocalParticles++;
} else if (block->getAABB().sqDistance(accessor.getPosition(idx)) < accessor.getInteractionRadius(idx)*accessor.getInteractionRadius(idx)) {
numShadowParticles++;
}
}
// count contacts here
const uint_t numContacts = 0;
BlockInfo blockInfo(numCells, numFluidCells, numNearBoundaryCells, numLocalParticles, numShadowParticles, numContacts, numberOfMESAPDSubCycles);
InfoCollectionPair infoCollectionEntry(block->getId(), blockInfo);
ic.insert( infoCollectionEntry );
for( auto nb = uint_t(0); nb < block->getNeighborhoodSize(); ++nb )
{
bs.sendBuffer( block->getNeighborProcess(nb) ) << infoCollectionEntry;
}
//note: is it necessary to add child blocks already into the info collection?
// here, we still have full geometrical information and can probably determine number of fluid and near boundary cells more easily
// however, the interesting (and most costly) blocks are never refined and thus their child infos is never needed
// see pe/amr/InfoCollection.cpp for an example
}
// size of buffer is unknown and changes with each send
bs.setReceiverInfoFromSendBufferState(false, true);
bs.sendAll();
for( auto recvIt = bs.begin(); recvIt != bs.end(); ++recvIt )
{
while( !recvIt.buffer().isEmpty() )
{
InfoCollectionPair val;
recvIt.buffer() >> val;
ic.insert(val);
}
}
}
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.numberOfLocalParticles, 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);
uint_t 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.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
blockInfo = combinedInfo;
}
}
} // 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 ParticlePresenceLevelDetermination.h
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "blockforest/BlockForest.h"
#include "lbm_mesapd_coupling/amr/InfoCollection.h"
namespace walberla {
namespace lbm_mesapd_coupling {
/*
* Class to determine the minimum level a block can be.
* For coupled LBM-PE simulations the following rules apply:
* - a moving particle will always remain on the finest block
* - a moving particle is not allowed to extend into an area with a coarser block
* - if no moving particle is present, the level can be as coarse as possible (restricted by the 2:1 rule)
* Therefore, if a particle, local or remote (due to particles that are larger than a block), is present on any of the
* neighboring blocks of a certain block, this block's target level is the finest level.
* This, together with a refinement checking frequency that depends on the maximum translational particle velocity,
* ensures the above given requirements.
*/
class ParticlePresenceLevelDetermination
{
public:
ParticlePresenceLevelDetermination( const shared_ptr<lbm_mesapd_coupling::InfoCollection> & infoCollection, uint_t finestLevel) :
infoCollection_( infoCollection ), finestLevel_( finestLevel)
{}
void operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
std::vector< const Block * > &, const BlockForest & /*forest*/ ) {
for (auto &minTargetLevel : minTargetLevels) {
uint_t currentLevelOfBlock = minTargetLevel.first->getLevel();
const uint_t numberOfParticlesInDirectNeighborhood = getNumberOfLocalAndShadowParticlesInNeighborhood(minTargetLevel.first);
uint_t targetLevelOfBlock = currentLevelOfBlock; //keep everything as it is
if ( numberOfParticlesInDirectNeighborhood > uint_t(0) )
{
// set block to finest level if there are particles nearby
targetLevelOfBlock = finestLevel_;
}
else
{
// block could coarsen since there are no particles nearby
if( currentLevelOfBlock > uint_t(0) )
targetLevelOfBlock = currentLevelOfBlock - uint_t(1);
}
WALBERLA_CHECK_LESS_EQUAL(std::abs(int_c(targetLevelOfBlock) - int_c(currentLevelOfBlock)), uint_t(1), "Only level difference of maximum 1 allowed!");
minTargetLevel.second = targetLevelOfBlock;
}
}
private:
uint_t getNumberOfLocalAndShadowParticlesInNeighborhood(const Block * block) {
auto numParticles = uint_t(0);
// add particles of current block
const auto infoIt = infoCollection_->find(block->getId());
WALBERLA_CHECK_UNEQUAL(infoIt, infoCollection_->end(), "Block with ID " << block->getId() << " not found in info collection!");
numParticles += infoIt->second.numberOfLocalParticles + infoIt->second.numberOfShadowParticles;
// add particles of all neighboring blocks
for(uint_t i = 0; i < block->getNeighborhoodSize(); ++i)
{
const BlockID &neighborBlockID = block->getNeighborId(i);
const auto infoItNeighbor = infoCollection_->find(neighborBlockID);
WALBERLA_CHECK_UNEQUAL(infoItNeighbor, infoCollection_->end(), "Neighbor block with ID " << neighborBlockID << " not found in info collection!");
numParticles += infoItNeighbor->second.numberOfLocalParticles + infoItNeighbor->second.numberOfShadowParticles;
}
return numParticles;
}
shared_ptr<lbm_mesapd_coupling::InfoCollection> infoCollection_;
uint_t finestLevel_;
};
} // 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 SubCyclingManager.cpp
//! \author Lukas Werner <lks.werner@fau.de>
//
//======================================================================================================================
#include "SubCyclingManager.h"
namespace walberla {
namespace lbm_mesapd_coupling {
SubCyclingManager::SubCyclingManager(size_t numberOfSubCycles, shared_ptr<WcTimingPool> timingPool)