From b81a635a1178247f75b71fcf71c107766197894c Mon Sep 17 00:00:00 2001
From: "armindamon.riess" <armindamon.riess@uzh.ch>
Date: Wed, 30 Nov 2022 10:26:57 +0100
Subject: [PATCH] correct force calculation, implement OpenMP
---
lib/nBodySim.cpp | 42 +++++++++++++++++++++++++++---------------
lib/nBodySim.hpp | 6 ++++--
2 files changed, 31 insertions(+), 17 deletions(-)
diff --git a/lib/nBodySim.cpp b/lib/nBodySim.cpp
index 89e1c0e..b681ddb 100644
--- a/lib/nBodySim.cpp
+++ b/lib/nBodySim.cpp
@@ -1,4 +1,6 @@
#include "nBodySim.hpp"
+#include <omp.h>
+#include <cmath>
#include <string>
#include <fstream>
#include <iomanip>
@@ -17,14 +19,14 @@ nBodySim::nBodySim(std::string& datafile) {
positions_ = new double[3*N_];
for (unsigned i=0; i < 3; ++i) {
for (unsigned j=0; j < N_; ++j) {
- file >> positions_[i*N_+j];
+ file >> positions_[i+3*j];
}
}
velocities_ = new double[3*N_];
for (unsigned i=0; i < 3; ++i) {
for (unsigned j=0; j < N_; ++j) {
- file >> velocities_[i*N_+j];
+ file >> velocities_[i+3*j];
}
}
@@ -39,10 +41,8 @@ nBodySim::nBodySim(std::string& datafile) {
}
forces_ = new double[3*N_];
- for (unsigned i=0; i < 3; ++i) {
- for (unsigned j=0; j < N_; ++j) {
- file >> forces_[i*N_+j];
- }
+ for (unsigned i=0; i < 3*N_; ++i) {
+ forces_[i] = 0.0;
}
}
@@ -56,29 +56,41 @@ nBodySim::~nBodySim() {
}
void nBodySim::calculateForces() {
+ #pragma omp parallel for
+ for (unsigned i=0; i < 3*N_; ++i) {
+ forces_[i] = 0.0;
+ }
// loop over each pair of particles and calculate the force between them
+ #pragma omp parallel for collapse(2)
for (unsigned i = 0; i < N_; i++) {
for (unsigned j = 0; j < N_; j++) {
if (i == j) continue;
// calculate distance between particles
- double dx = positions_[j] - positions_[i];
- double dy = positions_[j+1] - positions_[i+1];
- double dz = positions_[j+2] - positions_[i+2];
+ 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 r2 = dx*dx + dy*dy + dz*dz;
+ double r = std::sqrt(r2);
+ dx /= r;
+ dy /= r;
+ dz /= r;
// calculate force
- forces_[i] += masses_[j] * dx / (r2 + softening_[i]);
- forces_[i+1] += masses_[j] * dy / (r2 + softening_[i]);
- forces_[i+2] += masses_[j] * dz / (r2 + softening_[i]);
+ double mi = masses_[i];
+ double mj = masses_[j];
+ double s = softening_[i];
+ forces_[3*i] += mi*mj * dx / (r2 + s*s);
+ forces_[3*i+1] += mi*mj * dy / (r2 + s*s);
+ forces_[3*i+2] += mi*mj * dz / (r2 + s*s);
}
}
}
-void nBodySim::write2file(const double* array, std::string filename, unsigned N, unsigned dim) {
+void nBodySim::write2file(const double* array, std::string filename, unsigned dim) const {
// writes array to file, with N rows and dim columns
std::ofstream file(filename);
- for (unsigned i = 0; i < N; i++) {
+ for (unsigned i = 0; i < N_; i++) {
for (unsigned j = 0; j < dim; j++) {
- file << std::right << std::setw(16) << array[i+j] << " ";
+ file << std::right << std::setw(16) << array[dim*i+j] << " ";
}
file << std::endl;
}
diff --git a/lib/nBodySim.hpp b/lib/nBodySim.hpp
index 97c3d1a..c90d70d 100644
--- a/lib/nBodySim.hpp
+++ b/lib/nBodySim.hpp
@@ -10,9 +10,11 @@ public:
~nBodySim();
// methods
void calculateForces();
+ double* getMasses() { return masses_; }
+ double* getPositions() { return positions_; }
double* getForces() { return forces_; }
- void write2file(const double* array, std::string filename, unsigned N, unsigned dim);
- unsigned getN() { return N_; }
+ void write2file(const double* array, std::string filename, unsigned dim) const;
+ unsigned getN() const { return N_; }
private:
unsigned N_;
// 3d vectors are stored as a 1d array of length 3*N_ in the format x1, y1, z1, x2, y2, z2, ...
--
GitLab