The law that entropy always increases holds, I think, the supreme position among the laws of Nature. If someone points out to you that your pet theory of the universe is in disagreement with Maxwell’s equations — then so much the worse for Maxwell’s equations. If it is found to be contradicted by observation — well, these experimentalists do bungle things sometimes. But if your theory is found to be against the second law of thermodynamics I can give you no hope; there is nothing for it but to collapse in deepest humiliation.

— Sir Arthur Stanley Eddington

AS USUAL WORD PRESS’ LaTeX is annoying.  I need to post and go to bed.

Rather than continue to talk about global issues I’ll get back in the weeds this week and get into a technical conversation. I don’t really know where this is going, so this is a bit of stream of consciousness in thinking about a topic and developing a new idea. The inspiration for this came from my talk at the Multimat 2017 meeting, and considering how to fix problems I’ve seen with rarefactions. As a reminder, I had seen that some solvers produce solutions with small, but consistent violations of the second law of thermodynamics in their solution of expansions (e.g. rarefactions). Nothing catastrophic is observed, but it is a troubling failure from nominally robust solvers. This study was itself motivated by the observation of a systematic failure of these solvers to produce convergent solutions to very strong rarefactions, and examine what sort of character the solutions have under more ideal circumstances.

A few points are worth making about the solvers used and how they have been tested in the past. Mostly I’ve worked with finite volume codes, sort of the gold standard of production codes for fluid dynamics. These codes are very reliable, and well understood.  For the most part, the focus of test problems has been shock waves where bad methods can result in catastrophic instability for the codes. Rarefactions are far less studied and tested because they are generally benign, and don’t threaten the stability of the code. As a result, the rarefaction focused test problems are mostly missing. We do know that expansions can produce unphysical solutions for Eulerian codes at critical points (where the characteristic speeds go to zero, and numerical dissipation may vanish). Bad solutions can arise with strong rarefactions, but no one has pointed out that these solutions actually violate the second law of thermodynamics before. The end result is a relative ignorance about shortcomings of the code, and a potentially important outlet for improvement of the methods.

Von Neumann told Shannon to call his measure entropy, since “no one knows what entropy is, so in a debate you will always have the advantage.

― Jeremy Campbell

One of my suggestions about examining this problem is examining the solution to these problems with a wider variety of codes. This would include codes that do not use pure finite volume methods.  One example are methods based on flux differencing where the code can go to formally high-order accuracy for nonlinear problems. Control volume codes are limited to second-order accuracy and the leading nonlinear truncation error can produce the entropy-condition violating energy transfer in expansions, $C f_{uu} u_x u_{xx}$. For almost every control volume code these terms are dissipative in shock waves, thus providing additional stability to the codes in this dangerous configuration. The opposing reaction in expansions can go unnoticed because any imperfections in the solution are modulated by the physics of the problem. For this reason, the failing has gone completely unnoticed for decades.  A reasonable question to explore is whether codes based on different design principles exhibit the same problems, or produce solutions that satisfy the second law of thermodynamics more uniformly.

An important technique in defining flux difference schemes of high-order accuracy is flux splitting (more than second-order accuracy). The core idea is that approximating the fluxes to high order can produce higher formal accuracy then the variables. The question is does this produce solutions of a fundamentally different character with respect to entropy. Simply put, a flux splitting is a decomposition of the fluxes being differenced into negative and positive moving contributions. These fluxes are then differenced and then recomposed into the total flux. The splitting techniques add a directionality to the approximation needed for numerical stability associated with upwinding the approximation. The flux splitting techniques are closely related to Riemann solvers, but here to fore only include a small number of simple linearized Riemann solutions. I’d like to explore some greater generalization of this concept including flux splittings based on exact Riemann solvers.

The Riemann problem is the exact solution to the interaction of two discontinuous states described by hyperbolic equations. This analytic information can be used to develop numerical methods that encode this physically relevant information into the solution. In terms of numerical methods, the Riemann solution is a generalization of the principle of upwinding, where the physical direction of propagation is taken into account. The first person to describe this approach to numerical methods was SK Godunov in 1959. Godunov’s method was first-order accurate and used the exact solution to the Riemann problem. It was soon realized that one only needed to approximate the Riemann solution. This became a key development in the methods over time and allowed great progress. Over time it was realized that it also allowed great flexibility too.

