ghost task now also corrects smoothing length of particles that have strayed.

Former-commit-id: af8e0aa1402ab080c36fa3892b915febd5e83b9e

src/cell.h

@@ -67,6 +67,10 @@ struct cell {
     /* The tasks computing this cell's sorts. */
     struct task *sorts[14];
     
+    /* The tasks computing this cell's density. */
+    struct task *density[27];
+    int nr_density;
+    
     /* The ghost task to link density to interactions. */
     struct task *ghost;
     

src/engine.c

@@ -169,6 +169,14 @@ void engine_run ( struct engine *e , int sort_queues ) {
     int j, k;
     struct space *s = e->s;
     
+    /* Re-set the particle data. */
+    for ( k = 0 ; k < s->nr_parts ; k++ ) {
+        s->parts[k].wcount = 0.0f;
+        s->parts[k].wcount_dh = 0.0f;
+        s->parts[k].rho = 0.0f;
+        s->parts[k].rho_dh = 0.0f;
+        }
+    
     /* Run throught the tasks and get all the waits right. */
     for ( k = 0 ; k < s->nr_tasks ; k++ ) {
         s->tasks[k].done = 0;

src/part.h

@@ -72,6 +72,7 @@ struct part {
     /* Particle number density. */
     int icount;
     float wcount;
+    float wcount_dh;
     
     } __attribute__((aligned (32)));
     

src/runner.c

@@ -442,6 +442,7 @@ void runner_dosort ( struct runner *r , struct cell *c , int flags ) {
 void runner_doghost ( struct runner *r , struct cell *c ) {
 
     struct part *p;
+    struct cell *finger;
     int i, k, redo, count = c->count;
     int *pid;
     float ihg, ihg2;
@@ -478,14 +479,18 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
             ihg2 = ihg * ihg;
             p->rho *= ihg * ihg2;
             p->rho_dh *= ihg2 * ihg2;
+            
+            /* Update the smoothing length. */
+            p->h -= ( p->wcount + kernel_root - const_nwneigh ) / p->wcount_dh;
 
             /* Did we get the right number density? */
             if ( p->wcount + kernel_root > const_nwneigh + 1 ||
                  p->wcount + kernel_root < const_nwneigh - 1 ) {
-                printf( "runner_doghost: particle %i has bad wcount=%f.\n" , p->id , p->wcount + kernel_root );
+                printf( "runner_doghost: particle %i has bad wcount=%f.\n" , p->id , p->wcount + kernel_root ); fflush(stdout);
                 pid[redo] = pid[i];
                 redo += 1;
                 p->wcount = 0.0;
+                p->wcount_dh = 0.0;
                 p->rho = 0.0;
                 p->rho_dh = 0.0;
                 continue;
@@ -512,8 +517,40 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
             
         /* Re-set the counter for the next loop (potentially). */
         count = redo;
-        if ( count > 0 )
-            error( "Bad smoothing length, fixing this isn't implemented yet." );
+        if ( count > 0 ) {
+        
+            // error( "Bad smoothing length, fixing this isn't implemented yet." );
+            
+            /* Climb up the cell hierarchy. */
+            for ( finger = c ; finger != NULL ; finger = finger->parent ) {
+            
+                /* Run through this cell's density interactions. */
+                for ( k = 0 ; k < finger->nr_density ; k++ ) {
+                
+                    /* Self-interaction? */
+                    if ( finger->density[k]->type == task_type_self )
+                        runner_doself_subset_density( r , finger , c->parts , pid , count );
+                        
+                    /* Otherwise, pair interaction? */
+                    else if ( finger->density[k]->type == task_type_pair ) {
+                    
+                        /* Left or right? */
+                        if ( finger->density[k]->ci == finger )
+                            runner_dopair_subset_density( r , finger , c->parts , pid , count , finger->density[k]->cj );
+                        else
+                            runner_dopair_subset_density( r , finger , c->parts , pid , count , finger->density[k]->ci );
+                        
+                        }
+                
+                    /* Otherwise, sub interaction? */
+                    else if ( finger->density[k]->type == task_type_sub ) 
+                        runner_dosub_subset_density( r , finger->density[k]->ci , finger->density[k]->cj , c , c->parts , pid , count , finger->density[k]->flags );
+                
+                    }
+            
+                }
+        
+            }
             
         }
 

src/runner.h

@@ -29,6 +29,7 @@ enum {
     runner_timer_dopair_force,
     runner_timer_dosub_density,
     runner_timer_dosub_force,
+    runner_timer_dopair_subset,
     runner_timer_doghost,
     runner_timer_getpair,
     runner_timer_steal,

src/runner_iact.h

@@ -69,31 +69,35 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_density ( float r
 
     float r = sqrtf( r2 );
     float xi, xj;
-    float hg_inv;
+    float h_inv, hg_inv;
     float wi, wj, wi_dx, wj_dx;
     
     if ( r2 < hi*hi && pi != NULL ) {
         
-        hg_inv = kernel_igamma / hi;
+        h_inv = 1.0 / hi;
+        hg_inv = kernel_igamma * h_inv;
         xi = r * hg_inv;
         kernel_deval( xi , &wi , &wi_dx );
         
         pi->rho += pj->mass * wi;
         pi->rho_dh += -pj->mass * ( 3.0*wi + xi*wi_dx );
         pi->wcount += wi * ( 4.0f * M_PI / 3.0f * kernel_igamma3 );
+        pi->wcount_dh -= xi * h_inv * wi_dx * ( 4.0f * M_PI / 3.0f * kernel_igamma3 );
         pi->icount += 1;
         
         }
 
     if ( r2 < hj*hj && pj != NULL ) {
         
-        hg_inv = kernel_igamma / hj;
+        h_inv = 1.0 / hj;
+        hg_inv = kernel_igamma * h_inv;
         xj = r * hg_inv;
         kernel_deval( xj , &wj , &wj_dx );
         
         pj->rho += pi->mass * wj;
         pj->rho_dh += -pi->mass * ( 3.0*wj + xj*wj_dx );
         pj->wcount += wj * ( 4.0f * M_PI / 3.0f * kernel_igamma3 );
+        pj->wcount_dh -= xj * h_inv * wj_dx * ( 4.0f * M_PI / 3.0f * kernel_igamma3 );
         pj->icount += 1;
         
         }

src/space.c

@@ -129,6 +129,7 @@ void space_map_mkghosts ( struct cell *c , void *data ) {
 void space_map_clearnrtasks ( struct cell *c , void *data ) {
 
     c->nr_tasks = 0;
+    c->nr_density = 0;
 
     }
 
@@ -976,23 +977,43 @@ void space_maketasks ( struct space *s , int do_sort ) {
             }
         }
         
-    /* Count the number of tasks associated with each cell. */
+    /* Count the number of tasks associated with each cell and
+       store the density tasks in each cell. */
     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 )
+        if ( t->type == task_type_self ) {
             t->ci->nr_tasks += 1;
+            if ( t->subtype == task_subtype_density ) {
+                t->ci->density[ t->ci->nr_density ] = t;
+                t->ci->nr_density += 1;
+                }
+            }
         else if ( t->type == task_type_pair ) {
             t->ci->nr_tasks += 1;
             t->cj->nr_tasks += 1;
+            if ( t->subtype == task_subtype_density ) {
+                t->ci->density[ t->ci->nr_density ] = t;
+                t->ci->nr_density += 1;
+                t->cj->density[ t->cj->nr_density ] = t;
+                t->cj->nr_density += 1;
+                }
             }
         else if ( t->type == task_type_sub ) {
             t->ci->nr_tasks += 1;
             if ( t->cj != NULL )
                 t->cj->nr_tasks += 1;
+            if ( t->subtype == task_subtype_density ) {
+                t->ci->density[ t->ci->nr_density ] = t;
+                t->ci->nr_density += 1;
+                if ( t->cj != NULL ) {
+                    t->cj->density[ t->cj->nr_density ] = t;
+                    t->cj->nr_density += 1;
+                    }
+                }
             }
         }
-    
+        
     /* Append a ghost task to each cell. */
     space_map_cells( s , 1 , &space_map_mkghosts , s );