Merge branch 'check_interaction_radius' into 'master'

added asserts to ensure interactionRadius is set properly

See merge request walberla/walberla!439

apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp

@@ -309,6 +309,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
    // create bounding planes
    mesa_pd::data::Particle p0 = *ps->create(true);
    p0.setPosition(simulationDomain.minCorner());
+   p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
    p0.setOwner(mpi::MPIManager::instance()->rank());
    p0.setType(0);
@@ -317,6 +318,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
 
    mesa_pd::data::Particle p1 = *ps->create(true);
    p1.setPosition(simulationDomain.maxCorner());
+   p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
    p1.setOwner(mpi::MPIManager::instance()->rank());
    p1.setType(0);

apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp

@@ -443,6 +443,7 @@ int main( int argc, char **argv )
       // create two planes
       mesa_pd::data::Particle&& p0 = *ps->create(true);
       p0.setPosition(referenceVector);
+      p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
       p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
       p0.setOwner(mpi::MPIManager::instance()->rank());
       mesa_pd::data::particle_flags::set(p0.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
@@ -451,6 +452,7 @@ int main( int argc, char **argv )
 
       mesa_pd::data::Particle&& p1 = *ps->create(true);
       p1.setPosition(Vector3<real_t>(real_c(xSize),0,0));
+      p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
       p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
       p1.setOwner(mpi::MPIManager::instance()->rank());
       mesa_pd::data::particle_flags::set(p1.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);

apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp

@@ -119,11 +119,13 @@ int main( int argc, char ** argv )
    data::Particle&& p = *ps->create();
    p.setPosition(Vec3(0,0,2*radius));
    p.setLinearVelocity(Vec3(uTin, 0., -uNin));
+   p.setInteractionRadius(radius);
    p.setType(0);
 
    // create plane
    data::Particle&& p0 = *ps->create(true);
    p0.setPosition(Vec3(0,0,0));
+   p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p0.setShapeID(ss->create<data::HalfSpace>(Vector3<real_t>(0,0,1)));
    p0.setType(0);
    data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::INFINITE);

apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp

@@ -370,6 +370,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
    // create bounding planes
    mesa_pd::data::Particle&& p0 = *ps->create(true);
    p0.setPosition(simulationDomain.minCorner());
+   p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
    p0.setOwner(mpi::MPIManager::instance()->rank());
    p0.setType(0);
@@ -381,6 +382,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
       //only create top plane when no outflow BC should be set there
       mesa_pd::data::Particle&& p1 = *ps->create(true);
       p1.setPosition(simulationDomain.maxCorner());
+      p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
       p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
       p1.setOwner(mpi::MPIManager::instance()->rank());
       p1.setType(0);

apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp

@@ -290,6 +290,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
    // create bounding planes
    mesa_pd::data::Particle p0 = *ps->create(true);
    p0.setPosition(simulationDomain.minCorner());
+   p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
    p0.setOwner(mpi::MPIManager::instance()->rank());
    p0.setType(0);
@@ -298,6 +299,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
 
    mesa_pd::data::Particle p1 = *ps->create(true);
    p1.setPosition(simulationDomain.maxCorner());
+   p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
    p1.setOwner(mpi::MPIManager::instance()->rank());
    p1.setType(0);
@@ -306,6 +308,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
 
    mesa_pd::data::Particle p2 = *ps->create(true);
    p2.setPosition(simulationDomain.minCorner());
+   p2.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p2.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
    p2.setOwner(mpi::MPIManager::instance()->rank());
    p2.setType(0);
@@ -314,6 +317,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
 
    mesa_pd::data::Particle p3 = *ps->create(true);
    p3.setPosition(simulationDomain.maxCorner());
+   p3.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p3.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
    p3.setOwner(mpi::MPIManager::instance()->rank());
    p3.setType(0);
@@ -322,6 +326,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
 
    mesa_pd::data::Particle p4 = *ps->create(true);
    p4.setPosition(simulationDomain.minCorner());
+   p4.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p4.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,1,0) ));
    p4.setOwner(mpi::MPIManager::instance()->rank());
    p4.setType(0);
@@ -330,6 +335,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
 
    mesa_pd::data::Particle p5 = *ps->create(true);
    p5.setPosition(simulationDomain.maxCorner());
+   p5.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p5.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,-1,0) ));
    p5.setOwner(mpi::MPIManager::instance()->rank());
    p5.setType(0);

apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp

@@ -288,6 +288,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
 {
    mesa_pd::data::Particle p2 = *ps->create(true);
    p2.setPosition(simulationDomain.minCorner());
+   p2.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p2.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
    p2.setOwner(mpi::MPIManager::instance()->rank());
    p2.setType(0);
@@ -297,6 +298,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
 
    mesa_pd::data::Particle p3 = *ps->create(true);
    p3.setPosition(simulationDomain.maxCorner());
+   p3.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p3.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
    p3.setOwner(mpi::MPIManager::instance()->rank());
    p3.setType(0);
@@ -306,6 +308,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
 
    mesa_pd::data::Particle p4 = *ps->create(true);
    p4.setPosition(simulationDomain.minCorner());
+   p4.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p4.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,1,0) ));
    p4.setOwner(mpi::MPIManager::instance()->rank());
    p4.setType(0);
