Improved the way of computing the exact accelerations by removing the need of...

Improved the way of computing the exact accelerations by removing the need of recalculating exact forces that already exist

.gitignore

@@ -32,6 +32,7 @@ examples/*/*/*.hdf5
 examples/*/*/*.png
 examples/*/*/*.txt
 examples/*/*/used_parameters.yml
+examples/*/gravity_checks_*.dat
 
 tests/testPair
 tests/brute_force_standard.dat

examples/plot_gravity_checks.py

@@ -45,13 +45,25 @@ cols = ['b', 'g', 'r', 'm']
 step = int(sys.argv[1])
 
 # Find the files for the different expansion orders
-order_list = glob.glob("gravity_checks_step%d_order*.dat"%step)
+order_list = glob.glob("gravity_checks_swift_step%d_order*.dat"%step)
 num_order = len(order_list)
 
 # Get the multipole orders
 order = np.zeros(num_order)
 for i in range(num_order):
-    order[i] = int(order_list[i][26])
+    order[i] = int(order_list[i][32])
+
+# Read the exact accelerations first
+data = np.loadtxt('gravity_checks_exact_step%d.dat'%step)
+exact_ids = data[:,0]
+exact_pos = data[:,1:4]
+exact_a = data[:,4:7]
+# Sort stuff
+sort_index = np.argsort(exact_ids)
+exact_ids = exact_ids[sort_index]
+exact_pos = exact_pos[sort_index, :]
+exact_a = exact_a[sort_index, :]        
+exact_a_norm = np.sqrt(exact_a[:,0]**2 + exact_a[:,1]**2 + exact_a[:,2]**2)
     
 # Start the plot
 plt.figure()
@@ -72,6 +84,23 @@ if len(gadget2_file_list) != 0:
     gadget2_pos = gadget2_pos[sort_index, :]
     gadget2_a_exact = gadget2_a_exact[sort_index, :]
     gadget2_a_grav = gadget2_a_grav[sort_index, :]
+
+    # Cross-checks
+    if not np.array_equal(exact_ids, gadget2_ids):
+        print "Comparing different IDs !"
+
+    if np.max(np.abs(exact_pos - gadget2_pos)/np.abs(gadget2_pos)) > 1e-6:
+        print "Comparing different positions ! max difference:"
+        index = np.argmax(exact_pos[:,0]**2 + exact_pos[:,1]**2 + exact_pos[:,2]**2 - gadget2_pos[:,0]**2 - gadget2_pos[:,1]**2 - gadget2_pos[:,2]**2)
+        print "Gadget2 (id=%d):"%gadget2_ids[index], gadget2_pos[index,:], "exact (id=%d):"%exact_ids[index], exact_pos[index,:], "\n"
+
+    if np.max(np.abs(a_exact - gadget2_a_exact) / np.abs(gadget2_a_exact)) > 2e-6:
+        print "Comparing different exact accelerations ! max difference:"
+        index = np.argmax(exact_a[:,0]**2 + exact_a[:,1]**2 + exact_a[:,2]**2 - gadget2_a_exact[:,0]**2 - gadget2_a_exact[:,1]**2 - gadget2_a_exact[:,2]**2)
+        print "a_exact --- Gadget2:", gadget2_a_exact[index,:], "exact:", a_exact[index,:]
+        print "a_grav ---  Gadget2:", gadget2_a_grav[index,:], "exact:", a_grav[index,:],"\n"
+        print "pos ---     Gadget2: (id=%d):"%gadget2_ids[index], gadget2_pos[index,:], "exact (id=%d):"%ids[index], pos[index,:],"\n"
+
     
     # Compute the error norm
     diff = gadget2_a_exact - gadget2_a_grav
@@ -123,43 +152,33 @@ for i in range(num_order-1, -1, -1):
     data = np.loadtxt(order_list[i])
     ids = data[:,0]
     pos = data[:,1:4]
-    a_exact = data[:,4:7]
-    a_grav = data[:, 7:10]
+    a_grav = data[:, 4:7]
 
     # Sort stuff
     sort_index = np.argsort(ids)
     ids = ids[sort_index]
     pos = pos[sort_index, :]
-    a_exact = a_exact[sort_index, :]
-    a_grav = a_grav[sort_index, :]
+    a_grav = a_grav[sort_index, :]        
 
     # Cross-checks
-    if not np.array_equal(ids, gadget2_ids):
+    if not np.array_equal(exact_ids, ids):
         print "Comparing different IDs !"
 
-    if np.max(np.abs(pos - gadget2_pos)/np.abs(gadget2_pos)) > 1e-6:
+    if np.max(np.abs(exact_pos - pos)/np.abs(pos)) > 1e-6:
         print "Comparing different positions ! max difference:"
-        index = np.argmax(pos[:,0]**2 + pos[:,1]**2 + pos[:,2]**2 - gadget2_pos[:,0]**2 - gadget2_pos[:,1]**2 - gadget2_pos[:,2]**2)
-        print "Gadget2 (id=%d):"%gadget2_ids[index], gadget2_pos[index,:], "SWIFT (id=%d):"%ids[index], pos[index,:], "\n"
+        index = np.argmax(exact_pos[:,0]**2 + exact_pos[:,1]**2 + exact_pos[:,2]**2 - pos[:,0]**2 - pos[:,1]**2 - pos[:,2]**2)
+        print "SWIFT (id=%d):"%ids[index], pos[index,:], "exact (id=%d):"%exact_ids[index], exact_pos[index,:], "\n"
 
-    if np.max(np.abs(a_exact - gadget2_a_exact) / np.abs(gadget2_a_exact)) > 2e-6:
-        print "Comparing different exact accelerations ! max difference:"
-        index = np.argmax(a_exact[:,0]**2 + a_exact[:,1]**2 + a_exact[:,2]**2 - gadget2_a_exact[:,0]**2 - gadget2_a_exact[:,1]**2 - gadget2_a_exact[:,2]**2)
-        print "a_exact --- Gadget2:", gadget2_a_exact[index,:], "SWIFT:", a_exact[index,:]
-        print "a_grav ---  Gadget2:", gadget2_a_grav[index,:], "SWIFT:", a_grav[index,:],"\n"
-        print "pos ---     Gadget2: (id=%d):"%gadget2_ids[index], gadget2_pos[index,:], "SWIFT (id=%d):"%ids[index], pos[index,:],"\n"
-        
-        
+    
     # Compute the error norm
-    diff = a_exact - a_grav
+    diff = exact_a - a_grav
 
     norm_diff = np.sqrt(diff[:,0]**2 + diff[:,1]**2 + diff[:,2]**2)
-    norm_a = np.sqrt(a_exact[:,0]**2 + a_exact[:,1]**2 + a_exact[:,2]**2)
 
-    norm_error = norm_diff / norm_a
-    error_x = abs(diff[:,0]) / norm_a
-    error_y = abs(diff[:,1]) / norm_a
-    error_z = abs(diff[:,2]) / norm_a
+    norm_error = norm_diff / exact_a_norm
+    error_x = abs(diff[:,0]) / exact_a_norm
+    error_y = abs(diff[:,1]) / exact_a_norm
+    error_z = abs(diff[:,2]) / exact_a_norm
     
     # Bin the error
     norm_error_hist,_ = np.histogram(norm_error, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)

src/gravity.c

@@ -22,6 +22,7 @@
 
 /* Some standard headers. */
 #include <stdio.h>
