The code now implements Gadget-2 type viscosity in the form of a

Monaghan-1997 term with a Balsara-1995 switch.



Former-commit-id: e3be6134bf89c532f76bad1f933eade4ceb95fea

examples/SedovBlast/makeIC.py

@@ -63,7 +63,7 @@ for i in range(L):
             coords[index,0] = x
             coords[index,1] = y
             coords[index,2] = z
-            v[index,0] = 0.
+            v[index,0] = 0.1
             v[index,1] = 0.
             v[index,2] = 0.
             m[index] = mass

examples/test.c

@@ -959,6 +959,13 @@ int main ( int argc , char *argv[] ) {
         
         /* Take a step. */
         engine_step( &e , 0 );
+
+	if(j % 10 == 0)
+	  {
+	    char fileName[200];
+	    sprintf(fileName, "output_%05i.hdf5", j); 
+	    write_output(fileName, dim, parts, N, periodic);
+	  }
         
         /* Dump the first few particles. */
         for(k=0; k<10; ++k)

src/engine.c

@@ -296,13 +296,12 @@ void engine_step ( struct engine *e , int sort_queues ) {
 
             /* Scale the derivatives if they're freshly computed. */
             if ( p->dt <= dt_step ) {
-                p->u_dt *= p->POrho2;
                 p->h_dt *= p->h * 0.333333333f;
                 lcount += 1;
                 }
                 
             /* Update the particle's time step. */
-            p->dt = const_cfl * p->h / ( p->c + p->v_sig );
+            p->dt = const_cfl * p->h / (  p->v_sig );
 
             /* Update positions and energies at the half-step. */
             // #ifdef __SSE__

src/kernel.h

@@ -45,7 +45,7 @@ static float kernel_coeffs[ (kernel_degree + 1) * (kernel_ivals + 1) ] __attribu
       
       
 /**
- * @brief Computes the kernel and its derivative for a given distance x.
+ * @brief Computes the kernel and its derivative for a given distance x. Gives a sensible answer only if x<1.
  */
 
 __attribute__ ((always_inline)) INLINE static void kernel_deval ( float x , float *W , float *dW_dx ) {
@@ -65,7 +65,7 @@ __attribute__ ((always_inline)) INLINE static void kernel_deval ( float x , floa
 #ifdef VECTORIZE
 
 /**
- * @brief Computes the kernel and its derivative for a given distance x (Vectorized version)
+ * @brief Computes the kernel and its derivative for a given distance x (Vectorized version). Gives a sensible answer only if x<1.
  */
 
 __attribute__ ((always_inline)) INLINE static void kernel_deval_vec ( vector *x , vector *w , vector *dw_dx ) {
@@ -96,7 +96,7 @@ __attribute__ ((always_inline)) INLINE static void kernel_deval_vec ( vector *x
 #endif
 
 /**
- * @brief Computes the kernel for a given distance x
+ * @brief Computes the kernel for a given distance x. Gives a sensible answer only if x<1.
  */
 
 __attribute__ ((always_inline)) INLINE static void kernel_eval ( float x , float *W ) {

src/runner_iact.h

@@ -357,6 +357,7 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_force ( float r2
     float hi_inv, hi2_inv;
     float hj_inv, hj2_inv;
     float wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, w, dvdr;
+    float v_sig, omega_ij, Pi_ij;
     float f;
     int k;
     
@@ -373,32 +374,44 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_force ( float r2
     xj = r * hj_inv * kernel_igamma;
     kernel_deval( xj , &wj , &wj_dx );
     wj_dr = hj2_inv * hj2_inv * wj_dx;
-        
+                
+    /* Compute dv dot r. */
+    dvdr = ( pi->v[0] - pj->v[0] ) * dx[0] + ( pi->v[1] - pj->v[1] ) * dx[1] + ( pi->v[2] - pj->v[2] ) * dx[2];
+    dvdr *= ri;
+
+    /* Compute the relative velocity. (This is 0 if the particles move away from each other and negative otherwise) */
+    omega_ij = fminf(dvdr, 0.f);
+    
+    /* Compute signal velocity */
+    v_sig = pi->c + pj->c - 3.*omega_ij;
+
+    /* Compute viscosity tensor */
+    Pi_ij = -const_viscosity_alpha * v_sig * omega_ij / (pi->rho + pj->rho);
+
+    /* Apply Balsara switch */
+    Pi_ij *= (pi->Balsara + pj->Balsara);
+
     /* Get the common factor out. */
-    w = ri * ( pi->POrho2 * wi_dr + pj->POrho2 * wj_dr );
-        
+    w = ri * ( ( pi->POrho2 * wi_dr + pj->POrho2 * wj_dr ) - 0.25f * Pi_ij * ( wi_dr + wj_dr ) );
+
     /* 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;
         }
-        
-    /* Compute dv dot r. */
-    dvdr = ( pi->v[0] - pj->v[0] ) * dx[0] + ( pi->v[1] - pj->v[1] ) * dx[1] + ( pi->v[2] - pj->v[2] ) * dx[2];
-    dvdr *= ri;
-        
+                
     /* Get the time derivative for u. */
-    pi->u_dt += pj->mass * dvdr * wi_dr;
-    pj->u_dt += pi->mass * dvdr * wj_dr;
+    pi->u_dt += pi->POrho2 * pj->mass * dvdr * wi_dr + 0.125f * pj->mass * Pi_ij * dvdr * ( wi_dr + wj_dr );
+    pj->u_dt += pj->POrho2 * pi->mass * dvdr * wj_dr + 0.125f * pi->mass * Pi_ij * dvdr * ( wi_dr + wj_dr );
     
     /* Get the time derivative for h. */
     pi->h_dt -= pj->mass / pj->rho * dvdr * wi_dr;
     pj->h_dt -= pi->mass / pi->rho * dvdr * wj_dr;
     
     /* Update the signal velocity. */
-    pi->v_sig = fmaxf( pi->v_sig , pj->c - 3*dvdr );
-    pj->v_sig = fmaxf( pj->v_sig , pi->c - 3*dvdr );
+    pi->v_sig = fmaxf( pi->v_sig , v_sig );
+    pj->v_sig = fmaxf( pj->v_sig , v_sig );
     
     #ifdef HIST
     if ( hi > hj )

theory/latex/swift.tex

@@ -579,7 +579,7 @@ The viscosity tensor $\Pi_{ij}$ is given by:
 
 \begin{eqnarray*}
  \Pi_{ij} &=& -\alpha \frac{\left(c_i + c_j - 3w_{ij}\right)w_{ij}}{\rho_i + \rho_j} \\
- w_{ij} &=& \min\left(0, \frac{(\vec{v}_i - \vec{v}_j)\cdot(\vec{r}_i - \vec{r}_j}{|r_{ij}|}\right)
+ w_{ij} &=& \min\left(0, \frac{(\vec{v}_i - \vec{v}_j)\cdot(\vec{r}_i - \vec{r}_j)}{|r_{ij}|}\right)
 \end{eqnarray*}
 
 The $w_{ij}$ term can be re-used to express the signal velocity \ref{eq:sigvel}. The viscosity parameter $\alpha$