diff --git a/TODO.md b/TODO.md index 2cceccc..bf3f426 100644 --- a/TODO.md +++ b/TODO.md @@ -5,16 +5,4 @@ TODO * Memory bug when reading HDF5? x and v not allocated. -* Get rid of all global variables. Including config. - * Break main() into smaller chunks; operations that are timed should become independent functions. - -* Dynamically allocate the big arrays. - -* Create a key-value config file reader to replace the external parameter files. - -* Remove all ifdef blocks, options should be selected in the config file. - -* Optional HDF5 output. - -* In the config file let choose number of digits for ASCII output and 32 or 64-bit floats for HDF5. diff --git a/phigrape.cpp b/phigrape.cpp index 0d5b4dd..e889174 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -83,15 +83,28 @@ private: std::vector acc_loc, jrk_loc; }; -void calc_ext_grav(std::vector &external_gravity_components, int n, const std::vector &x, const std::vector &v, double *pot, double3 *acc, double3* jrk) -// TODO should just be a class that has this pointer array as a member -{ - std::fill(pot, pot+n, 0.); - for (auto component : external_gravity_components) { - if (component->is_active) - component->apply(n, x, v, pot, acc, jrk); +class Calc_ext_grav { +public: + void add_component(External_gravity &component) + { + components.push_back(&component); } -} + bool any_active() + { + for (auto component : components) + if (component->is_active) return true; + return false; + } + void operator()(int n, const std::vector &x, const std::vector &v, double *pot, double3 *acc, double3* jrk) + { + for (auto component : components) { + if (component->is_active) + component->apply(n, x, v, pot, acc, jrk); + } + } +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[]) { @@ -351,25 +364,18 @@ int main(int argc, char *argv[]) normalization_length = 1000/config.unit_length; normalization_velocity = 1.52484071426404437233e+01*sqrt(config.unit_length/config.unit_mass); } + Calc_ext_grav calc_ext_grav; Plummer ext_bulge(config.ext_m_bulge*normalization_mass, config.ext_b_bulge*normalization_length); + calc_ext_grav.add_component(ext_bulge); Miyamoto_Nagai ext_disk(config.ext_m_disk*normalization_mass, config.ext_a_disk*normalization_length, config.ext_b_disk*normalization_length); + calc_ext_grav.add_component(ext_disk); Plummer ext_halo_plummer(config.ext_m_halo_plummer*normalization_mass, config.ext_b_halo_plummer*normalization_length); + calc_ext_grav.add_component(ext_halo_plummer); Logarithmic_halo ext_log_halo(config.ext_log_halo_v*normalization_velocity, config.ext_log_halo_r*normalization_length); + calc_ext_grav.add_component(ext_log_halo); Dehnen ext_dehnen(config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma); - - std::vector external_gravity_components; - external_gravity_components.push_back(&ext_bulge); - external_gravity_components.push_back(&ext_disk); - external_gravity_components.push_back(&ext_halo_plummer); - external_gravity_components.push_back(&ext_log_halo); - external_gravity_components.push_back(&ext_dehnen); - bool has_external_gravity = false; - for (auto component : external_gravity_components) { - if (component->is_active) { - has_external_gravity = true; - break; - } - } + calc_ext_grav.add_component(ext_dehnen); + bool has_external_gravity = calc_ext_grav.any_active(); if (has_external_gravity) printf("External Potential: \n\n"); if (ext_bulge.is_active) printf("m_bulge = %.4E b_bulge = %.4E\n", config.ext_m_bulge*normalization_mass, config.ext_b_bulge*normalization_length); @@ -418,11 +424,10 @@ int main(int argc, char *argv[]) g6_set_tunit(51); g6_set_xunit(51); + bool grapite_active_search_flag = false; #ifdef ETICS grapite_set_dev_exec_threshold(config.grapite_dev_exec_threshold); - bool grapite_active_search_flag = config.grapite_active_search; -#else - bool grapite_active_search_flag = false; + grapite_active_search_flag = config.grapite_active_search; #endif int n_loc = N/n_proc; @@ -480,7 +485,8 @@ int main(int argc, char *argv[]) if (config.binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt[0], a[0], a[1], adot[0], adot[1]); } - calc_ext_grav(external_gravity_components, N, x, v, pot_ext, a, adot); + std::fill(pot_ext, pot_ext+N, 0.); + calc_ext_grav(N, x, v, pot_ext, a, adot); /* Energy control... */ if (myRank == rootRank) energy_contr(time_cur, timesteps, n_act_sum, calc_self_grav.g6_calls, N, m, x, v, pot, pot_ext); @@ -597,7 +603,8 @@ int main(int argc, char *argv[]) } /* Calculate gravity on active particles due to external forces */ - calc_ext_grav(external_gravity_components, n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new); + std::fill(pot_act_ext, pot_act_ext+n_act, 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 */ double min_dt = dt_max; @@ -706,7 +713,7 @@ int main(int argc, char *argv[]) // point. if (grapite_cep_index >= 0) { double3 xcm, vcm, xdc, vdc; - grapite_calc_center(N, m, (double(*)[3])x, (double(*)[3])v, xcm, vcm, xdc, vdc); + grapite_calc_center(N, m.data(), (double(*)[3])x.data(), (double(*)[3])v.data(), xcm, vcm, xdc, vdc); x[grapite_cep_index] = xdc; v[grapite_cep_index] = vdc; grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);