Stability Condition for Viscosity

[ltlancas]

I was recently looking at implementing a hyper-diffusion operator in Athena++ and was making sure that I understood the code as it is implemented for viscosity (the explicit method, not the super TimeStepping). I was particularly making sure that I understood where the time-step criteria were coming from. For the case of constant, isotropic kinematic viscosity, $\nu$, it seems from the hydro_diffusion.cpp file (specifically here) that the time step, $\Delta t$ is limited by
$$\Delta t < \frac{\Delta x^2}{2N_d\nu}$$
where $N_d$ is the number of dimensions and $\Delta x$ is the local grid spacing. This is then minimized over each dimension, (i.e. also $\Delta y$ and $\Delta z$).

Maybe I'm missing something but I think that this choice may be unstable, at least in certain situations. For example, considering the one-dimensional case, the above time-step criterion would be appropriate for the equation
$$\frac{\partial f}{\partial t} - \nu \frac{\partial^2 f}{\partial x^2} = 0$$
which was discretized on a grid with spacing $\Delta t$ and $\Delta x$.

However, the one-dimensional reduction of the parabolic part of the momentum equation is
$$\frac{\partial \left( \rho v_x\right)}{\partial t} - \frac{\partial}{\partial x} \left[\frac{4}{3}\rho\nu \frac{\partial v_x}{\partial x} \right] = 0$$
If I discretize the above and use the first-order finite difference scheme that is implemented in Athena in order to write the equation as a time-update and perform a von Neumann stability analysis, I get a stability condition of
$$\Delta t < \frac{3 \Delta x^2}{8 \nu}$$
This basically comes from the $4/3$ factor in the momentum equation which should be easy enough to see. I think the right condition for multiple dimensions would be to replace $N_d$ with $N_d + 1/3$, though I haven't really gone through it.

Anyway, does this make sense or am I missing something?

[ltlancas]

Actually, the averaging of densities (to get density at cell interfaces) should give you another factor of 1/2 so that the stability condition is
$$\Delta t < \frac{3\Delta x^2}{4\nu}$$
which means the choice in the code should be stable, sorry for the confusion. That said, I think the choice in the code may be too conservative?