In science if you know what you are doing you should not be doing it. In engineering if you do not know what you are doing you should not be doing it. Of course, you seldom, if ever, see the pure state.

– Richard Hamming

A simple Riemann solver can be defined by linearizing the problem, $f(u_l,u_r) = \frac{1}{2} \left[ f_l + f_r \right]$ $- R | \lambda | L$ $( u_r – u_l )$. The quantity $f$ is the nonlinear flux, $u_l, u_r$ are the states to the left and right of the interface. The dissipation is defined by the eigen-decomposition of the flux Jacobian, $\partial_u f = A = \partial_u f = R \lambda L$. This decomposition is contained of the right and left eigenvectors and the eigenvalues, $\lambda$. These eigenvalues are the characteristic velocities, which for gas dynamics are $u-c, u, u+c$ being the velocities and the sound speeds, $c$. This basic decomposition is the basis of flux splitting techniques.

The basic flux splitting takes the flux and decomposes it into right and left moving pieces, $f(u) = f(u)^- + f(u)^+$. One was to do this is choose a velocity, $\alpha > 0$, and create contributions where $f(u)^+ = \frac{1}{2} \left[ f(u) + \alpha u\right]$ and $f(u)^- = \frac{1}{2}\left[ f(u) - \alpha u\right]$. A simple choice of $\alpha = \Delta x/\Delta t$ creates the Lax-Friedrichs flux, the simplest (and most dissipative) Riemann solver. For the general linearized Riemann solver the flux splitting is $f(u)^+ = \frac{1}{2}\left[ f(u) + R |\lambda | L u\right]$ and $f(u)^- = \frac{1}{2}\left[ f(u) - R | \lambda | L u\right]$. The choice of the left and right states to evaluate the flux Jacobian defines the flux splitting. For example, if the states are evaluated using Roe’s recipe, we get the Roe flux splitting. If we evaluate the eigenvalues in a bounding fashion we can get the local Lax-Friedrichs method.

