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

Cleaned unused variables and fixed a bug with PN (also tested PN against old version)

Browse source
This commit is contained in:
Yohai Meiron 2020-04-07 22:40:54 -04:00
parent d9e3aea243
commit ebec3280f8
4 changed files with 41 additions and 51 deletions
Download patch file Download diff file Expand all files Collapse all files

2
.gitignore vendored
View file

@ -9,3 +9,5 @@
grapite-dev-exec-threshold grapite-dev-exec-threshold
phigrape phigrape
*.h5 *.h5
.*
.gitignore

2
io.cpp
View file

@ -81,7 +81,7 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub
H5Gget_info(group_id, &object_info); H5Gget_info(group_id, &object_info);
// Iterate over objects and add the number of each "step" group into a vector // Iterate over objects and add the number of each "step" group into a vector
std::vector<int> step_num_arr; std::vector<int> step_num_arr;
for (int i=0; i<object_info.nlinks; i++) { for (unsigned int i=0; i<object_info.nlinks; i++) {
char name_cstr[256]; char name_cstr[256];
H5Gget_objname_by_idx(group_id, i, name_cstr, 256); H5Gget_objname_by_idx(group_id, i, name_cstr, 256);
std::string name(name_cstr); std::string name(name_cstr);

4
n_bh.c
View file

@ -5,7 +5,7 @@
Last redaction : 2011.II.20. 13:00 Last redaction : 2011.II.20. 13:00
*/ */
int calc_force_n_BH(double m1, double xx1[], double vv1[], void calc_force_n_BH(double m1, double xx1[], double vv1[],
double m2, double xx2[], double vv2[], double m2, double xx2[], double vv2[],
double eps_BH, double eps_BH,
double *pot_n1, double a_n1[], double adot_n1[], double *pot_n1, double a_n1[], double adot_n1[],
@ -70,7 +70,5 @@ for(k=0;k<3;k++)
adot_n2[k] = m1*(v[k]/r3 - RP*x[k]/r4); adot_n2[k] = m1*(v[k]/r3 - RP*x[k]/r4);
} }
return(0);
} }
/***************************************************************************/ /***************************************************************************/

40
phigrape.cpp
View file

