diff --git a/Makefile b/Makefile index 6e4538f..06637d0 100644 --- a/Makefile +++ b/Makefile @@ -23,7 +23,7 @@ CPPFLAGS += -DHAS_HDF5 LIB += -lhdf5 -lz -ldl default: - $(MPICXX) $(CPPFLAGS) $(CXXFLAGS) $(INC) external.cpp io.cpp config.cpp phigrape.cpp -o $(EXECUTABLE) $(LIB) + $(MPICXX) $(CPPFLAGS) $(CXXFLAGS) $(INC) black_holes.cpp external.cpp io.cpp config.cpp phigrape.cpp -o $(EXECUTABLE) $(LIB) yebisu: CPPFLAGS := $(filter-out -DETICS, $(CPPFLAGS)) yebisu: default diff --git a/black_holes.cpp b/black_holes.cpp new file mode 100644 index 0000000..48685a6 --- /dev/null +++ b/black_holes.cpp @@ -0,0 +1,255 @@ +#include +#include +#include "black_holes.h" + + +// TODO do not include c files! +#define SQR(x) ((x)*(x)) +double L[3]; // needed in pn_bh_spin.c +#include "n_bh.c" +#include "pn_bh_spin.c" + +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 = 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; + + 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); +} + +// 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) {} +// void set_post_newtonian(const double c, const int pn_usage[7]) +// { +// this->c = c; +// std::copy(pn_usage, pn_usage+7, this->pn_usage); +// } +// void set_spins(const double spin1[3], const double spin2[3]) +// { +// 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_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, double m[], double3 x[], double3 v[], double pot[], double3 a[], double3 adot[], double dt[]); +// public: //TODO make private +// 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; +// }; + +void Black_hole_physics::adjust_softening( + double& pot1, double& pot2, + double3& acc1, double3& acc2, + double3& jrk1, double3& jrk2) +{ + 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; +} + +void Black_hole_physics::adjust_post_newtonian( + const double dt_bh, // pn_usage should be const + double3& acc1, double3& acc2, + double3& jrk1, double3& jrk2) +{ + // 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, + 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); + if (tmp == 505) exit(-1); // Very ugly way to terminate + + // NOTE we have these _corr variables accumulating the corrections before + // applying it. It's almost the same but different from applying each + // correction term in a loop. + double3 acc1_corr(0,0,0), acc2_corr(0,0,0), jrk1_corr(0,0,0), jrk2_corr(0,0,0); + for (int i=1; i<7; i++) { + acc1_corr += bbh_grav.a_pn1[i]; + acc2_corr += bbh_grav.a_pn2[i]; + jrk1_corr += bbh_grav.adot_pn1[i]; + jrk2_corr += bbh_grav.adot_pn2[i]; + } + acc1 += acc1_corr; + acc2 += acc2_corr; + jrk1 += jrk1_corr; + jrk2 += jrk2_corr; +} + +void Black_hole_physics::write_bh_data(double time_cur, double m[], double3 x[], double3 v[], double pot[], double3 a[], double3 adot[], 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) { + 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"); + } + fclose(out); +} + +void write_bh_nb_data(int nb, int smbh_count, double time_cur, int N, double m[], double3 x[], double3 v[]) +{ + //int nb = config->live_smbh_neighbor_number; + // TODO fix below!!! Should be class members. + int ind_sort[30000]; + double var_sort[30000]; + //TODO you don't want and probably don't need to allocate these here. Maybe just need size nb, or maybe use as private variables. + + auto out = fopen("bh_neighbors.dat", "a"); + + /* 1st BH */ + + for (int i_bh=0; i_bh < smbh_count; i_bh++) { + for (int i=0; i +#include "double3.h" + +struct Bbh_gravity { + double pot1, pot2; + double3 a1, a2, adot1, adot2, a_pn1[7], a_pn2[7], adot_pn1[7], adot_pn2[7]; + double spin1[3], spin2[3]; +}; + +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) {} + void set_post_newtonian(const double c, const int pn_usage[7]) + { + this->c = c; + std::copy(pn_usage, pn_usage+7, this->pn_usage); + } + void set_spins(const double spin1[3], const double spin2[3]) + { + 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_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, double m[], double3 x[], double3 v[], double pot[], double3 a[], double3 adot[], double dt[]); +public: //TODO make private + 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; +}; + +void write_bh_nb_data(int nb, int smbh_count, double time_cur, int N, double m[], double3 x[], double3 v[]); diff --git a/phigrape.cpp b/phigrape.cpp index 20bc500..86d8121 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -54,6 +54,7 @@ Version number : 19.04 Last redaction : 2019.04.16 12:55 *****************************************************************************/ #include "double3.h" +#include "black_holes.h" #include "config.h" #include "io.h" #include "external.h" @@ -99,7 +100,7 @@ Config *config; #define N_MAX (32*KB) -double L[3]; // needed in pn_bh_spin.c +// double L[3]; // needed in pn_bh_spin.c // Needed for things related to BHs #include "debug.h" @@ -131,11 +132,11 @@ double3 x_bh1, x_bh2, v_bh1, v_bh2; double pot_bh1, pot_bh2; double3 a_bh1, adot_bh1, a_bh2, adot_bh2; -#include "n_bh.c" +//#include "n_bh.c" -double3 a_pn1[7], adot_pn1[7], a_pn2[7], adot_pn2[7]; +//double3 a_pn1[7], adot_pn1[7], a_pn2[7], adot_pn2[7]; -#include "pn_bh_spin.c" +//#include "pn_bh_spin.c" #ifdef ETICS int grapite_cep_index; @@ -164,248 +165,6 @@ void get_CPU_time(double *time_real, double *time_user, double *time_syst) *time_user = *time_real; } - - -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 = 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; - - 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); -} - - - -class Black_hole_physics { -public: - Black_hole_physics() - : count(0) {} - Black_hole_physics(const double m1, const double m2, const int myRank, const int rootRank) - : m1(m1), m2(m2), count(2), myRank(myRank), rootRank(rootRank) {} - void set_post_newtonian(const double c, const int pn_usage[7]) - { - this->c = c; - std::copy(pn_usage, pn_usage+7, this->pn_usage); - } - void set_spins(const double spin1[3], const double spin2[3]) - { - std::copy(spin1, spin1+3, this->spin1); - std::copy(spin2, spin2+3, this->spin2); - } - void adjust_softening( - const double eps_old, const double eps_new, - const double3& x1, const double3& x2, - const double3& v1, const double3& v2, - 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& x1, double3& x2, - double3& v1, double3& v2, - double& pot1, double& pot2, - double3& acc1, double3& acc2, - double3& jrk1, double3& jrk2); -public: //TODO make private - double m1, m2; - int count; - int myRank, rootRank; - double c; - int pn_usage[7]; - double3 a_pn1[7], a_pn2[7], adot_pn1[7], adot_pn2[7]; - double spin1[3], spin2[3]; -}; - -void Black_hole_physics::adjust_softening( - const double eps_old, const double eps_new, - const double3& x1, const double3& x2, - const double3& v1, const double3& v2, - double& pot1, double& pot2, - double3& acc1, double3& acc2, - double3& jrk1, double3& jrk2) -{ - double pot_bh1, pot_bh2; - double3 a_bh1, a_bh2, adot_bh1, adot_bh2; - // calculate and "minus" the BH <-> BH softened pot, acc & jerk - two_body_gravity( - m1, x1, v1, - m2, x2, v2, - eps_old, - pot_bh1, a_bh1, adot_bh1, - pot_bh2, a_bh2, adot_bh2); - pot1 -= pot_bh1; - pot2 -= pot_bh2; - acc1 -= a_bh1; - acc2 -= a_bh2; - jrk1 -= adot_bh1; - jrk2 -= adot_bh2; - - // calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk - two_body_gravity( - m1, x1, v1, - m2, x2, v2, - eps_new, - pot_bh1, a_bh1, adot_bh1, - pot_bh2, a_bh2, adot_bh2); - pot1 += pot_bh1; - pot2 += pot_bh2; - acc1 += a_bh1; - acc2 += a_bh2; - jrk1 += adot_bh1; - jrk2 += adot_bh2; -} - -void Black_hole_physics::adjust_post_newtonian( - const double dt_bh, // pn_usage should be const - double3& x1, double3& x2, - double3& v1, double3& v2, - double& pot1, double& pot2, - double3& acc1, double3& acc2, - double3& jrk1, double3& jrk2) -{ - // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk - int tmp; - tmp = calc_force_pn_BH(m1, x1, v1, spin1, - m2, x2, v2, spin2, - c, dt_bh, pn_usage, - (double(*)[3])a_pn1, (double(*)[3])adot_pn1, - (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank); - if (tmp == 505) { - exit(-1); // Very ugly way to terminate - } - acc1 += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6]; - acc2 += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6]; - jrk1 += adot_pn1[1] + adot_pn1[2] + adot_pn1[3] + adot_pn1[4] + adot_pn1[5] + adot_pn1[6]; - jrk2 += adot_pn2[1] + adot_pn2[2] + adot_pn2[3] + adot_pn2[4] + adot_pn2[5] + adot_pn2[6]; -} - -void write_bh_data(double time_cur, double m[], double3 x[], double3 v[], double pot[], double3 a[], double3 adot[], double dt[]) -{ - if (config->live_smbh_count == 2) { - - auto out = fopen("bh.dat", "a"); - - for (int i=0; i < 2; i++) { - double3 *a_pn, *adot_pn, a_bh, adot_bh; - double pot_bh; - if (i==0) { - a_pn = a_pn1; - adot_pn = adot_pn1; - pot_bh = pot_bh1; - a_bh = a_bh1; - adot_bh = adot_bh1; - } else { - a_pn = a_pn2; - adot_pn = adot_pn2; - pot_bh = pot_bh2; - a_bh = a_bh2; - adot_bh = adot_bh2; - } - if (config->live_smbh_custom_eps >= 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 (config->binary_smbh_pn) { - 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"); - fclose(out); - } else if (config->live_smbh_count == 1) { - auto out = fopen("bh.dat", "a"); - double tmp_r = sqrt( SQR(x[0][0]) + SQR(x[0][1]) + SQR(x[0][2]) ); - double tmp_v = sqrt( SQR(v[0][0]) + SQR(v[0][1]) + SQR(v[0][2]) ); - double tmp_a = sqrt( SQR(a[0][0]) + SQR(a[0][1]) + SQR(a[0][2]) ); - double tmp_adot = sqrt( SQR(adot[0][0]) + SQR(adot[0][1]) + SQR(adot[0][2]) ); - 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], tmp_r, - v[0][0], v[0][1], v[0][2], tmp_v, - pot[0], - a[0][0], a[0][1], a[0][2], tmp_a, - adot[0][0], adot[0][1], adot[0][2], tmp_adot, - dt[0]); - fprintf(out,"\n"); - fclose(out); - } -} - -void write_bh_nb_data(double time_cur, int N, double m[], double3 x[], double3 v[]) -{ - int nb = config->live_smbh_neighbor_number; - int ind_sort[N_MAX]; - double var_sort[N_MAX]; - //TODO you don't want and probably don't need to allocate these here. Maybe just need size nb, or maybe use as private variables. - - auto out = fopen("bh_neighbors.dat", "a"); - - /* 1st BH */ - - for (int i_bh=0; i_bh < config->live_smbh_count; i_bh++) { - for (int i=0; i BH softened pot, acc & jerk - - calc_force_n_BH(m_bh1, x_bh1, v_bh1, - m_bh2, x_bh2, v_bh2, - eps, - &pot_bh1, a_bh1, adot_bh1, - &pot_bh2, a_bh2, adot_bh2); - - pot_act_new[i_bh1] -= pot_bh1; - pot_act_new[i_bh2] -= pot_bh2; - - a_act_new[i_bh1] -= a_bh1; - a_act_new[i_bh2] -= a_bh2; - - adot_act_new[i_bh1] -= adot_bh1; - adot_act_new[i_bh2] -= adot_bh2; - - // calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk - - calc_force_n_BH(m_bh1, x_bh1, v_bh1, - m_bh2, x_bh2, v_bh2, - eps_bh, - &pot_bh1, a_bh1, adot_bh1, - &pot_bh2, a_bh2, adot_bh2); - - pot_act_new[i_bh1] += pot_bh1; - pot_act_new[i_bh2] += pot_bh2; - - a_act_new[i_bh1] += a_bh1; - a_act_new[i_bh2] += a_bh2; - - adot_act_new[i_bh1] += adot_bh1; - adot_act_new[i_bh2] += adot_bh2; -} - int main(int argc, char *argv[]) { int skip_con=0; // skip_con should just be a static in the energy control function @@ -1030,6 +746,7 @@ int main(int argc, char *argv[]) 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.set_softening(eps, config->live_smbh_custom_eps); #ifdef ETICS grapite_read_particle_tags(N, config->grapite_mask_file_name.c_str(), myRank, n_loc); @@ -1075,35 +792,9 @@ int main(int argc, char *argv[]) MPI_Barrier(MPI_COMM_WORLD); if (config->live_smbh_count == 2) { - i_bh1 = 0; - i_bh2 = 1; - - if (config->live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(eps, config->live_smbh_custom_eps, x[0], x[1], v[0], v[1], pot[0], pot[1], a[0], a[1], adot[0], adot[1]); - - if (config->binary_smbh_pn) { - - // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk - int tmp; - tmp = calc_force_pn_BH(m[i_bh1], x[i_bh1], v[i_bh1], config->smbh1_spin.data(), - m[i_bh2], x[i_bh2], v[i_bh2], config->smbh2_spin.data(), - config->pn_c, dt[i_bh1], config->pn_usage.data(), - (double(*)[3])a_pn1, (double(*)[3])adot_pn1, - (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank); - - a[i_bh1] += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6]; - a[i_bh2] += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6]; - - adot[i_bh1] += adot_pn1[1] + adot_pn1[2] + adot_pn1[3] + adot_pn1[4] + adot_pn1[5] + adot_pn1[6]; - adot[i_bh2] += adot_pn2[1] + adot_pn2[2] + adot_pn2[3] + adot_pn2[4] + adot_pn2[5] + adot_pn2[6]; - - if (myRank == rootRank) { - if (tmp == 505) { - printf("PN RSDIST: %.8E \t %.8E \n", timesteps, time_cur); - fflush(stdout); - exit(-1); - } - } - } + black_hole_physics.set_xv(x[0], x[1], v[0], v[1]); + if (config->live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(pot[0], pot[1], a[0], a[1], adot[0], adot[1]); + if (config->binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt[0], a[0], a[1], adot[0], adot[1]); } calc_ext_grav(external_gravity_components, N, x, v, pot_ext, a, adot); @@ -1111,8 +802,6 @@ int main(int argc, char *argv[]) /* Wait to all processors to finish his works... */ MPI_Barrier(MPI_COMM_WORLD); - /* Wait to all processors to finish his works... */ - MPI_Barrier(MPI_COMM_WORLD); /* Energy control... */ if (myRank == rootRank) { @@ -1183,10 +872,10 @@ int main(int argc, char *argv[]) if (myRank == rootRank) { /* Write BH data... */ - if (config->live_smbh_output) 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, m, x, v, pot, a, adot, dt); /* Write BH NB data... */ - if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur, N, m, x, v); + if (config->live_smbh_neighbor_output) write_bh_nb_data(config->live_smbh_neighbor_number, config->live_smbh_count, time_cur, N, m, x, v); } /* if (myRank == rootRank) */ @@ -1361,39 +1050,9 @@ int main(int argc, char *argv[]) #endif if (config->live_smbh_count == 2) { - //if (config->live_smbh_custom_eps >= 0) adjust_bsmbh_softening(eps, config->live_smbh_custom_eps, i_bh1, i_bh2, m[0], m[1], x_act_new, v_act_new, pot_act_new, a_act_new, adot_act_new); - if (config->live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(eps, config->live_smbh_custom_eps, x_act_new[i_bh1], x_act_new[i_bh2], v_act_new[i_bh1], v_act_new[i_bh2], 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]); - - //TODO the below works but it interferes with the printing - /* if (config->binary_smbh_pn) black_hole_physics.adjust_post_newtonian( - dt[i_bh1], - x_act_new[i_bh1], x_act_new[i_bh2], - v_act_new[i_bh1], v_act_new[i_bh2], -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) { - // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk - int tmp; - tmp = calc_force_pn_BH(m[0], x_act_new[i_bh1], v_act_new[i_bh1], config->smbh1_spin.data(), - m[1], x_act_new[i_bh2], v_act_new[i_bh2], config->smbh2_spin.data(), - config->pn_c, dt[i_bh1], config->pn_usage.data(), - (double(*)[3])a_pn1, (double(*)[3])adot_pn1, - (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank); - - a_act_new[i_bh1] += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6]; - a_act_new[i_bh2] += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6]; - - adot_act_new[i_bh1] += adot_pn1[1] + adot_pn1[2] + adot_pn1[3] + adot_pn1[4] + adot_pn1[5] + adot_pn1[6]; - adot_act_new[i_bh2] += adot_pn2[1] + adot_pn2[2] + adot_pn2[3] + adot_pn2[4] + adot_pn2[5] + adot_pn2[6]; - - if (myRank == rootRank) { - if (tmp == 505) { - printf("PN RSDIST: TS = %.8E \t t = %.8E \n", timesteps, time_cur); - fflush(stdout); - exit(-1); - } - } - } + 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]); } calc_ext_grav(external_gravity_components, n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new); @@ -1508,10 +1167,10 @@ pot_act_new[i_bh1], pot_act_new[i_bh2], a_act_new[i_bh1], a_act_new[i_bh2], adot if (time_cur >= t_bh) { if (myRank == rootRank) { /* Write BH data... */ - if (config->live_smbh_output) 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, m, x, v, pot, a, adot, dt); /* Write BH NB data... */ - if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur, N, m, x, v); + if (config->live_smbh_neighbor_output) write_bh_nb_data(config->live_smbh_neighbor_number, config->live_smbh_count, time_cur, N, m, x, v); } /* if (myRank == rootRank) */