A spacetime smooth artificial viscosity method with wavelet noise indicator and shock collision scheme, Part 2: the 2D case
Abstract
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, spacetime 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 reactiondiffusion equations, which we call equations, whose solutions serve as spacetime smooth artificial viscosity indicators for shocks and contacts.
We develop a highorder numerical algorithm for gas dynamics in 2 which can accurately simulate the RayleighTaylor (RT) instability with KelvinHelmholtz (KH) rollup of the contact discontinuity, as well as shock collision and bounceback. 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 highfrequency 2 waveletbased noise detector together with an efficient and localized noise removal algorithm.
To test the methodology, we use a highly simplified WENObased 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 shockwall collision problem.
Contents
 1 Introduction
 2 The compressible Euler equations in 2
 3 The method in the 2 setting
 4 The anisotropic method for contact discontinuities
 5 The method for shockwall collision
 6 A waveletbased 2 noise detection and removal algorithm
 7 The numerical algorithm for the method
 8 The RayleighTaylor instability
 9 The Noh infinite strength shock problem
 10 The Sod circular explosion problem
 11 Concluding remarks
 A Numerical discretization and schemes
1 Introduction
This is the second in a twopart series of papers, in which we develop a highorder numerical algorithm, the method, to simulate compressible fluid flow with shock waves and contact discontinuities, as well as shockwall collision and bounceback. 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, multidimensional conservation laws such as the compressible Euler equations develop singularities in finitetime [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 RankineHugoniot conditions (see §2). The objective of a highorder numerical scheme is to produce a simulation which keeps the fronts sharp, while simultaneously providing highorder accuracy in smooth regions away from the front.
1.1 Numerical discretization
In §1.1 of [13], we described the tools necessary for designing highorder accurate numerical schemes in 1. In multi, similar tools are required to obtain nonoscillatory numerical schemes, but the multidimensional 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 higherorder numerical schemes and subsequently remove nonphysical 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 highorder 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 selfcontained development of the 1 method; for problems in one space dimension, the function is the solution to a forced scalar reactiondiffusion equation, and serves as a highly localized spacetime smooth indicator for the shock location. In 2, we again solve a scalar reactiondiffusion 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 timedependent 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 KelvinHelmholtz rollup without overly diffusing the mixing regions in such flows.
1.3 Shock collision and waveletbased noise removal
In part one [13], we consider the difficult problem of shockwall collision and bounceback 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 postcollision 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 timedependent artificial viscosity parameter naturally increases during shockwall collision and bounceback, allows for the addition of extra “wall viscosity” to the shock front. This suppresses postcollision noise and the wall heating error, and ensures that the solution retains highorder 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 wellknown issue. The first part of this work [13] provides a description of an efficient waveletbased noise detection and heat equationbased 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 RankineHugoniot 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 longtime evolution of contact discontinuities. A new 2 shockwall collision scheme is introduced in §5, and a waveletbased noise detection and localized heat equationbased 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 RayleighTaylor instability in §8, the Noh infinite strength shock problem in §9, the Sod circular explosion problem in §10, and a Sodtype shockwall 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
(1a)  
(1b) 
where the 4vector 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
(2) 
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
(3) 
We focus on solutions to the compressible Euler equations (1) that have a jump discontinuity across a timedependent curve
(4) 
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 RankineHugoniot 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 timedependent curves.
(5) 
In Fig.1, we see that captures both of the discontinuities, namely the shock front and the contact discontinuity, present in the solution.
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 spacetime 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
(6) 
where with , the grid spacings in the and directions, respectively, is the maximum wavespeed (3), and denotes the 2 Laplace operator. The spacetime smooth version of , which we denote by , is the solution to the scalar linear parabolic PDE
(7) 
where the function is a compression switch
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 RayleighTaylor 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.
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 spacetime smooth artificial viscosity: \cref@addtoresetequationparentequation
(8a)  
(8b)  
(8c)  
(8d)  
(8e) 
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
(9) 
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 spacetime smooth localizing function ensures that viscosity is added only at discontinuities, thereby ensuring that the solution retains highorder accuracy away from discontinuities, and, moreover, given that is a solution to a reactiondiffusion 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
longtime evolution of contact discontinuities in
the presence of RayleighTaylor (RT) instabilities which lead to KelvinHelmholtz (KH) rollup. We are particularly interested in the case that
across the contact discontinuity
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 Machnumber flows), instabilities that are generated by a small perturbation of the equilibrium interface position require a large number of timesteps to fully develop the KH rollup of the contact curve, and it is most often the case that fluid mixing (due to numerical dissipation) is present in this rollup 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 overdiffusion of the slight “hillsandvalleys” 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 ATHENA
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
(10a)  
(10b) 
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.
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 spacetime smoothing mechanism provided by the method to produce regularized versions of and , which we denote and , respectively: \cref@addtoresetequationparentequation
(11a)  
(11b) 
where is the operator defined in (6). We also define the function
(12) 
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 .
4.2 Directional artificial viscosity and the Euler system
We now consider the following Euler system for contact discontinuity evolution: \cref@addtoresetequationparentequation
(13a)  
(13b)  
(13c)  
(13d)  
(13e) 
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
(14) 
with .
We note that the artificial viscosity operator is a positive semidefinite 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 EulerLagrange 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 shockwall collision
We now present a simple extension of the 1 shockwall 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 shockwall collision and bounceback. This allows for the addition of extra “wall viscosity” during shockcollision, which results in the suppression of postcollision noise while maintaining highorder accuracy prior to collision.
The natural generalization to the twodimensional setting results in the 2 Euler scheme: \cref@addtoresetequationparentequation
(15a)  
(15b)  
(15c)  
(15d)  
(15e)  
(15f) 
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
(16) 
with the wall function defined by
Here, we assume that the shockwall 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 shockwall collision and bounceback.
6 A waveletbased 2 noise detection and removal algorithm
In this section, we extend the noise indicator algorithm presented in [13] to the twodimensional 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 equationbased 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.
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:

Zero mean:

“Quasiorthogonality” 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.
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 twodimensional setting is a plane, so a first attempt is to approximate by a plane in each of the subcells . However, this is not possible, since 3 points define a plane, whereas each of the subcells contains 4 cell center values.
Our solution is to divide into 8 regions as shown in Fig.7. Each of these regions contains precisely 3 cellcenter 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
(17) 
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
(18) 
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” highfrequency oscillation, namely a hat function (see Fig.8) centered in the domain with amplitude . The associated wavelet coefficient may then be calculated as
(19) 
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 nonzero 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 .
A localized heat equation with Dirichlet boundary conditions is then solved in each of the domains for a “denoised” solution , \cref@addtoresetequationparentequation
(20a)  
(20b)  
(20c) 
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