diff --git a/phigrape.cpp b/phigrape.cpp index 905655e..6c1e805 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -54,23 +54,36 @@ Version number : 19.04 Last redaction : 2019.04.16 12:55 *****************************************************************************/ +struct Parameters { + int live_smbh_count = 2; // ADD_BH1 or ADD_BH2 + double live_smbh_custom_eps = 0; // ADD_N_BH + bool live_smbh_output = true; // BH_OUT + bool live_smbh_neighbor_output = true; // BH_OUT_NB // NOTE needs `live_smbh_neighbor_number` + bool binary_smbh_influence_sphere_output = true; // BBH_INF // NOTE needs `binary_smbh_influence_radius_factor` + double binary_smbh_influence_radius_factor = 3.162277660168379497918067e+03; // R_INF + int live_smbh_neighbor_number = 10; // TODO make sure it's smaller than N? or just warn if it's not + bool binary_smbh_pn = true; // ADD_PN_BH + bool binary_smbh_spin = true; // ADD_SPIN_BH + bool dt_min_warning = true; // DT_MIN_WARNING +} parameters; + //#define NORM // Physical normalization -//#define ADD_BH1 // add the Single BH +// // // //#define ADD_BH1 // add the Single BH -//#define ADD_BH2 // add the Binary BH's -//#define ADD_N_BH // eps_BH = 0.0, but added only the Newtonian forces -//#define ADD_PN_BH // extra - added also the Post-Newton forces -//#define ADD_SPIN_BH // extra - added the SPIN for the BH's - DEFAULT !!! +// // // // #define ADD_BH2 // add the Binary BH's +// // // // #define ADD_N_BH // eps_BH = 0.0, but added only the Newtonian forces +// // // // #define ADD_PN_BH // extra - added also the Post-Newton forces +// // // // #define ADD_SPIN_BH // extra - added the SPIN for the BH's - DEFAULT !!! -//#define BH_OUT // extra output for BH's (live) -//#define BH_OUT_NB // extra output for the BH's neighbours (live) +// // // // #define BH_OUT // extra output for BH's (live) +// // // // #define BH_OUT_NB // extra output for the BH's neighbours (live) -//#define BBH_INF // BBH influence sphere... -//#define R_INF 10.0 // Factor for the influence sphere... if ( R < R_INF * DR_BBH ) -//#define R_INF2 (R_INF*R_INF) +// // // // #define BBH_INF // BBH influence sphere... +// // // // #define R_INF 3.162277660168379497918067e+03 // Factor for the influence sphere... if ( R < R_INF * DR_BBH ) +// // // // #define R_INF2 (R_INF*R_INF) -//#define DT_MIN_WARNING // dt < dt_min warning !!! +// // // // #define DT_MIN_WARNING // dt < dt_min warning !!! //#define BH_OUT_NB_EXT // extra output for the BH's neighbours (external) @@ -146,13 +159,7 @@ Last redaction : 2019.04.16 12:55 #include #include - -/* -double aaa; -double aaapars[5]; - -extern void qwerty_(double *aaa, double *aaapars); -*/ +#include #define G6_NPIPE 2048 #include "grape6.h" @@ -244,6 +251,13 @@ struct double3 { data[2] += a.data[2]; return *this; } + double3& operator-=(const double3& a) + { + data[0] -= a.data[0]; + data[1] -= a.data[1]; + data[2] -= a.data[2]; + return *this; + } double3& operator/=(const double& c) { data[0] /= c; @@ -268,6 +282,11 @@ double3 operator*(const double3& a, const double& c) return double3(a.data[0]*c, a.data[1]*c, a.data[2]*c); } +double3 operator/(const double3& a, const double& c) +{ + return double3(a.data[0]/c, a.data[1]/c, a.data[2]/c); +} + double3 operator+(const double3& a, const double3& b) { return double3(a.data[0]+b.data[0], a.data[1]+b.data[1], a.data[2]+b.data[2]); @@ -306,13 +325,21 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, tmp_a_bh_pn0, tmp_a_bh_pn1, tmp_a_bh_pn2, tmp_a_bh_pn2_5, tmp_a_bh_pn3, tmp_a_bh_pn3_5, tmp_a_bh_spin, tmp_adot_bh_pn0, tmp_adot_bh_pn1, tmp_adot_bh_pn2, tmp_adot_bh_pn2_5, tmp_adot_bh_pn3, tmp_adot_bh_pn3_5, tmp_adot_bh_spin; + +double L[3]; // needed in pn_bh_spin.c +// Needed for things related to BHs +#include "debug.h" +int ind_sort[N_MAX]; +double var_sort[N_MAX]; + + double3 xcm, vcm, mom, xdc, vdc, a2, a3, a2dot1; char processor_name[MPI_MAX_PROCESSOR_NAME], inp_fname[30], - out_fname[30], dbg_fname[30]; + out_fname[30]; /* global variables */ @@ -344,7 +371,7 @@ double3 x_act[N_MAX], v_act[N_MAX], a_act_tmp[N_MAX], adot_act_tmp[N_MAX], a_act_tmp_loc[N_MAX], adot_act_tmp_loc[N_MAX]; -FILE *inp, *out, *tmp_file, *dbg; +FILE *inp, *out; double CPU_time_real0, CPU_time_user0, CPU_time_syst0; double CPU_time_real, CPU_time_user, CPU_time_syst; @@ -375,8 +402,6 @@ double3 x_i[G6_NPIPE], v_i[G6_NPIPE], a_i[G6_NPIPE], jerk_i[G6_NPIPE]; int new_tunit=51, new_xunit=51; -int aflag=1, jflag=1, pflag=1; - double ti=0.0; double3 a2by18, a1by6, aby2; @@ -395,8 +420,8 @@ double eps_BH=0.0; #ifdef EXTPOT_GAL double m_bulge, a_bulge, b_bulge, m_disk, a_disk, b_disk, - m_halo, a_halo, b_halo, - x2_ij, y2_ij, z2_ij, + m_halo, a_halo, b_halo, // NOTE there are other "halo" variables in EXTPOT_GAL_LOG + x2_ij, y2_ij, z2_ij, // NOTE this variable exist also in EXTPOT_GAL_LOG r_tmp, r2_tmp, z_tmp, z2_tmp; #endif @@ -413,7 +438,9 @@ double v_halo, r_halo, #ifdef EXTPOT_BH double r2, rv_ij, - m_bh, b_bh, eps_bh; + m_bh, b_bh, eps_bh; // NOTE different from eps_BH + // NOTE there are other m_bh and b_bh defined outside this ifdef block + // NOTE this is a mess. Looks like eps_bh is never used, and the m_bh from outside the block is never used. #endif double pot_ext[N_MAX], pot_act_ext[N_MAX]; // for all EXTPOT @@ -428,18 +455,23 @@ double m_bh, m_bh1, m_bh2, b_bh, x_ij, y_ij, z_ij, vx_ij, vy_ij, vz_ij, rv_ij; -#ifdef BBH_INF +// #ifdef BBH_INF int inf_event[N_MAX]; -double x_bbhc[3], v_bbhc[3], DR2, tmp_r2; +double DR2, tmp_r2; +double3 x_bbhc, v_bbhc; double DV2, EB, SEMI_a, SEMI_a2; -#endif +// #endif // BBH_INF -#ifdef ADD_N_BH +// #ifdef ADD_N_BH -double x_bh1[3], x_bh2[3], v_bh1[3], v_bh2[3]; +// double x_bh1[3], x_bh2[3], v_bh1[3], v_bh2[3]; +double3 x_bh1, x_bh2, v_bh1, v_bh2; + +// double pot_bh1, a_bh1[3], adot_bh1[3], +// pot_bh2, a_bh2[3], adot_bh2[3]; +double pot_bh1, pot_bh2; +double3 a_bh1, adot_bh1, a_bh2, adot_bh2; -double pot_bh1, a_bh1[3], adot_bh1[3], - pot_bh2, a_bh2[3], adot_bh2[3]; //double eps_BH = 0.0; @@ -479,16 +511,17 @@ adot_n2 [0,1,2] for the 2 BH return - 0 if everything OK */ -#endif +// #endif // ADD_N_BH -#ifdef ADD_PN_BH +// #ifdef ADD_PN_BH double C_NB = 477.12; int usedOrNot[7] = {1, 1, 1, 1, 0, 0, 0}; -double a_pn1[7][3], adot_pn1[7][3], - a_pn2[7][3], adot_pn2[7][3]; +// double a_pn1[7][3], adot_pn1[7][3], +// a_pn2[7][3], adot_pn2[7][3]; +double3 a_pn1[7], adot_pn1[7], a_pn2[7], adot_pn2[7]; double s_bh1[3] = {0.0, 0.0, 1.0}; double s_bh2[3] = {0.0, 0.0, 1.0}; @@ -535,7 +568,7 @@ return - 0 if everything OK - 505 if BH's separation < 4 x (RSwarch1 + RSwarch2) */ -#endif +// #endif // ADD_N_BH #ifdef ETICS double t_exp, dt_exp=ETICS_DTSCF; // t_exp is just the initial value @@ -624,386 +657,164 @@ fclose(out); void write_bh_data() { + if (parameters.live_smbh_count == 2) { -#ifdef ADD_BH2 + out = fopen("bh.dat","a"); -out = fopen("bh.dat","a"); + for (int i=0; i < 2; i++) { +// double (*a_pn)[3], (*adot_pn)[3], pot_bh, *a_bh, *adot_bh; + 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; + } -/* 1st BH */ + tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) ); + tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) ); -// i=N-2; -i=0; + tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) ); + tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) ); -tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) ); -tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) ); + if (parameters.live_smbh_custom_eps >= 0) { -tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) ); -tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) ); + tmp_a_bh = sqrt( SQR(a_bh[0]) + SQR(a_bh[1]) + SQR(a_bh[2]) ); + tmp_adot_bh = sqrt( SQR(adot_bh[0]) + SQR(adot_bh[1]) + SQR(adot_bh[2]) ); -#ifdef ADD_N_BH + if (parameters.binary_smbh_spin) { + tmp_a_bh_pn0 = sqrt( SQR(a_pn[0][0]) + SQR(a_pn[0][1]) + SQR(a_pn[0][2]) ); + tmp_a_bh_pn1 = sqrt( SQR(a_pn[1][0]) + SQR(a_pn[1][1]) + SQR(a_pn[1][2]) ); + tmp_a_bh_pn2 = sqrt( SQR(a_pn[2][0]) + SQR(a_pn[2][1]) + SQR(a_pn[2][2]) ); + tmp_a_bh_pn2_5 = sqrt( SQR(a_pn[3][0]) + SQR(a_pn[3][1]) + SQR(a_pn[3][2]) ); + tmp_a_bh_pn3 = sqrt( SQR(a_pn[4][0]) + SQR(a_pn[4][1]) + SQR(a_pn[4][2]) ); + tmp_a_bh_pn3_5 = sqrt( SQR(a_pn[5][0]) + SQR(a_pn[5][1]) + SQR(a_pn[5][2]) ); + tmp_a_bh_spin = sqrt( SQR(a_pn[6][0]) + SQR(a_pn[6][1]) + SQR(a_pn[6][2]) ); -tmp_a_bh = sqrt( SQR(a_bh1[0]) + SQR(a_bh1[1]) + SQR(a_bh1[2]) ); -tmp_adot_bh = sqrt( SQR(adot_bh1[0]) + SQR(adot_bh1[1]) + SQR(adot_bh1[2]) ); - -#ifdef ADD_PN_BH - -tmp_a_bh_pn0 = sqrt( SQR(a_pn1[0][0]) + SQR(a_pn1[0][1]) + SQR(a_pn1[0][2]) ); -tmp_a_bh_pn1 = sqrt( SQR(a_pn1[1][0]) + SQR(a_pn1[1][1]) + SQR(a_pn1[1][2]) ); -tmp_a_bh_pn2 = sqrt( SQR(a_pn1[2][0]) + SQR(a_pn1[2][1]) + SQR(a_pn1[2][2]) ); -tmp_a_bh_pn2_5 = sqrt( SQR(a_pn1[3][0]) + SQR(a_pn1[3][1]) + SQR(a_pn1[3][2]) ); -tmp_a_bh_pn3 = sqrt( SQR(a_pn1[4][0]) + SQR(a_pn1[4][1]) + SQR(a_pn1[4][2]) ); -tmp_a_bh_pn3_5 = sqrt( SQR(a_pn1[5][0]) + SQR(a_pn1[5][1]) + SQR(a_pn1[5][2]) ); -tmp_a_bh_spin = sqrt( SQR(a_pn1[6][0]) + SQR(a_pn1[6][1]) + SQR(a_pn1[6][2]) ); - -tmp_adot_bh_pn0 = sqrt( SQR(adot_pn1[0][0]) + SQR(adot_pn1[0][1]) + SQR(adot_pn1[0][2]) ); -tmp_adot_bh_pn1 = sqrt( SQR(adot_pn1[1][0]) + SQR(adot_pn1[1][1]) + SQR(adot_pn1[1][2]) ); -tmp_adot_bh_pn2 = sqrt( SQR(adot_pn1[2][0]) + SQR(adot_pn1[2][1]) + SQR(adot_pn1[2][2]) ); -tmp_adot_bh_pn2_5 = sqrt( SQR(adot_pn1[3][0]) + SQR(adot_pn1[3][1]) + SQR(adot_pn1[3][2]) ); -tmp_adot_bh_pn3 = sqrt( SQR(adot_pn1[4][0]) + SQR(adot_pn1[4][1]) + SQR(adot_pn1[4][2]) ); -tmp_adot_bh_pn3_5 = sqrt( SQR(adot_pn1[5][0]) + SQR(adot_pn1[5][1]) + SQR(adot_pn1[5][2]) ); -tmp_adot_bh_spin = sqrt( SQR(adot_pn1[6][0]) + SQR(adot_pn1[6][1]) + SQR(adot_pn1[6][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 \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t\t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \n", - time_cur, m[i], - x[i][0], x[i][1], x[i][2], tmp_r, - v[i][0], v[i][1], v[i][2], tmp_v, - pot[i], - a[i][0], a[i][1], a[i][2], tmp_a, - adot[i][0], adot[i][1], adot[i][2], tmp_adot, - dt[i], - pot_bh1, - a_bh1[0], a_bh1[1], a_bh1[2], tmp_a_bh, - adot_bh1[0], adot_bh1[1], adot_bh1[2], tmp_adot_bh, - a_pn1[0][0], a_pn1[0][1], a_pn1[0][2], tmp_a_bh_pn0, - adot_pn1[0][0], adot_pn1[0][1], adot_pn1[0][2], tmp_adot_bh_pn0, - a_pn1[1][0], a_pn1[1][1], a_pn1[1][2], tmp_a_bh_pn1, - adot_pn1[1][0], adot_pn1[1][1], adot_pn1[1][2], tmp_adot_bh_pn1, - a_pn1[2][0], a_pn1[2][1], a_pn1[2][2], tmp_a_bh_pn2, - adot_pn1[2][0], adot_pn1[2][1], adot_pn1[2][2], tmp_adot_bh_pn2, - a_pn1[3][0], a_pn1[3][1], a_pn1[3][2], tmp_a_bh_pn2_5, - adot_pn1[3][0], adot_pn1[3][1], adot_pn1[3][2], tmp_adot_bh_pn2_5, - a_pn1[4][0], a_pn1[4][1], a_pn1[4][2], tmp_a_bh_pn3, - adot_pn1[4][0], adot_pn1[4][1], adot_pn1[4][2], tmp_adot_bh_pn3, - a_pn1[5][0], a_pn1[5][1], a_pn1[5][2], tmp_a_bh_pn3_5, - adot_pn1[5][0], adot_pn1[5][1], adot_pn1[5][2], tmp_adot_bh_pn3_5, - a_pn1[6][0], a_pn1[6][1], a_pn1[6][2], tmp_a_bh_spin, - adot_pn1[6][0], adot_pn1[6][1], adot_pn1[6][2], tmp_adot_bh_spin); - -#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 \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \n", - time_cur, m[i], - x[i][0], x[i][1], x[i][2], tmp_r, - v[i][0], v[i][1], v[i][2], tmp_v, - pot[i], - a[i][0], a[i][1], a[i][2], tmp_a, - adot[i][0], adot[i][1], adot[i][2], tmp_adot, - dt[i], - pot_bh1, - a_bh1[0], a_bh1[1], a_bh1[2], tmp_a_bh, - adot_bh1[0], adot_bh1[1], adot_bh1[2], tmp_adot_bh); - -#endif // ADD_PN_BH - -#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], tmp_r, - v[i][0], v[i][1], v[i][2], tmp_v, - pot[i], - a[i][0], a[i][1], a[i][2], tmp_a, - adot[i][0], adot[i][1], adot[i][2], tmp_adot, - dt[i]); - -#endif // ADD_N_BH - -/* 2nd BH */ - -// i=N-1; -i=1; - -tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) ); -tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) ); - -tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) ); -tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) ); - -#ifdef ADD_N_BH - -tmp_a_bh = sqrt( SQR(a_bh2[0]) + SQR(a_bh2[1]) + SQR(a_bh2[2]) ); -tmp_adot_bh = sqrt( SQR(adot_bh2[0]) + SQR(adot_bh2[1]) + SQR(adot_bh2[2]) ); - -#ifdef ADD_PN_BH - -tmp_a_bh_pn0 = sqrt( SQR(a_pn2[0][0]) + SQR(a_pn2[0][1]) + SQR(a_pn2[0][2]) ); -tmp_a_bh_pn1 = sqrt( SQR(a_pn2[1][0]) + SQR(a_pn2[1][1]) + SQR(a_pn2[1][2]) ); -tmp_a_bh_pn2 = sqrt( SQR(a_pn2[2][0]) + SQR(a_pn2[2][1]) + SQR(a_pn2[2][2]) ); -tmp_a_bh_pn2_5 = sqrt( SQR(a_pn2[3][0]) + SQR(a_pn2[3][1]) + SQR(a_pn2[3][2]) ); -tmp_a_bh_pn3 = sqrt( SQR(a_pn2[4][0]) + SQR(a_pn2[4][1]) + SQR(a_pn2[4][2]) ); -tmp_a_bh_pn3_5 = sqrt( SQR(a_pn2[5][0]) + SQR(a_pn2[5][1]) + SQR(a_pn2[5][2]) ); -tmp_a_bh_spin = sqrt( SQR(a_pn2[6][0]) + SQR(a_pn2[6][1]) + SQR(a_pn2[6][2]) ); - -tmp_adot_bh_pn0 = sqrt( SQR(adot_pn2[0][0]) + SQR(adot_pn2[0][1]) + SQR(adot_pn2[0][2]) ); -tmp_adot_bh_pn1 = sqrt( SQR(adot_pn2[1][0]) + SQR(adot_pn2[1][1]) + SQR(adot_pn2[1][2]) ); -tmp_adot_bh_pn2 = sqrt( SQR(adot_pn2[2][0]) + SQR(adot_pn2[2][1]) + SQR(adot_pn2[2][2]) ); -tmp_adot_bh_pn2_5 = sqrt( SQR(adot_pn2[3][0]) + SQR(adot_pn2[3][1]) + SQR(adot_pn2[3][2]) ); -tmp_adot_bh_pn3 = sqrt( SQR(adot_pn2[4][0]) + SQR(adot_pn2[4][1]) + SQR(adot_pn2[4][2]) ); -tmp_adot_bh_pn3_5 = sqrt( SQR(adot_pn2[5][0]) + SQR(adot_pn2[5][1]) + SQR(adot_pn2[5][2]) ); -tmp_adot_bh_spin = sqrt( SQR(adot_pn2[6][0]) + SQR(adot_pn2[6][1]) + SQR(adot_pn2[6][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 \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t\t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \n", - time_cur, m[i], - x[i][0], x[i][1], x[i][2], tmp_r, - v[i][0], v[i][1], v[i][2], tmp_v, - pot[i], - a[i][0], a[i][1], a[i][2], tmp_a, - adot[i][0], adot[i][1], adot[i][2], tmp_adot, - dt[i], - pot_bh2, - a_bh2[0], a_bh2[1], a_bh2[2], tmp_a_bh, - adot_bh2[0], adot_bh2[1], adot_bh2[2], tmp_adot_bh, - a_pn2[0][0], a_pn2[0][1], a_pn2[0][2], tmp_a_bh_pn0, - adot_pn2[0][0], adot_pn2[0][1], adot_pn2[0][2], tmp_adot_bh_pn0, - a_pn2[1][0], a_pn2[1][1], a_pn2[1][2], tmp_a_bh_pn1, - adot_pn2[1][0], adot_pn2[1][1], adot_pn2[1][2], tmp_adot_bh_pn1, - a_pn2[2][0], a_pn2[2][1], a_pn2[2][2], tmp_a_bh_pn2, - adot_pn2[2][0], adot_pn2[2][1], adot_pn2[2][2], tmp_adot_bh_pn2, - a_pn2[3][0], a_pn2[3][1], a_pn2[3][2], tmp_a_bh_pn2_5, - adot_pn2[3][0], adot_pn2[3][1], adot_pn2[3][2], tmp_adot_bh_pn2_5, - a_pn2[4][0], a_pn2[4][1], a_pn2[4][2], tmp_a_bh_pn3, - adot_pn2[4][0], adot_pn2[4][1], adot_pn2[4][2], tmp_adot_bh_pn3, - a_pn2[5][0], a_pn2[5][1], a_pn2[5][2], tmp_a_bh_pn3_5, - adot_pn2[5][0], adot_pn2[5][1], adot_pn2[5][2], tmp_adot_bh_pn3_5, - a_pn2[6][0], a_pn2[6][1], a_pn2[6][2], tmp_a_bh_spin, - adot_pn2[6][0], adot_pn2[6][1], adot_pn2[6][2], tmp_adot_bh_spin); - -#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 \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \n", - time_cur, m[i], - x[i][0], x[i][1], x[i][2], tmp_r, - v[i][0], v[i][1], v[i][2], tmp_v, - pot[i], - a[i][0], a[i][1], a[i][2], tmp_a, - adot[i][0], adot[i][1], adot[i][2], tmp_adot, - dt[i], - pot_bh2, - a_bh2[0], a_bh2[1], a_bh2[2], tmp_a_bh, - adot_bh2[0], adot_bh2[1], adot_bh2[2], tmp_adot_bh); - -#endif // ADD_PN_BH - -#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], tmp_r, - v[i][0], v[i][1], v[i][2], tmp_v, - pot[i], - a[i][0], a[i][1], a[i][2], tmp_a, - adot[i][0], adot[i][1], adot[i][2], tmp_adot, - dt[i]); - -#endif // ADD_N_BH - -fprintf(out,"\n"); - -fclose(out); - -#else - -#ifdef ADD_BH1 - -out = fopen("bh.dat","a"); - -// i=N-1; -i=0; - -tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) ); -tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) ); - -tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) ); -tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][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[i], - x[i][0], x[i][1], x[i][2], tmp_r, - v[i][0], v[i][1], v[i][2], tmp_v, - pot[i], - a[i][0], a[i][1], a[i][2], tmp_a, - adot[i][0], adot[i][1], adot[i][2], tmp_adot, - dt[i]); - -fprintf(out,"\n"); - -fclose(out); - -#endif // ADD_BH1 - -#endif // ADD_BH2 + tmp_adot_bh_pn0 = sqrt( SQR(adot_pn[0][0]) + SQR(adot_pn[0][1]) + SQR(adot_pn[0][2]) ); + tmp_adot_bh_pn1 = sqrt( SQR(adot_pn[1][0]) + SQR(adot_pn[1][1]) + SQR(adot_pn[1][2]) ); + tmp_adot_bh_pn2 = sqrt( SQR(adot_pn[2][0]) + SQR(adot_pn[2][1]) + SQR(adot_pn[2][2]) ); + tmp_adot_bh_pn2_5 = sqrt( SQR(adot_pn[3][0]) + SQR(adot_pn[3][1]) + SQR(adot_pn[3][2]) ); + tmp_adot_bh_pn3 = sqrt( SQR(adot_pn[4][0]) + SQR(adot_pn[4][1]) + SQR(adot_pn[4][2]) ); + tmp_adot_bh_pn3_5 = sqrt( SQR(adot_pn[5][0]) + SQR(adot_pn[5][1]) + SQR(adot_pn[5][2]) ); + tmp_adot_bh_spin = sqrt( SQR(adot_pn[6][0]) + SQR(adot_pn[6][1]) + SQR(adot_pn[6][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 \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t\t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \n", + time_cur, m[i], + x[i][0], x[i][1], x[i][2], tmp_r, + v[i][0], v[i][1], v[i][2], tmp_v, + pot[i], + a[i][0], a[i][1], a[i][2], tmp_a, + adot[i][0], adot[i][1], adot[i][2], tmp_adot, + dt[i], + pot_bh, + a_bh[0], a_bh[1], a_bh[2], tmp_a_bh, + adot_bh[0], adot_bh[1], adot_bh[2], tmp_adot_bh, + a_pn[0][0], a_pn[0][1], a_pn[0][2], tmp_a_bh_pn0, + adot_pn[0][0], adot_pn[0][1], adot_pn[0][2], tmp_adot_bh_pn0, + a_pn[1][0], a_pn[1][1], a_pn[1][2], tmp_a_bh_pn1, + adot_pn[1][0], adot_pn[1][1], adot_pn[1][2], tmp_adot_bh_pn1, + a_pn[2][0], a_pn[2][1], a_pn[2][2], tmp_a_bh_pn2, + adot_pn[2][0], adot_pn[2][1], adot_pn[2][2], tmp_adot_bh_pn2, + a_pn[3][0], a_pn[3][1], a_pn[3][2], tmp_a_bh_pn2_5, + adot_pn[3][0], adot_pn[3][1], adot_pn[3][2], tmp_adot_bh_pn2_5, + a_pn[4][0], a_pn[4][1], a_pn[4][2], tmp_a_bh_pn3, + adot_pn[4][0], adot_pn[4][1], adot_pn[4][2], tmp_adot_bh_pn3, + a_pn[5][0], a_pn[5][1], a_pn[5][2], tmp_a_bh_pn3_5, + adot_pn[5][0], adot_pn[5][1], adot_pn[5][2], tmp_adot_bh_pn3_5, + a_pn[6][0], a_pn[6][1], a_pn[6][2], tmp_a_bh_spin, + adot_pn[6][0], adot_pn[6][1], adot_pn[6][2], tmp_adot_bh_spin); + } 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 \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \n", + time_cur, m[i], + x[i][0], x[i][1], x[i][2], tmp_r, + v[i][0], v[i][1], v[i][2], tmp_v, + pot[i], + a[i][0], a[i][1], a[i][2], tmp_a, + adot[i][0], adot[i][1], adot[i][2], tmp_adot, + dt[i], + pot_bh, + a_bh[0], a_bh[1], a_bh[2], tmp_a_bh, + adot_bh[0], adot_bh[1], adot_bh[2], tmp_adot_bh); + } + } 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], tmp_r, + v[i][0], v[i][1], v[i][2], tmp_v, + pot[i], + a[i][0], a[i][1], a[i][2], tmp_a, + adot[i][0], adot[i][1], adot[i][2], tmp_adot, + dt[i]); + } + } + fprintf(out,"\n"); + fclose(out); + } else if (parameters.live_smbh_count == 1) { + out = fopen("bh.dat","a"); + tmp_r = sqrt( SQR(x[0][0]) + SQR(x[0][1]) + SQR(x[0][2]) ); + tmp_v = sqrt( SQR(v[0][0]) + SQR(v[0][1]) + SQR(v[0][2]) ); + tmp_a = sqrt( SQR(a[0][0]) + SQR(a[0][1]) + SQR(a[0][2]) ); + 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); + } } -#ifdef BH_OUT_NB - void write_bh_nb_data() { + int i_bh, nb = parameters.live_smbh_neighbor_number; + double tmp, tmp_r, tmp_v; -#ifdef ADD_BH2 + out = fopen("bh_nb.dat", "a"); -int i_bh, nb = 10; -double tmp, tmp_r, tmp_v; + /* 1st BH */ -out = fopen("bh_nb.dat","a"); + i_bh = 0; -/* 1st BH */ + for (int i_bh=0; i_bh < parameters.live_smbh_count; i_bh++) { + for (i=0; i 0)) { + out = fopen("bh.dat", "w"); + fclose(out); + } + if ((parameters.live_smbh_neighbor_output) && (parameters.live_smbh_count > 0)) { + out = fopen("bh_nb.dat", "w"); + fclose(out); + } + if (parameters.binary_smbh_influence_sphere_output) { + out = fopen("bbh.inf","w"); + fclose(out); + } -#ifdef BH_OUT - out = fopen("bh.dat","w"); - fclose(out); -#endif - -#ifdef BH_OUT_NB - out = fopen("bh_nb.dat","w"); - fclose(out); -#endif - -#else - -#ifdef ADD_BH1 - -#ifdef BH_OUT - out = fopen("bh.dat","w"); - fclose(out); -#endif - -#ifdef BH_OUT_NB - out = fopen("bh_nb.dat","w"); - fclose(out); -#endif - -#endif // ADD_BH1 - -#endif // ADD_BH2 - -#ifdef BBH_INF - out = fopen("bbh.inf","w"); - fclose(out); -#endif } else { // if (diskstep == 0) } // if (diskstep == 0) @@ -2076,96 +1864,92 @@ int main(int argc, char *argv[]) /* Wait to all processors to finish his works... */ MPI_Barrier(MPI_COMM_WORLD); -#ifdef ADD_BH2 + if (parameters.live_smbh_count == 2) { + i_bh1 = 0; + i_bh2 = 1; - i_bh1 = 0; - i_bh2 = 1; + if (parameters.live_smbh_custom_eps >= 0) { -#ifdef ADD_N_BH + m_bh1 = m[i_bh1]; + m_bh2 = m[i_bh2]; - m_bh1 = m[i_bh1]; - m_bh2 = m[i_bh2]; + for (k=0;k<3;k++) { + x_bh1[k] = x[i_bh1][k]; + v_bh1[k] = v[i_bh1][k]; - for (k=0;k<3;k++) { - x_bh1[k] = x[i_bh1][k]; - v_bh1[k] = v[i_bh1][k]; + x_bh2[k] = x[i_bh2][k]; + v_bh2[k] = v[i_bh2][k]; + } - x_bh2[k] = x[i_bh2][k]; - v_bh2[k] = v[i_bh2][k]; - } + // calculate and "minus" the BH <-> BH _softened_ pot, acc & jerk - // calculate and "minus" the BH <-> BH _softened_ pot, acc & jerk + tmp_i = 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); - tmp_i = 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[i_bh1] -= pot_bh1; + pot[i_bh2] -= pot_bh2; - pot[i_bh1] -= pot_bh1; - pot[i_bh2] -= pot_bh2; + for (k=0;k<3;k++) { + a[i_bh1][k] -= a_bh1[k]; + a[i_bh2][k] -= a_bh2[k]; - for (k=0;k<3;k++) { - a[i_bh1][k] -= a_bh1[k]; - a[i_bh2][k] -= a_bh2[k]; + adot[i_bh1][k] -= adot_bh1[k]; + adot[i_bh2][k] -= adot_bh2[k]; + } - adot[i_bh1][k] -= adot_bh1[k]; - adot[i_bh2][k] -= adot_bh2[k]; - } + // calculate and "plus" the new BH <-> BH _unsoftened_ pot, acc, jerk - // calculate and "plus" the new BH <-> BH _unsoftened_ pot, acc, jerk + tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1, + m_bh2, x_bh2, v_bh2, + parameters.live_smbh_custom_eps, + &pot_bh1, a_bh1, adot_bh1, + &pot_bh2, a_bh2, adot_bh2); - tmp_i = 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[i_bh1] += pot_bh1; + pot[i_bh2] += pot_bh2; - pot[i_bh1] += pot_bh1; - pot[i_bh2] += pot_bh2; + for (k=0;k<3;k++) { + a[i_bh1][k] += a_bh1[k]; + a[i_bh2][k] += a_bh2[k]; - for (k=0;k<3;k++) { - a[i_bh1][k] += a_bh1[k]; - a[i_bh2][k] += a_bh2[k]; + adot[i_bh1][k] += adot_bh1[k]; + adot[i_bh2][k] += adot_bh2[k]; + } + } - adot[i_bh1][k] += adot_bh1[k]; - adot[i_bh2][k] += adot_bh2[k]; - } + if (parameters.binary_smbh_pn) { -#endif // ADD_N_BH + // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk -#ifdef ADD_PN_BH + dt_bh_tmp = dt[0]; - // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk + tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1, + m_bh2, x_bh2, v_bh2, s_bh2, + C_NB, dt_bh_tmp, usedOrNot, + (double(*)[3])a_pn1, (double(*)[3])adot_pn1, + (double(*)[3])a_pn2, (double(*)[3])adot_pn2); - dt_bh_tmp = dt[0]; + for (k=0;k<3;k++) { + a[i_bh1][k] += a_pn1[1][k] + a_pn1[2][k] + a_pn1[3][k] + a_pn1[4][k] + a_pn1[5][k] + a_pn1[6][k]; + a[i_bh2][k] += a_pn2[1][k] + a_pn2[2][k] + a_pn2[3][k] + a_pn2[4][k] + a_pn2[5][k] + a_pn2[6][k]; - tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1, - m_bh2, x_bh2, v_bh2, s_bh2, - C_NB, dt_bh_tmp, usedOrNot, - a_pn1, adot_pn1, - a_pn2, adot_pn2); + adot[i_bh1][k] += adot_pn1[1][k] + adot_pn1[2][k] + adot_pn1[3][k] + adot_pn1[4][k] + adot_pn1[5][k] + adot_pn1[6][k]; + adot[i_bh2][k] += adot_pn2[1][k] + adot_pn2[2][k] + adot_pn2[3][k] + adot_pn2[4][k] + adot_pn2[5][k] + adot_pn2[6][k]; + } - for (k=0;k<3;k++) { - a[i_bh1][k] += a_pn1[1][k] + a_pn1[2][k] + a_pn1[3][k] + a_pn1[4][k] + a_pn1[5][k] + a_pn1[6][k]; - a[i_bh2][k] += a_pn2[1][k] + a_pn2[2][k] + a_pn2[3][k] + a_pn2[4][k] + a_pn2[5][k] + a_pn2[6][k]; - - adot[i_bh1][k] += adot_pn1[1][k] + adot_pn1[2][k] + adot_pn1[3][k] + adot_pn1[4][k] + adot_pn1[5][k] + adot_pn1[6][k]; - adot[i_bh2][k] += adot_pn2[1][k] + adot_pn2[2][k] + adot_pn2[3][k] + adot_pn2[4][k] + adot_pn2[5][k] + adot_pn2[6][k]; - } - - if (myRank == rootRank) { - if (tmp_i == 505) { - printf("PN RSDIST: %.8E \t %.8E \n", timesteps, time_cur); - fflush(stdout); - exit(-1); + if (myRank == rootRank) { + if (tmp_i == 505) { + printf("PN RSDIST: %.8E \t %.8E \n", timesteps, time_cur); + fflush(stdout); + exit(-1); + } + } } } -#endif // ADD_PN_BH - -#endif // ADD_BH2 - #ifdef EXTPOT calc_ext_grav_zero(); #endif @@ -2236,47 +2020,19 @@ int main(int argc, char *argv[]) dt[i] = dt_tmp; -#ifdef DT_MIN_WARNING - if (myRank == 0) { + if (parameters.dt_min_warning && (myRank == 0)) { if (dt[i] == dt_min) { printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[i], ind[i]); fflush(stdout); } } -#endif } /* i */ -#ifdef ADD_BH2 - - /* define the min. dt over all the part. and set it also for the BH... */ - - min_dt = dt[0]; - - for (i=1; i 0) { + double min_dt = *std::min_element(dt, dt+N); + for (int i=0; i 0) { + i_bh1 = 0; + i_bh2 = 1; + } #endif #ifdef TIMING @@ -2522,91 +2280,78 @@ int main(int argc, char *argv[]) DT_ACT_REDUCE += (CPU_tmp_user - CPU_tmp_user0); #endif -#ifdef ADD_BH2 + if (parameters.live_smbh_count == 2) { + if (parameters.live_smbh_custom_eps >= 0) { + m_bh1 = m_act[i_bh1]; + m_bh2 = m_act[i_bh2]; -#ifdef ADD_N_BH + x_bh1 = x_act_new[i_bh1]; + v_bh1 = v_act_new[i_bh1]; - m_bh1 = m_act[i_bh1]; - m_bh2 = m_act[i_bh2]; + x_bh2 = x_act_new[i_bh2]; + v_bh2 = v_act_new[i_bh2]; - for (k=0;k<3;k++) { - x_bh1[k] = x_act_new[i_bh1][k]; - v_bh1[k] = v_act_new[i_bh1][k]; + // calculate and "minus" the BH <-> BH softened pot, acc & jerk - x_bh2[k] = x_act_new[i_bh2][k]; - v_bh2[k] = v_act_new[i_bh2][k]; - } + tmp_i = 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); - // calculate and "minus" the BH <-> BH softened pot, acc & jerk + pot_act_new[i_bh1] -= pot_bh1; + pot_act_new[i_bh2] -= pot_bh2; - tmp_i = 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); + a_act_new[i_bh1] -= a_bh1; + a_act_new[i_bh2] -= a_bh2; - pot_act_new[i_bh1] -= pot_bh1; - pot_act_new[i_bh2] -= pot_bh2; + adot_act_new[i_bh1] -= adot_bh1; + adot_act_new[i_bh2] -= adot_bh2; - for (k=0;k<3;k++) { - a_act_new[i_bh1][k] -= a_bh1[k]; - a_act_new[i_bh2][k] -= a_bh2[k]; + // calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk - adot_act_new[i_bh1][k] -= adot_bh1[k]; - adot_act_new[i_bh2][k] -= adot_bh2[k]; - } + tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1, + m_bh2, x_bh2, v_bh2, + parameters.live_smbh_custom_eps, + &pot_bh1, a_bh1, adot_bh1, + &pot_bh2, a_bh2, adot_bh2); - // calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk + pot_act_new[i_bh1] += pot_bh1; + pot_act_new[i_bh2] += pot_bh2; - tmp_i = 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); + a_act_new[i_bh1] += a_bh1; + a_act_new[i_bh2] += a_bh2; - pot_act_new[i_bh1] += pot_bh1; - pot_act_new[i_bh2] += pot_bh2; + adot_act_new[i_bh1] += adot_bh1; + adot_act_new[i_bh2] += adot_bh2; + } - for (k=0;k<3;k++) { - a_act_new[i_bh1][k] += a_bh1[k]; - a_act_new[i_bh2][k] += a_bh2[k]; + if (parameters.binary_smbh_pn) { + // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk - adot_act_new[i_bh1][k] += adot_bh1[k]; - adot_act_new[i_bh2][k] += adot_bh2[k]; - } + dt_bh_tmp = dt[0]; -#endif // ADD_N_BH + tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1, + m_bh2, x_bh2, v_bh2, s_bh2, + C_NB, dt_bh_tmp, usedOrNot, + (double(*)[3])a_pn1, (double(*)[3])adot_pn1, + (double(*)[3])a_pn2, (double(*)[3])adot_pn2); -#ifdef ADD_PN_BH + 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]; - // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk + 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]; - dt_bh_tmp = dt[0]; - - tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1, - m_bh2, x_bh2, v_bh2, s_bh2, - C_NB, dt_bh_tmp, usedOrNot, - a_pn1, adot_pn1, - a_pn2, adot_pn2); - - for (k=0;k<3;k++) { - a_act_new[i_bh1][k] += a_pn1[1][k] + a_pn1[2][k] + a_pn1[3][k] + a_pn1[4][k] + a_pn1[5][k] + a_pn1[6][k]; - a_act_new[i_bh2][k] += a_pn2[1][k] + a_pn2[2][k] + a_pn2[3][k] + a_pn2[4][k] + a_pn2[5][k] + a_pn2[6][k]; - - adot_act_new[i_bh1][k] += adot_pn1[1][k] + adot_pn1[2][k] + adot_pn1[3][k] + adot_pn1[4][k] + adot_pn1[5][k] + adot_pn1[6][k]; - adot_act_new[i_bh2][k] += adot_pn2[1][k] + adot_pn2[2][k] + adot_pn2[3][k] + adot_pn2[4][k] + adot_pn2[5][k] + adot_pn2[6][k]; - } - - if (myRank == rootRank) { - if (tmp_i == 505) { - printf("PN RSDIST: TS = %.8E \t t = %.8E \n", timesteps, time_cur); - fflush(stdout); - exit(-1); + if (myRank == rootRank) { + if (tmp_i == 505) { + printf("PN RSDIST: TS = %.8E \t t = %.8E \n", timesteps, time_cur); + fflush(stdout); + exit(-1); + } + } } } -#endif // ADD_PN_BH - -#endif // ADD_BH2 #ifdef EXTPOT calc_ext_grav(); @@ -2647,25 +2392,10 @@ int main(int argc, char *argv[]) a2dot1abs = a2dot1.norm(); a3dot1abs = a3.norm(); -#ifdef ADD_BH2 - if ((ind_act[i] == 0) || (ind_act[i] == 1)) { + if ((parameters.live_smbh_count > 0) && (ind_act[i] < parameters.live_smbh_count)) dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); - } else { + else dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); - } -#else - -#ifdef ADD_BH1 - if (ind_act[i] == 0) { - dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); - } else { - dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); - } -#endif // ADD_BH1 - - dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); - -#endif // ADD_BH2 if (dt_new < dt_min) dt_tmp = dt_min; @@ -2690,114 +2420,92 @@ int main(int argc, char *argv[]) a_act[i] = a_act_new[i]; adot_act[i] = adot_act_new[i]; -#ifdef DT_MIN_WARNING - if (myRank == 0) { + if (parameters.dt_min_warning && (myRank == 0)) { if (dt_act[i] == dt_min) { printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt[i], ind_act[i]); fflush(stdout); } } -#endif + } /* i */ /* define the min. dt over all the act. part. and set it also for the BH... */ -#ifdef ADD_BH2 - min_dt = dt_act[0]; + if (parameters.live_smbh_count > 0) { + double min_dt = *std::min_element(dt_act, dt_act+n_act); + if (parameters.live_smbh_count>=1) dt_act[i_bh1] = min_dt; + if (parameters.live_smbh_count==2) dt_act[i_bh2] = min_dt; + } - for (i=1; i= t_bh) { if (myRank == rootRank) { -#ifdef BH_OUT /* Write BH data... */ - write_bh_data(); -#endif + if (parameters.live_smbh_output) write_bh_data(); -#ifdef BH_OUT_NB /* Write BH NB data... */ - write_bh_nb_data(); -#endif + if (parameters.live_smbh_neighbor_output) write_bh_nb_data(); } /* if (myRank == rootRank) */ @@ -2921,8 +2625,6 @@ int main(int argc, char *argv[]) } /* if (myRank == rootRank) */ - /* possible coordinate & velocity limits for ALL particles !!! */ - #ifdef ETICS_CEP // We are /inside/ a control step, so all particles must be synchronized; we can safely calculate their density centre. The acceleration and jerk currently in the memory are for the predicted position of the CEP, by calling grapite_calc_center we "correct" the position and velocity, but not the gravity at that point. @@ -2960,7 +2662,7 @@ int main(int argc, char *argv[]) } /* while (time_cur < t_end) */ - /* close the local GRAPE's */ + /* close the local GRAPEs */ g6_close(clusterid);