An adaptive multiresolution discontinuous Galerkin method with artificial viscosity for scalar hyperbolic conservation laws in multidimensions

An adaptive multiresolution discontinuous Galerkin method with artificial viscosity for scalar hyperbolic conservation laws in multidimensions

Juntao Huang Department of Mathematics, Michigan State University, East Lansing, MI 48824, USA. E-mail:    Yingda Cheng Department of Mathematics, Department of Computational Mathematics, Science and Engineering, Michigan State University, East Lansing, MI 48824, USA. E-mail: Research is supported by NSF grants DMS-1453661 and DMS-1720023.

In this paper, we develop an adaptive multiresolution discontinuous Galerkin (DG) scheme for scalar hyperbolic conservation laws in multidimensions. Compared with previous work for linear hyperbolic equations [25, 26], a class of interpolatory multiwavelets are applied to efficiently compute the nonlinear integrals over elements and edges in DG schemes. The resulting algorithm, therefore can achieve similar computational complexity as the sparse grid DG method for smooth solutions. Theoretical and numerical studies are performed taking into consideration of accuracy and stability with regard to the choice of the interpolatory multiwavelets. Artificial viscosity is added to capture the shock and only acts on the leaf elements taking advantages of the multiresolution representation. Adaptivity is realized by auto error thresholding based on hierarchical surplus. Accuracy and robustness are demonstrated by several numerical tests.

Keywords: Compressible Euler equations; positivity-preserving; chemical reactions; production-destruction equations; finite difference

1 Introduction

In this paper, we develop an adaptive multiresolution discontinuous Galerkin (DG) method for scalar nonlinear conservation laws in multidimensional case:


with appropriate initial and boundary conditions. Here , is the unknown function, and is the physical flux. We assume in the paper, but the discussion can be easily generalized to arbitrary box-shaped domains.

The DG method is a class of finite element methods using discontinuous approximation space for the numerical solutions and the test functions. The Runge-Kutta DG scheme for hyperbolic equations became very popular due to its provable stability and convergence, excellent conservation properties and accommodation for adaptivity and parallel implementations. We refer readers to the review papers [16, 14] for details. To adapt the degrees of freedom according to the local behavior of the numerical solution, many kinds of a posteriori error estimates have been designed for the DG schemes for hyperbolic equations, see e.g. [7, 1, 45, 29, 32, 30]. On the other hand, by using multiresolution analysis (MRA), automatic adaptivity can be achieved and no additional a posteriori error indicator is needed. Such ideas have been used to accelerate the computations for conservation laws under finite difference or finite volume frameworks [28, 8, 18, 4, 17, 12] and were used as trouble cell indicators for DG methods [50]. In recent years, there have been interests in developing adaptive multiresolution DG schemes [11, 5]. In particular, multiresolution-based adaptive DG schemes for solving one dimensional scalar conservation laws were proposed by Müller et al. in [33] and further extended to multidimensional cases [22], compressible flows [35, 21] and shallow water equations [20, 36].

Another idea to utilize the computational advantages of the MRA framework is called the sparse grid method [10], which is a well-known tool to compute high-dimensional PDEs and stochastic differential equations. Based on the attractive features of DG methods for solving convection-dominated problems, in recent years, we initiated a line of research developing the (adaptive) sparse grid DG methods, including the work for elliptic equations [51], transport equations [25], reaction-diffusion equations [40] and Vlasov-Maxwell equations [48]. For smooth solutions, the schemes we constructed can successfully reduce the number of degrees of freedom (DoF) of unknown from to for -dimensional problem, where is the uniform mesh size in each dimension. Stability and conservation of standard DG methods can be maintained. Errors are only slightly deteriorated for smooth solutions. Adaptivity can be incorporated naturally to treat solutions with less smoothness or local structures.

