diff --git a/src/domain_decomposition/MapPointToPeriodicDomain.cpp b/src/domain_decomposition/MapPointToPeriodicDomain.cpp index 82f365e0e5d9bd733f7207eb407a7a7cadc57036..1c86e7364e0ca2622d9a540b1cb496f3aed42c3e 100644 --- a/src/domain_decomposition/MapPointToPeriodicDomain.cpp +++ b/src/domain_decomposition/MapPointToPeriodicDomain.cpp @@ -40,7 +40,7 @@ void mapPointToPeriodicDomain( const std::array< bool, 3 > & periodic, const AAB shift = std::max( shift, real_t(0) ); x = domain.xMax() - shift; if( isIdentical( x, domain.xMax() ) || x < domain.xMin() ) - x = domain.xMin(); + x = std::nextafter(domain.xMax(), domain.xMin()); } else if( x >= domain.xMax() ) { @@ -61,7 +61,7 @@ void mapPointToPeriodicDomain( const std::array< bool, 3 > & periodic, const AAB shift = std::max( shift, real_t(0) ); y = domain.yMax() - shift; if( isIdentical( y, domain.yMax() ) || y < domain.yMin() ) - y = domain.yMin(); + y = std::nextafter(domain.yMax(), domain.yMin()); } else if( y >= domain.yMax() ) { @@ -82,7 +82,7 @@ void mapPointToPeriodicDomain( const std::array< bool, 3 > & periodic, const AAB shift = std::max( shift, real_t(0) ); z = domain.zMax() - shift; if( isIdentical( z, domain.zMax() ) || z < domain.zMin() ) - z = domain.zMin(); + z = std::nextafter(domain.zMax(), domain.zMin()); } else if( z >= domain.zMax() ) { diff --git a/src/pe/BlockFunctions.h b/src/pe/BlockFunctions.h index c1e8e466d6fb075480b289dbdc545b54013a9d25..8987f473882ae9c6b3c9494a81ddd2e42ff6ffa3 100644 --- a/src/pe/BlockFunctions.h +++ b/src/pe/BlockFunctions.h @@ -45,12 +45,20 @@ bool hasNeighborOwner(const BlockT& block, const Owner& owner) * Looks through all neighboring blocks to find the one whose AABB contains \a pt. * Also checks if \a pt is located in the block itself. * Returns -1 if no valid block is found otherwise the process rank of the containing block is returned. + * + * \attention If periodic boundaries are used you have to make sure the point is mapped to the domain before calling this function! */ template <class BlockT> -Owner findContainingProcess(const BlockT& block, math::Vector3<real_t> pt) +Owner findContainingProcess(const BlockT& block, const math::Vector3<real_t> pt) { + WALBERLA_DEBUG_SECTION() + { + auto pt2 = pt; + block.getBlockStorage().mapToPeriodicDomain(pt2); + WALBERLA_ASSERT_EQUAL(pt, pt2); + } + if (block.getAABB().contains(pt)) return Owner(int_c(block.getProcess()), block.getId().getID()); - if (!block.getBlockStorage().getDomain().contains(pt)) block.getBlockStorage().mapToPeriodicDomain(pt); for( uint_t i = uint_t(0); i != block.getNeighborhoodSize(); ++i ) { if (block.getNeighborAABB(i).contains(pt)) return Owner(int_c(block.getNeighborProcess(i)), block.getNeighborId(i).getID()); diff --git a/src/pe/rigidbody/Owner.h b/src/pe/rigidbody/Owner.h index 53e15de8816e96d0e91bee393bb1242dce51c9c6..3295ca272ef4ee69ff581f3c6a061da8ff3e6877 100644 --- a/src/pe/rigidbody/Owner.h +++ b/src/pe/rigidbody/Owner.h @@ -40,7 +40,7 @@ struct Owner int rank_; //< rank of the owner of the shadow copy IBlockID::IDType blockID_; //< block id of the block the shadow copy is located in - Owner(): rank_(-1) {} + Owner(): rank_(-1), blockID_(0) {} Owner(const int rank, const IBlockID::IDType blockID) : rank_(rank), blockID_(blockID) {} }; //************************************************************************************************* diff --git a/src/pe/synchronization/SyncNextNeighbors.h b/src/pe/synchronization/SyncNextNeighbors.h index 63aa9a2081db232e11a8677cca8eec6253c7fd40..83ad74ac0a839491f1d8023696354b61fa7f3cba 100644 --- a/src/pe/synchronization/SyncNextNeighbors.h +++ b/src/pe/synchronization/SyncNextNeighbors.h @@ -41,6 +41,7 @@ #include "blockforest/BlockForest.h" #include "core/mpi/BufferSystem.h" #include "core/timing/TimingTree.h" +#include "domain_decomposition/MapPointToPeriodicDomain.h" namespace walberla { namespace pe { @@ -56,7 +57,16 @@ void generateSynchonizationMessages(mpi::BufferSystem& bs, const Block& block, B WALBERLA_LOG_DETAIL( "Assembling of body synchronization message starts..." ); // position update - for( auto body = localStorage.begin(); body != localStorage.end(); ) { + for( auto body = localStorage.begin(); body != localStorage.end(); ) + { + //correct position to make sure body is always inside the domain! + if (!body->isFixed()) + { + auto pos = body->getPosition(); + block.getBlockStorage().mapToPeriodicDomain(pos); + body->setPosition(pos); + } + bool isInsideDomain = domain.contains( body->getAABB(), -dx ); WALBERLA_ASSERT( !body->isRemote(), "Central body storage contains remote bodies." ); @@ -146,7 +156,6 @@ void generateSynchonizationMessages(mpi::BufferSystem& bs, const Block& block, B { // Body is no longer locally owned (body is about to be migrated). Owner owner( findContainingProcess( block, gpos ) ); - WALBERLA_ASSERT_UNEQUAL( owner.blockID_, me.blockID_, "Position is " << gpos ); WALBERLA_LOG_DETAIL( "Local body " << b->getSystemID() << " is no longer on process " << body->MPITrait.getOwner() << " but on process " << owner ); @@ -164,9 +173,10 @@ void generateSynchonizationMessages(mpi::BufferSystem& bs, const Block& block, B // of which we own a shadow copy in the next position update since (probably) we no longer require the body but // are still part of its registration list. continue; - } - else + } else { + WALBERLA_ASSERT_UNEQUAL( owner.blockID_, me.blockID_, "Position is " << gpos << "\nlocal Block is: " << block.getAABB() ); + // New owner found among neighbors. WALBERLA_ASSERT_UNEQUAL( owner.blockID_, block.getId().getID(), "Migration is restricted to neighboring blocks." ); @@ -189,6 +199,11 @@ void generateSynchonizationMessages(mpi::BufferSystem& bs, const Block& block, B b->setRemote( true ); // Move body to shadow copy storage. + { + auto pos = b->getPosition(); + correctBodyPosition(domain, block.getAABB().center(), pos); + b->setPosition(pos); + } shadowStorage.add( localStorage.release( body ) ); // Note: We cannot determine here whether we require the body since we do not have up to date shadow copies. diff --git a/src/pe/synchronization/SyncShadowOwners.h b/src/pe/synchronization/SyncShadowOwners.h index cc0400d11d50f78597a9be5fd53666c84f386ae0..49752ae1746f0eb5c77490e01a261da07731a3ee 100644 --- a/src/pe/synchronization/SyncShadowOwners.h +++ b/src/pe/synchronization/SyncShadowOwners.h @@ -81,6 +81,14 @@ void updateAndMigrate( BlockForest& forest, BlockDataID storageID, const bool sy { BodyID b (bodyIt.getBodyID()); + //correct position to make sure body is always inside the domain! + if (!b->isFixed()) + { + auto pos = b->getPosition(); + block.getBlockStorage().mapToPeriodicDomain(pos); + b->setPosition(pos); + } + if( !b->isCommunicating() && !syncNonCommunicatingBodies ) { ++bodyIt; continue; @@ -131,6 +139,11 @@ void updateAndMigrate( BlockForest& forest, BlockDataID storageID, const bool sy b->setRemote( true ); // Move body to shadow copy storage. + { + auto pos = b->getPosition(); + correctBodyPosition(forest.getDomain(), block.getAABB().center(), pos); + b->setPosition(pos); + } shadowStorage.add( localStorage.release( bodyIt ) ); b->MPITrait.deregisterShadowOwner( owner ); @@ -215,7 +228,7 @@ void checkAndResolveOverlap( BlockForest& forest, BlockDataID storageID, const r BodyStorage& shadowStorage = (*storage)[1]; const Owner me( int_c(block.getProcess()), block.getId().getID() ); -// const math::AABB& blkAABB = block->getAABB(); + // const math::AABB& blkAABB = block->getAABB(); for( auto bodyIt = localStorage.begin(); bodyIt != localStorage.end(); ++bodyIt) { @@ -238,10 +251,10 @@ void checkAndResolveOverlap( BlockForest& forest, BlockDataID storageID, const r if (b->MPITrait.getOwner() == nbProcess) continue; // dont send to owner!! if (b->MPITrait.getBlockState( nbProcess.blockID_ )) continue; // only send to neighbor which do not know this body -// WALBERLA_LOG_DEVEL("neighobur aabb: " << block.getNeighborAABB(nb)); -// WALBERLA_LOG_DEVEL("isInsideDomain: " << isInsideDomain); -// WALBERLA_LOG_DEVEL("body AABB: " << b->getAABB()); -// WALBERLA_LOG_DEVEL("neighbour AABB: " << block.getNeighborAABB(nb)); + // WALBERLA_LOG_DEVEL("neighobur aabb: " << block.getNeighborAABB(nb)); + // WALBERLA_LOG_DEVEL("isInsideDomain: " << isInsideDomain); + // WALBERLA_LOG_DEVEL("body AABB: " << b->getAABB()); + // WALBERLA_LOG_DEVEL("neighbour AABB: " << block.getNeighborAABB(nb)); if( (isInsideDomain ? block.getNeighborAABB(nb).intersects( b->getAABB(), dx ) : block.getBlockStorage().periodicIntersect(block.getNeighborAABB(nb), b->getAABB(), dx)) ) { diff --git a/tests/domain_decomposition/CMakeLists.txt b/tests/domain_decomposition/CMakeLists.txt index 2c7cae2624be9baca0d50259c941223c233aa4f9..56fd5c358a5d74fc28ee78b910a0a08b22ca4524 100644 --- a/tests/domain_decomposition/CMakeLists.txt +++ b/tests/domain_decomposition/CMakeLists.txt @@ -4,6 +4,9 @@ # ################################################################################################### +waLBerla_compile_test( NAME MapPointToPeriodicDomain FILES MapPointToPeriodicDomain.cpp DEPENDS core ) +waLBerla_execute_test( NAME MapPointToPeriodicDomain ) + waLBerla_compile_test( NAME PeriodicIntersect FILES PeriodicIntersect.cpp DEPENDS core blockforest ) waLBerla_execute_test( NAME PeriodicIntersect ) diff --git a/tests/domain_decomposition/MapPointToPeriodicDomain.cpp b/tests/domain_decomposition/MapPointToPeriodicDomain.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ee60a1144fb6495c955e2765c0c7dfc1444a433a --- /dev/null +++ b/tests/domain_decomposition/MapPointToPeriodicDomain.cpp @@ -0,0 +1,85 @@ +//====================================================================================================================== +// +// 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 MapToPeriodicDomain.cpp +//! \author Sebastian Eibl <sebastian.eibl@fau.de> +// +//====================================================================================================================== + +#include "core/debug/TestSubsystem.h" +#include "core/Environment.h" +#include "core/math/AABB.h" +#include "core/math/Vector3.h" +#include "domain_decomposition/MapPointToPeriodicDomain.h" + +namespace walberla { +namespace testing { + +int main (int argc, char** argv) +{ + using namespace walberla::domain_decomposition; + + debug::enterTestMode(); + + walberla::Environment walberlaEnv( argc, argv ); + WALBERLA_UNUSED(walberlaEnv); + + std::array<bool, 3> isPeriodic{{true, true, true}}; + math::Vector3<real_t> minCorner(real_t(2), real_t(3), real_t(4)); + math::Vector3<real_t> maxCorner(real_t(3), real_t(4), real_t(5)); + math::AABB domain(minCorner, maxCorner); + + //minCorner -> minCorner + math::Vector3<real_t> p = minCorner; + mapPointToPeriodicDomain( isPeriodic, domain, p); + WALBERLA_CHECK_IDENTICAL( p[0], minCorner[0] ); + WALBERLA_CHECK_IDENTICAL( p[1], minCorner[1] ); + WALBERLA_CHECK_IDENTICAL( p[2], minCorner[2] ); + + //maxCorner -> minCorner + p = maxCorner; + mapPointToPeriodicDomain( isPeriodic, domain, p); + WALBERLA_CHECK_IDENTICAL( p[0], minCorner[0] ); + WALBERLA_CHECK_IDENTICAL( p[1], minCorner[1] ); + WALBERLA_CHECK_IDENTICAL( p[2], minCorner[2] ); + + //minCorner - epsilon -> maxCorner - epsilon + p[0] = std::nextafter(domain.xMin(), real_t(0)); + p[1] = std::nextafter(domain.yMin(), real_t(0)); + p[2] = std::nextafter(domain.zMin(), real_t(0)); + mapPointToPeriodicDomain( isPeriodic, domain, p); + WALBERLA_CHECK_IDENTICAL( p[0], std::nextafter(domain.xMax(), domain.xMin()) ); + WALBERLA_CHECK_IDENTICAL( p[1], std::nextafter(domain.yMax(), domain.yMin()) ); + WALBERLA_CHECK_IDENTICAL( p[2], std::nextafter(domain.zMax(), domain.zMin()) ); + + //maxCorner - epsilon -> maxCorner - epsilon + p[0] = std::nextafter(domain.xMax(), domain.xMin()); + p[1] = std::nextafter(domain.yMax(), domain.yMin()); + p[2] = std::nextafter(domain.zMax(), domain.zMin()); + mapPointToPeriodicDomain( isPeriodic, domain, p); + WALBERLA_CHECK_IDENTICAL( p[0], std::nextafter(domain.xMax(), domain.xMin()) ); + WALBERLA_CHECK_IDENTICAL( p[1], std::nextafter(domain.yMax(), domain.yMin()) ); + WALBERLA_CHECK_IDENTICAL( p[2], std::nextafter(domain.zMax(), domain.zMin()) ); + + return EXIT_SUCCESS; +} + +} //namespace testing +} //namespace walberla + +int main( int argc, char** argv ) +{ + return walberla::testing::main(argc, argv); +} diff --git a/tests/pe/CMakeLists.txt b/tests/pe/CMakeLists.txt index 49002bdc955b67fb6b97361a3a7e38f50e592fa4..bde1a805f1ed64cb52993eba626109dee85ce94a 100644 --- a/tests/pe/CMakeLists.txt +++ b/tests/pe/CMakeLists.txt @@ -94,8 +94,9 @@ set_property( TEST PE_SERIALIZEDESERIALIZE04 PROPERTY DEPENDS PE_SERIALIZEDESERI set_property( TEST PE_SERIALIZEDESERIALIZE08 PROPERTY DEPENDS PE_SERIALIZEDESERIALIZE04 ) #serialize runs of tets to avoid i/o conflicts when running ctest with -jN endif() -waLBerla_compile_test( NAME PE_SHADOWCOPY FILES ShadowCopy.cpp DEPENDS core blockforest ) -waLBerla_execute_test( NAME PE_SHADOWCOPY ) +waLBerla_compile_test( NAME PE_SHADOWCOPY FILES ShadowCopy.cpp DEPENDS core blockforest domain_decomposition ) +waLBerla_execute_test( NAME PE_SHADOWCOPY_NN COMMAND $<TARGET_FILE:PE_SHADOWCOPY> ) +waLBerla_execute_test( NAME PE_SHADOWCOPY_SO COMMAND $<TARGET_FILE:PE_SHADOWCOPY> --syncShadowOwners ) waLBerla_compile_test( NAME PE_SIMPLECCD FILES SimpleCCD.cpp DEPENDS core ) waLBerla_execute_test( NAME PE_SIMPLECCD ) @@ -110,10 +111,14 @@ waLBerla_execute_test( NAME PE_SYNCHRONIZATION09 COMMAND $<TARGET_FILE:PE_SYNC waLBerla_execute_test( NAME PE_SYNCHRONIZATION27 COMMAND $<TARGET_FILE:PE_SYNCHRONIZATION> PROCESSES 27) waLBerla_compile_test( NAME PE_SYNCHRONIZATIONDELETE FILES SynchronizationDelete.cpp DEPENDS core ) -waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE01 COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> ) -waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE03 COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> PROCESSES 3 LABELS longrun) -waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE09 COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> PROCESSES 9 LABELS longrun) -waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE27 COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> PROCESSES 27) +waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE01_NN COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> ) +waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE03_NN COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> PROCESSES 3 LABELS longrun) +waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE09_NN COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> PROCESSES 9 LABELS longrun) +waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE27_NN COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> PROCESSES 27) +waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE01_SO COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> --syncShadowOwners ) +waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE03_SO COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> --syncShadowOwners PROCESSES 3 LABELS longrun) +waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE09_SO COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> --syncShadowOwners PROCESSES 9 LABELS longrun) +waLBerla_execute_test( NAME PE_SYNCHRONIZATIONDELETE27_SO COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONDELETE> --syncShadowOwners PROCESSES 27) waLBerla_compile_test( NAME PE_SYNCHRONIZATIONLARGEBODY FILES SynchronizationLargeBody.cpp DEPENDS core ) waLBerla_execute_test( NAME PE_SYNCHRONIZATIONLARGEBODY01 COMMAND $<TARGET_FILE:PE_SYNCHRONIZATIONLARGEBODY> ) diff --git a/tests/pe/ShadowCopy.cpp b/tests/pe/ShadowCopy.cpp index 8c5dba7c4ad95888979f9b6c88f6109fc763d475..4a1670316768363246310de09f0e3c4e463c0441 100644 --- a/tests/pe/ShadowCopy.cpp +++ b/tests/pe/ShadowCopy.cpp @@ -42,6 +42,19 @@ int main( int argc, char** argv ) walberla::debug::enterTestMode(); walberla::MPIManager::instance()->initializeMPI( &argc, &argv ); + bool syncShadowOwners = false; + for( int i = 1; i < argc; ++i ) + { + if( std::strcmp( argv[i], "--syncShadowOwners" ) == 0 ) syncShadowOwners = true; + } + if (syncShadowOwners) + { + WALBERLA_LOG_DEVEL("running with syncShadowOwners"); + } else + { + WALBERLA_LOG_DEVEL("running with syncNextNeighbour"); + } + shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>(); // create blocks @@ -64,7 +77,7 @@ int main( int argc, char** argv ) // logging::Logging::instance()->setFileLogLevel(logging::Logging::DETAIL); // logging::Logging::instance()->includeLoggingToFile("ShadowCopy"); - bool syncShadowOwners = false; + std::function<void(void)> syncCall; if (!syncShadowOwners) { @@ -95,6 +108,23 @@ int main( int argc, char** argv ) WALBERLA_CHECK_FLOAT_EQUAL( sp->getRadius(), real_t(1.2) ); destroyBodyBySID( *globalBodyStorage, forest->getBlockStorage(), storageID, sid ); + WALBERLA_LOG_PROGRESS_ON_ROOT( " *** SPHERE AT BLOCK EDGE *** "); + sp = pe::createSphere( + *globalBodyStorage, + forest->getBlockStorage(), + storageID, + 999999999, + Vec3(0,2,2), + real_c(1.2)); + sid = sp->getSystemID(); + syncCall(); + sp = static_cast<SphereID> (getBody( *globalBodyStorage, forest->getBlockStorage(), storageID, sid, StorageSelect::LOCAL )); + sp->setPosition(real_c(-1)*std::numeric_limits<real_t>::epsilon(),2,2); + syncCall(); + sp = static_cast<SphereID> (getBody( *globalBodyStorage, forest->getBlockStorage(), storageID, sid, StorageSelect::LOCAL )); + WALBERLA_CHECK_NOT_NULLPTR(sp); + destroyBodyBySID( *globalBodyStorage, forest->getBlockStorage(), storageID, sid ); + WALBERLA_LOG_PROGRESS_ON_ROOT( " *** UNION *** "); UnionT* un = createUnion< boost::tuple<Sphere> >( *globalBodyStorage, forest->getBlockStorage(), storageID, 0, Vec3(2,2,2) ); auto sp1 = createSphere(un, 10, Vec3(real_t(4.9),2,2), real_t(1)); diff --git a/tests/pe/SynchronizationDelete.cpp b/tests/pe/SynchronizationDelete.cpp index fb40078fd1e4b121b02684800cf22d7c7537a2c0..f3b04cec96c65b7156a9ca56cd0d555a3673deb4 100644 --- a/tests/pe/SynchronizationDelete.cpp +++ b/tests/pe/SynchronizationDelete.cpp @@ -61,6 +61,19 @@ int main( int argc, char ** argv ) // logging::Logging::instance()->setFileLogLevel( logging::Logging::DETAIL ); // logging::Logging::instance()->includeLoggingToFile("SyncLog"); + bool syncShadowOwners = false; + for( int i = 1; i < argc; ++i ) + { + if( std::strcmp( argv[i], "--syncShadowOwners" ) == 0 ) syncShadowOwners = true; + } + if (syncShadowOwners) + { + WALBERLA_LOG_DEVEL("running with syncShadowOwners"); + } else + { + WALBERLA_LOG_DEVEL("running with syncNextNeighbour"); + } + shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage> (); // create blocks @@ -109,8 +122,17 @@ int main( int argc, char ** argv ) } } + std::function<void(void)> syncCall; + if (!syncShadowOwners) + { + syncCall = std::bind( pe::syncNextNeighbors<BodyTuple>, std::ref(forest->getBlockForest()), storageID, static_cast<WcTimingTree*>(nullptr), real_c(0.0), false ); + } else + { + syncCall = std::bind( pe::syncShadowOwners<BodyTuple>, std::ref(forest->getBlockForest()), storageID, static_cast<WcTimingTree*>(nullptr), real_c(0.0), false ); + } + for (int i = 0; i < 500; ++i){ - syncNextNeighbors<BodyTuple>(forest->getBlockForest(), storageID); + syncCall(); integrate(*forest, storageID, real_c(0.1)); }