[WIP] Functor Optimization

src/ladds/CollisionFunctor.cpp

@@ -94,7 +94,30 @@ void CollisionFunctor::AoSFunctor(Particle &i, Particle &j, bool newton3) {
 }
 
 void CollisionFunctor::SoAFunctorSingle(autopas::SoAView<SoAArraysType> soa, bool newton3) {
-  SoAFunctorPair(soa, soa, newton3);
+  const auto soaNumParticles = soa.getNumParticles();
+  if (soaNumParticles == 0) return;
+
+  // get pointers to the SoA
+  const auto *const __restrict ownedStatePtr1 = soa.template begin<Particle::AttributeNames::ownershipState>();
+
+  // outer loop over SoA1
+  for (size_t i = 0; i < soaNumParticles; ++i) {
+    if (ownedStatePtr1[i] == autopas::OwnershipState::dummy) {
+      // If the i-th particle is a dummy, skip this loop iteration.
+      continue;
+    }
+
+    // inner loop over SoA2
+    // custom reduction for collision collections (vector)
+#pragma omp declare reduction(vecMerge:CollisionCollectionT \
+                              : omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end()))
+    // alias because OpenMP needs it
+    auto &thisCollisions = _threadData[autopas::autopas_get_thread_num()].collisions;
+#pragma omp simd reduction(vecMerge : thisCollisions)
+    for (size_t j = i + 1; j < soaNumParticles; ++j) {
+      SoAKernel(i, j, soa, soa, newton3);
+    }
+  }
 }
 
 void CollisionFunctor::SoAFunctorPair(autopas::SoAView<SoAArraysType> soa1,
@@ -118,6 +141,7 @@ void CollisionFunctor::SoAFunctorPair(autopas::SoAView<SoAArraysType> soa1,
                               : omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end()))
     // alias because OpenMP needs it
     auto &thisCollisions = _threadData[autopas::autopas_get_thread_num()].collisions;
+    const auto soa2NumParticles = soa2.getNumParticles();
 #pragma omp simd reduction(vecMerge : thisCollisions)
     for (size_t j = 0; j < soa2.getNumParticles(); ++j) {
       SoAKernel(i, j, soa1, soa2, newton3);
@@ -141,13 +165,9 @@ void CollisionFunctor::SoAFunctorVerlet(autopas::SoAView<SoAArraysType> soa,
   }
 }
 
