diff --git a/TODO.md b/TODO.md index 22d82ff..6effc24 100644 --- a/TODO.md +++ b/TODO.md @@ -3,4 +3,6 @@ TODO * Break main() into smaller chunks; operations that are timed should become independent functions. -* Const everything \ No newline at end of file +* Const everything + +* OpenMP \ No newline at end of file diff --git a/external.h b/external.h index 0702d09..ead7920 100644 --- a/external.h +++ b/external.h @@ -1,4 +1,5 @@ #pragma once +#include #include #include "double3.h" @@ -15,11 +16,17 @@ public: } } virtual void calc_gravity() = 0; - bool is_active; + virtual void print_info() {} + void set_name(const std::string &name) + { + this->name = name; + } + bool is_active = false; protected: double potential; double3 acceleration, jerk; double3 x, v; + std::string name = "ext"; void set_coordinates(double3 x, double3 v) { this->x = x; @@ -30,23 +37,38 @@ protected: class Plummer : public External_gravity { public: Plummer(double m, double b) : m(m), b(b) {is_active=(m>0);} - void calc_gravity(); + void calc_gravity() override; + void print_info() override + { + if (!is_active) return; + printf("m_%-5s = %.4E b_%-5s = %.4E\n", name.c_str(), m, name.c_str(), b); + } private: double m, b; }; class Miyamoto_Nagai : public External_gravity { public: - Miyamoto_Nagai(double m, double a, double b) : m(m), a(a), b(b) {is_active=(m>0);} + Miyamoto_Nagai(double m, double a, double b) : m(m), a(a), b(b) {is_active=(m>0); this->set_name("disk");} void calc_gravity(); + void print_info() override + { + if (!is_active) return; + printf("m_%-5s = %.4E a_%-5s = %.4E b_%-5s = %.4E\n", name.c_str(), m, name.c_str(), a, name.c_str(), b); + } private: double m, a, b; }; class Logarithmic_halo : public External_gravity { public: - Logarithmic_halo(double v_halo, double r_halo) : v2_halo(v_halo*v_halo), r2_halo(r_halo*r_halo) {is_active=(r_halo>0);} - void calc_gravity(); + Logarithmic_halo(double v_halo, double r_halo) : v2_halo(v_halo*v_halo), r2_halo(r_halo*r_halo) {is_active=(r_halo>0); this->set_name("halo");} + void calc_gravity() override; + void print_info() override + { + if (!is_active) return; + printf("v_%-4s = %.6E r_%-4s = %.4E\n", name.c_str(), sqrt(v2_halo), name.c_str(), sqrt(r2_halo)); + } private: double v2_halo, r2_halo; }; @@ -54,7 +76,12 @@ private: class Dehnen : public External_gravity { public: Dehnen(double m, double r, double gamma) : m(m), r(r), gamma(gamma) {is_active=(m>0);} - void calc_gravity(); + void calc_gravity() override; + void print_info() override + { + if (!is_active) return; + printf("m_%-5s = %.4E r_%-5s = %.4E g_%-5s = %.4E\n", name.c_str(), m, name.c_str(), r, name.c_str(), gamma); + } private: double m, r, gamma; }; diff --git a/phigrape.cpp b/phigrape.cpp index 1cab74e..b82cb12 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -85,6 +85,13 @@ public: component->apply(n, x, v, pot, acc, jrk); } } + void print_info() + { + for (auto component : components) { + component->print_info(); + } + fflush(stdout); + } private: std::vector components; }; @@ -249,13 +256,8 @@ inline double aarseth_step(const double eta, const double dt, const double3 &a, int main(int argc, char *argv[]) { - /* read the input parameters */ - const Config config("phigrape.conf"); - timer.start(); - double timesteps=0.0, n_act_sum=0.0; - /* INIT the rand() !!! */ srand(19640916); /* it is just my birthday :-) */ @@ -276,6 +278,7 @@ int main(int argc, char *argv[]) /* Print the Rank and the names of processors */ printf("Rank of the processor %03d on %s \n", myRank, processor_name); + const Config config("phigrape.conf"); Input_data input_data; if (is_hdf5(config.input_file_name)) { #ifndef HAS_HDF5 @@ -294,16 +297,20 @@ int main(int argc, char *argv[]) auto &x = input_data.x; auto &v = input_data.v; + double t_disk = config.dt_disk*(1+floor(time_cur/config.dt_disk)); + double t_contr = config.dt_contr*(1+floor(time_cur/config.dt_contr)); + double t_bh = config.dt_bh*(1+floor(time_cur/config.dt_bh)); + if (myRank == rootRank) { printf("\n"); printf("Begin the calculation of phi-GRAPE program on %03d processors\n", n_proc); printf("\n"); - printf("N = %07d \t eps = %.6E \n", N, config.eps); - printf("t_beg = %.6E \t t_end = %.6E \n", time_cur, config.t_end); - printf("dt_disk = %.6E \t dt_contr = %.6E \n", config.dt_disk, config.dt_contr); + printf("N = %07d \t eps = %.6E\n", N, config.eps); + printf("t_beg = %.6E \t t_end = %.6E\n", time_cur, config.t_end); + printf("dt_disk = %.6E \t dt_contr = %.6E\n", config.dt_disk, config.dt_contr); printf("dt_bh = %.6E \n", config.dt_bh); - printf("eta = %.6E \n", config.eta); - printf("\n"); + printf("eta = %.6E\n\n", config.eta); + printf("t_disk = %.6E t_contr = %.6E t_bh = %.6E\n\n", t_disk, t_contr, t_bh); if ((diskstep == 0) && (time_cur == 0)) { FILE *out = fopen("contr.dat", "w"); @@ -327,33 +334,18 @@ int main(int argc, char *argv[]) } Calc_ext_grav calc_ext_grav; Plummer ext_bulge(config.ext_m_bulge*normalization_mass, config.ext_b_bulge*normalization_length); + ext_bulge.set_name("bulge"); 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); + ext_halo_plummer.set_name("halo"); 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); 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); - 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\n", config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma); - fflush(stdout); - - double t_disk = config.dt_disk*(1+floor(time_cur/config.dt_disk)); - double t_contr = config.dt_contr*(1+floor(time_cur/config.dt_contr)); - double t_bh = config.dt_bh*(1+floor(time_cur/config.dt_bh)); - - if (myRank == rootRank) { - printf("t_disk = %.6E t_contr = %.6E t_bh = %.6E \n\n", t_disk, t_contr, t_bh); - fflush(stdout); - } /* if (myRank == rootRank) */ + if (myRank == rootRank) calc_ext_grav.print_info(); /* some local settings for G6a boards */ int clusterid, numGPU; @@ -439,6 +431,7 @@ int main(int argc, char *argv[]) std::vector pot_ext(N, 0.); calc_ext_grav(N, x, v, pot_ext, a, adot); + double timesteps=0, n_act_sum=0; /* 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);