However, the main bottleneck for the sparse grid DG scheme developed so far is that, it is mainly for “linear” equations. Here, “linear” refers to either linear variable coefficient equations with given coefficients or coefficients that have some specified dependence on the unknowns, e.g. Vlasov systems through self-consistent field. There remain significant challenges to extend the methods to truly nonlinear problems in an efficient manner. For example, previous work in the literature on adaptive multiresolution DG schemes resort to the finest scale for the actual time evolution for the nonlinear terms. Therefore, the computational cost is proportional to the number of cells on the finest level, i.e. operations, and the reduced DoF in the solution representation is not realized in the actual computation. For nonlinear equations, there is only limited literature on collocation or finite difference based sparse grid methods [23], and the order of accuracy of the schemes is low. Sparse grid combination methods work for nonlinear problems, but they are less flexible in terms of adaptivity [39].

The objective of this work is to develop an adaptive multiresolution DG method to solve the nonlinear equation (1) efficiently. We would like to construct a scheme that can recover the computational complexity of the sparse grid method for smooth solutions. To represent a nonlinear function, it is natural to consider interpolation or collocation methods [31, 2]. This is achieved by using sparse grid collocation methods introduced in [49], which gives a framework to design adaptive sparse grid collocation onto arbitrary high order piecewise polynomial spaces. We approximate the nonlinear terms in the semi-discrete DG scheme for (1) by a linear combination of collocation basis functions up to required order of accuracy. We analyze the accuracy of the DG scheme with interpolation following the approach in [13, 34]. We compare different ways of sparse grid collocation methods and find that there exists some weak instability for Lagrange interpolation when solving (1). This motivates us to apply Hermite interpolation, which is more stable than the Lagrange interpolation [27]. With the Hermite multiwavelet interpolations, we can obtain satisfactory numerical results which coincide with our local truncation error analysis.

Another challenge we address in the paper is how to capture the shock and entropy solutions to (1). There are two approaches in the literature. The first one is to apply limiters to control spurious oscillations and at the same time maintain accuracy in smooth regions, e.g., the minmod-type limiter [15], the moment-based limiter [9] and WENO limiter [43]. However, it is quite difficult to impose limiters in the sparse grid DG methods, due to the global feature of the basis functions. Also a preliminary calculation from us shows that the piecewise constant sparse grid DG method in multidimensions is not monotone. This motivates us to use the second approach, which is to add artificial viscosity, see e.g. [6, 42, 24, 37, 38]. The idea is to add a diffusion term in the equation where the diffusion coefficient vanishes in the smooth region and becomes non-zero near the shock. This can be achieved by techniques such as entropy production [24] or local smoothness indicator [42]. We add an artificial viscosity term following the approach in [6]. Based on the estimate of the magnitudes of coefficients of hierarchical basis functions in [25, 26], we propose a smoothness indicator, which is built upon the inherent MRA and can automatically pick out the discontinuous regions. To improve the computational efficiency of our scheme, the implicit-explicit (IMEX) Runge-Kutta time integration is applied, where the nonlinear convection term is treated explicitly and the linear diffusion term is computed implicitly.

The rest of this paper is organized as follows. In Section 2, we review MRA associated with two sets of basis functions, i.e., the Alpert’s multiwavelets [3] and the interpolatory multiwavelets [49]. The adaptive multiresolution DG scheme is constructed in Section 3 using both sets of multiwavelets. The numerical performance is validated by linear advection equations, Burgers’ equations and KPP problems in Section 4. We conclude the paper in Section 5. The appendix collects the explicit formulas of the interpolatory multiwavelets used in this paper.

2 MRA and multiwavelets

In this section, we review MRA associated with piecewise polynomial space. We will start with Alpert’s multiwavelets [3] and then review the interpolatory multiwavelets [49].

2.1 Alpert’s multiwavelets

We first review MRA achieved by Alpert’s basis functions in one dimension [3]. We define a set of nested grids, where the -th level grid consists of uniform cells

for For notational convenience, we also denote The usual piecewise polynomial space of degree at most on the -th level grid for is denoted by


Then, we have the nested structure

We can now define the multiwavelet subspace , as the orthogonal complement of in with respect to the inner product on , i.e.,

For notational convenience, we let , which is the standard polynomial space of degree up to on . Therefore, we have .

Now we define a set of orthonormal basis associated with the space . The case of mesh level is trivial. We use the normalized shifted Legendre polynomials in and denote the basis by for . When , the orthonormal bases in are presented in [3] and denoted by

