test_125 now checks all particles against the similar brute-force call

src/hydro/Gadget2/hydro_part.h

@@ -61,46 +61,46 @@ struct part {
   /* Derivative of the density with respect to smoothing length. */
   float rho_dh;
 
-  // union {
+  union {
 
-  struct {
+    struct {
 
-    /* Number of neighbours. */
-    float wcount;
+      /* Number of neighbours. */
+      float wcount;
 
-    /* Number of neighbours spatial derivative. */
-    float wcount_dh;
+      /* Number of neighbours spatial derivative. */
+      float wcount_dh;
 
-    /* Particle velocity curl. */
-    float rot_v[3];
+      /* Particle velocity curl. */
+      float rot_v[3];
 
-    /* Particle velocity divergence. */
-    float div_v;
+      /* Particle velocity divergence. */
+      float div_v;
 
-  } density;
+    } density;
 
-  struct {
+    struct {
 
-    /* Balsara switch */
-    float balsara;
+      /* Balsara switch */
+      float balsara;
 
-    /* Pressure over density squared (including drho/dh term) */
-    float P_over_rho2;
+      /* Pressure over density squared (including drho/dh term) */
+      float P_over_rho2;
 
-    /* Particle sound speed. */
-    float soundspeed;
+      /* Particle sound speed. */
+      float soundspeed;
 
-    /* Signal velocity. */
-    float v_sig;
+      /* Signal velocity. */
+      float v_sig;
 
-    /* Entropy time derivative */
-    float entropy_dt;
+      /* Entropy time derivative */
+      float entropy_dt;
 
-    /* Time derivative of the smoothing length */
-    float h_dt;
+      /* Time derivative of the smoothing length */
+      float h_dt;
 
-  } force;
-  //  };
+    } force;
+  };
 
   /* Particle ID. */
   long long id;

src/tools.c

@@ -240,6 +240,70 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
   }
 }
 
+void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
+
+  float r2, hi, hj, hig2, hjg2, dx[3];
+  struct part *pi, *pj;
+
+  /* Implements a double-for loop and checks every interaction */
+  for (int i = 0; i < ci->count; ++i) {
+
+    pi = &ci->parts[i];
+    hi = pi->h;
+    hig2 = hi * hi * kernel_gamma2;
+
+    for (int j = 0; j < cj->count; ++j) {
+
+      pj = &cj->parts[j];
+      hj = pj->h;
+      hjg2 = hj * hj * kernel_gamma2;
+
+      /* Pairwise distance */
+      r2 = 0.0f;
+      for (int k = 0; k < 3; k++) {
+        dx[k] = ci->parts[i].x[k] - cj->parts[j].x[k];
+        r2 += dx[k] * dx[k];
+      }
+
+      /* Hit or miss? */
+      if (r2 < hig2 || r2 < hjg2) {
+
+        /* Interact */
+        runner_iact_nonsym_force(r2, dx, hi, hj, pi, pj);
+      }
+    }
+  }
+
+  /* Reverse double-for loop and checks every interaction */
+  for (int j = 0; j < cj->count; ++j) {
+
+    pj = &cj->parts[j];
+    hj = pj->h;
+    hjg2 = hj * hj * kernel_gamma2;
+
+    for (int i = 0; i < ci->count; ++i) {
+
+      pi = &ci->parts[i];
+      hi = pi->h;
+      hig2 = hi * hi * kernel_gamma2;
+
+      /* Pairwise distance */
+      r2 = 0.0f;
+      for (int k = 0; k < 3; k++) {
+        dx[k] = cj->parts[j].x[k] - ci->parts[i].x[k];
+        r2 += dx[k] * dx[k];
+      }
+
+      /* Hit or miss? */
+      if (r2 < hjg2 || r2 < hig2) {
+
+        /* Interact */
+        runner_iact_nonsym_force(r2, dx, hj, pi->h, pj, pi);
+      }
+    }
+  }
+}
+
 void self_all_density(struct runner *r, struct cell *ci) {
   float r2, hi, hj, hig2, hjg2, dxi[3];  //, dxj[3];
   struct part *pi, *pj;
@@ -287,6 +351,42 @@ void self_all_density(struct runner *r, struct cell *ci) {
   }
 }
 
