Improved external gravity calls style

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.

phigrape.cpp

@@ -83,15 +83,28 @@ private:
     std::vector<double3> acc_loc, jrk_loc;
 };
 
-void calc_ext_grav(std::vector<External_gravity*> &external_gravity_components, int n, const std::vector<double3> &x, const std::vector<double3> &v, double *pot, double3 *acc, double3* jrk)
+class Calc_ext_grav {
-// TODO should just be a class that has this pointer array as a member
+public:
-{
+    void add_component(External_gravity &component)
-    std::fill(pot, pot+n, 0.);
+    {
-    for (auto component : external_gravity_components) {
+        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<double3> &x, const std::vector<double3> &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<External_gravity*> 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<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &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);
-
+    calc_ext_grav.add_component(ext_dehnen);
-    std::vector<External_gravity*> external_gravity_components;
+    bool has_external_gravity = calc_ext_grav.any_active();
-    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;
-        }
-    }
 
     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;
+    grapite_active_search_flag = config.grapite_active_search;
-#else
-    bool grapite_active_search_flag = false;
 #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]);