Commit 8bb2c3c4 authored by Christoph Rettinger's avatar Christoph Rettinger

Added hydrodynamic force notifications, introduced them to tests and...

Added hydrodynamic force notifications, introduced them to tests and benchmarks, minor documentation and include fixes
parent c98cac51
Pipeline #22381 passed with stages
in 230 minutes and 41 seconds
......@@ -70,7 +70,7 @@
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
#include "mesa_pd/vtk/ParticleVtkOutput.h"
#include "timeloop/SweepTimeloop.h"
......@@ -207,8 +207,11 @@ public:
size_t idx = ac_->uidToIdx(sphereUid_);
if( idx != ac_->getInvalidIdx())
{
force = ac_->getHydrodynamicForce(idx);
torque = ac_->getHydrodynamicTorque(idx);
if(!mesa_pd::data::particle_flags::isSet(ac_->getFlags(idx),mesa_pd::data::particle_flags::GHOST))
{
force = ac_->getHydrodynamicForce(idx);
torque = ac_->getHydrodynamicTorque(idx);
}
}
WALBERLA_MPI_SECTION()
......@@ -549,9 +552,11 @@ int main( int argc, char **argv )
mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
syncNextNeighborFunc(*ps, *rpdDomain, overlap);
};
syncCall();
mesa_pd::mpi::ReduceProperty reduceProperty;
///////////////////////
// ADD DATA TO BLOCKS //
////////////////////////
......@@ -714,13 +719,16 @@ int main( int argc, char **argv )
// perform a single simulation step
timeloop.singleStep( timeloopTiming );
// sync hydrodynamic force to local particle
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
// average force
if( timestep == 0 )
{
lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueForAveragingKernel;
ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
ps->forEachParticle(false, mesa_pd::kernel::SelectLocal(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
}
ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor, averageHydrodynamicForceTorque, *accessor );
ps->forEachParticle(false, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque, *accessor );
// evaluation
timeloopTiming["Logging"].start();
......
......@@ -52,9 +52,6 @@
#include "lbm_mesapd_coupling/mapping/ParticleMapping.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/MovingParticleMapping.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/boundary/CurvedLinear.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/Reconstructor.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/ExtrapolationDirectionFinder.h"
#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/PdfReconstructionManager.h"
#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/LubricationCorrectionKernel.h"
......@@ -69,13 +66,9 @@
#include "mesa_pd/data/shape/Sphere.h"
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
#include "mesa_pd/kernel/SpringDashpot.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ContactFilter.h"
#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
#include "mesa_pd/vtk/ParticleVtkOutput.h"
#include "timeloop/SweepTimeloop.h"
......@@ -536,11 +529,8 @@ int main( int argc, char **argv )
syncCall();
mesa_pd::mpi::ReduceProperty reduceProperty;
lbm_mesapd_coupling::ResetHydrodynamicForceTorqueKernel resetHydrodynamicForceTorque;
real_t lubricationCutOffDistanceNormal = real_t(2) / real_t(3);
real_t lubricationCutOffDistanceTangentialTranslational = real_t(0.5);
real_t lubricationCutOffDistanceTangentialRotational = real_t(0.5);
......@@ -672,7 +662,6 @@ int main( int argc, char **argv )
{
if (contact_filter(acd.getIdx1(), acd.getIdx2(), *accessor, acd.getContactPoint(), *rpdDomain))
{
//WALBERLA_LOG_INFO(acd.getIdx1() << " " << acd.getIdx2() << " " << acd.getContactNormal() << " " << acd.getPenetrationDepth());
double_cast(acd.getIdx1(), acd.getIdx2(), *accessor, lubricationCorrectionKernel, *accessor, acd.getContactNormal(), acd.getPenetrationDepth());
}
}
......
......@@ -82,7 +82,6 @@
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/kernel/NonLinearSpringDashpot.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
......@@ -90,6 +89,7 @@
#include "mesa_pd/mpi/ReduceContactHistory.h"
#include "mesa_pd/mpi/ContactFilter.h"
#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
#include "mesa_pd/vtk/ParticleVtkOutput.h"
#include "timeloop/SweepTimeloop.h"
......@@ -818,11 +818,6 @@ int main( int argc, char **argv )
BlockDataID omegaBulkFieldID = field::addToStorage<ScalarField_T>( blocks, "omega bulk field", omegaBulk, field::fzyx, uint_t(0) );
// create the lattice model
// TRT
//LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( real_t(1) / relaxationTime ) );
// generated MRT
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
real_t lambda_d = lbm::collision_model::TRT::lambda_d( omega, magicNumber );
#ifdef WALBERLA_BUILD_WITH_CODEGEN
......@@ -868,11 +863,7 @@ int main( int argc, char **argv )
//linearCollisionResponse.setFrictionCoefficientStatic(0,0,frictionCoeff); // not used in this test case
linearCollisionResponse.setFrictionCoefficientDynamic(0,0,frictionCoeff);
// nonlinear model for ACTM
mesa_pd::mpi::ReduceProperty reduceProperty;
mesa_pd::mpi::ReduceContactHistory reduceAndSwapContactHistory;
// set up coupling functionality
......@@ -1057,15 +1048,17 @@ int main( int argc, char **argv )
{
// perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions
timeloop.singleStep( timeloopTiming );
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
if( averageForceTorqueOverTwoTimeSteps )
{
if( i == 0 )
{
lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueForAveragingKernel;
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
}
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, averageHydrodynamicForceTorque, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque, *accessor );
}
Vector3<real_t> hydForce(real_t(0));
......@@ -1165,7 +1158,7 @@ int main( int argc, char **argv )
// add hydrodynamic force
lbm_mesapd_coupling::AddHydrodynamicInteractionKernel addHydrodynamicInteraction;
ps->forEachParticle(useOpenMP, lbm_mesapd_coupling::RegularParticlesSelector(), *accessor, addHydrodynamicInteraction, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addHydrodynamicInteraction, *accessor );
hydForce = getForce(sphereUid,*accessor) - lubForce - collisionForce;
WALBERLA_ASSERT(!std::isnan(hydForce[0]) && !std::isnan(hydForce[1]) && !std::isnan(hydForce[2]), "Found nan value in hydrodynamic force = " << hydForce);
......
......@@ -74,12 +74,12 @@
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/SpringDashpot.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ContactFilter.h"
#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
#include "mesa_pd/vtk/ParticleVtkOutput.h"
#include "timeloop/SweepTimeloop.h"
......@@ -216,8 +216,8 @@ public:
{
pos = ac_->getPosition(idx);
transVel = ac_->getLinearVelocity(idx);
hydForce = ac_->getHydrodynamicForce(idx);
}
hydForce = ac_->getHydrodynamicForce(idx);
}
WALBERLA_MPI_SECTION()
......@@ -643,7 +643,6 @@ int main( int argc, char **argv )
mesa_pd::kernel::VelocityVerletWithShapePreForceUpdate vvIntegratorPreForce(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::VelocityVerletWithShapePostForceUpdate vvIntegratorPostForce(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::SpringDashpot collisionResponse(1);
mesa_pd::mpi::ReduceProperty reduceProperty;
// set up coupling functionality
......@@ -811,15 +810,17 @@ int main( int argc, char **argv )
timeloop.singleStep( timeloopTiming );
timeloopTiming["RPD"].start();
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
if( averageForceTorqueOverTwoTimeSteps )
{
if( i == 0 )
{
lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueForAveragingKernel;
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
}
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, averageHydrodynamicForceTorque, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque, *accessor );
}
for(auto subCycle = uint_t(0); subCycle < numRPDSubCycles; ++subCycle )
......
......@@ -78,12 +78,11 @@
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/SpringDashpot.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ContactFilter.h"
#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
#include "mesa_pd/vtk/ParticleVtkOutput.h"
#include "timeloop/SweepTimeloop.h"
......@@ -219,14 +218,13 @@ public:
if(!mesa_pd::data::particle_flags::isSet( ac_->getFlags(idx), mesa_pd::data::particle_flags::GHOST))
{
pos = ac_->getPosition(idx)[2];
hydForce = ac_->getHydrodynamicForce(idx)[2];
}
hydForce = ac_->getHydrodynamicForce(idx)[2];
}
WALBERLA_MPI_SECTION()
{
mpi::allReduceInplace( pos, mpi::SUM );
mpi::allReduceInplace( hydForce, mpi::SUM );
}
......@@ -601,7 +599,6 @@ int main( int argc, char **argv )
mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::SpringDashpot collisionResponse(1);
mesa_pd::mpi::ReduceProperty reduceProperty;
// set up coupling functionality
......@@ -768,15 +765,17 @@ int main( int argc, char **argv )
timeloop.singleStep( timeloopTiming );
timeloopTiming["RPD"].start();
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
if( averageForceTorqueOverTwoTimeSteps )
{
if( i == 0 )
{
lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueForAveragingKernel;
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
}
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, averageHydrodynamicForceTorque, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque, *accessor );
}
reduceProperty.operator()<mesa_pd::ForceTorqueNotification>(*ps);
......
......@@ -85,11 +85,13 @@
#include "mesa_pd/kernel/NonLinearSpringDashpot.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h"
#include "mesa_pd/mpi/ClearNextNeighborSync.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
#include "mesa_pd/mpi/ContactFilter.h"
#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
#include "mesa_pd/vtk/ParticleVtkOutput.h"
#include "timeloop/SweepTimeloop.h"
......@@ -539,7 +541,7 @@ real_t getAverageDensityInSystem(const shared_ptr<StructuredBlockStorage> & bloc
//////////
//*******************************************************************************************************************
/*!\brief PHYSICAL Test case of a sphere settling and colliding with a wall submerged in a viscous fluid.
/*!\brief PHYSICAL test case of a sphere settling and colliding with a wall submerged in a viscous fluid.
*
* The trajectory of the bouncing sphere is logged and compared to reference experiments.
*
......@@ -802,13 +804,28 @@ int main( int argc, char **argv )
"Unmatching domain decomposition in direction " << i << "!");
}
auto blocks = blockforest::createUniformBlockGrid( numberOfBlocksPerDirection[0], numberOfBlocksPerDirection[1], numberOfBlocksPerDirection[2],
shared_ptr< StructuredBlockForest > blocks;
if( readFromCheckPointFile )
{
WALBERLA_LOG_INFO_ON_ROOT("Reading block forest from file!");
blocks = blockforest::createUniformBlockGrid(checkPointFileName+"_forest.txt", cellsPerBlockPerDirection[0], cellsPerBlockPerDirection[1], cellsPerBlockPerDirection[2], false);
}
else
{
blocks = blockforest::createUniformBlockGrid( numberOfBlocksPerDirection[0], numberOfBlocksPerDirection[1], numberOfBlocksPerDirection[2],
cellsPerBlockPerDirection[0], cellsPerBlockPerDirection[1], cellsPerBlockPerDirection[2], dx,
0, false, false,
true, true, false, //periodicity
false );
WALBERLA_LOG_INFO_ON_ROOT("Domain decomposition:");
if(writeCheckPointFile)
{
WALBERLA_LOG_INFO_ON_ROOT("Writing block forest to file!");
blocks->getBlockForest().saveToFile(checkPointFileName+"_forest.txt");
}
}
WALBERLA_LOG_INFO_ON_ROOT("Domain partitioning:");
WALBERLA_LOG_INFO_ON_ROOT(" - blocks per direction = " << numberOfBlocksPerDirection );
WALBERLA_LOG_INFO_ON_ROOT(" - cells per block = " << cellsPerBlockPerDirection );
......@@ -827,6 +844,9 @@ int main( int argc, char **argv )
//init data structures
auto ps = walberla::make_shared<mesa_pd::data::ParticleStorage>(1);
auto ss = walberla::make_shared<mesa_pd::data::ShapeStorage>();
using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
BlockDataID particleStorageID;
if(readFromCheckPointFile)
{
......@@ -834,12 +854,10 @@ int main( int argc, char **argv )
particleStorageID = blocks->loadBlockData( checkPointFileName + "_mesa.txt", mesa_pd::domain::createBlockForestDataHandling(ps), "Particle Storage" );
} else {
particleStorageID = blocks->addBlockData(mesa_pd::domain::createBlockForestDataHandling(ps), "Particle Storage");
mesa_pd::mpi::ClearNextNeighborSync CNNS;
CNNS(*accessor);
}
auto ss = walberla::make_shared<mesa_pd::data::ShapeStorage>();
using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
// bounding planes
createPlaneSetup(ps,ss,blocks->getDomain(), applyOutflowBCAtTop);
......@@ -1183,15 +1201,17 @@ int main( int argc, char **argv )
{
// perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions
timeloop.singleStep( timeloopTiming );
reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
if( averageForceTorqueOverTwoTimeSteps )
{
if( i == 0 )
{
lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueForAveragingKernel;
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
}
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *accessor, averageHydrodynamicForceTorque, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque, *accessor );
}
real_t hydForce(real_t(0));
......
......@@ -52,10 +52,10 @@ if __name__ == '__main__':
ps.add_property("dw", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
# Properties for lbm_mesapd_coupling:
ps.add_property("hydrodynamicForce", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
ps.add_property("hydrodynamicTorque", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
ps.add_property("oldHydrodynamicForce", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
ps.add_property("oldHydrodynamicTorque", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
ps.add_property("hydrodynamicForce", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
ps.add_property("hydrodynamicTorque", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
ps.add_property("oldHydrodynamicForce", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
ps.add_property("oldHydrodynamicTorque", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
ch = mpd.add(data.ContactHistory())
ch.add_property("tangentialSpringDisplacement", "walberla::mesa_pd::Vec3", defValue="real_t(0)")
......@@ -111,6 +111,9 @@ if __name__ == '__main__':
ftn = mpd.add(mpi.PropertyNotification('ForceTorqueNotification'))
ftn.add_property('force', 'mesa_pd::Vec3', 'Vec3(real_t(0))')
ftn.add_property('torque', 'mesa_pd::Vec3', 'Vec3(real_t(0))')
hftn = mpd.add(mpi.PropertyNotification('HydrodynamicForceTorqueNotification'))
hftn.add_property('hydrodynamicForce', 'mesa_pd::Vec3', 'Vec3(real_t(0))')
hftn.add_property('hydrodynamicTorque', 'mesa_pd::Vec3', 'Vec3(real_t(0))')
hfn = mpd.add(mpi.PropertyNotification('HeatFluxNotification'))
hfn.add_property('heatFlux', 'real_t', 'real_t(0)')
mpd.add(mpi.ReduceContactHistory())
......
......@@ -30,7 +30,8 @@ namespace lbm_mesapd_coupling {
/*
* Kernel that adds the current hydrodynamic forces and torques onto the particles as forces and torques
*
* Should usually be carried out on local and ghost particles.
* When reducing hyd. force/torque before, this should usually be carried out only on local particles. (recommended)
* If not, it must be carried out on local and ghost particles.
*/
class AddHydrodynamicInteractionKernel
{
......
......@@ -29,9 +29,10 @@ namespace lbm_mesapd_coupling {
/*
* Kernel that averages the force/torque over two time steps in a rolling average fashion.
* This is said to damp oscillations of the interaction force/torque.
* See Ladd - " Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation", 1994, p.
* See Ladd - "Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation", 1994, p.
*
* Should usually be carried out on local and ghost particles.
* When reducing hyd. force/torque before, this should usually be carried out only on local particles. (recommended)
* If not, it must be carried out on local and ghost particles.
*/
class AverageHydrodynamicForceTorqueKernel
{
......
......@@ -30,8 +30,6 @@ namespace lbm_mesapd_coupling {
* Kernel that initializes the old hydrodynamic force/torque property of a particle with the currently set one.
* This should be used when starting the simulation (from anew or from checkpoint) and after load balancing.
* Only then, the following averaging kernel (AverageHydrodynamicForceTorqueKernel) applies the correct amount of force.
*
* Should usually be carried out on local and ghost particles.
*/
class InitializeHydrodynamicForceTorqueForAveragingKernel
{
......
......@@ -30,8 +30,6 @@ namespace lbm_mesapd_coupling {
/*
* Kernel that resets the values of hydrodynamicForce and hydrodynamicTorque currently stored to zero
*
* Should usually be carried out on local and ghost particles.
*/
class ResetHydrodynamicForceTorqueKernel
{
......
//======================================================================================================================
//
// 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 HydrodynamicForceTorqueNotification.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/ParticleStorage.h>
#include <mesa_pd/mpi/notifications/NotificationType.h>
#include <mesa_pd/mpi/notifications/reset.h>
#include <core/mpi/Datatype.h>
#include <core/mpi/RecvBuffer.h>
#include <core/mpi/SendBuffer.h>
namespace walberla {
namespace mesa_pd {
/**
* Trasmits force and torque information.
*/
class HydrodynamicForceTorqueNotification
{
public:
struct Parameters
{
id_t uid_;
mesa_pd::Vec3 hydrodynamicForce_;
mesa_pd::Vec3 hydrodynamicTorque_;
};
inline explicit HydrodynamicForceTorqueNotification( const data::Particle& p ) : p_(p) {}
const data::Particle& p_;
};
template <>
void reset<HydrodynamicForceTorqueNotification>(data::Particle& p)
{
p.setHydrodynamicForce( Vec3(real_t(0)) );
p.setHydrodynamicTorque( Vec3(real_t(0)) );
}
void reduce(data::Particle&& p, const HydrodynamicForceTorqueNotification::Parameters& objparam)
{
p.getHydrodynamicForceRef() += objparam.hydrodynamicForce_;
p.getHydrodynamicTorqueRef() += objparam.hydrodynamicTorque_;
}
void update(data::Particle&& p, const HydrodynamicForceTorqueNotification::Parameters& objparam)
{
p.setHydrodynamicForce( objparam.hydrodynamicForce_ );
p.setHydrodynamicTorque( objparam.hydrodynamicTorque_ );
}
} // namespace mesa_pd
} // namespace walberla
//======================================================================================================================
//
// Send/Recv Buffer Serialization Specialization
//
//======================================================================================================================
namespace walberla {
namespace mpi {
template< typename T, // Element type of SendBuffer
typename G> // Growth policy of SendBuffer
mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, const mesa_pd::HydrodynamicForceTorqueNotification& obj )
{
buf.addDebugMarker( "pn" );
buf << obj.p_.getUid();
buf << obj.p_.getHydrodynamicForce();
buf << obj.p_.getHydrodynamicTorque();
return buf;
}
template< typename T> // Element type of RecvBuffer
mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd::HydrodynamicForceTorqueNotification::Parameters& objparam )
{
buf.readDebugMarker( "pn" );
buf >> objparam.uid_;
buf >> objparam.hydrodynamicForce_;
buf >> objparam.hydrodynamicTorque_;
return buf;
}
template< >
struct BufferSizeTrait< mesa_pd::HydrodynamicForceTorqueNotification > {
static const bool constantSize = true;
static const uint_t size = BufferSizeTrait<id_t>::size +
BufferSizeTrait<mesa_pd::Vec3>::size +
BufferSizeTrait<mesa_pd::Vec3>::size +
mpi::BUFFER_DEBUG_OVERHEAD;
};
} // mpi
} // walberla
\ No newline at end of file
......@@ -148,6 +148,10 @@ void ParseMessage::operator()(int sender,
pIt->setGhostOwners(objparam.ghostOwners_);
pIt->setOldForce(objparam.oldForce_);
pIt->setOldTorque(objparam.oldTorque_);
pIt->setHydrodynamicForce(objparam.hydrodynamicForce_);
pIt->setHydrodynamicTorque(objparam.hydrodynamicTorque_);
pIt->setOldHydrodynamicForce(objparam.oldHydrodynamicForce_);
pIt->setOldHydrodynamicTorque(objparam.oldHydrodynamicTorque_);
WALBERLA_LOG_DETAIL( "Processed PARTICLE_MIGRATION_NOTIFICATION." );
......
......@@ -64,6 +64,10 @@ public:
uint_t type {0};
std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory> oldContactHistory {};
walberla::real_t temperature {real_t(0)};
walberla::mesa_pd::Vec3 hydrodynamicForce {real_t(0)};
walberla::mesa_pd::Vec3 hydrodynamicTorque {real_t(0)};
walberla::mesa_pd::Vec3 oldHydrodynamicForce {real_t(0)};
walberla::mesa_pd::Vec3 oldHydrodynamicTorque {real_t(0)};
};
inline explicit ParticleCopyNotification( const data::Particle& particle ) : particle_(particle) {}
......@@ -91,6 +95,10 @@ inline data::ParticleStorage::iterator createNewParticle(data::ParticleStorage&
pIt->setType(data.type);
pIt->setOldContactHistory(data.oldContactHistory);
pIt->setTemperature(data.temperature);
pIt->setHydrodynamicForce(data.hydrodynamicForce);
pIt->setHydrodynamicTorque(data.hydrodynamicTorque);
pIt->setOldHydrodynamicForce(data.oldHydrodynamicForce);
pIt->setOldHydrodynamicTorque(data.oldHydrodynamicTorque);
return pIt;
}
......@@ -133,6 +141,10 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
buf << obj.particle_.getType();
buf << obj.particle_.getOldContactHistory();
buf << obj.particle_.getTemperature();
buf << obj.particle_.getHydrodynamicForce();
buf << obj.particle_.getHydrodynamicTorque();
buf << obj.particle_.getOldHydrodynamicForce();
buf << obj.particle_.getOldHydrodynamicTorque();
return buf;
}
......@@ -156,6 +168,10 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
buf >> objparam.type;
buf >> objparam.oldContactHistory;
buf >> objparam.temperature;
buf >> objparam.hydrodynamicForce;
buf >> objparam.hydrodynamicTorque;
buf >> objparam.oldHydrodynamicForce;
buf >> objparam.oldHydrodynamicTorque;
return buf;
}
......
......@@ -48,6 +48,10 @@ public:
std::unordered_set<walberla::mpi::MPIRank> ghostOwners_ {};
walberla::mesa_pd::Vec3 oldForce_ {real_t(0)};
walberla::mesa_pd::Vec3 oldTorque_ {real_t(0)};
walberla::mesa_pd::Vec3 hydrodynamicForce_ {real_t(0)};
walberla::mesa_pd::Vec3 hydrodynamicTorque_ {real_t(0)};
walberla::mesa_pd::Vec3 oldHydrodynamicForce_ {real_t(0)};
walberla::mesa_pd::Vec3 oldHydrodynamicTorque_ {real_t(0)};
};
inline explicit ParticleMigrationNotification( const data::Particle& particle ) : particle_(particle) {}
......@@ -81,6 +85,10 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
buf << obj.particle_.getGhostOwners();
buf << obj.particle_.getOldForce();
buf << obj.particle_.getOldTorque();
buf << obj.particle_.getHydrodynamicForce();
buf << obj.particle_.getHydrodynamicTorque();
buf << obj.particle_.getOldHydrodynamicForce();
buf << obj.particle_.getOldHydrodynamicTorque();
return buf;
}
......@@ -92,6 +100,10 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
buf >> objparam.ghostOwners_;
buf >> objparam.oldForce_;
buf >> objparam.oldTorque_;
buf >> objparam.hydrodynamicForce_;
buf >> objparam.hydrodynamicTorque_;
buf >> objparam.oldHydrodynamicForce_;
buf >> objparam.oldHydrodynamicTorque_