From 388478352e36d9eae3a1eb09227cfd97dcd118fd Mon Sep 17 00:00:00 2001 From: Kristjan Kongas Date: Sat, 12 Oct 2019 19:23:10 +0300 Subject: [PATCH] Avoid checking collisions and computing forces twice between pairs of particles --- Lib/Integrators.h | 12 +++++++----- Lib/Universe.cpp | 38 ++++++++++++++++++++++++++------------ Lib/Universe.h | 11 ++++++----- Tests/IntegratorTest.cpp | 2 +- 4 files changed, 40 insertions(+), 23 deletions(-) diff --git a/Lib/Integrators.h b/Lib/Integrators.h index 19919f4..89e52cb 100644 --- a/Lib/Integrators.h +++ b/Lib/Integrators.h @@ -29,7 +29,8 @@ void advanceEuler(IntegrableState &x, const Differetiator &diff, double dT) { diff.prepareDifferentiation(x); static IntegrableState k1; - diff.derivative(k1, x); + static IntegrableState derivativeCache; // Only for performance reasons + diff.derivative(k1, derivativeCache, x); k1 *= dT; x += k1; } @@ -49,11 +50,12 @@ void advanceRungeKutta4(IntegrableState &x, const Differetiator &diff, double dT static IntegrableState xInitial, xAdditive; static IntegrableState k1, k2, k3, k4; + static IntegrableState derivativeCache; // Only for performance reasons xInitial = x; // Compute k1 and add to x - diff.derivative(k1, x); + diff.derivative(k1, derivativeCache, x); k1 *= dT; xAdditive = k1; @@ -63,7 +65,7 @@ void advanceRungeKutta4(IntegrableState &x, const Differetiator &diff, double dT // Compute k2 and add to x k1 *= .5; k1 += xInitial; - diff.derivative(k2, k1); + diff.derivative(k2, derivativeCache, k1); k2 *= dT; xAdditive = k2; xAdditive *= 2. / 6.; @@ -72,7 +74,7 @@ void advanceRungeKutta4(IntegrableState &x, const Differetiator &diff, double dT // Compute k3 and add to x k2 *= .5; k2 += xInitial; - diff.derivative(k3, k2); + diff.derivative(k3, derivativeCache, k2); k3 *= dT; xAdditive = k3; xAdditive *= 2. / 6.; @@ -80,7 +82,7 @@ void advanceRungeKutta4(IntegrableState &x, const Differetiator &diff, double dT // Compute k4 and add to x k3 += xInitial; - diff.derivative(k4, k3); + diff.derivative(k4, derivativeCache, k3); k4 *= dT; xAdditive = k4; xAdditive *= 1. / 6.; diff --git a/Lib/Universe.cpp b/Lib/Universe.cpp index 6172230..8ff2d40 100644 --- a/Lib/Universe.cpp +++ b/Lib/Universe.cpp @@ -134,10 +134,13 @@ void UniverseDifferentiator::prepareDifferentiation(UniverseState &state) const state.prepareDifferentiation(); } -void UniverseDifferentiator::derivative(UniverseState &der, UniverseState &state) const { +void UniverseDifferentiator::derivative(UniverseState &der, UniverseState &derivativeCache, UniverseState &state) const { + // Using derivative cache as another accumulator for forces to avoid data race + initForces(der, state); - computeForces(der, state); - forcesToAccel(der); + initForces(derivativeCache, state); + computeForces(der, derivativeCache, state); + forcesToAccel(der, derivativeCache); } void UniverseDifferentiator::initForces(UniverseState &der, const UniverseState &state) const { @@ -154,17 +157,20 @@ void UniverseDifferentiator::initForces(UniverseState &der, const UniverseState } } -void UniverseDifferentiator::computeForces(UniverseState &der, const UniverseState &state) const { +void UniverseDifferentiator::computeForces(UniverseState &der, UniverseState &derivativeCache, const UniverseState &state) const { assert(! state.state.empty()); size_t total = state.state.size() * state.state[0].size(); - cv::parallel_for_(cv::Range{ 0, (int) total }, UniverseDifferentiator::ParallelForces{ *this, der, state }); + cv::parallel_for_(cv::Range{ 0, (int) total }, + UniverseDifferentiator::ParallelForces{ *this, der, derivativeCache, state }); } -void UniverseDifferentiator::forcesToAccel(UniverseState &der) const { +void UniverseDifferentiator::forcesToAccel(UniverseState &der, const UniverseState &derivativeCache) const { for(size_t y = 0; y < der.state.size(); ++y) for(size_t x = 0; x < der.state[y].size(); ++x) for(size_t i = 0; i < der.state[y][x].size(); ++i) { auto &pDer = der.state[y][x][i]; + auto &pDerCache = derivativeCache.state[y][x][i]; + pDer.v += pDerCache.v; pDer.v *= 1. / pDer.type->getMass(); } } @@ -175,8 +181,9 @@ double UniverseDifferentiator::boundForce(double overEdge) const { } -UniverseDifferentiator::ParallelForces::ParallelForces(const UniverseDifferentiator &diff, UniverseState &der, const UniverseState &state) - : diff(diff), der(der), state(state) { +UniverseDifferentiator::ParallelForces::ParallelForces(const UniverseDifferentiator &diff, + UniverseState &der, UniverseState &derivativeCache, const UniverseState &state) + : diff(diff), der(der), derivativeCache(derivativeCache), state(state) { } void UniverseDifferentiator::ParallelForces::operator()(const cv::Range& range) const { @@ -189,15 +196,22 @@ void UniverseDifferentiator::ParallelForces::operator()(const cv::Range& range) for(int i0 = 0; i0 < (int) state.state[y0][x0].size(); ++i0) { const auto &pState0 = state.state[y0][x0][i0]; auto &pDer0 = der.state[y0][x0][i0]; - for(int y1 = std::max(0, y0 - 1); y1 < std::min((int) state.state.size(), y0 + 2); ++y1) - for(int x1 = std::max(0, x0 - 1); x1 < std::min((int) state.state[y1].size(), x0 + 2); ++x1) - for(int i1 = 0; i1 < (int) state.state[y1][x1].size(); ++i1) { - if(y0 == y1 && x0 == x1 && i0 == i1) continue; + // The if conditions avoid computing the same force twice + for(int y1 = std::max(0, y0 - 1); y1 < std::min((int) state.state.size(), y0 + 2); ++y1) { + if (y1 < y0) continue; + for (int x1 = std::max(0, x0 - 1); x1 < std::min((int) state.state[y1].size(), x0 + 2); ++x1) { + if(y1 == y0 && x1 < x0) continue; + for (int i1 = 0; i1 < (int) state.state[y1][x1].size(); ++i1) { + if (y0 == y1 && x0 == x1 && i0 <= i1) continue; const auto &pState1 = state.state[y1][x1][i1]; + auto &pDer1 = derivativeCache.state[y1][x1][i1]; Vector2D f = pState0.computeForce(pState1); pDer0.v += f; + pDer1.v -= f; } + } + } pDer0.v.x += diff.boundForce(-pState0.pos.x); pDer0.v.x -= diff.boundForce(pState0.pos.x - diff.config.sizeX); diff --git a/Lib/Universe.h b/Lib/Universe.h index 7e3e390..97ac39e 100644 --- a/Lib/Universe.h +++ b/Lib/Universe.h @@ -51,25 +51,26 @@ struct UniverseDifferentiator { UniverseDifferentiator(const UniverseConfig &config, const std::vector &_types); void prepareDifferentiation(UniverseState &state) const; // Has to be called once before every iteration - void derivative(UniverseState &der, UniverseState &state) const; + void derivative(UniverseState &der, UniverseState &derivativeCache, UniverseState &state) const; private: void initForces(UniverseState &der, const UniverseState &state) const; - void computeForces(UniverseState &der, const UniverseState &state) const; - void forcesToAccel(UniverseState &der) const; + void computeForces(UniverseState &der, UniverseState &derivativeCache, const UniverseState &state) const; + void forcesToAccel(UniverseState &der, const UniverseState &derivativeCache) const; double boundForce(double overEdge) const; public: class ParallelForces : public cv::ParallelLoopBody { public: - ParallelForces(const UniverseDifferentiator &diff, UniverseState &der, const UniverseState &state); + ParallelForces(const UniverseDifferentiator &diff, + UniverseState &der, UniverseState &derivativeCache, const UniverseState &state); void operator()(const cv::Range& range) const; ParallelForces& operator=(const ParallelForces &); private: const UniverseDifferentiator &diff; - UniverseState &der; + UniverseState &der, &derivativeCache; const UniverseState &state; }; }; diff --git a/Tests/IntegratorTest.cpp b/Tests/IntegratorTest.cpp index bad022d..91df220 100644 --- a/Tests/IntegratorTest.cpp +++ b/Tests/IntegratorTest.cpp @@ -11,7 +11,7 @@ class Linear { void prepareDifferentiation(double) const { } - void derivative(double &der, double x) const { + void derivative(double &der, double, double x) const { der = x; } };