diff --git a/black_holes.cpp b/black_holes.cpp index 7790a3b..2e6f5f7 100644 --- a/black_holes.cpp +++ b/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& 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 &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(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) { - // 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); } diff --git a/black_holes.h b/black_holes.h index 121ae7a..3d79184 100644 --- a/black_holes.h +++ b/black_holes.h @@ -3,6 +3,17 @@ #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]; @@ -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& 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 &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(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); public: //TODO make private - double m1, m2; + /////////////std::vector 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; @@ -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 db40311..86e704e 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) || (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)) diff --git a/phigrape.cpp b/phigrape.cpp index 4dd4865..f3d8546 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -435,11 +435,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 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 +446,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 pot_ext(N, 0.); @@ -532,7 +536,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 +551,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= 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 +628,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 +654,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);