Work in progress

src/engine.c

@@ -2110,11 +2110,11 @@ void engine_init_particles(struct engine *e) {
   if (e->policy & engine_policy_self_gravity) {
 
     mask |= 1 << task_type_grav_up;
-    mask |= 1 << task_type_self;
-    mask |= 1 << task_type_pair;
-    mask |= 1 << task_type_sub;
+    /* mask |= 1 << task_type_self; */
+    /* mask |= 1 << task_type_pair; */
+    /* mask |= 1 << task_type_sub; */
 
-    submask |= 1 << task_subtype_grav;
+    /* submask |= 1 << task_subtype_grav; */
   }
 
   /* Add the tasks corresponding to external gravity to the masks */
@@ -2335,20 +2335,20 @@ void engine_step(struct engine *e) {
   fclose(file);
 
   /* Check the gravity accelerations */
-  /* struct gpart *temp = malloc(s->nr_gparts * sizeof(struct gpart)); */
-  /* memcpy(temp, s->gparts, s->nr_gparts * sizeof(struct gpart)); */
-  /* gravity_n2(temp, s->nr_gparts); */
-  /* file = fopen("grav_brute.dat", "w"); */
-  /* for(size_t k = 0; k < s->nr_gparts; ++k) { */
-  /*   fprintf(file, "%lld %f %f %f %f %f %f\n", temp[k].id, */
-  /* 	    temp[k].x[0], temp[k].x[1], temp[k].x[2], */
-  /* 	    temp[k].a_grav[0], temp[k].a_grav[1], temp[k].a_grav[2]); */
-  /* } */
-  /* fclose(file); */
+  struct gpart *temp = malloc(s->nr_gparts * sizeof(struct gpart));
+  memcpy(temp, s->gparts, s->nr_gparts * sizeof(struct gpart));
+  gravity_n2(temp, s->nr_gparts);
+  file = fopen("grav_brute.dat", "w");
+  for(size_t k = 0; k < s->nr_gparts; ++k) {
+    fprintf(file, "%lld %f %f %f %f %f %f\n", temp[k].id,
+  	    temp[k].x[0], temp[k].x[1], temp[k].x[2],
+  	    temp[k].a_grav[0], temp[k].a_grav[1], temp[k].a_grav[2]);
+  }
+  fclose(file);
 
-  /* free(temp); */
+  free(temp);
 
-  float mass = 0.f;
+  double mass = 0.f;
   for (int i = 0; i < s->nr_cells; ++i) mass += s->cells[i].multipole.mass;
   message("Total mass: %f", mass);
 

src/gravity/Default/gravity_part.h

@@ -47,7 +47,7 @@ struct gpart {
   /* float tx; */
   /* float tv; */
 
-  float mass_interacted;
+  double mass_interacted;
   
   /* Anonymous union for id/part. */
   union {

src/runner_doiact_grav.h

@@ -25,7 +25,7 @@
 #include "clocks.h"
 #include "part.h"
 
-#define ICHECK 4934
+#define ICHECK 1
 
 /**
  * @brief Compute the recursive upward sweep, i.e. construct the
@@ -125,8 +125,9 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
     if (gp->ti_end > ti_current) continue;
 
     if (gp->id == -ICHECK)
-      message("id=%lld mass=%f CoM=[%f %f %f]\n", gp->id, multi.mass,
-              multi.CoM[0], multi.CoM[1], multi.CoM[2]);
+      message("id=%lld x=[%f %f %f] mass=%f CoM=[ %f %f %f ] size= %f\n", gp->id,
+	      gp->x[0], gp->x[1], gp->x[2], multi.mass,
+              multi.CoM[0], multi.CoM[1], multi.CoM[2], cj->h[0]);
 
     /* Compute the pairwise distance. */
     const float dx[3] = {multi.CoM[0] - gp->x[0],   // x
@@ -210,6 +211,23 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
   /* Anything to do here? */
   // if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
 
+  for (int pid = 0; pid < gcount_i; pid++) {
+
+    struct gpart *restrict gp = &gparts_i[pid];
+  if (gp->id == -ICHECK)
+    message("id=%lld size=%f\n", gp->id, cj->h[0]);
+
+  }
+
+  for (int pid = 0; pid < gcount_j; pid++) {
+
+    struct gpart *restrict gp = &gparts_j[pid];
+  if (gp->id == -ICHECK)
+    message("id=%lld size=%f\n", gp->id, ci->h[0]);
+
+  }
+
+  
   /* Loop over all particles in ci... */
   for (int pid = 0; pid < gcount_i; pid++) {
 
@@ -371,23 +389,23 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
     /*     cj->loc[2], cj->h[0]); */
   }
 
-  for (int i = 0; i < gcount_i; i++) {
+  /* for (int i = 0; i < gcount_i; i++) { */
 
-    struct gpart *const gp = &ci->gparts[i];
+  /*   struct gpart *const gp = &ci->gparts[i]; */
 
-    if (gp->id == -ICHECK)
-      message("id=%lld a=[%f %f %f]\n", gp->id, gp->a_grav[0], gp->a_grav[1],
-              gp->a_grav[2]);
-  }
+  /*   if (gp->id == -ICHECK) */
+  /*     message("id=%lld a=[%f %f %f]\n", gp->id, gp->a_grav[0], gp->a_grav[1], */
+  /*             gp->a_grav[2]); */
+  /* } */
 
-  for (int i = 0; i < gcount_j; i++) {
+  /* for (int i = 0; i < gcount_j; i++) { */
 
-    struct gpart *const gp = &cj->gparts[i];
+  /*   struct gpart *const gp = &cj->gparts[i]; */
 
-    if (gp->id == -ICHECK)
-      message("id=%lld a=[%f %f %f]\n", gp->id, gp->a_grav[0], gp->a_grav[1],
-              gp->a_grav[2]);
-  }
+  /*   if (gp->id == -ICHECK) */
+  /*     message("id=%lld a=[%f %f %f]\n", gp->id, gp->a_grav[0], gp->a_grav[1], */
+  /*             gp->a_grav[2]); */
+  /* } */
 
   /* Are both cells split ? */
   if (ci->split && cj->split) {

src/scheduler.c

@@ -172,8 +172,8 @@ void scheduler_splittasks(struct scheduler *s) {
       if (ci->split) {
 
         /* Make a sub? */
-        if (scheduler_dosub && (ci->count * ci->count < space_subsize ||
-                                ci->gcount * ci->gcount < space_subsize)) {
+        if (scheduler_dosub && (ci->count * ci->count < space_subsize)) {
+	  //                                ci->gcount * ci->gcount < space_subsize)) {
 
           /* convert to a self-subtask. */
           t->type = task_type_sub;
@@ -238,8 +238,8 @@ void scheduler_splittasks(struct scheduler *s) {
 
         /* Replace by a single sub-task? */
         if (scheduler_dosub &&
-            (ci->count * cj->count * sid_scale[sid] < space_subsize ||
-             ci->gcount * cj->gcount * sid_scale[sid] < space_subsize) &&
+            (ci->count * cj->count * sid_scale[sid] < space_subsize) &&// ||
+	    //ci->gcount * cj->gcount * sid_scale[sid] < space_subsize) &&
             sid != 0 && sid != 2 && sid != 6 && sid != 8) {
 
           /* Make this task a sub task. */