The construction follows a repeated Gram-Schmidt process and the explicit expression of the multiwavelet basis functions are provided in Table 1 in [3]. Note that such multiwavelet bases retain the orthonormal property of wavelet bases for different mesh levels, i.e.,


and the support of is in .

Multidimensional case when follows from a tensor-product approach. First we recall some basic notations. For a multi-index , where denotes the set of nonnegative integers, the and norms are defined as

The component-wise arithmetic operations and relational operations are defined as

By making use of the multi-index notation, we denote by the mesh level in a multivariate sense. We define the tensor-product mesh grid and the corresponding mesh size Based on the grid , we denote as an elementary cell, and

as the tensor-product piecewise polynomial space, where represents the collection of polynomials of degree up to in each dimension on cell . If we use equal mesh refinement of size in each coordinate direction, the grid and space will be denoted by and , respectively.

Based on a tensor-product construction, the multidimensional increment space can be defined as

Therefore, the standard tensor-product piecewise polynomial space on can be written as


while the sparse grid approximation space in [51] is


The dimension of scales as [51], which is significantly less than that of with exponential dependence on . The approximation results for are discussed in [51, 25], which has a stronger smoothness requirement than the traditional space. In this paper, we will not require the numerical solution to be in , but rather in and to be chosen adaptively similar to [26].

Finally, we define the basis functions in multidimensions as


for , and . The orthonormality of the bases can be established by (3).

2.2 Interpolatory multiwavelets

Alpert’s multiwavelets and the space are constructed so that they correspond to the difference of the projection on adjacent levels. The idea of the sparse grid collocation basis proposed in [49] is to switch the operator to be interpolation on nested grids. Below, we will outline the construction. Denote the set of interpolation points in the interval at mesh level 0 by . Here, the number of points in is . Then the interpolation points at mesh level , can be obtained correspondingly as

We require the points to be nested, i.e.


to save computational cost. This can be achieved by requiring , and then one can deduce (7) easily.

Given the nodes, we define the basis functions on the -th level grid as Lagrange or Hermite interpolation polynomials of degree which satisfy the property:

for and . It is easy to see that The constants will be specified later on in the paper. With the basis function at mesh level 0, we can define basis function at mesh level :

which is a complete basis set for

Next, we introduce the hierarchical representations. Define and for , then we have the decomposition

Denote the points in by . Then the points in for can be represented by

For notational convenience, we let The increment function space for is introduced as a function space that satisfies


and is defined through the multiwavelets that satisfies

for and . Then is given by

where . For completeness, we list the basis functions used in this paper in the appendix.

The construction above has close connection with interpolation operators. For a given function , we define as the standard Hermite interpolation on and have the representation

Clearly, The algorithm converting between the point values and the derivatives to hierarchical coefficients is given in [49], and by a standard argument in fast wavelet transform, can be performed in flops.

The multidimensional construction follows similar lines as in Section 2.1. We let


while the sparse grid approximation space is

Note that the construction by Alpert’s multiwavelet and the interpolatory multiwavelet gives the same sparse grid space. Finally, the interpolation operator in multidimension :

where the multidimensional basis functions are defined in the same approach as (6) by tensor products. If the space is switched from to some subset of e.g. the sparse grid space or some other subset of that is dynamically chosen, the interpolation operator can be defined accordingly, taking only multiwavelet basis functions that belong to that space.

3 Adaptive multiresolution DG evolution algorithm

In this section, we will describe the adaptive multiresolution DG scheme for (1). We will first introduce the DG scheme with multiresolution interpolation. The accuracy requirement for the interpolation operator is studied by local truncation error analysis. We then describe the adaptive strategy. Finally, the artificial viscosity is introduced based on the estimate of the coefficients of the hierarchical basis functions.

3.1 DG scheme with multiresolution interpolation

First, we review some basis notations about meshes. Let be the maximum mesh level and be the collection of all elementary cell , , . Define be the union of all the interfaces for all the elements in Here, for simplicity, we formulate the scheme with periodic boundary conditions, while we keep in mind other boundary conditions can be treated in the DG framework as well.

The semi-discrete DG scheme for the scalar conservation law reads as [15]


Here, is the numerical solution and is the test function. The numerical flux is taken to be the global Lax-Friedrichs flux:


where and the maximum is taken over the whole domain. belong to the same function space If we recover the standard (or full grid) DG method. If we obtain the sparse grid DG method. In this paper, we will take as a subset of that is chosen adaptively as outlined in Section 3.2.

In DG methods, the integrals over elements and edges are often approximated by numerical quadrature rules on each cell [13]. However, in sparse grid DG method, this naive approach would result in computational cost that is proportional to the number of fundamental elements, i.e., , and is still subject to the curse of dimensionality. To evaluate the integrals over elements and edges more efficiently with a cost proportional to the DoF of the underlying finite element space, we interpolate the nonlinear function by using the multiresolution Lagrange (or Hermite) interpolation basis functions introduced in Section 2.2. Therefore, the semi-discrete DG scheme with interpolation is


where is a multiresolution interpolation operator onto some finite element space with the same multiresolution structure as but of polynomial degree . The choice of will be specified later, which plays important roles in numerical stability and accuracy. Note that the numerical flux is only defined at edges, thus it remains to clarify the meaning of the interpolation . Since we use the global Lax-Friedrichs flux (10), we have

due to the linearity of the interpolation operator . Therefore, we only need to obtain the interpolation and then read the value on two sides of the edges to obtain and . Now, we discuss about numerical implementation. First, we read the (derivative) values of , which is a linear combination of Alpert’s basis functions at the chosen interpolation points. Second, we calculate the (derivative) values of at these interpolation points. Last, we transfer the (derivative) values to coefficients of interpolation basis, by using the algorithm introduced in [49]. At this point, the numerical integrations can be performed through a fast matrix-vector product as in [46]. We remark that the computational cost does not increase too much compared to the multiresolution DG schemes for linear equations introduced in [26]. The cost of the transformation from the (derivative) values to hierarchical coefficients is only linearly dependent on the dimension [49].

Now we discuss the choice of To preserve the accuracy of the original DG scheme (9), it is required that the interpolation operator reaches certain accuracy. Following [13], we rewrite the weak formulation (11) in the ODE form as


where is an operator onto which is a discrete approximation of and satisfies


To illustrate the ideas, we only consider the full grid or sparse grid DG methods, i.e. or For adaptive methods, similar intuitive arguments can be made, but rigorous proof is much harder. Using similar error estimates techniques in [13, 34], we have the following proposition on local truncation error:

Proposition 3.1 (Accuracy of semi-discrete DG scheme with interpolation).

Assume that the DG finite element space (standard or sparse) has polynomials up to degree if the interpolation operator in (11) has the accuracy of (standard) or (sparse) for sufficiently smooth functions, then the truncation error of the semi-discrete DG scheme with interpolation (11) is of order (standard) or (sparse). To be more precise, for sufficiently smooth function , the standard DG with interpolation (11) has the truncation error:


and the sparse grid DG with interpolation (11) has the truncation error:


Here, the constant may depend on the solution, but does not depend on .


To save space, we only show the proof for full grid DG space . Similar technique also applies to the sparse grid DG space using projection error estimates in [25].

We denote the standard projection operator onto the standard DG finite element space by , then




The estimate for is trivial using projection properties:


To estimate we consider any test function in DG space, and obtain

Here we use the multiplicative trace inequality and the inverse inequality, see e.g. Lemma 2.1 and Lemma 2.3 in [34]. We take to be in the inequality above and have

and eventually arrive at


Combining (18) and (17), we have the estimate for the truncation error (14). ∎

From the proposition above, we find that, one has to apply the interpolation in which the function space has polynomials of one degree higher than the original DG function space, if one would like to preserve the order of accuracy for the original (standard or sparse grid) DG scheme, i.e. we shall require For example, if we take quadratic polynomials for the DG space, then it is required to apply interpolation operator (Lagrange or Hermite interpolation) to treat the nonlinear terms. From our numerical test, it seems that it is not a necessary condition for the standard DG method, but it is necessary for the sparse grid DG method.