+void self_all_force(struct runner *r, struct cell *ci) {
+  float r2, hi, hj, hig2, hjg2, dxi[3];  //, dxj[3];
+  struct part *pi, *pj;
+
+  /* Implements a double-for loop and checks every interaction */
+  for (int i = 0; i < ci->count; ++i) {
+
+    pi = &ci->parts[i];
+    hi = pi->h;
+    hig2 = hi * hi * kernel_gamma2;
+
+    for (int j = i + 1; j < ci->count; ++j) {
+
+      pj = &ci->parts[j];
+      hj = pj->h;
+      hjg2 = hj * hj * kernel_gamma2;
+
+      if (pi == pj) continue;
+
+      /* Pairwise distance */
+      r2 = 0.0f;
+      for (int k = 0; k < 3; k++) {
+        dxi[k] = ci->parts[i].x[k] - ci->parts[j].x[k];
+        r2 += dxi[k] * dxi[k];
+      }
+
+      /* Hit or miss? */
+      if (r2 < hig2 || r2 < hjg2) {
+
+        /* Interact */
+        runner_iact_force(r2, dxi, hi, hj, pi, pj);
+      }
+    }
+  }
+}
+
 void pairs_single_grav(double *dim, long long int pid,
                        struct gpart *restrict gparts, const struct part *parts,
                        int N, int periodic) {

src/tools.h

@@ -40,6 +40,8 @@ void pairs_single_density(double *dim, long long int pid,
 
 void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj);
 void self_all_density(struct runner *r, struct cell *ci);
+void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj);
+void self_all_force(struct runner *r, struct cell *ci);
 
 void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic);
 

tests/difffloat.py

@@ -89,6 +89,7 @@ for i in range(n_lines_to_check):
         abs_diff = abs(data1[i,j] - data2[i,j])
 
         sum = abs(data1[i,j] + data2[i,j])
+        if sum < 1e8: continue
         if sum > 0:
             rel_diff = abs(data1[i,j] - data2[i,j]) / sum
         else:

tests/test125cells.sh.in

 #!/bin/bash
-rm brute_force_125_standard.dat swift_dopair_125_standard.dat
-
-./test125cells -p 4 -r 1 -d 0 -f standard
-
-#python @srcdir@/difffloat.py brute_force_125_standard.dat swift_dopair_125_standard.dat @srcdir@/tolerance.dat 4
-
+for v in {0..3}
+do
+    for p in {0..2}
+    do
+	rm brute_force_125_standard.dat swift_dopair_125_standard.dat
+	./test125cells -n 6 -r 1 -v $v -p $p -f standard
+	python @srcdir@/difffloat.py brute_force_125_standard.dat swift_dopair_125_standard.dat @srcdir@/tolerance_125.dat 6
+	echo ""
+    done
+done
+	
 exit $?

tests/test27cells.c

