A space-time smooth artificial viscosity method with wavelet noise indicator and shock collision scheme, Part 2: the 2-D case

A space-time smooth artificial viscosity method with wavelet noise indicator and shock collision scheme, Part 2: the 2-D case


This is the second part to our companion paper [13]. Herein, we generalize to two space dimensions the -method developed in [14, 13] for adding localized, space-time smooth artificial viscosity to nonlinear systems of conservation laws that propagate shock waves, rarefaction waves, and contact discontinuities. For gas dynamics, the -method couples the Euler equations to scalar reaction-diffusion equations, which we call -equations, whose solutions serve as space-time smooth artificial viscosity indicators for shocks and contacts.

We develop a high-order numerical algorithm for gas dynamics in 2- which can accurately simulate the Rayleigh-Taylor (RT) instability with Kelvin-Helmholtz (KH) roll-up of the contact discontinuity, as well as shock collision and bounce-back. Solutions to our -equations not only indicate the location of the shocks and contacts, but also track the geometry of the evolving fronts. This allows us to implement both directionally isotropic and anisotropic artificial viscosity schemes, the latter adding diffusion only in directions tangential to the evolving front.

We additionally produce a novel shock collision indicator function, which naturally activates during shock collision, and then smoothly deactivates. Moreover, we implement a high-frequency 2- wavelet-based noise detector together with an efficient and localized noise removal algorithm.

To test the methodology, we use a highly simplified WENO-based discretization scheme, devoid of any approximate Riemann solvers. We provide numerical results for some classical 2- test problems, including the RT problem, the Noh problem, a circular explosion problem from the Liska & Wendroff [9] review paper, and a shock-wall collision problem.


1 Introduction

This is the second in a two-part series of papers, in which we develop a high-order numerical algorithm, the -method, to simulate compressible fluid flow with shock waves and contact discontinuities, as well as shock-wall collision and bounce-back. In the first part [13], we developed our scheme in one space dimension. This second part is devoted to the more geometric problem of simulating shocks, contacts, and collisions in two space dimensions.

Compressible fluid flow is an example of a system of nonlinear conservation laws, which in two space dimensions, can be written as

where denotes spatial coordinates on the plane, denotes time, is a vector of conserved quantities, and and denote the horizontal and vertical flux functions, respectively.

Even with smooth initial conditions, multi-dimensional conservation laws such as the compressible Euler equations develop singularities in finite-time [18, 3] and, in general, solutions consist of propagating shock waves and contact discontinuities. In two space dimensions, shocks and contacts produce curves of discontinuities, also known as fronts, which propagate according to the Rankine-Hugoniot conditions (see §2). The objective of a high-order numerical scheme is to produce a simulation which keeps the fronts sharp, while simultaneously providing high-order accuracy in smooth regions away from the front.

1.1 Numerical discretization

In §1.1 of [13], we described the tools necessary for designing high-order accurate numerical schemes in 1-. In multi-, similar tools are required to obtain non-oscillatory numerical schemes, but the multi-dimensional analogues are generally limited by mesh considerations. For structured grids (such as products of uniform 1- grids), dimensional splitting is commonly used, decomposing the problem into a sequence of 1- problems. This technique is quite successful, but stringent mesh requirements prohibits its use on complex domains. Moreover, applications to PDE outside of variants of the Euler equations may be somewhat limited. For further discussion of the limitations of dimensional splitting, we refer the reader to Crandall & Majda [4], and Jiang & Tadmor [7]. For unstructured grids, dimensional splitting is not available and alternative approaches must be employed, necessitated by the lack of multi- Riemann solvers. WENO schemes on unstructured triangular grids have been developed in Hu & Shu [5], but using simplified methods, which employ reduced characteristic decompositions, can lead to a loss of monotonicity and stability.

Algorithms that explicitly introduce diffusion provide a simple way to stabilize higher-order numerical schemes and subsequently remove non-physical oscillations near shocks. We refer the reader to the introductory sections in [14, 13] for a review of the classical artificial viscosity method [23].

In this paper, we develop a stable high-order 2- numerical scheme that does not use approximate Riemann solvers or characteristic decompositions, but instead relies upon a 2- generalization of the -method.

1.2 Localization to shocks and contacts and capturing the geometry of the front

