changed definition of h, now consistent with literature.

Former-commit-id: 13a176e745a0a5e7fcbec354a0d90c4cbef96d48

examples/SedovBlast/makeIC.py

@@ -67,7 +67,7 @@ for i in range(L):
             v[index,1] = 0.
             v[index,2] = 0.
             m[index] = mass
-            h[index] = 2.251 * boxSize / L
+            h[index] = 2.251 / 2 * boxSize / L
             u[index] = internalEnergy
             ids[index] = index
             if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 2.01 * boxSize/L:

examples/SedovBlast/rdf.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ #                    Matthieu Schaller (matthieu.schaller@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/>.
+ # 
+ ##############################################################################
+
+import h5py
+import random
+import sys
+import math
+from numpy import *
+
+# Reads the HDF5 output of SWIFT and generates a radial density profile
+# of the different physical values.
+
+# Input values?
+if len(sys.argv) < 3 :
+    print "Useage: " , sys.argv[0] , " <filename> <nr. bins>"
+    exit()
+    
+# Get the input arguments
+fileName = sys.argv[1];
+nr_bins = int( sys.argv[2] );
+
+
+# Open the file
+fileName = sys.argv[1];
+file = h5py.File( fileName , 'r' )
+
+# Get the space dimensions.
+grp = file[ "/Header" ]
+boxsize = grp.attrs[ 'BoxSize' ]
+boxsize = boxsize[0]
+
+# Get the particle data
+grp = file.get( '/PartType0' )
+ds = grp.get( 'Coordinates' )
+coords = ds[...]
+ds = grp.get( 'Velocities' )
+v = ds[...]
+ds = grp.get( 'Masses' )
+m = ds[...]
+ds = grp.get( 'SmoothingLength' )
+h = ds[...]
+ds = grp.get( 'InternalEnergy' )
+u = ds[...]
+ds = grp.get( 'ParticleIDs' )
+ids = ds[...]
+ds = grp.get( 'Density' )
+rho = ds[...]
+
+# Set the origin to be the middle of the box
+origin = [ boxsize / 2 , boxsize / 2 , boxsize / 2 ]
+
+# Get the maximum radius
+r_max = boxsize / 2
+
+# Init the bins
+nr_parts = coords.shape[0]
+bins_v = zeros( nr_bins )
+bins_m = zeros( nr_bins )
+bins_h = zeros( nr_bins )
+bins_u = zeros( nr_bins )
+bins_rho = zeros( nr_bins )
+bins_count = zeros( nr_bins )
+bins_P = zeros( nr_bins )
+
+# Loop over the particles and fill the bins.
+for i in range( nr_parts ):
+
+    # Get the box index.
+    r = coords[i] - origin
+    r = sqrt( r[0]*r[0] + r[1]*r[1] + r[2]*r[2] )
+    ind = floor( r / r_max * nr_bins )
+    
+    # Update the bins
+    if ( ind < nr_bins ):
+        bins_count[ind] += 1
+        bins_v[ind] += sqrt( v[i,0]*v[i,0] + v[i,1]*v[i,1] + v[i,2]*v[i,2] )
+        bins_m[ind] += m[i]
+        bins_h[ind] += h[i]
+        bins_u[ind] += u[i]
+        bins_rho[ind] += rho[i]
+        bins_P[ind] += (2.0/3)*u[i]*rho[i]
+    
+# Loop over the bins and dump them
+print "# bucket left right count v m h u rho"
+for i in range( nr_bins ):
+
+    # Normalize by the bin volume.
+    r_left = r_max * i / nr_bins
+    r_right = r_max * (i+1) / nr_bins
+    vol = 4/3*math.pi*(r_right*r_right*r_right - r_left*r_left*r_left)
+    ivol = 1.0 / vol
+
+    print "%i %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e" % \
+        ( i , r_left , r_right , \
+          bins_count[i] * ivol , \
+          bins_v[i] / bins_count[i] , \
+          bins_m[i] * ivol , \
+          bins_h[i] / bins_count[i] , \
+          bins_u[i] / bins_count[i] , \
+          bins_rho[i] / bins_count[i] ,
+          bins_P[i] / bins_count[i] )
+    
+    

examples/SedovBlast/solution.py

@@ -33,7 +33,7 @@ P_0 = 1.e-5         # Background Pressure
 E_0 = 1.e2          # Energy of the explosion
 gamma = 5./3.     # Gas polytropic index
 
