Commit 0c6c9b84 authored by Christoph Rettinger's avatar Christoph Rettinger

Merge branch 'mr_coupling_updates' into 'master'

Used new functionalities and updated coupling

See merge request !254
parents 5c808551 8bb2c3c4
Pipeline #22396 passed with stages
in 498 minutes and 17 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();
......
......@@ -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):
......
......@@ -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)
......
......@@ -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::