Add content on performance benchmarks.

benchmarks/cpp/main.cpp

@@ -0,0 +1,195 @@
+#include <iostream>
+#include <math.h>
+#include <numeric>
+#include <vector>
+
+typedef double real;
+using namespace std;
+
+// Class Star, contains the properties:
+// mass (m)
+// position (r)
+// velocity (v)
+// accelerations (a and a0)
+
+class Star {
+public:
+  real m;
+  vector<real> r;
+  vector<real> v;
+  vector<real> a, a0;
+  // Default constructor
+  Star() {
+    r.assign(3, 0);
+    v.assign(3, 0);
+    a.assign(3, 0);
+    a0.assign(3, 0);
+  }
+  // Detailed constructor
+  Star(real mass, vector<real> pos, vector<real> vel) {
+    m = mass;
+    r = pos;
+    v = vel;
+  }
+  // Print function (overloaded << operator)
+  friend ostream &operator<<(ostream &so, const Star &si) {
+    so << si.m << " " << si.r[0] << " " << si.r[1] << " " << si.r[2] << " "
+       << si.v[0] << " " << si.v[1] << " " << si.v[2] << endl;
+    return so;
+  }
+};
+// Star cluster based on the Star class
+// A star cluster contains a number of stars
+// stored in the vector S
+//
+class Cluster : public Star {
+protected:
+public:
+  vector<Star> s;
+  Cluster() : Star() {}
+  // Computes the acceleration of each star in the cluster
+  void acceleration() {
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si)
+      si->a.assign(3, 0);
+    // For each star
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
+      vector<real> rij(3);
+      real init = 0.0;
+      // For each remaining star
+      for (vector<Star>::iterator sj = s.begin(); sj != s.end(); ++sj) {
+        if (si != sj) {
+          // Distance difference between the two stars
+          for (int i = 0; i != 3; ++i)
+            rij[i] = si->r[i] - sj->r[i];
+          // Sum of the dot product
+          real RdotR = inner_product(rij.begin(), rij.end(), rij.begin(), init);
+          real apre = 1. / sqrt(RdotR * RdotR * RdotR);
+          // Update accelerations
+          for (int i = 0; i != 3; ++i) {
+            si->a[i] -= sj->m * apre * rij[i];
+          }
+        } // end for
+      }   // si != sj
+    }     // end for
+  }       // end acceleration
+
+  // Update positions
+  void updatePositions(real dt) {
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
+      // Update the positions, based on the calculated accelerations and
+      // velocities
+      si->a0 = si->a;
+      for (int i = 0; i != 3; ++i) // for each axis (x/y/z)
+        si->r[i] += dt * si->v[i] + 0.5 * dt * dt * si->a0[i];
+    }
+  }
+
+  // Update velocities based on previous and new accelerations
+  void updateVelocities(real dt) {
+    // Update the velocities based on the previous and old accelerations
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
+      for (int i = 0; i != 3; ++i)
+        si->v[i] += 0.5 * dt * (si->a0[i] + si->a[i]);
+      si->a0 = si->a;
+    }
+  }
+
+  // Compute the energy of the system,
+  // contains an expensive O(N^2) part which can be moved to the acceleration
+  // part where this is already calculated
+  vector<real> energies() {
+    real init = 0;
+    vector<real> E(3), rij(3);
+    E.assign(3, 0);
+
+    // Kinetic energy
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si)
+      E[1] += 0.5 * si->m *
+              inner_product(si->v.begin(), si->v.end(), si->v.begin(), init);
+
+    // Potential energy
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
+      for (vector<Star>::iterator sj = si + 1; sj != s.end(); ++sj) {
+        for (int i = 0; i != 3; ++i)
+          rij[i] = si->r[i] - sj->r[i];
+        E[2] -= si->m * sj->m /
+                sqrt(inner_product(rij.begin(), rij.end(), rij.begin(), init));
+      }
+    }
+    E[0] = E[1] + E[2];
+    return E;
+  }
+
+  // Print function
+  friend ostream &operator<<(ostream &so, Cluster &cl) {
+    for (vector<Star>::iterator si = cl.s.begin(); si != cl.s.end(); ++si)
+      so << *si;
+    return so;
+  }
+};
+
+int main(int argc, char* argv[]) {
+
+  Cluster cl;
+  real m;
+  int dummy;
+  vector<real> r(3), v(3);
+
+  // Read input data from the command line (makeplummer | dumbp)
+  do {
+    cin >> dummy;
+    cin >> m;
+    for (int i = 0; i != 3; ++i)
+      cin >> r[i];
+    for (int i = 0; i != 3; ++i)
+      cin >> v[i];
+    cl.s.push_back(Star(m, r, v));
+  } while (!cin.eof());
+
+  // Remove the last one
+  cl.s.pop_back();
+
+  // Compute initial energu of the system
+  vector<real> E(3), E0(3);
+  E0 = cl.energies();
+  cerr << "Energies: " << E0[0] << " " << E0[1] << " " << E0[2] << endl;
+
+  // Start time, end time and simulation step
+  real t = 0.0;
+  real tend;
+
+  if (argc > 1)
+    tend = strtod(argv[1], NULL);
+  else
+    tend = 10.0;
+
+  real dt = 1e-3;
+  int k = 0;
+
+  // Initialize the accelerations
+  cl.acceleration();
+
+  // Start the main loop
+  while (t < tend) {
+    // Update positions based on velocities and accelerations
+    cl.updatePositions(dt);
+
+    // Get new accelerations
+    cl.acceleration();
+
+    // Update velocities
+    cl.updateVelocities(dt);
+
+    t += dt;
+    k += 1;
+    if (k % 100 == 0) {
+      E = cl.energies();
+      cout << "t= " << t << " E= " << E[0] << " " << E[1] << " " << E[2]
+           << " dE/E = " << (E[0] - E0[0]) / E0[0] << endl;
+    }
+  } // end while
+
+  cout << "number time steps: " << k << endl;
+
+  return 0;
+} // end program

