No man needs a vacation so much as the man who has just had one.

― Elbert Hubbard

A couple of things before I jump in: I was on vacation early this week, and by vacation I mean, vacation, nothing from work was happening, it was Paris and that’s too good to waste working, and the title of the post is partially right, its actually a situation that is much worse than that.

Among the most widely held and accepted facts for solving a hyperbolic PDE explicitly is the time step estimate, and it is not good enough. Not good enough to be relied upon for production work. That’s a pretty stunning thing to say you might think, but its actually pretty damn obvious. Moreover this basic fact also means that the same bounding estimates aren’t good enough to insure proper entropy production either since they are based on the same basic thing. In short order we have two fundamental pieces of the whole solution technology that can easily fall completely short of what is needed for production work.

It is probably holding progress back as more forgiving methods will be favored in the face of this shortcoming. What is more stunning is how bad the situation can actually be.

If you don’t already know, what is the accepted and standard approach? The time step size or dissipation is proportional to the characteristic speeds in the problem and then other specifics of a given method. The characteristic speeds are the absolute key to everything. For simple one dimensional gas dynamics, the characteristics are $u\pm c$ and $u$ where $u$ is the fluid velocity and $c$ is the sound speed. A simple bound can be used as $\lambda_{\max}=|u|+c$. The famous CFL or Courant condition gives a time step of $\Delta t = \min(A \Delta x/\lambda_{\max})$, where $A$ is a positive constant typically less than one. Then you’re off to the races and computing solutions to gas dynamics.

For dissipation purposes you look for the largest local wave speeds for a number of simple Riemann solvers, or other dissipation mechanism. This can be done locally or globally. If you choose the wave speed large enough, you are certain to meet the entropy condition.

The problem is that this way of doing things is as pervasive as it is wrong. The issue is that it ignores nonlinearities that can create far larger wave speeds. This happens on challenging problems and at the most challenging times for codes. Moreover the worse situation for this entire aspect of the technology happens under seemingly innocuous circumstances, a rarefaction wave, albeit it becomes acute for very strong rarefactions. The worst case is the flow into a vacuum state where the issue can become profoundly bad. Previously I would account for the nonlinearity of sound speeds for shock waves by adding a term similar to an artificial viscosity to the sound speed to account for this. This is only proportional to local velocity differences, and does not account for pressure effects. This approach looks like this for an ideal gas, $C = Co + \frac{\gamma + 1}{2}|\Delta u|$. This helps the situation a good deal and makes the time step selection or the wave speed estimate better, but far worse things can happen (the vacuum state issue mentioned above).

Here is the problem that causes this issue to persist. To really get your arms around this issue requires a very expensive solution to the Riemann problem, the exact Riemann solver. Moreover to get the bounding wave speed in general requires a nonlinear iteration using Newton’s method, and checking my own code 15 iterations of the solver. No one is going to do this to control time steps or estimate wave speeds.

Yuck!

I will postulate that this issue can be dealt with approximately. A first step would be to include the simplest nonlinear Riemann solver to estimate things; this uses the acoustic impedances to provide an estimate of the interior state of the Riemann fan. Here one puts the problem in a Lagrangian frame and solve the approximate wave curves, on the left side $-\rho_l C_l (u_* - u_l) = p_* - p_l$ and on the right side $\rho_r C_r$ $(u_r -u_*)$ $= p_r - p_*$ solving for $u_*$ and $p_*$. Then use the jump conditions for density, $1/\rho_* - 1/\rho = (u_l-u_*)/\rho_l C_l$ to find the interior densities (similar to the right). Then compute the interior sound speeds to get the bounds. The problem is that for a strong shock or rarefaction this approximation comes up very short, very short indeed.

The way to do better is to use a quadratic approximation for the wave speeds where the acoustic impedance changes to account for the strength of the interaction. I used these before to come up with a better approximate Riemann solver. The approximation is largely the same as before except, $\rho C_o \rightarrow \rho C_o + \rho s (u_* - u_o)$. Now you solve a quadratic problem, which is still closed form, but you have to deal with an unphysical root (which is straightforward using physical principles). For the strong rarefaction this still doesn’t work very well because the wave curve has a local minimum at a much lower velocity than the vacuum velocity.

I think this can be solved, but I need to work out the details and troubleshoot. Aside from this the details are similar to the description above.

Rider, William J. “An adaptive Riemann solver using a two-shock approximation.” Computers & fluids 28.6 (1999): 741-777.