As we noted above, part one [13] provides a self-contained development of the 1- -method; for problems in one space dimension, the function is the solution to a forced scalar reaction-diffusion equation, and serves as a highly localized space-time smooth indicator for the shock location. In 2-, we again solve a scalar reaction-diffusion equation for on the plane; in this instance, the function not only serves as a localized indicator for both shocks and contacts, but is also able to accurately represent the geometry of the evolving fronts. Having this geometry allows us to define time-dependent normal and tangent vectors to the fronts, which, in turn, enables the construction of artificial viscosity methods that add diffusion only in certain directions rather than uniformly across the mesh. As we shall explain below, this type of anisotropic artificial viscosity scheme is essential in accurately capturing Kelvin-Helmholtz roll-up without overly diffusing the mixing regions in such flows.

1.3 Shock collision and wavelet-based noise removal

In part one [13], we consider the difficult problem of shock-wall collision and bounce-back in 1-. In particular, when a shock wave collides with and reflects off of a fixed boundary, spurious oscillations develop behind the reflected shock. In addition to these post-collision oscillations, most schemes produce solutions that exhibit the phenomenon of anomalous wall heating [15]. A novel modification to the -method in 1- [13], wherein the time-dependent artificial viscosity parameter naturally increases during shock-wall collision and bounce-back, allows for the addition of extra “wall viscosity” to the shock front. This suppresses post-collision noise and the wall heating error, and ensures that the solution retains high-order accuracy away from the shock and prior to collision. In §5 of this paper, we present a generalization of the 1- shock collision scheme to the 2- setting, and apply the method to a circular explosion problem in §10.2.

The occurrence of high frequency noise in numerical solutions to PDE is a well-known issue. The first part of this work [13] provides a description of an efficient wavelet-based noise detection and heat equation-based noise removal algorithm for 1- gas dynamics simulations. Error analysis and convergence tests in [13] show that the noise detection and removal algorithm decreases the errors of numerical solutions and improves the rate of convergence. In §6 of this paper, we present the natural extension of this 1- algorithm to the 2- case, with applications to the RT problem in §8 and the Noh problem in §9.

1.4 Outline of the paper

In §2, we introduce the 2- compressible Euler equations and the Rankine-Hugoniot jump conditions relating the speed of propagation of curves of discontinuity with the jump discontinuities in the conservative variables. In §3, we present the isotropic -method, and in §4, the anisotropic -method, the latter designed for the long-time evolution of contact discontinuities. A new 2- shock-wall collision scheme is introduced in §5, and a wavelet-based noise detection and localized heat equation-based noise removal procedure is developed in §6. Details about the numerical methods implemented are provided in §7 and Appendix A. We then apply the methods to the Rayleigh-Taylor instability in §8, the Noh infinite strength shock problem in §9, the Sod circular explosion problem in §10, and a Sod-type shock-wall collision problem in §10.2

2 The compressible Euler equations in 2-

The compressible Euler equations in a 2- rectangular domain and a time interval are given in conservation law form as \cref@addtoresetequationparentequation


where the 4-vector and the flux functions and are defined as

Here, the vector denotes the two Cartesian coordinates, and are the density and energy fields, respectively, and are the velocities in the -direction and -direction, respectively, and we use the notation to denote the velocity vector field . The pressure function is defined by the equation of state


where and is the adiabatic constant. The system (1) represents the conservation of mass, momentum, and energy:

and defines a hyperbolic system, in the sense that both and have real eigenvalues; in particular, the four eigenvalues of are

and the four eigenvalues of are

The eigenvalues and correspond to sound waves, while the repeated eigenvalues , and correspond to vorticity and entropy waves. We define the maximum wave speed as


We focus on solutions to the compressible Euler equations (1) that have a jump discontinuity across a time-dependent curve


where for each index , represents either a shock front or a contact discontinuity. In the case that is a closed curve, we define to be the value of inside of and to be the value of outside of . We then set

Let and be the unit normal and tangent vectors to , respectively. The Rankine-Hugoniot conditions relate the speed of propagation of the curve of discontinuity with the jump discontinuity in the variables via the relation

If is a shock,

while if is a contact discontinuity, then

For the problems we consider in this paper, for both shocks and contacts.

3 The -method in the 2- setting

3.1 Smoothly localizing the curves of discontinuity and tracking the geometry

As noted above, the jump discontinuities of the density function describe the location of both the shock and contact fronts; as such, it is natural to use as an indicator for these time-dependent curves.1 Consequently, we track discontinuities by considering the quantity


In Fig.1, we see that captures both of the discontinuities, namely the shock front and the contact discontinuity, present in the solution.

Figure 1: The function for the Sod explosion problem at time .

