diff --git a/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp b/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
index 77524bf2da21e491336e7ada054fb3ac76483747..93da4ad7e225d7523530e4e5decf3ee014480fa0 100644
--- a/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
@@ -55,6 +55,7 @@
#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
#include "lbm_mesapd_coupling/DataTypes.h"
#include "lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
+#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/OmegaBulkAdaption.h"
@@ -714,6 +715,11 @@ int main( int argc, char **argv )
timeloop.singleStep( timeloopTiming );
// 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::SelectAll(), *accessor, averageHydrodynamicForceTorque, *accessor );
// evaluation
diff --git a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
index 12ef798b27eb670f4b32135f656be9217999a197..0e1aa9653c7687c7087f95479a721b7c76999880 100644
--- a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
@@ -65,6 +65,7 @@
#include "lbm_mesapd_coupling/DataTypes.h"
#include "lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h"
+#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/LubricationCorrectionKernel.h"
#include "lbm_mesapd_coupling/utility/OmegaBulkAdaption.h"
@@ -1057,9 +1058,13 @@ int main( int argc, char **argv )
// perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions
timeloop.singleStep( timeloopTiming );
-
- if( averageForceTorqueOverTwoTimeSteps && i!= 0)
+ 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::SelectAll(), *accessor, averageHydrodynamicForceTorque, *accessor );
}
diff --git a/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp b/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
index b8a2b2c1929518665c6dc6f84fbe2b6d78558afc..d0d186c16092a948f047f1ac6844af2a0052e3c8 100644
--- a/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
@@ -58,6 +58,7 @@
#include "lbm_mesapd_coupling/DataTypes.h"
#include "lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h"
+#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/LubricationCorrectionKernel.h"
#include "lbm_mesapd_coupling/utility/OmegaBulkAdaption.h"
@@ -372,7 +373,7 @@ int main( int argc, char **argv )
//numerical parameters
uint_t numberOfCellsInHorizontalDirection = uint_t(135);
- bool averageForceTorqueOverTwoTimSteps = true;
+ bool averageForceTorqueOverTwoTimeSteps = true;
bool conserveMomentum = false;
uint_t numRPDSubCycles = uint_t(1);
bool useVelocityVerlet = true;
@@ -393,7 +394,7 @@ int main( int argc, char **argv )
if( std::strcmp( argv[i], "--fluidType" ) == 0 ) { fluidType = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--numRPDSubCycles" ) == 0 ) { numRPDSubCycles = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--resolution" ) == 0 ) { numberOfCellsInHorizontalDirection = uint_c( std::atof( argv[++i] ) ); continue; }
- if( std::strcmp( argv[i], "--noForceAveraging" ) == 0 ) { averageForceTorqueOverTwoTimSteps = false; continue; }
+ if( std::strcmp( argv[i], "--noForceAveraging" ) == 0 ) { averageForceTorqueOverTwoTimeSteps = false; continue; }
if( std::strcmp( argv[i], "--conserveMomentum" ) == 0 ) { conserveMomentum = true; continue; }
if( std::strcmp( argv[i], "--baseFolder" ) == 0 ) { baseFolder = argv[++i]; continue; }
if( std::strcmp( argv[i], "--useEulerIntegrator" ) == 0 ) { useVelocityVerlet = false; continue; }
@@ -811,8 +812,13 @@ int main( int argc, char **argv )
timeloopTiming["RPD"].start();
- if( averageForceTorqueOverTwoTimSteps && i!= 0)
+ 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::SelectAll(), *accessor, averageHydrodynamicForceTorque, *accessor );
}
diff --git a/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp b/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
index 644a887df3aeb44096f2f267934bb9729035a89d..eac6530d09e73028be048bc494c7cda60b8bd1b7 100644
--- a/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
@@ -63,6 +63,7 @@
#include "lbm_mesapd_coupling/DataTypes.h"
#include "lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h"
+#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/OmegaBulkAdaption.h"
@@ -393,7 +394,7 @@ int main( int argc, char **argv )
bool vtkOutputAtEnd = true;
//numerical parameters
- bool averageForceTorqueOverTwoTimSteps = true;
+ bool averageForceTorqueOverTwoTimeSteps = true;
bool conserveMomentum = false;
uint_t numRPDSubCycles = uint_t(1);
std::string reconstructorType = "Grad"; // Eq, EAN, Ext, Grad
@@ -424,7 +425,7 @@ int main( int argc, char **argv )
if( std::strcmp( argv[i], "--logDensity" ) == 0 ) { logDensity = true; continue; }
if( std::strcmp( argv[i], "--vtkIOFreq" ) == 0 ) { vtkIOFreq = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--numRPDSubCycles" ) == 0 ) { numRPDSubCycles = uint_c( std::atof( argv[++i] ) ); continue; }
- if( std::strcmp( argv[i], "--noForceAveraging" ) == 0 ) { averageForceTorqueOverTwoTimSteps = false; continue; }
+ if( std::strcmp( argv[i], "--noForceAveraging" ) == 0 ) { averageForceTorqueOverTwoTimeSteps = false; continue; }
if( std::strcmp( argv[i], "--conserveMomentum" ) == 0 ) { conserveMomentum = true; continue; }
if( std::strcmp( argv[i], "--baseFolder" ) == 0 ) { baseFolder = argv[++i]; continue; }
if( std::strcmp( argv[i], "--reconstructorType" ) == 0 ) { reconstructorType = argv[++i]; continue; }
@@ -768,8 +769,13 @@ int main( int argc, char **argv )
timeloopTiming["RPD"].start();
- if( averageForceTorqueOverTwoTimSteps && i!= 0)
+ 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::SelectAll(), *accessor, averageHydrodynamicForceTorque, *accessor );
}
diff --git a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
index 0890c4ef497d9b840374a0cb4fe75bea2d4ff7c7..594cfd4325af59244525e1b570c6c4aa9891e7bb 100644
--- a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
@@ -68,6 +68,7 @@
#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/LubricationCorrectionKernel.h"
#include "lbm_mesapd_coupling/utility/OmegaBulkAdaption.h"
+#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
#include "mesa_pd/collision_detection/AnalyticContactDetection.h"
#include "mesa_pd/data/ParticleAccessorWithShape.h"
@@ -1183,9 +1184,13 @@ int main( int argc, char **argv )
// perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions
timeloop.singleStep( timeloopTiming );
-
- if( averageForceTorqueOverTwoTimeSteps && i!= 0)
+ 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::SelectAll(), *accessor, averageHydrodynamicForceTorque, *accessor );
}
diff --git a/src/lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h b/src/lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h
new file mode 100644
index 0000000000000000000000000000000000000000..07feb955e7f05f266db8a4cd89a562e6887a8406
--- /dev/null
+++ b/src/lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h
@@ -0,0 +1,54 @@
+//======================================================================================================================
+//
+// 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 InitializeHydrodynamicForceTorqueForAveragingKernel.h
+//! \ingroup lbm_mesapd_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "mesa_pd/data/IAccessor.h"
+
+namespace walberla {
+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
+{
+
+public:
+
+ InitializeHydrodynamicForceTorqueForAveragingKernel( ) = default;
+
+ template< typename ParticleAccessor_T >
+ void operator()(const size_t idx, ParticleAccessor_T& ac) const
+ {
+ static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
+
+ ac.setOldHydrodynamicForce( idx, ac.getHydrodynamicForce(idx) );
+ ac.setOldHydrodynamicTorque( idx, ac.getHydrodynamicTorque(idx) );
+ }
+};
+
+} // namespace lbm_mesapd_coupling
+} // namespace walberla
diff --git a/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp b/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
index ef8db99fa2e19a58c8c14ac96973722e18cff824..b3db226100e6bf087b210879bfff2af439e693d2 100644
--- a/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
+++ b/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
@@ -55,6 +55,7 @@
#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/PdfReconstructionManager.h"
#include "lbm_mesapd_coupling/utility/AddForceOnParticlesKernel.h"
#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
+#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
#include "lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
#include "lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h"
#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
@@ -361,7 +362,7 @@ int main( int argc, char **argv )
//numerical parameters
uint_t numberOfCellsInHorizontalDirection = uint_t(135);
- bool averageForceTorqueOverTwoTimSteps = true;
+ bool averageForceTorqueOverTwoTimeSteps = true;
bool conserveMomentum = false;
uint_t numRPDSubCycles = uint_t(1);
bool useVelocityVerlet = false;
@@ -377,7 +378,7 @@ int main( int argc, char **argv )
if( std::strcmp( argv[i], "--fluidType" ) == 0 ) { fluidType = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--numRPDSubCycles" ) == 0 ) { numRPDSubCycles = uint_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--resolution" ) == 0 ) { numberOfCellsInHorizontalDirection = uint_c( std::atof( argv[++i] ) ); continue; }
- if( std::strcmp( argv[i], "--noForceAveraging" ) == 0 ) { averageForceTorqueOverTwoTimSteps = false; continue; }
+ if( std::strcmp( argv[i], "--noForceAveraging" ) == 0 ) { averageForceTorqueOverTwoTimeSteps = false; continue; }
if( std::strcmp( argv[i], "--baseFolder" ) == 0 ) { baseFolder = argv[++i]; continue; }
if( std::strcmp( argv[i], "--useVV" ) == 0 ) { useVelocityVerlet = true; continue; }
if( std::strcmp( argv[i], "--useEANReconstructor" ) == 0 ) { useEANReconstructor = true; continue; }
@@ -742,8 +743,13 @@ int main( int argc, char **argv )
timeloopTiming["RPD"].start();
- if( averageForceTorqueOverTwoTimSteps && i!= 0)
+ if( averageForceTorqueOverTwoTimeSteps )
{
+ if( i == 0 )
+ {
+ lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel initializeHydrodynamicForceTorqueForAveragingKernel;
+ ps->forEachParticle(useOpenMP, sphereSelector, *accessor, initializeHydrodynamicForceTorqueForAveragingKernel, *accessor );
+ }
ps->forEachParticle(useOpenMP, sphereSelector, *accessor, averageHydrodynamicForceTorque, *accessor );
}