yohai/phigrape
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Prettied up the black hole output and removed some redundant variables

Browse source
This commit is contained in:
Yohai Meiron 2020-03-21 12:46:26 -04:00
parent 8f984b5158
commit 9b900eff6e
1 changed files with 70 additions and 206 deletions
Download patch file Download diff file Expand all files Collapse all files

274
phigrape.cpp
View file

@ -56,49 +56,8 @@ Last redaction : 2019.04.16 12:55
#include "config.h" #include "config.h"
Config *config; 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 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 // external potential (BH? or galactic?)
//#define EXTPOT_BH // BH - usually NB units //#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_DEH // Dehnen Galactic - usually physical units
//#define EXTPOT_GAL_LOG // Log 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 #ifdef ETICS
#include "grapite.h" #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. // 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; a_bh = a_bh2;
adot_bh = adot_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) { 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 ",
double tmp_a_bh = sqrt( SQR(a_bh[0]) + SQR(a_bh[1]) + SQR(a_bh[2]) ); time_cur, m[i],
double tmp_adot_bh = sqrt( SQR(adot_bh[0]) + SQR(adot_bh[1]) + SQR(adot_bh[2]) ); 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) { 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]) ); fprintf(out, "\t");
double tmp_a_bh_pn1 = sqrt( SQR(a_pn[1][0]) + SQR(a_pn[1][1]) + SQR(a_pn[1][2]) ); for (int pn_idx=0; pn_idx < 7; pn_idx++) {
double tmp_a_bh_pn2 = sqrt( SQR(a_pn[2][0]) + SQR(a_pn[2][1]) + SQR(a_pn[2][2]) ); 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());
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, "\n");
} else { } 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", 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], time_cur, m[i],
x[i][0], x[i][1], x[i][2], tmp_r, x[i][0], x[i][1], x[i][2], x[i].norm(),
v[i][0], v[i][1], v[i][2], tmp_v, v[i][0], v[i][1], v[i][2], v[i].norm(),
pot[i], pot[i],
a[i][0], a[i][1], a[i][2], tmp_a, a[i][0], a[i][1], a[i][2], a[i].norm(),
adot[i][0], adot[i][1], adot[i][2], tmp_adot, adot[i][0], adot[i][1], adot[i][2], adot[i].norm(),
dt[i]); dt[i]);
} }
} }
fprintf(out,"\n"); fprintf(out, "\n");
fclose(out); fclose(out);
} else if (config->live_smbh_count == 1) { } else if (config->live_smbh_count == 1) {
auto out = fopen("bh.dat", "a"); 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); get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif #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); g6_set_ti(clusterid, t);
int ni = n_act; // TODO why is this needed? 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++; g6_calls++;
} /* i */ } /* i */
/* Store the new value of the local partial force etc... */
// // // // // // // // for (int i=0; i<n_act; i++) {
// // // // // // // // int iii = ind_act[i];
// // // // // // // // pot_act_tmp_loc[iii] = pot_act_tmp[i];
// // // // // // // // a_act_tmp_loc[iii] = a_act_tmp[i];
// // // // // // // // adot_act_tmp_loc[iii] = adot_act_tmp[i];
// // // // // // // // } /* i */
#ifdef TIMING #ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst); get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_ACT_GRAV += (CPU_tmp_user - CPU_tmp_user0); DT_ACT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
@ -1342,7 +1224,6 @@ void energy_contr(const double time_cur, const double timesteps, const double n_
double E_tot = E_pot + E_kin + E_pot_ext; double E_tot = E_pot + E_kin + E_pot_ext;
double E_tot_corr = E_tot /* + E_corr */; // TODO what is this? double E_tot_corr = E_tot /* + E_corr */; // TODO what is this?
if (timesteps == 0.0) { if (timesteps == 0.0) {
/* initialize the E_tot_0 + etc... */ /* initialize the E_tot_0 + etc... */
@ -1404,27 +1285,23 @@ void energy_contr(const double time_cur, const double timesteps, const double n_
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank,
nj, diskstep=0, power, jjj, iii,
skip_con=0, tmp_i;
// These variables are moved from global double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank, dt_bh, t_bh=0.0, dt_bh_tmp,
i, j, k, nj, diskstep=0, power, jjj, iii, t_end, time_cur, dt_min, dt_max, min_t, min_t_loc, dt_new,
skip_con=0, tmp_i; eta_s, eta, eta_bh,
E_tot_0, E_tot_corr_0, E_tot_corr_sd_0,
double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, rcm_sum=0.0, vcm_sum=0.0,
dt_bh, t_bh=0.0, dt_bh_tmp, eps=0.0, eps2,
t_end, time_cur, dt_min, dt_max, min_t, min_t_loc, dt_new, a2_mod, adot2_mod,
eta_s, eta, eta_bh, dt_tmp, dt2half, dt3over6, dt4over24, dt5over120,
E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, dtinv, dt2inv, dt3inv,
rcm_sum=0.0, vcm_sum=0.0, a1abs, adot1abs, a2dot1abs, a3dot1abs,
eps=0.0, eps2, timesteps=0.0, n_act_sum=0.0, n_act_distr[N_MAX], g6_calls=0.0, g6_calls_sum=0.0,
a2_mod, adot2_mod, tmp;
dt_tmp, dt2half, dt3over6, dt4over24, dt5over120,
dtinv, dt2inv, dt3inv,
a1abs, adot1abs, a2dot1abs, a3dot1abs,
timesteps=0.0, n_act_sum=0.0, n_act_distr[N_MAX], g6_calls=0.0, g6_calls_sum=0.0,
tmp;
double3 xcm, vcm, mom, double3 xcm, vcm, mom,
xdc, vdc, xdc, vdc,
@ -1499,7 +1376,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
eta = config->eta; eta = config->eta;
strcpy(inp_fname, config->input_file_name.c_str()); strcpy(inp_fname, config->input_file_name.c_str());
if (config->binary_smbh_influence_sphere_output) for (i=0; i<N; i++) inf_event[i] = 0; // WARNING N wasn't set yet! if (config->binary_smbh_influence_sphere_output) for (int i=0; i<N; i++) inf_event[i] = 0; // WARNING N wasn't set yet!
/* /*
eps : Plummer softening parameter (can be even 0) eps : Plummer softening parameter (can be even 0)
@ -1724,10 +1601,8 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
fflush(stdout); fflush(stdout);
} /* if (myRank == rootRank) */ } /* if (myRank == rootRank) */
for (i=0; i<N; i++) { std::fill(t, t+N, time_cur);
t[i] = time_cur; std::fill(dt, dt+N, dt_min);
dt[i] = dt_min;
}
n_loc = N/n_proc; n_loc = N/n_proc;
@ -1786,12 +1661,11 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
/* load the nj particles to the G6 */ /* load the nj particles to the G6 */
for (j=0; j<nj; j++) { for (int j=0; j<nj; j++) {
a2by18 = {0, 0, 0}; a2by18 = {0, 0, 0};
a1by6 = {0, 0, 0}; a1by6 = {0, 0, 0};
aby2 = {0, 0, 0}; aby2 = {0, 0, 0};
g6_set_j_particle(clusterid, j, ind_loc[j], t_loc[j], dt_loc[j], m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]); g6_set_j_particle(clusterid, j, ind_loc[j], t_loc[j], dt_loc[j], m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]);
} /* j */ } /* j */
#ifdef ETICS #ifdef ETICS
@ -1826,7 +1700,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
n_act = N; n_act = N;
for (i=0; i<n_act; i++) { for (int i=0; i<n_act; i++) {
ind_act[i] = ind[i]; ind_act[i] = ind[i];
iii = ind_act[i]; iii = ind_act[i];
m_act[i] = m[iii]; m_act[i] = m[iii];
@ -1868,13 +1742,11 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
m_bh1 = m[i_bh1]; m_bh1 = m[i_bh1];
m_bh2 = m[i_bh2]; m_bh2 = m[i_bh2];
for (k=0;k<3;k++) { x_bh1 = x[i_bh1];
x_bh1[k] = x[i_bh1][k]; v_bh1 = v[i_bh1];
v_bh1[k] = v[i_bh1][k];
x_bh2[k] = x[i_bh2][k]; x_bh2 = x[i_bh2];
v_bh2[k] = v[i_bh2][k]; v_bh2 = v[i_bh2];
}
// calculate and "minus" the BH <-> BH _softened_ pot, acc & jerk // calculate and "minus" the BH <-> 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_bh1] -= pot_bh1;
pot[i_bh2] -= pot_bh2; pot[i_bh2] -= pot_bh2;
for (k=0;k<3;k++) { a[i_bh1] -= a_bh1;
a[i_bh1][k] -= a_bh1[k]; a[i_bh2] -= a_bh2;
a[i_bh2][k] -= a_bh2[k];
adot[i_bh1][k] -= adot_bh1[k]; adot[i_bh1] -= adot_bh1;
adot[i_bh2][k] -= adot_bh2[k]; adot[i_bh2] -= adot_bh2;
}
// calculate and "plus" the new BH <-> BH _unsoftened_ pot, acc, jerk // 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_bh1] += pot_bh1;
pot[i_bh2] += pot_bh2; pot[i_bh2] += pot_bh2;
for (k=0;k<3;k++) { a[i_bh1] += a_bh1;
a[i_bh1][k] += a_bh1[k]; a[i_bh2] += a_bh2;
a[i_bh2][k] += a_bh2[k];
adot[i_bh1][k] += adot_bh1[k]; adot[i_bh1] += adot_bh1;
adot[i_bh2][k] += adot_bh2[k]; adot[i_bh2] += adot_bh2;
}
} }
if (config->binary_smbh_pn) { 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_pn1, (double(*)[3])adot_pn1,
(double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank); (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank);
for (k=0;k<3;k++) { a[i_bh1] += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6];
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] += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6];
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_bh1] += adot_pn1[1] + adot_pn1[2] + adot_pn1[3] + adot_pn1[4] + adot_pn1[5] + adot_pn1[6];
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_bh2] += adot_pn2[1] + adot_pn2[2] + adot_pn2[3] + adot_pn2[4] + adot_pn2[5] + adot_pn2[6];
}
if (myRank == rootRank) { if (myRank == rootRank) {
if (tmp_i == 505) { 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 */ /* Define initial timestep for all particles on all nodes */
for (i=0; i<N; i++) { for (int i=0; i<N; i++) {
a2_mod = a[i][0]*a[i][0]+a[i][1]*a[i][1]+a[i][2]*a[i][2]; a2_mod = a[i][0]*a[i][0]+a[i][1]*a[i][1]+a[i][2]*a[i][2];
adot2_mod = adot[i][0]*adot[i][0]+adot[i][1]*adot[i][1]+adot[i][2]*adot[i][2]; adot2_mod = adot[i][0]*adot[i][0]+adot[i][1]*adot[i][1]+adot[i][2]*adot[i][2];
@ -2021,7 +1887,6 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
fflush(stdout); fflush(stdout);
} }
} }
} /* i */ } /* i */
if (config->live_smbh_count > 0) { if (config->live_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 */ /* load the nj particles to the G6 */
for (j=0; j<nj; j++) { for (int j=0; j<nj; j++) {
a2by18 = {0, 0, 0}; a2by18 = {0, 0, 0};
a1by6 = adot_loc[j]*(1./6.); a1by6 = adot_loc[j]*(1./6.);
aby2 = a_loc[j]*0.5; aby2 = a_loc[j]*0.5;
@ -2071,7 +1936,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
timesteps = 0.0; timesteps = 0.0;
n_act_sum = 0.0; n_act_sum = 0.0;
for (i=1; i<N+1; i++) { for (int i=1; i<N+1; i++) { // TODO why like this?
n_act_distr[i-1] = 0.0; n_act_distr[i-1] = 0.0;
} }
@ -2113,7 +1978,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
#else #else
min_t_loc = t[0]+dt[0]; min_t_loc = t[0]+dt[0];
for (j=0; j<n_loc; j++) { for (int j=0; j<n_loc; j++) {
jjj = j + myRank*n_loc; jjj = j + myRank*n_loc;
tmp = t[jjj] + dt[jjj]; tmp = t[jjj] + dt[jjj];
if (tmp < min_t_loc) min_t_loc = tmp; if (tmp < min_t_loc) min_t_loc = tmp;
@ -2156,7 +2021,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
#else #else
n_act = 0; n_act = 0;
for (i=0; i<N; i++) { for (int i=0; i<N; i++) {
if (t[i]+dt[i] == min_t ) { if (t[i]+dt[i] == min_t ) {
ind_act[n_act] = i; ind_act[n_act] = i;
n_act++; n_act++;
@ -2207,7 +2072,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
#ifdef TIMING #ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif #endif
for (i=0; i<n_act; i++) { for (int i=0; i<n_act; i++) {
iii = ind_act[i]; iii = ind_act[i];
m_act[i] = m[iii]; m_act[i] = m[iii];
@ -2240,7 +2105,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif #endif
for (i=0; i<n_act; i++) { for (int i=0; i<n_act; i++) {
dt_tmp = min_t - t_act[i]; dt_tmp = min_t - t_act[i];
dt2half = 0.5*dt_tmp*dt_tmp; dt2half = 0.5*dt_tmp*dt_tmp;
dt3over6 = (1./3.)*dt_tmp*dt2half; dt3over6 = (1./3.)*dt_tmp*dt2half;
@ -2363,7 +2228,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif #endif
for (i=0; i<n_act; i++) { for (int i=0; i<n_act; i++) {
dt_tmp = min_t - t_act[i]; dt_tmp = min_t - t_act[i];
dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0; dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0;
@ -2462,7 +2327,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB; SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB;
SEMI_a2 = SQR(SEMI_a); SEMI_a2 = SQR(SEMI_a);
for (i=2; i<n_act; i++) { for (int i=2; i<n_act; i++) {
tmp_r2 = SQR(x_act[i][0] - x_bbhc[0]) + SQR(x_act[i][1] - x_bbhc[1]) + SQR(x_act[i][2] - x_bbhc[2]); tmp_r2 = SQR(x_act[i][0] - x_bbhc[0]) + SQR(x_act[i][1] - x_bbhc[1]) + SQR(x_act[i][2] - x_bbhc[2]);
@ -2509,7 +2374,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
/* Return back the new coordinates + etc... of active part. to the global data... */ /* Return back the new coordinates + etc... of active part. to the global data... */
for (i=0; i<n_act; i++) { for (int i=0; i<n_act; i++) {
iii = ind_act[i]; iii = ind_act[i];
@ -2542,7 +2407,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif #endif
for (j=0; j<n_act; j++) { for (int j=0; j<n_act; j++) { // TODO would be nicer to use i instead of j here
#ifdef ETICS_CEP #ifdef ETICS_CEP
if (ind_act[j] == grapite_cep_index) grapite_update_cep(t_act[j], x_act[j], v_act[j], a_act[j], adot_act[j]); // All ranks should do it. if (ind_act[j] == grapite_cep_index) grapite_update_cep(t_act[j], x_act[j], v_act[j], a_act[j], adot_act[j]); // All ranks should do it.
#endif #endif
@ -2695,5 +2560,4 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
MPI_Finalize(); MPI_Finalize();
return 0; return 0;
} }