-t = 0.15           # Time of the solution
+t = 0.2           # Time of the solution
 
 N = 1000          # Number of radial points
 R_max = 3.        # Maximal radius

examples/SodShock/rhox.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ #                    Matthieu Schaller (matthieu.schaller@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/>.
+ # 
+ ##############################################################################
+
+import h5py
+import random
+import sys
+import math
+from numpy import *
+
+# Reads the HDF5 output of SWIFT and generates a radial density profile
+# of the different physical values.
+
+# Input values?
+if len(sys.argv) < 3 :
+    print "Useage: " , sys.argv[0] , " <filename> <nr. bins>"
+    exit()
+    
+# Get the input arguments
+fileName = sys.argv[1];
+nr_bins = int( sys.argv[2] );
+
+
+# Open the file
+fileName = sys.argv[1];
+file = h5py.File( fileName , 'r' )
+
+# Get the space dimensions.
+grp = file[ "/Header" ]
+boxsize = grp.attrs[ 'BoxSize' ]
+boxsize = boxsize[0]
+
+# Get the particle data
+grp = file.get( '/PartType0' )
+ds = grp.get( 'Coordinates' )
+coords = ds[...]
+ds = grp.get( 'Velocities' )
+v = ds[...]
+# ds = grp.get( 'Mass' )
+# m = ds[...]
+ds = grp.get( 'SmoothingLength' )
+h = ds[...]
+ds = grp.get( 'InternalEnergy' )
+u = ds[...]
+ds = grp.get( 'ParticleIDs' )
+ids = ds[...]
+ds = grp.get( 'Density' )
+rho = ds[...]
+
+# Get the maximum radius
+r_max = boxsize
+
+# Init the bins
+nr_parts = coords.shape[0]
+bins_v = zeros( nr_bins )
+bins_m = zeros( nr_bins )
+bins_h = zeros( nr_bins )
+bins_u = zeros( nr_bins )
+bins_rho = zeros( nr_bins )
+bins_count = zeros( nr_bins )
+bins_P = zeros( nr_bins )
+
+# Loop over the particles and fill the bins.
+for i in range( nr_parts ):
+
+    # Get the box index.
+    r = coords[i,0]
+    ind = floor( r / r_max * nr_bins )
+    
+    # Update the bins
+    bins_count[ind] += 1
+    bins_v[ind] += v[i,0] # sqrt( v[i,0]*v[i,0] + v[i,1]*v[i,1] + v[i,2]*v[i,2] )
+    # bins_m[ind] += m[i]
+    bins_h[ind] += h[i]
+    bins_u[ind] += u[i]
+    bins_rho[ind] += rho[i]
+    bins_P[ind] += (2.0/3)*u[i]*rho[i]
+    
+# Loop over the bins and dump them
+print "# bucket left right count v m h u rho"
+for i in range( nr_bins ):
+
+    # Normalize by the bin volume.
+    r_left = r_max * i / nr_bins
+    r_right = r_max * (i+1) / nr_bins
+    vol = 4/3*math.pi*(r_right*r_right*r_right - r_left*r_left*r_left)
+    ivol = 1.0 / vol
+
+    print "%i %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e" % \
+        ( i , r_left , r_right , \
+          bins_count[i] * ivol , \
+          bins_v[i] / bins_count[i] , \
+          bins_m[i] * ivol , \
+          bins_h[i] / bins_count[i] , \
+          bins_u[i] / bins_count[i] , \
+          bins_rho[i] / bins_count[i] ,
+          bins_P[i] / bins_count[i] )
+    
+    

examples/test.c

