From c917016e8566931d0d70a2fce23c5bf4a0323f20 Mon Sep 17 00:00:00 2001
From: Lukas Werner <lks.werner@fau.de>
Date: Fri, 16 Oct 2020 09:06:05 +0200
Subject: [PATCH] Adaptive grid refinement for LBM and MESA-PD
---
.../refinement/CurlBasedLevelDetermination.h | 178 +++++++++++++++++
src/lbm/refinement/all.h | 1 +
src/lbm_mesapd_coupling/amr/BlockInfo.h | 87 ++++++++
src/lbm_mesapd_coupling/amr/InfoCollection.h | 185 ++++++++++++++++++
.../ParticlePresenceLevelDetermination.h | 101 ++++++++++
.../utility/SubCyclingManager.cpp | 93 +++++++++
.../utility/SubCyclingManager.h | 83 ++++++++
7 files changed, 728 insertions(+)
create mode 100644 src/lbm/refinement/CurlBasedLevelDetermination.h
create mode 100644 src/lbm_mesapd_coupling/amr/BlockInfo.h
create mode 100644 src/lbm_mesapd_coupling/amr/InfoCollection.h
create mode 100644 src/lbm_mesapd_coupling/amr/level_determination/ParticlePresenceLevelDetermination.h
create mode 100644 src/lbm_mesapd_coupling/utility/SubCyclingManager.cpp
create mode 100644 src/lbm_mesapd_coupling/utility/SubCyclingManager.h
diff --git a/src/lbm/refinement/CurlBasedLevelDetermination.h b/src/lbm/refinement/CurlBasedLevelDetermination.h
new file mode 100644
index 000000000..d384449fb
--- /dev/null
+++ b/src/lbm/refinement/CurlBasedLevelDetermination.h
@@ -0,0 +1,178 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm/refinement/all.h b/src/lbm/refinement/all.h
index a63a1ba38..e55e8edb1 100644
--- a/src/lbm/refinement/all.h
+++ b/src/lbm/refinement/all.h
@@ -30,3 +30,4 @@
#include "TimeStepPdfPackInfo.h"
#include "TimeTracker.h"
#include "VorticityBasedLevelDetermination.h"
+#include "CurlBasedLevelDetermination.h"
\ No newline at end of file
diff --git a/src/lbm_mesapd_coupling/amr/BlockInfo.h b/src/lbm_mesapd_coupling/amr/BlockInfo.h
new file mode 100644
index 000000000..6b3e0706a
--- /dev/null
+++ b/src/lbm_mesapd_coupling/amr/BlockInfo.h
@@ -0,0 +1,87 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm_mesapd_coupling/amr/InfoCollection.h b/src/lbm_mesapd_coupling/amr/InfoCollection.h
new file mode 100644
index 000000000..13f8682d5
--- /dev/null
+++ b/src/lbm_mesapd_coupling/amr/InfoCollection.h
@@ -0,0 +1,185 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm_mesapd_coupling/amr/level_determination/ParticlePresenceLevelDetermination.h b/src/lbm_mesapd_coupling/amr/level_determination/ParticlePresenceLevelDetermination.h
new file mode 100644
index 000000000..9ad15259e
--- /dev/null
+++ b/src/lbm_mesapd_coupling/amr/level_determination/ParticlePresenceLevelDetermination.h
@@ -0,0 +1,101 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/lbm_mesapd_coupling/utility/SubCyclingManager.cpp b/src/lbm_mesapd_coupling/utility/SubCyclingManager.cpp
new file mode 100644
index 000000000..ca5c531d3
--- /dev/null
+++ b/src/lbm_mesapd_coupling/utility/SubCyclingManager.cpp
@@ -0,0 +1,93 @@
+//======================================================================================================================
+//
+// 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)
+ : numberOfSubCycles_(numberOfSubCycles), timingPool_(timingPool), currentTimeStep_(0) {}
+
+void SubCyclingManager::operator()() {
+ executeBeforeFunctions();
+ for (size_t s = 0; s < numberOfSubCycles_; ++s) {
+ executeDuringFunctions();
+ }
+ executeAfterFunctions();
+ ++currentTimeStep_;
+}
+
+SubCyclingManager::FuncHandle
+SubCyclingManager::addFuncBeforeSubCycles(const VoidVoidFunc &f, const std::string &identifier) {
+ IdentifiedFunc idFunc(identifier, f);
+ beforeFunctions_.emplace_back(idFunc);
+ return beforeFunctions_.size() - 1;
+}
+
+SubCyclingManager::FuncHandle
+SubCyclingManager::addFuncDuringSubCycles(const VoidVoidFunc &f, const std::string &identifier) {
+ IdentifiedFunc idFunc(identifier, f);
+ duringFunctions_.emplace_back(idFunc);
+ return duringFunctions_.size() - 1;
+}
+
+SubCyclingManager::FuncHandle
+SubCyclingManager::addFuncAfterSubCycles(const VoidVoidFunc &f, const std::string &identifier) {
+ IdentifiedFunc idFunc(identifier, f);
+ afterFunctions_.emplace_back(idFunc);
+ return afterFunctions_.size() - 1;
+}
+
+inline void SubCyclingManager::executeBeforeFunctions() {
+ executeFunctions(beforeFunctions_);
+}
+
+inline void SubCyclingManager::executeDuringFunctions() {
+ executeFunctions(duringFunctions_);
+}
+
+inline void SubCyclingManager::executeAfterFunctions() {
+ executeFunctions(afterFunctions_);
+}
+
+inline void SubCyclingManager::executeFunctions(std::vector<IdentifiedFunc>& functions) {
+ for (const auto &f: functions) {
+ startTiming(f.first);
+ f.second();
+ stopTiming(f.first);
+ }
+}
+
+inline void SubCyclingManager::startTiming(const std::string & name){
+ if (timingPool_ != nullptr) {
+ (*timingPool_)[name].start();
+ }
+}
+
+inline void SubCyclingManager::stopTiming(const std::string & name){
+ if (timingPool_ != nullptr) {
+ (*timingPool_)[name].end();
+ }
+}
+
+}
+}
\ No newline at end of file
diff --git a/src/lbm_mesapd_coupling/utility/SubCyclingManager.h b/src/lbm_mesapd_coupling/utility/SubCyclingManager.h
new file mode 100644
index 000000000..cd748446b
--- /dev/null
+++ b/src/lbm_mesapd_coupling/utility/SubCyclingManager.h
@@ -0,0 +1,83 @@
+//======================================================================================================================
+//
+// 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.h
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#include "core/DataTypes.h"
+#include "core/logging/Logging.h"
+#include "core/timing/TimingPool.h"
+
+namespace walberla {
+namespace lbm_mesapd_coupling {
+
+// loosely based on Timeloop class
+
+class SubCyclingManager {
+public:
+ typedef std::function<void ()> VoidVoidFunc;
+
+ /*! \name Construction & Destruction */
+ //@{
+ explicit SubCyclingManager(size_t numberOfSubCycles, shared_ptr<WcTimingPool> timingPool = nullptr);
+
+ virtual ~SubCyclingManager() {}
+ //@}
+
+ /*! \name Registration Functions */
+ //@{
+ typedef size_t FuncHandle;
+
+ FuncHandle addFuncBeforeSubCycles(const VoidVoidFunc &f, const std::string &identifier = "Other");
+ FuncHandle addFuncDuringSubCycles(const VoidVoidFunc &f, const std::string &identifier = "Other");
+ FuncHandle addFuncAfterSubCycles(const VoidVoidFunc &f, const std::string &identifier = "Other");
+ //@}
+
+ /*! \name Execution Functions */
+ //@{
+ void operator()();
+ //@}
+
+ /*! \name Bookkeeping Functions */
+ //@{
+ uint_t getCurrentTimeStep() const { return currentTimeStep_; }
+ void setCurrentTimeStep(uint_t timestep) { currentTimeStep_ = timestep; }
+ //@}
+
+private:
+ typedef std::pair<std::string, VoidVoidFunc> IdentifiedFunc;
+
+ inline void executeBeforeFunctions();
+ inline void executeDuringFunctions();
+ inline void executeAfterFunctions();
+
+ inline void executeFunctions(std::vector<IdentifiedFunc>& functions);
+
+ inline void startTiming(const std::string & name);
+ inline void stopTiming(const std::string & name);
+
+ size_t numberOfSubCycles_;
+ shared_ptr<WcTimingPool> timingPool_;
+ uint_t currentTimeStep_;
+
+ std::vector<IdentifiedFunc> beforeFunctions_;
+ std::vector<IdentifiedFunc> duringFunctions_;
+ std::vector<IdentifiedFunc> afterFunctions_;
+};
+
+}
+}
\ No newline at end of file
--
GitLab