better parallelization and some tweaks to space_prepare, space_rebuild, space_maketasks.

Former-commit-id: 410fd574761013d435619b0c1469f3a7e182a9b5

src/atomic.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * 
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ * 
+ ******************************************************************************/
+
+
+
+/* Get the inlining right. */
+#ifndef INLINE
+# if __GNUC__ && !__GNUC_STDC_INLINE__
+#  define INLINE extern inline
+# else
+#  define INLINE inline
+# endif
+#endif
+    
+#define atomic_add(v,i) __sync_fetch_and_add( v , i )
+#define atomic_inc(v) atomic_add( v , 1 )

src/engine.c

@@ -238,7 +238,9 @@ void engine_map_kick_first ( struct cell *c , void *data ) {
             /* Integrate other values if this particle will not be updated. */
             /* Init fields for density calculation. */
             if ( pdt > dt_step ) {
-                rho = p->rho *= expf( -3.0f * h_dt / h * dt );
+                // rho = p->rho *= expf( -3.0f * h_dt / h * dt );
+                float w = -3.0f * h_dt / h * dt;
+                rho = p->rho *= 1.0f + w*( -1.0f + w*( 0.5f - 1.0f/6.0f*w ) );
                 p->force.POrho2 = u * ( const_gamma - 1.0f ) / ( rho + h * p->rho_dh / 3.0f );
                 }
             else {
@@ -285,11 +287,6 @@ void engine_map_kick_first ( struct cell *c , void *data ) {
     c->dx_max = dx_max;
     c->h2dx_max = h2dx_max;
     
-    /* Clean out the task pointers. */
-    // c->sorts[0] = NULL;
-    // c->nr_tasks = 0;
-    // c->nr_density = 0;
-    
     }
 
 
@@ -519,17 +516,25 @@ void engine_step ( struct engine *e , int sort_queues ) {
     e->step += 1;
     
     /* Does the time step need adjusting? */
-    if ( e->dt == 0 )
+    if ( e->dt == 0 ) {
         e->dt = e->dt_orig;
-    while ( dt_min < e->dt ) {
-        e->dt *= 0.5;
-        e->step *= 2;
-        printf( "engine_step: dt_min dropped below time step, adjusting to dt=%e.\n" , e->dt );
+        while ( dt_min < e->dt )
+            e->dt *= 0.5;
+        while ( dt_min > 2*e->dt )
+            e->dt *= 2.0;
+        printf( "engine_step: dt_min=%.3e, adjusting time step to dt=%e.\n" , dt_min , e->dt );
         }
-    while ( dt_min > 2*e->dt && (e->step & 1) == 0 ) {
-        e->dt *= 2.0;
-        e->step /= 2;
-        printf( "engine_step: dt_min is larger than twice the time step, adjusting to dt=%e.\n" , e->dt );
+    else {
+        while ( dt_min < e->dt ) {
+            e->dt *= 0.5;
+            e->step *= 2;
+            printf( "engine_step: dt_min dropped below time step, adjusting to dt=%e.\n" , e->dt );
+            }
+        while ( dt_min > 2*e->dt && (e->step & 1) == 0 ) {
+            e->dt *= 2.0;
+            e->step /= 2;
+            printf( "engine_step: dt_min is larger than twice the time step, adjusting to dt=%e.\n" , e->dt );
+            }
         }
     
     /* Set the system time. */

src/runner_doiact.h

@@ -1815,8 +1815,18 @@ void DOSUB1 ( struct runner *r , struct cell *ci , struct cell *cj , int sid ) {
             }
             
         /* Otherwise, compute the pair directly. */
-        else
+        else {
+        
+            /* Do any of the cells need to be sorted first? */
+            if ( !(ci->sorted & (1 << sid) ) )
+                runner_dosort( r , ci , (1 << sid) , 0 );
+            if ( !(cj->sorted & (1 << sid) ) )
+                runner_dosort( r , cj , (1 << sid) , 0 );
+        
+            /* Compute the interactions. */
             DOPAIR1( r , ci , cj );
+            
+            }
     
         } /* otherwise, pair interaction. */
     
@@ -2085,8 +2095,18 @@ void DOSUB2 ( struct runner *r , struct cell *ci , struct cell *cj , int sid ) {
             }
             
         /* Otherwise, compute the pair directly. */
-        else
+        else {
+        
+            /* Do any of the cells need to be sorted first? */
+            if ( !(ci->sorted & (1 << sid) ) )
+                runner_dosort( r , ci , (1 << sid) , 0 );
+            if ( !(cj->sorted & (1 << sid) ) )
+                runner_dosort( r , cj , (1 << sid) , 0 );
+        
+            /* Compute the interactions. */
             DOPAIR2( r , ci , cj );
+            
+            }
     
         } /* otherwise, pair interaction. */
     

src/space.h

@@ -59,7 +59,7 @@ struct space {
     double h[3], ih[3];
     
     /* The minimum and maximum cutoff radii. */
-    double h_min, h_max, cell_min;
+    double h_max, cell_min;
     
     /* Current time step for particles. */
     float dt_step;
@@ -92,7 +92,8 @@ struct space {
     
     /* The list of tasks. */
     struct task *tasks;
-    int nr_tasks, tasks_size;
+    int nr_tasks, tasks_next;
+    int tasks_size;
     int *tasks_ind;
     
     /* General-purpose lock for this space. */

src/swift.h

@@ -24,6 +24,7 @@
 #include "cycle.h"
 #include "timers.h"
 #include "const.h"
+#include "atomic.h"
 #include "lock.h"
 #include "task.h"
 #include "part.h"

src/task.c

@@ -75,7 +75,7 @@ void task_cleanunlock ( struct task *t , int type ) {
  * @param tb The #task that will be unlocked.
  */
  
-void task_rmunlock( struct task *ta , struct task *tb ) {
+void task_rmunlock ( struct task *ta , struct task *tb ) {
 
     int k;
     
@@ -103,7 +103,7 @@ void task_rmunlock( struct task *ta , struct task *tb ) {
  * the task @c tb is not in the unlocks of @c ta.
  */
  
-void task_rmunlock_blind( struct task *ta , struct task *tb ) {
+void task_rmunlock_blind ( struct task *ta , struct task *tb ) {
 
     int k;
     
@@ -128,7 +128,19 @@ void task_rmunlock_blind( struct task *ta , struct task *tb ) {
  * @param tb The #task that will be unlocked.
  */
  
-void task_addunlock( struct task *ta , struct task *tb ) {
+void task_addunlock ( struct task *ta , struct task *tb ) {
+
+    /* Add the lock atomically. */
+    ta->unlock_tasks[ __sync_fetch_and_add( &ta->nr_unlock_tasks, 1 ) ] = tb;
+
+    /* Check a posteriori if we did not overshoot. */
+    if ( ta->nr_unlock_tasks > task_maxunlock )
+        error( "Too many unlock_tasks in task." );
+        
+    }
+    
+
+void task_addunlock_old ( struct task *ta , struct task *tb ) {
 
     int k;
     
@@ -137,6 +149,7 @@ void task_addunlock( struct task *ta , struct task *tb ) {
     /* Check if ta already unlocks tb. */
     for ( k = 0 ; k < ta->nr_unlock_tasks ; k++ )
         if ( ta->unlock_tasks[k] == tb ) {
+            error( "Duplicate unlock." );
             lock_unlock_blind( &ta->lock );
             return;
             }