Removed remaining macros and unused includes

TODO.md

@@ -6,5 +6,3 @@ TODO
 * Break main() into smaller chunks; operations that are timed should become independent functions.
 
 * Const everything
-
-* Remove unneeded includes

config.cpp

@@ -190,12 +190,16 @@ Config::Config(std::string file_name)
     eta = get_parameter<double>(dictionary, "eta");
     input_file_name = get_parameter<std::string>(dictionary, "input_file_name", "data.con");
     devices_per_node = get_parameter<int>(dictionary, "devices_per_node", 0);
-    dt_min_warning = get_parameter<bool>(dictionary, "dt_min_warning", false);
+    dt_max_power = get_parameter<int>(dictionary, "dt_max_power", -3);
+    dt_min_power = get_parameter<int>(dictionary, "dt_min_power", -36);
+    eta_s_corr = get_parameter<double>(dictionary, "eta_s_corr", 4);
+    eta_bh_corr = get_parameter<double>(dictionary, "eta_bh_corr", 4);
 
     output_hdf5 = get_parameter<bool>(dictionary, "output_hdf5", false);
     output_hdf5_double_precision = get_parameter<bool>(dictionary, "output_hdf5_double_precision", true);
     output_ascii_precision = get_parameter<int>(dictionary, "output_ascii_precision", 10);
     output_extra_mode = get_parameter<int>(dictionary, "output_extra_mode", 10);
+    dt_min_warning = get_parameter<bool>(dictionary, "dt_min_warning", false);
 
     live_smbh_count = get_parameter<int>(dictionary, "live_smbh_count", 0);
     live_smbh_custom_eps = get_parameter<double>(dictionary, "live_smbh_custom_eps", -1);

config.h

@@ -15,12 +15,16 @@ public:
     double              eta;
     std::string         input_file_name;
     int                 devices_per_node;
-    bool                dt_min_warning;
+    int                 dt_max_power;
+    int                 dt_min_power;
+    double              eta_s_corr;
+    double              eta_bh_corr;
    
     bool                output_hdf5;
     bool                output_hdf5_double_precision;
     int                 output_ascii_precision;
     int                 output_extra_mode;
+    bool                dt_min_warning;
    
     int                 live_smbh_count;
     double              live_smbh_custom_eps;

phigrape.conf

@@ -31,7 +31,25 @@ eta =                                   0.01
 # processes on a machine with a single device, set the value to 1 and use the
 # mpirun utility (or whatever is used in your job scheduler) to launch as many
 # processes as you like.
-#devices_per_node =                      1
+devices_per_node =                      1
+
+
+# The power of 2 of the maximum timestep. For example, if dt_max_power=-3 and
+# the adaptive timestep algotirhm suggests anything bigger than 0.125, the
+# timestep will be 0.125. If the 2^dt_max_power > dt_contr, the maximum timestep
+# is dt_contr. [default: -3]
+#dt_max_power =                          -3
+
+# The power of 2 of the minimum timestep. [default: -36]
+#dt_min_power =                          -36
+
+# Eta correction factor for the first step (the timestep will be calculated with
+# eta divided by eta_s_corr) [default: 4]
+#eta_s_corr =                            4.0
+
+# Eta correction factor for the black holes (the timestep for the black holes
+# will be calculated with eta divided by eta_bh_corr). [default: 4]
+#eta_bh_corr =                           4.0
 
 
 ##########

phigrape.cpp

@@ -1,21 +1,7 @@
-#define ETA_S_CORR   4.0
-#define ETA_BH_CORR  4.0
-
-#define DTMAXPOWER  -3.0
-#define DTMINPOWER -36.0
-
 #include <algorithm>
 #include <chrono>
-#include <math.h>
 #include <mpi.h>
 #include <numeric>
-#include <stdexcept>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <sys/time.h>
-#include <sys/resource.h>
-#include <time.h>
 
 #include "black_holes.h"
 #include "config.h"
@@ -28,22 +14,20 @@
 #include "grapite.h"
 #endif
 
-//chrono::steady_clock::time_point walltime_start;
-
 namespace std::chrono {
-struct Timer {
+    struct Timer {
-    void start()
+        void start()
-    {
+        {
-        t_start = steady_clock::now();
+            t_start = steady_clock::now();
-    }
+        }
-    void stop()
+        void stop()
-    {
+        {
-        t_stop = steady_clock::now();
+            t_stop = steady_clock::now();
-        time = duration_cast<nanoseconds>(t_stop - t_start).count()*1E-9;
+            time = duration_cast<nanoseconds>(t_stop - t_start).count()*1E-9;
-    }
+        }
-    double time; // seconds
+        double time; // seconds
-    steady_clock::time_point t_start, t_stop;
+        steady_clock::time_point t_start, t_stop;
-};
+    };
 }
 std::chrono::Timer timer;
 
@@ -306,22 +290,15 @@ int main(int argc, char *argv[])
     else
         ascii_read(config.input_file_name, diskstep, N, time_cur, m, x, v);
 
