move a couple of functinos to a new class utils

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})

lib/nBodySim.cpp

 #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) {

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_;

lib/utils.cpp

+#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(); 
+}

lib/utils.hpp

+#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