Added #include <stdexcept>

black_holes.cpp

@@ -10,64 +10,39 @@ double L[3]; // needed in pn_bh_spin.c
 #undef SQR
 /* END legacy inclusion */
 
-void two_body_gravity(
-    const double m1, const double3 x1, const double3 v1,
-    const double m2, const double3 x2, const double3 v2,
-    const double eps,
-    double& pot1, double3& acc1, double3& jrk1,
-    double& pot2, double3& acc2, double3& jrk2)
+void two_body_gravity(const Particle_ref& j, const Particle_ref& i, const double eps, double& pot, double3& acc, double3& jrk)
 {
-    double3 dx = x1 - x2;
-    double3 dv = v1 - v2;
+    double3 dx = i.x - j.x;
+    double3 dv = i.v - j.v;
     double r2 = dx.norm2() + eps*eps;
     double r  = sqrt(r2);
     double r3 = r2*r;
     double r4 = r2*r2;
     double RP = 3*(dx*dv)/r;
 
-    pot1 = -m2/r;
-    pot2 = -m1/r;
-    acc1 = -m2*dx/r3;
-    acc2 =  m1*dx/r3;
-    jrk1 = -m2*(dv/r3 - RP*dx/r4);
-    jrk2 =  m1*(dv/r3 - RP*dx/r4);
+    pot = -j.m/r;
+    acc = -j.m*dx/r3;
+    jrk = -j.m*(dv/r3 - RP*dx/r4);
 }
 
-void Black_hole_physics::adjust_softening(
-    double& pot1, double& pot2,
-    double3& acc1, double3& acc2,
-    double3& jrk1, double3& jrk2)
+void Black_hole_physics::adjust_softening(const std::vector<Particle_ref>& particles)
 {
     if (eps_new < 0) return;
-    // calculate and "minus" the BH <-> BH softened pot, acc & jerk
-    two_body_gravity(
-        m1, x1, v1,
-        m2, x2, v2,
-        eps_old,
-        bbh_grav.pot1, bbh_grav.a1, bbh_grav.adot1,
-        bbh_grav.pot2, bbh_grav.a2, bbh_grav.adot2);
-    pot1 -= bbh_grav.pot1;
-    pot2 -= bbh_grav.pot2;
-    acc1 -= bbh_grav.a1;
-    acc2 -= bbh_grav.a2;
-    jrk1 -= bbh_grav.adot1;
-    jrk2 -= bbh_grav.adot2;
-
-    // calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk
-    two_body_gravity(
-        m1, x1, v1,
-        m2, x2, v2,
-        eps_new,
-        bbh_grav.pot1, bbh_grav.a1, bbh_grav.adot1,
-        bbh_grav.pot2, bbh_grav.a2, bbh_grav.adot2);
-    pot1 += bbh_grav.pot1;
-    pot2 += bbh_grav.pot2;
-    acc1 += bbh_grav.a1;
-    acc2 += bbh_grav.a2;
-    jrk1 += bbh_grav.adot1;
-    jrk2 += bbh_grav.adot2;
+    double pot_old, pot_new;
+    double3 acc_old, acc_new, jrk_old, jrk_new;
+    for (int j = 0; j < particles.size(); j++) {
+        for (int i = 0; i < particles.size(); i++) {
+            if (i == j) continue;
+            two_body_gravity(particles[j], particles[i], eps_old, pot_old, acc_old, jrk_old);
+            two_body_gravity(particles[j], particles[i], eps_new, pot_new, acc_new, jrk_new);
+            particles[i].pot += pot_new - pot_old;
+            particles[i].acc += acc_new - acc_old;
+            particles[i].jrk += jrk_new - jrk_old;
+        }
+    }
 }
 