@@ -315,6 +318,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
 
    mesa_pd::data::Particle p5 = *ps->create(true);
    p5.setPosition(simulationDomain.maxCorner());
+   p5.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p5.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,-1,0) ));
    p5.setOwner(mpi::MPIManager::instance()->rank());
    p5.setType(0);

apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp

@@ -367,6 +367,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
    // create bounding planes
    mesa_pd::data::Particle&& p0 = *ps->create(true);
    p0.setPosition(simulationDomain.minCorner());
+   p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
    p0.setOwner(mpi::MPIManager::instance()->rank());
    p0.setType(0);
@@ -378,6 +379,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
       //only create top plane when no outflow BC should be set there
       mesa_pd::data::Particle&& p1 = *ps->create(true);
       p1.setPosition(simulationDomain.maxCorner());
+      p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
       p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
       p1.setOwner(mpi::MPIManager::instance()->rank());
       p1.setType(0);

apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadDistribution.cpp

@@ -389,6 +389,7 @@ void createPlane( const shared_ptr<mesa_pd::data::ParticleStorage> & ps, const s
 {
    mesa_pd::data::Particle&& p0 = *ps->create(true);
    p0.setPosition(position);
+   p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( normal.getNormalized() ));
    p0.setOwner(mpi::MPIManager::instance()->rank());
    p0.setType(0);

apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadEvaluation.cpp

@@ -460,6 +460,7 @@ int main( int argc, char **argv )
 
    mesa_pd::data::Particle&& p0 = *ps->create(true);
    p0.setPosition(generationDomain.minCorner());
+   p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
    p0.setOwner(mpi::MPIManager::instance()->rank());
    p0.setType(0);
@@ -468,6 +469,7 @@ int main( int argc, char **argv )
 
    mesa_pd::data::Particle&& p1 = *ps->create(true);
    p1.setPosition(generationDomain.maxCorner());
+   p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
    p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
    p1.setOwner(mpi::MPIManager::instance()->rank());
    p1.setType(0);

apps/benchmarks/GranularGas/CreateParticles.cpp

@@ -30,6 +30,7 @@ data::ParticleStorage::iterator createPlane( data::ParticleStorage& ps,
 {
    auto p0              = ps.create(true);
    p0->getPositionRef() = pos;
+   p0->getInteractionRadiusRef() = std::numeric_limits<real_t>::infinity();
    p0->getShapeIDRef()  = ss.create<data::HalfSpace>( normal );
    p0->getOwnerRef()    = walberla::mpi::MPIManager::instance()->rank();
    p0->getTypeRef()     = 0;

python/mesa_pd/templates/kernel/InsertParticleIntoLinkedCells.templ.h

@@ -77,7 +77,9 @@ inline void InsertParticleIntoLinkedCells::operator()(const size_t p_idx, Access
       ac.setNextParticle(p_idx, lc.infiniteParticles_.exchange(int_c(p_idx)));
    } else
    {
+      WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
       {%- for dim in range(3) %}
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
       int hash{{dim}} = static_cast<int>(std::floor((ac.getPosition(p_idx)[{{dim}}] - minCorner[{{dim}}]) * lc.invCellDiameter_[{{dim}}]));
       {%- endfor %}
       {%- for dim in range(3) %}

python/mesa_pd/templates/kernel/InsertParticleIntoSparseLinkedCells.templ.h

@@ -77,7 +77,9 @@ inline void InsertParticleIntoSparseLinkedCells::operator()(const size_t p_idx,
       ac.setNextParticle(p_idx, lc.infiniteParticles_.exchange(int_c(p_idx)));
    } else
    {
+      WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
       {%- for dim in range(3) %}
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
       int hash{{dim}} = static_cast<int>(std::floor((ac.getPosition(p_idx)[{{dim}}] - minCorner[{{dim}}]) * lc.invCellDiameter_[{{dim}}]));
       {%- endfor %}
       {%- for dim in range(3) %}

python/mesa_pd/templates/mpi/SyncGhostOwners.templ.cpp

@@ -72,6 +72,7 @@ void SyncGhostOwners::updateAndMigrate( data::ParticleStorage& ps,
    std::set<walberla::mpi::MPIRank> recvRanks; // potential message senders
    for( auto pIt = ps.begin(); pIt != ps.end(); ++pIt)
    {
+      WALBERLA_ASSERT_GREATER(pIt->getInteractionRadius(), 0_r, "Did you forget to set the interaction radius?");
       if (isSet( pIt->getFlags(), GHOST))
       {
          if (!isSet( pIt->getFlags(), NON_COMMUNICATING) || syncNonCommunicatingBodies)

python/mesa_pd/templates/mpi/SyncNextNeighbors.templ.cpp

@@ -87,6 +87,8 @@ void SyncNextNeighbors::generateSynchronizationMessages(walberla::mpi::BufferSys
    // position update
    for( auto pIt = ps.begin(); pIt != ps.end(); )
    {
+      WALBERLA_ASSERT_GREATER(pIt->getInteractionRadius(), 0_r, "Did you forget to set the interaction radius?");
+
       //skip all ghost particles
       if (data::particle_flags::isSet( pIt->getFlags(), data::particle_flags::GHOST))
       {

src/mesa_pd/common/AABBConversion.h

@@ -32,6 +32,7 @@ namespace mesa_pd {
 
 math::AABB getAABBFromInteractionRadius(const Vector3<real_t> & pos, const real_t interactionRadius )
 {
+   WALBERLA_ASSERT_GREATER(interactionRadius, 0_r, "Did you forget to set the interaction radius?");
    return math::AABB( pos[0]-interactionRadius, pos[1]-interactionRadius, pos[2]-interactionRadius,
                       pos[0]+interactionRadius, pos[1]+interactionRadius, pos[2]+interactionRadius );
 }

src/mesa_pd/kernel/InsertParticleIntoLinkedCells.h

@@ -72,8 +72,12 @@ inline void InsertParticleIntoLinkedCells::operator()(const size_t p_idx, Access
       ac.setNextParticle(p_idx, lc.infiniteParticles_.exchange(int_c(p_idx)));
    } else
    {
+      WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
       int hash0 = static_cast<int>(std::floor((ac.getPosition(p_idx)[0] - minCorner[0]) * lc.invCellDiameter_[0]));
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
       int hash1 = static_cast<int>(std::floor((ac.getPosition(p_idx)[1] - minCorner[1]) * lc.invCellDiameter_[1]));
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
       int hash2 = static_cast<int>(std::floor((ac.getPosition(p_idx)[2] - minCorner[2]) * lc.invCellDiameter_[2]));
       if (hash0 < 0) hash0 = 0;
       if (hash0 >= lc.numCellsPerDim_[0]) hash0 = lc.numCellsPerDim_[0] - 1;

src/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h

@@ -72,8 +72,12 @@ inline void InsertParticleIntoSparseLinkedCells::operator()(const size_t p_idx,
       ac.setNextParticle(p_idx, lc.infiniteParticles_.exchange(int_c(p_idx)));
    } else
    {
+      WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
       int hash0 = static_cast<int>(std::floor((ac.getPosition(p_idx)[0] - minCorner[0]) * lc.invCellDiameter_[0]));
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
       int hash1 = static_cast<int>(std::floor((ac.getPosition(p_idx)[1] - minCorner[1]) * lc.invCellDiameter_[1]));
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
       int hash2 = static_cast<int>(std::floor((ac.getPosition(p_idx)[2] - minCorner[2]) * lc.invCellDiameter_[2]));
       if (hash0 < 0) hash0 = 0;
       if (hash0 >= lc.numCellsPerDim_[0]) hash0 = lc.numCellsPerDim_[0] - 1;

src/mesa_pd/mpi/SyncGhostOwners.cpp

@@ -72,6 +72,7 @@ void SyncGhostOwners::updateAndMigrate( data::ParticleStorage& ps,
    std::set<walberla::mpi::MPIRank> recvRanks; // potential message senders
    for( auto pIt = ps.begin(); pIt != ps.end(); ++pIt)
    {
+      WALBERLA_ASSERT_GREATER(pIt->getInteractionRadius(), 0_r, "Did you forget to set the interaction radius?");
       if (isSet( pIt->getFlags(), GHOST))
       {
          if (!isSet( pIt->getFlags(), NON_COMMUNICATING) || syncNonCommunicatingBodies)

src/mesa_pd/mpi/SyncNextNeighbors.cpp

@@ -87,6 +87,8 @@ void SyncNextNeighbors::generateSynchronizationMessages(walberla::mpi::BufferSys
    // position update
    for( auto pIt = ps.begin(); pIt != ps.end(); )
    {
+      WALBERLA_ASSERT_GREATER(pIt->getInteractionRadius(), 0_r, "Did you forget to set the interaction radius?");
+
       //skip all ghost particles
       if (data::particle_flags::isSet( pIt->getFlags(), data::particle_flags::GHOST))
       {

src/mesa_pd/mpi/SyncNextNeighborsBlockForest.cpp

@@ -111,6 +111,8 @@ void SyncNextNeighborsBlockForest::generateSynchronizationMessages(walberla::mpi
    // position update
    for( auto pIt = ps.begin(); pIt != ps.end(); )
    {
+      WALBERLA_ASSERT_GREATER(pIt->getInteractionRadius(), 0_r, "Did you forget to set the interaction radius?");
+
       //skip all ghost particles
       if (data::particle_flags::isSet( pIt->getFlags(), data::particle_flags::GHOST))
       {