diff --git a/black_holes.cpp b/black_holes.cpp index 48685a6..0db9a13 100644 --- a/black_holes.cpp +++ b/black_holes.cpp @@ -2,12 +2,13 @@ #include #include "black_holes.h" - -// TODO do not include c files! -#define SQR(x) ((x)*(x)) +/* BEGIN legacy inclusion */ +// I'm not going to touch this C file +#define SQR(x) ((x)*(x)) double L[3]; // needed in pn_bh_spin.c -#include "n_bh.c" #include "pn_bh_spin.c" +#undef SQR +/* END legacy inclusion */ void two_body_gravity( const double m1, const double3 x1, const double3 v1, @@ -32,55 +33,6 @@ void two_body_gravity( 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, @@ -215,22 +167,12 @@ void Black_hole_physics::write_bh_data(double time_cur, double m[], double3 x[], fclose(out); } -void write_bh_nb_data(int nb, int smbh_count, double time_cur, int N, double m[], double3 x[], double3 v[]) +void Write_bh_nb_data::operator()(double time_cur) { - //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 +#include #include "double3.h" struct Bbh_gravity { @@ -56,4 +58,49 @@ public: //TODO make private 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[]); +//void write_bh_nb_data(int nb, int smbh_count, double time_cur, int N, double m[], double3 x[], double3 v[]); + +class Write_bh_nb_data { +public: + Write_bh_nb_data(int nb, int smbh_count, int N, double *m, double3 *x, double3 *v) + : nb(nb), smbh_count(smbh_count), N(N), m(m), x(x), v(v) + { + ind_sort.resize(N); + var_sort.resize(N); + out = fopen("bh_neighbors.dat", "w"); + } + ~Write_bh_nb_data() + { + fclose(out); + } + void operator()(double time_cur); +private: + int nb, smbh_count, N; + double *m; + double3 *x, *v; + std::vector ind_sort; + std::vector var_sort; + FILE *out; +}; + +class Binary_smbh_influence_sphere_output { +public: + Binary_smbh_influence_sphere_output(double factor, int N, double *m, double3 *x, double3 *v, double *pot, double *dt) + : factor(factor), m(m), x(x), v(v), pot(pot), dt(dt) + { + inf_event.assign(N, 0); + out = fopen("bbh_inf.dat", "w"); + } + + ~Binary_smbh_influence_sphere_output() + { + fclose(out); + } + void operator()(int ind_act[], int n_act, double timesteps, double time_cur); +private: + double factor; + double *m, *pot, *dt; + double3 *x, *v; + std::vector inf_event; + FILE *out; +}; \ No newline at end of file diff --git a/n_bh.c b/n_bh.c deleted file mode 100644 index 67223d9..0000000 --- a/n_bh.c +++ /dev/null @@ -1,74 +0,0 @@ -/***************************************************************************/ -/* - Coded by : Peter Berczik - Version number : 1.0 - Last redaction : 2011.II.20. 13:00 -*/ - -void calc_force_n_BH(double m1, double xx1[], double vv1[], - double m2, double xx2[], double vv2[], - double eps_BH, - double *pot_n1, double a_n1[], double adot_n1[], - double *pot_n2, double a_n2[], double adot_n2[]) -{ - -/* - INPUT - -m1 - mass of the 1 BH -xx1[0,1,2] - coordinate of the 1 BH -vv1[0,1,2] - velocity of the 1 BH - -m2 - mass of the 2 BH -xx2[0,1,2] - coordinate of the 2 BH -vv2[0,1,2] - velocity of the 2 BH - -eps_BH - force softening, can be even exactly 0.0 ! - - OUTPUT - -pot_n1 for the 1 BH -a_n1 [0,1,2] for the 1 BH -adot_n1 [0,1,2] for the 1 BH - -pot_n2 for the 2 BH -a_n2 [0,1,2] for the 2 BH -adot_n2 [0,1,2] for the 2 BH -*/ - - -int k; - -double r, r2, r3, r4, RP, x[3], v[3]; - - -for(k=0;k<3;k++) - { - x[k] = xx1[k] - xx2[k]; - v[k] = vv1[k] - vv2[k]; - } - -r2 = SQR(x[0]) + SQR(x[1]) + SQR(x[2]) + SQR(eps_BH); - -r = sqrt(r2); -r3 = r2*r; -r4 = r2*r2; - -RP = 3.0*(x[0]*v[0]+x[1]*v[1]+x[2]*v[2])/r; - -// Newton pot + acc & jerks - -*pot_n1 = -m2/r; -*pot_n2 = -m1/r; - -for(k=0;k<3;k++) - { - a_n1[k] = -m2*x[k]/r3; - a_n2[k] = m1*x[k]/r3; - - adot_n1[k] = -m2*(v[k]/r3 - RP*x[k]/r4); - adot_n2[k] = m1*(v[k]/r3 - RP*x[k]/r4); - } - -} -/***************************************************************************/ diff --git a/phigrape.cpp b/phigrape.cpp index 86d8121..da5b8d9 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -124,20 +124,6 @@ double DT_ACT_REDUCE; #endif - -/* external potential... */ - -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" - -//double3 a_pn1[7], adot_pn1[7], a_pn2[7], adot_pn2[7]; - -//#include "pn_bh_spin.c" - #ifdef ETICS int grapite_cep_index; #endif @@ -452,76 +438,6 @@ inline double aarseth_step(const double eta, const double dt_tmp, const double3 return sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); } -void binary_smbh_influence_sphere_output(int i_bh1, int i_bh2, int ind_act[], int n_act, double timesteps, double time_cur, double factor, int inf_event[], Source_particle_list *source_particle_list) -{ - //TODO !!IMPORTANT!! only open the file IF THERE IS ANYTHING TO WRITE!!! - //TODO inf_event to be static or something? - auto out = fopen("bbh_inf.dat", "a"); - - double m_bh1 = source_particle_list->m[0]; - double m_bh2 = source_particle_list->m[1]; - double3 x_bh1 = source_particle_list->x[0]; - double3 x_bh2 = source_particle_list->x[1]; - double3 v_bh1 = source_particle_list->v[0]; - double3 v_bh2 = source_particle_list->v[1]; - - double3 x_bbhc = (m_bh1*x_bh1 + m_bh2*x_bh2)/(m_bh1 + m_bh2); - double3 v_bbhc = (m_bh1*v_bh1 + m_bh2*v_bh2)/(m_bh1 + m_bh2); - - double DR2 = (x_bh1 - x_bh2).norm2(); - double DV2 = (v_bh1 - v_bh2).norm2(); - double EB = -(m_bh1 + m_bh2) / sqrt(DR2) + 0.5 * DV2; - double SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB; - double SEMI_a2 = SQR(SEMI_a); - - for (int i=0; ipot[0]; - double& pot_bh2 = source_particle_list->pot[1]; - double& m_act = source_particle_list->m[j_act]; - double3& x_act = source_particle_list->x[j_act]; - double3& v_act = source_particle_list->v[j_act]; - double& dt_act = source_particle_list->dt[j_act]; - double& pot_act = source_particle_list->pot[j_act]; - double tmp_r2 = (x_act - x_bbhc).norm2(); - - if (tmp_r2 < SEMI_a2*SQR(factor)) { - - if (inf_event[j_act] == 0) { - - fprintf(out,"INF1 %.6E %.16E %07d %07d %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E \n", - timesteps, time_cur, i, j_act, - sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2], - m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_bh1, - m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_bh2, - sqrt(tmp_r2), - m_act, x_act[0], x_act[1], x_act[2], v_act[0], v_act[1], v_act[2], pot_act, - dt_act); - - inf_event[j_act] = 1; - } - - } else { - if (inf_event[j_act] == 1) { - fprintf(out,"INF2 %.6E %.16E %07d %07d %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E \n", - timesteps, time_cur, i, ind_act[i], - sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2], - m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_bh1, - m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_bh2, - sqrt(tmp_r2), - m_act, x_act[0], x_act[1], x_act[2], v_act[0], v_act[1], v_act[2], pot_act, - dt_act); - } - - inf_event[j_act] = 0; - - } /* if (tmp_r2 < DR2*R_INF2) */ - - } /* i */ - fclose(out); -} int main(int argc, char *argv[]) { @@ -617,8 +533,6 @@ int main(int argc, char *argv[]) double dt_bh = config->dt_bh; if (myRank == rootRank) { - if (config->binary_smbh_influence_sphere_output) for (int i=0; ibinary_smbh_influence_sphere_output) { - out = fopen("bbh_inf.dat", "w"); - fclose(out); - } + // if (config->binary_smbh_influence_sphere_output) { + // out = fopen("bbh_inf.dat", "w"); + // fclose(out); + // } } get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0); @@ -748,6 +662,10 @@ int main(int argc, char *argv[]) } black_hole_physics.set_softening(eps, config->live_smbh_custom_eps); + 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); + #ifdef ETICS grapite_read_particle_tags(N, config->grapite_mask_file_name.c_str(), myRank, n_loc); grapite_set_dt_exp(config->dt_scf); @@ -875,7 +793,7 @@ int main(int argc, char *argv[]) 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(config->live_smbh_neighbor_number, config->live_smbh_count, time_cur, N, m, x, v); + if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur); } /* if (myRank == rootRank) */ @@ -1122,7 +1040,7 @@ int main(int argc, char *argv[]) if (config->binary_smbh_influence_sphere_output && (myRank == rootRank)) { //TODO clean up this mass. We don't want to have all these _act arrays; they are only needed for this lame function. - binary_smbh_influence_sphere_output(i_bh1, i_bh2, ind_act, n_act, timesteps, time_cur, config->binary_smbh_influence_radius_factor, inf_event, &source_particle_list); + binary_smbh_influence_sphere_output(ind_act, n_act, timesteps, time_cur); } #ifdef TIMING @@ -1170,7 +1088,7 @@ int main(int argc, char *argv[]) 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(config->live_smbh_neighbor_number, config->live_smbh_count, time_cur, N, m, x, v); + if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur); } /* if (myRank == rootRank) */