+#include <unistd.h>
 
 /* This object's header. */
 #include "gravity.h"
@@ -30,6 +31,50 @@
 #include "active.h"
 #include "error.h"
 
+/**
+ * @brief Checks whether the file containing the exact accelerations for
+ * the current choice of parameters already exists.
+ *
+ * @param e The #engine.
+ */
+int gravity_exact_force_file_exits(const struct engine *e) {
+
+  /* File name */
+  char file_name[100];
+  sprintf(file_name, "gravity_checks_exact_step%d.dat", e->step);
+
+  /* Does the file exist ? */
+  if (access(file_name, R_OK | W_OK) == 0) {
+
+    /* Let's check whether the header matches the parameters of this run */
+    FILE *file = fopen(file_name, "r");
+    if (!file) error("Problem reading gravity_check file");
+
+    char line[100];
+    char dummy1[10], dummy2[10];
+    double epsilon, newton_G;
+    int N;
+    /* Reads file header */
+    if (fgets(line, 100, file) != line) error("Problem reading title");
+    if (fgets(line, 100, file) != line) error("Problem reading G");
+    sscanf(line, "%s %s %le", dummy1, dummy2, &newton_G);
+    if (fgets(line, 100, file) != line) error("Problem reading N");
+    sscanf(line, "%s %s %d", dummy1, dummy2, &N);
+    if (fgets(line, 100, file) != line) error("Problem reading epsilon");
+    sscanf(line, "%s %s %le", dummy1, dummy2, &epsilon);
+    fclose(file);
+
+    /* Check whether it matches the current parameters */
+    if (N == SWIFT_GRAVITY_FORCE_CHECKS &&
+        (fabs(epsilon - e->gravity_properties->epsilon) / epsilon < 1e-5) &&
+        (fabs(newton_G - e->physical_constants->const_newton_G) / newton_G <
+         1e-5)) {
+      return 1;
+    }
+  }
+  return 0;
+}
+
 /**
  * @brief Run a brute-force gravity calculation for a subset of particles.
  *
@@ -47,6 +92,13 @@ void gravity_exact_force_compute(struct space *s, const struct engine *e) {
   const double const_G = e->physical_constants->const_newton_G;
   int counter = 0;
 
+  /* Let's start by checking whether we already computed these forces */
+  if (gravity_exact_force_file_exits(e)) {
+    message("Exact accelerations already computed. Skipping calculation.");
+    return;
+  }
+
+  /* No matching file present ? Do it then */
   for (size_t i = 0; i < s->nr_gparts; ++i) {
 
     struct gpart *gpi = &s->gparts[i];
@@ -140,25 +192,22 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
 
-  int counter = 0;
-
-  /* Some accumulators */
-  float err_rel_max = 0.f;
-  float err_rel_min = FLT_MAX;
-  float err_rel_mean = 0.f;
-  float err_rel_mean2 = 0.f;
-  float err_rel_std = 0.f;
-
-  char file_name[100];
-  sprintf(file_name, "gravity_checks_step%d_order%d.dat", e->step,
+  /* File name */
+  char file_name_swift[100];
+  sprintf(file_name_swift, "gravity_checks_swift_step%d_order%d.dat", e->step,
           SELF_GRAVITY_MULTIPOLE_ORDER);
-  FILE *file = fopen(file_name, "w");
-  fprintf(file, "# Gravity accuracy test G = %16.8e\n",
-          e->physical_constants->const_newton_G);
-  fprintf(file, "# %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s\n", "id",
-          "pos[0]", "pos[1]", "pos[2]", "a_exact[0]", "a_exact[1]",
-          "a_exact[2]", "a_grav[0]", "a_grav[1]", "a_grav[2]");
 
+  /* Creare files and write header */
+  FILE *file_swift = fopen(file_name_swift, "w");
+  fprintf(file_swift, "# Gravity accuracy test - SWIFT FORCES\n");
+  fprintf(file_swift, "# G= %16.8e\n", e->physical_constants->const_newton_G);
+  fprintf(file_swift, "# N= %d\n", SWIFT_GRAVITY_FORCE_CHECKS);
+  fprintf(file_swift, "# epsilon=%16.8e\n", e->gravity_properties->epsilon);
+  fprintf(file_swift, "# theta=%16.8e\n", e->gravity_properties->theta_crit);
+  fprintf(file_swift, "# %16s %16s %16s %16s %16s %16s %16s\n", "id", "pos[0]",
+          "pos[1]", "pos[2]", "a_swift[0]", "a_swift[1]", "a_swift[2]");
+
+  /* Output particle SWIFT accelerations  */
   for (size_t i = 0; i < s->nr_gparts; ++i) {
 
     struct gpart *gpi = &s->gparts[i];
@@ -167,64 +216,53 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
     if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
         gpart_is_starting(gpi, e)) {
 
-      fprintf(file,
-              "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e "
-              "%16.8e \n",
+      fprintf(file_swift, "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e \n",
               gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
-              gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2],
               gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2]);
-
-      const float diff[3] = {gpi->a_grav[0] - gpi->a_grav_exact[0],
-                             gpi->a_grav[1] - gpi->a_grav_exact[1],
-                             gpi->a_grav[2] - gpi->a_grav_exact[2]};
-
-      const float diff_norm =
-          sqrtf(diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]);
-      const float a_norm = sqrtf(gpi->a_grav_exact[0] * gpi->a_grav_exact[0] +
-                                 gpi->a_grav_exact[1] * gpi->a_grav_exact[1] +
-                                 gpi->a_grav_exact[2] * gpi->a_grav_exact[2]);
-
-      /* Compute relative error */
-      const float err_rel = diff_norm / a_norm;
-
-      /* Check that we are not above tolerance */
-      if (err_rel > rel_tol) {
-        message(
-            "Error too large !"
-            "gp->a_grav=[%3.6e %3.6e %3.6e] "
-            "gp->a_exact=[%3.6e %3.6e %3.6e], "
-            "gp->num_interacted=%lld, err=%f",
-            gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2],
-            gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2], gpi->num_interacted,
-            err_rel);
-
-        continue;
-      }
-
-      /* Construct some statistics */
-      err_rel_max = max(err_rel_max, fabsf(err_rel));
-      err_rel_min = min(err_rel_min, fabsf(err_rel));
-      err_rel_mean += err_rel;
-      err_rel_mean2 += err_rel * err_rel;
-      counter++;
     }
   }
 
