From 30115a72e952cb55f3de958cb16f6e9dc6fb09e5 Mon Sep 17 00:00:00 2001
From: "armindamon.riess" <armindamon.riess@uzh.ch>
Date: Wed, 30 Nov 2022 20:26:21 +0100
Subject: [PATCH] further improvement to computation time
---
lib/nBodySim.cpp | 50 +++++++++++++++++++++++++++++-------------------
lib/nBodySim.hpp | 2 +-
2 files changed, 31 insertions(+), 21 deletions(-)
diff --git a/lib/nBodySim.cpp b/lib/nBodySim.cpp
index adb2980..563cecd 100644
--- a/lib/nBodySim.cpp
+++ b/lib/nBodySim.cpp
@@ -58,45 +58,55 @@ nBodySim::~nBodySim() {
void nBodySim::runSimulation(double dt, unsigned nSteps) {
for (unsigned i=0; i < nSteps; ++i) {
calculateForces();
- updatePositions(dt);
+ doTimeStep(dt);
}
}
-///*
+
void nBodySim::calculateForces() {
#pragma omp parallel for
for (unsigned i=0; i < 3*nParticles_; ++i) {
forces_[i] = 0.0;
}
// loop over each pair of particles and calculate the force between them
- #pragma omp parallel for
+ #pragma omp parallel for schedule(guided)
for (unsigned i = 0; i < nParticles_; i++) {
+ double mi = masses_[i];
+ double s = softening_[i];
+ double px = positions_[3*i];
+ double py = positions_[3*i+1];
+ double pz = positions_[3*i+2];
+ double fx = 0.0;
+ double fy = 0.0;
+ double fz = 0.0;
for (unsigned j = i+1; j < nParticles_; j++) {
+ // if (j <= i) continue;
// calculate distance between particles
- double dx = positions_[3*j] - positions_[3*i];
- double dy = positions_[3*j+1] - positions_[3*i+1];
- double dz = positions_[3*j+2] - positions_[3*i+2];
+ double dx = positions_[3*j] - px;
+ double dy = positions_[3*j+1] - py;
+ double dz = positions_[3*j+2] - pz;
double r2 = dx*dx + dy*dy + dz*dz;
double r = std::sqrt(r2);
dx /= r;
dy /= r;
dz /= r;
- // calculate force
- double mi = masses_[i];
+ // calculate forces
double mj = masses_[j];
- double s = softening_[i];
- double fx = mi*mj * dx / (r2 + s*s);
- double fy = mi*mj * dy / (r2 + s*s);
- double fz = mi*mj * dz / (r2 + s*s);
- forces_[3*i] += fx;
- forces_[3*i+1] += fy;
- forces_[3*i+2] += fz;
- forces_[3*j] -= fx;
- forces_[3*j+1] -= fy;
- forces_[3*j+2] -= fz;
+ double fxij = mi*mj * dx / (r2 + s*s);
+ double fyij = mi*mj * dy / (r2 + s*s);
+ double fzij = mi*mj * dz / (r2 + s*s);
+ fx += fxij;
+ fy += fyij;
+ fz += fzij;
+ forces_[3*j] -= fxij;
+ forces_[3*j+1] -= fyij;
+ forces_[3*j+2] -= fzij;
}
+ forces_[3*i] += fx;
+ forces_[3*i+1] += fy;
+ forces_[3*i+2] += fz;
}
}
-//*/
+
/*
void nBodySim::calculateForces() {
#pragma omp parallel for
@@ -129,7 +139,7 @@ void nBodySim::calculateForces() {
}
*/
-void nBodySim::updatePositions(double dt) {
+void nBodySim::doTimeStep(double dt) {
#pragma omp parallel for
for (unsigned i=0; i < nParticles_; ++i) {
for (unsigned j=0; j < 3; ++j) {
diff --git a/lib/nBodySim.hpp b/lib/nBodySim.hpp
index df1fa4e..d296f15 100644
--- a/lib/nBodySim.hpp
+++ b/lib/nBodySim.hpp
@@ -14,7 +14,7 @@ public:
// loops over all pairs of particles and calculates forces
void calculateForces();
// updates positions and velocities using the calculated forces
- void updatePositions(double dt);
+ void doTimeStep(double dt);
// writes data to file, every row is a particle
void write2file(const double* array, std::string filename, unsigned dim) const;
--
GitLab