Skip to content

Commit

Permalink
Avoid checking collisions and computing forces twice between pairs of…
Browse files Browse the repository at this point in the history 
… particles
  • Loading branch information
@kongaskristjan
kongaskristjan committed Oct 12, 2019
1 parent 3b44952 commit 3884783
Show file tree
Hide file tree
Showing 4 changed files with 40 additions and 23 deletions.
12 changes: 7 additions & 5 deletions Lib/Integrators.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand All @@ -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;
Expand All @@ -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.;
Expand All @@ -72,15 +74,15 @@ 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.;
x += xAdditive;

// Compute k4 and add to x
k3 += xInitial;
diff.derivative(k4, k3);
diff.derivative(k4, derivativeCache, k3);
k4 *= dT;
xAdditive = k4;
xAdditive *= 1. / 6.;
Expand Down
38 changes: 26 additions & 12 deletions Lib/Universe.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand All @@ -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();
}
}
Expand All @@ -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 {
Expand All @@ -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);
Expand Down
11 changes: 6 additions & 5 deletions Lib/Universe.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,25 +51,26 @@ struct UniverseDifferentiator {

UniverseDifferentiator(const UniverseConfig &config, const std::vector<ParticleType> &_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;
};
};
Expand Down
2 changes: 1 addition & 1 deletion Tests/IntegratorTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
};
Expand Down

0 comments on commit 3884783

Please sign in to comment.