+#pragma omp declare simd
 void CollisionFunctor::SoAKernel(
     size_t i, size_t j, autopas::SoAView<SoAArraysType> &soa1, autopas::SoAView<SoAArraysType> &soa2, bool newton3) {
-  // TODO: as soon as this exception is removed / the SoAFunctor properly implemented
-  // remove the GTEST_SKIP in CollisionFunctorIntegrationTest::testAutoPasAlgorithm!
-  throw std::runtime_error(
-      "SoA kernel not up to date with AoS Kernel as it lacks the new linear interpolation distance.");
-
   // get pointers to the SoAs
   const auto *const __restrict x1ptr = soa1.template begin<Particle::AttributeNames::posX>();
   const auto *const __restrict y1ptr = soa1.template begin<Particle::AttributeNames::posY>();
@@ -156,6 +176,13 @@ void CollisionFunctor::SoAKernel(
   const auto *const __restrict y2ptr = soa2.template begin<Particle::AttributeNames::posY>();
   const auto *const __restrict z2ptr = soa2.template begin<Particle::AttributeNames::posZ>();
 
+  const auto *const __restrict vx1ptr = soa1.template begin<Particle::AttributeNames::velX>();
+  const auto *const __restrict vy1ptr = soa1.template begin<Particle::AttributeNames::velY>();
+  const auto *const __restrict vz1ptr = soa1.template begin<Particle::AttributeNames::velZ>();
+  const auto *const __restrict vx2ptr = soa2.template begin<Particle::AttributeNames::velX>();
+  const auto *const __restrict vy2ptr = soa2.template begin<Particle::AttributeNames::velY>();
+  const auto *const __restrict vz2ptr = soa2.template begin<Particle::AttributeNames::velZ>();
+
   const auto *const __restrict ownedStatePtr2 = soa2.template begin<Particle::AttributeNames::ownershipState>();
 
   const auto *const __restrict ptr1ptr = soa1.template begin<Particle::AttributeNames::ptr>();
@@ -164,7 +191,7 @@ void CollisionFunctor::SoAKernel(
   const auto *const __restrict id1ptr = soa1.template begin<Particle::AttributeNames::id>();
   const auto *const __restrict id2ptr = soa2.template begin<Particle::AttributeNames::id>();
 
-  if (ownedStatePtr2[j] == autopas::OwnershipState::dummy or id1ptr[i] == id2ptr[j]) {
+  if (ownedStatePtr2[j] == autopas::OwnershipState::dummy) {
     return;
   }
 
@@ -182,6 +209,52 @@ void CollisionFunctor::SoAKernel(
     return;
   }
 
+  // Get old time step position
+  const auto old_r_i_x = x1ptr[i] - vx1ptr[i] * _dt;
+  const auto old_r_i_y = y1ptr[i] - vy1ptr[i] * _dt;
+  const auto old_r_i_z = z1ptr[i] - vz1ptr[i] * _dt;
+  const auto old_r_j_x = x2ptr[j] - vx2ptr[j] * _dt;
+  const auto old_r_j_y = y2ptr[j] - vy2ptr[j] * _dt;
+  const auto old_r_j_z = z2ptr[j] - vz2ptr[j] * _dt;
+
+  // Compute nominator dot products
+  const auto vi_ri = vx1ptr[i] * old_r_i_x + vy1ptr[i] * old_r_i_y + vz1ptr[i] * old_r_i_z;
+  const auto vi_rj = vx1ptr[i] * old_r_j_x + vy1ptr[i] * old_r_j_y + vz1ptr[i] * old_r_j_z;
+  const auto vj_ri = vx2ptr[j] * old_r_i_x + vy2ptr[j] * old_r_i_y + vz2ptr[j] * old_r_i_z;
+  const auto vj_rj = vx2ptr[j] * old_r_j_x + vy2ptr[j] * old_r_j_y + vz2ptr[j] * old_r_j_z;
+
+  const auto nominator = vi_rj + vj_ri - vi_ri - vj_rj;
+
+  // Compute denominator dot products
+  const auto two_vi_vj = 2.0 * (vx1ptr[i] * vx2ptr[j] + vy1ptr[i] * vy2ptr[j] + vz1ptr[i] * vz2ptr[j]);
+  const auto vi_square = vx1ptr[i] * vx1ptr[i] + vy1ptr[i] * vy1ptr[i] + vz1ptr[i] * vz1ptr[i];
+  const auto vj_square = vx2ptr[j] * vx2ptr[j] + vy2ptr[j] * vy2ptr[j] + vz2ptr[j] * vz2ptr[j];
+
+  const auto denominator = vi_square + vj_square - two_vi_vj;
+
+  // Compute t at minimal distance
+  auto t = nominator / denominator;
+
+  // If in the past, minimum is at t = 0
+  // Else If in future timesteps, minimum for this is at t= _dt
+  t = std::clamp(t, 0., _dt);
+
+  // Compute actual distance by propagating along the line to t
+  const auto p1_x = old_r_i_x + vx1ptr[i] * t;
+  const auto p1_y = old_r_i_y + vy1ptr[i] * t;
+  const auto p1_z = old_r_i_z + vz1ptr[i] * t;
+  const auto p2_x = old_r_j_x + vx2ptr[j] * t;
+  const auto p2_y = old_r_j_y + vy2ptr[j] * t;
+  const auto p2_z = old_r_j_z + vz2ptr[j] * t;
+
+  const auto dr_lines_x = p1_x - p2_x;
+  const auto dr_lines_y = p1_y - p2_y;
+  const auto dr_lines_z = p1_z - p2_z;
+
+  const auto distanceSquare_lines = dr_lines_x * dr_lines_x + dr_lines_y * dr_lines_y + dr_lines_z * dr_lines_z;
+
+  if (distanceSquare_lines > _minorCutoffSquare) return;
+
   // store pointers to colliding pair
   if (id1ptr[i] < id2ptr[j]) {
     _threadData[autopas::autopas_get_thread_num()].collisions.emplace_back(ptr1ptr[i], ptr2ptr[j], dr2);

src/ladds/CollisionFunctor.h

@@ -31,16 +31,22 @@ class CollisionFunctor final : public autopas::Functor<Particle, CollisionFuncto
   }
 
   [[nodiscard]] bool allowsNonNewton3() final {
-    return true;
+    // as we are only interested in any interaction [i,j] it makes no sense to also look at [j,i]
+    return false;
   }
 
   [[nodiscard]] constexpr static auto getNeededAttr() {
-    return std::array<typename Particle::AttributeNames, 6>{Particle::AttributeNames::ptr,
-                                                            Particle::AttributeNames::id,
-                                                            Particle::AttributeNames::ownershipState,
-                                                            Particle::AttributeNames::posX,
-                                                            Particle::AttributeNames::posY,
-                                                            Particle::AttributeNames::posZ};
+    return std::array<typename Particle::AttributeNames, 9>{
+        Particle::AttributeNames::ptr,
+        Particle::AttributeNames::id,
+        Particle::AttributeNames::ownershipState,
+        Particle::AttributeNames::posX,
+        Particle::AttributeNames::posY,
+        Particle::AttributeNames::posZ,
+        Particle::AttributeNames::velX,
+        Particle::AttributeNames::velY,
+        Particle::AttributeNames::velZ,
+    };
   }
 
   [[nodiscard]] constexpr static auto getNeededAttr(std::false_type) {

src/ladds/Simulation.h

@@ -81,8 +81,8 @@ class Simulation {
    * @return Tuple of the collisions and whether AutoPas is currently in tuning mode.
    */
   std::tuple<CollisionFunctor::CollisionCollectionT, bool> collisionDetection(AutoPas_t &autopas,
-                                                            double deltaT,
-                                                            double conjunctionThreshold);
+                                                                              double deltaT,
+                                                                              double conjunctionThreshold);
 
   /**
    * Updates the configuration with the latest AutoPas configuration and writes it to a new YAML file.

src/ladds/particle/Constellation.cpp

@@ -15,7 +15,7 @@
 std::mt19937 Constellation::generator{42};
 
 Constellation::Constellation(const std::string &constellation_data_str, size_t interval, double altitudeDeviation)
-    : interval(interval), altitudeDeviation(altitudeDeviation), distribution(0., altitudeDeviation) {
+    : interval(interval), distribution(0., altitudeDeviation) {
   // split the 3 comma seperated arguments
   auto seperator1 = constellation_data_str.find(',', 0);
   auto seperator2 = constellation_data_str.find(',', seperator1 + 1);

src/ladds/particle/Constellation.h

@@ -141,12 +141,6 @@ class Constellation {
    */
   int planesDeployed = 0;
 
-  /**
-   * deviation parameter of the normal distribution that determines the deviation
-   * of the satellites base altitude
-   */
-  double altitudeDeviation;
-
   /**
    * seeded/deterministic random number generator used to add noise to the
    * altitudes of satellites

src/ladds/particle/Particle.h

@@ -24,7 +24,7 @@ class Particle final : public autopas::ParticleFP64 {
   /**
    * Enums used as ids for accessing and creating a dynamically sized SoA.
    */
-  enum AttributeNames : int { ptr, id, posX, posY, posZ, forceX, forceY, forceZ, ownershipState };
+  enum AttributeNames : int { ptr, id, posX, posY, posZ, velX, velY, velZ, forceX, forceY, forceZ, ownershipState };
 
   /**
    * The type for the SoA storage.
@@ -34,6 +34,9 @@ class Particle final : public autopas::ParticleFP64 {
                                                          double /*x*/,
                                                          double /*y*/,
                                                          double /*z*/,
+                                                         double /*vx*/,
+                                                         double /*vy*/,
+                                                         double /*vz*/,
                                                          double /*fx*/,
                                                          double /*fy*/,
                                                          double /*fz*/,
@@ -56,6 +59,12 @@ class Particle final : public autopas::ParticleFP64 {
       return getR()[1];
     } else if constexpr (attribute == AttributeNames::posZ) {
       return getR()[2];
+    } else if constexpr (attribute == AttributeNames::velX) {
+      return getV()[0];
+    } else if constexpr (attribute == AttributeNames::velY) {
+      return getV()[1];
+    } else if constexpr (attribute == AttributeNames::velZ) {
+      return getV()[2];
     } else if constexpr (attribute == AttributeNames::forceX) {
       return getF()[0];
     } else if constexpr (attribute == AttributeNames::forceY) {

tests/testladds/CollisionFunctorIntegrationTest.cpp

@@ -71,13 +71,12 @@ TEST_P(CollisionFunctorIntegrationTest, testAutoPasAlgorithm) {
 
   const auto &[traversal, dataLayout, newton3, cellSizeFactor] = GetParam();
 
-  // This is currently necessary until we implement the SoA functor
-  if (dataLayout == autopas::DataLayoutOption::soa) {
-    GTEST_SKIP_("SoAFunctor currently not implemented!");
-  }
-
   CollisionFunctor functor(_cutoff, 10.0, 0.1 * _cutoff);
 
+  if (not functor.allowsNonNewton3() and newton3 == autopas::Newton3Option::disabled) {
+    GTEST_SKIP_("Functor does not support Newton3==disabled!");
+  }
+
   // configure the AutoPas container
   autopas::AutoPas<Particle> autopas;
   // allow all container options since the traversal determines it uniquely

tests/testladds/CollisionFunctorTest.cpp

@@ -32,12 +32,9 @@ TEST(CollisionFunctorTest, ThreeParticles) {
 
   CollisionFunctor collisionFunctor(cutoff, 10.0, cutoff);
 
-  for (auto &di : debris) {
-    for (auto &dj : debris) {
-      if (di == dj) {
-        continue;
-      }
-      collisionFunctor.AoSFunctor(di, dj, newton3);
+  for (size_t i = 0; i < debris.size(); ++i) {
+    for (size_t j = i + 1; j < debris.size(); ++j) {
+      collisionFunctor.AoSFunctor(debris[i], debris[j], newton3);
     }
   }
 
@@ -70,12 +67,9 @@ TEST_P(CollisionFunctorTest, LinearInterpolationTest) {
   debris.emplace_back(x1, v1, 0.);
   debris.emplace_back(x2, v2, 1.);
 
-  for (auto &di : debris) {
-    for (auto &dj : debris) {
-      if (di == dj) {
-        continue;
-      }
-      collisionFunctor.AoSFunctor(di, dj, newton3);
+  for (size_t i = 0; i < debris.size(); ++i) {
+    for (size_t j = i + 1; j < debris.size(); ++j) {
+      collisionFunctor.AoSFunctor(debris[i], debris[j], newton3);
     }
   }