diff --git a/black_holes.cpp b/black_holes.cpp index 2e6f5f7..7790a3b 100644 --- a/black_holes.cpp +++ b/black_holes.cpp @@ -10,39 +10,64 @@ double L[3]; // needed in pn_bh_spin.c #undef SQR /* END legacy inclusion */ -void two_body_gravity(const Particle_ref& j, const Particle_ref& i, const double eps, double& pot, double3& acc, double3& jrk) +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) { - double3 dx = i.x - j.x; - double3 dv = i.v - j.v; + double3 dx = x1 - x2; + double3 dv = v1 - v2; double r2 = dx.norm2() + eps*eps; double r = sqrt(r2); double r3 = r2*r; double r4 = r2*r2; double RP = 3*(dx*dv)/r; - pot = -j.m/r; - acc = -j.m*dx/r3; - jrk = -j.m*(dv/r3 - RP*dx/r4); + 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); } -void Black_hole_physics::adjust_softening(const std::vector& particles) +void Black_hole_physics::adjust_softening( + double& pot1, double& pot2, + double3& acc1, double3& acc2, + double3& jrk1, double3& jrk2) { if (eps_new < 0) return; - 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; - } - } + // 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; } -#if 0 void Black_hole_physics::adjust_post_newtonian( const double dt_bh, // pn_usage should be const double3& acc1, double3& acc2, @@ -51,8 +76,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(masses[0], x1, v1, bbh_grav.spin1, - masses[1], x2, v2, bbh_grav.spin2, + tmp = calc_force_pn_BH(m1, x1, v1, bbh_grav.spin1, + m2, 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); @@ -73,22 +98,72 @@ void Black_hole_physics::adjust_post_newtonian( jrk1 += jrk1_corr; jrk2 += jrk2_corr; } -#endif -void Black_hole_physics::write_bh_data(const double time_cur, const int count, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &a, const std::vector &adot, const std::vector &dt) +void Black_hole_physics::write_bh_data(double time_cur, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &a, const std::vector &adot, const std::vector &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"); - for (int i = 0; i < count; i++) { + 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) { 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]); + 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"); } - fprintf(out, "\n"); fclose(out); } diff --git a/black_holes.h b/black_holes.h index 3d79184..121ae7a 100644 --- a/black_holes.h +++ b/black_holes.h @@ -3,17 +3,6 @@ #include #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]; @@ -24,8 +13,8 @@ class Black_hole_physics { public: Black_hole_physics() : count(0), c(0) {} - Black_hole_physics(const int count, const int myRank, const int rootRank) - : count(count), c(0), myRank(myRank), rootRank(rootRank) {} + 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) {} void set_post_newtonian(const double c, const int pn_usage[7]) { this->c = c; @@ -36,23 +25,34 @@ 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(const std::vector& particles); + void adjust_softening( + double& pot1, double& pot2, + double3& acc1, double3& acc2, + double3& jrk1, double3& jrk2); 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(const double time_cur, const int count, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &a, const std::vector &adot, const std::vector &dt); + void write_bh_data(double time_cur, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &a, const std::vector &adot, const std::vector &dt); public: //TODO make private - /////////////std::vector masses; + double m1, m2; 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; @@ -101,4 +101,4 @@ private: const std::vector &x, &v; std::vector inf_event; FILE *out; -}; +}; \ No newline at end of file diff --git a/config.cpp b/config.cpp index 86e704e..db40311 100644 --- a/config.cpp +++ b/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) - throw std::runtime_error("The number of live black holes (live_smbh_count) has to be greater than or equals to zero"); + 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 (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)) diff --git a/phigrape.cpp b/phigrape.cpp index 52ec687..4dd4865 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -2,7 +2,6 @@ #include #include #include -#include #include "black_holes.h" #include "config.h" @@ -436,10 +435,11 @@ int main(int argc, char *argv[]) calc_self_grav(time_cur, N, ind, x, v, pot, a, adot); Black_hole_physics black_hole_physics; - std::vector 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_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]); 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,18 +447,13 @@ int main(int argc, char *argv[]) double eps = config.eps; #endif black_hole_physics.set_softening(eps, config.live_smbh_custom_eps); - 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); + black_hole_physics.adjust_softening(pot[0], pot[1], a[0], a[1], adot[0], adot[1]); } } 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 pot_ext(N, 0.); @@ -537,7 +532,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, config.live_smbh_count, m, x, v, pot, a, adot, dt); + if (config.live_smbh_output) black_hole_physics.write_bh_data(time_cur, m, x, v, pot, a, adot, dt); if (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur); } /* if (myRank == rootRank) */ @@ -552,23 +547,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 */ - smbh_list.clear(); + int i_bh1=0, i_bh2=1; #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= 2) { - if (config.live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(smbh_list); - #if 0 + 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.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 */ @@ -629,7 +623,10 @@ int main(int argc, char *argv[]) } /* i */ /* define the min. dt over all the act. part. and set it also for the BH... */ - for (int i=0; i < config.live_smbh_count; i++) dt[i] = min_dt; + 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; + } if (config.binary_smbh_influence_sphere_output && (myRank == rootRank)) binary_smbh_influence_sphere_output(ind_act, n_act, timesteps, time_cur); @@ -655,7 +652,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, config.live_smbh_count, m, x, v, pot, a, adot, dt); + if (config.live_smbh_output) black_hole_physics.write_bh_data(time_cur, m, x, v, pot, a, adot, dt); /* Write BH NB data... */ if (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur);