Boundary violations are deeply experienced.

― David W. Earle

The post yesterday is closely connected to this one.

Yesterday I talked about a general purpose form for “limiters” applied to polynomial reconstructions, $p \left( \theta \right)$. For methods of this type the reconstruction is important because it defines the value of the function where there is not data, at cell-edges or faces where fluxes are computed, $u_{j\pm1/2}=p\left(\pm1/2\right)$. For a method in conservation form the edge-face values are essential to updating the solution. As I will discuss the value of these values and mindful modifications of the values are severely under-utilized in implementing methods. Various methods can be improved in terms of accuracy, resolution, dissipation and stability through the techniques discussed in this post.

All of these techniques will utilize an important observation for methods of this type; one can swap edge-face values without impacting the accuracy of the method. Consider the reconstructed edge value $u_{j+1/2}$ that can (usually) be two values with a reconstruction being applied to $p\left(\theta_j \right)$ and $p\left(\theta_{j+1} \right)$. If both reconstructions have equal accuracy we could trade these values, or replace one value by the other without impacting accuracy. Other properties like stability, dissipation or phase error are effected by this change, but formal accuracy is preserved.

One of the key aspects of the use of these edge values is a model for what an upwind method looks like applied to the data. Consider a nonlinear flux $f\left(u\right)$ where upwinding will be applied to the edge values from two neighboring cells to determine the physically admissible flux, $f\left(u\right)_{j+1/2}$. An upwind approximation will use edge values from cells $j$ and $j+1$. The simplest and most useful approximation of the flux is $f_{j+1/2} =$ $\frac{1}{2} \left( f_{j,j+1/2} + f_{j+1,j+1/2} \right) - \frac{1}{2} R \lvert\lambda\rvert$ $L \left( u_{j+1,j+1/2} - u_{j,j+1/2} \right)$ where $\partial_u f = R \lambda L$ is an eigen-decomposition of the flux Jacobian. The numerical dissipation is proportional to the eigenvalues $\lambda$ and the jump at the cell edge $u_{j+1,j+1/2} - u_{j,j+1/2}$. Generally speaking the physical direction of dissipation is defined by the jump from the cell-center (average) values, $u_{j+1}$ $latex– u_j$. The extrapolated edge values can (often) change the sign of the jump meaning that the approximation is locally anti-diffusive. A picture of faces is worth showing so one can understand what each swap looks like and envision the impact.

Perhaps the first application of this idea was the contact steepener defined by Colella and Woodward in their PPM method. Material interfaces also known as contact discontinuities are infamously diffused by Eulerian finite volume schemes even ones as high-resolution as PPM. The method introduced with PPM removes most of the diffusion usually induced by the method. The method works by locating valid interfaces and in the cells with the interface swapping the interface values at the edges. Often a material interface will be dissipative and swapping the edge values will be anti-diffusive. It is essential to apply a nonlinear stability mechanism to the swapped interface values and the polynomial reconstruction so that the anti-dissipative effects are not destabilizing. A second word of caution is the ill-advised nature of applying this technique at a nonlinear discontinuity such as a shock or rarefaction. For those flow features this technique would invite disaster.

A second, far safer application of edge swapping would be increasing the dissipation. In cases where the high-order local dissipation is anti-dissipative because the extrapolated values at the edges have a differently signed jump across the interface, swapping the edge values would make this dissipative. This condition has recently been dubbed as “entropy stable”. The condition can be computed by checking whether the jump in edge values $p(-1/2)_{j+1}$ $- p(1/2)_j$ has the same sign as the jump in cell-centered values $u_{j+1} - u_j$. This should be particularly important in stably and physically moving nonlinear structures on the grid.  A concern is that too much dissipation is created as the figure shows. In this case the jump is entropy stable, but larger than the first-order method. This could actually be unstable especially if a Lax-Friedrichs flux is used because it is more dissipative than upwinding.

The penultimate approach would be to use these techniques to remove dissipation entirely from the edge. This is enabled by simply setting the two edge values to be equal in value. Any convex average of the edge values would achieve this as well as simply over-writing one of the edge values by the other. Like the anti-dissipative methods this would need to be applied with care for the stability of the overall integration procedure. In several cases the resulting approximation is actually higher order accurate. For example by applying this to first-order upwinding, the approximation is promoted to second-order central differencing. This happens for odd-ordered upstream centered methods, i.e., the third-order upwind methods is promoted to fourth-order centered edge values, the fifth-order upwind to sixth-order centered, and so forth.

The final approach to discuss would simply take one of the edge values and completely replace it by the other. The most obvious way to achieve this would be to declare that the information at the edge is unambigiously one-directional, and use that direction to choose the edge values. In this case one would short-circuit the Riemann solver and would only be stable for super-sonic flow. Of course, the opposite would produce a demonstrably anti-diffusive effect and might be useful (and dangerous too).

Don’t live a normal life by default, push the boundaries of your potential.

Colella, Phillip, and Paul R. Woodward. “The piecewise parabolic method (PPM) for gas-dynamical simulations.” Journal of computational physics 54.1 (1984): 174-201.

Yang, Huanan. “An artificial compression method for ENO schemes: the slope modification method.” Journal of Computational Physics 89.1 (1990): 125-160.

Balsara, Dinshaw S., and Chi-Wang Shu. “Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy.”Journal of Computational Physics 160.2 (2000): 405-452.

Fjordholm, Ulrik S., Siddhartha Mishra, and Eitan Tadmor. “Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws.” SIAM Journal on Numerical Analysis 50.2 (2012): 544-573.

Banks, Jeffrey William, and Jeffrey Alan Furst Hittinger. “A new class of nonlinear finite-volume methods for Vlasov simulation.” Plasma Science, IEEE Transactions on 38.9 (2010): 2198-2207.

Dumbser, Michael, and Olindo Zanotti. “Very high order P N P M schemes on unstructured meshes for the resistive relativistic MHD equations.” Journal of Computational Physics 228.18 (2009): 6991-7006.

Munz, Claus-Dieter, et al. “Enhanced Accuracy for Finite-Volume and Discontinuous Galerkin Schemes via Non-intrusive Corrections.” Recent Developments in the Numerics of Nonlinear Hyperbolic Conservation Laws. Springer Berlin Heidelberg, 2013. 267-282.