Although the function is able to localize artificial viscosity to a curve of discontinuity, its use in artificial viscosity operators often serves to produce an oscillatory solution [14, 13]. This is due to the rough nature of the localizing function . Consequently, we first produce a space-time smooth variant of through the use of the -method. We define the natural 2- generalization of the 1- -equation as follows: we first define the operator as


where with , the grid spacings in the and directions, respectively, is the maximum wave-speed (3), and denotes the 2- Laplace operator. The space-time smooth version of , which we denote by , is the solution to the scalar linear parabolic PDE


where the function is a compression switch2 that ensures that vanishes in regions of expansion, where there are no discontinuities in the solution. The parameters and in (6) control the support and smoothness of the solution , respectively [14, 13].

The function provides not only the location of the shock and contact fronts, but also a good approximation to the geometry of the front. Specifically, the function is sufficiently localized so as to provide the shape of the evolving front. In Fig.2, we show results of the evolution of a contact discontinuity associated to the Rayleigh-Taylor instability (which is discussed in §8). As can be seen, and hence the function track the material interface accurately but the function is very rough, particularly in directions tangential to the contact front. The function , by contrast, is smooth and exhibits no oscillatory behavior, but is still highly localized to the contact discontinuity.

Figure 2: Comparison of and for the Rayleigh-Taylor instability at various times. The top row (a,b,c) are surface plots of , and the bottom row (d,e,f) are surface plots of the corresponding .

3.2 The isotropic -method and the Euler- system

We begin with a very natural generalization of the 1- -method to the 2- setting. In particular, we first consider the 2- Euler- system with isotropic space-time smooth artificial viscosity: \cref@addtoresetequationparentequation


where the pressure is given by the equation of state (2), is the operator defined in (6), is the forcing function defined in (5), and the artificial viscosity parameters are given by


with a constant.

We refer to the system (8) as the isotropic -method, because artificial viscosity is added uniformly in all directions (although clearly still localized to the fronts via the use of ). We remark that the particular form of artificial viscosity used in the momentum equations (8b) and (8c) as well as the energy equation (8d) ensure that total energy remains conserved, and that continues to evolve as the total energy function. For simplicity, suppose that periodic boundary conditions are enforced. On the one hand, integrating the energy equation (8d) yields . On the other hand, multiplying (8b) and (8c) by and , respectively, integrating over the domain , summing the resulting quantities and utilizing the conservation of mass equation (8a), the energy equation (8d), and the equation of state (2) yields

This shows that the velocity and the pressure adjust accordingly to maintain the relation (2), and that the Euler- system conserves the total energy.

The space-time smooth localizing function ensures that viscosity is added only at discontinuities, thereby ensuring that the solution retains high-order accuracy away from discontinuities, and, moreover, given that is a solution to a reaction-diffusion equation, it is smooth in both space and time. See [14] for the analysis of the solutions .

4 The anisotropic -method for contact discontinuities

We next consider the anisotropic -method, specifically designed for the long-time evolution of contact discontinuities in the presence of Rayleigh-Taylor (RT) instabilities which lead to Kelvin-Helmholtz (KH) roll-up. We are particularly interested in the case that across the contact discontinuity3, for which we also have that and .

Conservation of mass can be written as

and near an evolving front with tangent and normal vectors given by and , respectively, the divergence of the velocity vector field is given by

Across a contact curve, remains smooth, while can become extremely oscillatory due to a combination of the discontinuity of across the contact curve together with interpolation error of the contact curve onto a fixed grid (particularly in specifying the initial data). For simulations that require a great deal of time steps, it is important to add artificial viscosity in the tangential directions, but not in the normal directions.

In classical RT problems (and particularly for low Mach-number flows), instabilities that are generated by a small perturbation of the equilibrium interface position require a large number of time-steps to fully develop the KH roll-up of the contact curve, and it is most often the case that fluid mixing (due to numerical dissipation) is present in this roll-up region, so that the amplitude of the density and the gradient of the density can be significantly smaller than their maximum values.

The objective of our anisotropic -method is to add diffusion only in directions that are tangent to the contact discontinuity, while adding no artificial diffusion in directions that are normal to the contact curve. In doing so, we can maintain a very sharp interface, and prevent over-diffusion of the slight “hills-and-valleys” which arise in the KH mixing zones. Moreover, when generating the RT instability of a small interface perturbation on a uniform rectangular grid (rather than using a velocity perturbation as done by ATHENA4), spurious tangential spikes can form in the velocity fields and along the contact curve; these spikes, in turn, generate small-scale numerical KH structures which contaminate the solution (this is discussed further in §8.1). Consequently, it is necessary to remove these spikes while maintaining a sharp interface; this may be accomplished through the use of anisotropic diffusion.

