A Local Extrapolation Method for Hyperbolic Conservation Laws. I. The ENO Underlying Schemes

In this paper, we introduce a local extrapolation method (LEM) for the essentially non-oscillatory (ENO) schemes solving nonlinear hyperbolic conservation laws. The method extrapolates the numerical flux of the underlying scheme so that it keeps conservativity. We use a minmod type limiter to avoid spurious oscillations. We propose a new balancing technique that preserves the symmetry of a symmetric wave that works well for a wide range of CFL numbers. We also introduce two artificial compression procedures to the LEM which yield sharp resolutions of contact discontinuities. Numerical examples are presented to illustrate the performance of the method.     

Analysis of a high-order finite difference detector for discontinuities

High-order finite difference discontinuity detectors are essential for the location of discontinuities on discretized functions, especially in the application of high-order numerical methods for high-speed compressible flows for shock detection. The detectors are used mainly for switching between numerical schemes in regions of discontinuity to include artificial dissipation and avoid spurious oscillations. In this work a discontinuity detector is analysed by the construction of a piecewise polynomial function that incorporates jump discontinuities present on the function or its derivatives (up to third order) and the discussion on the selection of a cut-off value required by the detector. The detector function is also compared with other discontinuity detectors through numerical examples.     

A Shock-Capturing Algorithm for the Differential Equations of Dilation and Erosion

Dilation and erosion are the fundamental operations in morphological image processing. Algorithms that exploit the formulation of these processes in terms of partial differential equations offer advantages for non-digitally scalable structuring elements and allow sub-pixel accuracy. However, the widely-used schemes from the literature suffer from significant blurring at discontinuities. We address this problem by developing a novel, flux corrected transport (FCT) type algorithm for morphological dilation/erosion with a flat disc. It uses the viscosity form of an upwind scheme in order to quantify the undesired diffusive effects. In a subsequent corrector step we compensate for these artifacts by means of a stabilised inverse diffusion process that requires a specific nonlinear multidimensional formulation. We prove a discrete maximum–minimum principle in this multidimensional framework. Our experiments show that the method gives a very sharp resolution of moving fronts, and it approximates rotation invariance very well.     

Efficient numerical approximation of compressible multi-material flow for unstructured meshes

We consider the numerical approximation of multi-dimensional multi-material flows. This is a difficult topic related to the numerical smearing of contact discontinuities (or material interfaces or slip lines). Any Eulerian scheme will produce an artificial mixing zone. In this artificial mixture, the computation of thermodynamical variables (pressure, sound speed, temperature, …) is difficult to achieve correctly. For the stiff cases considered in this paper, with solid-liquid-gas interfaces, small errors on the thermodynamical variables lead to the blow up of the computation in many cases. In this paper, we review and explain this problem, then we propose solutions for 1D and 2D flows. Contrarily to front-tracking techniques, our schemes use the same formulation everywhere on the mesh. We provide several examples with the stiffened gas equation of state (EOS): inert materials (water-air), and chemically reactive materials (solid explosive-Plexiglass-air).     

An Antidiffusive Entropy Scheme for Monotone Scalar Conservation Laws

In a recent work J. Sci. Comput. 16, 479–524 (2001), B. Després and F. Lagoutière introduced a new approach to derive numerical schemes for hyperbolic conservation laws. Its most important feature is the ability to perform an exact resolution for a single traveling discontinuity. However their scheme is not entropy satisfying and can keep nonentropic discontinuities. The purpose of our work is, starting from the previous one, to introduce a new class of schemes for monotone scalar conservation laws, that satisfy an entropy inequality, while still resolving exactly the single traveling shocks or contact discontinuities. We show that it is then possible to have an excellent resolution of rarefaction waves, and also to avoid the undesirable staircase effect. In practice, our numerical experiments show second-order accuracy.     

Convection of Concentrated Vortices and Passive Scalars as Solitary Waves

A new version of a computational method, Vorticity Confinement, is described. Vorticity Confinement has been shown to efficiently treat thin features in multi-dimensional incompressible fluid flow, such as vortices and streams of passive scalars, and to convect them over long distances with no spreading due to numerical errors. Outside the features, where the flow is irrotational or the scalar vanishes, the method automatically reduces to conventional discretized finite difference fluid dynamic equations. The features are treated as a type of weak solution and, within the features, a nonlinear difference equation, as opposed to finite difference equation, is solved that does not necessarily represent a Taylor expansion discretization of a simple partial differential equation (PDE). The approach is similar to artificial compression and shock capturing schemes, where conservation laws are satisfied across discontinuities. For the features, the result of this conservation is that integral quantities such as total amplitude and centroid motion are accurately computed. Basically, the features are treated as multi-dimensional nonlinear discrete solitary waves that live on the computational lattice. These obey a confinement relation that is a generalization to multiple dimensions of 1-D discontinuity capturing schemes. A major point is that the method involves a discretization of a rotationally invariant operator, rather than a composition of separate 1-D operators, as in conventional discontinuity capturing schemes. The main objective of this paper is to introduce a new formulation of Vorticity Confinement that, compared to the original formulation, is simpler, allows more detailed analysis, and exactly conserves momentum for vortical flow. First, a short critique of conventional methods for these problems is given. The basic new method is then described. Finally, analysis of the new method and initial results are presented.     

A new all-speed flux scheme for the Euler equations

Since being proposed, the HLLEM-type schemes have been widely used because they are with high discontinuity resolutions and can be easily applied to the other system of hyperbolic conservation law. In this paper, we conduct theoretical analyses on the HLLE-type schemes’ performances at low speeds. By realizing that the excessive numerical dissipations corresponding to the velocity-difference terms of the momentum equations make these schemes incapable of obtaining physical solutions at low speeds, we adopt the function $g$ to control such dissipation. Also, we borrow the HLLEMS scheme’s construction and damp the shear waves in the vicinity of the shock to avoid the shock anomaly’s appearance. The moving contact discontinuity case and the Sod shock tube case show that the HLLEMS-AS scheme we propose in this paper can capture contact discontinuities and shocks as sharply as HLLEMS scheme. The Quirk’s odd–even test case and the hypersonic inviscid flow over a cylinder case demonstrate that HLLEMS-AS is robust against the shock anomaly. The inviscid low-speed flow around the NACA0012 airfoil case indicates that HLLEMS-AS is with a high resolution at low speeds. The turbulent flow past a backward facing step case demonstrates the shear wave capturing ability of the HLLEMS-AS scheme. These properties suggest that HLLEMS-AS is promising to be widely used in both cases of low speed and high speed.     

Multiplicative Operator Splittings in Nonlinear Diffusion: From Spatial Splitting to Multiple Timesteps

Operator splitting is a powerful concept used in many diversed fields of applied mathematics for the design of effective numerical schemes. Following the success of the additive operator splitting (AOS) in performing an efficient nonlinear diffusion filtering on digital images, we analyze the possibility of using multiplicative operator splittings to process images from different perspectives.We start by examining the potential of using fractional step methods to design a multiplicative operator splitting as an alternative to AOS schemes. By means of a Strang splitting, we attempt to use numerical schemes that are known to be more accurate in linear diffusion processes and apply them on images. Initially we implement the Crank-Nicolson and DuFort-Frankel schemes to diffuse noisy signals in one dimension and devise a simple extrapolation that enables the Crank-Nicolson to be used with high accuracy on these signals. We then combine the Crank-Nicolson in 1D with various multiplicative operator splittings to process images. Based on these ideas we obtain some interesting results. However, from the practical standpoint, due to the computational expenses associated with these schemes and the questionable benefits in applying them to perform nonlinear diffusion filtering when using long timesteps, we conclude that AOS schemes are simple and efficient compared to these alternatives.We then examine the potential utility of using multiple timestep methods combined with AOS schemes, as means to expedite the diffusion process. These methods were developed for molecular dynamics applications and are used efficiently in biomolecular simulations. The idea is to split the forces exerted on atoms into different classes according to their behavior in time, and assign longer timesteps to nonlocal, slowly-varying forces such as the Coulomb and van der Waals interactions, whereas the local forces like bond and angle are treated with smaller timesteps. Multiple timestep integrators can be derived from the Trotter factorization, a decomposition that bears a strong resemblance to a Strang splitting. Both formulations decompose the time propagator into trilateral products to construct multiplicative operator splittings which are second order in time, with the possibility of extending the factorization to higher order expansions. While a Strang splitting is a decomposition across spatial dimensions, where each dimension is subsequently treated with a fractional step, the multiple timestep method is a decomposition across scales. Thus, multiple timestep methods are a realization of the multiplicative operator splitting idea. For certain nonlinear diffusion coefficients with favorable properties, we show that a simple multiple timestep method can improve the diffusion process.     

Formulation of Kinetic Energy Preserving Conservative Schemes for Gas Dynamics and Direct Numerical Simulation of One-Dimensional Viscous Compressible Flow in a Shock Tube Using Entropy and Kinetic Energy Preserving Schemes

This paper follows up on the author’s recent paper “The Construction of Discretely Conservative Finite Volume Schemes that also Globally Conserve Energy or Enthalpy”. In the case of the gas dynamics equations the previous formulation leads to an entropy preserving (EP) scheme. It is shown in the present paper that it is also possible to construct the flux of a conservative finite volume scheme to produce a kinetic energy preserving (KEP) scheme which exactly satisfies the global conservation law for kinetic energy. A proof is presented for three dimensional discretization on arbitrary grids. Both the EP and KEP schemes have been applied to the direct numerical simulation of one-dimensional viscous flow in a shock tube. The computations verify that both schemes can be used to simulate flows with shock waves and contact discontinuities without the introduction of any artificial diffusion. The KEP scheme performed better in the tests.     

Shock Capturing Artificial Dissipation for High-Order Finite Difference Schemes

We study 2nd-, 4th-, 6th- and 8th-order accurate finite difference schemes approximating systems of conservation laws. Our goal is to utilize the high order of accuracy of the schemes for approximating complicated flow structures and add suitable diffusion operators to capture shocks. We choose appropriate viscosity terms and prove non-linear entropy stability. In the scalar case, entropy stability enables us to prove convergence to the unique entropy solution. Moreover, a limiter function that localizes the effect of the dissipation around discontinuities is derived. The resulting scheme is entropy stable for systems, and also converges to the entropy solution in the scalar case. We present a number of numerical experiments in order to demonstrate the robustness and accuracy of our scheme. The set of examples consists of a moving shock solution to the Burgers’ equation, a solution to the Euler equations that consists of a rarefaction and two contact discontinuities and a shock/entropy wave solution to the Euler equations (Shu’s test problem). Furthermore, we use the limited scheme to compute the solution to the linear advection equation and demonstrate that the limiter quickly vanishes for smooth flows and design/high-order of accuracy is retained. The numerical results in all experiments were very good. We observe a remarkable gain in accuracy when the order of the scheme is increased.     

A flux-coordinate-splitting technique for flows with shocks and contact discontinuities

The two-dimensional gasdynamic equations are solved everywhere in the flow field except in regions surrounding the contact discontinuites. A flux-vector-splitting (FVS) technique is applied to the Euler equations so that the directions of propagation of the signals and hence the shocks in the flow can be correctly captured. The split flux equations are solved using conventional second-order-accurate finite difference methods. In the regions surrounding the contact discontinuities, the gasdynamic equations are split into a set of one-dimensional equations. These are transformed in such a way that the density does not appear explicitly in the spatial derivatives of the resultant equations, which are of the Langrangian form. The equations are then solved using second-order-accurate finite difference schemes and numerical smearing of the contact discontinuities is avoided because the dependent variables are continuous across the discontinuities. Consequently, both shocks and contact discontinuities in a two-dimensional gasdynamic flow are accurately resolved. This flux-coordinate-splitting technique is used to calculate the gasdynamic flow in a shock tube, a converging cylindrical shock and the mixing of two supersonic streams. The results are compared with exact solutions and with those deduced from proven numerical techniques. Good correlations are obtained, especially in the sharp definition of contact discontinuities. Therefore, the proposed coordinate-splitting technique improves the resolution of contact discontinuities without affecting the overall calculations of the flow field. In view of this, the coordinate-splitting technique can also be used with other shock capturing techniques besides FVS to achieve the same results.     

Implementation of the CIP algorithm to magnetohydrodynamic simulations

An implementation of the Constrained Interpolation Profile (CIP) algorithm to magnetohydrodynamic (MHD) simulations is presented. First we transform the original momentum and magnetic induction equations to unfamiliar forms by introducing Elsässer variables [W.M. Elsässer, The hydromagnetic equations, Phys. Rev. (1950)]. In this formulation, while the compressional and pressure gradient terms remain as non-advective terms, the advective and magnetic stress terms are expressed in the form of an advection equation, which enables us to use the CIP algorithm. We have examined some 1D test problems using the code based on this formula. Linear Alfvén wave propagation tests reveal that the developed code is capable of solving any Alfvén wave propagation with only small numerical diffusion and phase errors up to k?h=2.5 (where ?h is the grid spacing). A shock tube test shows good agreement with a previous result with less numerical oscillation at the shock front and the contact discontinuity which are captured within a few grid points. Extension of the 1D code to the multi-dimensional case is straightforward. We have calculated the 3D nonlinear evolution of the Kelvin-Helmholtz instability (KHI) and compared the result with our previous study. We find that our new MHD code is capable of following the 3D turbulence excited by the KHI while retaining the solenoidal property of the magnetic field.     

A New Gradient Fidelity Term for Avoiding Staircasing Effect

Image denoising with some second order nonlinear PDEs often leads to a staircasing effect, which may produce undesirable blocky image. In this paper, we present a new gradient fidelity term and couple it with these PDEs to solve the problem. At first, we smooth the normal vector fields (i.e., the gradient fields) of the noisy image by total variation (TV) minimization and make the gradient of desirable image close to the smoothed normals, which is the idea of our gradient fidelity term. Then, we introduce the Euler-Lagrange equation of the gradient fidelity term into nonlinear diffusion PDEs for noise and staircasing removal. To speed up the computation of the vectorial TV minimization, the dual approach proposed by Bresson and Chan is employed. Some numerical experiments demonstrate that our gradient fidelity term can help to avoid the staircasing effect effectively, while preserving sharp discontinuities in images.     

High Accuracy Compact Schemes and Gibbs' Phenomenon

Compact difference schemes have been investigated for their ability to capture discontinuities. A new proposed scheme (Sengupta, Ganerwal and De (2003). J. Comp. Phys. 192(2), 677.) is compared with another from the literature Zhong (1998). J. Comp. Phys. 144, 622 that was developed for hypersonic transitional flows for their property related to spectral resolution and numerical stability. Solution of the linear convection equation is obtained that requires capturing discontinuities. We have also studied the performance of the new scheme in capturing discontinuous solution for the Burgers equation. A very simple but an effective method is proposed here in early diagnosis for evanescent discontinuities. At the discontinuity, we switch to a third order one-sided stencil, thereby retaining the high accuracy of solution. This produces solution with vastly reduced Gibbs' phenomenon of the solution. The essential causes behind Gibbs' phenomenon is also explained.     

Weak discontinuity annihilation in solidification modelling

Discontinuous behaviour provides substantial obstacles to the efficient application of mesh based numerical techniques. Accounting for strong discontinuities is presently of particular interest to the finite element research community with for example the development of cohesive and enriched elements to cater for material separation. Although strong discontinuities are of importance, of equal if not of greater interest and the focus in this paper, are weak discontinuities, which are present at any material change. A recent innovation for accounting for weak discontinuities has been the discovery of non-physical variables which are founded and defined using transport equations.This paper is concerned with the application of the non-physical approach to solidification modelling in the presence of more than one material discontinuity. A typical feature of the enthalpy-temperature response in solidification is discontinuities at phase transition temperatures as a consequence of phase change and latent heat release. In these circumstances, depending on the conditions that prevail, an element in a finite element mesh can have more than one discontinuity present.Presented in the paper is a methodology that can cater for multiple discontinuities. The non-physical approach permits the precise removal of weak discontinuities arising in the governing transport equations. In order to facilitate the application of the approach the finite element equations are presented in the form of weighted transport equations. The method utilises a non-physical form of enthalpy that possesses a remarkable source distribution like property at a discontinuity. It is demonstrated in the paper that it is through this property that multiple discontinuities can be exactly removed from an element so facilitating the use of continuous approximations.The new methodology is applied to a range of simple problems to provide an in-depth treatment and for ease of understanding to demonstrate the methods remarkable accuracy and stability.     

Diffusion-Inspired Shrinkage Functions and Stability Results for Wavelet Denoising

We study the connections between discrete one-dimensional schemes for nonlinear diffusion and shift-invariant Haar wavelet shrinkage. We show that one step of a (stabilised) explicit discretisation of nonlinear diffusion can be expressed in terms of wavelet shrinkage on a single spatial level. This equivalence allows a fruitful exchange of ideas between the two fields. In this paper we derive new wavelet shrinkage functions from existing diffusivity functions, and identify some previously used shrinkage functions as corresponding to well known diffusivities. We demonstrate experimentally that some of the diffusion-inspired shrinkage functions are among the best for translation-invariant multiscale wavelet denoising. Moreover, by transferring stability notions from diffusion filtering to wavelet shrinkage, we derive conditions on the shrinkage function that ensure that shift invariant single-level Haar wavelet shrinkage is maximum–minimum stable, monotonicity preserving, and variation diminishing.First online version published in June, 2005     

Vertex-Centered Linearity-Preserving Schemes for Nonlinear Parabolic Problems on Polygonal Grids

On arbitrary polygonal grids, a family of vertex-centered finite volume schemes are suggested for the numerical solution of the strongly nonlinear parabolic equations arising in radiation hydrodynamics and magnetohydrodynamics. We define the primary unknowns at the cell vertices and derive the schemes along the linearity-preserving approach. Since we adopt the same cell-centered diffusion coefficients as those in most existing finite volume schemes, it is required to introduce some auxiliary unknowns at the cell centers in the case of nonlinear diffusion coefficients. A second-order positivity-preserving algorithm is then suggested to interpolate these auxiliary unknowns via the primary ones. All the schemes lead to symmetric and positive definite linear systems and their stability can be rigorously analyzed under some standard and weak geometry assumptions. More interesting is that these vertex-centered schemes do not have the so-called numerical heat-barrier issue suffered by many existing cell-centered or hybrid schemes (Lipnikov et al. in J Comput Phys 305:111–126, 2016). Numerical experiments are also presented to show the efficiency and robustness of the schemes in simulating nonlinear parabolic problems.     

A High-Order Accurate Method for Frequency Domain Maxwell Equations with Discontinuous Coefficients

Maxwell equations contain a dielectric coefficient ɛ that describes the particular media. For homogeneous materials the dielectric coefficient is constant. There is a jump in this coefficient across the interface between differing media. This discontinuity can significantly reduce the order of accuracy of the numerical scheme. We present an analysis and implementation of a fourth order accurate algorithm for the solution of Maxwell equations with an interface between two media and so the dielectric coefficient is discontinuous. We approximate the discontinuous function by a continuous one either locally or in the entire domain. We study the one-dimensional system in frequency space. We only consider schemes that can be implemented for multidimensional problems both in the frequency and time domains.     

Locally Conservative Continuous Galerkin FEM for Pressure Equation in Two-Phase Flow Model in Subsurfaces

A typical two-phase model for subsurface flow couples the Darcy equation for pressure and a transport equation for saturation in a nonlinear manner. In this paper, we study a combined method consisting of continuous Galerkin finite element methods (CGFEMs) followed by a post-processing technique for Darcy equation and a nodal centered finite volume method (FVM) with upwind schemes for the saturation transport equation, in which the coupled nonlinear problem is solved in the framework of operator decomposition. The post-processing technique is applied to CGFEM solutions to obtain locally conservative fluxes which ensures accuracy and robustness of the FVM solver for the saturation transport equation. We applied both upwind scheme and upwind scheme with slope limiter for FVM on triangular meshes in order to eliminate the non-physical oscillations. Various numerical examples are presented to demonstrate the performance of the overall methodology.