Skip to content
Snippets Groups Projects
Commit b09aaf8b authored by Christoph Rettinger's avatar Christoph Rettinger
Browse files

Merge branch 'new_pe_coupling_functionalities' into 'master' 

New pe coupling functionalities

See merge request walberla/walberla!56
parents 6b14b8dd d9c623a7
Branches
Tags
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 580 additions and 440 deletions
apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
+ 41
− 153
View file @ b09aaf8b
......@@ -237,13 +237,14 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
return handling;
}
class ForceEval
/*
* Functionality for continuous logging of sphere properties during initial (resting) simulation
*/
class RestingSphereForceEvaluator
{
public:
ForceEval( SweepTimeloop* timeloop,
const shared_ptr< StructuredBlockStorage > & blocks, const ConstBlockDataID & bodyStorageID,
uint_t averageFrequency, const std::string & basefolder ) :
RestingSphereForceEvaluator( SweepTimeloop* timeloop, const shared_ptr< StructuredBlockStorage > & blocks,
const ConstBlockDataID & bodyStorageID, uint_t averageFrequency, const std::string & basefolder ) :
timeloop_( timeloop ), blocks_( blocks ), bodyStorageID_( bodyStorageID ), averageFrequency_( averageFrequency )
{
// initially write header in file
......@@ -325,15 +326,17 @@ private:
};
class SphereEval
/*
* Functionality for continuous logging of sphere properties during moving simulation
*/
class MovingSpherePropertyEvaluator
{
public:
SphereEval( SweepTimeloop* timeloop,
const shared_ptr< StructuredBlockStorage > & blocks, const ConstBlockDataID & bodyStorageID,
const std::string & basefolder,
const Vector3<real_t> & u_infty, const real_t Galileo, const real_t GalileoSim,
const real_t gravity, const real_t viscosity, const real_t diameter, const real_t densityRatio,
const uint_t numLBMSubCycles ) :
MovingSpherePropertyEvaluator( SweepTimeloop* timeloop, const shared_ptr< StructuredBlockStorage > & blocks,
const ConstBlockDataID & bodyStorageID, const std::string & basefolder,
const Vector3<real_t> & u_infty, const real_t Galileo, const real_t GalileoSim,
const real_t gravity, const real_t viscosity, const real_t diameter,
const real_t densityRatio, const uint_t numLBMSubCycles ) :
timeloop_( timeloop ), blocks_( blocks ), bodyStorageID_( bodyStorageID ),
u_infty_( u_infty ), gravity_( gravity ), viscosity_( viscosity ),
diameter_( diameter ), densityRatio_( densityRatio ), numLBMSubCycles_( numLBMSubCycles )
......@@ -494,95 +497,9 @@ private:
};
class ResetForce
{
public:
ResetForce( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & bodyStorageID )
: blocks_( blocks ), bodyStorageID_( bodyStorageID ) { }
void operator()()
{
for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin< pe::Sphere >( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end<pe::Sphere>(); ++bodyIt )
{
bodyIt->resetForceAndTorque();
}
}
}
private:
const shared_ptr< StructuredBlockStorage > blocks_;
const BlockDataID bodyStorageID_;
};
class GravitationalForce
{
public:
GravitationalForce( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & bodyStorageID,
real_t force )
: blocks_( blocks ), bodyStorageID_( bodyStorageID ), force_( force ) { }
void operator()()
{
for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
{
for( auto bodyIt = pe::LocalBodyIterator::begin< pe::Sphere >( *blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end<pe::Sphere>(); ++bodyIt )
{
bodyIt->addForce ( real_t(0), real_t(0), force_ );
}
}
}
private:
const shared_ptr< StructuredBlockStorage > blocks_;
const BlockDataID bodyStorageID_;
const real_t force_;
};
/*
* When using N LBM steps and one (larger) pe step in a single simulation step, the force applied on the pe bodies in each LBM sep has to be scaled
* by a factor of 1/N before running the pe simulation. This corresponds to an averaging of the force and torque over the N LBM steps and is said to
* damp oscillations. Usually, N = 2.
* See Ladd - "Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation", 1994, p. 302
*/
class AverageForce
{
public:
AverageForce( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & bodyStorageID,
const uint_t lbmSubCycles )
: blocks_( blocks ), bodyStorageID_( bodyStorageID ), invLbmSubCycles_( real_t(1) / real_c( lbmSubCycles ) ) { }
void operator()()
{
Vector3<real_t> force( real_t(0) );
Vector3<real_t> torque( 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 )
{
force = invLbmSubCycles_ * bodyIt->getForce();
torque = invLbmSubCycles_ * bodyIt->getTorque();
bodyIt->resetForceAndTorque();
bodyIt->setForce ( force );
bodyIt->setTorque( torque );
}
}
}
private:
const shared_ptr< StructuredBlockStorage > blocks_;
const BlockDataID bodyStorageID_;
real_t invLbmSubCycles_;
};
/*
* Result evaluation
* Result evaluation for the written VTK files
*
* input: vel (fluid velocity), u_p (particle velocity), u_infty (inflow velocity)
*
* needed data:
......@@ -598,45 +515,37 @@ class AverageForce
* vel_r_perp = vel_r * e_p_perp
* vel_r_Hz_perp = vel_r * e_p_Hz_perp
*/
class LogVTKInfo
class VTKInfoLogger
{
public:
LogVTKInfo( SweepTimeloop* timeloop,
const shared_ptr< StructuredBlockStorage > & blocks, const ConstBlockDataID & bodyStorageID,
const std::string & baseFolder,
const Vector3<real_t> u_infty, const real_t gravity,
const real_t diameter, const real_t densityRatio ) :
VTKInfoLogger( SweepTimeloop* timeloop, const shared_ptr< StructuredBlockStorage > & blocks,
const ConstBlockDataID & bodyStorageID, const std::string & baseFolder,
const Vector3<real_t> u_infty ) :
timeloop_( timeloop ), blocks_( blocks ), bodyStorageID_( bodyStorageID ), baseFolder_( baseFolder ), u_infty_( u_infty )
{
u_p_r_sqr_mag_ = real_t(0);
u_ref_ = std::sqrt( std::fabs(densityRatio - real_t(1)) * gravity * diameter );
}
{ }
void operator()()
{
Vector3<real_t> transVel( real_t(0) );
Vector3<real_t> u_p( real_t(0) );
for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
{
for( auto bodyIt = pe::LocalBodyIterator::begin<pe::Sphere>(*blockIt, bodyStorageID_ ); bodyIt != pe::LocalBodyIterator::end<pe::Sphere>(); ++bodyIt )
{
transVel = bodyIt->getLinearVel();
u_p = bodyIt->getLinearVel();
}
}
WALBERLA_MPI_SECTION()
{
mpi::allReduceInplace( transVel[0], mpi::SUM );
mpi::allReduceInplace( transVel[1], mpi::SUM );
mpi::allReduceInplace( transVel[2], mpi::SUM );
mpi::allReduceInplace( u_p[0], mpi::SUM );
mpi::allReduceInplace( u_p[1], mpi::SUM );
mpi::allReduceInplace( u_p[2], mpi::SUM );
}
// store particle velocity
u_p_ = transVel;
Vector3<real_t> u_p_r = u_p_ - u_infty_;
u_p_r_sqr_mag_ = u_p_r.sqrLength();
Vector3<real_t> u_p_r = u_p - u_infty_;
real_t u_p_r_sqr_mag = u_p_r.sqrLength();
real_t u_p_H = std::sqrt( u_p_r[0] * u_p_r[0] + u_p_r[1] * u_p_r[1]);
......@@ -660,8 +569,8 @@ public:
file.open( filename.c_str() );
file.precision(8);
file << "u_p = " << u_p_ << "\n";
file << "u_p_r_2 = " << u_p_r_sqr_mag_ << "\n\n";
file << "u_p = " << u_p << "\n";
file << "u_p_r_2 = " << u_p_r_sqr_mag << "\n\n";
file << "e_{pH} = " << e_p_H << "\n";
file << "e_{pparallel} = " << e_p_parallel << "\n";
......@@ -670,22 +579,6 @@ public:
file.close();
}
}
Vector3<real_t> getParticleVel() const
{
return u_p_;
}
real_t getURef() const
{
return u_ref_;
}
real_t getUPRSqrMag() const
{
return u_p_r_sqr_mag_;
}
private:
......@@ -698,10 +591,6 @@ private:
std::string baseFolder_;
Vector3<real_t> u_infty_;
Vector3< real_t > u_p_;
real_t u_p_r_sqr_mag_;
real_t u_ref_;
};
class PDFFieldCopy
......@@ -869,7 +758,7 @@ int main( int argc, char **argv )
///////////////////////////
const int numProcs = MPIManager::instance()->numProcesses();
if( numProcs < 16 ) WALBERLA_ABORT("At least 16 MPI ranks have to be used due to horizontal periodicity requirements!");
WALBERLA_CHECK(numProcs % 16 != 0, "An integer multiple of 16 MPI ranks has to be used due to horizontal periodicity and domain decomposition requirements!");
uint_t XBlocks, YBlocks, ZBlocks;
if( numProcs >= 192 )
......@@ -1130,11 +1019,11 @@ int main( int argc, char **argv )
}
// evaluate the drag force acting on the sphere...
shared_ptr< ForceEval > forceEval = walberla::make_shared< ForceEval >( &timeloopInit, blocks, bodyStorageID, averageFrequency, basefolder );
timeloopInit.addFuncAfterTimeStep( SharedFunctor< ForceEval >( forceEval ), "Evaluating drag force" );
shared_ptr< RestingSphereForceEvaluator > forceEval = walberla::make_shared< RestingSphereForceEvaluator >( &timeloopInit, blocks, bodyStorageID, averageFrequency, basefolder );
timeloopInit.addFuncAfterTimeStep( SharedFunctor< RestingSphereForceEvaluator >( forceEval ), "Evaluating drag force" );
// ...then reset the force
timeloopInit.addFuncAfterTimeStep( ResetForce( blocks, bodyStorageID ), "Resetting force on body");
timeloopInit.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesResetter( blocks, bodyStorageID ), "Resetting force on body");
if( vtkIOInit )
{
......@@ -1300,8 +1189,7 @@ int main( int argc, char **argv )
}
if( numLBMSubCycles != uint_t(1) )
timeloop.addFuncAfterTimeStep( AverageForce( blocks, bodyStorageID, numLBMSubCycles ), "Force averaging for several LBM steps" );
timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesScaler( blocks, bodyStorageID, real_t(1) / real_c(numLBMSubCycles) ), "Force averaging for several LBM steps" );
}else{
......@@ -1341,15 +1229,16 @@ int main( int argc, char **argv )
}
if( numLBMSubCycles != uint_t(1) )
timeloop.addFuncAfterTimeStep( AverageForce( blocks, bodyStorageID, numLBMSubCycles ), "Force averaging for several LBM steps" );
timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesScaler( blocks, bodyStorageID, real_t(1) / real_c(numLBMSubCycles) ), "Force averaging for several LBM steps" );
}
// add gravity
timeloop.addFuncAfterTimeStep( GravitationalForce( blocks, bodyStorageID, - gravity * ( densityRatio - real_t(1) ) * diameter * diameter * diameter * math::PI / real_t(6) ), "Add gravitational force" );
Vector3<real_t> extForcesOnSphere( real_t(0), real_t(0), - gravity * ( densityRatio - real_t(1) ) * diameter * diameter * diameter * math::PI / real_t(6));
timeloop.addFuncAfterTimeStep( pe_coupling::ForceOnBodiesAdder( blocks, bodyStorageID, extForcesOnSphere ), "Add external forces (gravity and buoyancy)" );
// evaluate the sphere properties
timeloop.addFuncAfterTimeStep( SphereEval( &timeloop, blocks, bodyStorageID, basefolder, uInfty, Galileo, GalileoSim, gravity, viscosity, diameter, densityRatio, numLBMSubCycles ), "Evaluate sphere");
timeloop.addFuncAfterTimeStep( MovingSpherePropertyEvaluator( &timeloop, blocks, bodyStorageID, basefolder, uInfty, Galileo, GalileoSim, gravity, viscosity, diameter, densityRatio, numLBMSubCycles ), "Evaluate sphere");
// advance pe rigid body simulation
timeloop.addFuncAfterTimeStep( pe_coupling::TimeStep( blocks, bodyStorageID, cr, syncCall, real_c(numLBMSubCycles), numPeSubCycles ), "pe Time Step" );
......@@ -1363,8 +1252,7 @@ int main( int argc, char **argv )
pdfFieldVTK->addBeforeFunction( pdfGhostLayerSync );
// function to output plane infos for vtk output
shared_ptr< LogVTKInfo > logVTKInfo = walberla::make_shared< LogVTKInfo >( &timeloop, blocks, bodyStorageID, basefolder, uInfty, gravity, diameter, densityRatio );
pdfFieldVTK->addBeforeFunction( SharedFunctor< LogVTKInfo >( logVTKInfo ) );
pdfFieldVTK->addBeforeFunction(VTKInfoLogger( &timeloop, blocks, bodyStorageID, basefolder, uInfty ));
// create folder for log_vtk files to not pollute the basefolder
boost::filesystem::path tpath2( basefolder+"/log_vtk" );
......
src/pe_coupling/utility/BodiesForceTorqueContainer.h 0 → 100644
+ 142
− 0
View file @ b09aaf8b
//======================================================================================================================
//
// 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 BodiesForceTorqueContainer.h
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "core/math/Vector3.h"
#include "domain_decomposition/StructuredBlockStorage.h"
#include "pe/rigidbody/BodyIterators.h"
#include "pe/synchronization/SyncForces.h"
#include <map>
#include <vector>
namespace walberla {
namespace pe_coupling {
class BodiesForceTorqueContainer
{
public:
BodiesForceTorqueContainer( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & bodyStorageID )
: blockStorage_( blockStorage ), bodyStorageID_( bodyStorageID )
{ }
void operator()()
{
store();
}
void store()
{
// 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 bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
auto & f = bodyForceTorqueMap_[ 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();
}
}
}
}
void setOnBodies()
{
// owning process sets the force/torque on the bodies
for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
{
for( auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::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] );
}
}
}
void clear()
{
bodyForceTorqueMap_.clear();
}
void swap( BodiesForceTorqueContainer & other )
{
std::swap( this->bodyForceTorqueMap_, other.bodyForceTorqueMap_ );
}
private:
shared_ptr<StructuredBlockStorage> blockStorage_;
const BlockDataID bodyStorageID_;
std::map< walberla::id_t, std::vector<real_t> > bodyForceTorqueMap_;
};
class BodyContainerSwapper
{
public:
BodyContainerSwapper( const shared_ptr<BodiesForceTorqueContainer> & cont1,
const shared_ptr<BodiesForceTorqueContainer> & cont2 )
: cont1_( cont1 ), cont2_( cont2 )
{ }
void operator()()
{
cont1_->swap( *cont2_ );
}
private:
shared_ptr<BodiesForceTorqueContainer> cont1_;
shared_ptr<BodiesForceTorqueContainer> cont2_;
};
} // namespace pe_coupling
} // namespace walberla
src/pe_coupling/utility/ForceOnBodiesAdder.h 0 → 100644
+ 66
− 0
View file @ b09aaf8b
//======================================================================================================================
//
// 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 ForceOnBodiesAdder.h
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "core/math/Vector3.h"
#include "domain_decomposition/StructuredBlockStorage.h"
#include "pe/rigidbody/BodyIterators.h"
namespace walberla {
namespace pe_coupling {
class ForceOnBodiesAdder
{
public:
ForceOnBodiesAdder( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & bodyStorageID,
const Vector3<real_t> & force )
: blockStorage_( blockStorage ), bodyStorageID_( bodyStorageID ), force_( force )
{ }
// set a constant force on all (only local, to avoid force duplication) bodies
void operator()()
{
for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
{
for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
{
bodyIt->addForce ( force_ );
}
}
}
void updateForce( const Vector3<real_t> & newForce )
{
force_ = newForce;
}
private:
shared_ptr<StructuredBlockStorage> blockStorage_;
const BlockDataID bodyStorageID_;
Vector3<real_t> force_;
};
} // namespace pe_coupling
} // namespace walberla
src/pe_coupling/utility/ForceTorqueOnBodiesResetter.h 0 → 100644
+ 60
− 0
View file @ b09aaf8b
//======================================================================================================================
//
// 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 ForceTorqueOnBodiesResetter.h
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "core/math/Vector3.h"
#include "domain_decomposition/StructuredBlockStorage.h"
#include "pe/rigidbody/BodyIterators.h"
namespace walberla {
namespace pe_coupling {
class ForceTorqueOnBodiesResetter
{
public:
ForceTorqueOnBodiesResetter( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & bodyStorageID )
: blockStorage_( blockStorage ), bodyStorageID_( bodyStorageID )
{ }
// resets forces and torques on all (local and remote) bodies
void operator()()
{
for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
bodyIt->resetForceAndTorque();
}
}
}
private:
shared_ptr<StructuredBlockStorage> blockStorage_;
const BlockDataID bodyStorageID_;
};
} // namespace pe_coupling
} // namespace walberla
src/pe_coupling/utility/ForceTorqueOnBodiesScaler.h 0 → 100644
+ 72
− 0
View file @ b09aaf8b
//======================================================================================================================
//
// 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 ForceTorqueOnBodiesScaler.h
//! \ingroup pe_coupling
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "core/math/Vector3.h"
#include "domain_decomposition/StructuredBlockStorage.h"
#include "pe/rigidbody/BodyIterators.h"
namespace walberla {
namespace pe_coupling {
// scales force/torquew on all bodies (local and remote) by a constant scalar value
// can e.g. be used to average the force/torque over two time steps
class ForceTorqueOnBodiesScaler
{
public:
ForceTorqueOnBodiesScaler( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & bodyStorageID,
const real_t & scalingFactor )
: blockStorage_( blockStorage ), bodyStorageID_( bodyStorageID ), scalingFactor_( scalingFactor )
{ }
// resets forces and torques on all (local and remote) bodies
void operator()()
{
Vector3<real_t> force(real_t(0));
Vector3<real_t> torque(real_t(0));
for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
force = scalingFactor_ * bodyIt->getForce();
torque = scalingFactor_ * bodyIt->getTorque();
bodyIt->resetForceAndTorque();
bodyIt->setForce ( force );
bodyIt->setTorque( torque );
}
}
}
private:
shared_ptr<StructuredBlockStorage> blockStorage_;
const BlockDataID bodyStorageID_;
const real_t scalingFactor_;
};
} // namespace pe_coupling
} // namespace walberla
src/pe_coupling/utility/TimeStep.h
+ 101
− 83
View file @ b09aaf8b
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla is free software: you can
// 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
// 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
//
// 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/>.
//
......@@ -44,92 +44,110 @@ class TimeStep
{
public:
explicit TimeStep( const shared_ptr<StructuredBlockStorage> & blockStorage,
const BlockDataID & bodyStorageID,
pe::cr::ICR & collisionResponse,
boost::function<void (void)> synchronizeFunc,
const real_t timestep = real_t(1), const uint_t intermediateSteps = uint_c(1) )
: timestep_( timestep )
, intermediateSteps_( ( intermediateSteps == 0 ) ? uint_c(1) : intermediateSteps )
, blockStorage_( blockStorage )
, bodyStorageID_(bodyStorageID)
, collisionResponse_(collisionResponse)
, synchronizeFunc_(synchronizeFunc)
{}
void operator()()
{
if( intermediateSteps_ == 1 )
{
collisionResponse_.timestep( timestep_ );
synchronizeFunc_( );
}
else
{
// sum up all forces from shadow copies on local body (owner)
pe::reduceForces( blockStorage_->getBlockStorage(), bodyStorageID_ );
// send total forces to shadow owners
pe::distributeForces( blockStorage_->getBlockStorage(), bodyStorageID_ );
// generate map from all known bodies (process local) to total forces
std::map< walberla::id_t, std::vector< real_t > > forceMap;
for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
explicit TimeStep( const shared_ptr<StructuredBlockStorage> & blockStorage,
const BlockDataID & bodyStorageID,
pe::cr::ICR & collisionResponse,
boost::function<void (void)> synchronizeFunc,
const real_t timeStep = real_t(1), const uint_t intermediateSteps = uint_c(1) )
: timeStep_( timeStep )
, intermediateSteps_( ( intermediateSteps == 0 ) ? uint_c(1) : intermediateSteps )
, blockStorage_( blockStorage )
, bodyStorageID_(bodyStorageID)
, collisionResponse_(collisionResponse)
, synchronizeFunc_(synchronizeFunc)
{}
void operator()()
{
if( intermediateSteps_ == 1 )
{
collisionResponse_.timestep( timeStep_ );
synchronizeFunc_( );
}
else
{
// during the intermediate time steps of the collision response, the currently acting forces
// (interaction forces, gravitational force, ...) have to remain constant.
// 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
std::map< walberla::id_t, std::vector< real_t > > forceMap;
for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
{
// 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 owner ship (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 stores as well.
for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
auto & f = forceMap[ bodyIt->getSystemID() ];
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] );
if ( bodyIt->isRemote() )
{
bodyIt->resetForceAndTorque();
}
}
auto & f = forceMap[ 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();
}
}
// perform pe time steps
const real_t subTimestep = timestep_ / real_c( intermediateSteps_ );
for( uint_t i = 0; i != intermediateSteps_; ++i )
}
// perform pe time steps
const real_t subTimestep = timeStep_ / real_c( intermediateSteps_ );
for( uint_t i = 0; i != intermediateSteps_; ++i )
{
// in the first set forces on local bodies are already set by force synchronization
if( i != 0 )
{
// in the first set forces on local bodies are already set by force synchronization
if( i != 0 ) {
for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
{
for( auto body = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID_); body != pe::LocalBodyIterator::end(); ++body )
{
const auto & f = forceMap[ body->getSystemID() ];
body->addForce ( f[0], f[1], f[2] );
body->addTorque( f[3], f[4], f[5] );
}
}
}
collisionResponse_.timestep( subTimestep );
synchronizeFunc_( );
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 )
{
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] );
}
}
}
}
}
collisionResponse_.timestep( subTimestep );
synchronizeFunc_( );
}
}
}
protected:
const real_t timestep_;
const real_t timeStep_;
const uint_t intermediateSteps_;
const uint_t intermediateSteps_;
shared_ptr<StructuredBlockStorage> blockStorage_;
const BlockDataID & bodyStorageID_;
shared_ptr<StructuredBlockStorage> blockStorage_;
const BlockDataID & bodyStorageID_;
pe::cr::ICR & collisionResponse_;
boost::function<void (void)> synchronizeFunc_;
pe::cr::ICR & collisionResponse_;
boost::function<void (void)> synchronizeFunc_;
}; // class TimeStep
......
src/pe_coupling/utility/all.h
+ 4
− 0
View file @ b09aaf8b
......@@ -22,5 +22,9 @@
#pragma once
#include "BodiesForceTorqueContainer.h"
#include "ForceOnBodiesAdder.h"
#include "ForceTorqueOnBodiesResetter.h"
#include "ForceTorqueOnBodiesScaler.h"
#include "LubricationCorrection.h"
#include "TimeStep.h"
tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
+ 2
− 25
View file @ b09aaf8b
......@@ -58,6 +58,7 @@
#include "pe_coupling/mapping/all.h"
#include "pe_coupling/momentum_exchange_method/all.h"
#include "pe_coupling/utility/all.h"
#include <vector>
#include <iomanip>
......@@ -331,30 +332,6 @@ private:
};
class ResetForce
{
public:
ResetForce( const shared_ptr< StructuredBlockStorage > & blocks,
const BlockDataID & bodyStorageID )
: blocks_( blocks ), bodyStorageID_( bodyStorageID )
{ }
void operator()()
{
for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
bodyIt->resetForceAndTorque();
}
}
}
private:
shared_ptr< StructuredBlockStorage > blocks_;
const BlockDataID bodyStorageID_;
};
class PDFFieldCopy
{
public:
......@@ -583,7 +560,7 @@ int main( int argc, char **argv )
<< AfterFunction( SharedFunctor< DragForceEvaluator >( forceEval ), "drag force evaluation" );
// resetting force
timeloop.addFuncAfterTimeStep( ResetForce( blocks, bodyStorageID ), "reset force on sphere");
timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesResetter( blocks, bodyStorageID ), "reset force on sphere");
timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
......
tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
+ 2
− 25
View file @ b09aaf8b
......@@ -61,6 +61,7 @@
#include "pe_coupling/mapping/all.h"
#include "pe_coupling/momentum_exchange_method/all.h"
#include "pe_coupling/utility/all.h"
#include "vtk/all.h"
#include "field/vtk/all.h"
......@@ -405,30 +406,6 @@ class ForceEval
};
class ResetForce
{
public:
ResetForce( const shared_ptr< StructuredBlockStorage > & blocks,
const BlockDataID & bodyStorageID )
: blocks_( blocks ), bodyStorageID_( bodyStorageID )
{ }
void operator()()
{
for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
bodyIt->resetForceAndTorque();
}
}
}
private:
shared_ptr< StructuredBlockStorage > blocks_;
const BlockDataID bodyStorageID_;
};
//////////
// MAIN //
//////////
......@@ -596,7 +573,7 @@ int main( int argc, char **argv )
timeloop.addFuncAfterTimeStep( makeSharedFunctor( forceEval ), "drag force evaluation" );
// resetting force
timeloop.addFuncAfterTimeStep( ResetForce( blocks, bodyStorageID ), "reset force on sphere");
timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesResetter( blocks, bodyStorageID ), "reset force on sphere");
if( vtkIO )
{
......
tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
+ 2
− 25
View file @ b09aaf8b
......@@ -58,6 +58,7 @@
#include "pe_coupling/mapping/all.h"
#include "pe_coupling/momentum_exchange_method/all.h"
#include "pe_coupling/utility/all.h"
#include <vector>
#include <iomanip>
......@@ -276,30 +277,6 @@ class TorqueEval
};
class ResetForce
{
public:
ResetForce( const shared_ptr< StructuredBlockStorage > & blocks,
const BlockDataID & bodyStorageID )
: blocks_( blocks ), bodyStorageID_( bodyStorageID )
{ }
void operator()()
{
for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
bodyIt->resetForceAndTorque();
}
}
}
private:
shared_ptr< StructuredBlockStorage > blocks_;
const BlockDataID bodyStorageID_;
};
class PDFFieldCopy
{
public:
......@@ -539,7 +516,7 @@ int main( int argc, char **argv )
<< AfterFunction( SharedFunctor< TorqueEval >( torqueEval ), "torque evaluation" );
}
// resetting force
timeloop.addFuncAfterTimeStep( ResetForce( blocks, bodyStorageID ), "reset force on sphere");
timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesResetter( blocks, bodyStorageID ), "reset force on sphere");
timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment