From 8cba5b5e1d24f8c2eac7589630fb6919d145400d Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Wed, 25 Mar 2020 16:58:46 -0400 Subject: [PATCH] Got rid of more global variables --- config.cpp | 2 +- phigrape.cpp | 158 ++++++++++----------------------------------------- 2 files changed, 32 insertions(+), 128 deletions(-) diff --git a/config.cpp b/config.cpp index aac05c2..301e259 100644 --- a/config.cpp +++ b/config.cpp @@ -121,7 +121,7 @@ void Config::error_checking() throw std::runtime_error("Binary black hole influence sphere output (binary_smbh_influence_sphere_output=true) requires exactly two live black holes (live_smbh_count=2)"); if (pn_usage.size() != 7) throw std::runtime_error("PN usage array (pn_usage) must have exactly seven components"); - for (int i; i<7; i++) + for (int i=0; i<7; i++) if (!((pn_usage[i] == 0) || (pn_usage[i] == 1))) throw std::runtime_error("PN usage array (pn_usage) must have ones and zeros only"); if (ext_units_physical && ((norm_mass == 0) || (norm_length == 0))) diff --git a/phigrape.cpp b/phigrape.cpp index e4153a1..99d88e7 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -325,26 +325,6 @@ double r2, rv_ij, // 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 - - -int i_bh, i_bh1, i_bh2, - num_bh = 0, num_bh1 = 0, num_bh2 = 0; - -double m_bh, m_bh1, m_bh2, b_bh, - r, r2, - x_ij, y_ij, z_ij, - vx_ij, vy_ij, vz_ij, rv_ij; - -// #ifdef BBH_INF -int inf_event[N_MAX]; -double DR2, tmp_r2; -double3 x_bbhc, v_bbhc; -double DV2, EB, SEMI_a, SEMI_a2; -// #endif // BBH_INF - -// #ifdef ADD_N_BH - // double x_bh1[3], x_bh2[3], v_bh1[3], v_bh2[3]; double3 x_bh1, x_bh2, v_bh1, v_bh2; @@ -358,99 +338,11 @@ double3 a_bh1, adot_bh1, a_bh2, adot_bh2; #include "n_bh.c" -/* -int 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 - -return - 0 if everything OK -*/ - -// #endif // ADD_N_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]; 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}; #include "pn_bh_spin.c" -/* -int calc_force_pn_BH(double m1, double xx1[], double vv1[], double ss1[], - double m2, double xx2[], double vv2[], double ss2[], - double CCC_NB, double dt_bh, - int usedOrNot[], - double a_pn1[][3], double adot_pn1[][3], - double a_pn2[][3], double adot_pn2[][3]) -*/ - -/* - 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 -spin1[0,1,2] - normalized spin 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 -spin2[0,1,2] - normalized spin of the 2 BH - -CCC_NB - Speed of light "c" in internal units -dt_BH - timestep of the BH's, needed for the SPIN integration - -usedOrNot[PN0, PN1, PN2, PN2.5, PN3, PN3.5, SPIN] - different PN term usage: PN1, PN2, PN2.5, PN3, PN3.5, SPIN - 0 1 2 3 4 5 6 - - OUTPUT - -a_pn1 [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 1 BH -adot_pn1[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 1 BH - -a_pn2 [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 2 BH -adot_pn2[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 2 BH - -return - 0 if everything OK - - 505 if BH's separation < 4 x (RSwarch1 + RSwarch2) -*/ - -// #endif // ADD_N_BH - #ifdef ETICS double t_exp, dt_exp=ETICS_DTSCF; // t_exp is just the initial value #ifdef ETICS_CEP @@ -460,28 +352,25 @@ int grapite_cep_index; void get_CPU_time(double *time_real, double *time_user, double *time_syst) { -struct rusage xxx; -double sec_u, microsec_u, sec_s, microsec_s; -struct timeval tv; + struct rusage xxx; + double sec_u, microsec_u, sec_s, microsec_s; + struct timeval tv; -getrusage(RUSAGE_SELF,&xxx); + getrusage(RUSAGE_SELF,&xxx); -sec_u = xxx.ru_utime.tv_sec; -sec_s = xxx.ru_stime.tv_sec; + sec_u = xxx.ru_utime.tv_sec; + sec_s = xxx.ru_stime.tv_sec; -microsec_u = xxx.ru_utime.tv_usec; -microsec_s = xxx.ru_stime.tv_usec; + microsec_u = xxx.ru_utime.tv_usec; + microsec_s = xxx.ru_stime.tv_usec; -*time_user = sec_u + microsec_u * 1.0E-06; -*time_syst = sec_s + microsec_s * 1.0E-06; + *time_user = sec_u + microsec_u * 1.0E-06; + *time_syst = sec_s + microsec_s * 1.0E-06; -// *time_real = time(NULL); - -gettimeofday(&tv, NULL); -*time_real = tv.tv_sec + 1.0E-06 * tv.tv_usec; - -*time_user = *time_real; + gettimeofday(&tv, NULL); + *time_real = tv.tv_sec + 1.0E-06 * tv.tv_usec; + *time_user = *time_real; } void read_data(char inp_fname[], int *diskstep, int *N, double *time_cur, int ind[], double m[], double3 x[], double3 v[]) @@ -885,7 +774,7 @@ DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0); } -void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, double& rcm_sum, double& vcm_sum, double& E_tot_0, double& E_tot_corr_0, double& E_tot_corr_sd_0, int &skip_con, int N, double m[], double3 x[], double3 v[], double pot[]) +void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, double& rcm_sum, double& vcm_sum, double& E_tot_0, double& E_tot_corr_0, double& E_tot_corr_sd_0, int &skip_con, int N, double m[], double3 x[], double3 v[], double pot[], double pot_ext[]) { // TODO maybe use static variables here for the previous step energy? double E_pot = 0.0; @@ -1058,6 +947,21 @@ int main(int argc, char *argv[]) double C_NB; int usedOrNot[7]; + double pot_ext[N_MAX], pot_act_ext[N_MAX]; // for all EXTPOT + + int i_bh1, i_bh2; + + double m_bh1, m_bh2; + + int inf_event[N_MAX]; + double DR2, tmp_r2; + double3 x_bbhc, v_bbhc; + double DV2, EB, SEMI_a, SEMI_a2; + + double s_bh1[3] = {0.0, 0.0, 1.0}; + double s_bh2[3] = {0.0, 0.0, 1.0}; + + /* INIT the rand() !!! */ srand(19640916); /* it is just my birthday :-) */ @@ -1503,7 +1407,7 @@ int main(int argc, char *argv[]) /* Energy control... */ if (myRank == rootRank) { - energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot); + energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot, pot_ext); } /* if (myRank == rootRank) */ #ifdef ETICS_DUMP @@ -2130,7 +2034,7 @@ int main(int argc, char *argv[]) if (time_cur >= t_contr) { if (myRank == rootRank) { - energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot); + energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot, pot_ext); /* write cont data */