4.1 Calculating the normal and tangential directions to

The first task is to accurately compute a good approximation to the tangent vectors to any curve of discontinuity , defined in (4). For the problems we consider here, this may be accomplished by setting \cref@addtoresetequationparentequation


where and denote the unit vectors in the and directions, respectively. In Fig.3, we provide vector plots of calculated using (10), as well as surface plots of each component of , for the specific case of the RT instability.

(a) Vector plot of at
(b) Vector plot of at
(c) Surface plot of at
(d) Surface plot of at
Figure 3: Calculation of the tangent vector for the RT instability.

As shown in Fig.3, the functions and suffer from the same issue affecting the localizing function ; namely, a lack of smoothness in both space and time. Consequently, we utilize the space-time smoothing mechanism provided by the -method to produce regularized versions of and , which we denote and , respectively: \cref@addtoresetequationparentequation


where is the operator defined in (6). We also define the function


which can be used to produce a “normalized” tangent vector .

In Fig.4, we provide vector plots of the vector , as well as surface plots of the components and . These should be contrasted with the corresponding figures in Fig.3; we see that is much smoother than , while remaining localized to the discontinuity .

(a) Vector plot of at
(b) Vector plot of at
(c) Surface plot of at
(d) Surface plot of at
Figure 4: Calculation of the smoothed tangent vector for the RT instability.

4.2 Directional artificial viscosity and the Euler- system

We now consider the following Euler- system for contact discontinuity evolution: \cref@addtoresetequationparentequation


where is given by (5), and by (10), by (6), and we utilize the Einstein summation convention, where a repeated free index in the same term implies summation over all values of that index. We note that and , and we use the two notations interchangeably. The artificial viscosity parameter is defined by


with .

We note that the artificial viscosity operator is a positive semi-definite operator. The proof is as follows: taking the scalar product of with and integrating over yields

for some . Here, denotes the smoothed tangential derivative operator. Thus, the anisotropic artificial viscosity operator is obtained as the Euler-Lagrange extremum associated to the function . Just as in the case of isotropic artificial viscosity, our solutions to the Euler- system preserve total energy (see §3.2).

5 The -method for shock-wall collision

We now present a simple extension of the 1- shock-wall collision scheme (see §3 [13]). Recall that the main novelty of the 1- shock wall collision scheme is the use of a wall function which naturally activates during shock-wall collision and bounce-back. This allows for the addition of extra “wall viscosity” during shock-collision, which results in the suppression of post-collision noise while maintaining high-order accuracy prior to collision.

The natural generalization to the two-dimensional setting results in the 2- Euler-- scheme: \cref@addtoresetequationparentequation


where is the operator defined in (6), is the forcing function defined in (5), and is a compression switch. The artificial viscosity parameters are given by


with the wall function defined by

Here, we assume that the shock-wall collision occurs at the bottom boundary . A Neumann boundary condition for is enforced at the bottom boundary . As in the 1- case, this results in the smooth growth in time of the amplitude of as the shock approaches the wall, and allows for the addition of “wall viscosity” during shock-wall collision and bounce-back.

6 A wavelet-based 2- noise detection and removal algorithm

In this section, we extend the noise indicator algorithm presented in [13] to the two-dimensional setting, using the same ideas as in the 1- case. Again, we construct a family of wavelets and obtain a set of wavelet coefficients , found by calculating the inner product of the highest frequency wavelets with the noisy function. These wavelet coefficients will indicate the location of the noise, and we employ a localized heat equation-based solver to remove this noise.

6.1 Construction of the highest frequency wavelets

We first discretize our grid by assuming we have cells in the -direction, and cells in the -direction. Label the cell centers by , where and are given by

with and . We group the cells into 3x3 blocks of 9 cells each, and then define the highest frequency wavelet with support over the domain spanned by the cell centers in each of these blocks, as shown in Fig.5.

left boundary

right boundary

top boundary

bottom boundary

Figure 5: Grid setup for the construction of the wavelets in 2-. The dashed black curve denotes the boundary , while the solid black lines indicate cell edges. The black dots indicate cell centers, and the domains bounded by red lines indicate the support of each of the highest frequency wavelets.

This yields a set of highest frequency wavelets. We denote these wavelets by , for and , with each is supported in the rectangular domain .