-    double eps = config.eps;
-    double eta = config.eta;
-    double t_end = config.t_end;
-    double dt_disk = config.dt_disk;
-    double dt_contr = config.dt_contr;
-    double dt_bh = 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, eps);
+        printf("N       = %07d \t eps      = %.6E \n", N, config.eps);
-        printf("t_beg   = %.6E \t t_end    = %.6E \n", time_cur, t_end);
+        printf("t_beg   = %.6E \t t_end    = %.6E \n", time_cur, config.t_end);
-        printf("dt_disk = %.6E \t dt_contr = %.6E \n", dt_disk, dt_contr);
+        printf("dt_disk = %.6E \t dt_contr = %.6E \n", config.dt_disk, config.dt_contr);
-        printf("dt_bh   = %.6E \n", dt_bh);
+        printf("dt_bh   = %.6E \n", config.dt_bh);
-        printf("eta     = %.6E \n", eta);
+        printf("eta     = %.6E \n", config.eta);
         printf("\n");
 
         if ((diskstep == 0) && (time_cur == 0)) {
@@ -365,18 +342,12 @@ int main(int argc, char *argv[])
     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;
+    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));
-    const double dt_min = pow(2, DTMINPOWER);
+    double t_bh    = config.dt_bh*(1+floor(time_cur/config.dt_bh));
-    const double dt_max = std::min({dt_disk, dt_contr, pow(2, DTMAXPOWER)});
-
-    double t_disk  = dt_disk*(1+floor(time_cur/dt_disk));
-    double t_contr = dt_contr*(1+floor(time_cur/dt_contr));
-    double t_bh    = dt_bh*(1+floor(time_cur/dt_bh));
 
     if (myRank == rootRank) {
-        printf("t_disk = %.6E   t_contr = %.6E   t_bh = %.6E \n", t_disk, t_contr, t_bh);
+        printf("t_disk = %.6E   t_contr = %.6E   t_bh = %.6E \n\n", t_disk, t_contr, t_bh);
-        printf("\n");
         fflush(stdout);
     } /* if (myRank == rootRank) */
 
@@ -413,6 +384,7 @@ int main(int argc, char *argv[])
     grapite_set_t_exp(time_cur);
 #endif
 
+    const double dt_min = pow(2, config.dt_min_power);
     std::vector<int> ind(N);
     std::iota(begin(ind), end(ind), 0);
     /* load the nj particles to the G6 */
@@ -442,7 +414,7 @@ int main(int argc, char *argv[])
     std::vector<double> 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 calc_self_grav(N, n_loc, clusterid, npipe, config.eps);
     calc_self_grav(time_cur, N, ind, x, v, pot, a, adot);
 
     Black_hole_physics black_hole_physics;
@@ -458,7 +430,7 @@ int main(int argc, char *argv[])
         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);
+    black_hole_physics.set_softening(config.eps, config.live_smbh_custom_eps);
 
     std::vector<double> pot_ext(N, 0.);
     calc_ext_grav(N, x, v, pot_ext, a, adot);
@@ -482,13 +454,14 @@ int main(int argc, char *argv[])
     }
 #endif
 
+    const double dt_max = std::min({config.dt_disk, config.dt_contr, pow(2, config.dt_max_power)});
     std::vector<double> dt(N);
     /* Define initial timestep for all particles on all nodes */
     for (int j=0; j<N; j++) {
         double a2_mod = a[j].norm2();
         double adot2_mod = adot[j].norm2();
 
-        double dt_tmp, eta_s = eta/ETA_S_CORR;
+        double dt_tmp, eta_s = config.eta/config.eta_s_corr;
         if (adot2_mod==0) dt_tmp = eta_s; // That's weird, when will we have such a case?
         else              dt_tmp = eta_s*sqrt(a2_mod/adot2_mod);
 
@@ -538,7 +511,7 @@ int main(int argc, char *argv[])
     } /* if (myRank == rootRank) */
 
     /* The main integration loop */
-    while (time_cur <= t_end) {
+    while (time_cur <= config.t_end) {
 
         /* Define the minimal time and the active particles on all the nodes */
         double min_t = active_search.get_minimum_time(t, dt);
@@ -601,8 +574,8 @@ int main(int argc, char *argv[])
 
             //TODO make beautiful
             double eta_curr;
-            if ((config.live_smbh_count > 0) && (ind_act[i] < config.live_smbh_count)) eta_curr = eta_bh;
+            if ((config.live_smbh_count > 0) && (ind_act[i] < config.live_smbh_count)) eta_curr = config.eta/config.eta_bh_corr;
-            else eta_curr = eta;
+            else eta_curr = config.eta;
 
             double dt_new = aarseth_step(eta_curr, dt_tmp, a_act_new[i], adot_act_new[i], a2, a3);
 
@@ -675,7 +648,7 @@ int main(int argc, char *argv[])
 
             } /* if (myRank == rootRank) */
 
-            t_bh += dt_bh;
+            t_bh += config.dt_bh;
         } /* if (time_cur >= t_bh) */
 
         if (time_cur >= t_contr) {
@@ -702,7 +675,7 @@ int main(int argc, char *argv[])
             }
 #endif
 
-            t_contr += dt_contr;
+            t_contr += config.dt_contr;
         } /* if (time_cur >= t_contr) */
 
         if (time_cur >= t_disk) {
@@ -720,7 +693,7 @@ int main(int argc, char *argv[])
                 grapite_dump(out_fname, 2);
             }
 #endif
-            t_disk += dt_disk;
+            t_disk += config.dt_disk;
         } /* if (time_cur >= t_disk) */
     } /* while (time_cur < t_end) */