From d9ed8e5131b1048edffbfd666e31c87b8c335a18 Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Wed, 28 Oct 2020 19:18:37 -0400 Subject: [PATCH] Got rid of remaining manual allocations, moved some variable declarations to where they are needed --- TODO.md | 6 ++- black_holes.cpp | 2 +- black_holes.h | 8 ++- external.cpp | 1 - external.h | 2 +- io.cpp | 8 +-- io.h | 2 +- phigrape.cpp | 128 ++++++++++++++++++++---------------------------- 8 files changed, 67 insertions(+), 90 deletions(-) diff --git a/TODO.md b/TODO.md index bf3f426..2c926f8 100644 --- a/TODO.md +++ b/TODO.md @@ -1,8 +1,10 @@ TODO ==== -* Make the reamining arrays vectors or smart pointers. - * Memory bug when reading HDF5? x and v not allocated. * Break main() into smaller chunks; operations that are timed should become independent functions. + +* Const everything + +* Remove unneeded includes \ No newline at end of file diff --git a/black_holes.cpp b/black_holes.cpp index c14b970..7790a3b 100644 --- a/black_holes.cpp +++ b/black_holes.cpp @@ -99,7 +99,7 @@ void Black_hole_physics::adjust_post_newtonian( jrk2 += jrk2_corr; } -void Black_hole_physics::write_bh_data(double time_cur, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, double3 a[], double3 adot[], double dt[]) +void Black_hole_physics::write_bh_data(double time_cur, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &a, const std::vector &adot, const std::vector &dt) { // This function logs data on the black hole(s). It uses both external data // (the arguments to this function) and optionall internal data to this diff --git a/black_holes.h b/black_holes.h index 402a628..121ae7a 100644 --- a/black_holes.h +++ b/black_holes.h @@ -46,7 +46,7 @@ public: const double dt_bh, // pn_usage should be const double3& acc1, double3& acc2, double3& jrk1, double3& jrk2); - void write_bh_data(double time_cur, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, double3 a[], double3 adot[], double dt[]); + void write_bh_data(double time_cur, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &a, const std::vector &adot, const std::vector &dt); public: //TODO make private double m1, m2; int count; @@ -83,7 +83,7 @@ private: class Binary_smbh_influence_sphere_output { public: - Binary_smbh_influence_sphere_output(double factor, int N, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, double *dt) + Binary_smbh_influence_sphere_output(double factor, int N, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &dt) : factor(factor), m(m), x(x), v(v), pot(pot), dt(dt) { inf_event.assign(N, 0); @@ -97,9 +97,7 @@ public: void operator()(const std::vector& ind_act, int n_act, double timesteps, double time_cur); private: double factor; - const std::vector &pot; - const std::vector &m; - double /**pot,*/ *dt; + const std::vector &pot, &m, &dt; const std::vector &x, &v; std::vector inf_event; FILE *out; diff --git a/external.cpp b/external.cpp index 27c5f71..3e927ad 100644 --- a/external.cpp +++ b/external.cpp @@ -59,7 +59,6 @@ void Logarithmic_halo::calc_gravity() jerk = tmp * (rv_ij * x - r2_r2_halo * v); } - void Dehnen::calc_gravity() { diff --git a/external.h b/external.h index 08d989f..0702d09 100644 --- a/external.h +++ b/external.h @@ -4,7 +4,7 @@ class External_gravity { public: - void apply(const int n_act, const std::vector &x, const std::vector &v, double pot[], double3 a[], double3 adot[]) + void apply(const int n_act, const std::vector &x, const std::vector &v, std::vector &pot, std::vector &a, std::vector &adot) { for (int i=0; iset_coordinates(x[i], v[i]); diff --git a/io.cpp b/io.cpp index 6bc8e5f..d0a187d 100644 --- a/io.cpp +++ b/io.cpp @@ -146,7 +146,7 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub #endif } -void h5_write(const std::string file_name, const int step_num, const int N, const double t, std::vector &m, std::vector &x, std::vector &v, const std::vector &pot, const double3 *acc, const double3 *jrk, const int extra_mode=0, const bool use_double_precision=true) +void h5_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &acc, const std::vector &jrk, const int extra_mode=0, const bool use_double_precision=true) { #ifdef HAS_HDF5 hid_t file_id, group_id, attribute_id, dataspace_id; @@ -174,7 +174,7 @@ void h5_write(const std::string file_name, const int step_num, const int N, cons H5Sclose(dataspace_id); }; - write_dataset("MASS", 1, m.data()); + write_dataset("MASS", 1, (double*)m.data()); // casting away const... write_dataset("X", 2, (double*)x.data()); write_dataset("V", 2, (double*)v.data()); @@ -182,8 +182,8 @@ void h5_write(const std::string file_name, const int step_num, const int N, cons bool write_acc = (extra_mode >> 1) & 1; bool write_jrk = (extra_mode >> 2) & 1; if (write_pot) write_dataset("POT", 1, (double*)pot.data()); - if (write_acc) write_dataset("ACC", 2, (double*)acc); - if (write_jrk) write_dataset("JRK", 2, (double*)jrk); + if (write_acc) write_dataset("ACC", 2, (double*)acc.data()); + if (write_jrk) write_dataset("JRK", 2, (double*)jrk.data()); H5Gclose(group_id); H5Fclose(file_id); diff --git a/io.h b/io.h index 874ca27..7571592 100644 --- a/io.h +++ b/io.h @@ -12,5 +12,5 @@ void ascii_write(const std::string file_name, const int step_num, const int N, c void h5_read(const std::string file_name, int *step_num, int *N, double *t, double m[], double3 x[], double3 v[]); // In case the code is compiled without HDF5 support, the implementation of this function just throws an error -void h5_write(const std::string file_name, const int step_num, const int N, const double t, std::vector &m, std::vector &x, std::vector &v, const std::vector &pot, const double3 *acc, const double3 *jrk, const int extra_mode=0, const bool use_double_precision=true); +void h5_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &acc, const std::vector &jrk, const int extra_mode=0, const bool use_double_precision=true); // In case the code is compiled without HDF5 support, the implementation of this function just throws an error diff --git a/phigrape.cpp b/phigrape.cpp index e889174..bab4157 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -58,21 +58,20 @@ public: jrk_loc.resize(N); } void operator()(const double t, const int n_act, std::vector &ind_act, std::vector &x_act, std::vector &v_act, - std::vector& pot, double3 acc[], double3 jrk[]) + std::vector& pot, std::vector &acc, std::vector &jrk) { g6_set_ti(clusterid, t); for (int i=0; iis_active) return true; return false; } - void operator()(int n, const std::vector &x, const std::vector &v, double *pot, double3 *acc, double3* jrk) + void operator()(int n, const std::vector &x, const std::vector &v, std::vector &pot, std::vector &acc, std::vector &jrk) { for (auto component : components) { if (component->is_active) @@ -106,7 +105,7 @@ private: std::vector components; }; -void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, int N, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, double pot_ext[]) +void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, int N, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, const std::vector &pot_ext) { double E_pot = 0; for (int i=0; i &t, std::vector &dt) { double min_t_loc, min_t; #ifdef ETICS @@ -197,7 +195,7 @@ public: MPI_Allreduce(&min_t_loc, &min_t, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); return min_t; } - void get_active_indices(const double min_t, const double t[], const double dt[], int ind_act[], int& n_act) + void get_active_indices(const double min_t, const std::vector &t, const std::vector &dt, std::vector &ind_act, int &n_act) { #ifdef ETICS if (grapite_active_search_flag) { @@ -226,11 +224,11 @@ public: } private: int myRank, n_proc, n_loc, N; - int *ind_act_loc; + std::vector ind_act_loc; bool grapite_active_search_flag; }; -inline void calc_high_derivatives(const double dt_tmp, const double3 a_old, const double3 a_new, const double3 a1_old, const double3 a1_new, double3& a2, double3& a3) +inline void calc_high_derivatives(const double dt_tmp, const double3 &a_old, const double3 &a_new, const double3 &a1_old, const double3 &a1_new, double3 &a2, double3 &a3) { double dtinv = 1/dt_tmp; double dt2inv = dtinv*dtinv; @@ -244,7 +242,7 @@ inline void calc_high_derivatives(const double dt_tmp, const double3 a_old, cons a3 = 12*a0mia1*dt3inv + 6*ad0plad1*dt2inv; } -inline void corrector(const double dt_tmp, const double3 a2, const double3 a3, double3& x, double3& v) +inline void corrector(const double dt_tmp, const double3 &a2, const double3 &a3, double3 &x, double3 &v) { double dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0; double dt4over24 = dt3over6*dt_tmp/4.0; @@ -254,7 +252,7 @@ inline void corrector(const double dt_tmp, const double3 a2, const double3 a3, d v += dt3over6*a2 + dt4over24*a3; } -inline double aarseth_step(const double eta, const double dt, const double3 a, const double3 a1, const double3 a2, const double3 a3) +inline double aarseth_step(const double eta, const double dt, const double3 &a, const double3 &a1, const double3 &a2, const double3 &a3) { double a1abs = a.norm(); double adot1abs = a1.norm(); @@ -274,8 +272,6 @@ int main(int argc, char *argv[]) double timesteps=0.0, n_act_sum=0.0; - double3 zeros = {0, 0, 0}; // Dummy; can't really be const because of the GRAPE interface. - /* INIT the rand() !!! */ srand(19640916); /* it is just my birthday :-) */ @@ -298,8 +294,6 @@ int main(int argc, char *argv[]) int diskstep, N; double time_cur; - // The memory for m, x, and v is allocated inside h5_read or ascii_read - //double *m; std::vector m; std::vector x, v; if (is_hdf5(config.input_file_name)) { @@ -312,20 +306,6 @@ int main(int argc, char *argv[]) else ascii_read(config.input_file_name, diskstep, N, time_cur, m, x, v); - std::vector ind(N); - std::iota(begin(ind), end(ind), 0); - double3 *a = new double3[N], *adot = new double3[N]; - std::vector pot(N); - double *pot_ext = new double[N], *t = new double[N], *dt = new double[N]; - - /* data for active particles */ - int n_act; - std::vector ind_act(N); - std::vector pot_act_new(N); - std::vector x_act_new(N), v_act_new(N); - double *pot_act_ext = new double[N]; - double3 *a_act_new = new double3[N], *adot_act_new = new double3[N]; - double eps = config.eps; double eta = config.eta; double t_end = config.t_end; @@ -382,8 +362,7 @@ int main(int argc, char *argv[]) if (ext_disk.is_active) printf("m_disk = %.4E a_disk = %.4E b_disk = %.4E\n", config.ext_m_disk*normalization_mass, config.ext_a_disk*normalization_length, config.ext_b_disk*normalization_length); if (ext_halo_plummer.is_active) printf("m_halo = %.4E b_halo = %.4E\n", config.ext_m_halo_plummer*normalization_mass, config.ext_b_halo_plummer*normalization_length); if (ext_log_halo.is_active) printf("v_halo = %.6E r_halo = %.6E \n", config.ext_log_halo_v*normalization_velocity, config.ext_log_halo_r*normalization_length); - if (ext_dehnen.is_active) printf("m_ext = %.6E r_ext = %.6E \t g_ext = %.3E \n", config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma); - printf("\n"); + if (ext_dehnen.is_active) printf("m_ext = %.6E r_ext = %.6E \t g_ext = %.3E \n\n", config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma); fflush(stdout); double eta_bh = eta/ETA_BH_CORR; @@ -401,9 +380,6 @@ int main(int argc, char *argv[]) fflush(stdout); } /* if (myRank == rootRank) */ - std::fill(t, t+N, time_cur); - std::fill(dt, dt+N, dt_min); - /* some local settings for G6a boards */ int clusterid, numGPU; if (config.devices_per_node==0) { @@ -431,33 +407,19 @@ int main(int argc, char *argv[]) #endif int n_loc = N/n_proc; - Calc_self_grav calc_self_grav(N, n_loc, clusterid, npipe, eps); - Active_search active_search(myRank, n_proc, n_loc, N, grapite_active_search_flag); - Black_hole_physics black_hole_physics; - if (config.live_smbh_count == 1) - black_hole_physics = Black_hole_physics(m[0], 0, myRank, rootRank); - else if (config.live_smbh_count == 2) - black_hole_physics = Black_hole_physics(m[0], m[1], myRank, rootRank); - if (config.binary_smbh_pn) { - black_hole_physics.set_post_newtonian(config.pn_c, config.pn_usage.data()); - if (config.pn_usage[6]) black_hole_physics.set_spins(config.smbh1_spin.data(), config.smbh2_spin.data()); - } - black_hole_physics.set_softening(eps, config.live_smbh_custom_eps); - - Binary_smbh_influence_sphere_output binary_smbh_influence_sphere_output(config.binary_smbh_influence_radius_factor, N, m, x, v, pot, dt); - - Write_bh_nb_data write_bh_nb_data(config.live_smbh_neighbor_number, config.live_smbh_count, N, m, x, v); - #ifdef ETICS grapite_read_particle_tags(N, config.grapite_mask_file_name.c_str(), myRank, n_loc); grapite_set_dt_exp(config.dt_scf); grapite_set_t_exp(time_cur); #endif + std::vector ind(N); + std::iota(begin(ind), end(ind), 0); /* load the nj particles to the G6 */ + double3 zeros = {0, 0, 0}; // Dummy; can't really be const because of the GRAPE interface. for (int k=0; k a(N), adot(N); + std::vector pot(N); + /* define the all particles as a active on all the processors for the first time grav calc. */ + Calc_self_grav calc_self_grav(N, n_loc, clusterid, npipe, eps); calc_self_grav(time_cur, N, ind, x, v, pot, a, adot); - if (config.live_smbh_count == 2) { + Black_hole_physics black_hole_physics; + if (config.live_smbh_count == 1) + black_hole_physics = Black_hole_physics(m[0], 0, myRank, rootRank); + else if (config.live_smbh_count == 2) { + black_hole_physics = Black_hole_physics(m[0], m[1], myRank, rootRank); black_hole_physics.set_xv(x[0], x[1], v[0], v[1]); if (config.live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(pot[0], pot[1], a[0], a[1], adot[0], adot[1]); - if (config.binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt[0], a[0], a[1], adot[0], adot[1]); + if (config.binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt_min, a[0], a[1], adot[0], adot[1]); } + if (config.binary_smbh_pn) { + black_hole_physics.set_post_newtonian(config.pn_c, config.pn_usage.data()); + if (config.pn_usage[6]) black_hole_physics.set_spins(config.smbh1_spin.data(), config.smbh2_spin.data()); + } + black_hole_physics.set_softening(eps, config.live_smbh_custom_eps); - std::fill(pot_ext, pot_ext+N, 0.); + std::vector pot_ext(N, 0.); calc_ext_grav(N, x, v, pot_ext, a, adot); /* Energy control... */ @@ -507,6 +482,7 @@ int main(int argc, char *argv[]) } #endif + std::vector dt(N); /* Define initial timestep for all particles on all nodes */ for (int j=0; j 0) { - double min_dt = *std::min_element(dt, dt+N); + double min_dt = *std::min_element(begin(dt), end(dt)); for (int i=0; i ind_act(N); + std::vector x_act_new(N), v_act_new(N), a_act_new(N), adot_act_new(N); + std::vector t(N, time_cur), pot_act_new(N); + + // Functors for the main integration loop + Active_search active_search(myRank, n_proc, n_loc, N, grapite_active_search_flag); + Binary_smbh_influence_sphere_output binary_smbh_influence_sphere_output(config.binary_smbh_influence_radius_factor, N, m, x, v, pot, dt); + Write_bh_nb_data write_bh_nb_data(config.live_smbh_neighbor_number, config.live_smbh_count, N, m, x, v); + if (myRank == rootRank) { + if (config.live_smbh_output) black_hole_physics.write_bh_data(time_cur, m, x, v, pot, a, adot, dt); + if (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur); + } /* if (myRank == rootRank) */ + /* The main integration loop */ while (time_cur <= t_end) { @@ -563,7 +544,8 @@ int main(int argc, char *argv[]) double min_t = active_search.get_minimum_time(t, dt); /* Get indices of all particles that will be active in this bunch */ - active_search.get_active_indices(min_t, t, dt, ind_act.data(), n_act); + int n_act; + active_search.get_active_indices(min_t, t, dt, ind_act, n_act); /* Find the BH(s) indices in the active list */ int i_bh1=0, i_bh2=1; @@ -603,7 +585,7 @@ int main(int argc, char *argv[]) } /* Calculate gravity on active particles due to external forces */ - std::fill(pot_act_ext, pot_act_ext+n_act, 0.); + std::vector pot_act_ext(N, 0.); calc_ext_grav(n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new); /* correct the active particles positions etc... on all the nodes */ @@ -757,10 +739,6 @@ int main(int argc, char *argv[]) fflush(stdout); } /* if (myRank == rootRank) */ - delete[] a; delete[] adot; delete[] pot_ext; delete[] t; delete[] dt; delete[] a_act_new; delete[] adot_act_new; delete[] pot_act_ext; - /* Finalize the MPI work */ MPI_Finalize(); - - return 0; }