diff --git a/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp b/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
index 93da4ad7e225d7523530e4e5decf3ee014480fa0..ecf14716ca36f464bf00c482e0a83b8c63c41856 100644
--- a/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
@@ -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();
diff --git a/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py b/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
index bf127f8eb0ca4ddbaf74e5e950477a39eb8432cc..461291dbd48f83f17c49783740666f5cbbec5955 100644
--- a/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
+++ b/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
@@ -5,15 +5,26 @@ from lbmpy_walberla import generate_lattice_model
from pystencils_walberla import CodeGeneration
from lbmpy.creationfunctions import create_lb_collision_rule
-from lbmpy.moments import is_even, get_order
+from lbmpy.moments import is_even, get_order, MOMENT_SYMBOLS
from lbmpy.stencils import get_stencil
-from lbmpy.methods import mrt_orthogonal_modes_literature
with CodeGeneration() as ctx:
stencil = get_stencil("D3Q19", 'walberla')
omega = sp.symbols("omega_:%d" % len(stencil))
- moments = mrt_orthogonal_modes_literature(stencil, True, False)
+
+ x, y, z = MOMENT_SYMBOLS
+ one = sp.Rational(1, 1)
+ sq = x ** 2 + y ** 2 + z ** 2
+ moments = [
+ [one, x, y, z], # [0, 3, 5, 7]
+ [sq - 1], # [1]
+ [3 * sq ** 2 - 6 * sq + 1], # [2]
+ [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
+ [3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 11, 13, 14, 15]
+ [(2 * sq - 3) * (3 * x ** 2 - sq), (2 * sq - 3) * (y ** 2 - z ** 2)], # [10, 12]
+ [(y ** 2 - z ** 2) * x, (z ** 2 - x ** 2) * y, (x ** 2 - y ** 2) * z] # [16, 17, 18]
+ ]
method = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False, nested_moments=moments)
def modification_func(moment, eq, rate):
diff --git a/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py b/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
index d52d67dfc738cc6403f2b246d9cccb61a0ffd914..9293f976fc96c424cf685dc1d549e5bdaa0d7d61 100644
--- a/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
+++ b/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
@@ -5,10 +5,9 @@ from lbmpy_walberla import generate_lattice_model
from pystencils_walberla import CodeGeneration
from lbmpy.creationfunctions import create_lb_collision_rule
-from lbmpy.moments import is_even, get_order
+from lbmpy.moments import is_even, get_order, MOMENT_SYMBOLS
from lbmpy.stencils import get_stencil
from lbmpy.forcemodels import Luo
-from lbmpy.methods import mrt_orthogonal_modes_literature
with CodeGeneration() as ctx:
@@ -18,7 +17,19 @@ with CodeGeneration() as ctx:
stencil = get_stencil("D3Q19", 'walberla')
omega = sp.symbols("omega_:%d" % len(stencil))
- moments = mrt_orthogonal_modes_literature(stencil, True, False)
+ x, y, z = MOMENT_SYMBOLS
+ one = sp.Rational(1, 1)
+ sq = x ** 2 + y ** 2 + z ** 2
+ moments = [
+ [one, x, y, z], # [0, 3, 5, 7]
+ [sq - 1], # [1]
+ [3 * sq ** 2 - 6 * sq + 1], # [2]
+ [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
+ [3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 11, 13, 14, 15]
+ [(2 * sq - 3) * (3 * x ** 2 - sq), (2 * sq - 3) * (y ** 2 - z ** 2)], # [10, 12]
+ [(y ** 2 - z ** 2) * x, (z ** 2 - x ** 2) * y, (x ** 2 - y ** 2) * z] # [16, 17, 18]
+ ]
+
methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False,
force_model=forcemodel, nested_moments=moments)
diff --git a/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp b/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp
index 0be7fe660efa3adc5d4199bdb8f85e7fb9c00e0c..dcb7737fe5c8809802e22a3bd5802fd8321a2369 100644
--- a/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp
@@ -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());
}
}
diff --git a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
index 0e1aa9653c7687c7087f95479a721b7c76999880..dceb28ec8148f3c6d051fa83d37c965be57e9ace 100644
--- a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
@@ -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);
diff --git a/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp b/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
index d0d186c16092a948f047f1ac6844af2a0052e3c8..a0423c979d33650740333653cf361bc95948b0cf 100644
--- a/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
@@ -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 )
diff --git a/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp b/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
index eac6530d09e73028be048bc494c7cda60b8bd1b7..86dad96117e972522f3da3433ac725c294bc1794 100644
--- a/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
@@ -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);
diff --git a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
index 594cfd4325af59244525e1b570c6c4aa9891e7bb..676390ebe8fe1eb353b119d22a765ad61af19de0 100644
--- a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
@@ -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));
diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index 42572d41ffaaf44f7edb932c918678ae8d311b33..67da03b5162b24a3418dbc46995c3741f801a93c 100755
--- a/python/mesa_pd.py
+++ b/python/mesa_pd.py
@@ -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())
diff --git a/src/lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h b/src/lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h
index 4ca474cc0a24d59e9798bbfe4a387ac390a5da70..768b5a8eb280d4d4b6444b9182c22706d4b9bc08 100644
--- a/src/lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h
+++ b/src/lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h
@@ -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
{
diff --git a/src/lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h b/src/lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h
index 31c828e0d1c574d6b372048ff007ccaeefe931fc..e354a68b97c0cebc6aa479b964ecd7e565a08549 100644
--- a/src/lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h
+++ b/src/lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h
@@ -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
{
diff --git a/src/lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h b/src/lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h
index 07feb955e7f05f266db8a4cd89a562e6887a8406..4daa31e5a8d95be58f76b8e35f39c82f34044f99 100644
--- a/src/lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h
+++ b/src/lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h
@@ -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
{
diff --git a/src/lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h b/src/lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h
index 1e7dd22c42e296f9f98e694def5a457480764683..3cb84299928fc8b803f5e209af25191e23107abf 100644
--- a/src/lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h
+++ b/src/lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h
@@ -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
{
diff --git a/src/mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h b/src/mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h
new file mode 100644
index 0000000000000000000000000000000000000000..6155351601b2c0cd16673c6d865eb041152a77c4
--- /dev/null
+++ b/src/mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h
@@ -0,0 +1,121 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/mesa_pd/mpi/notifications/ParseMessage.h b/src/mesa_pd/mpi/notifications/ParseMessage.h
index 1794965c02c8e01f36be1149c41f61f9d55b7340..63aa0d542f4d969e6659e64507fb3119335c5473 100644
--- a/src/mesa_pd/mpi/notifications/ParseMessage.h
+++ b/src/mesa_pd/mpi/notifications/ParseMessage.h
@@ -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." );
diff --git a/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h b/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
index 1d9ca382cc3574306b470fea2cbe4b78b0c90cc8..2ac5383c17fdda618019c338ef36f47437da1cbf 100644
--- a/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
@@ -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;
}
diff --git a/src/mesa_pd/mpi/notifications/ParticleMigrationNotification.h b/src/mesa_pd/mpi/notifications/ParticleMigrationNotification.h
index 78246ac932580a03d5530c1a9e69f044f3d73662..836fa29ee8453dd1df578de69e6eccc2e00d2f31 100644
--- a/src/mesa_pd/mpi/notifications/ParticleMigrationNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleMigrationNotification.h
@@ -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_;
return buf;
}
diff --git a/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp b/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
index b3db226100e6bf087b210879bfff2af439e693d2..5aa9a51c2c8292815638f33b7e0916e9c6ae57e2 100644
--- a/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
+++ b/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
@@ -79,6 +79,7 @@
#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"
@@ -206,8 +207,8 @@ public:
{
pos = ac_->getPosition(idx);
transVel = ac_->getLinearVelocity(idx);
+ hydForce = ac_->getHydrodynamicForce(idx);
}
- hydForce = ac_->getHydrodynamicForce(idx);
}
WALBERLA_MPI_SECTION()
@@ -743,14 +744,17 @@ int main( int argc, char **argv )
timeloopTiming["RPD"].start();
+ // -> full hydrodynamic force/torque info is available on local particle
+ reduceProperty.operator()<mesa_pd::HydrodynamicForceTorqueNotification>(*ps);
+
if( averageForceTorqueOverTwoTimeSteps )
{
if( i == 0 )
{
lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueForAveragingKernel;
- ps->forEachParticle(useOpenMP, sphereSelector, *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
}
- ps->forEachParticle(useOpenMP, sphereSelector, *accessor, averageHydrodynamicForceTorque, *accessor );
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, averageHydrodynamicForceTorque, *accessor );
}
for(auto subCycle = uint_t(0); subCycle < numRPDSubCycles; ++subCycle )
@@ -762,7 +766,7 @@ int main( int argc, char **argv )
syncCall();
}
- ps->forEachParticle(useOpenMP, sphereSelector, *accessor, addHydrodynamicInteraction, *accessor );
+ ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addHydrodynamicInteraction, *accessor );
ps->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, addGravitationalForce, *accessor );
// lubrication correction