@@ -646,7 +646,7 @@ void kernel_dump ( int N ) {
     // float dw_dx4[4] __attribute__ ((aligned (16)));
 
     for ( k = 0 ; k <= N ; k++ ) {
-        x = ((float)k) / N * kernel_igamma;
+        x = ((float)k) / N;
         x4[3] = x4[2]; x4[2] = x4[1]; x4[1] = x4[0]; x4[0] = x;
         kernel_deval( x , &w , &dw_dx );
         // kernel_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
@@ -790,13 +790,13 @@ int main ( int argc , char *argv[] ) {
 
     /* Read particles and space information from (GADGET) IC */
     tic = getticks();
-    read_ic(ICfileName, dim, &parts, &N, &periodic);
+    read_ic( ICfileName , dim , &parts , &N , &periodic );
     printf( "main: reading particle properties took %.3f ms.\n" , ((double)(getticks() - tic)) / CPU_TPS * 1000 ); fflush(stdout);
     
     /* Apply h scaling */
     if(scaling != 1.0)
       for ( k = 0 ; k < N ; k++ )
-	parts[k].h *= scaling;
+	    parts[k].h *= scaling;
     
     /* Apply shift */
     if(shift[0] !=0 || shift[1] !=0 || shift[2] !=0 )

examples/test_single.c

@@ -64,15 +64,15 @@ int main ( int argc , char *argv[] ) {
 
     int k, N = 100;
     struct part p1, p2;
-    float r2, dx[3] = { 0.0f , 0.0f , 0.0f };
+    float x, w, dwdx, r2, dx[3] = { 0.0f , 0.0f , 0.0f };
     
     /* Init the particles. */
     for ( k = 0 ; k < 3 ; k++ ) {
         p1.a[k] = 0.0f; p1.v[k] = 0.0f; p1.x[k] = 0.0;
         p2.a[k] = 0.0f; p2.v[k] = 0.0f; p2.x[k] = 0.0;
         }
-    p1.rho = 1.0f; p1.mass = 9.7059e-4; p1.h = 0.222871287;
-    p2.rho = 1.0f; p2.mass = 9.7059e-4; p2.h = 0.222871287;
+    p1.rho = 1.0f; p1.mass = 9.7059e-4; p1.h = 0.222871287 / 2;
+    p2.rho = 1.0f; p2.mass = 9.7059e-4; p2.h = 0.222871287 / 2;
     p1.force.c = 0.0f; p1.force.balsara = 0.0f;
     p2.force.c = 0.0f; p2.force.balsara = 0.0f;
     p1.u = 1.e-5 / ((const_gamma - 1.)*p1.rho);
@@ -87,19 +87,32 @@ int main ( int argc , char *argv[] ) {
     for ( k = 1 ; k <= N ; k++ ) {
     
         /* Set the distance/radius. */
-        dx[0] = -((float)k)/N * fmaxf( p1.h , p2.h );
+        dx[0] = -((float)k)/N * fmaxf( p1.h , p2.h ) * kernel_gamma;
         r2 = dx[0]*dx[0];
         
         /* Clear the particle fields. */
-        p1.a[0] = 0.0f; p1.force.u_dt = 0.0f;
-        p2.a[0] = 0.0f; p2.force.u_dt = 0.0f;
+        /* p1.a[0] = 0.0f; p1.force.u_dt = 0.0f;
+        p2.a[0] = 0.0f; p2.force.u_dt = 0.0f; */
         
         /* Interact the particles. */
-        runner_iact_force( r2 , dx , p1.h , p2.h , &p1 , &p2 );
+        // runner_iact_force( r2 , dx , p1.h , p2.h , &p1 , &p2 );
+        
+        /* Clear the particle fields. */
+        p1.rho = 0.0f; p1.density.wcount = 0.0f;
+        p2.rho = 0.0f; p2.density.wcount = 0.0f;
+        
+        /* Interact the particles. */
+        runner_iact_density( r2 , dx , p1.h , p2.h , &p1 , &p2 );
+        
+        /* Evaluate just the kernel. */
+        x = fabsf( dx[0] ) / p1.h;
+        kernel_deval( x , &w , &dwdx );
         
         /* Output the results. */
-        printf( "%.3e %.3e %.3e %.3e %.3e\n" ,
-            -dx[0] , p1.a[0] , p1.force.u_dt , p2.a[0] , p2.force.u_dt );
+        printf( "%.3e %.3e %.3e %.3e %.3e %.3e %.3e\n" ,
+            // -dx[0] , p1.a[0] , p1.force.u_dt , p2.a[0] , p2.force.u_dt ,
+            -dx[0] , p1.rho , p1.density.wcount , p2.rho , p2.density.wcount ,
+            w , dwdx );
     
         } /* loop over radii. */
     

src/const.h

@@ -21,7 +21,7 @@
 
 /* Hydrodynamical constants. */
 #define const_gamma             (5.0f/3.0f)
-#define const_viscosity_alpha   0.8
+#define const_viscosity_alpha   0.8f
 
 /* Time integration constants. */
 #define const_cfl               0.3f

src/engine.c

@@ -229,7 +229,7 @@ void engine_map_kick_first ( struct cell *c , void *data ) {
                         (x[1] - x_old[1])*(x[1] - x_old[1]) +
                         (x[2] - x_old[2])*(x[2] - x_old[2]) );
             dx_max = fmaxf( dx_max , dx );
-            h2dx_max = fmaxf( h2dx_max , dx*2 + h );
+            h2dx_max = fmaxf( h2dx_max , dx*2 + h*kernel_gamma );
 
             /* Update positions and energies at the half-step. */
             p->v[0] = v[0] += dt * a[0];
@@ -362,12 +362,63 @@ void engine_collect_kick2 ( struct cell *c ) {
  * @brief Compute the force on a single particle brute-force.
  */
 
-void engine_single ( double *dim , long long int pid , struct part *__restrict__ parts , int N , int periodic ) {
+void engine_single_density ( double *dim , long long int pid , struct part *__restrict__ parts , int N , int periodic ) {
 
     int i, k;
     double r2, dx[3];
-    float fdx[3], ihg;
-    struct part p, p2;
+    float fdx[3], ih;
+    struct part p;
+    
+    /* Find "our" part. */
+    for ( k = 0 ; k < N && parts[k].id != pid ; k++ );
+    if ( k == N )
+        error( "Part not found." );
+    p = parts[k];
+    
+    /* Clear accumulators. */
+    ih = 1.0f / p.h;
+    p.rho = 0.0f; p.rho_dh = 0.0f;
+    p.density.wcount = 0.0f; p.density.wcount_dh = 0.0f;
+	p.density.div_v = 0.0;
+	for ( k=0 ; k < 3 ; k++)
+		p.density.curl_v[k] = 0.0;
+            
+    /* Loop over all particle pairs (force). */
+    for ( k = 0 ; k < N ; k++ ) {
+        if ( parts[k].id == p.id )
+            continue;
+        for ( i = 0 ; i < 3 ; i++ ) {
+            dx[i] = p.x[i] - parts[k].x[i];
+            if ( periodic ) {
+                if ( dx[i] < -dim[i]/2 )
+                    dx[i] += dim[i];
+                else if ( dx[i] > dim[i]/2 )
+                    dx[i] -= dim[i];
+                }
+            fdx[i] = dx[i];
+            }
+        r2 = fdx[0]*fdx[0] + fdx[1]*fdx[1] + fdx[2]*fdx[2];
+        if ( r2 < p.h*p.h*kernel_gamma2 ) {
+            runner_iact_nonsym_density( r2 , fdx , p.h , parts[k].h , &p , &parts[k] );
+            }
+        }
+        
+    /* Dump the result. */
+    p.rho = ih * ih * ih * ( p.rho + p.mass*kernel_root );
+    p.rho_dh = p.rho_dh * ih * ih * ih * ih;
+    p.density.wcount = ( p.density.wcount + kernel_root ) * ( 4.0f / 3.0 * M_PI * kernel_gamma3 );
+    printf( "pairs_single: part %lli (h=%e) has wcount=%e, rho=%e, rho_dh=%e.\n" , p.id , p.h , p.density.wcount , p.rho , p.rho_dh );
+    fflush(stdout);
+    
+    }
+
+
+void engine_single_force ( double *dim , long long int pid , struct part *__restrict__ parts , int N , int periodic ) {
+
+    int i, k;
+    double r2, dx[3];
+    float fdx[3];
+    struct part p;
     
     /* Find "our" part. */
     for ( k = 0 ; k < N && parts[k].id != pid ; k++ );
@@ -376,12 +427,12 @@ void engine_single ( double *dim , long long int pid , struct part *__restrict__
     p = parts[k];
     
     /* Clear accumulators. */
-    ihg = kernel_igamma / p.h;
     p.a[0] = 0.0f; p.a[1] = 0.0f; p.a[2] = 0.0f;
     p.force.u_dt = 0.0f; p.force.h_dt = 0.0f; p.force.v_sig = 0.0f;
             
     /* Loop over all particle pairs (force). */
-    for ( k = 0 ; k < N ; k++ ) {
+    // for ( k = 0 ; k < N ; k++ ) {
+    for ( k = N-1 ; k >= 0 ; k-- ) {
         if ( parts[k].id == p.id )
             continue;
         for ( i = 0 ; i < 3 ; i++ ) {
@@ -395,21 +446,19 @@ void engine_single ( double *dim , long long int pid , struct part *__restrict__
             fdx[i] = dx[i];
             }
         r2 = fdx[0]*fdx[0] + fdx[1]*fdx[1] + fdx[2]*fdx[2];
-        if ( r2 < p.h*p.h || r2 < parts[k].h*parts[k].h ) {
+        if ( r2 < p.h*p.h*kernel_gamma2 || r2 < parts[k].h*parts[k].h*kernel_gamma2 ) {
             p.a[0] = 0.0f; p.a[1] = 0.0f; p.a[2] = 0.0f;
             p.force.u_dt = 0.0f; p.force.h_dt = 0.0f; p.force.v_sig = 0.0f;
-            p2 = parts[k];
-            runner_iact_nonsym_force( r2 , fdx , p.h , p2.h , &p , &p2 );
-            printf( "pairs_simple: interacting particles %lli and %lli (rho=%.3e,r=%e): a=[%.3e,%.3e,%.3e], u_dt=%.3e, h_dt=%.3e, v_sig=%.3e.\n" ,
+            runner_iact_nonsym_force( r2 , fdx , p.h , parts[k].h , &p , &parts[k] );
+            /* printf( "pairs_simple: interacting particles %lli and %lli (rho=%.3e,r=%e): a=[%.3e,%.3e,%.3e], u_dt=%.3e, h_dt=%.3e, v_sig=%.3e.\n" ,
                 pid , p2.id , p2.rho , sqrtf(r2) ,
                 p.a[0] , p.a[1] , p.a[2] ,
-                p.force.u_dt , p.force.h_dt , p.force.v_sig );
+                p.force.u_dt , p.force.h_dt , p.force.v_sig ); */
             }
         }
         
     /* Dump the result. */
-    ihg = kernel_igamma / p.h;
-    printf( "pairs_single: part %lli (h=%e) has wcount=%f, rho=%f.\n" , p.id , p.h , p.density.wcount + 32.0/3 , ihg * ihg * ihg * ( p.rho + p.mass*kernel_root ) );
+    printf( "pairs_single: part %lli (h=%e) has a=[%.3e,%.3e,%.3e], udt=%e.\n" , p.id , p.h , p.a[0] , p.a[1] , p.a[2] , p.force.u_dt );
     fflush(stdout);
     
     }
@@ -472,6 +521,8 @@ void engine_step ( struct engine *e , int sort_queues ) {
             }
         }
         
+    // engine_single_density( e->s->dim , 6178 , e->s->parts , e->s->nr_parts , e->s->periodic );
+
     /* Start the clock. */
     TIMER_TIC_ND
     
@@ -488,11 +539,12 @@ void engine_step ( struct engine *e , int sort_queues ) {
     /* Stop the clock. */
     TIMER_TOC(timer_runners);
     
+    // engine_single_force( e->s->dim , 6178 , e->s->parts , e->s->nr_parts , e->s->periodic );
+    
     // for(k=0; k<10; ++k)
     //   printParticle(parts, k);
-    // engine_single( e->s->dim , 432626 , e->s->parts , e->s->nr_parts , e->s->periodic );
     // printParticle( parts , 432626 );
-    // printParticle( parts , 432628 );
+    // printParticle( e->s->parts , 6178 , e->s->nr_parts );
 
     /* Collect the cell data from the second kick. */
     for ( k = 0 ; k < e->s->nr_cells ; k++ ) {

src/kernel.h

@@ -32,16 +32,15 @@
 /* Coefficients for the kernel. */ 
 #define kernel_degree 3
 #define kernel_ivals 2
-#define kernel_gamma 0.5f
-#define kernel_igamma 2.0f
-#define kernel_igamma3 8.0f
-#define kernel_igamma4 16.0f
+#define kernel_gamma 2.0f
+#define kernel_gamma2 4.0f
+#define kernel_gamma3 8.0f
 static float kernel_coeffs[ (kernel_degree + 1) * (kernel_ivals + 1) ] __attribute__ ((aligned (16))) =
     { 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 ( kernel_coeffs[ kernel_degree ] )
-#define kernel_wroot ( 4.0/3.0*M_PI*kernel_igamma3*kernel_coeffs[ kernel_degree ] )
+#define kernel_wroot ( 4.0/3.0*M_PI*kernel_coeffs[ kernel_degree ] )
       
       
 /**

src/runner.c

@@ -335,7 +335,7 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
     struct cell *finger;
     int i, k, redo, count = c->count;
     int *pid;
-    float h, hg, ihg, ihg2, ihg4, h_corr, rho, rho_dh, u, fc;
+    float h, ih, ih2, ih4, h_corr, rho, wcount, rho_dh, wcount_dh, u, fc;
     float normDiv_v, normCurl_v;
     float dt_step = r->e->dt_step;
     TIMER_TIC
@@ -372,23 +372,21 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
             
   	            /* Some smoothing length multiples. */
 	            h = p->h;
-                hg = kernel_gamma * h;
-                ihg = 1.0f / hg;
-                ihg2 = ihg * ihg;
-                ihg4 = ihg2 * ihg2;
-
-                /* Adjust the computed weighted number of neighbours. */
-                p->density.wcount += kernel_wroot;
+                ih = 1.0f / h;
+                ih2 = ih * ih;
+                ih4 = ih2 * ih2;
 
 		        /* Final operation on the density. */
-                p->rho = rho = ihg * ihg2 * ( p->rho + p->mass*kernel_root );
-                p->rho_dh = rho_dh = p->rho_dh * ihg4;
+                p->rho = rho = ih * ih2 * ( p->rho + p->mass*kernel_root );
+                p->rho_dh = rho_dh = p->rho_dh * ih4;
+                wcount = ( p->density.wcount + kernel_root ) * ( 4.0f / 3.0 * M_PI * kernel_gamma3 );
+                wcount_dh = p->density.wcount_dh * ih * ( 4.0f / 3.0 * M_PI * kernel_gamma3 );
                     
                 /* Compute the smoothing length update (Newton step). */
-                if ( p->density.wcount_dh == 0.0f )
+                if ( wcount_dh == 0.0f )
                     h_corr = p->h;
                 else {
-                    h_corr = ( const_nwneigh - p->density.wcount ) / p->density.wcount_dh;
+                    h_corr = ( const_nwneigh - wcount ) / wcount_dh;
 
                     /* Truncate to the range [ -p->h/2 , p->h ]. */
                     h_corr = fminf( h_corr , h );
@@ -400,9 +398,9 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
                 cp->h = p->h;
 
                 /* Did we get the right number density? */
-                if ( p->density.wcount > const_nwneigh + const_delta_nwneigh ||
-                     p->density.wcount < const_nwneigh - const_delta_nwneigh ) {
-                    // printf( "runner_doghost: particle %lli (h=%e,depth=%i) has bad wcount=%.3f.\n" , p->id , p->h , c->depth , p->density.wcount ); fflush(stdout);
+                if ( wcount > const_nwneigh + const_delta_nwneigh ||
+                     wcount < const_nwneigh - const_delta_nwneigh ) {
+                    // printf( "runner_doghost: particle %lli (h=%e,depth=%i) has bad wcount=%.3f.\n" , p->id , p->h , c->depth , wcount ); fflush(stdout);
                     // p->h += ( p->density.wcount + kernel_root - const_nwneigh ) / p->density.wcount_dh;
                     pid[redo] = pid[i];
                     redo += 1;
@@ -417,8 +415,8 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
                     }
 
                 /* Pre-compute some stuff for the balsara switch. */
-		        normDiv_v = fabs( p->density.div_v / rho * ihg4 );
-		        normCurl_v = sqrtf( p->density.curl_v[0] * p->density.curl_v[0] + p->density.curl_v[1] * p->density.curl_v[1] + p->density.curl_v[2] * p->density.curl_v[2] ) / rho * ihg4;
+		        normDiv_v = fabs( p->density.div_v / rho * ih4 );
+		        normCurl_v = sqrtf( p->density.curl_v[0] * p->density.curl_v[0] + p->density.curl_v[1] * p->density.curl_v[1] + p->density.curl_v[2] * p->density.curl_v[2] ) / rho * ih4;
                 
                 /* As of here, particle force variables will be set. Do _NOT_
                    try to read any particle density variables! */
@@ -428,10 +426,10 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
                 p->force.c = fc = sqrtf( const_gamma * ( const_gamma - 1.0f ) * u );
 
                 /* Compute the P/Omega/rho2. */
-                p->force.POrho2 = u * ( const_gamma - 1.0f ) / ( rho + hg * rho_dh / 3.0f );
+                p->force.POrho2 = u * ( const_gamma - 1.0f ) / ( rho + h * rho_dh / 3.0f );
 
-		/* Balsara switch */
-		p->force.balsara = normDiv_v / ( normDiv_v + normCurl_v + 0.0001f * fc * ihg );
+		        /* Balsara switch */
+		        p->force.balsara = normDiv_v / ( normDiv_v + normCurl_v + 0.0001f * fc * ih );
                 
                 /* Reset the acceleration. */
                 for ( k = 0 ; k < 3 ; k++ )
@@ -543,7 +541,7 @@ void runner_dokick2 ( struct runner *r , struct cell *c ) {
             }
 
         /* Update the particle's time step. */
-        p->dt = pdt = const_cfl * kernel_gamma * h / ( p->force.v_sig );