diff --git a/src/ladds/CollisionFunctor.cpp b/src/ladds/CollisionFunctor.cpp index bab04e74..d6d7e4ef 100644 --- a/src/ladds/CollisionFunctor.cpp +++ b/src/ladds/CollisionFunctor.cpp @@ -94,7 +94,30 @@ void CollisionFunctor::AoSFunctor(Particle &i, Particle &j, bool newton3) { } void CollisionFunctor::SoAFunctorSingle(autopas::SoAView 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(); + + // 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 soa1, @@ -118,6 +141,7 @@ void CollisionFunctor::SoAFunctorPair(autopas::SoAView 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 soa, } } +#pragma omp declare simd void CollisionFunctor::SoAKernel( size_t i, size_t j, autopas::SoAView &soa1, autopas::SoAView &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(); const auto *const __restrict y1ptr = soa1.template begin(); @@ -156,6 +176,13 @@ void CollisionFunctor::SoAKernel( const auto *const __restrict y2ptr = soa2.template begin(); const auto *const __restrict z2ptr = soa2.template begin(); + const auto *const __restrict vx1ptr = soa1.template begin(); + const auto *const __restrict vy1ptr = soa1.template begin(); + const auto *const __restrict vz1ptr = soa1.template begin(); + const auto *const __restrict vx2ptr = soa2.template begin(); + const auto *const __restrict vy2ptr = soa2.template begin(); + const auto *const __restrict vz2ptr = soa2.template begin(); + const auto *const __restrict ownedStatePtr2 = soa2.template begin(); const auto *const __restrict ptr1ptr = soa1.template begin(); @@ -164,7 +191,7 @@ void CollisionFunctor::SoAKernel( const auto *const __restrict id1ptr = soa1.template begin(); const auto *const __restrict id2ptr = soa2.template begin(); - 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); diff --git a/src/ladds/CollisionFunctor.h b/src/ladds/CollisionFunctor.h index e8eb56cb..59699ce9 100644 --- a/src/ladds/CollisionFunctor.h +++ b/src/ladds/CollisionFunctor.h @@ -31,16 +31,22 @@ class CollisionFunctor final : public autopas::Functor{Particle::AttributeNames::ptr, - Particle::AttributeNames::id, - Particle::AttributeNames::ownershipState, - Particle::AttributeNames::posX, - Particle::AttributeNames::posY, - Particle::AttributeNames::posZ}; + return std::array{ + 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) { diff --git a/src/ladds/Simulation.h b/src/ladds/Simulation.h index d8590eaa..7c219024 100644 --- a/src/ladds/Simulation.h +++ b/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 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. diff --git a/src/ladds/particle/Constellation.cpp b/src/ladds/particle/Constellation.cpp index d1130751..ee35407a 100644 --- a/src/ladds/particle/Constellation.cpp +++ b/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); diff --git a/src/ladds/particle/Constellation.h b/src/ladds/particle/Constellation.h index 4f57eab5..37ec4bb5 100644 --- a/src/ladds/particle/Constellation.h +++ b/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 diff --git a/src/ladds/particle/Particle.h b/src/ladds/particle/Particle.h index e2fde120..6576e2da 100644 --- a/src/ladds/particle/Particle.h +++ b/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) { diff --git a/tests/testladds/CollisionFunctorIntegrationTest.cpp b/tests/testladds/CollisionFunctorIntegrationTest.cpp index 443bf0ac..fc7b5b54 100644 --- a/tests/testladds/CollisionFunctorIntegrationTest.cpp +++ b/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 autopas; // allow all container options since the traversal determines it uniquely diff --git a/tests/testladds/CollisionFunctorTest.cpp b/tests/testladds/CollisionFunctorTest.cpp index 9a67060b..74357a34 100644 --- a/tests/testladds/CollisionFunctorTest.cpp +++ b/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); } }