benchmarks/python/__init__.py

@@ -0,0 +1 @@
+

benchmarks/python/bench.py

@@ -0,0 +1,138 @@
+import sys
+import math
+import time
+from datetime import timedelta
+
+import numpy as np
+import pandas as pd
+
+from transonic import boost
+
+def load_input_data(path):
+    df = pd.read_csv(
+        path, names=["mass", "x", "y", "z", "vx", "vy", "vz"], delimiter=r"\s+"
+    )
+
+    masses = np.ascontiguousarray(df["mass"].values)
+    positions = np.ascontiguousarray(df.loc[:, ["x", "y", "z"]].values)
+    velocities = np.ascontiguousarray(df.loc[:, ["vx", "vy", "vz"]].values)
+
+    return masses, positions, velocities
+
+
+def advance_positions(positions, velocities, accelerations, time_step):
+    positions += time_step * velocities + 0.5 * accelerations * time_step ** 2
+
+
+def advance_velocities(velocities, accelerations, accelerations1, time_step):
+    velocities += 0.5 * time_step * (accelerations + accelerations1)
+
+
+def compute_distance(vec):
+    return math.sqrt(sum(vec ** 2))
+
+def compute_accelerations(accelerations, masses, positions):
+    number_of_particles = masses.size
+    for index_p0 in range(number_of_particles - 1):
+        position0 = positions[index_p0]
+        masses0 = masses[index_p0]
+        vector = np.empty(3)
+        for index_p1 in range(index_p0 + 1, number_of_particles):
+            masses1 = masses[index_p1]
+            position1 = positions[index_p1]
+            for i in range(3):
+                vector[i] = position0[i] - position1[i]
+
+            distance = sum(vector ** 2) * math.sqrt(sum(vector ** 2))
+            coef_m1 = masses0 / distance
+            coef_m2 = masses1 / distance
+            for i in range(3):
+                accelerations[index_p0][i] -= coef_m1 * vector[i]
+                accelerations[index_p1][i] += coef_m2 * vector[i]
+    return accelerations
+
+
+@boost
+def pythran_loop(
+    time_step: float,
+    nb_steps: int,
+    masses: "float[]",
+    positions: "float[:,:]",
+    velocities: "float[:,:]",
+):
+
+    accelerations = np.zeros_like(positions)
+    accelerations1 = np.zeros_like(positions)
+
+    accelerations = compute_accelerations(accelerations, masses, positions)
+
+    time = 0.0
+    energy0, _, _ = compute_energies(masses, positions, velocities)
+    energy_previous = energy0
+
+    for step in range(nb_steps):
+        advance_positions(positions, velocities, accelerations, time_step)
+        # swap acceleration arrays
+        accelerations, accelerations1 = accelerations1, accelerations
+        accelerations.fill(0)
+        compute_accelerations(accelerations, masses, positions)
+        advance_velocities(velocities, accelerations, accelerations1, time_step)
+        time += time_step
+
+        if not step % 100:
+            energy, _, _ = compute_energies(masses, positions, velocities)
+            # f-strings supported by Pythran>=0.9.8
+            #print(
+            #    f"t = {time_step * step:5.2f}, E = {energy:.7f}, "
+            #    f"dE/E = {(energy - energy_previous) / energy_previous:+.7f}"
+            #)
+            energy_previous = energy
+
+    return math.floor((10000000*energy)/10000000), math.floor((10000000*energy0)/10000000)
+    #return energy, energy0
+
+def compute_kinetic_energy(masses, velocities):
+    return 0.5 * np.sum(masses * np.sum(velocities ** 2, 1))
+
+
+def compute_potential_energy(masses, positions):
+    nb_particules = masses.size
+    pe = 0.0
+    for index_p0 in range(nb_particules - 1):
+        mass0 = masses[index_p0]
+        for index_p1 in range(index_p0 + 1, nb_particules):
+            mass1 = masses[index_p1]
+            vector = positions[index_p0] - positions[index_p1]
+            distance = compute_distance(vector)
+            pe -= (mass0 * mass1) / distance
+    return pe
+
+
+def compute_energies(masses, positions, velocities):
+    energy_kin = compute_kinetic_energy(masses, velocities)
+    energy_pot = compute_potential_energy(masses, positions)
+    return energy_kin + energy_pot, energy_kin, energy_pot
+
+
+if __name__ == "__main__":
+
+    try:
+        time_end = float(sys.argv[2])
+    except IndexError:
+        time_end = 10.0
+
+    time_step = 0.001
+    nb_steps = int(time_end / time_step) + 1
+
+    path_input = sys.argv[1]    
+    masses, positions, velocities = load_input_data(path_input)
+
+    start = time.time()
+    for i in range(5):
+        energy, energy0 = pythran_loop(time_step, nb_steps, masses, positions, velocities)
+    end = time.time()
+
+    print(f"Final dE/E = {(energy - energy0) / energy0:.7e}")
+    print(
+        f"{nb_steps} time steps run in {timedelta(seconds=end-start)}"
+    )