From c6d8a70ce0a1ef1302f9801b99fbde8abcc8f3ab Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Wed, 7 Aug 2019 12:13:14 +0200
Subject: [PATCH] [ADD] next neighbors sync optimized for blockforest
---
src/mesa_pd/kernel/AssocToBlock.h | 80 +++++
.../mpi/SyncNextNeighborsBlockForest.cpp | 297 ++++++++++++++++++
.../mpi/SyncNextNeighborsBlockForest.h | 80 +++++
tests/mesa_pd/CMakeLists.txt | 3 +
tests/mesa_pd/kernel/SyncNextNeighbors.cpp | 17 +-
.../kernel/SyncNextNeighborsBlockForest.cpp | 168 ++++++++++
6 files changed, 640 insertions(+), 5 deletions(-)
create mode 100644 src/mesa_pd/kernel/AssocToBlock.h
create mode 100644 src/mesa_pd/mpi/SyncNextNeighborsBlockForest.cpp
create mode 100644 src/mesa_pd/mpi/SyncNextNeighborsBlockForest.h
create mode 100644 tests/mesa_pd/kernel/SyncNextNeighborsBlockForest.cpp
diff --git a/src/mesa_pd/kernel/AssocToBlock.h b/src/mesa_pd/kernel/AssocToBlock.h
new file mode 100644
index 000000000..a795a3051
--- /dev/null
+++ b/src/mesa_pd/kernel/AssocToBlock.h
@@ -0,0 +1,80 @@
+//======================================================================================================================
+//
+// 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 AssocToBlock.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+
+#include <blockforest/BlockForest.h>
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+
+class AssocToBlock
+{
+public:
+ explicit AssocToBlock(const std::shared_ptr<BlockForest>& bf) : bf_(bf) {}
+
+ template <typename Accessor>
+ void operator()(const size_t i, Accessor& ac) const;
+private:
+ std::shared_ptr<BlockForest> bf_ = nullptr;
+};
+
+template <typename Accessor>
+inline void AssocToBlock::operator()(const size_t idx,
+ Accessor& ac) const
+{
+ blockforest::Block*& currentBlock = ac.getCurrentBlockRef(idx);
+
+ if (currentBlock != nullptr)
+ {
+ if (currentBlock->getAABB().contains(ac.getPosition(idx)))
+ {
+ return;
+ } else
+ {
+ currentBlock = nullptr;
+ }
+ }
+
+ for (auto& blk : bf_->getBlockMap())
+ {
+ if (blk.second->getAABB().contains(ac.getPosition(idx)))
+ {
+ currentBlock = blk.second.get();
+ return;
+ }
+ }
+
+ WALBERLA_CHECK_NOT_NULLPTR(currentBlock, ac.getPosition(idx) << "\n" << bf_->begin()->getAABB());
+}
+
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/src/mesa_pd/mpi/SyncNextNeighborsBlockForest.cpp b/src/mesa_pd/mpi/SyncNextNeighborsBlockForest.cpp
new file mode 100644
index 000000000..00f336f81
--- /dev/null
+++ b/src/mesa_pd/mpi/SyncNextNeighborsBlockForest.cpp
@@ -0,0 +1,297 @@
+//======================================================================================================================
+//
+// 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 SyncNextNeighborsBlockForest.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#include "SyncNextNeighborsBlockForest.h"
+#include <mesa_pd/domain/BlockForestDomain.h>
+
+#include <mesa_pd/mpi/RemoveAndNotify.h>
+
+namespace walberla {
+namespace mesa_pd {
+namespace mpi {
+
+void correctParticlePosition(Vec3& pt,
+ const Vec3& center,
+ const math::AABB& domain,
+ const std::array<bool, 3>& periodic)
+{
+ const Vec3 dis = pt - center;
+
+ if (periodic[0] && (-domain.xSize() * 0.5 > dis[0])) pt[0] += domain.xSize();
+ if (periodic[0] && (+domain.xSize() * 0.5 < dis[0])) pt[0] -= domain.xSize();
+
+ if (periodic[1] && (-domain.ySize() * 0.5 > dis[1])) pt[1] += domain.ySize();
+ if (periodic[1] && (+domain.ySize() * 0.5 < dis[1])) pt[1] -= domain.ySize();
+
+ if (periodic[2] && (-domain.zSize() * 0.5 > dis[2])) pt[2] += domain.zSize();
+ if (periodic[2] && (+domain.zSize() * 0.5 < dis[2])) pt[2] -= domain.zSize();
+}
+
+void SyncNextNeighborsBlockForest::operator()(data::ParticleStorage& ps,
+ const std::shared_ptr<blockforest::BlockForest>& bf,
+ const real_t dx) const
+{
+ if (numProcesses_ == 1) return;
+
+
+ for (auto& blk : bf->getBlockMap())
+ {
+ for( uint_t i = uint_t(0); i != blk.second->getNeighborhoodSize(); ++i )
+ {
+ auto nbProcessRank = blk.second->getNeighborProcess(i);
+ if (bs.sendBuffer(nbProcessRank).isEmpty())
+ {
+ // fill empty buffers with a dummy byte to force transmission
+ bs.sendBuffer(nbProcessRank) << walberla::uint8_c(0);
+ }
+ }
+ }
+
+ generateSynchronizationMessages(ps, bf, dx);
+
+ // size of buffer is unknown and changes with each send
+ bs.setReceiverInfoFromSendBufferState(false, true);
+ bs.sendAll();
+
+ // Receiving the updates for the remote rigid bodies from the connected processes
+ WALBERLA_LOG_DETAIL( "Parsing of particle synchronization response starts..." );
+ ParseMessage parseMessage;
+ domain::BlockForestDomain domain(bf);
+ for( auto it = bs.begin(); it != bs.end(); ++it )
+ {
+ walberla::uint8_t tmp;
+ it.buffer() >> tmp;
+ while( !it.buffer().isEmpty() )
+ {
+ parseMessage(it.rank(), it.buffer(), ps, domain);
+ }
+ }
+ WALBERLA_LOG_DETAIL( "Parsing of particle synchronization response ended." );
+}
+
+void SyncNextNeighborsBlockForest::generateSynchronizationMessages(data::ParticleStorage& ps,
+ const std::shared_ptr<blockforest::BlockForest>& bf,
+ const real_t dx) const
+{
+ const uint_t ownRank = uint_c(rank_);
+ std::array<bool, 3> periodic;
+ periodic[0] = bf->isPeriodic(0);
+ periodic[1] = bf->isPeriodic(1);
+ periodic[2] = bf->isPeriodic(2);
+
+ WALBERLA_LOG_DETAIL( "Assembling of particle synchronization message starts..." );
+
+ // position update
+ for( auto pIt = ps.begin(); pIt != ps.end(); )
+ {
+ //skip all ghost particles
+ if (data::particle_flags::isSet( pIt->getFlags(), data::particle_flags::GHOST))
+ {
+ ++pIt;
+ continue;
+ }
+
+ //skip all particles that do not communicate (create ghost particles) on other processes
+ if (data::particle_flags::isSet( pIt->getFlags(), data::particle_flags::NON_COMMUNICATING))
+ {
+ ++pIt;
+ continue;
+ }
+
+ auto& currentBlock = pIt->getCurrentBlockRef();
+ WALBERLA_CHECK_NOT_NULLPTR(currentBlock);
+ if (isInsideAABB(pIt->getPosition(), pIt->getInteractionRadius() + dx, currentBlock->getAABB()))
+ {
+ //no sync needed
+ //just delete ghost particles if there are any
+
+ for (const auto& ghostOwner : pIt->getGhostOwners() )
+ {
+ auto& buffer( bs.sendBuffer(static_cast<walberla::mpi::MPIRank>(ghostOwner)) );
+
+ WALBERLA_LOG_DETAIL( "Sending removal notification for particle " << pIt->getUid() << " to process " << ghostOwner );
+
+ packNotification(buffer, ParticleRemovalNotification( *pIt ));
+ }
+
+ pIt->getGhostOwnersRef().clear();
+
+ ++pIt;
+ continue;
+ }
+
+ //correct position to make sure particle is always inside the domain!
+ //everything is decided by the master particle therefore ghost particles are not touched
+ if (!data::particle_flags::isSet( pIt->getFlags(), data::particle_flags::FIXED) &&
+ !data::particle_flags::isSet( pIt->getFlags(), data::particle_flags::GHOST))
+ {
+ bf->mapToPeriodicDomain( pIt->getPositionRef() );
+ }
+
+ // Note: At this point we know that the particle was locally owned before the position update.
+ WALBERLA_CHECK_EQUAL(pIt->getOwner(), ownRank);
+
+ WALBERLA_LOG_DETAIL( "Processing local particle " << pIt->getUid() );
+
+ // Update nearest neighbor processes.
+ auto isInsideDomain = isInsideAABB(pIt->getPosition(), pIt->getInteractionRadius(), bf->getDomain());
+ std::vector<int> ranksAlreadyTreated{int_c(ownRank)};
+ for( uint_t nb = uint_t(0); nb < currentBlock->getNeighborhoodSize(); ++nb )
+ {
+ auto nbProcessRank = currentBlock->getNeighborProcess(nb);
+ if (std::find(ranksAlreadyTreated.begin(), ranksAlreadyTreated.end(), int_c(nbProcessRank)) != ranksAlreadyTreated.end())
+ {
+ continue;
+ }
+ auto nbAABB = currentBlock->getNeighborAABB(nb);
+ auto sqDistance = isInsideDomain
+ ? sqDistancePointToAABB(pIt->getPosition(), nbAABB)
+ : sqDistancePointToAABBPeriodic(pIt->getPosition(), nbAABB, bf->getDomain(), periodic);
+ auto tmp = pIt->getInteractionRadius() + dx;
+ if( sqDistance < tmp*tmp)
+ {
+ ranksAlreadyTreated.emplace_back(int_c(nbProcessRank));
+ auto ghostOwnerIt = std::find( pIt->getGhostOwners().begin(), pIt->getGhostOwners().end(), nbProcessRank );
+ if( ghostOwnerIt != pIt->getGhostOwners().end() )
+ {
+ // already a ghost there -> update
+ auto& buffer( bs.sendBuffer(nbProcessRank) );
+ WALBERLA_LOG_DETAIL( "Sending update notification for particle " << pIt->getUid() << " to process " << (nbProcessRank) );
+ packNotification(buffer, ParticleUpdateNotification( *pIt ));
+ } else
+ {
+ // no ghost there -> create ghost
+ auto& buffer( bs.sendBuffer(nbProcessRank) );
+ WALBERLA_LOG_DETAIL( "Sending shadow copy notification for particle " << pIt->getUid() << " to process " << (nbProcessRank) );
+ packNotification(buffer, ParticleCopyNotification( *pIt ));
+ pIt->getGhostOwnersRef().insert( int_c(nbProcessRank) );
+ }
+ }
+ }
+ for (auto ghostOwnerIt = pIt->getGhostOwnersRef().begin();
+ ghostOwnerIt != pIt->getGhostOwnersRef().end();
+ )
+ {
+ if (std::find(ranksAlreadyTreated.begin(),
+ ranksAlreadyTreated.end(),
+ int_c(*ghostOwnerIt)) == ranksAlreadyTreated.end())
+ {
+ // In case the rigid particle no longer intersects the remote process nor interacts with it but is registered,
+ // send removal notification.
+ auto& buffer( bs.sendBuffer(*ghostOwnerIt) );
+
+ WALBERLA_LOG_DETAIL( "Sending removal notification for particle " << pIt->getUid() << " to process " << *ghostOwnerIt );
+
+ packNotification(buffer, ParticleRemovalNotification( *pIt ));
+
+ ghostOwnerIt = pIt->getGhostOwnersRef().erase(ghostOwnerIt);
+
+ continue;
+ }
+ ++ghostOwnerIt;
+ }
+
+ //particle has left subdomain?
+ if (currentBlock->getAABB().contains(pIt->getPosition()))
+ {
+ // particle still is locally owned after position update.
+ WALBERLA_LOG_DETAIL( "Owner of particle " << pIt->getUid() << " is still process " << pIt->getOwner() );
+ } else
+ {
+ //find new owner
+ int ownerRank = -1;
+ for( uint_t i = uint_t(0); i != currentBlock->getNeighborhoodSize(); ++i )
+ {
+ if (currentBlock->getNeighborAABB(i).contains(pIt->getPosition()))
+ {
+ ownerRank = int_c(currentBlock->getNeighborProcess(i));
+ }
+ }
+
+ if( ownerRank != int_c(ownRank) )
+ {
+ WALBERLA_LOG_DETAIL( "Local particle " << pIt->getUid() << " is no longer on process " << ownRank << " but on process " << ownerRank );
+
+ if( ownerRank < 0 )
+ {
+ // No owner found: Outflow condition.
+ WALBERLA_LOG_DETAIL( "Sending deletion notifications for particle " << pIt->getUid() << " due to outflow." );
+
+ // Registered processes receive removal notification in the remove() routine.
+ pIt = removeAndNotify( bs, ps, pIt );
+
+ continue;
+ }
+
+ WALBERLA_LOG_DETAIL( "Sending migration notification for particle " << pIt->getUid() << " to process " << ownerRank << "." );
+ //WALBERLA_LOG_DETAIL( "Process registration list before migration: " << pIt->getGhostOwners() );
+
+ // Set new owner and transform to ghost particle
+ pIt->setOwner(ownerRank);
+ data::particle_flags::set( pIt->getFlagsRef(), data::particle_flags::GHOST );
+
+ // currently position is mapped to periodically to global domain,
+ // this might not be the correct position for a ghost particle
+ correctParticlePosition( pIt->getPositionRef(), currentBlock->getAABB().center(), bf->getDomain(), periodic );
+
+ // Correct registration list (exclude new owner and us - the old owner) and
+ // notify registered processes (except for new owner) of (remote) migration since they possess a ghost particle.
+ auto ownerIt = std::find( pIt->getGhostOwners().begin(), pIt->getGhostOwners().end(), ownerRank );
+ WALBERLA_CHECK_UNEQUAL(ownerIt, pIt->getGhostOwners().end(), "New owner has to be former ghost owner!" );
+
+ pIt->getGhostOwnersRef().erase( ownerIt );
+
+ for( auto ghostRank : pIt->getGhostOwners() )
+ {
+ auto& buffer( bs.sendBuffer(static_cast<walberla::mpi::MPIRank>(ghostRank)) );
+
+ WALBERLA_LOG_DETAIL( "Sending remote migration notification for particle " << pIt->getUid() <<
+ " to process " << ghostRank );
+
+ packNotification(buffer, ParticleRemoteMigrationNotification( *pIt, ownerRank ));
+ }
+
+ pIt->getGhostOwnersRef().insert( int_c(ownRank) );
+
+ // Send migration notification to new owner
+ auto& buffer( bs.sendBuffer(ownerRank) );
+ packNotification(buffer, ParticleMigrationNotification( *pIt ));
+
+ pIt->getGhostOwnersRef().clear();
+
+ continue;
+ }
+ }
+
+ ++pIt;
+ }
+
+ WALBERLA_LOG_DETAIL( "Assembling of particle synchronization message ended." );
+}
+
+} // namespace mpi
+} // namespace mesa_pd
+} // namespace walberla
diff --git a/src/mesa_pd/mpi/SyncNextNeighborsBlockForest.h b/src/mesa_pd/mpi/SyncNextNeighborsBlockForest.h
new file mode 100644
index 000000000..0f9b01453
--- /dev/null
+++ b/src/mesa_pd/mpi/SyncNextNeighborsBlockForest.h
@@ -0,0 +1,80 @@
+//======================================================================================================================
+//
+// 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 SyncNextNeighborsBlockForest.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/Flags.h>
+#include <mesa_pd/data/ParticleStorage.h>
+#include <mesa_pd/domain/IDomain.h>
+#include <mesa_pd/mpi/notifications/PackNotification.h>
+#include <mesa_pd/mpi/notifications/ParseMessage.h>
+#include <mesa_pd/mpi/notifications/ParticleCopyNotification.h>
+#include <mesa_pd/mpi/notifications/ParticleMigrationNotification.h>
+#include <mesa_pd/mpi/notifications/ParticleRemoteMigrationNotification.h>
+#include <mesa_pd/mpi/notifications/ParticleRemovalNotification.h>
+#include <mesa_pd/mpi/notifications/ParticleUpdateNotification.h>
+
+#include <blockforest/BlockForest.h>
+
+#include <core/mpi/BufferSystem.h>
+#include <core/logging/Logging.h>
+
+namespace walberla {
+namespace mesa_pd {
+namespace mpi {
+
+/**
+ * Kernel which updates all ghost particles.
+ *
+ * \ingroup mesa_pd_mpi
+ */
+class SyncNextNeighborsBlockForest
+{
+public:
+ void operator()(data::ParticleStorage& ps,
+ const std::shared_ptr<blockforest::BlockForest>& blockforest,
+ const real_t dx = real_t(0)) const;
+
+ int64_t getBytesSent() const { return bs.getBytesSent(); }
+ int64_t getBytesReceived() const { return bs.getBytesReceived(); }
+
+ int64_t getNumberOfSends() const { return bs.getNumberOfSends(); }
+ int64_t getNumberOfReceives() const { return bs.getNumberOfReceives(); }
+private:
+ void generateSynchronizationMessages(data::ParticleStorage& ps,
+ const std::shared_ptr<blockforest::BlockForest>& blockforest,
+ const real_t dx) const;
+
+ mutable walberla::mpi::BufferSystem bs = walberla::mpi::BufferSystem( walberla::mpi::MPIManager::instance()->comm() );
+
+ int numProcesses_ = walberla::mpi::MPIManager::instance()->numProcesses();
+ int rank_ = walberla::mpi::MPIManager::instance()->rank();
+};
+
+} // namespace mpi
+} // namespace mesa_pd
+} // namespace walberla
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index 69278f743..e8945e940 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -137,6 +137,9 @@ waLBerla_execute_test( NAME MESA_PD_Kernel_SyncGhostOwnersLarge PROCESSES 27 )
waLBerla_compile_test( NAME MESA_PD_Kernel_SyncNextNeighbors FILES kernel/SyncNextNeighbors.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_SyncNextNeighbors PROCESSES 27 )
+waLBerla_compile_test( NAME MESA_PD_Kernel_SyncNextNeighborsBlockForest FILES kernel/SyncNextNeighborsBlockForest.cpp DEPENDS core )
+waLBerla_execute_test( NAME MESA_PD_Kernel_SyncNextNeighborsBlockForest PROCESSES 27 )
+
waLBerla_compile_test( NAME MESA_PD_Kernel_TemperatureIntegration FILES kernel/TemperatureIntegration.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_TemperatureIntegration )
diff --git a/tests/mesa_pd/kernel/SyncNextNeighbors.cpp b/tests/mesa_pd/kernel/SyncNextNeighbors.cpp
index 93a7d6768..99869faf7 100644
--- a/tests/mesa_pd/kernel/SyncNextNeighbors.cpp
+++ b/tests/mesa_pd/kernel/SyncNextNeighbors.cpp
@@ -64,14 +64,21 @@ int main( int argc, char ** argv )
walberla::mpi::MPIManager::instance()->useWorldComm();
//logging::Logging::instance()->setStreamLogLevel(logging::Logging::DETAIL);
-// logging::Logging::instance()->includeLoggingToFile("MESA_PD_Kernel_SyncNextNeighbor");
-// logging::Logging::instance()->setFileLogLevel(logging::Logging::DETAIL);
+ logging::Logging::instance()->includeLoggingToFile("MESA_PD_Kernel_SyncNextNeighbor");
+ logging::Logging::instance()->setFileLogLevel(logging::Logging::DETAIL);
//init domain partitioning
auto forest = blockforest::createBlockForest( AABB(-15,-15,-15,15,15,15), // simulation domain
Vector3<uint_t>(3,3,3), // blocks in each direction
- Vector3<bool>(true, true, true) // periodicity
+ Vector3<bool>(true, true, true), // periodicity
+ 27, // number of processes
+ 1 // initial refinement
);
+ math::AABB localDomain = forest->begin()->getAABB();
+ for (auto& iBlk : *forest)
+ {
+ localDomain.merge(iBlk.getAABB());
+ }
domain::BlockForestDomain domain(forest);
std::array< bool, 3 > periodic;
periodic[0] = forest->isPeriodic(0);
@@ -118,10 +125,10 @@ int main( int argc, char ** argv )
SNN(ps, domain);
//check
- if (sqDistancePointToAABBPeriodic(pos, forest->begin()->getAABB(), forest->getDomain(), periodic) <= radius * radius)
+ if (sqDistancePointToAABBPeriodic(pos, localDomain, forest->getDomain(), periodic) <= radius * radius)
{
WALBERLA_CHECK_EQUAL(ps.size(), 1);
- if (forest->begin()->getAABB().contains(pos))
+ if (localDomain.contains(pos))
{
WALBERLA_CHECK(!data::particle_flags::isSet(ps.begin()->getFlags(), data::particle_flags::GHOST));
} else
diff --git a/tests/mesa_pd/kernel/SyncNextNeighborsBlockForest.cpp b/tests/mesa_pd/kernel/SyncNextNeighborsBlockForest.cpp
new file mode 100644
index 000000000..25a811757
--- /dev/null
+++ b/tests/mesa_pd/kernel/SyncNextNeighborsBlockForest.cpp
@@ -0,0 +1,168 @@
+//======================================================================================================================
+//
+// 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 SyncNextNeighborsBlockForest.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/data/ParticleStorage.h>
+#include <mesa_pd/domain/BlockForestDomain.h>
+#include <mesa_pd/kernel/AssocToBlock.h>
+#include <mesa_pd/kernel/ParticleSelector.h>
+#include <mesa_pd/mpi/SyncNextNeighborsBlockForest.h>
+
+#include <blockforest/BlockForest.h>
+#include <blockforest/Initialization.h>
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+#include <core/mpi/Reduce.h>
+
+#include <iostream>
+#include <memory>
+
+namespace walberla {
+namespace mesa_pd {
+
+const real_t radius = real_t(1);
+
+walberla::id_t createSphere(data::ParticleStorage& ps, domain::IDomain& domain)
+{
+ walberla::id_t uid = 0;
+ auto owned = domain.isContainedInProcessSubdomain( uint_c(walberla::mpi::MPIManager::instance()->rank()), Vec3(0,0,0) );
+ if (owned)
+ {
+ data::Particle&& p = *ps.create();
+ p.getPositionRef() = Vec3(0,0,0);
+ p.getInteractionRadiusRef() = radius;
+ p.getRotationRef() = Rot3(Quat());
+ p.getLinearVelocityRef() = Vec3(1,2,3);
+ p.getAngularVelocityRef() = Vec3(4,5,6);
+ p.getOwnerRef() = walberla::mpi::MPIManager::instance()->rank();
+ uid = p.getUid();
+ WALBERLA_LOG_DETAIL("SPHERE CREATED");
+ }
+
+ walberla::mpi::allReduceInplace(uid, walberla::mpi::SUM);
+ return uid;
+}
+
+int main( int argc, char ** argv )
+{
+ Environment env(argc, argv);
+ WALBERLA_UNUSED(env);
+ walberla::mpi::MPIManager::instance()->useWorldComm();
+
+ //logging::Logging::instance()->setStreamLogLevel(logging::Logging::DETAIL);
+ logging::Logging::instance()->includeLoggingToFile("MESA_PD_Kernel_SyncNextNeighbor");
+ logging::Logging::instance()->setFileLogLevel(logging::Logging::DETAIL);
+
+ //init domain partitioning
+ auto forest = blockforest::createBlockForest( AABB(-15,-15,-15,15,15,15), // simulation domain
+ Vector3<uint_t>(3,3,3), // blocks in each direction
+ Vector3<bool>(true, true, true), // periodicity
+ 27, // number of processes
+ 1 // initial refinement
+ );
+ //checking blocks distribution
+ WALBERLA_CHECK_EQUAL(forest->size(), 8);
+
+ math::AABB localDomain = forest->begin()->getAABB();
+ for (auto& iBlk : *forest)
+ {
+ localDomain.merge(iBlk.getAABB());
+ }
+
+ domain::BlockForestDomain domain(forest);
+ std::array< bool, 3 > periodic;
+ periodic[0] = forest->isPeriodic(0);
+ periodic[1] = forest->isPeriodic(1);
+ periodic[2] = forest->isPeriodic(2);
+
+ //init data structures
+ auto ps = std::make_shared<data::ParticleStorage>(100);
+ data::ParticleAccessor ac(ps);
+
+ //initialize particle
+ auto uid = createSphere(*ps, domain);
+ WALBERLA_LOG_DEVEL_ON_ROOT("uid: " << uid);
+
+ //init kernels
+ kernel::AssocToBlock assoc(forest);
+ mpi::SyncNextNeighborsBlockForest SNN;
+
+ std::vector<real_t> deltas { real_t(0.1),
+ real_t(4.9),
+ real_t(5.1),
+ real_t(9.9),
+ real_t(10.1),
+ real_t(14.9),
+ real_t(-14.9),
+ real_t(-10.1),
+ real_t(-9.9),
+ real_t(-5.1),
+ real_t(-4.9),
+ real_t(-0.1)};
+
+ for (auto delta : deltas)
+ {
+ WALBERLA_LOG_DEVEL_ON_ROOT(delta);
+
+ ps->forEachParticle(false, kernel::SelectLocal(), ac, assoc, ac);
+
+ auto pos = Vec3(1,-1,1) * delta;
+ WALBERLA_LOG_DETAIL("checking position: " << pos);
+ // owner moves particle to new position
+ auto pIt = ps->find(uid);
+ if (pIt != ps->end())
+ {
+ if (!data::particle_flags::isSet(pIt->getFlags(), data::particle_flags::GHOST))
+ {
+ pIt->setPosition(pos);
+ }
+ }
+
+ //sync
+ SNN(*ps, forest);
+
+ //check
+ if (sqDistancePointToAABBPeriodic(pos, localDomain, forest->getDomain(), periodic) <= radius * radius)
+ {
+ WALBERLA_CHECK_EQUAL(ps->size(), 1);
+ if (localDomain.contains(pos))
+ {
+ WALBERLA_CHECK(!data::particle_flags::isSet(ps->begin()->getFlags(), data::particle_flags::GHOST));
+ } else
+ {
+ WALBERLA_CHECK(data::particle_flags::isSet(ps->begin()->getFlags(), data::particle_flags::GHOST));
+ }
+ } else
+ {
+ WALBERLA_CHECK_EQUAL(ps->size(), 0);
+ }
+ }
+
+
+ return EXIT_SUCCESS;
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
+int main( int argc, char ** argv )
+{
+ return walberla::mesa_pd::main(argc, argv);
+}
--
GitLab