Another approach to generating a flux splitting does not use the variables in the expression of splitting, and only uses the fluxes. In this case the expressions are developed in terms of the sign of the eigenvalues/characteristic velocity. The splitting then works only as a scaling by eigenvector decomposition of the flux Jacobian. The expressions are somewhat simplified, as $f^+ = \frac{1}{2} \left(f + R \mbox{sign}(\lambda) L \right) f$ and $f^- = \frac{1}{2} ( f – R \mbox{sign(\lambda) L ) f$. We not in passing that the smooth or soft version of the sign function might be extremely useful in this type of splitting and introducing a continuously differentiable function (https://wjrider.wordpress.com/2017/03/24/smoothed-operators/ https://wjrider.wordpress.com/2017/03/24/smoothed-operators/). By the same token, the absolute value used in the usual flux splitting approach could also be smoothed to similar effect. We need to take care in our choices to assure that the accuracy of the resulting numerical method is not negatively impacted. We get into some very big problems when we want to generalize to other Riemann solvers. Examples of these solvers are the HLL family of solvers, and the most classical Riemann solver, the exact solver or close approximations to that approach (e.g., a single iteration of the nonlinear Newton’s method used in the exact solver). How can these important methods be utilized in flux splitting methods? For very strong wave interactions these classes of methods are extremely valuable and not presently possible to be used effectively in flux splitting.

Nature not only suggests to us problems, she suggests their solution.

—Henri Poincare´

We can start with the simpler case of the HLL type of flux, which has an algebraic description. The HLL flux is defined using the space-time diagram by integrating the equations to derive a flux. The simplest form of the flux uses bounds for the wave speeds and neglecting all of the structure inside the Riemann fan resulting in a simple closed form expression for the flux, $f_{lr} = \left[a_r f_l – a_l f_r + a_l a_r \left( u_r – u_l \right)\right]/(a_r – a_l)$. The flux is quite simple, but dependent on the estimates for the smallest and largest wave speeds in the system. The left wave speed, $a_l$ is the smallest wave speed and needs to be bounded at zero (i.e., it is negative). The right most wave speed is $a_r$ and is bounded below by zero. The HLL flux has the benefit of reducing to simple upwind flux for the system if all the wave speeds are either negative or positive. For a flux splitting we need to take this apart into negative and positive moving pieces for the purposes of splitting nearby fluxes as we did with the Roe, or flavors of Lax-Friedrichs.

The flux splitting can be defined almost by inspection. The positive flux is $f^+ = (a_r f – a_l a_r u) / (a_r – a_l)$. The negative flux is $f^- = (- a_l f + a_l a_r u) / (a_r – a_l)$. This is a wonderfully simple result, and meets all the basic requirements for a flux splitting. Unfortunately, the HLL flux is extremely dissipative, thus lacking some degree of practical utility. Still we expect this flux splitting to be quite robust especially for strong waves with the proviso that the wave speed estimates bound the physical wave speeds. This is a much more delicate estimate than usually recognized. The case of a reflected wave can produce wave speeds nonlinearly that exceed the wave speeds in the initial data.

The harder case is the class of exact Riemann solvers that are defined algorithmically and do not have a closed form. After using an exact Riemann solver we do have a set of initial left and right states, and the resolved state at the centering point $x/t=0$. If we desire a flux splitting, it needs to be defined in terms of these variables. The trick in this endeavor is choosing an algebraic structure to help produce a workable flux splitting technique. We build upon the experience of the HLL flux partially because we can incorporate the knowledge arising from the exact solution into the algebraic structure to good effect. In particular, the nature of the one-sided differencing can be reproduced effectively. This requires the wave speed bounds to use the interior states of the solution.

The exact flux is different than the HLL flux, and this will be defined by changing the dissipation vector in the flux. Our chosen structure is a flux defined by $f_{lr} = \(a_r f_l – a_l f_r – D \left( u_r – u_l \right)) / (a_r – a_l)$. If we can derive the form for $D$ our work will be done. The positive flux is $f^+ = (a_r f + D u) / (a_r – a_l})$. The negative flux is $f^- = (- a_l f – D u) / (a_r – a_l)$. Now we just have a little bit of algebra to arrive at our final expression. The math is nice and straightforward, $D = (a_r f_l – a_l f_r – (a_r – a_l) f_{lr} ) / (u_r – u_l)$. A couple comments are needed at this point. When the states become equal, the solver becomes ill defined, $u_l = u_r$. Fortunately, this is exactly where the linearized flux splitting approaches or HLL would be ideal.

The secret to being wrong isn’t to avoid being wrong! The secret is being willing to be wrong. The secret is realizing that wrong isn’t fatal.

― Seth Godin

Godunov, S. K. “A finite difference method for the computation of discontinuous solutions of the equations of fluid dynamics.” Sbornik: Mathematics 47, no. 8-9 (1959): 357-393.

Van Leer, Bram. “Flux-vector splitting for the Euler equations.” In Eighth international conference on numerical methods in fluid dynamics, pp. 507-512. Springer Berlin/Heidelberg, 1982.

Harten, Amiram, Peter D. Lax, and Bram Van Leer. “On upstream differencing and Godunov-type schemes for hyperbolic conservation laws.” In Upwind and High-Resolution Schemes, pp. 53-79. Springer Berlin Heidelberg, 1997

Shu, Chi-Wang, and Stanley Osher. “Efficient implementation of essentially non-oscillatory shock-capturing schemes, II.” Journal of Computational Physics 83, no. 1 (1989): 32-78..

Jiang, Guang-Shan, and Chi-Wang Shu. “Efficient implementation of weighted ENO schemes.” Journal of computational physics126, no. 1 (1996): 202-228.

Quirk, James J. “A contribution to the great Riemann solver debate.” International Journal for Numerical Methods in Fluids 18, no. 6 (1994): 555-574.