diff --git a/phigrape.cpp b/phigrape.cpp index 1b89885..5f52247 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -56,49 +56,8 @@ Last redaction : 2019.04.16 12:55 #include "config.h" Config *config; -// struct Parameters { -// double eps = 1E-4; -// double t_end = 1; -// double dt_disk = 0.125; -// double dt_contr = 0.125; -// double dt_bh = 0.125; -// double eta = 0.01; -// std::string input_file_name = "data.con"; -// bool dt_min_warning = true; // DT_MIN_WARNING -// 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` -// 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_influence_sphere_output = true; // BBH_INF // NOTE needs `binary_smbh_influence_radius_factor` -// double binary_smbh_influence_radius_factor = 3.162277660168379497918067e+03; // R_INF -// } parameters; - //#define NORM // Physical normalization -// // // //#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 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 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 BH_OUT_NB_EXT // extra output for the BH's neighbours (external) - -//#define STARDESTR // disruption of stars by BH tidal forces - -//#define STARDESTR_EXT // disruption of stars by external BH tidal forces - //#define EXTPOT // external potential (BH? or galactic?) //#define EXTPOT_BH // BH - usually NB units @@ -107,21 +66,6 @@ Config *config; //#define EXTPOT_GAL_DEH // Dehnen Galactic - usually physical units //#define EXTPOT_GAL_LOG // Log Galactic - usually physical units -//#define EXTPOT_SC // SC extra POT for Bek test runs... -//#define DATAFILE eff0.05.tsf10.00.tab -//#define M_SC_DIM 100001 - -//#define CMCORR // CM correction in the zero step and in every dt_contr - -//#define LIMITS // for "0" mass particles !!! -//#define R_LIMITS 1.0E+03 // for "0" mass particles !!! - -//#define LIMITS_NEW // for "0" mass particles !!! - -//#define LIMITS_ALL_BEG // for ALL particles at the beginning... -//#define LIMITS_ALL // for ALL particles -//#define R_LIMITS_ALL 1.0E+03 // for ALL particles - #ifdef ETICS #include "grapite.h" // why do we need CEP as a compilaion flag... just have it always on when ETICS is on. IF there is no CEP, there should be a graceful skipping of those operations. @@ -592,85 +536,37 @@ void write_bh_data(double time_cur, double m[], double3 x[], double3 v[], double a_bh = a_bh2; adot_bh = adot_bh2; } - - double tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) ); - double tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) ); - - double tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) ); - double tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) ); - if (config->live_smbh_custom_eps >= 0) { - - double tmp_a_bh = sqrt( SQR(a_bh[0]) + SQR(a_bh[1]) + SQR(a_bh[2]) ); - double tmp_adot_bh = sqrt( SQR(adot_bh[0]) + SQR(adot_bh[1]) + SQR(adot_bh[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 ", + 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) { - double tmp_a_bh_pn0 = sqrt( SQR(a_pn[0][0]) + SQR(a_pn[0][1]) + SQR(a_pn[0][2]) ); - double tmp_a_bh_pn1 = sqrt( SQR(a_pn[1][0]) + SQR(a_pn[1][1]) + SQR(a_pn[1][2]) ); - double tmp_a_bh_pn2 = sqrt( SQR(a_pn[2][0]) + SQR(a_pn[2][1]) + SQR(a_pn[2][2]) ); - double tmp_a_bh_pn2_5 = sqrt( SQR(a_pn[3][0]) + SQR(a_pn[3][1]) + SQR(a_pn[3][2]) ); - double tmp_a_bh_pn3 = sqrt( SQR(a_pn[4][0]) + SQR(a_pn[4][1]) + SQR(a_pn[4][2]) ); - double tmp_a_bh_pn3_5 = sqrt( SQR(a_pn[5][0]) + SQR(a_pn[5][1]) + SQR(a_pn[5][2]) ); - double tmp_a_bh_spin = sqrt( SQR(a_pn[6][0]) + SQR(a_pn[6][1]) + SQR(a_pn[6][2]) ); - - double tmp_adot_bh_pn0 = sqrt( SQR(adot_pn[0][0]) + SQR(adot_pn[0][1]) + SQR(adot_pn[0][2]) ); - double tmp_adot_bh_pn1 = sqrt( SQR(adot_pn[1][0]) + SQR(adot_pn[1][1]) + SQR(adot_pn[1][2]) ); - double tmp_adot_bh_pn2 = sqrt( SQR(adot_pn[2][0]) + SQR(adot_pn[2][1]) + SQR(adot_pn[2][2]) ); - double tmp_adot_bh_pn2_5 = sqrt( SQR(adot_pn[3][0]) + SQR(adot_pn[3][1]) + SQR(adot_pn[3][2]) ); - double tmp_adot_bh_pn3 = sqrt( SQR(adot_pn[4][0]) + SQR(adot_pn[4][1]) + SQR(adot_pn[4][2]) ); - double tmp_adot_bh_pn3_5 = sqrt( SQR(adot_pn[5][0]) + SQR(adot_pn[5][1]) + SQR(adot_pn[5][2]) ); - double 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); + 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], tmp_r, - v[i][0], v[i][1], v[i][2], tmp_v, + 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], tmp_a, - adot[i][0], adot[i][1], adot[i][2], tmp_adot, + 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"); + fprintf(out, "\n"); fclose(out); } else if (config->live_smbh_count == 1) { auto out = fopen("bh.dat", "a"); @@ -991,11 +887,6 @@ void calc_self_grav(double t, double eps2, double &g6_calls, int n_loc, get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); #endif -// // // // // // // /* bad variables... */ -// // // // // // // double pot_act_tmp[N_MAX]; -// // // // // // // double3 a_act_tmp[N_MAX]; -// // // // // // // double3 adot_act_tmp[N_MAX]; - g6_set_ti(clusterid, t); int ni = n_act; // TODO why is this needed? @@ -1014,15 +905,6 @@ void calc_self_grav(double t, double eps2, double &g6_calls, int n_loc, g6_calls++; } /* i */ - /* Store the new value of the local partial force etc... */ - -// // // // // // // // for (int i=0; ieta; strcpy(inp_fname, config->input_file_name.c_str()); - if (config->binary_smbh_influence_sphere_output) for (i=0; ibinary_smbh_influence_sphere_output) for (int i=0; i BH _softened_ pot, acc & jerk @@ -1887,13 +1759,11 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, 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]; + a[i_bh1] -= a_bh1; + a[i_bh2] -= a_bh2; - adot[i_bh1][k] -= adot_bh1[k]; - adot[i_bh2][k] -= adot_bh2[k]; - } + adot[i_bh1] -= adot_bh1; + adot[i_bh2] -= adot_bh2; // calculate and "plus" the new BH <-> BH _unsoftened_ pot, acc, jerk @@ -1906,13 +1776,11 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, 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]; + a[i_bh1] += a_bh1; + a[i_bh2] += a_bh2; - adot[i_bh1][k] += adot_bh1[k]; - adot[i_bh2][k] += adot_bh2[k]; - } + adot[i_bh1] += adot_bh1; + adot[i_bh2] += adot_bh2; } if (config->binary_smbh_pn) { @@ -1927,13 +1795,11 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, (double(*)[3])a_pn1, (double(*)[3])adot_pn1, (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank); - 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]; + 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][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]; - } + 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_i == 505) { @@ -1997,7 +1863,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, /* Define initial timestep for all particles on all nodes */ - for (i=0; ilive_smbh_count > 0) { @@ -2044,7 +1909,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, /* load the nj particles to the G6 */ - for (j=0; j