

― E.B. White

Yesterday I looked at a simple question regarding Riemann solvers. The conclusion of this brief investigation was that more examination is warranted. A large part of the impetus for the question comes from a recent research emphasis on positivity preserving discretizations. This means that quantities that are physically positive like density and energy are garuanteed to be positive in numerical solutions.

These have always focused on the variables being solved for. The fluxes used to build the solution procedure have been examined as producing physically-admissible solutions. I noticed that in least one case, the momentum flux, the flux should be positive-definite and its numerical approximation might not be. I’ll provide a little background on the thought process that got me to the question. The point I’ll make is much more exotic and esoteric in nature.

Back when I was in Los Alamos I did work on the nature of the truncation error using a technique called “modified equations”. This technique uses the venerated Taylor series expansion to describe the order of the approximation error in a numerical method. Unlike many analysis methods, which are limited to linear equations, the modified equation method gives results for nonlinear equations. One of things I noticed in the process of analysis explains a common problem with upwind differencing: entropy violating solutions to expansion waves often called rarefaction shocks.apjs307031f1_lr

For normal fluid dynamics shocks occur upon compression, and rarefactions occur upon expansion. If an expanding flow has a shock wave it is unphysical. In cases where wave speed used in upwinding goes through zero, the dissipation in upwinding goes to zero as the dissipation is directly proportional to the wave speed. This happens at “sonic points” where the velocity is equal to the sound speed, u\pm c=0.

We can see what happens using the modified equation approach for upwinding. Take a general discretization for upwind differing in conservation form, \Delta x \partial_x f(u) \approx f_{j+1/2} - f_{j-1/2}. We use the upwind approximation f_{j+1/2} = \frac{1}{2} \left( f(u)_j + f_{j+1} \right) - \frac{1}{2} \left| \partial_u f\right|\left( u_{j+1}-u_j \right). We plug all of this into the equations and expand u in a Taylor series, u(x+j\Delta x) = u(x) + \Delta x \partial_x u + \frac{1}{2} (\Delta x)^2 \partial_{xx}u +\ldots.

When we plunge into the math, and simplify we get some really interesting results, \partial_x f(u) \approx \partial_x\left[ f(u) +\frac{\Delta x}{2}\left| \partial_u f\right| \partial_x u +\frac{(\Delta x)^2}{6}\left( \partial_{uu}f (u_x)^2 + \partial_u f u_{xx}\right)\right]. Here is the key to rarefaction shocks; when the wavespeed, \partial_u f is near zero, the dissipation is actually governed by a higher order term, \partial_{uu} f (\partial_x u)^2.It is also notable that the dissipation aids the upwind dissipation for compressions, and thus shock waves. Anti-dissipation is a shock wave would be utterly catastrophic.

Usually fluids are convex, \partial_{uu} f>0 thus when u_x>0 the term in question is anti-dissipative. This is intrinsically unphysical. The dissipation from upwinding from the lower order term proportional to \Delta x is not large enough to offset the anti-dissipation. This produces the conditions needed for a rarefactions shock wave. I’ve worried that these effects can consistently cause problems in solutions and undermine the entropy satisfying solutions. These terms were not considered in the basic theory of upwinding introduced by Harten, Hyman and Lax using modified equation analysis [HHL]. What makes matters worse is that the anti-dissipative terms will dominate in expansions when we take the upwind approximation to second-order.

The most complicated skill is to be simple.

― Dejan Stojanovic

[HLL] Harten, Amiram, James M. Hyman, Peter D. Lax, and Barbara Keyfitz. “On finite‐difference approximations and entropy conditions for shocks.” Communications on Pure and Applied Mathematics 29, no. 3 (1976): 297-322.