@@ -249,7 +249,7 @@ void runner_doself1_density(struct runner *r, struct cell *ci);
 /* And go... */
 int main(int argc, char *argv[]) {
   size_t runs = 0, particles = 0;
-  double h = 1.2348, size = 1., rho = 1.;
+  double h = 1.23485, size = 1., rho = 1.;
   double perturbation = 0.;
   char outputFileNameExtension[200] = "";
   char outputFileName[200] = "";
@@ -266,7 +266,7 @@ int main(int argc, char *argv[]) {
   srand(0);
 
   char c;
-  while ((c = getopt(argc, argv, "m:s:h:p:r:t:d:f:v:")) != -1) {
+  while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:")) != -1) {
     switch (c) {
       case 'h':
         sscanf(optarg, "%lf", &h);
@@ -274,7 +274,7 @@ int main(int argc, char *argv[]) {
       case 's':
         sscanf(optarg, "%lf", &size);
         break;
-      case 'p':
+      case 'n':
         sscanf(optarg, "%zu", &particles);
         break;
       case 'r':
@@ -300,9 +300,10 @@ int main(int argc, char *argv[]) {
 
   if (h < 0 || particles == 0 || runs == 0) {
     printf(
-        "\nUsage: %s -p PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
-        "\nGenerates a cell pair, filled with particles on a Cartesian grid."
-        "\nThese are then interacted using runner_dopair1_density."
+        "\nUsage: %s -n PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
+        "\nGenerates 27 cells, filled with particles on a Cartesian grid."
+        "\nThese are then interacted using runner_dopair1_density() and "
+        "runner_doself1_density()."
         "\n\nOptions:"
         "\n-h DISTANCE=1.2348 - Smoothing length in units of <x>"
         "\n-m rho             - Physical density in the cell"
@@ -316,13 +317,14 @@ int main(int argc, char *argv[]) {
   }
 
   /* Help users... */
+  message("Adiabatic index: ga = %f", hydro_gamma);
   message("Smoothing length: h = %f", h * size);
   message("Kernel:               %s", kernel_name);
-  message("Neighbour target: N = %f",
-          h * h * h * 4.0 * M_PI * kernel_gamma3 / 3.0);
+  message("Neighbour target: N = %f", h * h * h * kernel_norm);
   message("Density target: rho = %f", rho);
   message("div_v target:   div = %f", vel == 2 ? 3.f : 0.f);
   message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f);
+
   printf("\n");
 
   /* Build the infrastructure */

tests/test27cells.sh.in

 #!/bin/bash
-rm brute_force_27_standard.dat swift_dopair_27_standard.dat
 
-./test27cells -p 6 -r 1 -d 0 -f standard
-
-python @srcdir@/difffloat.py brute_force_27_standard.dat swift_dopair_27_standard.dat @srcdir@/tolerance.dat 6
+for v in {0..3}
+do
+    rm brute_force_27_standard.dat swift_dopair_27_standard.dat
+    ./test27cells -n 6 -r 1 -d 0 -f standard -v $v
+    python @srcdir@/difffloat.py brute_force_27_standard.dat swift_dopair_27_standard.dat @srcdir@/tolerance_27.dat 6
+    echo ""
+done
 
 exit $?

tests/test27cellsPerturbed.sh.in

 #!/bin/bash
 rm brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat
 
-./test27cells -p 6 -r 1 -d 0.1 -f perturbed
+./test27cells -n 6 -r 1 -d 0.1 -f perturbed
 
-python @srcdir@/difffloat.py brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat @srcdir@/tolerance.dat 6
+python @srcdir@/difffloat.py brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat @srcdir@/tolerance_27.dat 6
 
 exit $?

tests/testPair.sh.in

@@ -3,6 +3,6 @@ rm brute_force_standard.dat swift_dopair_standard.dat
 
 ./testPair -p 6 -r 1 -d 0 -f standard
 
-python @srcdir@/difffloat.py brute_force_standard.dat swift_dopair_standard.dat @srcdir@/tolerance.dat
+python @srcdir@/difffloat.py brute_force_standard.dat swift_dopair_standard.dat @srcdir@/tolerance_pair.dat
 
 exit $?

tests/testPairPerturbed.sh.in

@@ -3,6 +3,6 @@ rm brute_force_perturbed.dat swift_dopair_perturbed.dat
 
 ./testPair -p 6 -r 1 -d 0.1 -f perturbed
 
-python @srcdir@/difffloat.py brute_force_perturbed.dat swift_dopair_perturbed.dat @srcdir@/tolerance.dat
+python @srcdir@/difffloat.py brute_force_perturbed.dat swift_dopair_perturbed.dat @srcdir@/tolerance_pair.dat
 
 exit $?

tests/tolerance_125.dat

+#   ID    pos_x    pos_y    pos_z      v_x      v_y      v_z        h      rho    div_v        S        u        P        c      a_x      a_y      a_z     h_dt    v_sig    dS/dt    du/dt
+    0	  1e-5	   1e-5	    1e-5       1e-5	1e-5	 1e-5	    1e-5   1e-5	  1e-5	       1e-5	1e-5	 1e-5	  1e-5	 1e-5	  1e-5	   1e-5	   1e-5	   1e-5	    1e-5     1e-5
+    0	  1e-5	   1e-5	    1e-5       1e-5	1e-5	 1e-5	    1e-5   1e-5	  1e-5	       1e-5	1e-5	 1e-5	  1e-5	 1e-5	  1e-5	   1e-5	   1e-5	   1e-5	    1e-5     1e-5

tests/tolerance_pair.dat

+#   ID      pos_x      pos_y      pos_z        v_x        v_y        v_z           rho        rho_dh        wcount     wcount_dh         div_v       curl_vx       curl_vy       curl_vz
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-5	      1e-5	    2e-5       3e-2		 1e-5	     1e-5	   1e-5		 1e-5
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-5	      1.2e-5	    1e-5       1e-2		 1e-4	     1e-4	   1e-4		 1e-4