fix density computation, include self.

Former-commit-id: e2337231bd3a93920a628055854316bdf7d240d0

src/Makefile.am

@@ -20,11 +20,11 @@
 AUTOMAKE_OPTIONS=gnu
 
 # Add the debug flag to the whole thing
-AM_CFLAGS = -g -O3 -Wall -Werror -ffast-math -fstrict-aliasing -ftree-vectorize \
-    -funroll-loops $(SIMD_FLAGS) $(OPENMP_CFLAGS) \
-    -DTIMER -DCOUNTER -DCPU_TPS=2.67e9
-# AM_CFLAGS = -Wall -Werror $(OPENMP_CFLAGS) \
+# AM_CFLAGS = -g -O3 -Wall -Werror -ffast-math -fstrict-aliasing -ftree-vectorize \
+#     -funroll-loops $(SIMD_FLAGS) $(OPENMP_CFLAGS) \
 #     -DTIMER -DCOUNTER -DCPU_TPS=2.67e9
+AM_CFLAGS = -Wall -Werror $(OPENMP_CFLAGS) \
+    -DTIMER -DCOUNTER -DCPU_TPS=2.67e9
 
 # Assign a "safe" version number
 AM_LDFLAGS = $(LAPACK_LIBS) $(BLAS_LIBS) $(HDF5_LDFLAGS) -version-info 0:0:0

src/engine.c

@@ -208,16 +208,16 @@ void engine_step ( struct engine *e , int sort_queues ) {
         u_bar[k] = p->u + hdt * p->u_dt;
         
         /* Move the particles with the velocitie at the half-step. */
-        p->x[0] += dt * v_bar[3*k+0];
-        p->x[1] += dt * v_bar[3*k+1];
-        p->x[2] += dt * v_bar[3*k+2];
+        // p->x[0] += dt * v_bar[3*k+0];
+        // p->x[1] += dt * v_bar[3*k+1];
+        // p->x[2] += dt * v_bar[3*k+2];
         
         /* Update positions and energies at the half-step. */
-        p->v[0] += dt * p->a[0];
-        p->v[1] += dt * p->a[1];
-        p->v[2] += dt * p->a[2];
-        p->u *= expf( p->u_dt / p->u * dt );
-        p->h *= expf( -1.0f * p->h_dt / p->h * dt );
+        // p->v[0] += dt * p->a[0];
+        // p->v[1] += dt * p->a[1];
+        // p->v[2] += dt * p->a[2];
+        // p->u *= expf( p->u_dt / p->u * dt );
+        // p->h *= expf( -1.0f * p->h_dt / p->h * dt );
         
         /* Integrate other values if this particle will not be updated. */
         if ( p->dt > dt_max ) {

src/ic.c

@@ -333,7 +333,6 @@ void read_ic ( char* fileName, double dim[3], struct part **parts,  int* N, int*
   readArray(h_grp, "ParticleIDs", ULONGLONG, *N, 1, *parts, id, COMPULSORY);
   readArray(h_grp, "TimeStep", FLOAT, *N, 1, *parts, dt, OPTIONAL);
   readArray(h_grp, "Acceleration", FLOAT, *N, 3, *parts, a, OPTIONAL);
-  readArray(h_grp, "TimeStep", FLOAT, *N, 1, *parts, dt, OPTIONAL );
   readArray(h_grp, "Density", FLOAT, *N, 1, *parts, rho, OPTIONAL );
   readArray(h_grp, "Acceleration", FLOAT, *N, 3, *parts, a, OPTIONAL );
 

src/runner.c

@@ -360,16 +360,17 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
                 /* Adjust the computed rho. */
                 ihg = kernel_igamma / p->h;
                 ihg2 = ihg * ihg;
-                p->rho *= ihg * ihg2;
+                p->rho = ihg * ihg2 * ( p->rho + p->mass*kernel_root );
                 p->rho_dh *= ihg2 * ihg2;
+                p->wcount += kernel_wroot;
 
                 /* Update the smoothing length. */
-                p->h -= ( p->wcount + kernel_root - const_nwneigh ) / p->wcount_dh;
+                p->h -= ( p->wcount - 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 %lli (h=%e,depth=%i) has bad wcount=%f.\n" , p->id , p->h , c->depth , p->wcount + kernel_root ); fflush(stdout);
+                if ( p->wcount > const_nwneigh + 1 ||
+                     p->wcount < const_nwneigh - 1 ) {
+                    printf( "runner_doghost: particle %lli (h=%e,depth=%i) has bad wcount=%f.\n" , p->id , p->h , c->depth , p->wcount ); fflush(stdout);
                     // p->h += ( p->wcount + kernel_root - const_nwneigh ) / p->wcount_dh;
                     pid[redo] = pid[i];
                     redo += 1;

src/runner_iact.h

@@ -68,7 +68,8 @@ static float kernel_coeffs[ (kernel_degree + 1) * (kernel_ivals + 1) ] __attribu
     { 3.0/4.0*M_1_PI , -3.0/2.0*M_1_PI , 0.0 , M_1_PI , 
       -0.25*M_1_PI , 3.0/2.0*M_1_PI , -3.0*M_1_PI , M_2_PI , 
       0.0 , 0.0 , 0.0 , 0.0 };
-#define kernel_root ( 4.0/3.0*M_PI*kernel_igamma3*kernel_coeffs[ kernel_degree ] )
+#define kernel_root ( kernel_coeffs[ kernel_degree ] )
+#define kernel_wroot ( 4.0/3.0*M_PI*kernel_igamma3*kernel_coeffs[ kernel_degree ] )
       
       
 /**
@@ -148,7 +149,7 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_density ( float r
         kernel_deval( xi , &wi , &wi_dx );
         
         pi->rho += pj->mass * wi;
-        pi->rho_dh += -pj->mass * ( 3.0*wi + xi*wi_dx );
+        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;
@@ -353,8 +354,8 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_force ( float r2
     /* Use the force, Luke! */
     for ( k = 0 ; k < 3 ; k++ ) {
         f = dx[k] * w;
-        pi->a[k] += pj->mass * f;
-        pj->a[k] += -pi->mass * f;
+        pi->a[k] -= pj->mass * f;
+        pj->a[k] += pi->mass * f;
         }
         
     /* Compute dv dot r. */
@@ -486,8 +487,8 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_vec_force ( float
         pi[k]->h_dt += pih_dt.f[k];
         pj[k]->h_dt += pjh_dt.f[k];
         for ( j = 0 ; j < 3 ; j++ ) {
-            pi[k]->a[j] += pia[j].f[k];
-            pj[k]->a[j] -= pja[j].f[k];
+            pi[k]->a[j] -= pia[j].f[k];
+            pj[k]->a[j] += pja[j].f[k];
             }
         }
         
@@ -532,7 +533,7 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_nonsym_force ( fl
     /* Use the force, Luke! */
     for ( k = 0 ; k < 3 ; k++ ) {
         f = dx[k] * w;
-        pi->a[k] += pj->mass * f;
+        pi->a[k] -= pj->mass * f;
         }
         
     /* Compute dv dot r. */
@@ -646,7 +647,7 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_nonsym_vec_force
         pi[k]->u_dt += piu_dt.f[k];
         pi[k]->h_dt += pih_dt.f[k];
         for ( j = 0 ; j < 3 ; j++ )
-            pi[k]->a[j] += pia[j].f[k];
+            pi[k]->a[j] -= pia[j].f[k];
         }
 
 #else