Commit 2d079319 authored by Christoph Rettinger's avatar Christoph Rettinger
Browse files

Merge branch 'pecoupling_utils' into 'master'

Implementation changes pecoupling module

See merge request walberla/walberla!118
parents ddd2c7b7 cd174b62
......@@ -21,13 +21,13 @@
#pragma once
#include "blockforest/StructuredBlockForest.h"
#include "core/math/Vector3.h"
#include "domain_decomposition/StructuredBlockStorage.h"
#include "pe/rigidbody/BodyIterators.h"
#include "pe/synchronization/SyncForces.h"
#include <map>
#include <vector>
#include <array>
namespace walberla {
namespace pe_coupling {
......@@ -36,9 +36,14 @@ class BodiesForceTorqueContainer
{
public:
BodiesForceTorqueContainer( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & bodyStorageID )
: blockStorage_( blockStorage ), bodyStorageID_( bodyStorageID )
{ }
typedef std::map< walberla::id_t, std::array<real_t,6> > ForceTorqueStorage_T;
BodiesForceTorqueContainer( const shared_ptr<StructuredBlockForest> & blockForest, const BlockDataID & bodyStorageID )
: blockForest_( blockForest ), bodyStorageID_( bodyStorageID )
{
// has to be added to the forest (not the storage) to register correctly
bodyForceTorqueStorageID_ = blockForest->addBlockData(make_shared<blockforest::AlwaysCreateBlockDataHandling<ForceTorqueStorage_T> >(), "BodiesForceTorqueContainer");
}
void operator()()
{
......@@ -50,38 +55,18 @@ public:
// clear map
clear();
// sum up all forces/torques from shadow copies on local body (owner)
pe::reduceForces( blockStorage_->getBlockStorage(), bodyStorageID_ );
// send total forces/torques to shadow owners
pe::distributeForces( blockStorage_->getBlockStorage(), bodyStorageID_ );
// (re-)build map
for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
{
auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
auto & f = bodyForceTorqueMap_[ bodyIt->getSystemID() ];
auto & f = (*bodyForceTorqueStorage)[ bodyIt->getSystemID() ];
// only add if body has not been added already before (from another block)
if( f.empty() )
{
const auto & force = bodyIt->getForce();
f.push_back( force[0] );
f.push_back( force[1] );
f.push_back( force[2] );
const auto & torque = bodyIt->getTorque();
f.push_back( torque[0] );
f.push_back( torque[1] );
f.push_back( torque[2] );
}
// reset of force/torque on remote bodies necessary to erase the multiple occurrences of forces/torques on bodies
// (due to call to distributeForces() before)
if ( bodyIt->isRemote() )
{
bodyIt->resetForceAndTorque();
}
const auto & force = bodyIt->getForce();
const auto & torque = bodyIt->getTorque();
f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
}
}
......@@ -89,42 +74,52 @@ public:
void setOnBodies()
{
// owning process sets the force/torque on the bodies
for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
// set the force/torque stored in the block-local map onto all bodies
for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
{
for( auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
const auto &f = bodyForceTorqueMap_[bodyIt->getSystemID()];
WALBERLA_ASSERT( !f.empty(), "When attempting to set force/torque on local body " << bodyIt->getSystemID() << " at position " << bodyIt->getPosition() << ", body was not found in map!");
bodyIt->addForce ( f[0], f[1], f[2] );
bodyIt->addTorque( f[3], f[4], f[5] );
const auto f = bodyForceTorqueStorage->find( bodyIt->getSystemID() );
if( f != bodyForceTorqueStorage->end() )
{
const auto & ftValues = f->second;
bodyIt->addForce ( ftValues[0], ftValues[1], ftValues[2] );
bodyIt->addTorque( ftValues[3], ftValues[4], ftValues[5] );
}
// else: new body has arrived that was not known before
}
}
}
void clear()
{
bodyForceTorqueMap_.clear();
for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
{
auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
bodyForceTorqueStorage->clear();
}
}
void swap( BodiesForceTorqueContainer & other )
{
std::swap( this->bodyForceTorqueMap_, other.bodyForceTorqueMap_ );
std::swap( bodyForceTorqueStorageID_, other.bodyForceTorqueStorageID_);
}
private:
shared_ptr<StructuredBlockStorage> blockStorage_;
shared_ptr<StructuredBlockStorage> blockForest_;
const BlockDataID bodyStorageID_;
std::map< walberla::id_t, std::vector<real_t> > bodyForceTorqueMap_;
BlockDataID bodyForceTorqueStorageID_;
};
class BodyContainerSwapper
{
public:
BodyContainerSwapper( const shared_ptr<BodiesForceTorqueContainer> & cont1,
const shared_ptr<BodiesForceTorqueContainer> & cont2 )
BodyContainerSwapper( const shared_ptr<BodiesForceTorqueContainer> & cont1, const shared_ptr<BodiesForceTorqueContainer> & cont2 )
: cont1_( cont1 ), cont2_( cont2 )
{ }
......
//======================================================================================================================
//
// 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 BodySelectorFunctions.cpp
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#include "BodySelectorFunctions.h"
#include "pe/rigidbody/RigidBody.h"
namespace walberla {
namespace pe_coupling {
bool selectAllBodies(pe::BodyID /*bodyID*/)
{
return true;
}
bool selectRegularBodies(pe::BodyID bodyID)
{
return !bodyID->hasInfiniteMass() && !bodyID->isGlobal();
}
bool selectFixedBodies(pe::BodyID bodyID)
{
return bodyID->hasInfiniteMass() && !bodyID->isGlobal();
}
bool selectGlobalBodies(pe::BodyID bodyID)
{
return bodyID->isGlobal();
}
} // namespace pe_coupling
} // namespace walberla
......@@ -21,29 +21,18 @@
#pragma once
#include "pe/Types.h"
#include "pe/rigidbody/RigidBody.h"
namespace walberla {
namespace pe_coupling {
bool selectAllBodies(pe::BodyID /*bodyID*/)
{
return true;
}
bool selectRegularBodies(pe::BodyID bodyID)
{
return !bodyID->hasInfiniteMass() && !bodyID->isGlobal();
}
bool selectFixedBodies(pe::BodyID bodyID)
{
return bodyID->hasInfiniteMass() && !bodyID->isGlobal();
}
bool selectGlobalBodies(pe::BodyID bodyID)
{
return bodyID->isGlobal();
}
bool selectAllBodies(pe::BodyID /*bodyID*/);
bool selectRegularBodies(pe::BodyID bodyID);
bool selectFixedBodies(pe::BodyID bodyID);
bool selectGlobalBodies(pe::BodyID bodyID);
} // namespace pe_coupling
} // namespace walberla
......@@ -35,7 +35,7 @@
#include <functional>
#include <map>
#include <array>
namespace walberla {
namespace pe_coupling {
......@@ -88,43 +88,25 @@ public:
// Since they are reset after the call to collision response, they have to be stored explicitly before.
// Then they are set again after each intermediate step.
// sum up all forces/torques from shadow copies on local body (owner)
pe::reduceForces( blockStorage_->getBlockStorage(), bodyStorageID_ );
// send total forces/torques to shadow owners
pe::distributeForces( blockStorage_->getBlockStorage(), bodyStorageID_ );
// generate map from all known bodies (process local) to total forces/torques
std::map< walberla::id_t, std::vector< real_t > > forceMap;
// this has to be done on a block-local basis, since the same body could reside on several blocks from this process
typedef domain_decomposition::IBlockID::IDType BlockID_T;
std::map< BlockID_T, std::map< walberla::id_t, std::array< real_t, 6 > > > forceTorqueMap;
for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
{
BlockID_T blockID = blockIt->getId().getID();
auto& blockLocalForceTorqueMap = forceTorqueMap[blockID];
// iterate over local and remote bodies and store force/torque in map
// Remote bodies are required since during the then following collision response time steps,
// bodies can change ownership (i.e. a now remote body could become a local body ).
// Since only the owning process re-sets the force/torque later, remote body values have to be stored as well.
for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
auto & f = forceMap[ bodyIt->getSystemID() ];
auto & f = blockLocalForceTorqueMap[ bodyIt->getSystemID() ];
// only add if body has not been added already before (from another block)
if( f.empty() )
{
const auto & force = bodyIt->getForce();
f.push_back( force[0] );
f.push_back( force[1] );
f.push_back( force[2] );
const auto & torque = bodyIt->getTorque();
f.push_back( torque[0] );
f.push_back( torque[1] );
f.push_back( torque[2] );
}
const auto & force = bodyIt->getForce();
const auto & torque = bodyIt->getTorque();
// reset of force/torque on remote bodies necessary to erase the multiple occurrences of forces/torques on bodies
// (due to call to distributeForces() before)
if ( bodyIt->isRemote() )
{
bodyIt->resetForceAndTorque();
}
f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
}
}
......@@ -133,18 +115,26 @@ public:
for( uint_t i = 0; i != numberOfSubIterations_; ++i )
{
// in the first iteration, forces on local bodies are already set by force synchronization
// in the first iteration, forces are already set
if( i != 0 )
{
for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
{
// only owning process sets force/torque on bodies
for( auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
BlockID_T blockID = blockIt->getId().getID();
auto& blockLocalForceTorqueMap = forceTorqueMap[blockID];
// re-set stored force/torque on bodies
for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
const auto & f = forceMap[ bodyIt->getSystemID() ];
WALBERLA_ASSERT( !f.empty(), "When attempting to set force/torque on local body " << bodyIt->getSystemID() << " at position " << bodyIt->getPosition() << ", body was not found in map!");
bodyIt->addForce ( f[0], f[1], f[2] );
bodyIt->addTorque( f[3], f[4], f[5] );
const auto f = blockLocalForceTorqueMap.find( bodyIt->getSystemID() );
if( f != blockLocalForceTorqueMap.end() )
{
const auto & ftValues = f->second;
bodyIt->addForce ( ftValues[0], ftValues[1], ftValues[2] );
bodyIt->addTorque( ftValues[3], ftValues[4], ftValues[5] );
}
}
}
}
......
......@@ -156,3 +156,15 @@ waLBerla_execute_test( NAME SphereWallCollisionBehaviorDPMFuncTest COMMAND $<TAR
waLBerla_compile_test( FILES discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp DEPENDS blockforest pe timeloop )
waLBerla_execute_test( NAME HinderedSettlingDynamicsDPMFuncTest COMMAND $<TARGET_FILE:HinderedSettlingDynamicsDPM> --funcTest PROCESSES 4 LABELS longrun CONFIGURATIONS RelWithDbgInfo )
###################################################################################################
# Utility tests
###################################################################################################
waLBerla_compile_test( FILES utility/BodiesForceTorqueContainerTest.cpp DEPENDS blockforest pe timeloop )
waLBerla_execute_test( NAME BodiesForceTorqueContainerTest COMMAND $<TARGET_FILE:BodiesForceTorqueContainerTest> PROCESSES 1 )
waLBerla_execute_test( NAME BodiesForceTorqueContainerParallelTest COMMAND $<TARGET_FILE:BodiesForceTorqueContainerTest> PROCESSES 3 )
waLBerla_compile_test( FILES utility/PeSubCyclingTest.cpp DEPENDS blockforest pe timeloop )
waLBerla_execute_test( NAME PeSubCyclingTest COMMAND $<TARGET_FILE:PeSubCyclingTest> PROCESSES 1 )
waLBerla_execute_test( NAME PeSubCyclingParallelTest COMMAND $<TARGET_FILE:PeSubCyclingTest> PROCESSES 3 )
\ 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 BodiesForceTorqueContainerTest.cpp
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#include "blockforest/Initialization.h"
#include "core/DataTypes.h"
#include "core/Environment.h"
#include "core/debug/TestSubsystem.h"
#include "core/math/all.h"
#include "pe/basic.h"
#include "pe/utility/DestroyBody.h"
#include <pe_coupling/utility/all.h>
namespace force_torque_container_test
{
///////////
// USING //
///////////
using namespace walberla;
typedef boost::tuple<pe::Sphere> BodyTypeTuple ;
/*!\brief Test cases for the force torque container provided by the coupling module
*
* Spheres at different positions are created and force(s) and torque(s) are applied.
* Then, they are store in the container and reset on the body.
* Then, in some cases, the sphere is moved to cross block borders, and the stored forces/torques are reapplied by the container again.
* The obtained force/torque values are compared against the originally set (and thus expected) values.
*/
//////////
// MAIN //
//////////
int main( int argc, char **argv )
{
debug::enterTestMode();
mpi::Environment env( argc, argv );
// uncomment to have logging
//logging::Logging::instance()->setLogLevel(logging::Logging::LogLevel::DETAIL);
const real_t dx = real_t(1);
const real_t radius = real_t(5);
///////////////////////////
// DATA STRUCTURES SETUP //
///////////////////////////
Vector3<uint_t> blocksPerDirection(uint_t(3), uint_t(1), uint_t(1));
Vector3<uint_t> cellsPerBlock(uint_t(20), uint_t(20), uint_t(20));
Vector3<bool> periodicity(true, false, false);
// create fully periodic domain with refined blocks
auto blocks = blockforest::createUniformBlockGrid( blocksPerDirection[0], blocksPerDirection[1], blocksPerDirection[2],
cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2],
dx,
0, false, false,
periodicity[0], periodicity[1], periodicity[2],
false );
// pe body storage
pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "Storage");
auto sphereMaterialID = pe::createMaterial( "sphereMat", real_t(1) , real_t(0.3), real_t(0.2), real_t(0.2), real_t(0.24), real_t(200), real_t(200), real_t(0), real_t(0) );
// pe coupling
const real_t overlap = real_t( 1.5 ) * dx;
std::function<void(void)> syncCall = std::bind( pe::syncNextNeighbors<BodyTypeTuple>, boost::ref(blocks->getBlockForest()), bodyStorageID, static_cast<WcTimingTree*>(NULL), overlap, false );
// sphere positions for test scenarios
Vector3<real_t> positionInsideBlock(real_t(10), real_t(10), real_t(10));
Vector3<real_t> positionAtBlockBorder(real_t(19.5), real_t(10), real_t(10));
Vector3<real_t> positionAtBlockBorderUpdated(real_t(20.5), real_t(10), real_t(10));
Vector3<real_t> positionAtBlockBorder2(real_t(20) + radius + overlap - real_t(0.5), real_t(10), real_t(10));
Vector3<real_t> positionAtBlockBorderUpdated2(real_t(20) + radius + overlap + real_t(0.5), real_t(10), real_t(10));
Vector3<real_t> testForce(real_t(2), real_t(1), real_t(0));
Vector3<real_t> torqueOffset = Vector3<real_t>(real_t(1), real_t(0), real_t(0));
pe_coupling::ForceTorqueOnBodiesResetter resetter(blocks, bodyStorageID);
shared_ptr<pe_coupling::BodiesForceTorqueContainer> container1 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
shared_ptr<pe_coupling::BodiesForceTorqueContainer> container2 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
pe_coupling::BodyContainerSwapper containerSwapper(container1, container2);
//////////////////
// Inside block //
//////////////////
{
std::string testIdentifier("Test: sphere inside block");
WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
positionInsideBlock, radius, sphereMaterialID, false, true, false);
syncCall();
uint_t applicationCounter( 0 );
for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
++applicationCounter;
auto pos = bodyIt->getPosition();
bodyIt->addForceAtPos(testForce, pos+torqueOffset);
}
}
mpi::allReduceInplace(applicationCounter, mpi::SUM);
container1->store();
resetter();
containerSwapper();
container2->setOnBodies();
Vector3<real_t> expectedForce = applicationCounter * testForce;
Vector3<real_t> expectedTorque = applicationCounter * ( torqueOffset % testForce );
WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting force: " << expectedForce);
WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting torque: " << expectedTorque);
Vector3<real_t> actingForce(real_t(0));
Vector3<real_t> actingTorque(real_t(0));
for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
actingForce += bodyIt->getForce();
actingTorque += bodyIt->getTorque();
}
}
mpi::allReduceInplace(actingForce[0], mpi::SUM);
mpi::allReduceInplace(actingForce[1], mpi::SUM);
mpi::allReduceInplace(actingForce[2], mpi::SUM);
mpi::allReduceInplace(actingTorque[0], mpi::SUM);
mpi::allReduceInplace(actingTorque[1], mpi::SUM);
mpi::allReduceInplace(actingTorque[2], mpi::SUM);
WALBERLA_ROOT_SECTION()
{
WALBERLA_CHECK_FLOAT_EQUAL(actingForce[0], expectedForce[0], "Mismatch in force0");
WALBERLA_CHECK_FLOAT_EQUAL(actingForce[1], expectedForce[1], "Mismatch in force1");
WALBERLA_CHECK_FLOAT_EQUAL(actingForce[2], expectedForce[2], "Mismatch in force2");
WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[0], expectedTorque[0], "Mismatch in torque0");
WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[1], expectedTorque[1], "Mismatch in torque1");
WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[2], expectedTorque[2], "Mismatch in torque2");
}
// clean up
for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
bodyIt->markForDeletion();
}
}
syncCall();
container1->clear();
container2->clear();
WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
}
/////////////////////
// At block border //
/////////////////////
{
std::string testIdentifier("Test: sphere at block border");
WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
positionAtBlockBorder, radius, sphereMaterialID, false, true, false);
syncCall();
uint_t applicationCounter( 0 );
for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
++applicationCounter;
auto pos = bodyIt->getPosition();
bodyIt->addForceAtPos(testForce, pos+torqueOffset);
}
}
mpi::allReduceInplace(applicationCounter, mpi::SUM);
container1->store();
resetter();
containerSwapper();
container2->setOnBodies();
Vector3<real_t> expectedForce = applicationCounter * testForce;
Vector3<real_t> expectedTorque = applicationCounter * ( torqueOffset % testForce );
WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting force: " << expectedForce);
WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting torque: " << expectedTorque);