Entropy viscosity (EV) method for solving nonlinear conservation laws

Why EV?

Spectral methods are well known for efficient performance: They only need a relative low number of degrees of freedom to be able to approximate smooth solutions for a given accuracy compared to finite difference codes. However, spectral methods are ill-suited to approximate non smooth problems. This is because the Gibbs and Runge phenomena start to spoil the convergence rate and accuracy in these cases. Nevertheless, problems of interest often posses discontinuous solutions or solutions with insufficient analyticity, e.g. non smooth pressure distribution of matter that indicates the surface of a star (contact discontinuity).

Over the years different approaches have been developed that aim to circumvent the convergence problems of spectral methods. Among those methods is the EV method (guermond2008?) which is the topic of this project.

How does it work?

Consider a 1D conservation law of the form \[ \partial_{t} u + \partial_{x} f(u) = 0, \] where \(f(u)\) is the flux function. If \(f(u) = \frac{1}{2}u^2\), then one obtains Burgers’ equation. This equation is well known to develop discontinuous solutions, e.g. when evolving a sine wave over a periodic grid.

The EV method is a modification of the much earlier established artificial viscosity (AV) approach. The idea of the AV method is to augment a given conservation law with an additional term that contains a second order spatial derivative times a small constant, e.g. \[ \partial_{t} u + \partial_{x} f(u) = \partial_x (\nu \partial_x u), \] \(\nu\) is called viscosity, because this additional term mimics a diffusive contribution. The addition of this diffusive term causes the solution to smooth out over time. If \(\nu \rightarrow 0\) then the initial equation is recovered.

To approximate the unmodified PDE as good as possible, the AV and EV method model the viscosity \(\nu\) as a nonlinear function depending on \(u\). The difference between AV and EV is how \(\nu\) is constructed. Focusing on the EV method, \(\nu\) is constructed based on the so-called entropy inequality \[ \partial_{t} E(u) + f'(u) \partial_{x} E(u) \leq 0, \] where \(E\) is termed entropy (in reference to physical problems). The above inequality arises as a necessary condition for the selection of a unique solution from studying the weak form of a 1D conservation law, cf. (leveque1992?). The essence of this inequality is that \(E_t + f' E_x\) should be only non zero in regions with very steep gradients and discontinuities. This makes it a good candidate for a shock detector, because it flags potential troubled cells. Constructing \(\nu\) based on the entropy inequality should then lead to a smoothing of the solution only in those regions where it should be necessary.

Numerical implementation

Implementing the EV method is straight forward: Firstly, one has to code a convection-diffusion problem. This is described in some text books, cf. (hesthaven2007?) for a DG implementation. Given a working convection-diffusion problem, one can apply the EV method by simply computing the entropy viscosity \(\nu\) before evaluating the time derivative of \(u\).

According to (guermond2008?), (zingan2013?), in the 1D case there are two parameters that need to be adjusted. One is used to limited the viscosity to a maximum and to avoid too extensive diffusion. Another one is used to adjust the sensitivity of the method.

First results

Below follow plots of results of test calculations that show that the method does work and improve simulation results.

Pseudospectral method

Figure shows the time evolution of sine wave initial data computed with a Pseudospectral approximation and the EV method turned on. This is a test case to check if the entropy viscosity does not introduce undesired diffusion for smooth problems. The plotted solution qualitatively agrees with the one computed without EV (and also the analytic solution).

Figure shows the time evolution of sine wave initial data computed with a Pseudospectral approximation and the EV method turned off. As one can see the solution is dominated by oscillations. In figure the result of the same calculation but with enabled EV is plotted. Clearly, in this case EV does its job and avoids spurious oscillations. This result also qualitatively agrees with the analytic solution of the 1D Burgers’ equation (TODO provide reference).

Discontinuous Galerkin method

Figure shows the time evolution of a theta function initial data computed with a DG approximation. The sharp edged rectangle shows the initial data and the oscillating line corresponds to the evolved data. The analytic solution to this problem would agree with the initial data after it has been evolved for \(T=1\), which corresponds to one full revolution in the periodic domain. However, we observe that if EV is not turned, undesired oscillations are present. Comparing this result with figures , , we see that EV does its job. The difference between the last two plots is the difference in resolution. Increasing the polynomial order \(N\) and the number of grid elements \(K\) shows that the shocks are resolved sharper.

Figures , show the time evolution of sine wave initial date computed with a DG approximation. Both plots show a rather sharped edged curve as the final state at \(T = 0.6\). These results differ from the ones discussed in the previous section, because there the time evolution has only been carried out till \(T = 0.25\), which is when the shock is fully developed. The addition of EV allows to evolve data even if shocks have already been developed. The presented solution is indeed the correct one, cf. (yu2018?) (they show results for \(T = 0.3\) for a flux function \(f(u) = u^2\)). When comparing the two plots directly one notices that some oscillations remain for the case where \(N\) is even and \(K\) is odd. Such differences are also observed in calculations without EV and are closely connected with the parity of \(N, K\), e.g. with the grid point placement and cell division of the computational domain. This needs further investigation.

Summary and current status

I have briefly outlined the motivation and basic idea that underlies the EV method. I have then presented first tests and results that confirm that the EV method does work and gives considerable improvement for problems with discontinuous data.

The above results were produced using a Julia package evm.jl I wrote. At the time of writing these notes the package is capable of solving - 1D advection equation, - 1D Burgers’ equation, - general (nonlinear) conservation laws (not tested) with periodic boundary conditions and the EV method.

What is next?

  • Quantify convergence rates for 1D problems.
  • Solve (special relativistic) Euler’s equation in 1D.
  • Make comparison with Yorgos’ finite difference results.
  • Implement EV method in BAMPS:
    • Study 3D problems (TOV star).
  • Investigate alternative artificial viscosity approaches
    • Averaged modal decay model, (klockner2011?),
    • Other methods discussed in (yu2017?).