diff --git a/CMakeLists.txt b/CMakeLists.txt
index a6401c510279b036255b9befed570129e031b3e9..52f2c169733a92d2eadea0cd423311f6b1de6b62 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -8,11 +8,13 @@ set(CMAKE_CXX_FLAGS "-O3 -Wall")
set(HEADERS lib/nBodySim.hpp
lib/node.hpp
lib/tree.hpp
+ lib/utils.hpp
)
set(SOURCES lib/nBodySim.cpp
lib/node.cpp
lib/tree.cpp
+ lib/utils.cpp
)
add_library(nbodyLib ${HEADERS} ${SOURCES})
diff --git a/lib/nBodySim.cpp b/lib/nBodySim.cpp
index 6e5bda6e4c219eecfd44685dac76b903f255d316..79523614bea1ae14a3dd1093f4d53505eb5edb70 100644
--- a/lib/nBodySim.cpp
+++ b/lib/nBodySim.cpp
@@ -1,6 +1,4 @@
#include "nBodySim.hpp"
-#include "tree.hpp"
-#include "utils.hpp"
#include <omp.h>
#include <cmath>
#include <string>
@@ -52,6 +50,8 @@ nBodySim::nBodySim(std::string datafile) {
Vector center{0.0, 0.0, 0.0};
tree_ = new Tree(mass_, positions_, forces_, softening_, nParticles_, maxCoord, center);
+
+ utils_ = Utilities(nParticles_, positions_, velocities_, forces_, mass_, softening_, tree_);
}
nBodySim::~nBodySim() {
@@ -73,28 +73,28 @@ void nBodySim::runSimulation(double dt, unsigned nSteps) {
for (unsigned i=0; i < nSteps; ++i) {
write2file(positions_, positionsFile, 3, skip);
- energyFile << calculateTotalEnergy() << std::endl;
+ energyFile << utils_.calculateTotalEnergy() << std::endl;
- // Vector COM = getCenterOfMass();
+ // Vector COM = utils_.getCenterOfMass();
// Integrate with leapfrog
auto start = std::chrono::high_resolution_clock::now();
kick(dt/2.0);
drift(dt);
treeUpdateForces();
- // double dev = verifyForces();
+ // double dev = utils_.verifyForces();
kick(dt/2.0);
updateTree();
auto stop = std::chrono::high_resolution_clock::now();
auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(stop - start);
std::cout << std::right << "Step " << std::setw(5) << i
- // << ": MFM = " << std::setw(12) << calculateMeanForceMagnitude()
+ // << ": MFM = " << std::setw(12) << utils_.calculateMeanForceMagnitude()
// << ", DEV = " << std::setw(12) << dev
- // << ", MCD = " << std::setw(8) << calculateMeanCenterDistance()
- // << ", MID = " << std::setw(8) << calculateMeanInterparticleDistance()
+ // << ", MCD = " << std::setw(8) << utils_.calculateMeanCenterDistance()
+ // << ", MID = " << std::setw(8) << utils_.calculateMeanInterparticleDistance()
// << ", COM = " << "(" << COM[0] << ", " << COM[1] << ", " << COM[2] << ")"
- // << ", E = " << std::setw(12) << calculateTotalEnergy()
+ // << ", E = " << std::setw(12) << utils_.calculateTotalEnergy()
<< ", T = " << std::setw(8) << duration.count()/1000.0 << " s"
<< std::endl;
}
@@ -157,44 +157,6 @@ void nBodySim::treeUpdateForces() {
}
}
-double nBodySim::verifyForces() const {
- double avgdev = 0.0;
- // loop over each pair of particles and calculate the force between them
- #pragma omp parallel for // schedule(guided)
- for (unsigned i = 0; i < nParticles_; ++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 = 0; j < nParticles_; ++j) {
- if (j == i) continue;
- // calculate distance between particles
- double dx = positions_[3*j+0] - 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);
-
- // calculate forces
- double factor = mass_*mass_ / (r2 + softening_*softening_);
- double fxij = factor * dx / r;
- double fyij = factor * dy / r;
- double fzij = factor * dz / r;
- fx += fxij;
- fy += fyij;
- fz += fzij;
- }
- double f = std::sqrt(fx*fx + fy*fy + fz*fz);
- double ft = std::sqrt(forces_[3*i]*forces_[3*i] + forces_[3*i+1]*forces_[3*i+1] + forces_[3*i+2]*forces_[3*i+2]);
- double df = std::abs(f - ft);
- avgdev += df;
- }
- avgdev /= nParticles_;
- return avgdev;
-}
-
void nBodySim::kick(double dt) {
#pragma omp parallel for
for (unsigned i=0; i < nParticles_; ++i) {
@@ -217,92 +179,6 @@ void nBodySim::updateTree() {
tree_->update();
}
-double nBodySim::calculateMeanInterparticleDistance() const {
- double meanInterparticleDistance = 0.0;
- // loop over each pair of particles and calculate the distance between them
- #pragma omp parallel for reduction(+:meanInterparticleDistance) // schedule(guided, 64)
- for (unsigned i = 0; i < nParticles_; i++) {
- double px = positions_[3*i];
- double py = positions_[3*i+1];
- double pz = positions_[3*i+2];
- double localMeanInterparticleDistance = 0.0;
- for (unsigned j = i+1; j < nParticles_; j++) {
- // calculate distance between particles
- 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);
- localMeanInterparticleDistance += r;
- }
- meanInterparticleDistance += localMeanInterparticleDistance;
- }
- meanInterparticleDistance /= (nParticles_*(nParticles_-1)/2);
- return meanInterparticleDistance;
-}
-
-double nBodySim::calculateMeanForceMagnitude() const {
- double forceMagnitude = 0.0;
- #pragma omp parallel for reduction(+:forceMagnitude)
- for (unsigned i = 0; i < nParticles_; i++) {
- double fx = forces_[3*i];
- double fy = forces_[3*i+1];
- double fz = forces_[3*i+2];
- double f2 = fx*fx + fy*fy + fz*fz;
- double f = std::sqrt(f2);
- forceMagnitude += f;
- }
- forceMagnitude /= nParticles_;
- return forceMagnitude;
-}
-
-double nBodySim::calculateMeanCenterDistance() const {
- double avgDist = 0.0;
- #pragma omp parallel for reduction(+:avgDist)
- for (unsigned i = 0; i < nParticles_; i++) {
- double x = positions_[3*i];
- double y = positions_[3*i+1];
- double z = positions_[3*i+2];
- double d2 = x*x + y*y + z*z;
- double d = std::sqrt(d2);
- avgDist += d;
- }
- avgDist /= nParticles_;
- return avgDist;
-}
-
-double nBodySim::calculateTotalEnergy() const {
- double totalEnergy = 0.0;
- #pragma omp parallel for reduction(+:totalEnergy)
- for (unsigned i = 0; i < nParticles_; i++) {
- double px = positions_[3*i];
- double py = positions_[3*i+1];
- double pz = positions_[3*i+2];
- double vx = velocities_[3*i];
- double vy = velocities_[3*i+1];
- double vz = velocities_[3*i+2];
- double v2 = vx*vx + vy*vy + vz*vz;
- double kineticEnergy = 0.5 * mass_ * v2;
- double potentialEnergy = 0.0;
- for (unsigned j = 0; j < nParticles_; ++j) {
- if (j == i) continue;
- // calculate distance between particles
- 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);
- potentialEnergy -= mass_*mass_ / r;
- }
- totalEnergy += kineticEnergy + potentialEnergy;
- }
- return totalEnergy;
-}
-
-Vector nBodySim::getCenterOfMass() const {
- return tree_->getCOM();
-}
-
void nBodySim::write2file(const double* array, std::ofstream& file, unsigned dim, unsigned skip) const {
// writes array to file, with N rows and dim columns
for (unsigned i = 0; i < nParticles_; i+=skip) {
diff --git a/lib/nBodySim.hpp b/lib/nBodySim.hpp
index 4b6be7dc9ba9416045eb5966d43a5e956297a489..5f887dbb5c72883a555ae1394d4b3bc15ec3571c 100644
--- a/lib/nBodySim.hpp
+++ b/lib/nBodySim.hpp
@@ -2,6 +2,7 @@
#define NBODYSIM_HPP
#include "tree.hpp"
+#include "utils.hpp"
#include <string>
class nBodySim {
@@ -18,8 +19,6 @@ public:
void updateForces();
// Use treecode and multipole expansion to calculate forces
void treeUpdateForces();
- // Loops over all pairs of particles and calculates forces, checks if forces are correct
- double verifyForces() const;
// Leapfrog: update velocities
void kick(double dt);
@@ -27,18 +26,7 @@ public:
void drift(double dt);
// Update tree
void updateTree();
-
- // Loop over all pairs of particles and calculate mean interparticle distance
- double calculateMeanInterparticleDistance() const;
- // Loop over all force vectors and calculate mean force magnitude
- double calculateMeanForceMagnitude() const;
- // loop over every particle and calculate the mean distance from the center
- double calculateMeanCenterDistance() const;
- // Loop over all particles and calculate the total energy
- double calculateTotalEnergy() const;
- // get center of mass
- Vector getCenterOfMass() const;
-
+
// Writes data to file, every row is a particle
void write2file(const double* array, std::ofstream& file, unsigned dim, unsigned skip) const;
// Write current state of simulation to file
@@ -53,6 +41,8 @@ public:
unsigned getNParticles() const { return nParticles_; }
private:
+ Utilities utils_;
+
Tree* tree_;
unsigned nParticles_;
diff --git a/lib/utils.cpp b/lib/utils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a9c8ac7aff268868b3c78afb5c1b73838aaf8e92
--- /dev/null
+++ b/lib/utils.cpp
@@ -0,0 +1,146 @@
+#include "utils.hpp"
+
+#include <cmath>
+#include <omp.h>
+
+Utilities::Utilities() {
+ nParticles_ = 0;
+ positions_ = nullptr;
+ velocities_ = nullptr;
+ forces_ = nullptr;
+ mass_ = 0.0;
+ tree_ = nullptr;
+}
+
+Utilities::Utilities(unsigned nParticles, double* positions, double* velocities, double* forces, double mass, double softening, Tree* tree) {
+ nParticles_ = nParticles;
+ positions_ = positions;
+ velocities_ = velocities;
+ forces_ = forces;
+ mass_ = mass;
+ tree_ = tree;
+}
+
+double Utilities::verifyForces() const {
+ double avgdev = 0.0;
+ // loop over each pair of particles and calculate the force between them
+ #pragma omp parallel for // schedule(guided)
+ for (unsigned i = 0; i < nParticles_; ++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 = 0; j < nParticles_; ++j) {
+ if (j == i) continue;
+ // calculate distance between particles
+ double dx = positions_[3*j+0] - 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);
+
+ // calculate forces
+ double factor = mass_*mass_ / (r2 + softening_*softening_);
+ double fxij = factor * dx / r;
+ double fyij = factor * dy / r;
+ double fzij = factor * dz / r;
+ fx += fxij;
+ fy += fyij;
+ fz += fzij;
+ }
+ double f = std::sqrt(fx*fx + fy*fy + fz*fz);
+ double ft = std::sqrt(forces_[3*i]*forces_[3*i] + forces_[3*i+1]*forces_[3*i+1] + forces_[3*i+2]*forces_[3*i+2]);
+ double df = std::abs(f - ft);
+ avgdev += df;
+ }
+ avgdev /= nParticles_;
+ return avgdev;
+}
+
+double Utilities::calculateMeanInterparticleDistance() const {
+ double meanInterparticleDistance = 0.0;
+ // loop over each pair of particles and calculate the distance between them
+ #pragma omp parallel for reduction(+:meanInterparticleDistance) // schedule(guided, 64)
+ for (unsigned i = 0; i < nParticles_; i++) {
+ double px = positions_[3*i];
+ double py = positions_[3*i+1];
+ double pz = positions_[3*i+2];
+ double localMeanInterparticleDistance = 0.0;
+ for (unsigned j = i+1; j < nParticles_; j++) {
+ // calculate distance between particles
+ 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);
+ localMeanInterparticleDistance += r;
+ }
+ meanInterparticleDistance += localMeanInterparticleDistance;
+ }
+ meanInterparticleDistance /= (nParticles_*(nParticles_-1)/2);
+ return meanInterparticleDistance;
+}
+
+double Utilities::calculateMeanForceMagnitude() const {
+ double forceMagnitude = 0.0;
+ #pragma omp parallel for reduction(+:forceMagnitude)
+ for (unsigned i = 0; i < nParticles_; i++) {
+ double fx = forces_[3*i];
+ double fy = forces_[3*i+1];
+ double fz = forces_[3*i+2];
+ double f2 = fx*fx + fy*fy + fz*fz;
+ double f = std::sqrt(f2);
+ forceMagnitude += f;
+ }
+ forceMagnitude /= nParticles_;
+ return forceMagnitude;
+}
+
+double Utilities::calculateMeanCenterDistance() const {
+ double avgDist = 0.0;
+ #pragma omp parallel for reduction(+:avgDist)
+ for (unsigned i = 0; i < nParticles_; i++) {
+ double x = positions_[3*i];
+ double y = positions_[3*i+1];
+ double z = positions_[3*i+2];
+ double d2 = x*x + y*y + z*z;
+ double d = std::sqrt(d2);
+ avgDist += d;
+ }
+ avgDist /= nParticles_;
+ return avgDist;
+}
+
+double Utilities::calculateTotalEnergy() const {
+ double totalEnergy = 0.0;
+ #pragma omp parallel for reduction(+:totalEnergy)
+ for (unsigned i = 0; i < nParticles_; i++) {
+ double px = positions_[3*i];
+ double py = positions_[3*i+1];
+ double pz = positions_[3*i+2];
+ double vx = velocities_[3*i];
+ double vy = velocities_[3*i+1];
+ double vz = velocities_[3*i+2];
+ double v2 = vx*vx + vy*vy + vz*vz;
+ double kineticEnergy = 0.5 * mass_ * v2;
+ double potentialEnergy = 0.0;
+ for (unsigned j = 0; j < nParticles_; ++j) {
+ if (j == i) continue;
+ // calculate distance between particles
+ 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);
+ potentialEnergy -= mass_*mass_ / r;
+ }
+ totalEnergy += kineticEnergy + potentialEnergy;
+ }
+ return totalEnergy;
+}
+
+Vector Utilities::getCenterOfMass() const {
+ return tree_->getCOM();
+}
diff --git a/lib/utils.hpp b/lib/utils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8dec913d32c90d19147a03cee1319079aec300e5
--- /dev/null
+++ b/lib/utils.hpp
@@ -0,0 +1,35 @@
+#ifndef UTILS_HPP
+#define UTILS_HPP
+
+#include "tree.hpp"
+
+class Utilities {
+public:
+ // Constructor
+ Utilities();
+ Utilities(unsigned nParticles, double* positions, double* velocities, double* forces, double mass, double softening, Tree* tree);
+
+ // Loops over all pairs of particles and calculates forces, checks if forces are correct
+ double verifyForces() const;
+ // Loop over all pairs of particles and calculate mean interparticle distance
+ double calculateMeanInterparticleDistance() const;
+ // Loop over all force vectors and calculate mean force magnitude
+ double calculateMeanForceMagnitude() const;
+ // Loop over every particle and calculate the mean distance from the center
+ double calculateMeanCenterDistance() const;
+ // Loop over all particles and calculate the total energy
+ double calculateTotalEnergy() const;
+ // Get center of mass
+ Vector getCenterOfMass() const;
+
+private:
+ unsigned nParticles_;
+ double* positions_;
+ double* velocities_;
+ double* forces_;
+ double mass_;
+ double softening_;
+ Tree* tree_;
+};
+
+#endif /* UTILS_HPP */
\ No newline at end of file