In Proposition 3.1, we only estimate the truncation error, and this is far from rigorous error estimate that takes into account stability. In numerical experiments, we observe that the standard DG is stable with the Lagrange interpolation. However, the sparse grid DG with Lagrange interpolation is weakly unstable and will blow up with very fine mesh for polynomials of high degrees (see the numerical results in Table 2 and Table 3 in Section 4). With Hermite interpolation, the sparse grid DG scheme is more stable and produce satisfactory convergence rate (see Table 4 in Section 4). Actually, for standard DG with quadrature rules applied in each element, if the truncation error satisfies the required order of accuracy, then the convergence and error estimate can be guanranteed [34]. However, it is not true for the sparse grid DG method from our numerical experiments. This indicates that the standard DG method is more stable than the sparse grid DG method in this sense. We also remark that, since the interpolation operator introduced here is global but not local, the approach in [34] would probably fail to obtain the rigorous error estimate here. We will leave the detailed analysis as future work.

3.2 Adaptivity

In this section, we review the adaptive procedure introduced in [26] to determine the space The method is very similar to those in [26], except that two sets of basis functions are involved and they are adaptively chosen at the same time.

In the adaptive DG algorithm, we specify the maximum mesh level and an accuracy threshold . Details on the optimal choice of can be found for example in [33]. The same adaptive multiresolution projection method in [26] is applied here as the numerical initial condition for DG schemes. The error indicator using norm is used. The details are omitted and we refer readers to Algorithm 1 in [26].

The scheme is implemented by hash table as the underlying data structure. We now introduce the concepts of child, parent and leaf elements. If an element with satisfies the condition that there exists an integer such that and , where denotes the unit vector in the direction, and the support of is within that of , then is called a child element of . Accordingly, element is called a parent element of . If an element does not have its child element in the hash table, then we call it a leaf element.

The time evolution consists of four steps. The first step is the prediction step, which means given the hash table that stores the numerical solution at time step and the associated leaf table , we need to predict the location where the details becomes significant at the next time step , then add more elements in order to capture the fine structures. We solve for from to using a cheap solver, e.g. the forward Euler discretization. Here, the interpolation operator is determined by accuracy requirement, and has the same multiresolution structure as determined by the hash table corresponding to the numerical solution . The predicted solution at is denoted by . Note that to save cost, that the artificial viscosity term as introduced in Section 3.3 does not need to be included in the prediction step.

The second step is the refinement step according to the predicted solution . We traverse the hash table and if an element satisfies the refinement criteria


where denotes the hierarchical coefficient corresponding to the basis i.e. (19) indicates that such an element becomes significant at the next time step, then we need to refine the mesh by adding its children elements to . The detailed procedure is described as follows. For a child element of , if it has been already added to , i.e. , we do nothing; if not, we add the element to and let the associated detail coefficients . Moreover, we need to increase the number of children by one for all elements that has as its child element and remove the parent elements of from the leaf table if they have been added. Finally, we obtain a larger hash table and the associated approximation space and the leaf table .

Then, based on the updated hash table , we evolve the numerical solution by the DG formulation with space . Namely, we solve for from to , to generate the precoarsened solution , by using the the accurate solver with artificial viscosity in Section 3.3. Here, the interpolation operator should determined by the updated hash table . Note that in the artificial viscosity we fix to be such that the matrix for the diffusion term only need to be resembled one time in each time step.

The last step is to coarsen by removing elements that become insignificant at time level The hash table that stores the numerical solution is recursively coarsened by the following procedure. The leaf table is traversed, and if an element satisfies the coarsening criterion


where is a prescribed error constant, then we remove the element from both table and , and set the associated coefficients . For each of its parent elements in table , we decrease the number of children by one. If the number becomes zero, i.e, the element has no child any more, then it is added to the leaf table accordingly. Repeat the coarsening procedure until no element can be removed from the table . By removing only the leaf element at each time, we avoid generating “holes” in the hash table. The output of this coarsening procedure are the updated hash table and leaf table, denoted by and respectively, and the compressed numerical solution . In practice, is chosen to be smaller than for safety. In the simulations presented in this paper, we use .