@ -103,7 +103,7 @@ Config *config;
#define KB 1024 #define KB 1024
#define MB (KB*KB) #define MB (KB*KB)
#define N_MAX (6*MB) #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 // Needed for things related to BHs
@ -336,11 +336,12 @@ void calc_self_grav(double t, double eps2, double &g6_calls, int n_loc,
void calc_ext_grav(std::vector<External_gravity*> &external_gravity_components, int n_act, double3 *x_act_new, double3 *v_act_new, double *pot_act_ext, double3 *a_act_new, double3* adot_act_new) void calc_ext_grav(std::vector<External_gravity*> &external_gravity_components, int n_act, double3 *x_act_new, double3 *v_act_new, double *pot_act_ext, double3 *a_act_new, double3* adot_act_new)
{ {
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
/* Define the external potential for all active particles on all nodes */ /* Define the external potential for all active particles on all nodes */
int ni = n_act; // TODO redundant?
std::fill(pot_act_ext, pot_act_ext+n_act, 0.); std::fill(pot_act_ext, pot_act_ext+n_act, 0.);
for (auto component : external_gravity_components) { for (auto component : external_gravity_components) {
@ -348,10 +349,6 @@ int ni = n_act; // TODO redundant?
component->apply(n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new); component->apply(n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new);
} }
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
/* For simple Plummer potential... */ /* For simple Plummer potential... */
#ifdef TIMING #ifdef TIMING
@ -679,8 +676,7 @@ void adjust_bsmbh_softening(const double eps, const double eps_bh, const int i_b
// calculate and "minus" the BH <-> BH softened pot, acc & jerk // calculate and "minus" the BH <-> BH softened pot, acc & jerk
double tmp_i; calc_force_n_BH(m_bh1, x_bh1, v_bh1,
tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
m_bh2, x_bh2, v_bh2, m_bh2, x_bh2, v_bh2,
eps, eps,
&pot_bh1, a_bh1, adot_bh1, &pot_bh1, a_bh1, adot_bh1,
@ -697,7 +693,7 @@ void adjust_bsmbh_softening(const double eps, const double eps_bh, const int i_b
// 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, calc_force_n_BH(m_bh1, x_bh1, v_bh1,
m_bh2, x_bh2, v_bh2, m_bh2, x_bh2, v_bh2,
eps_bh, eps_bh,
&pot_bh1, a_bh1, adot_bh1, &pot_bh1, a_bh1, adot_bh1,
@ -711,28 +707,24 @@ void adjust_bsmbh_softening(const double eps, const double eps_bh, const int i_b
adot_act_new[i_bh1] += adot_bh1; adot_act_new[i_bh1] += adot_bh1;
adot_act_new[i_bh2] += adot_bh2; adot_act_new[i_bh2] += adot_bh2;
} }
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank, int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank,
nj, diskstep=0, power, jjj, iii, diskstep=0, power, jjj, iii,
skip_con=0, tmp_i; skip_con=0, tmp_i;
double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
dt_bh, t_bh=0.0, dt_bh_tmp, dt_bh, t_bh=0.0, dt_bh_tmp,
t_end, time_cur, dt_min, dt_max, min_t, min_t_loc, dt_new, t_end, time_cur, dt_min, dt_max, min_t, min_t_loc,
eta_s, eta, eta_bh, eta_s, eta, eta_bh,
E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0,
rcm_sum=0.0, vcm_sum=0.0, rcm_sum=0.0, vcm_sum=0.0,
eps=0.0, eps2, eps=0.0, eps2,
a2_mod, adot2_mod, a2_mod, adot2_mod,
dt_tmp, dt2half, dt3over6, dt4over24, dt5over120, dt_tmp, dt2half, dt3over6,
dtinv, dt2inv, dt3inv, timesteps=0.0, n_act_sum=0.0, n_act_distr[N_MAX], g6_calls=0.0, g6_calls_sum=0.0;
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,
@ -785,9 +777,7 @@ int main(int argc, char *argv[])
double m_bh1, m_bh2; double m_bh1, m_bh2;
int inf_event[N_MAX]; int inf_event[N_MAX];
double DR2, tmp_r2;
double3 x_bbhc, v_bbhc; double3 x_bbhc, v_bbhc;
double DV2, EB, SEMI_a, SEMI_a2;
double s_bh1[3] = {0.0, 0.0, 1.0}; double s_bh1[3] = {0.0, 0.0, 1.0};
double s_bh2[3] = {0.0, 0.0, 1.0}; double s_bh2[3] = {0.0, 0.0, 1.0};
@ -823,7 +813,7 @@ int main(int argc, char *argv[])
if (is_hdf5(config->input_file_name)) { if (is_hdf5(config->input_file_name)) {
#ifndef HAS_HDF5 #ifndef HAS_HDF5
fprintf(stderr, "ERROR: input file is in HDF5 format, but the code was compiled without HDF5 support\n") fprintf(stderr, "ERROR: input file is in HDF5 format, but the code was compiled without HDF5 support\n");
return -1; return -1;
#endif #endif
h5_read(config->input_file_name, &diskstep, &N, &time_cur, m, x, v); h5_read(config->input_file_name, &diskstep, &N, &time_cur, m, x, v);
@ -1097,7 +1087,7 @@ int main(int argc, char *argv[])
// 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, calc_force_n_BH(m_bh1, x_bh1, v_bh1,
m_bh2, x_bh2, v_bh2, m_bh2, x_bh2, v_bh2,
eps, eps,
&pot_bh1, a_bh1, adot_bh1, &pot_bh1, a_bh1, adot_bh1,
@ -1114,7 +1104,7 @@ int main(int argc, char *argv[])
// 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, calc_force_n_BH(m_bh1, x_bh1, v_bh1,
m_bh2, x_bh2, v_bh2, m_bh2, x_bh2, v_bh2,
config->live_smbh_custom_eps, config->live_smbh_custom_eps,
&pot_bh1, a_bh1, adot_bh1, &pot_bh1, a_bh1, adot_bh1,
@ -1427,8 +1417,8 @@ int main(int argc, char *argv[])
dt_bh_tmp = dt[0]; dt_bh_tmp = dt[0];
tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1, tmp_i = calc_force_pn_BH(m_bh1, x_act_new[i_bh1], v_act_new[i_bh1], s_bh1,
m_bh2, x_bh2, v_bh2, s_bh2, m_bh2, x_act_new[i_bh2], v_act_new[i_bh2], s_bh2,
C_NB, dt_bh_tmp, usedOrNot, C_NB, dt_bh_tmp, usedOrNot,
(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);