We now have to fix a form for the 2- wavelet. The two key properties that are required of the wavelet family are:

  1. Zero mean:

  2. “Quasi-orthogonality” of the form:

    so that each of the highest frequency wavelets is orthogonal to every other wavelet of the same frequency.

We recall from [13] that each member of the 1- wavelet family is oscillatory and orthogonal to all linear functions. Consequently, we design the 2- wavelet family to also have these two properties. For simplicity, assume that a highest frequency wavelet is supported in the domain . The highest frequency wavelets are then obtained by translation of to the domain . Our wavelets take the form shown in Fig.6.

(a) the highest frequency wavelet
(b) cross sectional view of along
Figure 6: The highest frequency wavelet in 2-. We assume that the support of is .

The exact formula for the highest frequency wavelet that is supported in the domain is

where and are given by

Here, the values of and are chosen so that satisfies and the zero mean condition , and can be calculated as

The highest frequency wavelets are then obtained by translation of to the domain .

Now, given a function defined at the cell centers , we wish to calculate the inner product of with each of the highest frequency wavelets. The first step is to approximate the function over , the support of . In 1-, a function is approximated over the interval , the support of a 1- highest frequency wavelet, by linearly interpolating between the cell center values (see §4 in [13]). The analogue of a line in the two-dimensional setting is a plane, so a first attempt is to approximate by a plane in each of the sub-cells . However, this is not possible, since 3 points define a plane, whereas each of the sub-cells contains 4 cell center values.

Figure 7: Dividing into 8 regions.

Our solution is to divide into 8 regions as shown in Fig.7. Each of these regions contains precisely 3 cell-center values, and these 3 values define a plane. Thus, we can approximate by a plane in each of these 8 regions, and define a piecewise linear function such that . We may then approximate the -th wavelet coefficient as


where . One can verify that the inner product of with an arbitrary plane is identically zero i.e. if for some constants and , then for every and . This is analogous to the result in 1- that each of the highest frequency wavelets is orthogonal to linear functions.

6.2 Algorithm

Given a noisy function , we now present an algorithm for first detecting and then subsequently removing the noise. The algorithm is almost identical to that presented for the 1- case (see §4 in [13]).

Noise detection in the presence of discontinuities

We first calculate, using formula (17), the wavelet coefficients, each associated with one of the highest frequency wavelets. The coefficients that are largest in magnitude should indicate the location of the noise. However, as in the 1- case, it is possible that the wavelet coefficients that are largest in magnitude are actually indicating the location of the curve of discontinuity . Since the -method is taking care of artificial diffusion in this region, one needs to manually “turn off” the noise detection in a small region surrounding the shock curve . There are numerous ways to do this, but we employ one of the simplest methods, namely to turn off the noise detection if


where is some value between 0 and 1. A typical range of values for is .

With noise detection deactivated in the region surrounding , the wavelet coefficients that are largest in magnitude now indicate the location of the noise. We now “turn on” a noise detector function in the domain if the associated wavelet coefficient satisfies . The constant is the wavelet coefficient obtained from a “typical” high-frequency oscillation, namely a hat function (see Fig.8) centered in the domain with amplitude . The associated wavelet coefficient may then be calculated as

(a) reference oscillation: a hat function
(b) cross sectional view along
Figure 8: The 2- reference oscillation: a hat function with amplitude .

Noise removal with a localized heat equation

The noise removal process in the 2- case is identical to that in 1-. We first construct the domain , given by the union of all domains such that the noise detector function is non-zero in . We write as the union of its connected subsets , and then define the domains as the domain extended by one cell in each outward in each direction (see Fig.9). For example, if , then .

Figure 9: Construction of the domains and . The squares bounded by the red lines indicate the support of each of the highest frequency wavelets . The shaded green regions indicate where the noise indicator algorithm detects noise, and represent the domains . The hatched regions indicate the extension of each to . The domains are the regions bounded by the solid black lines.

A localized heat equation with Dirichlet boundary conditions is then solved in each of the domains for a “de-noised” solution , \cref@addtoresetequationparentequation


while for and . The time is a fictitious time, introduced solely for the diffusion mechanism, while is a small constant, which we refer to as the noise removal viscosity. Equation (20b) is the initial condition, and (20c) is a Dirichlet boundary condition ensuring continuity of over the domain .

However, as in the 1- case, it is not necessary to explicitly construct the domains . Instead, one can use the noise detector function and solve a modified heat equation with Dirichlet boundary conditions, given by \cref@addtoresetequationparentequation