+  message("Written SWIFT accelerations in file '%s'.", file_name_swift);
+
   /* Be nice */
-  fclose(file);
+  fclose(file_swift);
 
-  /* Final operation on the stats */
-  if (counter > 0) {
-    err_rel_mean /= counter;
-    err_rel_mean2 /= counter;
-    err_rel_std = sqrtf(err_rel_mean2 - err_rel_mean * err_rel_mean);
-  }
+  if (!gravity_exact_force_file_exits(e)) {
+
+    char file_name_exact[100];
+    sprintf(file_name_exact, "gravity_checks_exact_step%d.dat", e->step);
+
+    FILE *file_exact = fopen(file_name_exact, "w");
+    fprintf(file_exact, "# Gravity accuracy test - EXACT FORCES\n");
+    fprintf(file_exact, "# G= %16.8e\n", e->physical_constants->const_newton_G);
+    fprintf(file_exact, "# N= %d\n", SWIFT_GRAVITY_FORCE_CHECKS);
+    fprintf(file_exact, "# epsilon=%16.8e\n", e->gravity_properties->epsilon);
+    fprintf(file_exact, "# theta=%16.8e\n", e->gravity_properties->theta_crit);
+    fprintf(file_exact, "# %16s %16s %16s %16s %16s %16s %16s\n", "id",
+            "pos[0]", "pos[1]", "pos[2]", "a_exact[0]", "a_exact[1]",
+            "a_exact[2]");
 
-  /* Report on the findings */
-  message("Checked gravity for %d gparts.", counter);
-  message("Error on |a_grav|: min=%e max=%e mean=%e std=%e", err_rel_min,
-          err_rel_max, err_rel_mean, err_rel_std);
+    /* Output particle exact accelerations  */
+    for (size_t i = 0; i < s->nr_gparts; ++i) {
 
+      struct gpart *gpi = &s->gparts[i];
+
+      /* Is the particle was active and part of the subset to be tested ? */
+      if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
+          gpart_is_starting(gpi, e)) {
+
+        fprintf(
+            file_exact, "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e \n",
+            gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
+            gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2]);
+      }
+    }
+
+    message("Written exact accelerations in file '%s'.", file_name_exact);
+
+    /* Be nice */
+    fclose(file_exact);
+  }
 #else
   error("Gravity checking function called without the corresponding flag.");
 #endif