+#if 0
 void Black_hole_physics::adjust_post_newtonian(
     const double dt_bh, // pn_usage should be const
     double3& acc1, double3& acc2,
@@ -76,8 +51,8 @@ void Black_hole_physics::adjust_post_newtonian(
     // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk
     // TODO maybe have the PN terms as local variables here?
     int tmp;
-    tmp = calc_force_pn_BH(m1, x1, v1, bbh_grav.spin1,
-                           m2, x2, v2, bbh_grav.spin2,
+    tmp = calc_force_pn_BH(masses[0], x1, v1, bbh_grav.spin1,
+                           masses[1], x2, v2, bbh_grav.spin2,
                            c,  dt_bh, pn_usage,
                            (double(*)[3])bbh_grav.a_pn1, (double(*)[3])bbh_grav.adot_pn1,
                            (double(*)[3])bbh_grav.a_pn2, (double(*)[3])bbh_grav.adot_pn2, myRank, rootRank);
@@ -98,72 +73,22 @@ void Black_hole_physics::adjust_post_newtonian(
     jrk1 += jrk1_corr;
     jrk2 += jrk2_corr;
 }
+#endif
 
-void Black_hole_physics::write_bh_data(double time_cur, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &a, const std::vector<double3> &adot, const std::vector<double> &dt)
+void Black_hole_physics::write_bh_data(const double time_cur, const int count, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &a, const std::vector<double3> &adot, const std::vector<double> &dt)
 {
-    // This function logs data on the black hole(s). It uses both external data
-    // (the arguments to this function) and optionall internal data to this
-    // object (the most recently calculated force and post-Newtonian terms).
     auto out = fopen("bh.dat", "a");
-    if (count == 2) {
-        for (int i=0; i < 2; i++) {
-            double3 *a_pn, *adot_pn, a_bh, adot_bh;
-            double pot_bh;
-            if (i==0) {
-                a_pn    = bbh_grav.a_pn1;
-                adot_pn = bbh_grav.adot_pn1;
-                pot_bh  = bbh_grav.pot1;
-                a_bh    = bbh_grav.a1;
-                adot_bh = bbh_grav.adot1;
-            } else {
-                a_pn    = bbh_grav.a_pn2;
-                adot_pn = bbh_grav.adot_pn2;
-                pot_bh  = bbh_grav.pot2;
-                a_bh    = bbh_grav.a2;
-                adot_bh = bbh_grav.adot2;
-            }
-            if (eps_new >= 0) {
-                fprintf(out, "%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E ",
-                        time_cur, m[i],
-                        x[i][0], x[i][1], x[i][2], x[i].norm(),
-                        v[i][0], v[i][1], v[i][2], v[i].norm(),
-                        pot[i],
-                        a[i][0],    a[i][1],    a[i][2],    a[i].norm(),
-                        adot[i][0], adot[i][1], adot[i][2], adot[i].norm(),
-                        dt[i],
-                        pot_bh,
-                        a_bh[0],    a_bh[1],    a_bh[2],    a_bh.norm(),
-                        adot_bh[0], adot_bh[1], adot_bh[2], adot_bh.norm());
-                if (c > 0) {
-                    fprintf(out, "\t");
-                    for (int pn_idx=0; pn_idx < 7; pn_idx++) {
-                        fprintf(out, "\t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E ", a_pn[pn_idx][0], a_pn[pn_idx][1], a_pn[pn_idx][2], a_pn[pn_idx].norm(), adot_pn[pn_idx][0], adot_pn[pn_idx][1], adot_pn[pn_idx][2], adot_pn[pn_idx].norm());
-                    }
-                }
-                fprintf(out, "\n");
-            } else {
-                fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \n",
-                        time_cur, m[i],
-                        x[i][0], x[i][1], x[i][2], x[i].norm(),
-                        v[i][0], v[i][1], v[i][2], v[i].norm(),
-                        pot[i],
-                        a[i][0],    a[i][1],    a[i][2],    a[i].norm(),
-                        adot[i][0], adot[i][1], adot[i][2], adot[i].norm(),
-                        dt[i]);
-            }
-        }
-        fprintf(out, "\n");
-    } else if (count == 1) {
+    for (int i = 0; i < count; i++) {
         fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \n",
-                time_cur, m[0],
-                x[0][0], x[0][1], x[0][2], x[0].norm(),
-                v[0][0], v[0][1], v[0][2], v[0].norm(),
-                pot[0],
-                a[0][0],    a[0][1],    a[0][2],    a[0].norm(),
-                adot[0][0], adot[0][1], adot[0][2], adot[0].norm(),
-                dt[0]);
-        fprintf(out,"\n");
+            time_cur, m[i],
+            x[i][0], x[i][1], x[i][2], x[i].norm(),
+            v[i][0], v[i][1], v[i][2], v[i].norm(),
+            pot[i],
+            a[i][0],    a[i][1],    a[i][2],    a[i].norm(),
+            adot[i][0], adot[i][1], adot[i][2], adot[i].norm(),
+            dt[i]);
     }
+    fprintf(out, "\n");
     fclose(out);
 }
 

black_holes.h

@@ -3,6 +3,17 @@
 #include <vector>
 #include "double3.h"
 
+struct Particle_ref {
+    Particle_ref(double& m, double3& x, double3& v, double& pot, double3& acc, double3& jrk) :
+        m(m), x(x), v(v), pot(pot), acc(acc), jrk(jrk) {}
+    const double&  m;
+    const double3& x;
+    const double3& v;
+    double&        pot;
+    double3&       acc;
+    double3&       jrk;
+};
+
 struct Bbh_gravity {
     double pot1, pot2;
     double3 a1, a2, adot1, adot2, a_pn1[7], a_pn2[7], adot_pn1[7], adot_pn2[7];
@@ -13,8 +24,8 @@ class Black_hole_physics {
 public:
     Black_hole_physics()
         : count(0), c(0) {}
-    Black_hole_physics(const double m1, const double m2, const int myRank, const int rootRank)
-        : m1(m1), m2(m2), count(2), c(0), myRank(myRank), rootRank(rootRank) {}
+    Black_hole_physics(const int count, const int myRank, const int rootRank)
+        : count(count), c(0), myRank(myRank), rootRank(rootRank) {}
     void set_post_newtonian(const double c, const int pn_usage[7])
     {
         this->c = c;
@@ -25,34 +36,23 @@ public:
         std::copy(spin1, spin1+3, this->bbh_grav.spin1);
         std::copy(spin2, spin2+3, this->bbh_grav.spin2);
     }
-    void set_xv(const double3& x1, const double3& x2, const double3& v1, const double3& v2)
-    {
-        this->x1 = x1;
-        this->x2 = x2;
-        this->v1 = v1;
-        this->v2 = v2;
-    }
     void set_softening(const double eps_old, const double eps_new)
     {
         this->eps_old = eps_old;
         this->eps_new = eps_new;
     }
-    void adjust_softening(
-        double& pot1, double& pot2,
-        double3& acc1, double3& acc2,
-        double3& jrk1, double3& jrk2);
+    void adjust_softening(const std::vector<Particle_ref>& particles);
 
     void adjust_post_newtonian(
         const double dt_bh, // pn_usage should be const
         double3& acc1, double3& acc2,
         double3& jrk1, double3& jrk2);
-    void write_bh_data(double time_cur, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &a, const std::vector<double3> &adot, const std::vector<double> &dt);
+    void write_bh_data(const double time_cur, const int count, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &a, const std::vector<double3> &adot, const std::vector<double> &dt);
 public: //TODO make private
-    double m1, m2;
+    /////////////std::vector<double> masses;
     int count;
     int myRank, rootRank;
     double eps_old, eps_new;
-    double3 x1, v1, x2, v2;
     double c;
     int pn_usage[7];
     Bbh_gravity bbh_grav;

config.cpp

@@ -133,8 +133,8 @@ void Config::error_checking()
     if (output_hdf5)
         throw std::runtime_error("HDF5 output format (output_hdf5=true) requires the code to be compiled with HAS_HDF5");
 #endif
-    if ((live_smbh_count < 0) || (live_smbh_count > 2))
-        throw std::runtime_error("The number of live black holes (live_smbh_count) can be 0, 1, or 2");
+    if (live_smbh_count < 0)
+        throw std::runtime_error("The number of live black holes (live_smbh_count) has to be greater than or equals to zero");
     if (binary_smbh_pn && (live_smbh_count!=2))
         throw std::runtime_error("Post-Newtonian gravity (binary_smbh_pn=true) requires live_smbh_count=2");
     if (binary_smbh_pn && (pn_c <= 0))

phigrape.cpp

@@ -2,6 +2,7 @@
 #include <chrono>
 #include <mpi.h>
 #include <numeric>
+#include <stdexcept>
 
 #include "black_holes.h"
 #include "config.h"
@@ -435,11 +436,10 @@ int main(int argc, char *argv[])
     calc_self_grav(time_cur, N, ind, x, v, pot, a, adot);
 
     Black_hole_physics black_hole_physics;
-    if (config.live_smbh_count == 1)
-        black_hole_physics = Black_hole_physics(m[0], 0, myRank, rootRank);
-    else if (config.live_smbh_count == 2) {
-        black_hole_physics = Black_hole_physics(m[0], m[1], myRank, rootRank);
-        black_hole_physics.set_xv(x[0], x[1], v[0], v[1]);
+    std::vector<Particle_ref> smbh_list;
+    if (config.live_smbh_count >= 1)
+        black_hole_physics = Black_hole_physics(config.live_smbh_count, myRank, rootRank);
+    else if (config.live_smbh_count >= 2) {
         if (config.live_smbh_custom_eps >= 0) {
 #ifdef ETICS
             double eps = (config.grapite_smbh_star_eps >= 0)?config.grapite_smbh_star_eps:config.eps;
@@ -447,13 +447,18 @@ int main(int argc, char *argv[])
             double eps = config.eps;
 #endif
             black_hole_physics.set_softening(eps, config.live_smbh_custom_eps);
-            black_hole_physics.adjust_softening(pot[0], pot[1], a[0], a[1], adot[0], adot[1]);
+            for (int i = 0; i < config.live_smbh_count; i++)
+                smbh_list.emplace_back(m[i], x[i], v[i], pot[i], a[i], adot[i]);
+            black_hole_physics.adjust_softening(smbh_list);
         }
     }
     if (config.binary_smbh_pn) {
+        throw std::runtime_error("This is the triple+ SMBH version, it cannot do PN yet!");
+        #if 0
         black_hole_physics.set_post_newtonian(config.pn_c, config.pn_usage.data());
         if (config.pn_usage[6]) black_hole_physics.set_spins(config.smbh1_spin.data(), config.smbh2_spin.data());
         black_hole_physics.adjust_post_newtonian(dt_min, a[0], a[1], adot[0], adot[1]);
+        #endif
     }
 
     std::vector<double> pot_ext(N, 0.);
@@ -532,7 +537,7 @@ int main(int argc, char *argv[])
     Binary_smbh_influence_sphere_output binary_smbh_influence_sphere_output(config.binary_smbh_influence_radius_factor, N, m, x, v, pot, dt);
     Write_bh_nb_data write_bh_nb_data(config.live_smbh_neighbor_number, config.live_smbh_count, N, m, x, v);
     if (myRank == rootRank) {
-        if (config.live_smbh_output) black_hole_physics.write_bh_data(time_cur, m, x, v, pot, a, adot, dt);
+        if (config.live_smbh_output) black_hole_physics.write_bh_data(time_cur, config.live_smbh_count, m, x, v, pot, a, adot, dt);
         if (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur);
     } /* if (myRank == rootRank) */
 
@@ -547,23 +552,23 @@ int main(int argc, char *argv[])
         active_search.get_active_indices(min_t, t, dt, ind_act, n_act);
 
         /* Find the BH(s) indices in the active list */
-        int i_bh1=0, i_bh2=1;
+        smbh_list.clear();
 #ifdef ETICS
         /* Unlike with the simple active search, with GPU accelerated GRAPite
         active search, the list of active indices is not sorted. */
         int n_bh = config.live_smbh_count;
         if (config.grapite_active_search && (n_bh>0)) {
             int act_def_grapite_bh_count = 0;
-            int i_bh[n_bh];
             for (int i=0; i<n_act; i++) {
                 if (ind_act[i]<n_bh) {
-                    i_bh[ind_act[i]] = i;
-                    if (++act_def_grapite_bh_count == n_bh) break;
+                    smbh_list.emplace_back(m[ind_act[i]], x_act_new[i], v_act_new[i], pot_act_new[i], a_act_new[i], adot_act_new[i]);
+                    if (act_def_grapite_bh_count++ == n_bh) break;
                 }
             }
-            i_bh1 = i_bh[0];
-            if (n_bh == 2) i_bh2 = i_bh[1];
         }
+#else
+        for (int i = 0; i < config.live_smbh_count; i++)
+            smbh_list.emplace_back(m[ind_act[i]], x_act_new[i], v_act_new[i], pot_act_new[i], a_act_new[i], adot_act_new[i]);
 #endif
 
         /* predict the active particles positions etc... on all the nodes */
@@ -572,10 +577,11 @@ int main(int argc, char *argv[])
         /* Calculate gravity on active particles */
         calc_self_grav(min_t, n_act, ind_act, x_act_new, v_act_new, pot_act_new, a_act_new, adot_act_new);
 
-        if (config.live_smbh_count == 2) {
-            black_hole_physics.set_xv(x_act_new[i_bh1], x_act_new[i_bh2], v_act_new[i_bh1], v_act_new[i_bh2]);
-            if (config.live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(pot_act_new[i_bh1], pot_act_new[i_bh2], a_act_new[i_bh1], a_act_new[i_bh2], adot_act_new[i_bh1], adot_act_new[i_bh2]);
+        if (config.live_smbh_count >= 2) {
+            if (config.live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(smbh_list);
+            #if 0
             if (config.binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt[i_bh1], a_act_new[i_bh1], a_act_new[i_bh2], adot_act_new[i_bh1], adot_act_new[i_bh2]);
+            #endif
         }
 
         /* Calculate gravity on active particles due to external forces */
@@ -623,10 +629,7 @@ int main(int argc, char *argv[])
         } /* i */
 
         /* define the min. dt over all the act. part. and set it also for the BH... */
-        if (config.live_smbh_count > 0) {
-            if (config.live_smbh_count>=1) dt[0] = min_dt;
-            if (config.live_smbh_count==2) dt[1] = min_dt;
-        }
+        for (int i=0; i < config.live_smbh_count; i++) dt[i] = min_dt;
 
         if (config.binary_smbh_influence_sphere_output && (myRank == rootRank))
             binary_smbh_influence_sphere_output(ind_act, n_act, timesteps, time_cur);
@@ -652,7 +655,7 @@ int main(int argc, char *argv[])
         if (time_cur >= t_bh) {
             if (myRank == rootRank) {
                 /* Write BH data... */
-                if (config.live_smbh_output) black_hole_physics.write_bh_data(time_cur, m, x, v, pot, a, adot, dt);
+                if (config.live_smbh_output) black_hole_physics.write_bh_data(time_cur, config.live_smbh_count, m, x, v, pot, a, adot, dt);
 
                 /* Write BH NB data... */
                 if (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur);