diff --git a/examples/SedovBlast/makeIC.py b/examples/SedovBlast/makeIC.py
index 6a142334ed314407e8cf9884440992e1d02a9164..723b70396b55e832d67dd134b2610d2e21bbff3b 100644
--- a/examples/SedovBlast/makeIC.py
+++ b/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
diff --git a/examples/test.c b/examples/test.c
index 57669065b088801601a7f10a1d0eb8464c8ba737..9adf85d901cb94967969f8080d021ca17271cb87 100644
--- a/examples/test.c
+++ b/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)
diff --git a/src/engine.c b/src/engine.c
index 0020dfb2879cd6c38d4544a695de130809f4337e..cf9f564c61d1b59b1448b2dcae193805fb3615ad 100644
--- a/src/engine.c
+++ b/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__
diff --git a/src/kernel.h b/src/kernel.h
index 493eb8f910bf1324e71e7eb4abffe96c5ce34215..09781ee080263eb6e69694cdff30a8afd54e5b0b 100644
--- a/src/kernel.h
+++ b/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 ) {
diff --git a/src/runner_iact.h b/src/runner_iact.h
index 43192c952d7c7749a8ba06f0c55fb1ac4105df98..107b3a1f02845bbeb9803d237259fb7fe558bcee 100644
--- a/src/runner_iact.h
+++ b/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 )
diff --git a/theory/latex/swift.tex b/theory/latex/swift.tex
index 184943a21aef5401236255a4bcb87296bb33a33a..f489f5896b60a6848e0a0b7a6966587f1402b2ed 100755
--- a/theory/latex/swift.tex
+++ b/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$