3.3 Artificial viscosity

For capturing shock, we add artificial viscosity following the approach in [6] and arrive at the semi-discrete formulation


where is the artificial viscosity. The artificial viscosity is piecewise constant in each element and depends on . It is only imposed in the leaf element (since the sharp gradient and shock will only appear in the leaf element) and is determined in the following approach:

where and are constants chosen empirically. In the numerical experiments, we typically take and . The parameters and are given as


In the regions where the solutions are smooth, should be the same order as . If the solution is discontinuous, should be much larger than .

This smoothness indicator is motivated by the estimate of the coefficients of hierarchical basis functions in [25, 26]. It is shown in [25, 26] that for a function ,


where and is a constant independent of mesh level . Therefore, by assuming that , we can obtain a local estimate in each element on mesh level : for any index ,


Therefore, for sufficiently smooth functions, the coefficients should decay like

Remark 3.1.

There are still many problems to be explored on the artificial viscosity. The first one is the specific form of the artificial viscosity term. Here, for simplicity, we only add an artificial viscosity term in (3.3). One may also add a physical diffusion term and then discrete it using local DG [42] or interior penalty DG [37]. The second issue is how to choose the optimal parameters and in the artificial viscosity to obtain a sharp shock profile. The artificial neural network introduced in [44, 19] might be helpful for this problem. We will explore these subjects in future work.

The diffusion coefficient is of order for trouble cells and zero for normal cells. Thus, the explicit time integration in both convection and diffusion terms in (3.3) will yield CFL condition . For hyperbolic problems with DG discretizations using polynomials of degree and upwind numerical flux and a stage explicit RK method of order , the CFL constant is around [16]. However, for solving diffusion equation with local DG discretization with polynomials and alternating numerical flux, the CFL constant is around 0.0555 for , 0.0169 for , 0.0063 for , and 0.003 for , if coupled with explicit Runge-Kutta methods of the corresponding order111The CFL constants are provided by Chi-Wang Shu from Brown University in personal communications., which is much smaller than the CFL constant for convection terms, especially for polynomials of high degrees. If the alternating numerical flux is replaced by the central flux for the diffusion equation, the CFL constant is slightly larger but still much smaller than the CFL constant for the convection part: 0.125 for , 0.0384 for , 0.0158 for and 0.0083 for .

To obtain better computational efficiency, we avoid explicit time integrations and apply the IMEX time discretizations where the convection term is treated explicitly and the diffusion term implicitly. Here, we only present the third-order IMEX method introduced in [41], which will be coupled with the DG space of quadratic polynomials. The explicit part is the same with the explicit third-order strong stability preserving (SSP) Runge-Kutta method [47] and the implicit part has four stages. To be precise, for the ODE systems:


where denotes the non-stiff term (convection parts) and the stiff term (diffusion parts). The IMEX scheme for (26) reads as


with the stage and the parameters


The other parameters not listed above are zero.

By using the IMEX time integrator, the time step restriction remains the same as determined by the convection term. Note that the artificial viscosity is determined by and will keep unchanged in the middle stages of time evolution from to . Therefore, the matrix for the diffusion term only need to be resembled one time in each time step. Also, we only need to solve a linear system in which the coefficient matrix is symmetric positive definite and also sparse (there exist only a small portion of elements with non-zero viscosity). In the computation, we apply the conjugate gradient method to solve this linear system. We also remark that, for smooth solutions, this scheme will reduce to the explicit time integrations when coupled with the semi-discrete DG scheme with artificial viscosity (3.3), since the artificial viscosity will automatically vanish and then IMEX scheme (27) reduces to the third-order SSP RK method.

4 Numerical results

In this section, we perform numerical experiments to validate the accuracy and robustness of our scheme. The computational domain is for 1D and for 2D. Periodic boundary condition is imposed. When testing accuracy for smooth solutions, we apply the TVD Runge-Kutta time discretizations [47]: second-order RK method for the piecewise linear finite element space () and third-order RK method for the quadratic () and cubic () finite element space. When testing the capability for capturing discontinuous solutions, we use the quadratic finite element space (