several bug fixes.

Former-commit-id: d0161ae876e742cd8a08b6af01a461704a3cfb59

src/cell.h

@@ -38,7 +38,7 @@ struct cell {
     double h[3];
     
     /* Max radii in this cell. */
-    double r_max;
+    double h_max;
     
     /* The depth of this cell in the tree. */
     int depth, split;

src/engine.c

@@ -55,6 +55,10 @@
  * @brief Sort the tasks in topological order over all queues.
  *
  * @param e The #engine.
+ *
+ * TODO: Return the indices tid as these are the tasks sorted according
+ * to their ranks. They can then be dropped into the queues in order
+ * of these indices.
  */
  
 void engine_ranktasks ( struct engine *e ) {

src/runner.c

@@ -485,6 +485,9 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
                 printf( "runner_doghost: particle %i has bad wcount=%f.\n" , p->id , p->wcount + kernel_root );
                 pid[redo] = pid[i];
                 redo += 1;
+                p->wcount = 0.0;
+                p->rho = 0.0;
+                p->rho_dh = 0.0;
                 continue;
                 }
 

src/runner_doiact.h

@@ -99,7 +99,7 @@ void DOPAIR_NAIVE ( struct runner *r , struct cell *ci , struct cell *cj ) {
             pj = &parts_j[ pjd ];
         
             /* Compute the pairwise distance. */
-            r2 = 0.0;
+            r2 = 0.0f;
             for ( k = 0 ; k < 3 ; k++ ) {
                 dx[k] = pix[k] - pj->x[k];
                 r2 += dx[k]*dx[k];
@@ -117,7 +117,7 @@ void DOPAIR_NAIVE ( struct runner *r , struct cell *ci , struct cell *cj ) {
         } /* loop over the parts in ci. */
         
     #ifdef TIMER_VERBOSE
-        printf( "runner_dopair_naive[%02i]: %i/%i parts at depth %i (r_max=%.3f/%.3f) took %.3f ms.\n" , r->id , count_i , count_j , ci->depth , ci->r_max , cj->r_max , ((double)TIMER_TOC(TIMER_DOPAIR)) / CPU_TPS * 1000 );
+        printf( "runner_dopair_naive[%02i]: %i/%i parts at depth %i (r_max=%.3f/%.3f) took %.3f ms.\n" , r->id , count_i , count_j , ci->depth , ci->h_max , cj->h_max , ((double)TIMER_TOC(TIMER_DOPAIR)) / CPU_TPS * 1000 );
     #else
         TIMER_TOC(TIMER_DOPAIR);
     #endif
@@ -144,7 +144,8 @@ void DOPAIR ( struct runner *r , struct cell *ci , struct cell *cj ) {
     struct part *pi, *pj, *parts_i, *parts_j;
     double pix[3], pjx[3], di, dj;
     float dx[3], hi, hi2, hj, hj2, r2;
-    double hi_max, hj_max, di_max, dj_min;
+    double hi_max, hj_max;
+    double di_max, dj_min;
     int count_i, count_j;
     TIMER_TIC
     
@@ -172,9 +173,17 @@ void DOPAIR ( struct runner *r , struct cell *ci , struct cell *cj ) {
     for ( rshift = 0.0 , k = 0 ; k < 3 ; k++ )
         rshift += shift[k]*runner_shift[ 3*sid + k ];
         
-    /* printf( "runner_dopair: doing pair [ %g %g %g ]/[ %g %g %g ] with %i/%i parts and shift = [ %g %g %g ].\n" ,
-        ci->loc[0] , ci->loc[1] , ci->loc[2] , cj->loc[0] , cj->loc[1] , cj->loc[2] ,
-        ci->count , cj->count , shift[0] , shift[1] , shift[2] ); fflush(stdout); */
+    /* for ( k = 0 ; k < ci->count ; k++ )
+        if ( ci->parts[k].id == 561590 )
+            break;
+    if ( k == ci->count )
+        for ( k = 0 ; k < cj->count ; k++ )
+            if ( cj->parts[k].id == 561590 )
+                break;
+    if ( k < cj->count )
+        printf( "runner_dopair: doing pair [ %g %g %g ]/[ %g %g %g ] with %i/%i parts, h_max=%g/%g, and shift = [ %g %g %g ] (rshift=%g).\n" ,
+            ci->loc[0] , ci->loc[1] , ci->loc[2] , cj->loc[0] , cj->loc[1] , cj->loc[2] ,
+            ci->count , cj->count , ci->h_max , cj->h_max , shift[0] , shift[1] , shift[2] , rshift ); fflush(stdout); */
     /* for ( hi = 0 , k = 0 ; k < ci->count ; k++ )
         hi += ci->parts[k].r;
     for ( hj = 0 , k = 0 ; k < cj->count ; k++ )
@@ -186,7 +195,7 @@ void DOPAIR ( struct runner *r , struct cell *ci , struct cell *cj ) {
     sort_j = &cj->sort[ sid*(cj->count + 1) ];
     
     /* Get some other useful values. */
-    hi_max = ci->r_max - rshift; hj_max = cj->r_max - rshift;
+    hi_max = ci->h_max - rshift; hj_max = cj->h_max;
     count_i = ci->count; count_j = cj->count;
     parts_i = ci->parts; parts_j = cj->parts;
     di_max = sort_i[count_i-1].d - rshift;
@@ -216,7 +225,7 @@ void DOPAIR ( struct runner *r , struct cell *ci , struct cell *cj ) {
             pj = &parts_j[ sort_j[pjd].i ];
         
             /* Compute the pairwise distance. */
-            r2 = 0.0;
+            r2 = 0.0f;
             for ( k = 0 ; k < 3 ; k++ ) {
                 dx[k] = pix[k] - pj->x[k];
                 r2 += dx[k]*dx[k];
@@ -237,7 +246,7 @@ void DOPAIR ( struct runner *r , struct cell *ci , struct cell *cj ) {
     tic = getticks(); */
 
     /* Loop over the parts in cj. */
-    for ( pjd = 0 ; pjd < count_j && sort_j[pjd].d - hj_max < di_max ; pjd++ ) {
+    for ( pjd = 0 ; pjd < count_j && ( 1 || sort_j[pjd].d - hj_max < di_max ) ; pjd++ ) {
     
         /* Get a hold of the jth part in cj. */
         pj = &parts_j[ sort_j[ pjd ].i ];
@@ -257,7 +266,7 @@ void DOPAIR ( struct runner *r , struct cell *ci , struct cell *cj ) {
             pi = &parts_i[ sort_i[pid].i ];
             
             /* Compute the pairwise distance. */
-            r2 = 0.0;
+            r2 = 0.0f;
             for ( k = 0 ; k < 3 ; k++ ) {
                 dx[k] = pi->x[k] - pjx[k];
                 r2 += dx[k]*dx[k];
@@ -275,7 +284,7 @@ void DOPAIR ( struct runner *r , struct cell *ci , struct cell *cj ) {
         } /* loop over the parts in ci. */
 
     #ifdef TIMER_VERBOSE
-        printf( "runner_dopair[%02i]: %i/%i parts at depth %i (r_max=%.3f/%.3f, h=%.3f) took %.3f ms.\n" , r->id , count_i , count_j , ci->depth , ci->r_max , cj->r_max , fmax(ci->h[0],fmax(ci->h[1],ci->h[2])) , ((double)(TIMER_TOC(TIMER_DOPAIR))) / CPU_TPS * 1000 );
+        printf( "runner_dopair[%02i]: %i/%i parts at depth %i (r_max=%.3f/%.3f, h=%.3f) took %.3f ms.\n" , r->id , count_i , count_j , ci->depth , ci->h_max , cj->h_max , fmax(ci->h[0],fmax(ci->h[1],ci->h[2])) , ((double)(TIMER_TOC(TIMER_DOPAIR))) / CPU_TPS * 1000 );
     #else
         TIMER_TOC(TIMER_DOPAIR);
     #endif
@@ -320,7 +329,7 @@ void DOSELF ( struct runner *r , struct cell *c ) {
             pj = &parts[pjd];
         
             /* Compute the pairwise distance. */
-            r2 = 0.0;
+            r2 = 0.0f;
             for ( k = 0 ; k < 3 ; k++ ) {
                 dx[k] = pix[k] - pj->x[k];
                 r2 += dx[k]*dx[k];
@@ -330,7 +339,7 @@ void DOSELF ( struct runner *r , struct cell *c ) {
             if ( r2 < hi2 || r2 < pj->h*pj->h ) {
             
                 IACT( r2 , dx , hi , pj->h , pi , pj );
-            
+                
                 }
         
             } /* loop over all other particles. */

src/space.c

@@ -41,9 +41,6 @@
 /* Error macro. */
 #define error(s) { printf( "%s:%s:%i: %s\n" , __FILE__ , __FUNCTION__ , __LINE__ , s ); abort(); }
 
-/* Convert cell location to ID. */
-#define cell_getid( cdim , i , j , k ) ( (int)(k) + (cdim)[2]*( (int)(j) + (cdim)[1]*(int)(i) ) )
-
 /* Split size. */
 int space_splitsize = space_splitsize_default;
 
@@ -95,6 +92,10 @@ void space_map_clearsort ( struct cell *c , void *data ) {
 
 /**
  * @brief Mapping function to append a ghost task to each cell.
+ *
+ * Looks for the super cell, e.g. the highest-level cell above each
+ * cell for which a pair is defined. All ghosts below this cell will
+ * depend on the ghost of their parents (sounds spooky, but it isn't).
  */
 
 void space_map_mkghosts ( struct cell *c , void *data ) {
@@ -109,7 +110,7 @@ void space_map_mkghosts ( struct cell *c , void *data ) {
         if ( finger->nr_tasks > 0 )
             c->super = finger;
             
-    /* Make the ghost task, if needed. */
+    /* Make the ghost task */
     if ( c->super != c || c->nr_tasks > 0 )
         c->ghost = space_addtask( s , task_type_ghost , task_subtype_none , 0 , 0 , c , NULL , NULL , 0 , NULL , 0 );
 
@@ -359,7 +360,7 @@ void space_splittasks ( struct space *s ) {
         hj = fmax( cj->h[0] , fmax( cj->h[1] , cj->h[2] ) );
             
         /* Should this task be split-up? */
-        if ( ci->split && cj->split && ci->r_max < hi/2 && cj->r_max < hj/2 ) {
+        if ( ci->split && cj->split && ci->h_max < hi/2 && cj->h_max < hj/2 ) {
         
             /* Get the relative distance between the pairs, wrapping. */
             for ( k = 0 ; k < 3 ; k++ ) {
@@ -800,7 +801,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
     int i, j, k, ii, jj, kk, iii, jjj, kkk, cid, cjd;
     int *cdim = s->cdim;
     int nr_tasks_old = s->nr_tasks;
-    struct task *t, *t2;
+    struct task *t , *t2;
     int pts[7][8] = { { -1 , 12 , 10 , 9 , 4 , 3 , 1 , 0 } ,
                       { -1 , -1 , 11 , 10 , 5 , 4 , 2 , 1 } ,
                       { -1 , -1 , -1 , 12 , 7 , 6 , 4 , 3 } , 
@@ -968,7 +969,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
                     if ( t->flags & ( 1 << i ) ) {
                         for ( j = 0 ; j < 8 ; j++ )
                             if ( t->ci->progeny[j] != NULL )
-                                task_rmunlock( t->ci->progeny[j]->sorts[i] , t );
+                                task_rmunlock_blind( t->ci->progeny[j]->sorts[i] , t );
                         t->ci->sorts[i] = NULL;
                         }
             t->type = task_type_none;
@@ -976,7 +977,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
         }
         
     /* Count the number of tasks associated with each cell. */
-    space_map_cells( s , 1 , &space_map_mkghosts , NULL );
+    space_map_cells( s , 1 , &space_map_clearnrtasks , NULL );
     for ( k = 0 ; k < s->nr_tasks ; k++ ) {
         t = &s->tasks[k];
         if ( t->type == task_type_self )
@@ -1010,8 +1011,6 @@ void space_maketasks ( struct space *s , int do_sort ) {
             
         /* Otherwise, pair interaction? */
         else if ( t->type == task_type_pair && t->subtype == task_subtype_density ) {
-            if ( t->ci->ghost == NULL )
-                printf( "space_maketasks: cell at %lx has no ghost!\n" , (unsigned long int)t->ci );
             task_addunlock( t , t->ci->super->ghost );
             task_addunlock( t , t->cj->super->ghost );
             t2 = space_addtask( s , task_type_pair , task_subtype_force , 0 , 0 , t->ci , t->cj , NULL , 0 , NULL , 0 );
@@ -1069,7 +1068,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
 void space_split ( struct space *s , struct cell *c ) {
 
     int k, count;
-    double h, h_limit, r_max = 0.0;
+    double h, h_limit, h_max = 0.0;
     struct cell *temp;
     
     /* Check the depth. */
@@ -1084,10 +1083,10 @@ void space_split ( struct space *s , struct cell *c ) {
         h = c->parts[k].h;
         if ( h <= h_limit )
             count += 1;
-        if ( h > r_max )
-            r_max = h;
+        if ( h > h_max )
+            h_max = h;
         }
-    c->r_max = r_max;
+    c->h_max = h_max;
             
     /* Split or let it be? */
     if ( count > c->count*space_splitratio && c->count > space_splitsize ) {
@@ -1113,7 +1112,7 @@ void space_split ( struct space *s , struct cell *c ) {
                 temp->loc[2] += temp->h[2];
             temp->depth = c->depth + 1;
             temp->split = 0;
-            temp->r_max = 0.0;
+            temp->h_max = 0.0;
             temp->parent = c;
             c->progeny[k] = temp;
             }
@@ -1237,7 +1236,7 @@ void space_init ( struct space *s , double dim[3] , struct part *parts , int N ,
     if ( h_cells < h_max )
         h_cells = h_max;
     for ( k = 0 ; k < 3 ; k++ ) {
-        cdim[k] = ceil( dim[k] / h_cells );
+        cdim[k] = floor( dim[k] / h_cells );
         h[k] = dim[k] / cdim[k];
         ih[k] = 1.0 / h[k];
         }

src/space.h

@@ -24,10 +24,13 @@
 #define space_maxdepth                  10
 #define space_cellallocchunk            1000
 #define space_splitratio                0.875
-#define space_splitsize_default         800
+#define space_splitsize_default         400
 #define space_dosub                     1
 
 
+/* Convert cell location to ID. */
+#define cell_getid( cdim , i , j , k ) ( (int)(k) + (cdim)[2]*( (int)(j) + (cdim)[1]*(int)(i) ) )
+
 /* Split size. */
 extern int space_splitsize;
 

src/task.c

@@ -66,6 +66,30 @@ void task_rmunlock( struct task *ta , struct task *tb ) {
     }
     
 
+/**
+ * @brief Remove an unlock_task from the given task.
+ *
+ * @param ta The unlocking #task.
+ * @param tb The #task that will be unlocked.
+ *
+ * Differs from #task_rmunlock in that it will not fail if
+ * the task @c tb is not in the unlocks of @c ta.
+ */
+ 
+void task_rmunlock_blind( struct task *ta , struct task *tb ) {
+
+    int k;
+    
+    for ( k = 0 ; k < ta->nr_unlock_tasks ; k++ )
+        if ( ta->unlock_tasks[k] == tb ) {
+            ta->nr_unlock_tasks -= 1;
+            ta->unlock_tasks[k] = ta->unlock_tasks[ ta->nr_unlock_tasks ];
+            return;
+            }
+
+    }
+    
+
 /**
  * @brief Add an unlock_task to the given task.
  *

src/task.h

@@ -62,4 +62,5 @@ struct task {
 
 /* Function prototypes. */
 void task_rmunlock( struct task *ta , struct task *tb );
+void task_rmunlock_blind( struct task *ta , struct task *tb );
 void task_addunlock( struct task *ta , struct task *tb );